I have ecological data, and a physiological response variable. I'm suspicious there are circadian rhythms at play, but the current COSINOR models are a bit heavy for what I'm trying to do. What is the best way to handle day of the year in a glm (if there is one!) that accounts for its cyclical nature (ie 1 = 366)?
There may be a better way to do it, but you can transform your JDay variable to create pairs of sin and cos harmonics, which become a Fourier series, i.e.
s1 = sin(2 * pi * JDay / 366)
c1 = cos(2 * pi * JDay / 366)
The result should be that s1 are c1 cyclic from -1 to 1 as they make their way around the unit circle (in y and x coordinates respectively).
To add them into a model, just add in your linear term of your variable of interest and interact your variable of interest with each of the harmonics, which will give you a coefficient for the linear and each interaction term, and you can recreate the temporally-varying function by subbing these coefficients back into your model formula.
So if your model formula was (assuming cov1 is what we think changes throughout the year) (sorry for the mess...):
y ~ a + (beta1 * cov1) + (beta2 * cov1 * s1) + (beta3 * cov1 * c1) + ... etc
Then you can get your temporally varying coefficient function by subbing those three coefficients (betas) into:
beta1 +
beta2 * sin( 2 * pi * JDay / 366) +
beta3 * cos(2 * pi * JDay / 366)
Which will give you a sequence of coefficients throughout the year (which will be cyclic; 1 = 366).
It should look something like the image that I attached (from using 4 harmonics), which was an animal resource selection coefficient which varied throughout the day. If there is no temporal change to the coefficient then it should just be flat throughout the year. You will just have a single function which is the mean, and to get the uncertainty estimates like I have there you have to bootstrap by taking beta1s, beta2s, etc from the normal distribution given to you by the glm (i.e. the mean and sd of each coefficient in the model output), and for each set of betas you can create a new function. I did this for 100 different sets of betas to create 100 functions. A bit fiddly but not so bad to iterate in R.
To get more rapidly varying responses you can just keep adding in harmonics in pairs of sin and cos whilst decreasing the period (theoretically your results should get more accurate as you add more harmonics, as your function can approximate more functions, but may not be necessary, try just two first).
s2 = sin(4 * pi * JDay / 366)
c2 = cos(4 * pi * JDay / 366)
Hopefully that makes at least some sense. It's not difficult to implement once you get your head around thinking about it, but feel free to message or email me to get some more info (or some messy code).
- Badges
- Science topic
- Similar topics
- Geoscience
- Asia
- India