Field guide to linear mixed-effects models

Mixed-effects models in R, for observations that are not independent. For the Kalahari Field Course.
Author

Arpat Ozgul

Last updated

October 19, 2026

This guide was prepared for the Kalahari Field Course (EEE331) of the University of Zurich.

library(tidyverse)
library(lme4)
library(lmerTest)   # adds p-values to lmer summaries

The idea

In the previous practical, we learned to think in terms of one modelling framework:

glm(
  Response ~ Predictors,
  family = ...
)

Sections with red titles are extras They go a little beyond what you need for the main practical, so don’t worry if you do not follow every detail.

Now we add one more question:

Are my observations independent?

If they are, a GLM may be appropriate.

If they are not, we may need a mixed-effects model.

We don’t need a completely new type of model. We simply extend the model we already know.

Response ~ Predictors + (1 | Group)

Our biological question

Imagine we want to know:

Does sparrow weaver nest height affect nest temperature?

We visit 10 trees. Each tree has 6 nests, and we measure the height and temperature of every nest once.

So we have:

  • 60 nest measurements
  • but only 10 trees

The data below are simulated for teaching. Importantly, they contain a hidden tree-level driver: canopy cover differs among trees, and denser canopy makes trees cooler. Because average nest height also differs among trees, canopy cover creates an apparent relationship between nest height and temperature when tree structure is ignored.

We will initially pretend that canopy cover was not measured, and discover what the random effect does with this hidden variation. Later we will add canopy cover back into the model.

## Average nest height differs among trees
tree_mean_height <- seq(2, 5, length.out = 10)

## Nest heights also vary within each tree
within_tree_height <- c(-0.60, -0.36, -0.12,
                         0.12,  0.36,  0.60)

## A tree-level variable: canopy cover (%)
canopy_cover <- seq(20, 80, length.out = 10)

## Trees differ slightly for additional reasons we did not measure
tree_noise <- c(0.25, -0.20, 0.10, -0.15, 0.20,
               -0.10,  0.15, -0.25, 0.05, -0.05)

## Denser canopy makes trees cooler
tree_effect <- 5 - 0.10 * canopy_cover + tree_noise

## Small nest-to-nest variation
nest_noise <- c(0.20, -0.30, 0.10,
               -0.20,  0.25, -0.05)

nest_data <- expand_grid(
  Tree = 1:10,
  Nest = 1:6
) %>%
  mutate(
    Height = tree_mean_height[Tree] + within_tree_height[Nest],
    CanopyCover = canopy_cover[Tree],
    Temperature = 25 + tree_effect[Tree] + nest_noise[Nest],
    Tree = factor(Tree)
  )

nest_data
# A tibble: 60 x 5
   Tree   Nest Height CanopyCover Temperature
   <fct> <int>  <dbl>       <dbl>       <dbl>
 1 1         1   1.4         20          28.4
 2 1         2   1.64        20          28.0
 3 1         3   1.88        20          28.4
 4 1         4   2.12        20          28.0
 5 1         5   2.36        20          28.5
 6 1         6   2.6         20          28.2
 7 2         1   1.73        26.7        27.3
 8 2         2   1.97        26.7        26.8
 9 2         3   2.21        26.7        27.2
10 2         4   2.45        26.7        26.9
# i 50 more rows

Are our observations independent?

At first glance, we have 60 observations.

But are they really 60 independent observations?

But the 60 nests are not all equivalent independent replicates.

Nests from the same tree share:

  • canopy cover
  • local shade
  • tree architecture
  • exposure
  • local microclimate

So six nests on Tree 1 are not the same as six nests on six different trees.

This is non-independence. If we ignore it and pretend all 60 nests are independent, we create pseudoreplication.

Step 1: Fit the model as if all nests were independent

From the previous practical, we know how to fit a Gaussian GLM:

model_naive <- glm(
  Temperature ~ Height,
  family = gaussian,
  data = nest_data
)

summary(model_naive)
Call:
glm(formula = Temperature ~ Height, family = gaussian, data = nest_data)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  31.0223     0.3751   82.71   <2e-16 ***
Height       -1.7207     0.1027  -16.75   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.6865766)

    Null deviance: 232.500  on 59  degrees of freedom
Residual deviance:  39.821  on 58  degrees of freedom
AIC: 151.68

Number of Fisher Scoring iterations: 2

The model suggests a very strong relationship between nest height and temperature.

Before interpreting it, look at the data:

