Worked solutions — Mixed Models & Model Selection

Exercise 1 — Streams: complete pooling vs no pooling vs partial pooling

(a) Equations and parameter counts.

• lm1: $y_{ij} = \beta_0 + \varepsilon_{ij}$, $\varepsilon \sim \mathcal{N}(0, \sigma^2)$. Two parameters: $\beta_0, \sigma^2$.
• lm2: $y_{ij} = \alpha_j + \varepsilon_{ij}$ with $j = 1, \ldots, 100$ and a free $\alpha_j$ per stream. 101 parameters (100 means + $\sigma^2$).
• lm3: $y_{ij} = \beta_0 + b_{0,j} + \varepsilon_{ij}$, $b_{0,j} \sim \mathcal{N}(0, \sigma_\text{stream}^2)$. Just 3 parameters regardless of the number of streams ($\beta_0, \sigma_\text{stream}^2, \sigma^2$).

(b) Boxplot of residuals from lm1 by stream. Entire streams sit systematically above or below zero. That is the visual fingerprint of a violated independence assumption — the three replicate observations within a stream are far more similar to each other than to observations in other streams.

boxplot(split(resid(lm1), Streams$Stream), ylab="Residual", xlab="Stream")
abline(0, 0, lwd=3)

(c) Shrinkage in action. Comparing coef(lm2) with coef(lm3) stream by stream, every lm3 estimate is shifted slightly toward the grand mean (63.26). The shift is 1.49% of the deviation for every stream — the same percentage for all 100, not larger for the extreme ones. The design is balanced (exactly 3 observations per stream), and the shrinkage weight depends only on $n_j$:

It is small because among-stream variance (459) dwarfs the uncertainty in a 3-observation mean (≈6.9): each stream's own data are precise relative to how much streams truly differ, so the model barely pulls them. Stream 11 (the lowest, 8.92) moves to 9.73; Stream 43 (the highest, 111.69) moves to 110.97.

# Compare side by side
data.frame(
 noPool   = coef(lm(Density ~ factor(Stream) - 1, data=Streams)),
 partPool = coef(lm3)$Stream[ , "(Intercept)"]
)

(d) Variance components. Approximate output:

Random effects:
 Formula: ~1 | Stream
 (Intercept) Residual
StdDev: 21.43 4.56

intervals(lm3)
 lower est. upper
sd((Intercept)) 18.60 21.43 24.68
sigma 3.6 4.6 5.8

ICC = $21.43^2 / (21.43^2 + 4.56^2) \approx 0.957$. 97% of the variance is between streams; only 3% is the within-stream replicate noise. Ignoring the random effect (using lm1) would give standard errors that are far too small — the effective sample size is closer to 18 (streams) than to 54 (observations).

Exercise 2 — Length–weight random-intercept model

(a) Equation. $W_{ij} = \beta_0 + b_{0,j} + \beta_1 L_{ij} + \varepsilon_{ij}$, $b_{0,j} \sim \mathcal{N}(0, \sigma_S^2)$, $\varepsilon_{ij} \sim \mathcal{N}(0, \sigma^2)$. Biologically: $\beta_1$ is the population-average slope of weight on length; $\beta_0$ is the population-average intercept; $b_{0,j}$ is how much subject $j$'s expected weight at $L=0$ deviates from the population mean (a "scale" effect of subject identity).

(b) Fixed effects from summary(m1). From the current data file ($n = 100$, 10 subjects):

Fixed effects: Weight ~ Length
                  Value  Std.Error  DF   t-value
(Intercept)     -0.0267     0.0956  89    -0.279
Length           3.0172     0.0296  89   101.975

Population-average slope $\approx 3.02$ per unit length, with intervals(m1) giving (2.958, 3.076) — nowhere near zero. Note the Length/Weight columns are on a transformed/index scale (ranges $\approx 0.7$–3.0 and 1.4–9.4), so the slope is in those units; do not report it in grams. The intercept itself is biologically peculiar (length = 0) — don't over-interpret it.

(c) Among-subject SD. $\approx 0.211$. So two subjects of the same length are expected to differ in weight by about $\pm 0.21$ (1 SD), or roughly $\pm 0.42$ (2 SD covers most of the population), on the weight scale.

(d) Within-subject residual SD. ~0.21. ICC = $0.211^2 / (0.211^2 + 0.206^2) \approx 0.51$. Even after accounting for length, about half the remaining variation is between subjects — individual fish sit persistently above or below the population line, so repeated measures on a fish are strongly non-independent. Grading note: a student reporting ICC near 0.5 and concluding the random intercept is doing real work is right; the covariate does not explain away the among-subject variation here.

