Data Structures & Their Hierarchy

What do we mean, though, by this hierarchical structure notion and why does it matter? Basic statistical methods including the classic 1-way ANOVA, simple linear and logistic regression, and the \(\chi^2\) test (along with many others!) assume independence of sampling. However, different types of hierarchical structures violate this assumption. The hierarchy can take many forms, but one can imagine such examples as pupils in the same classrooms, patients seen by the same doctor, or plants in the same field, may not necessarily be treated as independent. That is, if one were to assess the skill of a large number of pupils across many schools in the taking of some standardized test, one can imagine that students in the same class may not function as independent measurements (perhaps via some shared aspect of having the same teacher), and further, that students in the same school may also be thought of as having some correlation in their scores, maybe through experiencing shared school-level teaching curricula or resources. Performing inference on all the students across classrooms and schools without accounting for this hierarchical structure will often lead to dangerously incorrect results. For further discussion of this notion and how it can be treated, Hagadorn and Shaffer (2019) have a nice discussion of this, though their terminology is slightly different than ours.

To deal with this challenge, researchers typically turn to an approach we call a linear mixed-effects model, which is an extension of a standard linear model as described above, typically containing only fixed effects. The key addition is the notion of a random effect, defined below, which allows for the accounting of within-group correlation, while also facilitating inference about both group-specific and population-level effects.

To define a linear mixed-effects model, we can introduce a new concept, known as a random effect (also known as a varying effect). However, we immediately find ourselves in definition jail with respect to what, specifically, a random effect means relative to a fixed effect.

Random Effects vs. Fixed Effects

To start, we’ll develop some intuition for what a random effect is and why we might want to use one.Imagine we’re conducting a study investigating the effect of soil pH on crop yield. We have a number of experimental plots, which are nested within farms across different regions. Since each farm contributes multiple plots, we end up with a hierarchical data structure where observations within farms cannot be assumed independent. The set-up, visually, may look something like this:

graph TD
    Farm1[Farm 1]:::farm1
    Farm2[Farm 2]:::farm2
    Farm3[Farm 3]:::farm3
    
    Farm1 --> P1[Plot 1]:::plot
    Farm1 --> P2[Plot 2]:::plot
    Farm1 --> P3[Plot 3]:::plot
    Farm1 --> P4[Plot 4]:::plot
    Farm1 --> P5[Plot 5]:::plot
    
    Farm2 --> P6[Plot 6]:::plot
    Farm2 --> P7[Plot 7]:::plot
    Farm2 --> P8[Plot 8]:::plot
    Farm2 --> P9[Plot 9]:::plot
    Farm2 --> P10[Plot 10]:::plot
    
    Farm3 --> P11[Plot 11]:::plot
    Farm3 --> P12[Plot 12]:::plot
    Farm3 --> P13[Plot 13]:::plot
    Farm3 --> P14[Plot 14]:::plot
    Farm3 --> P15[Plot 15]:::plot
    
    classDef farm1 fill:#DEB887,stroke:#8B4513,stroke-width:5px,color:#000
    classDef farm2 fill:#F4A460,stroke:#A0522D,stroke-width:5px,color:#000
    classDef farm3 fill:#D2B48C,stroke:#8B7355,stroke-width:5px,color:#000
    classDef plot fill:#FFFFFF,stroke:#000000,stroke-width:4px,color:#000

Typically a random effect is thought of as a mechanism for capturing otherwise “unexplained” variability in your dataset that, importantly, is not central to your analysis. This is in opposition to a fixed effect which captures variation that you do in fact want to understand as an explanatory variable in your data. In this example, we can consider the farms as not being central to our analysis, and simply adding components of variation to our estimates that we may not care about, but that nonetheless we need to account for. We use the definition of Kreft and De (1998) and Gelman (2005) and define the two things with respect to some grouping variable of interest. In our example, that grouping variable is the farm. A fixed effect is constant across groups (i.e. farms) and random effects vary. In our example, a random effect across each farm \(i\), we can define as \(\alpha_i\), and the fixed effect (slope) corresponds to \(\beta\), which manifests as parallel lines for different groups. This is the model \[y_{i} = \beta + \alpha_i\] where \(y\) is our variable of interest (crop yield).

