[R-sig-ME] marginality principle for mixed effects models
Douglas Bates
bates at stat.wisc.edu
Mon Apr 6 22:36:21 CEST 2009
On Mon, Apr 6, 2009 at 1:57 PM, Ben Bolker <bolker at ufl.edu> wrote:
>
> I think you could also say
>
> data(Oats,package="nlme")
>
> without loading the whole nlme package.
Indeed. However, the versions of the data sets in MEMSS are slightly
different from the nlme versions. The versions in nlme with all
inherit from the groupedData class. The versions in MEMSS are data
frames. I now feel that the groupedData class was a bit too baroque
and that it is better to stick with something simpler like a data
frame. Every once in a while some of the characteristics of the
groupedData class, like converting grouping factors to ordered
factors, causes unexpected results. That is why the groupedData class
was not propagated to the lme4 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