(e) Random-slope test.

m_slope <- lmer(Weight ~ Length + (1 + Length | Subject), data=LenW)
anova(m1b, m_slope, refit=FALSE) # REML LRT on random structure

Typical output shows a small LRT statistic (often $\chi^2 < 1$ on 2 df), $p$-value > 0.5, AND a singular-fit warning indicating the slope variance is at the boundary. The data don't support a random slope. We use refit = FALSE because we're comparing models that differ in random structure only, with identical fixed effects — REML log-likelihoods are valid for that comparison.

Exercise 3 — Von Bertalanffy nonlinear mixed model

(a) Why expect $L_\infty$ to vary? Asymptotic body size is heritable and depends on individual conditions across life history — food, water temperature, parasites, genetics. Forcing all individuals to share the same $L_\infty$ is the same mistake as forcing all streams to share the same mean density.

(b) Diagnostic. boxplot(split(residuals(m_fixed), AgeLen$Subject)) shows whole subjects sitting consistently above or below the zero line — the within-subject residuals are not independent. Same diagnostic as the streams example.

(c) Residual SE comparison.

summary(m_fixed)$sigma # ~ 9.5 (made-up but typical)
summary(m_nlme)$sigma # ~ 3.5
VarCorr(m_nlme) # SD(Linf | Subject) ~ 8.5

The within-subject residual SE of m_nlme (≈ 3.5) is much smaller than the nls residual SE (≈ 9.5) because the random effect absorbed the among-subject variation in $L_\infty$. By the variance partition $\sigma^2_\text{nls} \approx \sigma_{L_\infty}^2 + \sigma^2_\text{within}$, the within-subject residual has to be smaller — that is the whole point of accounting for the random effect.

(d) What does m_nlme$coefficients$fixed[1] + m_nlme$coefficients$random$Subject produce? The BLUP of each subject's $L_\infty$: the population mean plus the partially-pooled subject-specific deviation. Equivalent to the BLUPs you'd get from coef(m_nlme) for the $L_\infty$ component.

Compared with running nls separately per subject, the BLUPs:

• Stabilize estimates for subjects with sparse age coverage by borrowing strength from the rest.
• Use a Normal prior on $L_\infty$ implicitly (the random-effects distribution).
• Give consistent uncertainty bounds and let you do further inference (predict for a new individual, etc.).

Exercise 4 — Trophic position with crossed random effects

(a) Random vs fixed for Lake and Fish_Species. Both are exchangeable groupings of observations (multiple fish per lake; multiple fish per species). If the goal is to generalize the length-trophic-position relationship to other lakes or other species, random effects are appropriate. Counter-argument: only 3 species and 6 lakes is on the edge of identifiability for a random-effect variance; species-specific feeding ecologies (S1 vs S2 vs S3) are biologically distinct enough that we may care about them by name — a case for fixed effects on species and possibly on lake.

(b) Bivariate Normal interpretation.

Each species draws both an intercept deviation and a slope deviation from this distribution; they can be correlated.

(c) AICc table and interpretation of M8.

aictab(list(M1=M1,M2=M2,M3=M3,M4=M4,M5=M5,M6=M6,M7=M7,M8=M8))
#    K   AICc  Delta_AICc  AICcWt
# M8 7   36.33       0.00    0.84
# M2 9   39.64       3.31    0.16
# M1 5   82.81      46.48    0.00
# M7 7   86.44      50.11    0.00
# M5 6  273.54     237.21    0.00
# M3 4  281.69     245.36    0.00
# M4 4  463.04     426.71    0.00
# M6 6  467.21     430.88    0.00
# Ordering: M8 < M2 < M1 < M7 < M5 < M3 < M4 < M6

summary(M8)
Random effects:
 Groups Name Variance Std.Dev. Corr
 Lake (Intercept) 0.205 0.453
 Fish_Species (Intercept) 0.867 0.931
 Z_Length 0.025 0.157 1.00
 Residual 0.050 0.224
Fixed effects:
 Estimate Std.Error
(Intercept) -0.0009 0.569
Z_Length 0.422 0.092

