[R-sig-ME] marginality principle for mixed effects models
Ben Bolker
bolker at ufl.edu
Mon Apr 6 20:57:10 CEST 2009
I think you could also say
data(Oats,package="nlme")
without loading the whole nlme package.
Douglas Bates wrote:
> A minor technical point. It is not a good idea to have the lme4 and
> nlme packages attached at the same time. They both define generic
> functions fixef, ranef, etc. and these definitions conflict.
>
> There is a package called MEMSS that provides the data sets from nlme
> without conflicts.
>
> On Fri, Apr 3, 2009 at 9:17 PM, Steven McKinney <smckinney at bccrc.ca> wrote:
>>
>> Dear list,
>>
>> I'm working on mapping all my old-school fixed effects linear
>> modeling teachings into the fabulous new(er) mixed effects world.
>>
>> I've been looking at the "Oats" data used in the nlme library.
>> Here's a set of commands that build successively richer models,
>> and tests to evaluate differences due to variety.
>>
>> library("nlme")
>> library("lme4")
>> plot(Oats, inner = TRUE)
>> print(om1 <- lmer(yield ~ nitro + (1 | Block), data = Oats))
>> print(om2 <- lmer(yield ~ nitro + Variety + (1 | Block), data = Oats))
>> print(om3 <- lmer(yield ~ nitro + (1 | Block/Variety), data = Oats))
>> print(om4 <- lmer(yield ~ nitro + Variety + (1 | Block/Variety), data = Oats))
>>
>> anova(om4, om3) ## Does this violate the principle of marginality? (Testing main effect in presence of interaction?)
>> anova(om4, om2) ## Would this be the first test of variety?
>> anova(om4, om1) ## Would this not be the omnibus test for variety?
>> anova(om2, om1)
>>
>> If I wanted to assess the effect of variety, it seems to me
>> the first test I'd want to look at is model om4 versus model om2.
>> This it seems to me would be in the spirit of testing an interaction
>> term with both main effects in the model in the old fixed effects
>> world.
>>
>> An omnibus test of the importance of variety would seem to me
>> to be model om4 versus om1.
>>
>> Testing om4 versus om3 seems to me to violate the marginality
>> principle (testing a main effect in the presence of an interaction
>> involving that main effect). Or is there something different in
>> the mixed effects world - does this marginality principle hold
>> for this scenario? The plots and all the other tests seem to
>> strongly suggest that there are important differences due to variety,
>> but this test suggests otherwise. This test does not seem
>> appropriate.
>>
>> Any comments? Is this paradigm mapping reasonable?
>>
>>
>>
>>> library("nlme")
>>> library("lme4")
>>> plot(Oats, inner = TRUE) # Lattice plot of data
>>> print(om1 <- lmer(yield ~ nitro + (1 | Block), data = Oats))
>> Linear mixed-effects model fit by REML
>> Formula: yield ~ nitro + (1 | Block)
>> Data: Oats
>> AIC BIC logLik MLdeviance REMLdeviance
>> 610.7 617.5 -302.4 616.4 604.7
>> Random effects:
>> Groups Name Variance Std.Dev.
>> Block (Intercept) 243.34 15.599
>> Residual 254.99 15.968
>> number of obs: 72, groups: Block, 6
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) 81.872 7.104 11.524
>> nitro 73.667 8.416 8.753
>>
>> Correlation of Fixed Effects:
>> (Intr)
>> nitro -0.355
>>
>>
>>> print(om2 <- lmer(yield ~ nitro + Variety + (1 | Block), data = Oats))
>> Linear mixed-effects model fit by REML
>> Formula: yield ~ nitro + Variety + (1 | Block)
>> Data: Oats
>> AIC BIC logLik MLdeviance REMLdeviance
>> 601.6 613 -295.8 608.8 591.6
>> Random effects:
>> Groups Name Variance Std.Dev.
>> Block (Intercept) 245.03 15.653
>> Residual 234.73 15.321
>> number of obs: 72, groups: Block, 6
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) 81.872 7.069 11.582
>> nitro 73.667 8.075 9.123
>> Variety1 2.646 2.211 1.196
>> Variety2 -3.174 1.277 -2.486
>>
>> Correlation of Fixed Effects:
>> (Intr) nitro Varty1
>> nitro -0.343
>> Variety1 0.000 0.000
>> Variety2 0.000 0.000 0.000
>>
>>
>>> print(om3 <- lmer(yield ~ nitro + (1 | Block/Variety), data = Oats))
>> Linear mixed-effects model fit by REML
>> Formula: yield ~ nitro + (1 | Block/Variety)
>> Data: Oats
>> AIC BIC logLik MLdeviance REMLdeviance
>> 601 610.1 -296.5 604.3 593
>> Random effects:
>> Groups Name Variance Std.Dev.
>> Variety:Block (Intercept) 121.33 11.015
>> Block (Intercept) 210.28 14.501
>> Residual 165.51 12.865
>> number of obs: 72, groups: Variety:Block, 18; Block, 6
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) 81.872 6.944 11.79
>> nitro 73.667 6.781 10.86
>>
>> Correlation of Fixed Effects:
>> (Intr)
>> nitro -0.293
>>
>>
>>> print(om4 <- lmer(yield ~ nitro + Variety + (1 | Block/Variety), data = Oats))
>> Linear mixed-effects model fit by REML
>> Formula: yield ~ nitro + Variety + (1 | Block/Variety)
>> Data: Oats
>> AIC BIC logLik MLdeviance REMLdeviance
>> 594.5 608.1 -291.2 601.3 582.5
>> Random effects:
>> Groups Name Variance Std.Dev.
>> Variety:Block (Intercept) 108.94 10.437
>> Block (Intercept) 214.57 14.648
>> Residual 165.55 12.867
>> number of obs: 72, groups: Variety:Block, 18; Block, 6
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) 81.872 6.946 11.786
>> nitro 73.667 6.781 10.863
>> Variety1 2.646 3.539 0.748
>> Variety2 -3.174 2.043 -1.553
>>
>> Correlation of Fixed Effects:
>> (Intr) nitro Varty1
>> nitro -0.293
>> Variety1 0.000 0.000
>> Variety2 0.000 0.000 0.000
>>
>>
>>
>>
>>> anova(om4, om3) ## Does this violate the principle of marginality? (Testing main effect in presence of interaction?)
>> Data: Oats
>> Models:
>> om3: yield ~ nitro + (1 | Block/Variety)
>> om4: yield ~ nitro + Variety + (1 | Block/Variety)
>> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
>> om3 4 612.30 621.41 -302.15
>> om4 6 613.28 626.94 -300.64 3.0197 2 0.2209
>>
>>
>>> anova(om4, om2) ## Would this be the first test for effect of variety?
>> Data: Oats
>> Models:
>> om2: yield ~ nitro + Variety + (1 | Block)
>> om4: yield ~ nitro + Variety + (1 | Block/Variety)
>> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
>> om2 5 618.85 630.23 -304.42
>> om4 6 613.28 626.94 -300.64 7.5629 1 0.005958 **
>> ---
>> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>
>>
>>> anova(om4, om1) ## Would this not be the omnibus test for variety?
>> Data: Oats
>> Models:
>> om1: yield ~ nitro + (1 | Block)
>> om4: yield ~ nitro + Variety + (1 | Block/Variety)
>> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
>> om1 3 622.40 629.23 -308.20
>> om4 6 613.28 626.94 -300.64 15.114 3 0.001722 **
>> ---
>> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>
>>
>>> anova(om2, om1)
>> Data: Oats
>> Models:
>> om1: yield ~ nitro + (1 | Block)
>> om2: yield ~ nitro + Variety + (1 | Block)
>> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
>> om1 3 622.40 629.23 -308.20
>> om2 5 618.85 630.23 -304.42 7.5512 2 0.02292 *
>> ---
>> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>
>>
>> Steven McKinney, Ph.D.
>>
>> Statistician
>> Molecular Oncology and Breast Cancer Program
>> British Columbia Cancer Research Centre
>>
>> email: smckinney +at+ bccrc +dot+ ca
>>
>> tel: 604-675-8000 x7561
>>
>> BCCRC
>> Molecular Oncology
>> 675 West 10th Ave, Floor 4
>> Vancouver B.C.
>> V5Z 1L3
>> Canada
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc
More information about the R-sig-mixed-models
mailing list