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: 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 — you're 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 — 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 18 streams. Two of those streams (call them Stream 11 and Stream 12) only had three observations each, and one of those three happened to be unusually low. If you fit a separate intercept per stream (no pooling), those two streams will look extreme. If you fit one grand mean and ignore stream (complete pooling), you'll 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 — 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; the data dominate. A stream with 3 observations gets pulled a lot; we don't trust the data as much.
- 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.
The shrinkage formula is not magic — it's 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 — 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. AIC tries to predict CV behavior from in-sample likelihood. When they agree, that's reassuring. When they disagree, trust CV for prediction questions.
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$ or the AIC of the best model.
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 — 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 — "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.