[R-sig-ME] marginality principle for mixed effects models
Douglas Bates
bates at stat.wisc.edu
Mon Apr 6 20:51:47 CEST 2009
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
>
More information about the R-sig-mixed-models
mailing list