Mixed effects models in the morning, model selection in the afternoon

(a) Let $y_{ij}$ be density measurement $i$ in stream $j$. lm1: $y_{ij} = \beta_0 + \varepsilon_{ij}$ — two parameters ($\beta_0, \sigma^2$). lm2: $y_{ij} = \alpha_j + \varepsilon_{ij}$ with a free mean $\alpha_j$ for each of the $j=1,\dots,100$ streams — 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)$ — only 3 parameters ($\beta_0, \sigma_\text{stream}^2, \sigma^2$), regardless of how many streams there are.

(b) Residuals from lm1 cluster strongly by stream — entire streams are shifted above or below zero. This violates the independence assumption of OLS. Standard errors from lm1 are too small (anti-conservative).

(c) Stream 11 has the lowest mean (8.9) and Stream 43 the highest (111.7); every stream has exactly 3 observations. The partially-pooled estimates from lm3 sit slightly closer to the grand mean of 63.26 — Stream 11 moves 8.9 → 9.7, Stream 43 moves 111.7 → 111.0.

What shrinkage does here is simple: every stream is pulled toward the grand mean by the same percentage of its distance from that mean — about 1.5%. Because it is the same fraction for all of them, a stream that starts far from the grand mean moves more in absolute terms — Stream 11 sits about 54 units below the mean and is nudged ~0.8, Stream 43 about 48 units above and is nudged ~0.7 — but no stream is shrunk by a larger fraction than any other. Each estimate keeps roughly 98.5% of its own stream average and borrows just ~1.5% from the overall mean.

Two things set how big that pull is: how noisy each stream's own average is, and how different the streams genuinely are. Here the streams differ a lot relative to the noise in a 3-observation average, so the model mostly trusts each stream's own data and barely shrinks — and because every stream has the same number of observations (3), every stream gets the same small pull. Shrinkage grows when a group has few or noisy observations, or when the groups are nearly alike; it fades toward nothing when each group's own data are already precise.

(d) Approximate output:

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

Fixed effects: Density ~ 1
 Value Std.Error DF t-value
(Intercept) 63.26 2.16 200 29.3

intervals(lm3)$reStruct$Stream
 lower est. upper
sd((Intercept)) 18.60 21.43 24.68
intervals(lm3)$sigma
 lower est. upper
 4.14 4.56 5.03

ICC = $21.43^2 / (21.43^2 + 4.56^2) \approx 459 / 480 \approx 0.957$. Almost all variation in fish density is between streams; within a stream, the three replicate measurements are very similar. That is exactly the situation in which ignoring the random effect would give you wildly over-confident inference.

(a) $W_{ij} = \beta_0 + b_{0,j} + \beta_1 L_{ij} + \varepsilon_{ij}$, with $b_{0,j}\sim\mathcal{N}(0,\sigma_S^2)$ and $\varepsilon_{ij}\sim\mathcal{N}(0,\sigma^2)$. $\beta_1$ is the population-average slope of weight on length, shared across subjects. $b_{0,j}$ is the subject's intercept deviation from the population mean.

(b) Approximate output: $\hat\beta_0 \approx -0.03$, $\hat\beta_1 \approx 3.02$ — i.e., a one-unit increase in the length variable corresponds to about 3.0 units of weight on average across subjects. This file's Length and Weight columns are on a transformed/index scale (ranges $\approx 0.7$–3.0 and 1.4–9.4), not raw centimetres and grams, so read the slope in those units rather than attaching a gram figure to it. CIs from intervals(m1) exclude zero comfortably.

(c) Among-subject SD on the intercept is around 0.211; so subjects of the same length differ in expected weight by roughly $\pm 0.21$ (1 SD) or $\pm 0.42$ (2 SD) on the weight scale, even after accounting for the length effect.

(d) Within-subject residual SD around 0.21. ICC = $0.211^2 / (0.211^2 + 0.206^2) \approx 0.045 / 0.087 \approx 0.51$. Even after the length covariate, about half the remaining variation is between subjects: individual fish are persistently heavier or lighter than the population line predicts, so repeated measurements on one fish are far from independent. That is exactly why the random intercept is needed here.

