An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20090108/a1be605c/attachment-0001.pl>
cosinor analysis
2 messages · Anne Berger, Charles C. Berry
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