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 [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-ecology |
Free forum by Nabble | Edit this page |