(e) The REML LRT (use refit = FALSE because we are comparing random-effect structures, not fixed effects, and REML likelihoods are valid for that comparison). With this dataset, the random-slope model usually shows little improvement — the LRT $p$-value is large, and the slope variance is near zero (singular fit warning). Conclusion: a random intercept suffices.

(a) Asymptotic body size depends on genetics, food availability, and individual life history. Pooling all fish into a single $L_\infty$ ignores real biological heterogeneity.

(b) Residuals from the fixed-effects nls form coherent vertical clusters by subject — entire subjects are systematically above or below zero. This is the same diagnostic as in the streams example: independence is violated.

(c) The mixed-model within-subject residual SE is much smaller than the nls residual SE, because the random effect on $L_\infty$ absorbed the among-subject variability that nls was forced to dump into residuals. Total variance is partitioned: $\sigma^2_\text{total} \approx \sigma_{L_\infty}^2 + \sigma^2_\text{within}$.

(d) It returns the BLUP (best linear unbiased predictor) of each subject's $L_\infty$: population mean plus the partially-pooled random deviation. Compared with running a separate nls per subject, BLUPs use information from other subjects to stabilize estimates for subjects with sparse or noisy data. Subjects with only a few age-points get pulled hard toward the population mean; subjects with many measurements stay near their own MLE.

(a) Both are grouping variables along which observations are nested (multiple fish per species, multiple fish per lake). If you'd be happy to generalize the length–TP relationship to other lakes or other species, treating them as random captures that intent. The counter-argument: with only 3 species and 6 lakes, you have very few "draws" from each random-effects population — the among-group variance is hard to estimate, and a fixed-effects model with explicit species/lake terms (and possibly an interaction) is more defensible.

(b) The model assumes:

i.e., each species has a (intercept, slope) pair drawn from a bivariate Normal centered at the fixed effects, with its own SDs and a correlation $\rho$. Species with high intercepts may also have steep slopes (or shallow, depending on $\rho$).

(c) Best AICc is typically M8 (random intercept & slope on species, random intercept on lake). Approximate variance components:

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

Most slope variation lives across species — different feeding ecologies. The answer to "or across lakes?" is in what M8 does not contain: the winning model has a random slope on species but only a random intercept on lake, so there is no lake slope term in the VarCorr block to interpret. That absence is the result, not an oversight. Two things confirm it. First, the ranking: M2 is M8 plus a lake slope, and it is 3.31 AICc worse — the extra slope does not pay for its parameters. Second, if you fit M2 anyway and look, its lake slope SD is 0.024, against a species slope SD of 0.157 — nearly an order of magnitude smaller, and next to the lake intercept SD of 0.45. So lake shifts the trophic-position baseline up or down but barely bends the length–TP relationship, while species changes both. Within-fish residual (0.22) is small relative to species variation. Note the perfect-correlation warning: the slope–intercept correlation comes back as +1.00, a classic singular fit. 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. With only 3 species this is not surprising. Treat it as a warning that the random-slope structure is over-specified for this dataset, not as a real biological correlation.

Caveat on the AICc table: aictab() warns that ranking models with different fixed effects needs ML, not REML. Here the fixed part (Z_Length) is identical across M1–M8 and only the random structure differs, so the REML comparison is exactly the case REML is for — the warning can be noted and set aside.

(d) With few levels, the variance components are estimated with so little information that singular fits and unstable correlations are common — and the inferential machinery (Wald CIs, LRTs) is poorly calibrated. Switch to fixed effects for species when (i) you have ≤ 3–4 levels, (ii) you care about each level by name (S1 vs S2 vs S3 is biologically meaningful), or (iii) you cannot defend a "population of species" sampling story.

(a) We're comparing models that differ in their fixed-effect structure. REML log-likelihoods are not comparable across different fixed-effect specifications. anova() on two REML lmer fits will (by default) refit them with ML before computing the LRT, but it's better practice to fit them with REML = FALSE from the start so you can see what you're comparing.

