[R] Conservative "ANOVA tables" in lmer

Fabian Scheipl f.abian at gmx.net
Fri Sep 8 14:25:05 CEST 2006


Dear list,

I have written functions to perform simulation-based tests of 
HO: Var(Random Effects)=0, beta_foo=0 
in linear mixed models based on the exact distribution of the LR- and
Restricted LR-test statistic (see research presented in
Crainiceanu, Ruppert: "LRT in LMM with 1 variance component",
J. R. Statist. Soc. B (2004), 66, Part 1, pp. 165-185 )
-they are about 15-20 times faster than the parametric bootstrap.

At the moment, the exact distributions are only easily simulated for the
case of 1 single variance component/random effect and i.i.d. errors; feasible approximations for the "multivariate" case are currently
being investigated and will be implemented soon.

the syntax looks something like this:

#begin code:

data(sleepstudy)
summary(sleepstudy)  #Effect of sleep deprivation on reaction time
xyplot(Reaction~Days|Subject, data=sleepstudy)
m<-lmer(Reaction~Days+(Days-1|Subject),data=sleepstudy) 
#random slopes, but no random intercept  
#doesna make sense, but it's just an example
summary(m)

#test for individual heterogeneity based on RLRT
#(No restrictions on fixed effects under H0)
#HO: lambda=Var(RandomSlopes)/Var(error)==0 <==> Var(RandomSlopes)==0

t3<-RLRT1SimTest(m, lambda0=0, seed=5, nsim=10000)

#will produce output:
#HO: lambda = 0 ; p-value = 0 
# observed lambda = 0.06259639 

#test for influence of Days based on LRT 
#(restriction on fixed efects: beta_Days==0)

m0<-lm(Reaction~1,data=sleepstudy)
t4<-LRT1SimTest(m, m0, seed=10, nsim=10000)

#will produce output:
#Model under HO:  Reaction ~ 1 ;
#Model under HA:  Reaction ~ Days + (Days - 1 | Subject) ;
# p-value = 0 
# observed lambda = 0.06259639 

#end code 
                      
If you are interested in using these functions i'll be glad
to send them to you-
be aware, however, that you can only use them for testing
"1 Random Effect" vs. "no Random Effect" in a model with i.i.d. errors!!

The plan is to put them in a package beginning next year and
use them as a basis for an (exact) anova.lmer() method.

Greetings,
Fabian Scheipl


--



More information about the R-help mailing list