[R-sig-ME] account for temporal correlation

[Ben Bolker]

matteo dossena <matteo.dossena at ...> writes:

> 
> Dear list,
> 
> I have a data frame like this:
> 
set.seed(5)
mydata <- data.frame(var = rnorm(100,20,1),
                      temp = sin(sort(rep(c(1:10),10))),
                      subj = as.factor(rep(c(1:10),5)),
                      time = sort(rep(c(1:10),10)),
                      trt = rep(c("A","B"), 50))
> 
> I need to model the response of var as a function of temp*trt
> and to do so I'm using the following model:
>
library(nlme)
model <- lme(var~temp*trt,random=~1|subj,mydata)
> 
> however, since i have repeated measurement of the same subject,
> i.e. I measured var in subj1 at time1,2,3..
> I must consider the non independence of the residuals.
> moreover, temp is also a function of time, but i'm not sure how
> to include this in my model.
> 
> I'm following the approach in Zuur et al 2009, so I checked for
> temporal auto-correlation using the
> function afc()
> In fact the residuals follow the temporal patter of temperature with time.
> However, here I'm not interested in the temporal dependence of temperature
> and consequently the effect of
> this on var.
> Instead what i need to do is to account for the
> correlation between consecutive measures (made on the same
> subject) in the error term of the model. 
> 
> here is my attempt to do so:
> 

model1 <- lme(var~temp*trt,random=~1|subj,
                correlation=corAR1(form=~1|subj),mydata)

model1$modelStruct$corStruct

Correlation structure of class corAR1 representing
        Phi 
-0.05565362 

  You shouldn't be nesting time within subject.  'subject' is all the grouping
you need here.

intervals(model1)

 gives (approximate!!) CI for the correlation structure parameter
(-0.27,0.77) in this case

  Of course, in this case we don't expect to find anything interesting
because these are simulated data without any correlation built in.

  These plots are useful.

plot(ACF(model),alpha=0.05)
plot(ACF(model1),alpha=0.05)  ## should be ALMOST identical to the one above
 ## taking correlation into account:
plot(ACF(model1,resType="normalized"),alpha=0.05)