Mixed models & model selection: the graded hand-in

Reading a mixed model someone else fitted

Part (a) needs no computer. Most of the mixed models you meet in your career will arrive as printed output in a paper or a collaborator's email, and reading one is a different skill from fitting one.

A colleague weighs nestlings in 18 broods and asks whether chicks in larger broods are lighter. They send you this and nothing else:

Linear mixed model fit by REML ['lmerMod']
Formula: mass ~ broodsize + (1 | brood)

Random effects:
 Groups   Name        Variance Std.Dev.
 brood    (Intercept) 7.218    2.687
 Residual             2.249    1.500
Number of obs: 74, groups:  brood, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  22.2315     2.0053  11.087
broodsize    -0.6663     0.4560  -1.461

1. From that output alone: (i) how much of the variation in chick mass is between broods rather than within them? Compute the ICC and say it in a sentence a field biologist would use. (ii) What is the estimated effect of one extra chick in the brood, and how confident would you be that it is real? (iii) The colleague concludes "brood size doesn't matter." Give one reason that conclusion might be premature, given the design they used.

Now back to the length–weight data from the walkthrough, where you already have m1 and m1b fitted:

TheData <- read.csv(here::here("Day3_MixedModels_ModelSelection", "data", "length_weight.csv"))
LenW <- data.frame(TheData); LenW$Subject <- as.factor(LenW$Subject)

m1  <- lme(Weight ~ Length, data=LenW, random = ~ 1 | Subject, method="REML")
m1b <- lmer(Weight ~ Length + (1 | Subject), data=LenW)   # same model, lme4 syntax

2. Report the within-subject residual SD from m1, and compute the ICC, $\tau^2/(\tau^2 + \sigma^2)$, where $\tau$ is the among-subject SD and $\sigma$ is the residual SD. Then use those variance components to answer a design question.

   A new study can afford exactly 60 measurements. It could measure 6 groups 10 times each, or 30 groups twice each. Its goal is the population mean of the response, the average across all groups rather than the mean of any one group, so here “better” means estimating that population mean more precisely, with a smaller standard error. For $g$ groups each measured $n$ times, that standard error is

   $$\text{SE} = \sqrt{\dfrac{\tau^2}{g} + \dfrac{\sigma^2}{g\,n}}.$$
   1. Compute the SE for both allocations using your length–weight estimates of $\tau$ and $\sigma$.
   2. Repeat for a second, unrelated system where the variation splits differently, with $\tau = 12$ and $\sigma = 20$ (for example, fish density among streams). Its ICC is $12^2/(12^2 + 20^2) = 0.26$, well below the length–weight ICC. Working both cases shows whether the best design depends on the ICC.
   3. In each case, which allocation is better, and by what percentage does it shrink the SE? Explain the result using the two terms under the square root: which one gets smaller when the 60 measurements are spread over more groups, which one stays the same, and how the ICC decides how much there is to gain.
3. Fit a random-slope model lmer(Weight ~ Length + (1 + Length | Subject)) and compare it to m1b with anova(m1b, m_slope, refit = FALSE). Is the random slope supported? Explain why refit = FALSE is correct here given that refitting would use ML.

If you leave the fish out of the model, where does the fish-to-fish variation go?

Fish differ from one another, and they also get measured imperfectly. Those are two separate sources of variation, and a mixed model can tell them apart. A fixed-effects fit has only one place to put variation, the residual, so it has to blend them.

Here you build the data yourself, so you know exactly how big each source really is, and then watch what each fit reports.

library(nlme)
set.seed(77)

nfish <- 12; ages <- 1:8
Kap <- 0.25; T0 <- 0
sd_Linf <- 12    # fish really do differ: SD of 12 in asymptotic length
sd_obs  <- 3     # and each measurement is off by about 3