In this guide we’ll go over the formalism behind this kind of approach, when it’s useful, and how to implement a linear mixed-effects model in common programming languages.

We’ll also focus mainly on linear mixed-effects models with random intercepts Figure 2, but will discuss the existence briefly of random slopes, and thus, random effects models. The formulation will be consistent with the classic frequentist approach & language. We assume you’re familiar with the representation of fixed effects and simple linear regression in Figure 2.

Theoretical Intuition

[THINK BACK TO INTUITION FROM BEFORE]

This is exactly the use case for something like linear mixed-effects models. Such an approach allows us to account for within-group correlation while doing inference about both group-specific and population-level effects. To formalize this, say our simple linear model is

where we have measurements of crop yield, \(y_{ij}\) in plots \(i\), nested within farm \(j\), as well as some soil pH levels, \(x_{ij}\). It follows that we may want to understand if there is some effect of the fixed effect soil pH on the crop yield. That is, we might want to do some form of standard linear regression to learn about our \(\beta\) values. But, we’re faced with the possibility that each farm has characteristics specific to itself that may influence the crop yield. This could be a variety of different features, but often we can assume that some or all of them are unobserved.

Therefore, we can use a random effect to help control for this unobserved dependence.

Random Intercepts Model

Now that we have a foundation of a random effect, let us consider them in the context where there is also a fixed effect that we wish to estimate. The model, as we wrote it before, will now include our single fixed effect of interest, and our random effect of farm.

\[Y_{ij} = \beta_0 + \beta_1 X_{ij} + \alpha_j + \epsilon_{ij}\]

This is to be interpreted, in plain english, as a simple linear mixed-effects model with a fixed effect of soil pH, and a random effect on farm. In this formulation, we are assuming that our random effect takes the form of a random intercept. Which is to say, we assume the relationship between soil pH and yield to operates with the same form in each place, but that there are likely differences in soil pH generally between farm locations. So the per-unit change of \(X_{ij}\) leads to the same per-unit change in \(Y_{ij}\), and therefore we can estimate one population-level \(\beta\) value. Let’s visualize this:

In Figure 3 we can clearly see that there is a discernible relationship between points of the same colour, but, at least looking at it, we can see that the the relationship between the two variables are probably linear. We can preview what will come out of our analysis by putting some lines through each set of points.

Show the code

# we'll show how to read in this data example later on!
ggplot(data = field_data) +
    geom_point(aes(x = soil_pH, y = yield, colour = farm_id)) +
    geom_smooth(
        aes(
            x = soil_pH, y = yield,
            group = farm_id, colour = farm_id
        ),
        method = lm
    ) +
    theme_bw() +
    labs(x = "Soil pH", y = "Crop Yield") +
    scale_colour_manual("Farm ID", values = MoMAColors::moma.colors("OKeeffe"))

There are some obvious differences here in Figure 4, and it appears as though perhaps the intercepts for the different farms might be different. Having done this brief visualization, let’s look at how we can actually estimate our model.

Formalism

So that we do not lose generality, let us consider the model where there are two random variables: \(\mathcal{Y}\), the \(n\)-dimensional vector of responses, and \(\mathcal{B}\), the \(q\)-dimensional vector of random effects. We of course observe some value of \(\mathcal{Y}\), which we note as \(\mathbf{y}\), but based on how we defined a random effect as a random variable, we definitionally do not observe the value \(\mathbf{b}\) of \(\mathcal{B}\).