ggplot(nest_data, aes(Height, Temperature)) +
  geom_point() +
  geom_smooth(method = "lm")
`geom_smooth()` using formula = 'y ~ x'

It certainly looks as if higher nests are cooler.

If we stopped here, we might conclude:

Nest height strongly affects temperature.

But now ask the key question:

Are these observations independent?

Step 2: Look at the sampling structure

Colour the same observations by tree:

ggplot(nest_data,
       aes(Height, Temperature,
           colour = Tree)) +
  geom_point()

Now the pattern looks quite different.

Within each tree, there is almost no relationship between nest height and temperature.

Instead, some trees are generally warmer and some trees are generally cooler.

The apparent overall relationship arises because:

  • cooler trees tend to have higher nests
  • warmer trees tend to have lower nests

The naïve GLM cannot distinguish those two sources of variation.

Step 3: Add tree as a random effect

We now tell the model that nests from the same tree belong together.

model_mixed <- lmer(
  Temperature ~ Height + (1 | Tree),
  data = nest_data
)

summary(model_mixed)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Temperature ~ Height + (1 | Tree)
   Data: nest_data

REML criterion at convergence: 53.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3837 -0.8860  0.1233  0.8655  1.1620 

Random effects:
 Groups   Name        Variance Std.Dev.
 Tree     (Intercept) 4.21562  2.0532  
 Residual             0.05006  0.2237  
Number of obs: 60, groups:  Tree, 10

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) 25.034839   0.694681 11.434117  36.038 3.92e-13 ***
Height      -0.009954   0.070094 49.875114  -0.142    0.888    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
       (Intr)
Height -0.353

The important new term is:

(1 | Tree)

This means:

Allow each tree to have its own baseline temperature.

The model still asks whether Height affects temperature, but it no longer pretends that all nests are independent.

What changed?

Compare the estimated height effect:

tibble(
  Model = c("Naive GLM", "Mixed model"),
  HeightEstimate = c(
    coef(model_naive)["Height"],
    fixef(model_mixed)["Height"]
  )
)
# A tibble: 2 x 2
  Model       HeightEstimate
  <chr>                <dbl>
1 Naive GLM         -1.72   
2 Mixed model       -0.00995

And compare the model summaries:

summary(model_naive)
Call:
glm(formula = Temperature ~ Height, family = gaussian, data = nest_data)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  31.0223     0.3751   82.71   <2e-16 ***
Height       -1.7207     0.1027  -16.75   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.6865766)

    Null deviance: 232.500  on 59  degrees of freedom
Residual deviance:  39.821  on 58  degrees of freedom
AIC: 151.68

Number of Fisher Scoring iterations: 2
summary(model_mixed)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Temperature ~ Height + (1 | Tree)
   Data: nest_data

REML criterion at convergence: 53.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3837 -0.8860  0.1233  0.8655  1.1620 

Random effects:
 Groups   Name        Variance Std.Dev.
 Tree     (Intercept) 4.21562  2.0532  
 Residual             0.05006  0.2237  
Number of obs: 60, groups:  Tree, 10

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) 25.034839   0.694681 11.434117  36.038 3.92e-13 ***
Height      -0.009954   0.070094 49.875114  -0.142    0.888    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
       (Intr)
Height -0.353

In the naïve GLM, height appears to have a very strong effect.

Once we account for tree identity, the height effect is essentially gone.

The important lesson is:

The 60 observations gave us lots of measurements, but they did not give us 60 independent trees.

Why did the result change?

The naïve GLM only sees one overall pattern:

Taller nests are cooler.

The mixed model separates this into:

  • differences among trees
  • differences within trees

In this example, almost all of the apparent relationship occurs among trees. Once we account for tree identity, there is very little evidence that nest height itself affects temperature.

Fixed effects and random effects

A useful way to think about them is:

Type Role in the model Example
Fixed effect A variable whose effect we want to estimate Nest height, canopy cover
Random effect A grouping structure that allows observations from the same group to be correlated Tree

So here:

Temperature ~ Height + (1 | Tree)
  • Height is a fixed effect.
  • Tree is the random effect.

The random effect is not just a statistical correction. It represents unexplained differences among trees that make nests from the same tree more similar to one another.

What does the random effect tell us?

We can inspect the variation among trees:

VarCorr(model_mixed)
 Groups   Name        Std.Dev.
 Tree     (Intercept) 2.05320 
 Residual             0.22374

The tree-level variance tells us how much the baseline temperature differs among trees, after accounting for nest height.

