Trigonometric regression

[Marc Jarczok]

Hi stata users,

i´m using sine and cosine terms as predictors in modeling periodic time series.
My aim is to extract the MESOR, Amplitude and Acrophase of an individual 24hour recording.

For those interested in the contend / meaning of the data: The time series consists of RMSSD values (an expression of ones heart rate variability) of a single individual collapsed to ever 5.35minute intervals but with gaps (e.g. if the recording went not very well during a 5.35 min interval). The recording started at 10:37 for 24hours.


So I ran the following syntax on my 288 data points (4 missing resulting in 284 analyzed segments)

.         reg rmssd  sintimeA costimeA
 
      Source |       SS       df       MS              Number of obs =     284
-------------+------------------------------           F(  2,   281) =   38.24
       Model |  30413.6508     2  15206.8254           Prob > F      =  0.0000
    Residual |  111745.603   281  397.671186           R-squared     =  0.2139
-------------+------------------------------           Adj R-squared =  0.2083
       Total |  142159.254   283   502.32952           Root MSE      =  19.942
 
------------------------------------------------------------------------------
       rmssd |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    sintimeA |  -5.174174   1.769928    -2.92   0.004    -8.658175   -1.690173
    costimeA |  -13.42496   1.616504    -8.30   0.000    -16.60696   -10.24297
       _cons |   29.10342   1.201788    24.22   0.000     26.73777    31.46907
------------------------------------------------------------------------------
 
.         predict yhatA ,   xb
 
.         gen mesorA=_b[_cons]
 
.         gen ampA = sqrt(_b[costimeA]^2 + _b[sintimeA]^2)
 
.         gen acrA = atan(-_b[sintimeA] / _b[costimeA])
 
.         list ampA mesorA acrA in 1/1
 
     +---------------------------------+
     |     ampA     mesorA        acrA |
     |---------------------------------|
  1. | 14.38755   29.10342   -.3678697 |
     +---------------------------------+

But when inspecting the results in a graph using the "twoway function" graph command one can see that the acrophase of the long-dashed curve (As is) is not corresponding to the yhat curve (solid line) (although MESOR and Amplitude do so).
In my example here, the appropriate acrophase value would be -10.6 (as opposed to -.036787) to match the yhat curve (short dashed line labeled "Expected" bc i expected the function to fit the yhat line).
What is it that I am missing here? Do i have an error in calculating the Acrophase value?

In addition, one can see that the data contain ultradian rhythms, which can be captured by e.g. doubling the amount of cycles.
But how would the amplitude & acrophase be calculated in that case?

Thanks for your help & time!


Marc

[Clyde Schechter]

I think there is something wrong with your expectations about acrophase. If you want it in radians, as your formula based on the atan() function suggests, it will always be between -pi and +pi, so 10.6 is not a possible value in any case. If you want it in time units, then you have to take the value you got from your formula and transform it into the time metric. Presumably in calculating your sine and cosine terms you transformed clock time to something like 2pi*(time-start_time)/24 first. So you have to apply the inverse of that transformation to your result to get it back to the time metric.

Turning to your coding practices (an aside): it rarely makes sense to create variables to only hold constants. So unless your -gen- statements are simply the first stage in some calculations that will introduce variation in later steps, you probably should store mesorA, ampA, and acrA as scalars rather than variables.

[Nick Cox]

I think you have a stray minus sign in there.

I get from differentiating a curve like a + s sin(theta) + c cos(theta) that tangent of phase = s/c and given that 24 hours appear to correspond to 2 pi radians, your peak is at about

 
. di (12 * _pi) *  atan(5.174174/13.42496)
13.868364

hours, which matches the graph. Here I cancelled the minus signs by eye.

[Marc Jarczok]

Dear Clyde, Dear Nick, thank you so much for your response.
@ Clyde: The example is from a data processing loop where several thousand recordings are processed and exported as new datasets with one line per individual holding the id and the newly generated vars. I was aware that 10.6 is not a possible values in terms of radiants, but this is the values that makes the function correspond to the yhat function.

@ Nick: I´ll check & correct my formulas when i´m back at work. Can you tell me how to calculate Acrophase & Amplitude when using several cycles (as suggested in your SJ Article "In praise of trigonometric predictors " (2006)?

[Nick Cox]

I'd have to check the definitions of those quantities. It is easily enough to search for maxima numerically by applying summarize to fitted curves and looking up the times at which they occur.

[Marc Jarczok]

Dear all, after a pause i had time to go back to this topic. I found the error, it´s in the computation & expression of the acrophase and the time scale definition of sintimeA and costimeA -

I defined sintimeA and costimeA by 24hours
gen double sintimeA= sin( 2*_pi*datetime/24)
gen double costimeA= cos( 2*_pi*datetime/24)

instead of 24h my divisor should have been 288 (my amount of 5.35 minute intervals)
The acrophase needs to be converted from radiants to the according time unit i.e.:

gen acrA = ((atan(-_b[sintimeA] / _b[costimeA]) ) * 180/_pi ) * (288/360)


this gives the correct result. Now it´s pretty obvious to me :D

Thanks again for your help & time.
Best Marc