[R-sig-eco] Periodic regression

[Dixon, Philip M [STAT]]

Thiago,

I would not try to interpret each periodic component separately.  That linear regression model is just a convenient way to fit a model with a phase-shifted cosine curve:

log (mean index) = b0 + A*cos(c*time - theta)

The amplitude of the oscillation is A; the phase shift is theta.  Those can be estimated from the regression coefficients in the linear model:
log (mean index) = b0 + b1 cos(c*t) + b2 sin(c*t)

as:
A = sqrt(b1^2 + b2^2)
theta = atan(b2/b1)

Here's R code to implement these calculations:
tad.lm<-glm(index~(cos(2*pi*(hora/24)))+(sin(2*pi*(hora/24))), 
  family=quasipoisson)
tad.coef <- coef(tad.lm)[2:3]

A <- sqrt(sum(tad.coef^2))   # where tad.coef is only periodic components
phase <- atan2(tad.coef[2],tad.coef[1])
#  in radians, in the appropriate quadrant

maxhora <- phase*24/(2*pi)
maxhora <- ifelse(maxhora < 0, 24+maxhora, maxhora)
#  and in hours, shifting by 24 hr if negative

Standard errors are obtained by the delta method.
Confidence intervals are best obtained by bootstrapping.  You should think hard about what sort of bootstrap to use (observations, standardized residuals, parametric on the estimates, parametric on the observations).

Chester Bliss's pamphlet on Periodic regression in Biology and ?? is a good introduction, but remember that his computing resources is 60 years old.

The test of "no aggregation" that I suggested (one of many possibilities) = test of "no periodic components".  That is the test that both b1 = 0 and b2 = 0 in the linear model above.  Can use either a drop in quasi-deviance or a C Beta test on the estimates.  Asymptotically equivalent; usually similar but rarely identical in practical samples.  

R code is:

tad.lm<-glm(index~(cos(2*pi*(hora/24)))+(sin(2*pi*(hora/24))), 
  family=quasipoisson)
tad.lm0 <- glm(index~1, family=quasipoisson)

# F (drop in quasi deviance) test for both periodic coefficients = 0
anova(tad.lm0, tad.lm, test='F')

# C Beta test for both periodic coefficients = 0

tad.coef <- coef(tad.lm)[2:3]
tad.vc <- vcov(tad.lm)[2:3,2:3]
# extract coefficients and variance-covariance matrix for periodic terms
#  i.e. dropping the estimated intercept

npar <- length(tad.coef)
nobs <- length(hora)
# extract # parameters in test and # observations in data set

tad.f <- (t(tad.coef) %*% solve(tad.vc) %*% tad.coef) / npar
c(F = tad.f, probF=pf(tad.f, npar, nobs-(npar+1), lower=F))
 
Unless there are issues of general interest to the list, I suggest any further discussion be done off-list.

Best wishes,
Philip Dixon