The unconditional distribution of \(\mathcal{B}\) (that is, the distribution before we observe any data) is given by \(\mathcal{B}\sim\mathcal{N}(\mathbf{0},\Sigma_\theta)\), where \(\Sigma_\theta\) is a variance-covariance matrix that depends on our variance-component parameter, \(\theta\). The variance-component parameter \(\theta\) controls how much variation we expect to see in our random effects. In the simplest case of a random intercept model, \(\theta\) might be a single number that determines the standard deviation of the random intercepts. In our soil pH example, this would control how much we expect the baseline crop yield to vary from field to field. In more complex models, \(\theta\) can be a vector of multiple parameters that control not just the variance of different random effects, but also how they covary with each other. For instance, if we had both random intercepts and random slopes for pH effect by field, \(\theta\) would include parameters controlling the variance of intercepts, the variance of slopes, and the correlation between them.

In this case, we end up with a linear predictor \(\mathbf{X}\beta + \mathbf{Z}\mathbf{b}\), which is the conditional mean of \(\mathcal{Y}\), given \(\mathcal{B} = \mathbf{b}\). Since we are dealing with the linear mixed-effects model case for now, the conditional distribution \((\mathcal{Y}|\mathcal{B}=\mathbf{b})\) is given by \[ (\mathcal{Y}|\mathcal{B}=\mathbf{b}) \sim \mathcal{N}(\mathbf{X}\beta + \mathbf{Z}\mathbf{b}, \sigma^2 \mathbf{I}),\] where \(\mathbf{I}\) is the identity matrix.

To make the computational machinery work more efficiently, we introduce what are called spherical random effects, denoted \(\mathcal{U}\). The term “spherical” here refers to random effects that have been standardized so that they have no correlation with each other and all have the same variance. Think of this as similar to converting raw test scores to z-scores. These spherical random effects are related to our actual random effects through the transformation \[\mathcal{B}=\Lambda_\theta\mathcal{U},\] where \(\Lambda_\theta\) is called the relative covariance factor. This transformation essentially allows us to work with a simpler, standardized version of the problem during computation, then transform back to get results on the original scale. In our soil pH example, if we had random effects for different fields that were correlated (perhaps because nearby fields share soil characteristics), the spherical random effects would represent a transformed version where those correlations have been removed.

This is where the rubber meets the road computationally. We arrive at the challenge of computing the maximum likelihood estimation of the variance-component parameter \(\theta\), which we term \(\hat{\theta}\), that minimizes what is called the profiled log-likelihood. The actual value of \(\hat{\theta}\) is obtained through numerical optimization, the details of which will depend on the software package at hand. In the process of evaluating our profiled log-likelihood and finding this optimal \(\hat{\theta}\), we simultaneously obtain several other quantities: \(\hat{\beta}\), which is our estimate for the coefficients of the fixed effects, and \(\tilde{\mathbf{u}}_{\hat{\theta}}\), which is the conditional mode of the spherical random effects \(\mathcal{U}\) (that is, the value of \(\mathbf{u}\) that is most likely given the data we observed).

We then transform these conditional modes of the spherical random effects back to the original scale of our random effects using \[\tilde{\mathbf{b}}_{\hat{\theta}} = \Lambda_{\hat{\theta}}\tilde{\mathbf{u}}_{\hat{\theta}}.\] These values \(\tilde{\mathbf{b}}_{\hat{\theta}}\) represent the values of \(\mathbf{b}\) that maximize the density of the conditional distribution \(\mathcal{B}|\mathcal{Y} = \mathbf{y}\). In plain language, these are our best guesses for the random effect values given the data we actually observed. You will often see these called the best linear unbiased predictors (BLUPs) or conditional modes. In our soil pH example, these would be the estimated effect of each specific field on crop yield. For a linear mixed model where all distributions are Gaussian (normal), these conditional modes are also equal to the conditional means, which is a convenient mathematical property.

Estimation via (Restricted) Maximum Likelihood

When we use linear mixed-effects models, we often fit the models via a maximum likelihood estimation procedure, with the most common form of it being restricted maximum likelihood (REML). More details about this method along with some notes to pay attention to are in the appendix on REML

