Last chance! 50% off unlimited learning
Sale ends in
cosinor(angle,x=NULL,code=NULL,period=24,plot=FALSE,opti=FALSE)
circadian.mean(angle, hours=TRUE)
circular.mean(angle) #angles in radians
circadian.cor(angle,hours=TRUE) #angles in radians
circular.cor(angle) #angles in radians
circadian.linear.cor(angle,x,hours=TRUE)
circadian.cor
here. Just as the Pearson r is a ratio of covariance to the square root of the product of two variances, so is the circular correlation. The circular covariance of two circular vectors is defined as the average product of the sines of the deviations from the circular mean. The variance is thus the average squared sine of the angular deviations from the circular mean. Circular statistics are used for data that vary over a period (e.g., one day) or over directions (e.g., wind direction or bird flight). Jammalamadaka and Lund (2006) give a very good example of the use of circular statistics in calculating wind speed and direction. The code from CircStats and circular was adapted to allow for analysis of data from various studies of mood over the day.
circular.mean
and circular.cor
are just circadian.mean
and circadian.cor
but with input given in radians rather than hours.
The cosinor function will either iteratively fit cosines of the angle to the observed data (opti=TRUE) or use the circular by linear regression to estimate the best fitting phase angle. If cos.t <- cos(time) and sin.t = sin(time) (expressed in hours), then beta.c and beta.s may be found by regression and the phase is $sign(beta.c) * acos(beta.c/\sqrt(beta.c^2 + beta.s^2)) * 12/pi$
Simulations (see examples) suggest that with incomplete times, perhaps the optimization procedure yields slightly better fits with the correct phase than does the linear model, but the differences are very small. In the presence of noisey data, these advantages seem to reverse. The recommendation thus seems to be to use the linear model approach (the default).
time <- seq(1:24)
pure <- matrix(time,24,18)
pure <- cos((pure + col(pure))*pi/12)
matplot(pure,type="l",main="Pure circadian arousal rhythms",xlab="time of day",ylab="Arousal")
p <- cosinor(time,pure)
#set.seed(42)
noisy <- pure + rnorm(24*18)
n <- cosinor(time,noisy)
#small.pure <- pure[c(6:18),]
small.pure <- pure[c(8,11,14,17,20,23),]
#small.noisy <- noisy[c(6:18),]
small.noisy <- noisy[c(8,11,14,17,20,23),]
matplot(small.noisy,type="l",main="Noisy circadian arousal rhythms",
xlab="time of day",ylab="Arousal")
#sp <- cosinor(time[c(6:18)],small.pure) #linear fit
sp <- cosinor(time[c(8,11,14,17,20,23)],small.pure)
spo <- cosinor(time[c(8,11,14,17,20,23)],small.pure,opti=TRUE) #iterative fit
sn <- cosinor(time[c(8,11,14,17,20,23)],small.noisy) #linear
sno <- cosinor(time[c(8,11,14,17,20,23)],small.noisy,opti=TRUE) #iterative
sum.df <- data.frame(pure=p,noisy = n, small=sp,small.noise = sn,
small.opt=spo,small.noise.opt=sno)
round(sum.df,2)
round(circadian.cor(sum.df[,c(1,3,5,7,9,11)]),2) #compare alternatives
round(cor(sum.df[,c(2,4,6,8,10,12)]),2)
Run the code above in your browser using DataLab