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$. Because both columns are logged, this slope is the allometric exponent: weight scales as about the cube of length, $\text{weight} \propto \text{length}^{3.02}$, which is what you expect when a fish grows without changing shape, since mass tracks volume and volume tracks length cubed. A slope near 3 is the biologically meaningful reading, not "3 units of weight per unit of length". CIs from intervals(m1) exclude zero comfortably, though for an allometric exponent the interesting comparison is against 3, not against 0.

(c) Among-subject SD on the intercept is around 0.211, on the log-weight scale. A difference on a log scale is a ratio once you exponentiate: $e^{0.211} = 1.23$, so at any given length a one-SD subject is about 23% heavier than the population average, and a two-SD subject about 53% heavier. That is the interpretation to report; "±0.21" alone is not in units anyone can picture.

(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: $0.1060$ against $5.4168$ on this file, because the random effect on $L_\infty$ absorbed the among-subject variability that nls was forced to dump into residuals.

which is the nls residual variance of $29.34$. The additive partition is a random-intercept result. It holds when the random effect enters the mean with coefficient 1. $L_\infty$ does not, so you have to carry the attenuation factor through.

(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. In general, subjects with few or noisy age-points get pulled hard toward the population mean while data-rich subjects stay near their own MLE, but do not go looking for that contrast in fish_growth.csv: all 10 fish have exactly 10 age-points each, so the design is balanced and every subject receives the same amount of pull. Seeing sparse subjects shrink harder than dense ones requires an unbalanced dataset.

(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, and AIC agrees: 300.8 without the slope against 77.0 with it, a drop of 224 units. Both tools point the same way, which is the ordinary case.

(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.

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, so 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 A1. 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.

(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 A3 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.292,\ 0.463]$, 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     6  41.88  -14.94  42.59  1   6.7e-11 ***

AICcmodavg::aictab(list(intercepts_only=m_int, with_species_slope=m_slop))
#                    K   AICc  Delta_AICc  AICcWt
# with_species_slope 6  42.36    0.00      1.00
# intercepts_only    5  82.81   40.45      0.00

The LRT gives $\chi^2 = 42.6$ on 1 df, $p \approx 7\times10^{-11}$; AICc favors the random-slope model by $\Delta$AICc $= 40.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 A3, 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 7\times10^{-11}$) 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.