(b) The LRT statistic is twice the gain in maximized log-likelihood from adding the one fixed slope: $\Lambda = 2(\ell_\text{alt} - \ell_\text{null})$. Here the ML log-likelihoods are $\ell_\text{null} = -146.4$ and $\ell_\text{alt} = -33.5$, so

Body length is overwhelmingly supported — this is not a marginal effect. You can sanity-check the magnitude against the Wald view: the slope is $\hat\beta = 0.420$ with SE $0.019$, so $t = 22$ and $t^2 \approx 467$; the LRT and Wald statistics are the same test asymptotically and both are astronomically far into the tail (they differ here only because $t^2$ uses the curvature at the estimate while $\Lambda$ uses the actual likelihood drop). With $n = 180$ observations spanning the length range, a standardized-length slope of 0.42 SD of trophic position per 1 SD of length is a strong, well-estimated signal, not a whisper.

(c) REML gives less-biased variance components. The point estimates of the fixed effects barely change, but the standard errors are more honest under REML.

(a) Top-left (fitted vs Pearson residual): looking for randomness around zero with constant spread — funnel shape means heteroscedasticity. Top-right (covariate vs residual): looking for any remaining structure — curvature implies missing nonlinearity. Bottom-left (residual by species): all boxes centered on zero with similar spread — a species sitting away from zero means the random intercept didn't fully absorb the species effect. Bottom-right (Q-Q of Lake random intercepts): the assumption is that lake intercepts are Normal; serious departures from the line are a warning.

(b) Random effects have an explicit Normal assumption in their definition. Raw residuals are also assumed Normal, but their distribution is contaminated by the random-effect distribution and Pearson scaling — the Q-Q on random effects is a more direct check of the assumption.

(c) Options: log-transform the response; switch to a Gamma or Lognormal GLMM with a log link; allow a variance function in nlme::lme(weights = varPower()); or use a model with a built-in mean-variance relationship (Poisson/NB for counts, Beta for proportions).

The nonsense-variable model is within 2 AIC units of the correct model essentially 100% of the time — it always is. This is not luck; it falls straight out of the definition of AIC. Write $\text{AIC} = -2\ell + 2k$ for each model and subtract the simpler ($k$ parameters) from the more complex ($k+1$ parameters):

because $2(\ell_{\text{big}} - \ell_{\text{small}})$ is exactly the likelihood-ratio statistic $\chi^2_1$ for the extra term (twice the gain in maximized log-likelihood). This is the general identity worth keeping — adding one parameter moves AIC by $2 - \chi^2_1$: the penalty for a parameter is a flat $+2$, while the reward is a random $\chi^2_1$. Two consequences follow at once. Since $\chi^2_1 \ge 0$ always, $\Delta\text{AIC} \le 2$ always — a nested one-extra-parameter model is within 2 by construction, whatever the data say. And the junk model wins outright ($\Delta\text{AIC} < 0$) exactly when $\chi^2_1 > 2$, which for a true null happens with probability $P(\chi^2_1 > 2) = 0.157$ — about 16% of the time. The $\Delta\text{AIC} < 2$ rule does not screen out useless variables at all.

The lesson: the $\Delta\text{AIC} < 2$ rule is structurally biased toward including junk. Anderson & Burnham's own writing has been used selectively here; many authors (Arnold 2010 is the standard reference) push back on the "include any model within 2" practice and emphasize that the included variable's coefficient and CI have to do real work.

Approximate results across 500 reps (your seed will differ):

AIC and AICc retain a noise predictor in the top model about 40–45% of the time. You can predict that rate rather than only reporting it, using the $\Delta\text{AIC} = 2 - \chi^2_1$ identity from Exercise 7: a single noise predictor enters the AIC-best model exactly when it clears $\chi^2_1 > 2$ — probability $0.157$ each. With three near-independent noise predictors, the chance that at least one slips in is