Why restricted? Maximum likelihood estimation of variance components in mixed-effects models suffers from a well-known deficiency: the estimates are biased downward because the procedure fails to account for the degrees of freedom consumed in estimating the fixed effects, much as the sample variance requires division by \(n-1\) rather than \(n\) to correct for the estimation of the mean. Restricted maximum likelihood remedies this problem by basing inference on \(n-p\) error contrasts—linear combinations of the data that contain no information about the fixed effects—thereby properly accounting for the uncertainty introduced by estimating \(\beta\). While REML produces less biased estimates of variance components and is thus preferred for final model fitting, it cannot be used to compare models with different fixed effects structures, as the restricted likelihoods under different fixed effects specifications involve different transformations of the data and are therefore not comparable. Consequently, the recommended practice is to use ordinary maximum likelihood when selecting among competing fixed effects structures via likelihood ratio tests or information criteria, and then to refit the chosen model using REML to obtain optimal estimates of the variance components.

Random Intercepts Model

Recall our soil pH and crop yield example where we showed a quick plot in Figure 4 indicating that there were almost certainly differences between our farms in terms of their intercepts. So, we can simply fit a linear mixed-effects model, with the goal of estimating a (fixed) effect of soil pH on crop yield (i.e., what is the “slope” of the relationship), and a random effect on the intercept (i.e., what are the intercepts of each farm). Quickly, let’s set up our model.

In our soil pH example, our statistical model is \[Y_{ij} = \beta_0 + \beta_1 X_{ij} + \alpha_{j} + \epsilon_{ij}\] where \(Y_{ij}\) is the measurement of crop yield at plot \(i\) in farm \(j\), \(\beta_0\) is our population-level intercept, \(\beta_1\) is the to-be-estimated effect or coefficient of \(X_{ij}\), the measurement of soil pH at plot \(i\) in farm \(j\), \(\alpha_j\) is the random effect on the intercept for each group (farm), \(j\), and then \(\epsilon\) is our error term.

We need to translate this into code. In these languages, we use what’s called a formula to describe our statistical model. The general form of that formula is:

response ~ 1 + fixed_effect + (1 | random_intercept)

where we can add arbitrarily many fixed effects, e.g.,

response ~ 1 + fixed_effect_1 + fixed_effect_2 + (1 | random_intercept_1) + (1 | random_intercept_2)

Now, to some code! Let’s first get a sense of what our data actually look like for this example.

# some set up!
library(readr)
library(here)
# load in the data from our local location
field_data <- read_csv(
    here("./research-guides/mixed-effects-models/assets/data/crop_yield_data.csv"),
    show_col_types = FALSE
)
head(field_data)

# A tibble: 6 × 4
  farm_id plot_id soil_pH yield
    <dbl>   <dbl>   <dbl> <dbl>
1       1       1    7.06 143. 
2       1       2    6.58 106. 
3       1       3    6.32 122. 
4       1       4    6.39 138. 
5       1       5    6.41 112. 
6       1       6    6.19  92.7

# some set up!
using DataFrames, CSV, CategoricalArrays, Base.Filesystem
#  project_root = @__DIR__ # if we need to double check where the directory is pointing rn
# load the data from local location
field_data = CSV.read(
    joinpath(@__DIR__, "assets", "data", "crop_yield_data.csv"),
    DataFrame
)

150×4 DataFrame
 Row │ farm_id  plot_id  soil_pH  yield
     │ Int64    Int64    Float64  Float64
─────┼─────────────────────────────────────
   1 │       1        1  7.06464  143.0
   2 │       1        2  6.57797  106.061
   3 │       1        3  6.32341  122.438
   4 │       1        4  6.38518  137.577
   5 │       1        5  6.41172  112.132
   6 │       1        6  6.19462   92.6556
   7 │       1        7  7.32745  135.555
   8 │       1        8  5.90546  118.104
  ⋮  │    ⋮        ⋮        ⋮        ⋮
 144 │       5      144  5.40329   65.022
 145 │       5      145  5.96985   81.1347
 146 │       5      146  6.85269   77.3114
 147 │       5      147  6.1639    65.7406
 148 │       5      148  6.52858   64.3481
 149 │       5      149  6.86806   69.387
 150 │       5      150  6.73596   87.4023
                           135 rows omitted

