The R Book

Fixed Effects in Hierarchical Sampling

Given that the gender effect in our hierarchical sampling example on p. 639 was so large, and that gender makes a sensible fixed effect (it has informative factor levels: male and female), we might fit gender as a main effect. The important point to note is that when you want to compare models with different fixed effects using lme you must change the fitting method from the default REML to the optional maximum likelihood method="ML". This then allows you to use anova to compare lme models with different fixed effects:

model10<-lme(subject~gender,random=~1|town/district/street/family/gender,
                                          method="ML")
model11<-lme(subject~1,random=~1|town/district/street/family/gender,
                                     method="ML")
anova(model10,model11)

         Model   df       AIC        BIC      logLik     Test    L.Ratio  p-value
model10      1    8  3331.584   3368.218   -1657.792
model11      2    7  3352.221   3384.276   -1669.111   1 vs 2   22.63755  <.0001

It is clear that the model with gender as a fixed effect (model10) is vastly superior to the model with out any fixed effects (p < 0.0001). It has a much lower AIC, despite its extra parameter. The variance components have been little affected by fitting gender as a fixed effect, and the effect size of gender is given by:

summary(model10)
Fixed effects: subject ~ gender
                  Value   Std.Error    DF     t-value   p-value
(Intercept)    8.976328   0.6332402   360   14.175234         0
gendermale    -1.930773   0.3936610   179   -4.904659         0

You can see what the parameter values are by looking at the treatment means:

tapply(subject,gender,mean) ...