Linf_i <- rnorm(nfish, 100, sd_Linf)      # each fish's own Linf
sim <- do.call(rbind, lapply(1:nfish, function(i) data.frame(
  Subject = factor(i), Age = ages,
  Length  = Linf_i[i]*(1-exp(-Kap*(ages - T0))) + rnorm(length(ages), 0, sd_obs))))
gd <- groupedData(Length ~ Age | Subject, data = sim)

# Fit 1: one curve for all fish. No way to say that fish differ.
fit_nls  <- nls(Length ~ Linf*(1-exp(-Kappa*(Age-Tzero))), data = sim,
                start = c(Linf=100, Kappa=.2, Tzero=0))

# Fit 2: each fish gets its own Linf, drawn from a distribution.
fit_nlme <- nlme(Length ~ Linf*(1-exp(-Kappa*(Age-Tzero))),
                 fixed = Linf + Kappa + Tzero ~ 1, random = Linf ~ 1 | Subject,
                 data = gd, start = c(Linf=100, Kappa=.2, Tzero=0))

1. What does each fit say the measurement error is? Run summary(fit_nls)$sigma and VarCorr(fit_nlme). You built the data with a measurement error of 3 and a fish-to-fish SD of 12. Which fit recovers both numbers, and what does the other one report instead?
2. Work out where that wrong number from nls() came from. No calculus needed. Just look at the growth equation. Length is $L_\infty$ multiplied by the factor $(1 - e^{-\kappa a})$, and that factor is the fraction of its final size a fish has reached by age $a$.

   So two fish whose $L_\infty$ differ by 12 are not 12 apart at every age. At age $a$ they are $12 \times (\text{that fraction})$ apart, close together when young and far apart when old.

   1. Compute the fraction at each age: 1 - exp(-0.25 * 1:8). How far apart are two such fish at age 1? At age 8?
   2. At any age the fish are spread out by both that gap and the measurement error. Squares add, so the spread at age $a$ is sqrt((12*fraction)^2 + 3^2). Compute all eight.
   3. nls must describe all eight ages with one number. Take sqrt(mean(SD^2)) of your eight values and compare it with the $\sigma$ that nls reported in (a). Close?
3. What that does to the residual plot. Your eight numbers in (b) should grow with age. But nls used a single $\sigma$ for all of them, so its residuals must be too small at young ages and too large at old ones. Draw plot(fitted(fit_nls), resid(fit_nls)) and describe the shape.

   Here is the trap. That picture looks exactly like non-constant variance, and the usual fixes are a variance function or a log transform. Neither would help. Say in one or two sentences what is actually wrong with the model.

4. What you gain by separating the two. Compare the standard errors of Kappa and Tzero between the fits (summary(fit_nls)$coefficients and summary(fit_nlme)$tTable). Neither parameter was given a random effect, yet one fit gives much tighter values. Every standard error is scaled by how noisy the model thinks the data are. Use that to explain the difference.
5. How many fish do you need? Change one line of the simulation code, nfish, and leave everything else exactly as it is: sd_Linf stays 12, sd_obs stays 3, Kap and T0 stay as given. Those are the true values you are simulating from, so they do not change; what changes is how many fish the model gets to learn them from. Re-run the simulation block and refit fit_nlme with nfish <- 4, and again with nfish <- 20, and each time report the estimated fish-to-fish SD from VarCorr(fit_nlme). Compare each estimate with the value it is trying to recover, which is 12 every time. What does the pattern say about how many groups you need before a random effect is worth fitting?

Fish trophic position, and a random-effect structure the data cannot support

The QCBS dataset qcbs_w6_data.csv gives trophic position vs body length for three fish species (S1, S2, S3) sampled across six lakes (L1–L6). The biological question is whether trophic position increases with body length, and does that relationship differ by species and by lake?

This problem has a second purpose. One of the candidate models below asks the data for something it cannot give, and it is the model that wins the model-selection table. Learning to spot that is the point.

