[R] cosinor analysis
Charles C. Berry
cberry at tajo.ucsd.edu
Thu Jan 8 20:01:57 CET 2009
On Thu, 8 Jan 2009, Anne Berger wrote:
> Hallo,
> I didn´t found any facilities for Halbergs cosinor analysis in
R. This analysis is well known in the Chronobiology as the least
square approximation of time series using cosine function of known
period (in my case of 24hours-period). I tried to write a script but
crashed...
> Can you give me some advices, please!?
Anne,
I append a bunch of functions that implement the analysis of Y.L. Tong
(1976) Biometrics 32:85-94.
These are offered as is. To understand them, you will probably need to
refer to the Tong article. The notation I used in coding should come
close to matching that in Tong.
The following code should simulate 20 data for 20 subjects and fit a
cosinor function to each of them:
source("cosinor.R") ## file of all the functions below
X <- seq(0,24,by=2)
junk <- sim.cosinor(20,1,.1,1,.1,X,0,.1,.1, pi/12)
cosinor.lm.each( y ~ r.ij + s.ij, junk, ~ id )
HTH,
Chuck
###############################################################
### ###
### cosinor.R - functions to help in cosinor analysis ###
### ###
### Author: Charles C.Berry ###
### Date: June 12, 2002 ###
### Copyright: GPL Version 2 or higher ###
### ###
###############################################################
### following Y.L. Tong (1976) Biometrics 32:85-94
### and the notation and equation numbering therein
two.to.six <-
function( M.i, A.i, omega, t.ij, phi.i )
{
beta.i <- A.i * cos( phi.i )
gamma.i <- A.i * sin( omega * phi.i )
r.ij <- cos( omega * t.ij )
s.ij <- sin( omega * t.ij )
res <- list( M.i = M.i, beta.i=beta.i, gamma.i = gamma.i,
r.ij = r.ij, s.ij = s.ij )
attr(res,"omega") <- omega
res
}
six.to.two <-
function( M.i, beta.i, gamma.i, omega, r.ij, s.ij )
{
A.i <- sqrt( beta.i^2 + gamma.i^2 )
phi.i <- atan( gamma.i / beta.i )
t.ij <- acos( r.ij ) / omega
res <- list( M.i = M.i, A.i = A.i, phi.i = phi.i, t.ij = t.ij )
attr(res,"omega") <- omega
res
}
six.to.two.lm <-
function(fit,frame,omega)
{
cfs <- coef(fit)
M.i <- unname( cfs[ 1 ] )
beta.i <- unname( cfs["r.ij"] )
gamma.i <- unname( cfs["s.ij"] )
res <- six.to.two( M.i, beta.i, gamma.i, omega, frame$r.ij, frame$s.ij )
attr(res,"omega") <- omega
res
}
six.to.two.lsfit <-
function(fit, omega)
{
cfs <- coef(fit)
M.i <- unname( cfs[ 1, ] )
beta.i <- unname( cfs["r.ij",] )
gamma.i <- unname( cfs["s.ij",] )
A.i <- sqrt( beta.i^2 + gamma.i^2 )
phi.i <- atan( gamma.i / beta.i )
resss <- colSums( residuals( fit ) ^ 2)
totss <- colSums( fit$qr$qt[ -1, ]^2 )
regss <- totss - resss
n <- nrow( residuals( fit ) )
p <- 3
fstat <- (regss/ (p-1))/(resss/(n - p))
res <- list( M.i = M.i, A.i = A.i, phi.i = phi.i, fstat=fstat )
attr(res,"omega") <- omega
res
}
eq.two <-
function(M.i, A.i, phi.i, t.ij , omega)
{
## just expectation - no epsilon.ij used
M.i + A.i * cos( omega * t.ij - phi.i )
}
sim.two <-
function(M.mean, M.sigma, A.mean, A.sigma,
t.ij, phi.mean, phi.sigma, eps.sigma, omega, n=1 )
{
nct <-
if (length(dim(t.ij))==0) length(t.ij) else ncol(t.ij)
M.i <- rep( rnorm( n, M.mean, M.sigma ) , each = nct )
A.i <- rep( abs( rnorm( n, A.mean, A.sigma ) ), each = nct )
phi.i <- rep( rnorm( n, phi.mean, phi.sigma ) , each = nct )
epsilon.ij <- rnorm( n * nct, 0, eps.sigma )
res <- eq.two( M.i, A.i, phi.i, t.ij, omega ) + epsilon.ij
res
}
sim.cosinor <-
function(n, M.mean, M.sigma, A.mean, A.sigma,
t.ij, phi.mean, phi.sigma, eps.sigma, omega )
{
y <- sim.two( M.mean, M.sigma, A.mean, A.sigma,
t.ij, phi.mean, phi.sigma, eps.sigma, omega, n)
res <-
data.frame( y = y, t.ij = rep( t.ij, n ), id = rep( seq(n), each=length(t.ij)))
attr(res,"omega") <- omega
res
}
r.and.s <-
function(x,omega)
{
as.data.frame(two.to.six(0,0,omega, x$t.ij, 0)[c("r.ij","s.ij")])
}
cosinor.lm <-
function(formoola, frame ){
if ( !all( is.element( c("r.ij","s.ij"), colnames(frame)))) frame[, c("r.ij","s.ij")] <-
r.and.s( frame, attr(frame, "omega"))
fit <- lm(formoola, frame )
c( unlist( six.to.two.lm( fit, frame, attr( frame, "omega") )[ 1:3 ] ),
fstat=unname( summary(fit)$fstatistic["value"] ) )
}
cosinor.lm.each <-
function(formoola, frame, id ){
if ( !all( is.element( c("r.ij","s.ij"), colnames(frame)))) frame[, c("r.ij","s.ij")] <-
r.and.s( frame, attr(frame, "omega"))
id <- update(id, ~.-1 )
id.var <- all.vars( id )
id.expr <- parse(text=id.var)
clusters <- frame[,id.var]
res <- list()
for (i in unique(clusters)){
res[[i]] <- cosinor.lm( formoola, subset( frame, i == eval(id.expr) ) )
}
do.call("rbind" , res )
}
### EXAMPLE:
### X <- seq(0,24,by=2)
### junk <- sim.cosinor(20,1,.1,1,.1,X,0,.1,.1, pi/12)
###
### cosinor.lm.each( y ~ r.ij + s.ij, junk, ~ id )
###
cosinor.lsfit.each <-
function( ymat, rs.mat, omega )
{
fit <- lsfit( rs.mat, ymat )
do.call( "cbind", six.to.two.lsfit(fit, omega ))
}
########################
### end of cosinor.R ###
########################
> Thanks
> Anne Berger
> Institute of Zoo- and Wildlife Research, Berlin, Germany
> Swedish University of Agricultural Sciences, Dept. of Wildlife, Fish and Environmental Studies, Umeå, Sweden
>
> [[alternative HTML version deleted]]
>
>
Charles C. Berry (858) 534-2098
Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
More information about the R-help
mailing list