[R] Conservative "ANOVA tables" in lmer
Fabian Scheipl
f.abian at gmx.net
Fri Sep 8 16:24:30 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