dat <- read.csv(here::here("Day3_MixedModels_ModelSelection", "data", "qcbs_w6_data.csv"))
dat$Fish_Species <- as.factor(dat$Fish_Species)
dat$Lake         <- as.factor(dat$Lake)
dat$Z_Length     <- scale(dat$Fish_Length)[,1]   # z-score for numerical stability
dat$Z_TP         <- scale(dat$Trophic_Pos)[,1]

M1 <- lmer(Z_TP ~ Z_Length + (1|Fish_Species) + (1|Lake), data=dat, REML=TRUE)
M3 <- lmer(Z_TP ~ Z_Length + (1|Fish_Species), data=dat, REML=TRUE)
M4 <- lmer(Z_TP ~ Z_Length + (1|Lake), data=dat, REML=TRUE)
M7 <- lmer(Z_TP ~ Z_Length + (1|Fish_Species) + (1+Z_Length|Lake), data=dat, REML=TRUE)
M8 <- lmer(Z_TP ~ Z_Length + (1+Z_Length|Fish_Species) + (1|Lake), data=dat, REML=TRUE)

library(AICcmodavg)
aictab(list(M1=M1,M3=M3,M4=M4,M7=M7,M8=M8))

1. Why are both Lake and Fish_Species candidates for random effects rather than fixed? Give the argument that would push you the other way.
2. What two random effects estimate. The term (1 + Z_Length | Fish_Species) gives each species its own deviation from the overall intercept, $b_{0j}$, and from the overall slope, $b_{1j}$. The model does not estimate those three species’ deviations as free parameters. It estimates how such deviations vary across a population of species, and that variation is summarized by a variance–covariance matrix: $$\begin{pmatrix}\sigma_0^2 & \rho\,\sigma_0\sigma_1\\ \rho\,\sigma_0\sigma_1 & \sigma_1^2\end{pmatrix}$$ The diagonal holds the variance of the species intercepts, $\sigma_0^2$, and of the species slopes, $\sigma_1^2$. The off-diagonal is their covariance, built from $\rho$, the correlation between a species’ intercept and its slope. VarCorr(M8) reports $\sigma_0$ and $\sigma_1$ in its Std.Dev. column and $\rho$ in its Corr column. How many parameters does this matrix cost, and what are they? What would a positive $\rho$ mean for these fish?
3. Report the AICc table and name the winner. Then run isSingular() on all five fits:

   sapply(list(M1=M1,M3=M3,M4=M4,M7=M7,M8=M8), lme4::isSingular)

   Two of them come back TRUE, and one of those two is the model that just won. What did R print when you fitted them, and did you notice it at the time?

   Which parameter broke? Run VarCorr(M8) and look at the Corr column. The estimated correlation between species intercepts and species slopes is exactly $1.000$, not “close to 1” but sitting on the largest value a correlation can take. That is what a boundary (singular) fit means. It is the correlation that broke because M8 has only three species, and so only three (intercept, slope) pairs to learn it from. Three values are enough for a rough idea of how spread out the intercepts are, or the slopes. They say almost nothing about a correlation, because three points often fall close to a straight line even when there is no relationship at all, so the estimate can run to $\pm 1$ whatever the truth is. Report the correlation and the two species standard deviations, and say which of the three you would trust.

4. Ask for less: the || notation. Writing a double bar instead of a single one keeps both random effects but forces their correlation to zero, so the model estimates two parameters instead of three:

   # (1 + Z_Length | Fish_Species)   intercepts, slopes, AND their correlation  — 3 parameters
   # (1 + Z_Length || Fish_Species)  intercepts and slopes, correlation fixed at 0 — 2 parameters
   
   M8u <- lmer(Z_TP ~ Z_Length + (1+Z_Length||Fish_Species) + (1|Lake), data=dat, REML=TRUE)
   lme4::isSingular(M8u)
   VarCorr(M8u)

   The correlation is a genuinely estimated parameter, so dropping it lowers the parameter count by exactly one. You will see that as $K = 6$ for M8u against $K = 7$ for M8 in the AICc table, and as 1 degree of freedom in a likelihood-ratio test between them. Does M8u converge? Compare its two standard deviations with M8's. Has dropping the correlation cost you anything you actually wanted to know?
