A gentle introduction to mixed-effects models

In a previous post, I discussed how to use linear regression to analyse whether one or more variables influence an outcome. Linear regression is great, but it is not always the best solution; one of the extremely important assumptions of linear regression is that the observations are independent. However, there are a lot of situations when that does not apply. This post will explore how to use mixed-effects models in R to deal with this!

Experimental designs giving non-independent measures

Not interested in the theory? Jump directly to the example!

Two of the main types of experimental design where you end up with non-independent observations are repeated measures and nested designs.
In repeated measures design, you are taking measurements from the same experimental unit multiple times; a typical example is taking multiple blood samples from the same patient over the course of several visits. Because the measures always come from the same individual, they are not independent, and they are likely to be more correlated between each other than they are with measures from a different individual.

In nested design, you are taking measurements that have some hierarchical or nested structure. For example, you may take biopsies from various patients at different clinics, and then analyse multiple cells in each biopsy. The cells coming from the same biopsy come from the same piece of tissue, thus they are not independent and may have some higher degree of correlation than cells from different biopsies. Furthermore, patients at the same clinic may also show some higher degree of correlation because they are treated by the same medical staff, using the same procedures and maybe because they live in the same area. So we have a nesting of the type clinic > patient/biopsy > cells.

Why should I bother (replication and pseudoreplication)?

So… what’s the big deal with having non-independent samples? We will see some examples below, but the issue is that we risk creating models that do not correctly represent the real variability in our sample. For instance, think about the situation on the left in the figure above. If we model those as 12 independent samples, we are assuming that there are 11 degrees of freedom in our system, which is not true, as each measurement at time 2 will be somehow related to the respective measurement at time 1.

Replication is an essential part of experimental design. We perform a measurement multiple times on multiple experimental units to increase how precisely the estimations we perform on the samples approximate the parameters of the population.

By incorrectly increasing the degrees of freedom in our model, we generate inaccurate statistical models that can misrepresent reality. The process of (maliciously or not) inflating the degrees of freedom is called pseudoreplication, and it is a well-known issue in biological and other sciences (see for example Lazic 2010 or Lazic et al. 2018). By reporting having measured 12 independent samples instead of reporting 4 independent groups of 3 non-independent samples, and analysing your data according, you are falling into the pseudoreplication trap… and bad things can happen there!

A simple example

Let’s say we measure the blood level of substance X, over the course of 3 time points in 10 people. We want to know whether time has an effect on the levels of X. I created some example data to use for this post.

Let’s have a quick look at the data

ggplot(data, aes(x=Time, y=X)) +
geom_line(aes(col=Patient))

This plot tells us that time has no clear effect on the level of X, and that there is quite a bit of variability amongst patients in their X levels. This is a good example of repeated measures and one that will allow us to see the issues with pseudoreplication. So, let’s start to analyse these data using R!

The bad way

We start by creating a simple linear model to analyse our data such as:

$$X = \beta_0 + \beta_1 * Time + \epsilon$$

model <- lm(X ~ Time, data)
summary(model)