# MixedModels.jl requires any grouping variable (anything behind a `|`) to be
# categorical. lme4 coerces this for us, Julia does not, so we do it by hand
field_data.farm_id = categorical(field_data.farm_id)

150-element CategoricalArray{Int64,1,UInt32}:
 1
 1
 1
 1
 1
 1
 1
 1
 1
 1
 ⋮
 5
 5
 5
 5
 5
 5
 5
 5
 5

field_data

150×4 DataFrame
 Row │ farm_id  plot_id  soil_pH  yield
     │ Cat…     Int64    Float64  Float64
─────┼─────────────────────────────────────
   1 │ 1              1  7.06464  143.0
   2 │ 1              2  6.57797  106.061
   3 │ 1              3  6.32341  122.438
   4 │ 1              4  6.38518  137.577
   5 │ 1              5  6.41172  112.132
   6 │ 1              6  6.19462   92.6556
   7 │ 1              7  7.32745  135.555
   8 │ 1              8  5.90546  118.104
  ⋮  │    ⋮        ⋮        ⋮        ⋮
 144 │ 5            144  5.40329   65.022
 145 │ 5            145  5.96985   81.1347
 146 │ 5            146  6.85269   77.3114
 147 │ 5            147  6.1639    65.7406
 148 │ 5            148  6.52858   64.3481
 149 │ 5            149  6.86806   69.387
 150 │ 5            150  6.73596   87.4023
                           135 rows omitted

We’re dealing with clean data in this case, so we’re skipping some common data cleaning steps that may normally come here.

library(lme4, quietly = TRUE)
library(ggplot2)

# random intercepts model
model_ri <- lmer(
    yield ~ 1 + soil_pH + (1 | farm_id),
    data = field_data,
    REML = TRUE
)

# this shows the model summary
summary(model_ri)

Linear mixed model fit by REML ['lmerMod']
Formula: yield ~ 1 + soil_pH + (1 | farm_id)
   Data: field_data

REML criterion at convergence: 1184.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.77492 -0.67190  0.03397  0.78095  2.55810 

Random effects:
 Groups   Name        Variance Std.Dev.
 farm_id  (Intercept) 405.1    20.13   
 Residual             146.5    12.10   
Number of obs: 150, groups:  farm_id, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)   -7.407     15.921  -0.465
soil_pH       16.034      2.002   8.007

Correlation of Fixed Effects:
        (Intr)
soil_pH -0.823

using MixedModels, StatsBase, StatsModels

# random intercepts model (structure essentially identical)
model_ri = fit(
    MixedModel,
    @formula(yield ~ 1 + soil_pH + (1 | farm_id)),
    field_data,
    REML=true
) # model summary displayed automatically

Linear mixed model fit by REML
 yield ~ 1 + soil_pH + (1 | farm_id)
 REML criterion at convergence: 1184.4262456837705

Variance components:
            Column   Variance Std.Dev.
farm_id  (Intercept)  405.0656 20.1262
Residual              146.5110 12.1042
 Number of obs: 150; levels of grouping factors: 5

  Fixed-effects parameters:
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  -7.40672    15.9206   -0.47    0.6418
soil_pH      16.0336      2.00237   8.01    <1e-14
──────────────────────────────────────────────────

# get out the variance components
var_components <- as.data.frame(VarCorr(model_ri))
farm_var <- var_components$vcov[1]
residual_var <- var_components$vcov[2]

# calculate icc manually
icc <- farm_var / (farm_var + residual_var)
cat("Intraclass Correlation Coefficient:", round(icc, 3), "\n")

Intraclass Correlation Coefficient: 0.734 

# get out the variance components
vc = VarCorr(model_ri);
farm_var = vc.σρ.farm_id.σ[1]^2;
residual_var = vc.s^2;