5. The judgment call. Add M8u to the AICc table. M8 still comes out ahead by about 6 AICc units, and anova(M8u, M8, refit = FALSE) returns $p = 0.004$, so both criteria say the correlation parameter is "worth it". Yet that parameter is pinned at its boundary, and AICc charged M8 a full parameter for a correlation stuck at 1.000.

   Which model do you report, and why? Your answer should say what a likelihood-ratio test and an AICc difference can and cannot tell you about a parameter sitting on the edge of its allowed range.

6. Three species and six lakes is on the edge of what should be treated as random at all. When would you switch species to a fixed effect instead, and what would you lose by doing so?

Detection probability across regions, a binomial GLMM

You fit a binomial mixed model to detection/non-detection data (a camera-trap or PIT-array design; the algebra is identical). Cameras are nested within regions:

m <- glmer(Detect ~ scale(Covariate) + (1 | Region/Camera), data=d, family=binomial)
VarCorr(m)
# Groups        Name        Variance Std.Dev.
# Camera:Region (Intercept) 1.05     1.02
# Region        (Intercept) 0.46     0.68
# Fixed effects:
#                   Estimate Std.Error z
# (Intercept)       -1.4      0.30     -4.6
# scale(Covariate)   0.62     0.13      4.8

1. Interpret the output in words. (i) The fixed-effect estimate for scale(Covariate) is $0.62$ on the logit scale. In which direction does detection change as the covariate increases, and is the effect statistically clear (look at the estimate relative to its standard error)? (ii) Compare the two random-effect standard deviations, Region ($\sqrt{0.46}=0.68$) and Camera-within-Region ($\sqrt{1.05}=1.02$): at which level, among regions, or among cameras within a region, is there more unexplained variation in detectability, and what would that suggest biologically?
2. For a camera at mean covariate value with average region and camera intercepts, what is the expected detection probability? (Use plogis().)
3. Translate the covariate slope into an odds ratio per 1 SD, with an approximate 95% CI.
4. Why include the camera-level (rather than only the region-level) random intercept at all? What does dropping it do to your standard errors?

Nested counts and overdispersion with observation-level random effects

An agency flies aerial surveys that produce counts nested three deep: 4 management units × 6 transects per unit × 3 flights per transect = 72 transect-flights, each a single count. Counts vary among units, among transects within units, and among flights within transects.

1. Diagram the nesting, then write the equation for a Poisson GLMM with random intercepts at each level. With one count per flight, the flight-level term has exactly one observation per level. Recall from the walkthrough that in a Normal model such a term could not be separated from the residual error. A Poisson GLMM has no separate residual variance, though, so the term is identifiable there. It is the observation-level random effect (OLRE), and it is a legitimate, widely used way to handle overdispersion in Poisson counts. In one sentence, say what it adds to each count and why that lets the variance of the counts exceed their mean.
2. In lme4 syntax, what is the difference between (1 | unit/transect) and (1 | unit) + (1 | transect)? If transect IDs are reused across units (transect "1" exists in every unit but they are different transects), which do you use?
3. Now that you know about observation-level random effects, what are the two ways to address overdispersion in a mixed model for counts? For each, say what it changes about the model and how it lets the variance of the counts exceed their mean.

How often does each criterion pick noise as "best"?

Generate a response that is completely unrelated to its predictors, fit every subset with dredge(), and record how often the top model retains at least one noise variable under each criterion.

library(MuMIn)
options(na.action = "na.fail")
set.seed(126)
n <- 50; a <- 5; sdy <- 3
nummodels <- 500
pstor <- AICstor <- AICcstor <- BICstor <- rep(NA, nummodels)

