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

_______________________________________________
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

Reply via email to