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)