in line with the simulated 40–45% (the small excess is joint selection and the mild correlation among the fitted predictors). The same reusable move explains the two conservative criteria — each just raises the bar the extra variable must clear. The family-wise p-value rate is $1 - 0.95^3 \approx 0.14$. BIC replaces AIC's flat penalty of 2 with $\ln n \approx 3.9$ for $n = 50$, so a junk predictor now needs $\chi^2_1 > 3.9$ (probability $P(\chi^2_1 > 3.9) \approx 0.048$ each, hence $1 - (1 - 0.048)^3 \approx 0.14$ overall) — about as strict as the p-value screen, and matching the observed $\approx 0.16$. AIC's penalty of 2 is simply too small to keep junk out.

When would BIC backfire? When true effects are weak. BIC's heavier penalty buys Type-I safety at the cost of Type-II error — see the next exercise.

Approximate per-rep averages (out of 3 true effects):

BIC misses about 1.5 of the 3 true predictors on average — essentially tied with the p-value approach. AIC and AICc are noticeably more inclusive of weak effects (missing ~1.1–1.2), which is what you want when the science is about small effects. Lesson: there is no free lunch. Whatever Type-I behavior you buy with a heavier penalty, you pay back in Type-II.

(c) Bigger n. Re-run both harnesses — this one (three real effects, Type II) and Exercise 8 (pure noise, Type I) — at $n = 50$ and again at $n = 500$. Averaged over 500 datasets, out of 3 predictors:

Read the columns. AIC's Type-I rate barely moves with sample size (0.51 → 0.47): it keeps roughly half a noise predictor no matter how much data you give it, because its penalty ($2$ per parameter) never grows. BIC's Type-I rate collapses (0.20 → 0.05): its penalty is $\log(n)$ per parameter, which climbs with $n$, so junk gets squeezed out as data accumulate. Meanwhile every Type-II rate falls as $n$ grows — even the tiny $0.01\,x_3$ effect starts to clear the bar once there is enough data to see it. The Type-I/Type-II trade-off does not go away; it just slides: BIC keeps buying cleaner false-positive control at the price of the weakest real effects, and AIC keeps the weak effects at the price of some persistent false positives.

(a) Both criteria pick m3 — a pure-noise predictor. They agree, and they are both wrong:

     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

The two orderings are nearly identical (Spearman correlation 0.95). This is the first hint that CV is not an independent check on AIC.

(b) m3 beats the truth by $\Delta$AIC = 0.658 and by 0.036 in CV MSE — six tenths of one percent. But the interesting part is why the margin has that value. m3 is m0 with one extra junk predictor, so the identity from the core-concepts section applies exactly:

The entire margin is just the likelihood-ratio statistic for x3, which is $\chi^2_1 = 2.658$ ($p = 0.103$ — not close to significant). x3 wins for one reason and one reason only: its $\chi^2$ landed above 2. That is exactly the event you already met when the $\Delta\text{AIC} < 2$ rule was picked apart — $P(\chi^2_1 > 2) = 15.7\%$, the rate at which a junk variable takes the top spot outright. Exercise 10 is not a new phenomenon; it is that same 16% event, caught in the act.

So do not read this table by asking which models fall within 2 AIC of the winner. As you saw earlier, a nested junk model is within 2 by construction, so that question has a guaranteed answer and tells you nothing. Ask instead what x3's coefficient does: $-1.914$ with a 95% CI of $[-4.232,\ 0.403]$, spanning zero. The variable that "won" cannot even sign its own effect.

(c) Over 2000 simulated datasets the intercept-only model wins 58.4% of the time on LOOCV and 57.6% on AIC. That gap of 0.8 points is smaller than the Monte Carlo error on either estimate (±2.2), so the honest reading is that the two criteria perform identically. They select the same model 90% of the time. Both are fooled by noise in roughly 42% of datasets.

This is not a coincidence of this simulation. For a Gaussian linear model, leave-one-out CV and AIC are asymptotically equivalent — Stone (1977) proved they are the same criterion in the limit, and at $n=50$ they have already nearly converged. Holding data out does not buy you protection AIC lacks; it is a different route to the same number.