for(i in 1:nummodels) {
 y <- rnorm(n, mean=a, sd=sdy)
 x1 <- runif(n); x2 <- runif(n); x3 <- runif(n)
 d <- data.frame(y, x1, x2, x3)
 reg <- lm(y ~ ., data=d)
 pstor[i]    <- sum(summary(reg)$coefficients[2:4,4] < 0.05)
 AICstor[i]  <- sum(!is.na(dredge(reg, rank="AIC")[1, 2:4]))
 AICcstor[i] <- sum(!is.na(dredge(reg, rank="AICc")[1, 2:4]))
 BICstor[i]  <- sum(!is.na(dredge(reg, rank="BIC")[1, 2:4]))
}
data.frame(p = mean(pstor > 0), AIC = mean(AICstor > 0),
           AICc = mean(AICcstor > 0), BIC = mean(BICstor > 0))

1. Start from the walkthrough loop. You built this simulation in the afternoon walkthrough and found roughly p 0.16, AIC 0.45, AICc 0.39, BIC 0.16. Run it here and report the four proportions. With no real effect anywhere, the AICc-best model still contains at least one noise predictor about 40% of the time, and the walkthrough showed why.
2. How far behind is everything else? A first-place model means little if the other models are close behind it. dredge() lists every candidate model with a delta column, which is that model’s AICc minus the best model’s AICc, and the intercept-only model is the row with no predictors at all. Create gap <- null_d <- rep(NA, nummodels) before the loop, and replace the AICc line inside the loop with these four lines:

   dd          <- dredge(reg, rank = "AICc")
   AICcstor[i] <- sum(!is.na(dd[1, 2:4]))
   gap[i]      <- dd$AICc[2] - dd$AICc[1]                     # best vs runner-up
   null_d[i]   <- dd$delta[rowSums(!is.na(dd[, 2:4])) == 0]   # intercept-only model's delta

   Report the median gap, the fraction of runs in which the best model is more than 2 AICc units ahead of the runner-up, and the fraction in which the intercept-only model is within 2 AICc units of the best.
3. What a reader needs. A paper reports that the top AICc model includes x2. Using your numbers from (b), explain in two or three sentences why “ranked first” on its own is weak evidence that x2 matters, and name the numbers a reader would need in order to judge the claim.

When BIC misses a weak-but-real effect

Three predictors have known effect sizes, one strong, one modest, one weak. Ask which criterion keeps the weak one.

library(MuMIn); options(na.action="na.fail")
set.seed(11); n <- 60
strong <- runif(n); modest <- runif(n); weak <- runif(n)
y <- 0.8*strong - 0.4*modest + 0.1*weak + rnorm(n, 0, 0.5)
d <- data.frame(y, strong, modest, weak)
m <- lm(y ~ strong + modest + weak, data=d)
dredge(m, rank="AIC"); dredge(m, rank="AICc"); dredge(m, rank="BIC")

1. Run the single simulation. Which predictors does each criterion retain in its top model?
2. Wrap it in a 500-rep loop (resampling y and the predictors each time) and report the proportion of reps in which each criterion retains weak. Do the same for modest and strong.
3. Two criteria, two questions. AIC is efficient: it estimates how well a model will predict new observations, and for large samples it ranks models the same way leave-one-out cross-validation does. BIC is consistent: if the model that generated the data is among the candidates, BIC selects it with probability approaching 1 as the sample size grows, and a difference in BIC between two models approximates the log of a Bayes factor, a measure of the evidence for one model over the other.
   1. Which criterion answers “does including weak improve prediction of new observations?”, and which answers “do the data provide evidence that weak belongs in the model that generated them?” Justify each match.
   2. Using your retention rates from (b), say how far each answer can be trusted at $n = 60$. In particular, if BIC leaves weak out, is that evidence that weak has no effect?
