R/cosinor.R
In frenchja/psych: Procedures for Psychological, Psychometric, and Personality Research
#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)
}
"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.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
x <- cbind(angle,x)} #convert to radians
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) }
frenchja/psych documentation built on May 16, 2019, 2:49 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#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)
}
"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.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
x <- cbind(angle,x)} #convert to radians
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) }
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.