[R-sig-ME] LRT tests with lme4 and single random effects
Ned Dochtermann
ned.dochtermann at gmail.com
Thu Oct 20 00:59:45 CEST 2011
Hi All,
As part of looking at some issues with repeatability from a repeated
measures mixed-model analysis I've been asked by a colleague about
statistical power in regards to the "significance" of repeatability
estimates. In other packages (SAS and ASReml) this would be determined based
on the significance of the random effect. The problem I ran into of course
is that there isn't a ready-made way to test the significance of a random
factor in a single random effect model since there isn't a function that has
commensurate likelihoods with lme4. I know I could do this with lme and gls
but as I had been using lmer to look at the accuracy of repeatability
estimates, I'd prefer to stick with the lme4 library. I know mcmcsamp could
provide CI's around the random term's variance but that's not quite what
we're looking for.
The work around I came up with was to add in an uninformative (except to a
small degree due to chance) variable and include it as a random factor in a
separate model. The deviance for this model could then be used as a stand-in
for a model omitting the random factor.
My code to do this was:
icc<-0.3; n<-100; reps<-2
ind.0<-rnorm(n,0,sqrt(icc))
ind.e<-rnorm(n*reps,0,sqrt(1-icc))
ind<-gl(n,1)
ind.data<-data.frame(cbind(sample(c(1:n),n*reps,T),
rep(seq(n),reps),
matrix(ind.0[ind]+ind.e,ncol=1)))
col.names<-c("rnd","ind","y"); names(ind.data)<-col.names
mm1<-lmer(y~1+(1|ind),data=ind.data)
mmN<-lmer(y~1+(1|rnd),data=ind.data)
mm1.dev<-summary(mm1)@deviance[1]
mmN.dev<-summary(mmN)@deviance[1]
(lrt.p<-1-pchisq(abs(mm1.dev-mmN.dev),df=1))
If lrt.p is less than alpha this then suggests a "significant" random
effect. Due to the issues with p-values I suppose lrt.p should also be
divided by two. It seems like this should work, although since 1|rnd will by
chance reduce the deviance some small amount it might be conservative. Is
this a completely idiotic approach due to me missing something obvious to
statisticians that was not obvious to me as a user?
Sorry if this has come up, I did search but didn't find anything that quite
fit what I was looking for.
Ned
--
Ned Dochtermann
Department of Biology
University of Nevada, Reno
ned.dochtermann at gmail.com
http://wolfweb.unr.edu/homepage/mpeacock/Ned.Dochtermann/
http://www.researcherid.com/rid/A-7146-2010
--
More information about the R-sig-mixed-models
mailing list