# calculate icc manually
icc = farm_var / (farm_var + residual_var);
println("Intraclass Correlation Coefficient: ", round(icc, digits=3))

Intraclass Correlation Coefficient: 0.734

# residuals and fitted values
field_data$residuals <- residuals(model_ri)
field_data$fitted <- fitted(model_ri)

# diagnostic plots
par(mfrow = c(2, 2))

# Q-Q plot for normality of residuals
qqnorm(field_data$residuals, main = "Q-Q Plot of Residuals")
qqline(field_data$residuals, col = "red")

# homoscedasticity check
plot(field_data$fitted, field_data$residuals,
    xlab = "Fitted Values", ylab = "Residuals",
    main = "Residuals vs Fitted Values"
)
abline(h = 0, col = "red", lty = 2)

# residuals
hist(field_data$residuals,
    breaks = 20,
    main = "Distribution of Residuals",
    xlab = "Residuals"
)

# residuals vs predictor
plot(field_data$soil_pH, field_data$residuals,
    xlab = "Soil pH", ylab = "Residuals",
    main = "Residuals vs Soil pH"
)
abline(h = 0, col = "red", lty = 2)

par(mfrow = c(1, 1))

using Plots, Distributions, StatsPlots

# extract residuals and fitted values
field_data.residuals = residuals(model_ri);
field_data.fitted = fitted(model_ri);

# Q-Q plot of residuals
p1 = plot(qqplot(Normal, field_data.residuals),
    xlabel="Theoretical Quantiles",
    ylabel="Sample Quantiles",
    title="Q-Q Plot of Residuals",
    legend=false);

# homoscedasticity 
p2 = scatter(field_data.fitted, field_data.residuals,
    xlabel="Fitted Values",
    ylabel="Residuals",
    title="Residuals vs Fitted Values",
    legend=false,
    alpha=0.6);
hline!([0], color=:red, linestyle=:dash);

# hist of residuals
p3 = histogram(field_data.residuals,
    xlabel="Residuals",
    ylabel="Frequency",
    title="Distribution of Residuals",
    legend=false,
    bins=20);

# residuals vs predictor
p4 = scatter(field_data.soil_pH, field_data.residuals,
    xlabel="Soil pH",
    ylabel="Residuals",
    title="Residuals vs Soil pH",
    legend=false,
    alpha=0.6);
hline!([0], color=:red, linestyle=:dash);

plot!(p1, p2, p3, p4, layout=(2, 2), size=(1100, 900))

# get random effects
random_effects <- ranef(model_ri)$farm_id
names(random_effects) <- "intercept"

# Q-Q plot on random effects
par(mfrow = c(1, 2))

qqnorm(random_effects$intercept,
    main = "Q-Q Plot of Random Intercepts"
)
qqline(random_effects$intercept, col = "red")

# hist of random effects
hist(random_effects$intercept,
    breaks = 10,
    main = "Distribution of Random Intercepts",
    xlab = "Random Intercept"
)

par(mfrow = c(1, 1))

# extract via ranef but then turn from Vector{Matrix{Float64}} to just a flat
# vector for easy plotting
random_effects = vec(reduce(vcat, ranef(model_ri)));
# Q-Q plot for random effects
p1 = plot(qqplot(Normal, random_effects),
    xlabel="Theoretical Quantiles",
    ylabel="Sample Quantiles",
    title="Q-Q Plot of Random Intercepts",
    legend=false);
# hist of random effects
p2 = histogram(random_effects,
    xlabel="Random Intercept",
    ylabel="Frequency",
    title="Distribution of Random Intercepts",
    legend=false,
    bins=10);

plot(p1, p2, layout=(1, 2), size=(800, 350))

# get a sequence of pH values for prediction
ph_range <- seq(min(field_data$soil_pH),
    max(field_data$soil_pH),
    length.out = 100
)

# Population-level (fixed effects) prediction
pop_pred <- predict(model_ri,
    newdata = data.frame(soil_pH = ph_range),
    re.form = NA
)

