[R-sig-ME] marginality principle for mixed effects models
Steven McKinney
smckinney at bccrc.ca
Sat Apr 4 04:17:47 CEST 2009
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
More information about the R-sig-mixed-models
mailing list