What we did today, in one paragraph
We tackled the fact that real ecological data are not independent. Fish are inside lakes, plots are inside sites, measurements are repeated on the same individual. Mixed effects models let us include these grouping structures as random effects: rather than estimating a separate parameter for each group, we estimate how much groups vary and let each group's estimate be pulled toward the population mean. In the afternoon, we faced the harder question. When we have several plausible models, how do we pick one? We compared the likelihood ratio test, AIC, AICc, BIC, and cross-validation, and ran simulations to see what each criterion does, and doesn't, protect us from. The honest answer at the end of the day is that there is no automated procedure that replaces thinking carefully about your science question.
The big-picture story
Every modeling decision so far this week has been about getting the data-generating story right. Day 1 gave us distributions to tell the story for one observation. Day 2 wired predictors into the mean via a link function. Day 3 adds the next layer: observations are organized into groups that share something we haven't measured. We need to model that shared-but-unmeasured component, or our standard errors will be too small and our discoveries will mostly be false.
The one sentence that unifies the day: independence is the assumption you should distrust by default, and mixed models are how you replace it with something more honest.
Fixed effects vs random effects, in English
The scenario we'll use (and yes, "badjuice" is made up). A new graduate student joins your lab to study razor clam abundance along 10 beaches on the Oregon coast. At each beach they record clam abundance and the local concentration of badjuice, a fictitious environmental toxin invented for this course (it is not a real chemical), and they repeat the survey for four years. The student wants to know: does badjuice affect clam abundance?
That one scenario already contains the whole distinction below: badjuice is the thing you want a number for, and beach is the grouping you must respect but don't care about one-by-one.
The technical distinction is often muddled. Here is the most useful way to think about it:
- A fixed effect is something whose specific value you want to know. "What is the effect of badjuice on clams?" Badjuice is a fixed effect because the slope is the answer to your question.
- A random effect is a grouping variable. Each level (beach 1, beach 2, ..., beach 10) is one draw from a larger pool of possible levels. You don't particularly need to know beach 7's specific intercept; you need to make sure that putting all 10 beaches' observations together didn't pretend they were independent when they weren't.
A useful test: ask yourself "if next year I had a brand new beach, would I want my model to be able to say anything useful about it?" If yes, beach is a random effect, and you are learning about beaches as a population. If you only care about these particular 10 beaches in this particular year, you could just use fixed effects for them.
A few rules of thumb that catch most decisions:
- Fewer than ~5 levels? Probably use a fixed effect, because the random effect variance is essentially unidentifiable.
- Treatment, sex, species (when there are 2 or 3)? Fixed.
- Stream, plot, site, year (when there are many)? Random.
- Repeated measurements on the same animal? Animal is a random effect (or you have a multilevel model with animal nested in higher groups).
What "shrinkage" actually does (the magic of partial pooling)
Imagine measuring fish density in 100 streams, three measurements in each, the Day 3 stream dataset. A few streams come out looking extreme (Stream 11 has the lowest average, 8.9; Stream 43 the highest, 111.7, against a grand mean of 63.3). If you fit a separate intercept per stream (no pooling), those streams keep their raw averages exactly as measured. If you fit one grand mean and ignore stream (complete pooling), you miss the fact that streams really are different.
A mixed model splits the difference. It says: "the streams probably aren't quite as different as their raw means suggest, so let me pull the extreme ones a little toward the overall average." That pull is shrinkage. The strength of the pull depends on:
- How much data each stream has. A stream with 30 observations gets pulled very little; its own data dominate. A stream with 3 observations gets pulled more; we trust its own average less.
- How variable streams are overall. If among-stream variance is huge, even small samples convince us streams really do differ. If it's small, the model assumes most variation is residual noise and pulls hard.
Both dials are pinned in the Day 3 stream data. Every stream has exactly the same three measurements, so the first dial cannot differ between streams. And streams differ enormously (among-stream SD 21.4 against a within-stream SD of 4.6), so the second dial is set to "barely pull at all". The result: every one of the 100 streams is pulled by the same 1.49% of its distance from the grand mean, keeping 98.5% of its own average. Stream 11 moves from 8.9 to 9.7; Stream 43 from 111.7 to 111.0. Nothing here shrinks harder than anything else. To see that you would need an unbalanced design, where some streams have far fewer measurements than others. The interactive shrinkage demo lets you move both dials yourself and watch strong shrinkage appear.
The shrinkage formula is not magic. It is exactly what an idealized Bayesian would do with a prior "fish density is around the population mean, with variance $\sigma_\text{stream}^2$" combined with the actual data for each stream. Mixed models are quietly Bayesian.
Why this is good science. Shrinkage gives more honest predictions for new streams. A separate-intercept (no-pooling) model overfits, predicting extreme values for new streams that won't be there. The complete-pooling model underfits, ignoring real variation. Partial pooling is the predictive sweet spot, and the principle is general (anywhere you have groups of related estimates, partial pooling tends to outperform either extreme).
Model selection in English
You'll often have several plausible models. The question "which one is best?" is actually three different questions in disguise:
| The question you're really asking | The right tool |
|---|
| Does adding this predictor explain significantly more variation than chance? | Likelihood ratio test (for nested models) |
| Of these models, which would be expected to predict best on new data? | AIC / AICc (in theory); cross-validation (in practice) |
| Of these models, which is most plausible given a (heavier) penalty for complexity? | BIC |
| Which predictor is causally important? | None of the above; this needs a study design |
The trap that has wrecked many graduate theses is to run all subsets of predictors through dredge() or stepwise selection, then write up the "best" model as if it were a normal regression. The coefficients in that best model are biased away from zero, the standard errors are too small, and the p-values are meaningless. This is true even when you use AIC instead of p-values. The simulations in today's lab show that AIC keeps junk predictors in the top model surprisingly often.
What AIC actually is
AIC, the Akaike Information Criterion, is
$\text{AIC} = -2\,\ell(\hat\theta) + 2k$
where $\ell(\hat\theta)$ is the maximized log-likelihood and $k$ is the number of parameters. Lower is better.
The most accurate plain-English description is this: AIC is an estimate of how far (in a precise statistical sense called Kullback–Leibler divergence) the model's predictive distribution is from the truth. Models with lower AIC are estimated to be closer to the truth. The 2-per-parameter penalty is what you'd derive by working out how much in-sample likelihood overstates out-of-sample predictive quality, on average, under some regularity assumptions.
What AIC is not:
- Not a p-value. There is no null hypothesis being tested.
- Not a measure of effect size or biological importance.
- Not a hard cutoff. The famous "$\Delta\text{AIC} < 2$" rule of thumb is actually quite weak. Adding a pure noise predictor costs a penalty of exactly 2 but buys back some fit, so the net change is $\Delta\text{AIC} = 2 - \chi^2_1$, always under 2 and about 1 on average. Noise variables are therefore "competitive" by construction.
- Not magic. AIC's variance is large in small samples (use AICc), and it tends to retain more complex models than necessary in noise simulations.
The corrected version, AICc, adds an extra penalty $2k(k+1)/(n-k-1)$ that matters when sample size is small relative to the number of parameters. Use AICc by default unless $n/k$ is large. BIC uses a $k \ln n$ penalty, which is larger than AIC's $2k$ as soon as $n \gtrsim 8$, so BIC is more conservative about retaining predictors, at the cost of missing weak real effects.
What cross-validation actually is
Cross-validation answers the question "how well does this model predict data it hasn't seen?" by literally hiding some data, fitting the model on what's left, and seeing how badly the predictions miss the hidden values. Repeat for every observation (leave-one-out CV) or every fold of, say, ten (10-fold CV), and average the prediction errors.
If a model is overfit, using random noise to make the training data look great, it will look bad on the held-out data. CV catches that directly.
But CV is not an independent second opinion on AIC. For a Gaussian linear model, leave-one-out CV and AIC are asymptotically the same criterion (Stone 1977), and at the sample sizes in this lab they have already nearly converged. In Exercise A3 the two rankings correlate at 0.95, they pick the same model about 90% of the time, and over 2000 simulated datasets they recover the truth at indistinguishable rates (58.4% vs 57.6%). Holding data out does not buy you protection that AIC lacks. It is a slower route to the same number. That equivalence is the reason AIC is everywhere: LOOCV costs you $n$ refits, AIC gives the same answer in closed form from one fit.
For any modeling project whose goal is "predict X for new units", cross-validation should appear in your methods. It is much more honest than reporting the in-sample $R^2$. Reporting AIC instead is not a lesser substitute (for Gaussian models AIC is the closed-form version of LOOCV), but neither one answers "is this effect real?" Both rank models for prediction; only the coefficient, its interval, and the study design speak to whether an effect exists.
Common pitfalls
"I'll just put everything as a random effect to be safe." Random effects need enough levels (≥ 5, ideally ≥ 8) to estimate the variance component. With 2–3 levels, the model is essentially unidentifiable and you'll see singular fit warnings. Use fixed effects there.
"My fixed effect went away when I added the random effect." That can happen when the variable you thought was the predictor was really just a marker for the group. (Imagine "lake pH" as a fixed effect when really all you have is one pH value per lake. Adding a lake random effect absorbs the pH information.) The fix is study design, not statistics.
Comparing REML log-likelihoods across different fixed effects. Don't. REML likelihoods only mean the same thing when the fixed-effect part is the same. For fixed-effect comparisons, refit with ML. For random-effect comparisons, REML is fine.
"$\Delta\text{AIC} = 1.7$ means model A and model B are equally good." Not really. The lab simulations show that a model containing a pure noise predictor sits within 2 AIC units of the truth essentially always, because adding one parameter can raise AIC by at most 2, so a junk variable is never ruled out by the "within 2" rule. Treat $\Delta\text{AIC} < 2$ as "the data don't strongly distinguish them," not as evidence that the extra predictor matters.
Running dredge() and reporting the resulting p-values. The selection process biases coefficients and shrinks standard errors. If you must use all-subsets selection, at minimum (i) report it as exploratory, (ii) validate on held-out data, and (iii) don't quote p-values from the selected model as if it had been pre-specified.
Singular fit warnings from lmer. Usually a sign that a random-effect variance is being estimated at or near zero, or a correlation is being estimated near $\pm 1$. Often resolved by simplifying the random structure (drop a random slope, drop a correlation), z-scoring predictors, or noting that you don't have enough levels to estimate that piece.
Looking ahead
Tomorrow (Day 4) we'll formalize what we've been quietly doing with maximum likelihood: the engine inside lm, glm, lmer, nls, and nlme. Once you see how MLE works, you'll see that almost everything in this course is the same idea applied to different distributions and link functions. We'll also start Bayesian inference, which gives a much more direct answer to many of the questions today's frequentist tools handle awkwardly, such as "what is the probability my effect is positive?", "what is the posterior on $\sigma_\text{stream}$?", and so on.
If today's mixed-model material didn't fully click, re-work the lab exercises and the problem set before tomorrow morning.