Mixed models & model selection: the graded hand-in

Which streams get pulled hardest toward the mean?

The walkthrough used a balanced design, with every stream sampled the same number of times. Real survey effort is never balanced, and the imbalance is exactly where partial pooling earns its keep. Here you simulate an unbalanced design, so you know the truth, and watch what the model does with it.

library(nlme)
set.seed(413)

nvisit <- c(rep(2, 8), rep(5, 8), rep(30, 8))   # 8 streams each at 2, 5, 30 visits
tau <- 12;  sig <- 20;  mu <- 80              # TRUE among-stream SD, within SD, grand mean

stream_mean <- rnorm(length(nvisit), mu, tau)      # each stream's real density
Streams <- do.call(rbind, lapply(seq_along(nvisit), function(i)
  data.frame(Stream = factor(i),
             Density = rnorm(nvisit[i], stream_mean[i], sig))))

m <- lme(Density ~ 1, random = ~ 1 | Stream, data = Streams)

nopool  <- tapply(Streams$Density, Streams$Stream, mean)   # each stream on its own
partial <- coef(m)[, 1]                                    # partially pooled
n_i     <- as.numeric(table(Streams$Stream))

1. Predict before you compute. Partial pooling keeps a fraction $B_i$ of each stream's own mean and pulls the rest toward the grand mean, where
   $B_i = \dfrac{\tau^2}{\tau^2 + \sigma^2/n_i}$
   Using the true $\tau = 12$ and $\sigma = 20$, work out $B_i$ by hand for a stream with 2 visits and for one with 30. Which group will be pulled further, and roughly how many times further? Write the numbers down before running anything else.
