Why Grouping Structure Matters
The pseudoreplication problem
Imagine you want to know whether nitrogen deposition increases stream algal biomass. You sample four streams — two receiving high deposition, two receiving low — and at each stream you collect ten samples per visit over six visits. That looks like \(4 \times 10 \times 6 = 240\) data points.
But those 240 observations are not 240 independent pieces of
information. Samples from the same stream on the same day share
water chemistry, discharge, and temperature. Visits to the same stream
over time share the same substrate and riparian canopy. Treating all 240
rows as independent (i.e., plugging them into lm()) gives
you wildly over-confident standard errors and p-values — this is
pseudoreplication.
📌 The mixed-effects solution
Mixed-effects models include both:
- Fixed effects — the predictors you care about and want to make inferences on (nitrogen deposition, temperature, treatment)
- Random effects — grouping factors that introduce non-independence (stream site, watershed, plot, individual, year)
Random effects absorb the variance attributable to groups, leaving the fixed-effect estimates clean and the standard errors honest.
Fixed versus random: how to decide
| Effect type | Examples | Question | Inference scope |
|---|---|---|---|
| Fixed | Nitrogen treatment (high vs. low) | Does N deposition change algal biomass? | Only these two levels |
| Fixed | Temperature (°C), canopy cover (%) | Does temperature predict DO? | Continuous gradient |
| Random | Stream site ID (sampled from many possible sites) | How much does biomass vary among sites? | Population of all streams |
| Random | Sampling year (drawn from many possible years) | How much does biomass vary among years? | Population of all years |
| Random | Individual plant or animal (sub-sampled repeatedly) | How much does the individual vary? | Population of all individuals |
✅ The practical rule
If you have fewer than ~5 levels and they exhaust all possibilities (e.g., two treatment groups, three land-use types), treat the factor as fixed. If you have many levels that represent a random draw from a larger population (sites, years, individuals, blocks), treat it as random.
Simulate a motivating dataset
We will use a single running example throughout this module: dissolved oxygen (DO, mg/L) measured at multiple stream monitoring stations as a function of water temperature (°C). Temperature is the fixed predictor; station is the grouping variable whose structure changes across sections.
set.seed(42)
n_stations <- 8
n_obs <- 12 # observations per station
# Station-level intercept offsets (random intercepts)
station_intercepts <- rnorm(n_stations, mean = 0, sd = 1.8)
# Station-level slope offsets (random slopes, used later)
station_slopes <- rnorm(n_stations, mean = 0, sd = 0.07)
stations <- paste0("S", 1:n_stations)
do_data <- lapply(seq_along(stations), function(i) {
temp <- runif(n_obs, 10, 24)
do <- (12.0 + station_intercepts[i]) +
(-0.35 + station_slopes[i]) * temp +
rnorm(n_obs, 0, 0.5)
data.frame(
station = stations[i],
temp_C = temp,
DO_mgl = do,
slope_offset = station_slopes[i] # stored for illustration
)
}) |>
bind_rows() |>
mutate(station = factor(station))
head(do_data)## station temp_C DO_mgl slope_offset
## 1 S1 15.43352 11.160632 0.1412897
## 2 S1 19.59238 10.985931 0.1412897
## 3 S1 10.05528 13.316682 0.1412897
## 4 S1 21.66083 9.731652 0.1412897
## 5 S1 10.10268 12.230557 0.1412897
## 6 S1 12.90723 10.892272 0.1412897ggplot(do_data, aes(temp_C, DO_mgl, colour = station, group = station)) +
geom_point(alpha = 0.5, size = 1.8) +
geom_smooth(method = "lm", se = FALSE, linewidth = 0.9, alpha = 0.85) +
scale_colour_brewer(palette = "Dark2", name = "Station") +
labs(x = "Water temperature (°C)", y = "Dissolved oxygen (mg/L)",
title = "Station-specific DO–temperature relationships",
subtitle = "Parallel slopes suggest random intercepts are sufficient; some slope variation is also visible") +
theme_minimal(base_size = 13) +
theme(
legend.position = "right",
plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
plot.subtitle = element_text(colour = "#7a7a72", size = 10)
)
Dissolved oxygen versus water temperature at eight stream monitoring stations. Each line is a separate OLS fit for one station. The negative temperature–DO relationship is consistent across sites, but the intercept (baseline DO level) varies markedly — some stations are consistently higher or lower than others. This station-level variation is exactly what random intercepts will absorb.