That equivalence is the best reason AIC is everywhere. LOOCV is the gold-standard estimate of out-of-sample prediction error, but it costs you $n$ refits; AIC gives you the same answer in closed form from a single fit. So when the goal is genuinely prediction — which model will forecast best on data it has not seen — AIC is doing exactly its job, and doing it cheaply, which is why it became the standard tool. The trouble in this exercise is not that AIC predicts badly; it is that we asked a prediction criterion an inference question — "which variable is real?" — that it was never built to answer. Much of the ecology literature slips here, reading a low-$\Delta$AIC predictor as an established effect. Keep the two jobs apart: AIC (≈ LOOCV) ranks models for prediction; coefficients, intervals, and study design tell you whether an effect is real.

(d) Not a defect, and no. With eight candidate models and $n=50$, three noise predictors will produce spurious in-sample fit some of the time, and any criterion that ranks models on fit to this one dataset will sometimes rank noise first. The criteria are answering the question they were asked — "which of these eight predicts best?" — and dutifully returning an answer even though the honest answer is "none of them, and they are all the same." Selection returns a winner whether or not one deserves to win.

What helps is not a cleverer criterion but a different question. Once a table is in front of you, the rank is the least informative thing in it — make the winning variable justify itself with its coefficient and interval, as in (b), where x3's CI comfortably contains zero (this is Arnold's 2010 point). And upstream of that, keep the candidate set small and tied to hypotheses you held before seeing the data: the 42% failure rate is a direct consequence of dredging eight models out of pure noise. Fewer, better-motivated models is the only move here that changes the failure rate rather than relabeling it.

Caveat: both criteria are about prediction. If your question is "is variable X causally important," neither AIC nor CV answers it — you need a causal frame and (often) a designed experiment.

(a) The two models share fixed effects but differ in random structure. We compare REML log-likelihoods because they are unbiased for variance components. Setting refit = FALSE prevents anova from refitting them under ML.

(b) Run them. Here the two agree emphatically — the species-level random slope is strongly supported both ways:

anova(m_int, m_slop, refit=FALSE)
#         npar   AIC    logLik  Chisq Df  Pr(>Chisq)
# m_int      5  82.47  -36.23
# m_slop     7  35.68  -10.84  50.79  2   9.4e-12 ***

AICcmodavg::aictab(list(intercepts_only=m_int, with_species_slope=m_slop))
#                    K   AICc  Delta_AICc  AICcWt
# with_species_slope 7  36.33    0.00      1.00
# intercepts_only    5  82.81   46.48      0.00

The LRT gives $\chi^2 = 50.8$ on 2 df, $p \approx 9\times10^{-12}$; AICc favors the random-slope model by $\Delta$AICc $= 46.5$ (Akaike weight $\approx 1.00$). This is easy territory — keep the slope. But notice what each statistic claimed, because they answer different questions and will not always agree. The LRT asks "is the larger model significantly better at maximizing likelihood?" — an inference question about whether the slope variance is really nonzero. AICc asks "is the larger model expected to predict better out of sample?" When they disagree (LRT says "not significant" but AICc still prefers the bigger model), let the question decide: for inference — does the random slope matter biologically? — lean on the LRT and the variance-component CI; for prediction, trust AICc (which, from Exercise 10, is ≈ LOOCV).

(c) One catch makes the LRT $p$-value untrustworthy in borderline cases. The null hypothesis is that the slope variance is exactly zero — but a variance cannot be negative, so zero sits on the very edge of the values the parameter is allowed to take. When you test a parameter at the boundary of its space, the usual $\chi^2$ reference distribution is the wrong one, and it makes the reported $p$-value too large (conservative). A genuinely important random slope can therefore come back looking "not significant."

Here it does not matter — the effect is so strong ($p \approx 9\times10^{-12}$) that no correction would change the verdict. But in a borderline case the boundary correction can flip your conclusion. The practical rule: for a boundary test — "is this variance zero?" — do not trust anova()'s $p$-value. Use RLRsim::exactRLRT() or a parametric bootstrap, which build the correct reference distribution by simulation and give you an honest $p$-value.