Almost all the slope variability lives across species (different feeding ecologies → different length-TP slopes). Lake variability mostly shifts the intercept (a lake-wide trophic-position baseline). Note that M8 has no lake slope at all — the VarCorr block lists Lake with only (Intercept), so "across lakes?" is answered by the term's absence rather than by a number in the table. The ranking backs this up: M2 is M8 plus a lake slope and lands 3.31 AICc worse, and M2's own lake slope SD is 0.024 versus a species slope SD of 0.157. Within-fish residual is small. Note the slope–intercept correlation of +1.00 for the species (intercept, slope) — a singular fit, unsurprising given only 3 species. That is a degenerate (boundary) estimate — lmer has run out of data to identify the correlation and has pinned it at the edge of its range. Treat it as a warning that the random-slope structure is over-specified for this dataset, not as a real biological correlation.

(d) Practical risk with few levels. Variance components are essentially unidentifiable, singular fits are common, and any inference about the variance (CIs, LRTs) is poorly calibrated. Switch to fixed effects when (i) ≤ 3–4 levels, (ii) you care about each level by name, or (iii) the "draw from a larger population" interpretation is not credible. In this dataset, a model with fixed effects for Species and an interaction Length × Species — with Lake either fixed or random depending on the science question — would be very defensible.

Exercise 5 — ML vs REML for fixed-effect comparison

(a) Why REML = FALSE? We are comparing models that differ in fixed effects. REML log-likelihoods are computed on data with the fixed effects partialed out — they are not comparable across different fixed-effect specifications. Use ML.

(b) LRT result.

anova(fit_no_x, fit_with_x)
# Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
# 4 ... ... ... ...
# 5 ... ... ... ... ~6-8 1 ~0.005-0.01

$\chi^2 \approx 6$–$8$ on 1 df, $p \approx 0.005$–$0.01$. Body length is supported as a predictor of trophic position even after accounting for lake and species.

(c) Refitting with REML for reporting. Variance components estimated by ML are biased downward (they ignore the df consumed by the fixed effects). REML corrects that bias. The fixed-effect point estimates barely change, but the random-effect SDs and the SEs that depend on them are more honest under REML.

Exercise 6 — Mixed-model diagnostics

(a) What each plot checks.

(b) Why Q-Q the random effects? The model explicitly assumes the random effects are Normal. The raw residuals are a mixture of within-group residuals and random-effect deviations; their distribution is harder to interpret. Q-Q-plotting the BLUPs directly checks the assumption that was made.

(c) If residual variance scales with fitted values. Options in rough order:

1. Try a log-transform of the response (if it's positive).
2. Fit a Gamma or Lognormal GLMM (glmmTMB, lme4::glmer with appropriate family).
3. Allow a variance function: nlme::lme(weights = varPower()) models residual SD as a power of fitted value.
4. If the response is a count or proportion, use a distribution with a built-in mean-variance relationship (Poisson/NB or Binomial/Beta-Binomial).

Exercise 7 — The "$\Delta\text{AIC} < 2$" misconception

(a) Single-run experience. Running the code a few times shows the model with the noise predictor is "competitive" every time — $\Delta\text{AIC}$ always lands below 2, usually somewhere in 0.5–2.0, and every so often it goes negative and the junk model wins outright. It never exceeds 2. Students who report seeing a value above 2 have a bug. If a student expects a large noise coefficient to raise $\Delta\text{AIC}$, correct that now: a larger $|\hat\beta_{x2}|$ means a larger likelihood-ratio statistic, which pushes $\Delta\text{AIC}$ down, not up.

(b) Formal simulation. Record the actual $\Delta\text{AIC}$, not just the indicator — the indicator hides the point.

set.seed(20260516); reps <- 1000
dAIC <- replicate(reps, {
 n <- 1000; x <- runif(n); x2 <- rnorm(n)
 y <- 2 + 3*x + rnorm(n, 0, 0.2)
 AIC(lm(y ~ x + x2)) - AIC(lm(y ~ x))
})
mean(dAIC < 2)   # 1
max(dAIC)        # 2 (as printed; strictly below — 2 - max = 7.0e-08)
mean(dAIC)       # 1.048022
mean(dAIC < 0)   # 0.141
quantile(dAIC, c(0, .25, .5, .75, 1))
#      0%      25%      50%      75%     100%
# -9.1541   0.7405   1.5848   1.9046   2.0000

The fraction is not "high," it is exactly 1.000 — 1000 out of 1000 — and that is not a property of this seed. Reason: adding one parameter costs a penalty of exactly $2\times 1$ and buys back the likelihood-ratio statistic for that parameter, so

and $\chi^2_1 \ge 0$ always. So $\Delta\text{AIC} < 2$ with probability 1. No simulation was ever needed; the algebra settles it. The simulation's only job is to make the ceiling visible: the maximum over 1000 reps sits $7\times10^{-8}$ below 2, exactly as $2 - \chi^2_1$ requires. The mean is $2 - E[\chi^2_1] = 2 - 1 = 1$, matching the observed 1.048.

The second number matters just as much: $\Delta\text{AIC} < 0$ in 14.1% of reps, i.e. the pure-noise model has the lower AIC and wins outright. That is $P(\chi^2_1 > 2) = 15.7\%$ in theory, recovered by simulation. A junk variable takes the top spot roughly one time in seven.

Exercise 8 — Type I error simulation

(a) Prediction. AIC will be the worst offender — its penalty (2 per parameter) is small enough that noise predictors easily slip into the top model.

(b) Numbers. Approximate results across 500 reps:

(From the lab's set.seed(126), 500 reps, $n = 50$ — matches the table in lab Exercise 8.) AIC and AICc keep a noise predictor in the top model about 40–45% of the time. BIC and the naive $p$-value column are both far lower and, on this seed, exactly tied at 0.162 — BIC because its $\ln n = \ln 50 \approx 3.9$ penalty is nearly double AIC's fixed 2, and $p$ because it is not doing model selection at all, just testing each coefficient once (family-wise rate $1 - 0.95^3 \approx 0.14$). BIC is not "the most conservative" here — it is about as strict as the $p$-value screen, and the real outlier is AIC's penalty of 2 being too small to keep junk out.

(c) Why BIC outperforms here. BIC's penalty grows with $n$, so for $n = 50$ it's ≈ 3.9 per parameter — almost twice AIC's. When effects are truly zero, this extra penalty drops noise predictors. Backfires when: real effects are weak. BIC will miss them (next exercise).

Exercise 9 — Type II error simulation (weak true effects)

(a) Average number of missed real effects (out of 3):

(From the lab's set.seed(2026), 500 reps, $n = 50$ — matches the table in lab Exercise 9.) BIC misses about 1.46 of the 3 true predictors on average — essentially tied with the p-value approach (1.51), just as it was in Exercise 8. AIC and AICc are the inclusive ones (missing ~1.11–1.19), which is what you want when the science is about small effects. The pattern across both exercises is the same: BIC and the $p$-value screen behave alike at this $n$, and AIC buys its extra power over both by paying in Type-I error.

(b) Implication for BIC. If your science cares about weak effects (subtle population trends, small treatment differences), BIC will under-detect them. Pick AIC/AICc — or, better, don't reduce the question to a single criterion. Look at the coefficient and its CI.

(c) Sample-size sweep. Rerunning with $n = 500$ should show all criteria detecting the strong effect (x1, true coefficient 1) reliably; AIC/AICc start picking up x2 (coefficient -0.1); only at much larger $n$ does anything reliably detect x3 (coefficient 0.01). BIC needs the largest $n$ before retaining each effect — its sample-size-dependent penalty trades off precisely against detection of weak effects.

Exercise 10 — LOOCV vs AIC for noise data

(a) Single seed (set.seed(126), n=50, y unrelated to x1, x2, x3).

cv <- sapply(forms, cvmse)
aic <- sapply(forms, function(f) AIC(glm(f, data=d)))
round(data.frame(LOOCV=cv, AIC=aic, dAIC=aic-min(aic))[order(cv), ], 3)

     LOOCV     AIC  dAIC
m3   5.995 233.059 0.000
m0   6.031 233.717 0.658
m13  6.052 233.392 0.333
m1   6.078 233.948 0.889
m23  6.205 235.059 2.000
m2   6.239 235.715 2.655
m123 6.286 235.389 2.330
m12  6.309 235.939 2.880

Both criteria pick m3 — a pure-noise predictor. They do not disagree; they fail together, and their orderings correlate at 0.95 (Spearman). Students expecting CV to rescue them from AIC's mistake will find that CV made the same mistake.

Note the scale: the intercept-only CV MSE is 6.031, not the $\sigma^2 = 9$ you might predict from rnorm(n, 5, 3). This particular sample drew var(y) = 5.910, well below the population value — a useful incidental reminder that at $n=50$ the realized variance wanders.

(b) is where the real lesson lives. m3 is m0 plus one junk predictor, so the identity from the core-concepts section pins the margin exactly: $\Delta\text{AIC} = 2 - \chi^2_1 = 2 - 2.658 = -0.658$. The whole "win" is the likelihood-ratio statistic for x3 ($\chi^2_1 = 2.658$, $p = 0.103$) happening to exceed 2 — the $P(\chi^2_1 > 2) = 15.7\%$ event already dissected in the $\Delta\text{AIC} < 2$ discussion. x3's coefficient is $-1.914$, 95% CI $[-4.232, 0.403]$, spanning zero.

Do not accept an answer that reads this table via "models within 2 AIC are competitive" — the lab has already shown that a nested junk model is within 2 by construction.

(c) 100-replication summary.

set.seed(1)
res <- replicate(100, {
 y <- rnorm(50, 5, 3); x1 <- runif(50); x2 <- runif(50); x3 <- runif(50)
 d <- data.frame(y, x1, x2, x3)
 forms <- list(m0 = y~1, m1=y~x1, m2=y~x2, m3=y~x3,
 m12=y~x1+x2, m13=y~x1+x3, m23=y~x2+x3, m123=y~x1+x2+x3)
 cv <- sapply(forms, function(f) cv.glm(d, glm(f, data=d))$delta[1])
 aic <- sapply(forms, function(f) AIC(glm(f, data=d)))
 c(cv_picks_truth = unname(which.min(cv) == 1),
 aic_picks_truth = unname(which.min(aic) == 1),
 both_agree = unname(which.min(cv) == which.min(aic)))
})
round(rowMeans(res), 2)

 cv_picks_truth aic_picks_truth      both_agree
           0.58            0.57            0.92

LOOCV 0.58, AIC 0.57. Over 2000 replications these settle to 58.4% and 57.6% — a gap of 0.8 points against a Monte Carlo error of ±2.2, i.e. no detectable difference. The two criteria pick the same model in over 90% of datasets, and each is fooled by noise about 42% of the time.

(The unname() calls only stop which.min's labels leaking into the row names; without them the output reads cv_picks_truth.m1, which looks like a result but is just the first replicate's winner.)

(d) Is this a defect in the criteria? No — and no better criterion fixes it. For a Gaussian linear model, leave-one-out CV and AIC are asymptotically equivalent (Stone 1977): they are the same criterion in the limit, and by $n=50$ they have essentially converged. Holding data out is not an independent safeguard against AIC's errors; it is a longer road to the same number. Checked from $n=30$ to $n=1000$, the two track each other within half a percentage point at every sample size.

With eight candidates and $n=50$, three noise predictors sometimes produce spurious in-sample fit, and any criterion that ranks models by fit to a single dataset will sometimes rank noise first. The criteria answer the question asked — "which of these eight predicts best?" — and return a winner even when the honest answer is "none, and they are indistinguishable."

Exercise 11 — LRT vs AICc for a random-slope decision

(a) Why refit = FALSE? Both models share the same fixed-effect specification but differ in random-effect structure. REML log-likelihoods are comparable in that setting and give less-biased variance components. refit = FALSE tells anova not to refit them under ML.

(b) Disagreement scenario. AICc and the REML LRT often agree, but not always. Each statistic answers a different question:

• LRT: "Is the larger model significantly better at maximizing likelihood (given the asymptotic null distribution)?"
• AICc: "Is the larger model expected to predict better, after a parsimony penalty?"

For an inference question about whether the random slope matters biologically, the LRT plus the CI on the slope variance is the cleanest evidence. For a prediction question, prefer AICc (and ideally cross-validation).

(c) Boundary-test conservatism. Testing whether a variance component is zero puts the null on the boundary of the parameter space (variances cannot be negative), so the classical $\chi^2$ null overstates uncertainty. The correction depends on how many parameters you added. Adding one variance component with no correlation term: the null is $\tfrac{1}{2}\chi^2_0 + \tfrac{1}{2}\chi^2_1$, and halving the reported $p$-value is right. Adding a random slope: you add a variance and a correlation, anova() reports Df = 2, and the null is $\tfrac{1}{2}\chi^2_1 + \tfrac{1}{2}\chi^2_2$ — halving is then wrong. Worked on the Exercise 2 / P2 model ($\chi^2 = 0.078$, Df 2): R gives $p = 0.962$, halving gives $0.481$, correct mixture gives $0.871$. Grading note: a student who reaches for "divide by 2" on the random-slope test has learned the rule without its condition. Practical fix: RLRsim::exactRLRT() or a parametric bootstrap.