The effect of chick

Until now, we ignored differences due to chicken; but these do exist and are substantial:

ggplot(data = ChickWeight,aes(x = Time,y = weight,colour = Chick,group = Chick)) + 
       geom_line() + facet_grid(~Diet) + theme(legend.position = 'none')

Including random effects in the model

Models with fixed effects and random effects other than error are called mixed models. We use the package lme4 for fitting mixed models.

In matrix notation a mixed model is

\(\mathbf{y} = \mathbf{Xb} + \mathbf{Zu} + \mathbf{e}\)

with \(\mathbf{u} \sim \mathbf{N}(0,\sigma^2_u)\) and \(\mathbf{e} \sim \mathbf{N}(0,\sigma^2_e)\). In this model, vector \(\mathbf{u}\) may be a vector of chickeffects.

Fitting the model in lme4

lme41 <- lmer(weight~Diet + Time + (1|Chick),ChickWeight)
lme41

## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ Diet + Time + (1 | Chick)
##    Data: ChickWeight
## REML criterion at convergence: 5584.004
## Random effects:
##  Groups   Name        Std.Dev.
##  Chick    (Intercept) 22.92   
##  Residual             28.27   
## Number of obs: 578, groups:  Chick, 50
## Fixed Effects:
## (Intercept)        Diet2        Diet3        Diet4         Time  
##      11.244       16.210       36.543       30.013        8.717

Checking differences

We can wonder if this model is an improvement over the ANCOVA model ancW2. We do this indirectly because the function can not handle an object of class and one of class lm.

ancW2 <- lm(weight~Diet + Time, ChickWeight)
AIC(logLik(ancW2))

## [1] 5789.607

AIC(logLik(lme41))

## [1] 5598.004

We cautiously conclude (this is no a hypothesis test) that including a random intercept for Chick improves the fit of our model.

Including a random regression coefficient

We could also think that not only the intercept but also the slope is dependent on the chick:

lme42 <- lmer(weight~Diet + Time + (Time|Chick), ChickWeight)
lme42

## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ Diet + Time + (Time | Chick)
##    Data: ChickWeight
## REML criterion at convergence: 4803.754
## Random effects:
##  Groups   Name        Std.Dev. Corr 
##  Chick    (Intercept) 12.40         
##           Time         3.76    -0.98
##  Residual             12.78         
## Number of obs: 578, groups:  Chick, 50
## Fixed Effects:
## (Intercept)        Diet2        Diet3        Diet4         Time  
##      26.356        2.839        2.004        9.255        8.444

Testing

Test whether model with random slopes and intercepts fits better than model with random intercepts only

anova(lme41, lme42)

## refitting model(s) with ML (instead of REML)

## Data: ChickWeight
## Models:
## lme41: weight ~ Diet + Time + (1 | Chick)
## lme42: weight ~ Diet + Time + (Time | Chick)
##       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## lme41  7 5619.2 5649.7 -2802.6   5605.2                             
## lme42  9 4834.1 4873.3 -2408.0   4816.1 789.12      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Return to main page

• Exercises mixed models
  • Solutions exercises mixed models
• Return to main page