[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