4. Now change the sample size. Run the loop from (b) again with n = 500, and also rerun the pure-noise loop from Problem 6 at n = 500. Fill in this table with the rate at which each criterion keeps the weak-but-real predictor (its power, one minus the Type II error) and the rate at which it keeps at least one pure-noise predictor (its Type I error):

                   keeps weak (real)      keeps noise (Type I)
                   n = 60    n = 500      n = 50    n = 500
   AIC             ____      ____         ____      ____
   AICc            ____      ____         ____      ____
   BIC             ____      ____         ____      ____

   One criterion's Type I error barely moves with $n$ and one collapses. Explain both from the size of the penalties, $2$ against $\ln n$, and then say what that means for a criterion's ability to find a weak effect as data accumulate. Which behaviour is the consistency described in (c), and which is the price AIC pays for being efficient? Which would you want when the goal is to identify the predictors in the model that generated the data, and which when the goal is to predict new observations?

AIC vs LOOCV on the trophic-position data

Return to qcbs_w6_data.csv from Problem 3, but now treat species and lake as ordinary fixed factors so we can rank candidate fixed-effect models by both an information criterion and cross-validation. The question is predictive. Which set of predictors best forecasts a fish's trophic position?

library(boot); library(AICcmodavg)
dat <- read.csv(here::here("Day3_MixedModels_ModelSelection", "data", "qcbs_w6_data.csv"))
dat$Fish_Species <- as.factor(dat$Fish_Species)
dat$Lake         <- as.factor(dat$Lake)
dat$Z_Length     <- scale(dat$Fish_Length)[,1]

forms <- list(
 null    = Trophic_Pos ~ 1,
 length  = Trophic_Pos ~ Z_Length,
 len_sp  = Trophic_Pos ~ Z_Length + Fish_Species,
 len_lk  = Trophic_Pos ~ Z_Length + Lake,
 full    = Trophic_Pos ~ Z_Length + Fish_Species + Lake,
 interact= Trophic_Pos ~ Z_Length * Fish_Species + Lake)

# LOOCV MSE (delta[1]) and AICc for each candidate
cv   <- sapply(forms, function(f) cv.glm(dat, glm(f, data=dat))$delta[1])
aicc <- sapply(forms, function(f) AICc(glm(f, data=dat)))
data.frame(cv_MSE = round(cv, 4), AICc = round(aicc, 1))

1. Report the table. Which model minimizes AICc, and which minimizes the LOOCV MSE? Do they agree?
2. Compute $\Delta$AICc for every model and the Akaike weights. Is there a single clearly-best model, or a cluster within $\Delta\text{AICc} < 2$?
3. Compare the best model's LOOCV MSE to var(dat$Trophic_Pos). Roughly what fraction of the response variance does the model's prediction error remove? Is the "best" model actually good in absolute terms, or just best of a weak field?
4. If AICc and LOOCV disagreed, which would you trust for a pure prediction task, and why?

Model averaging and asking the right question

A size-at-age analysis is fit with four competing candidate models. M1 has a fixed temperature effect, M2 a fixed regional effect, M3 both, and M0 is intercept-only. The AICc table:

Model   K   AICc     deltaAICc   weight
M3      5   1240.0   0.0         0.45
M1      4   1240.6   0.6         0.33
M2      4   1241.9   1.9         0.17
M0      3   1244.7   4.7         0.04

1. Write the Akaike-weight formula. Confirm the weights in the table from the $\Delta$AICc column, and list which models are "competitive" by the $\Delta\text{AICc} < 2$ rule.
2. Why can model-averaging coefficients mislead when the candidate models contain different, correlated predictors? Why is model-averaging predictions usually safer?
3. If the temperature coefficient in M3 is $\hat\beta_T = 0.03$ with SE 0.02, what does that say about the practical importance of temperature, regardless of which model "wins"? What single number would you also report so a reader can judge how much of the inference rides on one model?