Hi
Is there a way to investigate auto correlation in a model fit with lmer, similar to the ACF.lme function for lme models? I haven't been able to find any thread discussion on the topic and some advice would be much appreciated. I thought that maybe I could use the by-subject outputs from acf.fnc  (languageR package), but this doesn't work (see below).

library(nlme)
library(lme4)
library(languageR)

##consider this data set
nID<-20
nwithin<-10
prosw<-nID*nwithin

growth<-vector(length=prosw)
age<-seq(2,11, length=nwithin)
a<-rnorm(prosw,mean=2)
b<-120;c<-10;d<-5;e<- -2;rho<--0.9
for (i in 1:nID)
{
sf<-0.1*rnorm(1);sc<-2*rnorm(1);ss<-rnorm(1,mean=2);err<-0.1*rnorm(1)
for (j in 1:nwithin)
{
growth[nwithin*(i-1)+j]<-sf +  b + (c+sc)*a[nwithin*(i-1)+j] - (d+ss)*age[j] +sf*a[nwithin*(i-1)+j]+
                   err
err<-err*rho + 0.1*rnorm(1)
}
}
a<-seq((1980+sample(0:10,1)),length.out=10,by=1)
b<-seq((1980+sample(0:10,1)),length.out=10,by=1)
c<-seq((1980+sample(0:10,1)),length.out=10,by=1)
d<-seq((1980+sample(0:10,1)),length.out=10,by=1)
e<-seq((1980+sample(0:10,1)),length.out=10,by=1)
f<-seq((1980+sample(0:10,1)),length.out=10,by=1)
g<-seq((1980+sample(0:10,1)),length.out=10,by=1)
h<-seq((1980+sample(0:10,1)),length.out=10,by=1)
i<-seq((1980+sample(0:10,1)),length.out=10,by=1)
j<-seq((1980+sample(0:10,1)),length.out=10,by=1)
k<-seq((1980+sample(0:10,1)),length.out=10,by=1)
l<-seq((1980+sample(0:10,1)),length.out=10,by=1)
m<-seq((1980+sample(0:10,1)),length.out=10,by=1)
n<-seq((1980+sample(0:10,1)),length.out=10,by=1)
o<-seq((1980+sample(0:10,1)),length.out=10,by=1)
p<-seq((1980+sample(0:10,1)),length.out=10,by=1)
q<-seq((1980+sample(0:10,1)),length.out=10,by=1)
r<-seq((1980+sample(0:10,1)),length.out=10,by=1)
s<-seq((1980+sample(0:10,1)),length.out=10,by=1)
t<-seq((1980+sample(0:10,1)),length.out=10,by=1)
year<-c(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t)##sorry couldn't get for loop to work for this!
age<-rep(age, nID)
ID<-rep(seq(1,nID), each=nwithin)
growthdata<-data.frame(ID,age,year,growth)
growthdata$ID<-as.factor(growthdata$ID)

##I would like to fit a mixed model with crossed random effects of ID and year (not every level of ID is present in each year). It is highly likely that there is correlation between sequential measurements within ID. I can account for this using an auto-regressive correlation structure of the form correlation=corAR1(form=~age|ID) in lme, but due to the high number of ID and year levels in my full dataset, crossing ID with year is not possible. Therefore, I thought I'd try lmer and fit the model:

model.lmer<-lmer(growth~age+(age|ID)+(1|as.factor(year)),growthdata)

###  and see how if performed. I fit a similar model in lme to compare the outputs of ACF and acf.fnc using model residuals. I naively thought that averaging across levels of 'lag' in the acf.fnc output may work, but it obviously doesn't! My real data has a much higher level of autocorrelation.

model.lme<-lme(growth~age,random=~age|ID,growthdata)
ACF(model.lme)
lag         ACF
   0  1.00000000
   1 -0.16324159
   2 -0.05552011
   3 -0.16254222
   4 -0.10331816
   5 -0.06783792
   6 -0.03268318
   7  0.08196762
   8 -0.34495610
   9  0.17929713
growthdata$resid<-resid(model.lme)
y<-acf.fnc(growthdata,group='ID',time='age',x='resid',plot=F)
yy<-as.matrix(tapply(y$Acf,y$Lag,mean))
         [,1]
0  1.00000000
1 -0.12350304
2 -0.10730408
3 -0.04293020
4 -0.09515590
5 -0.07507144
6  0.01232782
7  0.01487862
8 -0.07217361
9 -0.01106817

Cheers

John







	[[alternative HTML version deleted]]