A larger value means that trees still differ strongly from one another for reasons not explained by the fixed effects in the model.

You can also inspect the estimated tree effects:

ranef(model_mixed)$Tree
   (Intercept)
1    3.2286794
2    2.1175297
3    1.7548987
4    0.8433540
5    0.5306242
6   -0.4308217
7   -0.8433540
8   -1.9046024
9   -2.2672334
10  -3.0290744

Some trees are estimated to be generally warmer than average and others cooler than average.

But this raises an important biological question:

Why do the trees differ?

Can a fixed effect explain the random-effect variance?

This is the key to understanding what the random effect was doing in our example.

Remember that when we created the simulated data, trees differed in canopy cover. Trees with denser canopy were cooler. Average nest height also differed among trees.

When we fitted:

Temperature ~ Height + (1 | Tree)

we deliberately left canopy cover out.

The model therefore did not know why Tree 1 was warmer than Tree 8. It simply allowed the trees to have different baseline temperatures through (1 | Tree).

In other words, the random effect absorbed tree-level variation caused partly by the hidden canopy-cover effect.

Now look at canopy cover:

ggplot(nest_data, aes(CanopyCover, Temperature)) +
  geom_point() +
  geom_smooth(method = "lm")
`geom_smooth()` using formula = 'y ~ x'

And add it as a fixed effect:

model_canopy <- lmer(
  Temperature ~ Height + CanopyCover + (1 | Tree),
  data = nest_data
)

summary(model_canopy)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Temperature ~ Height + CanopyCover + (1 | Tree)
   Data: nest_data

REML criterion at convergence: 18.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.60145 -0.84075  0.07006  0.82811  1.31978 

Random effects:
 Groups   Name        Variance Std.Dev.
 Tree     (Intercept) 0.02438  0.1561  
 Residual             0.04997  0.2235  
Number of obs: 60, groups:  Tree, 10

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) 30.083550   0.174718 11.333309 172.183   <2e-16 ***
Height       0.011905   0.070409 49.000003   0.169    0.866    
CanopyCover -0.102504   0.004617 34.730145 -22.203   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Height
Height      -0.403       
CanopyCover -0.246 -0.763

Compare the random-effect variance before and after adding canopy cover:

VarCorr(model_mixed)
 Groups   Name        Std.Dev.
 Tree     (Intercept) 2.05320 
 Residual             0.22374 
VarCorr(model_canopy)
 Groups   Name        Std.Dev.
 Tree     (Intercept) 0.15613 
 Residual             0.22354

The variance attributed to Tree becomes much smaller.

Why?

Because we have now explained a large part of why trees differ in temperature. Variation that previously appeared as unexplained differences among trees can now be attributed to canopy cover.

Random-effect variance is variation among groups that remains unexplained by the fixed effects in the model.

This also explains the misleading result from the naïve GLM. Nest height itself did not create the strong temperature pattern. Trees differed in canopy cover, canopy cover affected temperature, and average nest height also differed among trees. Ignoring the grouping structure mixed these sources of variation together.

When should I think about random effects?

Ask whether several observations share the same sampling unit.

If you sampled… A possible random effect
Several nests on the same tree Tree
Several trees within the same site Site
Repeated measurements on the same animal Individual
Repeated detections from the same camera Camera
Several plots within a larger block Block
Repeated observations in the same year Year

The principle is:

Observations that belong to the same group may be more similar to each other than observations from different groups.

What if the sampling design is nested?

Sometimes there is more than one level of grouping.

Imagine that we measure several nests on each tree, and several trees occur within each site:

Site
 └── Tree
      └── Nest measurements

Trees are nested within sites. We can write:

(1 | Site/Tree)

which expands to:

(1 | Site) + (1 | Site:Tree)

This allows sites to have different baseline temperatures, and trees within those sites to differ further.

Use nested random effects when one grouping unit genuinely occurs within another grouping unit.

A small but important notation point: if Tree is nested within Site, the order is (1 | Site/Tree), not (1 | Tree/Site).

Random intercepts vs. random slopes

So far we have used a random intercept:

(1 | Tree)

This allows each tree to have a different baseline temperature.

Sometimes we may also expect the effect of a predictor to differ among groups. For example:

(Height | Tree)

This is shorthand for:

(1 + Height | Tree)

It allows each tree to have its own baseline temperature and its own relationship between height and temperature.

This is a random-slope model. You probably will not need one for the field-course projects, but it is useful to recognise the notation.