Call:
lm(formula = X ~ Time, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.2700  -2.5912   0.9825   3.7788  10.7300 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   29.980      2.879  10.413 3.91e-11 ***
Time          -0.205      1.333  -0.154    0.879    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.96 on 28 degrees of freedom
Multiple R-squared:  0.0008443,	Adjusted R-squared:  -0.03484 
F-statistic: 0.02366 on 1 and 28 DF,  p-value: 0.8789

The summary of the model tells us that Time does not influence the levels of X however, you can see that the model is representing a system with 28 degrees of freedom, which we know is wrong! Indeed, what is happening is that we are considering the measurements as all independent.

We can see this if we plot the residuals by patient

In the definition of our linear model we assume $\epsilon \sim N(0, \sigma_{\epsilon}^{2})$, which, albeit being true here when looking at all the residuals (the orange bar in the figure), it is clearly not true when looking at residuals from the single individuals. What is happening is that we are artificially inflating the sample variability, by putting the “patient effect” inside the residuals.

Let’s see if there is a better way!

The useless way

Seen that the problem lies in the individual patients we may be tempted to use it as a factor.

model2 <- lm(X ~ Time + Patient, data)

Indeed, this solves the problem, and now the patient effect is not in the residuals anymore, however the model summary raises a different issue

Call:
lm(formula = X ~ Time + Patient, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4383 -1.0083  0.1025  0.9967  3.3617 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  36.9433     1.6598  22.258 4.52e-15 ***
Time         -0.2050     0.5082  -0.403 0.691176    
PatientB     -5.9333     1.8557  -3.197 0.004743 ** 
PatientC      1.4000     1.8557   0.754 0.459845    
PatientD    -16.1000     1.8557  -8.676 4.92e-08 ***
PatientE     -8.4000     1.8557  -4.527 0.000231 ***
PatientF     -6.5333     1.8557  -3.521 0.002286 ** 
PatientG     -8.8333     1.8557  -4.760 0.000136 ***
PatientH    -15.8667     1.8557  -8.550 6.16e-08 ***
PatientI     -4.2000     1.8557  -2.263 0.035521 *  
PatientJ     -5.1667     1.8557  -2.784 0.011824 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.273 on 19 degrees of freedom
Multiple R-squared:  0.9014,	Adjusted R-squared:  0.8495 
F-statistic: 17.37 on 10 and 19 DF,  p-value: 1.615e-07

Indeed, the residual standard error is about 1/3 of what it was before, however, what is happening here is that we are creating a lot of dummy variables to represent the patients. So our model is something like

$$X = \beta_0 + \beta_1 * Time + \beta_2 * PB + \beta_3 * PC + \ldots + \epsilon$$

Where

$$PB = \begin{cases}& 1\text{ if Patient} = B\\ & 0\text{ otherwise}\end{cases}$$

$$PC = \begin{cases}& 1\text{ if Patient} = C\\ & 0\text{ otherwise}\end{cases}$$

And so on…

So this model is not very useful, as it only models those 10 patients and cannot extrapolate to a patient 11 unless you estimate a new variable for them, which will need measuring X for that patient. Essentially, we have created a model that tells us what we already know and cannot do anything else.

The proper way, using mixed-effects models!

So, we finally get to do things the proper way, using mixed-effects models!

The idea is that we don’t use the patient as a predictor or, in mixed-effect lingo as a fixed effect, but rather like a random effect, a random variable, independent from the residuals, which influences our outcome but that we are not interested in determining. In other terms

$$X_{kij} = \beta_0 + \beta_1 * Time + b + epsilon$$

Where we have 10 random effects $b \sim N(0, \sigma_b^{2})$, each of which is shared by the 3 measurements from each patient, plus 30 residuals $\epsilon \sim N(0, \sigma_{\epsilon}^{2})$, one per observation.

Unfortunately, lm cannot handle mixed-effect models, so we need to use the lme function in the nlme package.

library(nlme)

model3 <- lme(X ~ Time, random=~1|Patient, data)

The syntax of lme is essentially the same as lm, plus a random term, which I will explain below.

summary(model3)

Linear mixed-effects model fit by REML
 Data: data 
       AIC      BIC    logLik
  166.4679 171.7967 -79.23394

Random effects:
 Formula: ~1 | Patient
        (Intercept) Residual
StdDev:    5.611054 2.272793

Fixed effects: X ~ Time 
             Value Std.Error DF   t-value p-value
(Intercept) 29.980  2.086551 19 14.368206  0.0000
Time        -0.205  0.508212 19 -0.403375  0.6912
 Correlation: 
     (Intr)
Time -0.487

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.47025716 -0.48710670 -0.02399607  0.31374056  1.54213668 

Number of Observations: 30
Number of Groups: 10 

R has correctly identified that we have 30 observations but only 10 experimental units! The rest of the summary is pretty much the same as for lm, however instead of getting $R^2$ as a goodness-of-fit measure we get three measures, AIC, BIC and log-likelihood. Briefly, lower values of AIC or BIC, or higher value of log-likelihood indicate a better fit.

Random intercepts and random slopes

So, what exactly is happening here? Let’s look at the estimated coefficients

> coef(model3)

  (Intercept)   Time
A    36.58225 -0.205
B    30.95659 -0.205
C    37.90966 -0.205
D    21.31711 -0.205
E    28.61783 -0.205
F    30.38770 -0.205
G    28.20697 -0.205
H    21.53834 -0.205
I    32.60004 -0.205
J    31.68350 -0.205

R is estimating a different intercept for each patient, whilst keeping the same slope for Time. We have created a random-intercept mixed-effects model.

Why does it only change the intercept? That is because we defined the random effects as ~1 | Patient. This tells R we want to estimate a different intercept (indicated by 1) for each patient. We can also estimate different slopes for the different fixed effects in our model, by specifying them in the random effects formula. For instance

model4 <- lme(X ~ Time, random=~Time|Patient, data)
coef(model4)

  (Intercept)        Time
A    34.98618  0.66094525
B    32.97207 -1.20892331
C    35.51697  1.07422047
D    21.09373 -0.17600178
E    26.48515  0.85413897
F    33.41293 -1.72205606
G    28.04105 -0.13863582
H    20.51766  0.22698304
I    35.62665 -1.70147529
J    31.14762  0.08080454

Similarly, if we wanted a random slope model, where we change only the slope and not the intercept we would use ~Time-1 | Patient.

Finally, let’s plot the model’s prediction in the three situations

ggplot(data, aes(Time, X)) +
  geom_point(aes(col=Patient)) +
  geom_line(aes(x = Time, y = predict(model3), col=Patient))

So, which one should I choose? Random intercept, random slope, or both? It depends on your particular problem, and this is where domain-specific knowledge comes into play. Which model makes more sense in biological terms?

Mixed-effects models are a much more complex topic and there is much more you can do with them, two good reference books are “Mixed Effects Models in S and S-PLUS” by Douglas Bates and José Pinheiro and “Extending the Linear Model with R” by Julian Faraway.

And that is all for today, hope you enjoyed this post!