#functions for doing diurnal rhythm analyses #Heavily adapted from the circular package # #a function to estimate diurnal phase of mood data #the input is a data frame or matrix with #time of measurement (in 24 hour clock) #and then the mood measures (1 or many) #Version of October 22, 2008 #seriously revised April 12, 2009 # #modified February 14, 2015 to adjust the call to mean(x[1]) #and to make the grouping function actually work #find the best fitting phase (in hours) #cleaned up March 9, 2015 to allow a more natural calling sequence #modifed Jan - April, 2016 to make cleaner code #Added cosinor.period April 21, 2016 to iteratively fit period as an option "cosinor.period" <- function(angle,x=NULL,code=NULL,data = NULL, hours=TRUE,period=seq(23,26,1),plot=FALSE,opti=FALSE,na.rm=TRUE) { #first, organize the data in terms of the input if(!is.null(data)) { if(is.matrix(data)) data <- data.frame(data) if(is.character(angle)) angle <- which(colnames(data) == angle) if(!is.null(code)) { if(is.character(code)) codeloc <- which(colnames(data) ==code) x <- data[,c(angle,x,codeloc)] } else {x <- data[,c(angle,x)]} angle <- x[1] x <- x[-1] } else { if (is.null(x) && is.null(code)) {angle <- data.frame(angle) x <- angle angle<- angle[,1] } else {x <- data.frame(x) x <- cbind(angle,x) angle <- x[1] x <- x[-1] } } xdata <- x #we need to save this to do iterative fits old.angle <- angle per <- period fit.period <- list() for (i in 1:length(per)) { period <- per[i] if(hours) { angle <- old.angle*2*pi/period x <- cbind(angle,xdata)} #convert to radians nvar <- dim(xdata)[2] -1 if(is.null(code)) { fit <- cosinor1(angle,x[-1],period=period,opti=opti,na.rm=na.rm) #if there is a code (i.e., subject id), then do the fits for each separate subject m.resp <- mean(x[,1]) s.resp <- sd(x[,1])} else { #x <- angle fit.list <- by(x,x[code],function(x) cosinor1(angle=x[1],x=x[-c(1,which(colnames(x)== code))],period=period,opti=opti,na.rm=na.rm)) #this is the case if code is specified ncases <- length(fit.list) fit <- matrix(unlist(fit.list),nrow=ncases,byrow=TRUE) colnames(fit) <- c(paste(colnames(x)[-c(1,which(colnames(x)== code))], "phase",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "fit",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "amp",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "sd",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "mean",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "intercept",sep=".")) rownames(fit) <- names(fit.list) } fit.period[[i]] <- list(fit) } x <- NA #just to avoid being told that x is a global #now, find for each variable and for each subject, that value of of fit which is maximized, and then what is the ncols <- 6 * length(x) fit.m <- matrix(unlist(fit.period),nrow=ncases,byrow=FALSE) #the fits are every nvar * 6 elements starting at nvar + 1 maxfit <- per np <- length(per) fits <- cbind(matrix(NA,nrow=ncases,ncol=nvar),fit) for (j in 1:ncases) { #do it for each subject for (i in 1:nvar) {#do it for each variable for (p in 1:np) {#find the fits for all periods #maxfit[p] <- fit.m[j,(p-1) * nvar * 6 + nvar+1] maxfit[p] <- fit.period[[p]][[1]][j,i+nvar] } max.period <- which.max(maxfit) fits[j,i] <- per[max.period] fits[j,i+nvar] <- fit.period[[max.period]][[1]][j,i] fits[j,i+2*nvar] <- fit.period[[max.period]][[1]][j,i+nvar] fits[j,i+3*nvar] <- fit.period[[max.period]][[1]][j,i + 2*nvar] fits[j,i+4*nvar] <- fit.period[[max.period]][[1]][j,i+3* nvar] fits[j,i+5*nvar] <- fit.period[[max.period]][[1]][j,i + 4*nvar] fits[j,i+6*nvar] <- fit.period[[max.period]][[1]][j,i+5* nvar] } } return(fits) } #revised 9/15/18 to handle radians correctly and to handle character names for variables "circadian.phase" <- "cosinor" <- function(angle,x=NULL,code=NULL,data = NULL, hours=TRUE,period=24,plot=FALSE,opti=FALSE,na.rm=TRUE) { #first, organize the data in terms of the input if(!is.null(data)) { if(is.matrix(data)) data <- data.frame(data) if(is.character(angle)) angle <- which(colnames(data) == angle) if(is.character(x)) x <- which(colnames(data) ==x) if(!is.null(code)) { if(is.character(code)) codeloc <- which(colnames(data) ==code) x <- data[,c(angle,x,codeloc)] } else {x <- data[,c(angle,x)]} angle <- x[1] x <- x[-1] } else { if (is.null(x) && is.null(code)) {angle <- data.frame(angle) x <- angle angle<- angle[,1] } else {x <- data.frame(x) x <- cbind(angle,x) angle <- x[1] x <- x[-1] } } if(hours) { angle <- angle*2*pi/period } #convert to radians x <- cbind(angle,x) nvar <- dim(x)[2] if(is.null(code)) { fit <- cosinor1(angle,x[-1],period=period,opti=opti,na.rm=na.rm) #if there is a code (i.e., subject id), then do the fits for each separate subject m.resp <- mean(x[,1]) s.resp <- sd(x[,1]) if(plot) {#curve(cos((*x-fit[1,1])*s.resp+m.resp)*pi/12),add=TRUE) } #this draws the first fitted variable }} else {#x <- angle fit.list <- by(x,x[code],function(x) cosinor1(angle=x[1],x=x[-c(1,which(colnames(x)== code))],period=period,opti=opti,na.rm=na.rm)) #this is the case if code is specified ncases <- length(fit.list) fit<- matrix(unlist(fit.list),nrow=ncases,byrow=TRUE) colnames(fit) <- c(paste(colnames(x)[-c(1,which(colnames(x)== code))], "phase",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "fit",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "amp",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "sd",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "mean",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "intercept",sep=".")) rownames(fit) <- names(fit.list) } x <- NA #just to avoid being told that x is a global return(fit) } # cosinor1 actually does the work # it either uses a fitting function (optimize) from core R # or calls a linear regression fit "cosinor1" <- function(angle,x,period=24,opti=FALSE,na.rm=TRUE) { response <- x n.var <- dim(x)[2] if(is.null(n.var)) n.var <-1 fit <- matrix(NaN,nrow=n.var,ncol=6) for (i in 1:n.var) { if(opti) {fits <- optimize(f=phaser,c(0,24),time=angle,response=x[,i],period=period,maximum=TRUE) #iterative fit fit[i,1] <- fits$maximum fit[i,2] <- fits$objective } else {fits <- cosinor.lm2 (angle,x[,i],period=period,na.rm=na.rm) #simple linear regression based upon sine and cosine of time fit[i,1] <- fits[[1]] #this is the acrophase fit[i,2] <- fits[[2]] #this is the correlation of the fit with the data fit[i,3] <- fits[[3]] #this is the amplitude fit[i,4] <- fits[[4]] #The standard deviation of the observed scores fit[i,5] <- fits[[5]] #the mean of the observed scores fit[i,6] <- fits[[6]] #the predicted value of the observed score at the intercept } } colnames(fit) <- c("phase","fit","amplitude","sd","mean","intercept") rownames(fit) <- colnames(x) return(fit) } "phaser" <- function(phase,time,response,period) { #this is used in the iterative fit procedure phaser <- cor(cos(((time-phase)*2*pi)/period),response,use="pairwise")} #the alternative to the iterative procedure is simple linear regression of the cosine and sine "cosinor.lm2" <- function(angle,y,period=24,na.rm=na.rm) { p2 <- period/2 cos.t <- cos(angle) #angle is already in radians! sin.t <- sin(angle) dat.df <- data.frame(iv1=cos.t,iv2=sin.t,y) cor.dat <- cor(dat.df,use="pairwise") beta1 <- (cor.dat[3,1] - cor.dat[3,2] * cor.dat[1,2])/(cor.dat[1,1]-cor.dat[1,2]^2) beta2 <- (cor.dat[3,2] - cor.dat[3,1] * cor.dat[1,2])/(cor.dat[1,1]-cor.dat[1,2]^2) #note, these are standardized beta weights # phase <- ( sign(beta2) *acos( beta1/sqrt(beta1^2 + beta2^2)))*p2/pi #this is in hours phase <- atan2(beta2,beta1) #in radians intercept <- cos(phase) #the value at time 0 phase <- phase *p2/pi #convert to hours phase <- phase %% period r <- cor(cos(angle-phase*pi/p2),y,use="pairwise") sdy <- sd(y,na.rm=na.rm) meany <- mean(y,na.rm=na.rm) #amp <- r *sdy/.7223 amp <- sqrt(beta1^2 + beta2^2) #see Chow 2009 among others -- note we are finding the standardized amp intercept <- intercept * amp * sdy + meany #amp <- r * sd(y,na.rm=TRUE)/sd(cos(angle[,1]),na.rm=TRUE) #R <- sqrt(cor.dat[3,1]*beta1 + cor.dat[3,2]*beta2) #these are identical fit <- list(phase=phase,R=r,amp=amp,sd=sdy,mean=meany,intercept) return(fit) } "cosinor.plot" <- function(angle,x=NULL,data = NULL,IDloc=NULL,ID=NULL,hours=TRUE,period=24,na.rm=TRUE,ylim=NULL,ylab="observed",xlab="Time (double plotted)",main="Cosine fit",add=FALSE,multi=FALSE,typ="l",...) { if(!multi) {main <- paste("ID = ",ID," ",x)} if(!is.null(data)) { if(is.matrix(data)) data <- data.frame(data) # if(is.character(angle)) angle <- which(colnames(data) == angle) if(!is.null(IDloc)) { x <- data[data[,IDloc]==ID,c(angle,x)] angle <- x[,1,drop=FALSE] } else {x <- data[,c(angle,x)] angle <- x[,1,drop=FALSE] main <- c(main," ",IDloc)} } else { if (is.null(x) && is.null(data)) {x <- data.frame(x) x <- angle angle<- angle[1] } else {x <- data.frame(x) x <- cbind(angle,x) angle <- x[1] x <- x[-1] } } if(hours) { angle <- angle*2*pi/24 } xx <- 1:96 fit <- cosinor1(angle, x = x[2], period = period, na.rm = na.rm) m.resp <- mean(x[, 2], na.rm = TRUE) s.resp <- sd(x[, 2], na.rm = TRUE) sd.time <- sd(cos(angle[, 1]), na.rm = TRUE) if(missing(ylim)) ylim=c(min(x[,2],(m.resp - fit[1,3]),na.rm=na.rm),max(x[,2],(m.resp + fit[1,3]),na.rm=na.rm)) if(!multi | !missing(main)){main <- paste(main," ",round(fit[1,1],2)) } else {main=main} # plot(xx/2,cos((xx/2-fit[1,1])*pi/12)*s.resp*fit[1,2]/.707+ m.resp,typ=typ,ylim=ylim,...) # plot(xx/2,cos((xx/2-fit[1,1])*pi/12)*fit[1,3]*s.resp+ #m.resp,typ="l",ylim=ylim,ylab=ylab,xlab=xlab,main=main,...) #plot the lines first if(!add) {plot(xx/2,cos((xx/2-fit[1,1])*pi/12)*s.resp*fit[1,2]/.707+ m.resp,typ=typ,ylim=ylim,ylab=ylab,xlab=xlab,main=main,...) } else { points(xx/2,cos((xx/2-fit[1,1])*pi/12)*s.resp*fit[1,2]/.707+ m.resp,typ=typ,ylim=ylim,ylab=ylab,xlab=xlab,main=main,...) } if(!multi) { points(c(x[,1],x[,1] + 24),rep(x[,2],2),...)} else {points(xx/2,cos((xx/2-fit[1,1])*pi/12)*s.resp*fit[1,2]/.707+ m.resp,typ="l",...)} #this draws the first fitted variable } #Added March 26, 2015 to do split half (first/second) reliabilities "circadian.reliability" <- function(angle,x=NULL,code=NULL,data = NULL,min=16, oddeven=FALSE,hours=TRUE,period=24,plot=FALSE,opti=FALSE,na.rm=TRUE) { cl <- match.call() if(!is.null(data)) { if(is.character(angle)) angle <- which(colnames(data) == angle) if(!is.null(code)) { if(is.character(code)) codeloc <- which(colnames(data) ==code) x <- data[,c(angle,x,codeloc)] } else {x <- data[,c(angle,x)]} angle <- x[1] x <- x[-1] } else { if (is.null(x) && is.null(code)) {x <- angle angle<- angle[,1] } else {x <- cbind(angle,x) angle <- x[1] x <- x[-1] } } if(hours) { angle <- angle*2*pi/period x <- cbind(angle,x) } n.obs <- dim(x)[1] if(is.null(code)) { fit <- cosinor.rel(angle,x,period=period,na.rm=na.rm) #if there is a code (i.e., subject id), then do the fits for each separate subject m.resp <- mean(x[,1]) s.resp <- sd(x[,1]) } else { fit.list <- by(x,x[,code],function(x) cosinor.rel(angle=x[1],x=x[-c(1,which(colnames(x)== code))],min=min,oddeven=oddeven,na.rm=na.rm)) #this is the case if code is specified ncases <- length(fit.list) fit <- matrix(unlist(fit.list),nrow=ncases,byrow=TRUE) colnames(fit) <- c(paste(colnames(x)[-c(1,which(colnames(x)== code))], "phase1",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "phase2",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "fit1",sep="."),paste(colnames(x)[-c(1,which(colnames(x)== code))], "fit2",sep=".")) rownames(fit) <- names(fit.list) } nvar <-ncol(fit)/4 r <- circadian.cor(fit[,1:(nvar*2)]) r.fit <- cor(fit[,(nvar*2+1):ncol(fit)],use="pairwise") splithalf <- split.fit <- rep(NA,nvar) for (i in 1:nvar) {splithalf[i] <- r[i,(nvar+i)] split.fit[i] <- r.fit[i,(nvar+i)]} rel <- splithalf * 2/(1+splithalf) fit.rel <- split.fit * 2/(1+split.fit) names(rel) <- paste(colnames(x)[-c(1,which(colnames(x)== code))]) names(fit.rel) <- paste(colnames(x)[-c(1,which(colnames(x)== code))]) # x <- NA #just to avoid being told that x is a global #now do the F test between the two splits split.F <- circadian.split.F(fit) result <- list(phase.rel=rel,fit.rel=fit.rel,split.F = split.F, splits=fit,Call=cl) class(result) <- c("psych","circadian","reliability") return(result) } # cosinor.rel actually does the work # or calls a linear regression fit "cosinor.rel" <- function(angle,x,min=16,oddeven=FALSE,period=24,na.rm=TRUE) { response <- x n.var <- dim(x)[2] n.obs <- dim(x)[1] if(is.null(n.var)) n.var <-1 fit <- matrix(NaN,nrow=n.var,ncol=4) if(n.obs >= min) { for (i in 1:n.var) { if(oddeven) {fits1 <- cosinor.lm2 (angle[seq(1,n.obs,2),1],x[seq(1,n.obs,2),i],na.rm=na.rm) fits2 <- cosinor.lm2 (angle[seq(2,n.obs,2),1],x[seq(2,n.obs,2),i],na.rm=na.rm) } else { fits1 <- cosinor.lm2 (angle[1:n.obs/2,1],x[1:n.obs/2,i],na.rm=na.rm) fits2 <- cosinor.lm2 (angle[(n.obs/2+1):n.obs,1],x[(n.obs/2+1):n.obs,i],na.rm=na.rm) #simple linear regression based upon sine and cosine of time} #simple linear regression based upon sin and cosine of time } fit[i,1] <- fits1[[1]] #this is the acrophase fit[i,3] <- fits1[[2]] #this is the correlation of the fit with the data fit[i,2] <- fits2[[1]] #this is the acrophase fit[i,4] <- fits2[[2]] #this is the correlation of the fit with the data } } colnames(fit) <- c("phase","phase2","fit","fit2") rownames(fit) <- colnames(x) return(fit) } "circadian.split.F" <- function(angle,hours=TRUE,na.rm=TRUE) { nvar <- ncol(angle)/4 stats1 <- circadian.stats(angle[,1:nvar]) stats2 <- circadian.stats(angle[,(nvar+1):(nvar*2)]) pool <- rbind(angle[,1:nvar],angle[,(nvar+1):(nvar*2)]) all <- circadian.stats(pool) allR <- all$n * all$R within <- matrix(c(stats1$n*stats1$R,stats2$n*stats2$R),ncol=2) rownames(within) <- rownames(all) ngroups <- 2 SSb <- rowSums(within) - allR SSw <- all$n - rowSums(within) dfw <- all$n - ngroups MSw <- SSw/dfw dfb = ngroups -1 MSb <- SSb/dfb F <- MSb/MSw prob <- 1-pf(F,dfb,dfw) F.df <- data.frame(SSb= SSb,dfb=dfb,MSb=MSb,SSw=SSw,dfw=dfw,MSw=MSw,F=F,prob=prob) result<- list(pooled =all,group1 =stats1, group2=stats2 ,F=F.df) class(result) <- c("psych","circadian") return(result) } ## # #find the mean phase of output from cosiner or any other circular data set #can find the mean phase of data in radians or hours (default) # "circadian.mean" <- function(angle,data=NULL,hours=TRUE,na.rm=TRUE) { if(!is.null(data)) angle <- data[,angle] if(hours) { angle <- angle*2*pi/24 } x <- cos(angle) y <- sin(angle) if (is.vector(angle)) { mx <- mean(x,na.rm=na.rm) my <- mean(y,na.rm=na.rm) } else { mx <- colMeans(x,na.rm=na.rm) my <- colMeans(y,na.rm=na.rm) } mean.angle <- sign(my) * acos((mx)/sqrt(mx^2+my^2)) # mean.angle <- atan(my/mx) #according to circular stats, but the other form is clearer if (hours) {mean.angle <- mean.angle*24/(2*pi) mean.angle[mean.angle <= 0] <- mean.angle[mean.angle<=0] + 24} return(mean.angle) } "circadian.sd" <- function(angle,data=NULL,hours=TRUE,na.rm=TRUE) { if(!is.null(data)) angle <- data[,angle] if(hours) { angle <- angle*2*pi/24 } nvar <- dim(angle)[2] if(is.null(nvar)) nvar <- 1 x <- cos(angle) y <- sin(angle) if(nvar > 1) { mx <- colSums(x,na.rm=na.rm) my <- colSums(y,na.rm=na.rm) n.obs <- colSums(!is.na(angle))} else { mx <- sum(x,na.rm=na.rm) my <- sum(y,na.rm=na.rm) n.obs <- sum(!is.na(angle))} R <- sqrt(mx^2+my^2)/n.obs mean.angle <- sign(my) * acos((mx/n.obs)/R) Rvar <- 1 - R sd <- sqrt(-2 * log(R)) #for (i in 1:nvar) {#the logic is that larger deviations are weighted more, up to the sin(theta) # var[i] <- sum(sin(angle[,i] -mean.angle[i])^2 ,na.rm=na.rm) } #n.obs <- colSums(!is.na(angle)) ##but these are in radians! if(hours) {#sd <- sd * 24/(pi*2) Rvar <- Rvar * 24/(pi*2)} return(list(Rvar=Rvar,sd =sd,R= R)) } "circadian.stats" <- function(angle,data=NULL,hours=TRUE,na.rm=TRUE) { cl <- match.call() means <- circadian.mean(angle=angle,data=data,hours=hours,na.rm=na.rm) csd <- circadian.sd(angle=angle,hours=hours,na.rm=na.rm) if(!is.null(data)) angle <- data[,angle] if(length(means)>1 ) { n.obs <- colSums(!is.na(angle))} else {n.obs <- sum(!is.na(angle)) } R <- csd$R if(hours) {sd <- csd$sd*24/(2*pi)} else {sd <- csd$sd} z <- n.obs * R^2 p <- exp(-z) result<- data.frame(n=n.obs,mean=means,sd=sd,R,z=z,p=p) #result <- list(n=n.obs,mean=means,sd=sd,R,z=z,p=p,call=cl) class(result) <- c("psych","circadian","data.frame") return(result) } "circadian.F" <- function(angle,group,data=NULL,hours=TRUE,na.rm=TRUE) { if(!is.null(data)) {angle <- data[,angle] group <- data[,group]} stats <- by(angle,group,circadian.stats) all <- circadian.stats(angle) allR <- all$n * all$R nR <- lapply(stats,function(x) x$n * x$R) ngroups <- length(nR) within <- matrix(unlist(nR),ncol=ngroups) rownames(within) <- rownames(all) SSb <- rowSums(within) - allR SSw <- all$n - rowSums(within) dfw <- all$n - ngroups MSw <- SSw/dfw dfb = ngroups -1 MSb <- SSb/dfb F <- MSb/MSw prob <- 1-pf(F,dfb,dfw) F.df <- data.frame(SSb= SSb,dfb=dfb,MSb=MSb,SSw=SSw,dfw=dfw,MSw=MSw,F=F,prob=prob) result<- list(pooled =all,bygroup = stats,F=F.df) class(result) <- c("psych","circadian") return(result) } print_circadian <- function(x,short=TRUE,digits=2) { if(!is.null(x$Call)) {cat("Call: ") print(x$Call)} cat("\nCircadian Statistics :\n") if(!is.null(x$F)) { cat("\nCircadian F test comparing groups :\n") print(round(x$F,digits)) if(short) cat("\n To see the pooled and group statistics, print with the short=FALSE option") } if(!is.null(x$pooled) && !short) { cat("\nThe pooled circadian statistics :\n") print( x$pooled)} if(!is.null(x$bygroup) && !short) {cat("\nThe circadian statistics by group:\n") print(x$bygroup)} #if(!is.null(x$result)) print(round(x$result,digits)) if(!is.null(x$phase.rel)) { cat("\nSplit half reliabilities are split half correlations adjusted for test length\n") x.df <- data.frame(phase=x$phase.rel,fits=x$fit.rel) print(round(x.df,digits)) } if(is.data.frame(x)) {class(x) <- "data.frame" print(round(x,digits=digits)) } } ## The circular correlation matrix of phase data #adapted from the circStats package #with various modifications for the study of mood data #one change is not use atan but rather use cosine over length # "circadian.cor" <- function(angle,data=NULL,hours=TRUE,na.rm=TRUE) { if(!is.null(data)) angle <- data[,angle] if(hours) { angle <- angle*2*pi/24 } nvar <- dim(angle)[2] correl <- diag(nvar) x <- cos(angle) y <- sin(angle) mx <- colMeans(x,na.rm=na.rm) my <- colMeans(y,na.rm=na.rm) mean.angle <- sign(my) * acos((mx)/sqrt(mx^2+my^2)) for (i in 1:nvar) {#the logic is that larger deviations are weighted more, up to the sin(theta) for (j in 1:i) {covar <- sum(sin(angle[,i] -mean.angle[i]) *sin(angle[,j] -mean.angle[j]),na.rm=na.rm) correl[i,j] <- correl[j,i] <- covar} } var <- diag(correl)/colSums(!is.na(angle)) sd <- diag(sqrt(1/diag(correl))) correl <- sd %*% correl %*% sd colnames(correl) <- rownames(correl) <- colnames(angle) return(correl) } #to find the correlation of a linear variable (e.g., a personality trait) with a circular one (e.g., phase) "circadian.linear.cor" <- function(angle,x=NULL,data=NULL,hours=TRUE) { if(!is.null(data)) angle <- data[,angle] if(hours) { angle <- angle*2*pi/24 } if(is.null(x)) {x <- angle[2:dim(angle)[2]] angle <- angle[1]} cos.angle <- cos(angle) sin.angle <- sin(angle) cor.cos <- cor(cos.angle,x,use="pairwise") cor.sin <- cor(sin.angle,x,use="pairwise") if(!is.vector(angle)) {cor.cs <- diag(cor(cos.angle,sin.angle,use="pairwise"))} else {cor.cs <- cor(cos.angle,sin.angle,use="pairwise")} R <- sqrt((cor.cos^2 + cor.sin^2 - 2 * cor.cos * cor.sin * cor.cs)/(1-cor.cs^2))*sign(cor.cos) return(R) } "circadian.plot" <- function(angle,x=NULL,hours=TRUE,title="Polar plot") { if(hours) { angle <- angle*2*pi/24 } x1 <- cos(angle) * x y1 <- sin(angle) * x maxx <- max(x) plot(x1,y1,axes=FALSE,xlab="",ylab="",xlim=c(-maxx,maxx),ylim=c(-maxx,maxx),asp=1) segments <- 51 angles <- (0:segments) * 2 * pi/segments unit.circle <- cbind(cos(angles), sin(angles)) points(unit.circle*maxx,typ="l") text(maxx,0,"24",pos=4) text(-maxx,0,"12",pos=2) text(0,maxx,"6",pos=3) text(0,-maxx,"18",pos=1) } "circular.mean" <- function(angle,na.rm=TRUE) { x <- cos(angle) y <- sin(angle) if (is.vector(angle)) { mx <- mean(x,na.rm=na.rm) my <- mean(y,na.rm=na.rm) } else { mx <- colMeans(x,na.rm=na.rm) my <- colMeans(y,na.rm=na.rm) } mean.angle <- sign(my) * acos((mx)/sqrt(mx^2+my^2)) #mean.angle <- atan(my/mx) #according to circular stats, but the other form is clearer return(mean.angle) } "circular.cor" <- function(angle,na.rm=TRUE) { nvar <- dim(angle)[2] correl <- diag(nvar) x <- cos(angle) y <- sin(angle) mx <- colMeans(x,na.rm=na.rm) my <- colMeans(y,na.rm=na.rm) mean.angle <- sign(my) * acos((mx)/sqrt(mx^2+my^2)) for (i in 1:nvar) {#the logic is that larger deviations are weighted more, up to the sin(theta) for (j in 1:i) {covar <- sum(sin(angle[,i] -mean.angle[i]) *sin(angle[,j] -mean.angle[j])) correl[i,j] <- correl[j,i] <- covar} } var <- diag(correl) sd <- diag(sqrt(1/diag(correl))) correl <- sd %*% correl %*% sd colnames(correl) <- rownames(correl) <- colnames(angle) return(list(correl,var)) } #deprecated "circadian.linear.cor.2" <- function(angle,x,hours=TRUE) { if(hours) { angle <- angle*2*pi/24 } cos.angle <- cos(angle) sin.angle <- sin(angle) cor.cos <- cor(cos.angle,x,use="pairwise") cor.sin <- cor(sin.angle,x,use="pairwise") if(!is.vector(angle)) {cor.cs <- diag(cor(cos.angle,sin.angle))} else {cor.cs <- cor(cos.angle,sin.angle)} R <- sqrt((cor.cos^2 + cor.sin^2 - 2 * cor.cos * cor.sin * cor.cs)/(1-cor.cs^2)) return(R) }