# plot
ggplot(field_data, aes(x = soil_pH, y = yield)) +
    # put the individual farms
    geom_point(aes(color = factor(farm_id)), alpha = 0.6) +
    # now farm-specific lines
    geom_line(aes(y = fitted, group = farm_id, color = factor(farm_id)),
        linewidth = 0.8, alpha = 0.7
    ) +
    # now the trend overall
    geom_line(
        data = data.frame(soil_pH = ph_range, yield = pop_pred),
        aes(x = soil_pH, y = yield),
        color = "black", linewidth = 1.5, linetype = "dashed"
    ) +
    labs(
        x = "Soil pH",
        y = "Crop Yield (kg/ha)",
        title = "Random Intercepts Model: Crop Yield vs Soil pH",
        subtitle = "Dashed line: population average; Colored lines: farm-specific predictions",
        color = "Farm"
    ) +
    theme_minimal() +
    theme(legend.position = "right")

# get a sequence of pH values for prediction
ph_range = collect(range(minimum(field_data.soil_pH),
    maximum(field_data.soil_pH),
    length=100));
# Population-level (fixed effects) prediction
pop_pred_data = DataFrame(
    soil_pH=ph_range,
    farm_id=repeat(["Farm_01"], Base.length(ph_range)),
    yield=zeros(Base.length(ph_range))
);
pop_pred = predict(model_ri, pop_pred_data, new_re_levels=:population);
# plot
p = scatter(field_data.soil_pH, field_data.yield,
    group=field_data.farm_id,
    xlabel="Soil pH",
    ylabel="Crop Yield (kg/ha)",
    title="Random Intercepts Model: Crop Yield vs Soil pH",
    alpha=0.6,
    markersize=4,
    legend=:outerright,
    label=permutedims(unique(field_data.farm_id)));

# add farm-specific fitted lines
for farm in unique(field_data.farm_id)
    farm_data = field_data[field_data.farm_id.==farm, :]
    plot!(farm_data.soil_pH, farm_data.fitted,
        label="",
        linewidth=2,
        alpha=0.7)
end

# add population-level trend finally
plot!(ph_range, pop_pred,
    color=:black,
    linewidth=2,
    linestyle=:dash,
    label="Population Average",
    size=(1100, 900))

When should I use a mixed model instead of traditional regression?

Use mixed models when your data has a natural grouping or clustering structure and observations within groups are likely correlated. Examples include repeated measures on the same subjects, students within schools, or spatial clusters. You cannot just ignore this variation. If you ignore this structure, your standard errors will be wrong and your inferences invalid!

I’m working in R/Julia and want to know more about some of the code I’m trying to use

You should check out the excellent documentation that exits for both resources. The lme4 package has the docs of course, but it also has a paper attached to it (Bates et al. 2014) which also may be helpful. The Julia MixedModels.jl docs maintained by Alday and Bates (2025) are excellent.

How do I decide what to include as random effects?

An age old question! Often, depending on what kind of inference you want to make, the same effects could be considered either fixed OR random. Start by identifying the grouping structure in your data. Include random intercepts for grouping factors where you expect baseline differences across groups. Be parsimonious: more complex random effects structures require more data to estimate reliably. A common rule of thumb is to have at least 5-10 groups per random effect parameter.

I have a random effect on my intercept, but I think that within my groups, the relationship between my fixed effect and my response will be different. What should I do?

You are no longer talking about a fixed effect! That, referencing Figure 2, you are now interested in a purely random effects model. Proceed with the fitting of the model the same way you would otherwise in most cases, after consulting a suitable text.

What’s the difference between ML and REML estimation?

REML (Restricted Maximum Likelihood) is generally preferred because it produces unbiased estimates of variance components. However, you can’t use it for model comparison. Use ML when comparing models with different fixed effects structures using likelihood ratio tests, then refit the final model with REML. For comparing models with different random effects, either ML or REML is appropriate.

What if I don’t have normal data (e.g. binomial, count data, etc)?