2. Now check. Extract $\hat\tau$ and $\hat\sigma$ with VarCorr(m), compute $B_i$ for every stream, and tabulate the mean $|$nopool − partial$|$ separately for the 2-visit, 5-visit and 30-visit streams. Did your prediction hold?
3. Explain in two or three sentences why the formula behaves this way. What is $\sigma^2/n_i$ measuring, and what happens to $B_i$ as $n_i$ grows?
4. One of the 2-visit streams will have an unusually extreme raw mean, simply because two observations are not many. Find it, report its no-pooling and partially-pooled estimates, and say which you would put in a report and why.
5. A colleague proposes dropping every stream with fewer than 5 visits, "because those means are unreliable." Using your answer to (b), explain what partial pooling already does about unreliable means, and what dropping the streams would cost.

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("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 and compute the ICC. Contrast it with the ICC implied by Problem 1's true values ($\tau = 12$, $\sigma = 20$) and explain what the difference tells you about the two designs.
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.

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 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? Re-run the whole thing with nfish <- 4, then with nfish <- 20, refitting fit_nlme each time, and report the estimated fish-to-fish SD. It was built as 12 in every case. 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("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 does the term (1 + Z_Length | Fish_Species) assume about the population of slopes across species? Write out the implied bivariate Normal for the (intercept, slope) pair, and count how many parameters that costs. Two variances plus what?
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?
4. Find the parameter that 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 precisely on the largest value a correlation is allowed to take. That is what boundary (singular) fit means.

   Why did the correlation break rather than one of the two variances? You have three species, so three (intercept, slope) pairs, which is three points. Rather than reason about it in the abstract, watch what three points can and cannot tell you. Both simulations below use data with a known answer:

   set.seed(4)
   
   # Correlation from 3 points, when the TRUE correlation is zero
   r3 <- replicate(20000, cor(rnorm(3), rnorm(3)))
   median(abs(r3))        # the typical |r| you get from three points
   mean(abs(r3) > 0.9)     # how often it looks nearly perfect
   
   # Standard deviation from 3 points, when the TRUE sd is 1
   s3 <- replicate(20000, sd(rnorm(3)))
   median(s3)
   mean(abs(s3 - 1) < 0.25)  # how often it lands within 25% of the truth

   1. Run the first block. Every one of those 20,000 correlations came from data with no relationship at all. What is the typical $|r|$, and how often does it come back above 0.9?
   2. Run the second. How far off is a standard deviation estimated from three points, and how often is it roughly right?
   3. Now put the two together. In two or three sentences, say why M8's correlation was the parameter that failed. Which of the two quantities can three groups say something useful about, and which one will happily return an extreme answer no matter what the truth is?

   If you want to see it directly, plot one of the samples: plot(rnorm(3), rnorm(3)). Three points often look like a line.

5. 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)

   Does it converge? Compare its two standard deviations with M8's. Has dropping the correlation cost you anything you actually wanted to know?
6. 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.

   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.

7. 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 the overdispersion check

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, what happens to the flight-level term, and what is the standard trick for using it anyway?
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. For the Poisson GLMM, how do you check for overdispersion? Write the Pearson-dispersion calculation and state the threshold that should worry you.
4. If the data are overdispersed, name two concrete fixes and say what each does.

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. Bring that loop across. There is nothing new to discover in re-running it, and the rest of this problem extends it rather than repeating it.
2. Now apply the rule the lecture warns about. Taking the AICc-best model is one decision rule; a stricter one is only accept the winner if it beats the runner-up by more than 2 AICc units, and otherwise decline to choose. Modify the loop to record both whether the top model contains a predictor and the gap dd$AICc[2] - dd$AICc[1]:

   dd  <- dredge(lm(y ~ ., data=d), rank="AICc")
   noise_top[i] <- sum(!is.na(dd[1, 2:4]))
   gap[i]       <- dd$AICc[2] - dd$AICc[1]

   Report three numbers: how often the winner is more than 2 units clear; among those runs, how often it still contains a noise predictor; and the overall fraction of runs in which the strict rule keeps a noise predictor.
3. The strict rule cuts the overall false-positive rate a long way. Explain how it does that. Look at how often it declines to pick anything at all, and at the noise rate conditional on a clear winner existing. Is the rule choosing better models, or mostly refusing to choose? What does the median gap between best and runner-up tell you about how often a clear winner exists when nothing is real?
4. Is AICc uniquely bad, or is this just Type I error? Every selection rule keeps noise sometimes, so the 2-unit rule only means something next to a benchmark. Repeat the gap calculation for BIC (dredge(reg, rank="BIC"), gap db$BIC[2] - db$BIC[1]) and record the same three numbers, plus the plain p < 0.05 rate you already have. Fill in this table:

                         picks noise   winner >2 clear   noise GIVEN clear
   p < 0.05 (any of 3)      ____%            --                 --
   AICc best                ____%          ____%              ____%
   BIC  best                ____%          ____%              ____%
   
   median gap to runner-up:  AICc ____    BIC ____

5. Two things in that table need explaining. (i) BIC declares a clear winner far more often than AICc does. Work out why from the size of the two penalties at $n = 50$, and check your reasoning against the median gaps. (ii) BIC is wrong less often when it commits, yet its overall rate of "confident and wrong" is higher than AICc's. Explain how both can be true at once, and say which of the two numbers you would want to know before trusting a published model-selection result.

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. Which criterion would you use if the science question is "is the weak predictor worth studying further?", and which if it is "does the weak signal generalize out of sample?" Justify each choice.
4. Now change the sample size. Run the loop from (b) again with n = 500, and also rerun the pure-noise loop from Problem 7 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 of the two behaviours would you want from a criterion used for inference, and which for prediction?

AIC vs LOOCV on the trophic-position data

Return to qcbs_w6_data.csv from Problem 4, 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("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? State one reason AICc could mislead on a dataset this size.

Model averaging, evidence ratios, 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. The biologist's actual question is "does temperature affect size-at-age?" Re-frame the table to answer that question directly. Compute the summed weight of models containing temperature vs those without, and interpret it as an evidence ratio rather than "which model wins."
4. 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?