Do not confuse this with nesting. (Height | Tree) does not mean that Height is nested within Tree; it means that the effect of Height is allowed to vary among trees.

Testing and reporting the results

As in the first practical, fitting the model and testing a biological effect are separate steps.

For these Gaussian mixed models, the lmerTest package gives approximate p-values for the fixed effects in the model summary:

summary(model_mixed)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Temperature ~ Height + (1 | Tree)
   Data: nest_data

REML criterion at convergence: 53.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3837 -0.8860  0.1233  0.8655  1.1620 

Random effects:
 Groups   Name        Variance Std.Dev.
 Tree     (Intercept) 4.21562  2.0532  
 Residual             0.05006  0.2237  
Number of obs: 60, groups:  Tree, 10

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept) 25.034839   0.694681 11.434117  36.038 3.92e-13 ***
Height      -0.009954   0.070094 49.875114  -0.142    0.888    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
       (Intr)
Height -0.353

For the height effect, the important part is the Height row under Fixed effects. It gives the estimated effect, its standard error, a t statistic, and a p-value.

The main reporting rule is the same as before:

Describe the biological result first, then give the statistical evidence.

For example, we can extract the actual result from our model:

How would I report this?

After accounting for differences among trees, nest height was not associated with nest temperature (linear mixed-effects model: β = -0.01 ± 0.07 SE, t~49.9~ = -0.14, p = 0.888).

Notice that we do not simply write “the effect was non-significant”. We say what biological relationship was tested and then report the statistical evidence.

Now consider the model in which we measured canopy cover and added it as a fixed effect:

How would I report this?

Trees with denser canopy were cooler (linear mixed-effects model: β = -0.1 ± 0 SE, t~34.7~ = -22.2, p < 0.001). Adding canopy cover also substantially reduced the unexplained variance among trees.

That last sentence is biologically important: the fixed effect helps explain why the groups differed.

Comparing mixed models with AIC

We can also compare competing biological models using AIC, just as in the first practical.

There is one important technical point for lmer(). By default, lmer() fits models using REML (REML = TRUE). If two models differ in their fixed effects, their REML-based AIC values should not be compared directly. Refit them using maximum likelihood:

model_mixed_ML <- lmer(
  Temperature ~ Height + (1 | Tree),
  data = nest_data,
  REML = FALSE
)

model_canopy_ML <- lmer(
  Temperature ~ Height + CanopyCover + (1 | Tree),
  data = nest_data,
  REML = FALSE
)

AIC(model_mixed_ML, model_canopy_ML)
                df      AIC
model_mixed_ML   4 58.80402
model_canopy_ML  5 10.93046

The model with the lower AIC has more support, after balancing model fit against model complexity.

How would I report this?

The model including canopy cover was better supported than the model containing nest height alone (ΔAIC = 47.9), suggesting that canopy cover explains important variation in nest temperature among trees.

For this course, the main thing to remember is simple:

p-values  -> evidence for particular fixed effects
AIC       -> comparison of competing models

If you compare lmer() models with different fixed effects using AIC, fit them with REML = FALSE.

Adding a random effect itself is not a significance-testing trick. It is about representing the sampling design correctly.

A brief note about diagnostics

As with the GLMs in the first practical, a fitted mixed model should be checked before we trust it.

For a simple Gaussian mixed model, start by looking at the residuals:

plot(model_mixed)

qqnorm(resid(model_mixed))
qqline(resid(model_mixed))

With more complicated generalized mixed models (glmer()), especially Binomial or Poisson models, packages such as DHARMa can be useful.

The roadmap

The roadmap from the first practical still works. We have simply expanded the third question.

Question Why does it matter?
What is my response? Chooses the model family.
What are my predictors? Defines the fixed effects in the model.
Are the observations independent? Determines whether random effects may be needed.
How will I compare models? p-values or AIC, depending on the question.

The basic progression is therefore:

Response ~ Predictors

and, when observations belong to groups:

Response ~ Predictors + (1 | Group)

Final message

In the first practical, we learned to think about models rather than separate statistical tests.

Here we added one more idea: the model should also represent how the data were collected.

If observations from the same tree, animal, camera, site, or other sampling unit are not independent, we can represent that grouping with a random effect:

Response ~ Predictors + (1 | Group)

The random effect allows groups to differ from one another. If we later measure a biological variable that explains those differences, adding it as a fixed effect may reduce the unexplained random-effect variance.

And remember:

More measurements within the same groups give us more information about those groups, but they do not create new independent groups.