Then you’re looking for a generalized linear mixed-effects model (GLMM). The concepts underpining a GLMM are very similar with respect to the effects and the intuition, but require some extra care. Fitting them can be challenging, though many excellent resources exist to help you fit them out of the box. Bolker et al. (2009) is a very accessible introduction to the topic.

How do I assess model fit for mixed models?

Unlike OLS regression, \(R^2\) is more complex for mixed models. Use conditional \(R^2\) (variance explained by fixed and random effects) and marginal \(R^2\) (variance explained by fixed effects only). Information criteria (AIC, BIC) are useful for model comparison. For validation, consider cross-validation at the group level.

Can I include both spatial and temporal correlation structures?

Yes, but this requires careful consideration. Packages like nlme in R allow you to specify complex correlation structures within groups. However, computational demands increase substantially. For spatiotemporal data, consider whether the spatial correlation is adequately captured by random effects or whether you need explicit correlation structures.

What assumptions are required to fit this kind of model?

Many of the assumptions of a linear mixed effect model are the same as a standard linear model, which are covered in the OLS regression research guide. Briefly, the most important assumptions are:

1. Linearity: The relationships between the predictors and the response variable are linear
2. Independence of residuals: Conditional on the fixed and random effects, the residuals are independent of each other.
3. Normality of residuals: The residuals are normally distributed.
4. Constant variance of residuals (Homoscedasticity): The spread of the residuals is consistent across all levels of the predictors.
5. No multicollinearity: The independent variables are not strongly correlated with one another.
6. Normal distribution of random effects: The random effects are assumed to be normally distributed

Checking these (mostly) occurs after the model is fit or at the experimental design stage, so we’ll focus exclusively on the assumptions that can be checked once a model has been fit or during the fitting process.

I don’t know if my effect should be a fixed effect or a random effect.

This is challenge often faced by researchers and the technicalities scope beyond this guide. This can be answered at multiple levels of complexity. We encourage any readers to use the answer (or references) that answers the question to a sufficient level of detail that they are comfortable accounting for. There is genuinely no “wrong” approach here, more just additional levels of complexity that one can incorporate into their thinking and modeling.

Beginner:

1. If you are trying to measure a “population-level” effect, then treat the thing as a fixed effect. If it represents variation that should be explained so you can focus on your effect of interest, treat it as a random effect.

Intermediate:

2. Use a Hausman test to look at the difference between the results when you use a fixed vs. random effects model and whether or not they’re significant. If what you care about is a model that best describes your data under the set of constraints you’re working with (i.e. you only have information about a subset of variables on hand), that can help you pick the designation (however see Fielding (2004), Bell, Fairbrother, and Jones (2019), and #3)

Advanced:

3. If you’re an econometrics person, resist the urge to use a Hausman test and read Bell, Fairbrother, and Jones (2019). Fixed effects vs. random effects isn’t really about whether predictors correlate with unobserved effects, it’s about whether you treat level-2 entities as unrelated fixed dummies versus draws from a distribution. Treating anything as a random draw almost always allows for richer possibilities, and we would encourage the use of the Hausman test to start (and end!) with testing whether within and between effects differ. If they do, we encourage using a within-between (REWB or Mundlak model) which gives you everything a fixed-effect model but with the added ability to estimate between effects, include time-invariant variables, and add random slopes (Mundlak 1978).

Some Suggestions for Core Textbooks:

• Gelman, Andrew, and Jennifer Hill. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press, 2006.
• Pinheiro, José C., and Douglas M. Bates. Mixed-Effects Models in S and S-PLUS. Springer, 2000.
• Fitzmaurice, Garrett M., Nan M. Laird, and James H. Ware. Applied Longitudinal Analysis. 2nd ed. Wiley, 2011.
• McCulloch, Charles E., and Shayle R. Searle. Generalized, Linear, and Mixed Models. Wiley, 2001.
• Gelman, Andrew, et al. Bayesian data analysis. Chapman and Hall/CRC, 1995.(an updated version to February 2025 exists in pdf form)

Appendix A - REML