Chapter 14: Multilevel Modelling

Welcome to the online content for Chapter 14!

As always, I’ll assume that you’ve already read up to this chapter of the book and worked through the online content for the previous chapters. If not, please do that first.

As always, click the ‘Run Code’ buttons below to execute the R code. Remember to wait until they say ‘Run Code’ before you press them. And be careful to run these boxes in order if later boxes depend on you having done other things previously.

Ignoring clustering

We began the chapter by ignoring clustering, so let’s reproduce those analyses here.

First, let’s read in the data set:

We have 20 people, with their locations, ages and heights.

Let’s plot height against age:

We see the same pattern of points that I showed in the chapter.

Now, let’s ignore the clustering into West and East and just do an ordinary regression of height on age. We use the lm function to do this, just as we did in Chapter 8.

I’ve called the model ‘model1’, because we’re going to produce a sequence of different models, which we’ll want to compare, so it’s helpful to number them.

We see a positive value of 0.8206 for the slope of the line, which corresponds to how it looks in the plot underneath. This is the effect of age on location. The values for adjusted \(R^2\) and the \(p\) value of .007 match those in the chapter.

We can do separate regressions for the two locations by using the subset function:

The values for adjusted \(R^2\) and the \(p\) values match those given in the chapter for each location.

Null single-level model

To run multilevel models in R, we’ll use a package called lme4. So, let’s install that:

To run the null single-level model, we don’t actually need lme4. We do it just with lm:

There’s no slope here (it’s forced to be zero), and the estimate of the intercept is just the mean height of 164.9765 cm.

Null 2-level model

From now on, we need to use the lmer function in lme4:

The vertical line | in the formula tells R that location is the random effect.

The REML=FALSE argument is a technicality to do with the method used for estimating the parameters. I don’t discuss this in the book.

We can see AIC and BIC values, which are information criteria (described in Chapter 10). They’re similar to deviance, and, for all of these, smaller is better. AIC and BIC are larger than the deviance because they penalise models for the parameters that they include.

Under ‘Random effects’, we get an estimate of the between-location (Level 2) variance in height (9.450) and the within-location between-person (Level 1) residual variance (2.622). These match the values stated in the chapter.

If we express the first of these as a fraction of the total, we get

This tells us that 78% of all of the variance in height can be ascribed to differences in location.

To see if Model 5 is significantly better than Model 4 (the null single-level model), we can do a chi-squared test by using the anova function, which compares the fit of two models:

Both the AIC and the BIC decrease, and we can see that the deviance has dropped from 106.576 to 83.259, which is good. (Less deviance means that the model fits better.)

The difference is the 23.317, which is a chi-squared statistic on 1 degree of freedom, giving a \(p\) value of 1.374e-06, which means \(1.374\times10^{-6}\) . For such a small number (a probability of about 1 in a million), we write \(p<.001\). This means that there’s a significant improvement from Model 4 to Model 5.

Random-intercepts model

Here, the code is:

We get a negative value for the slope for age (-0.9083), because increasing age predicts lower heights.

We can compare this with Model 5:

As before, both AIC and BIC decrease, which is good. The deviance drops by 15.354, with one additional parameter, and the chi-squared test (on 1 degree of freedom again) gives a \(p\) value of 8.911e-05, which is also \(<.001\). So, again, we have a significant improvement in model fit, this time from Model 5 to Model 6.

Random-slopes model

Finally, we’ll try allowing the slopes to vary:

We get a warning message “boundary (singular) fit: see help(‘isSingular’)”, which we’ll ignore. This just means that the model fits, but the random effects are very small, as we’ll see.

The slope for age is -0.9114, which is very similar to the slope for age that we had in Model 6.

Let’s compare Model 7 with Model 6:

This time, both AIC and BIC increase, which is not good. And the deviance drops by only 0.0066. This time, we have lost 2 more degrees of freedom, so we have to test 0.0066 against a chi-squared distribution with 2 degrees of freedom. The \(p\) value is 0.9967, so there’s no significant improvement to the model.

We’d probably choose to stick with Model 6 as out best model, and ignore Model 7.