Mixed models & model selection — the graded hand-in

Streams: complete pooling vs no pooling vs partial pooling

The file stream_density.csv holds fish-density measurements from 100 streams, three replicate observations per stream. Load it with read.csv():

library(nlme); library(lme4)

TheData <- read.csv("data/stream_density.csv")
Streams <- data.frame(TheData)
Streams$Stream <- factor(Streams$Stream)

lm1 <- lm(Density ~ 1, data=Streams)                 # complete pooling
lm2 <- lm(Density ~ Stream - 1, data=Streams)        # no pooling
lm3 <- lme(fixed = Density ~ 1, random = ~ 1 | Stream, data=Streams) # partial pooling

1. Write out the equation each of the three models fits, and state how many parameters each estimates. Which model would you use if you wanted to make inferences about fish density in other Oregon streams not in this dataset, and why?
2. Plot residuals from lm1 by stream with boxplot(split(resid(lm1), Streams$Stream)). What do you see, and what does it tell you about the validity of lm1?
3. Compare coef(lm2) (the 100 no-pooling stream means) with coef(lm3) (the 100 partially-pooled estimates). Pick the most extreme stream and describe how far it shrank toward the overall mean, and why.
4. From summary(lm3) and intervals(lm3), report the among-stream SD (with its 95% CI), the within-stream residual SD, and compute the ICC. Interpret the ICC in one sentence.

Length–weight: random intercepts with a fixed-effect covariate

length_weight.csv has length and weight for individuals grouped by subject. The biological question: is there a length–weight relationship, and does the intercept vary by subject?

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

1. Write out the equation the model fits and say what each parameter means biologically. Report the fixed-effect intercept and slope: what is the group-mean length–weight relationship?
2. From summary(m1) / VarCorr(m1b), report the among-subject SD of the intercept. Translate it into "subjects of the same length differ by ± X units in expected weight" (the columns are on an index scale, so keep the answer in the data's own units).
3. Report the within-subject residual SD and compute the ICC. Contrast it with the ICC from Problem 1 and explain the difference.
4. 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.

Von Bertalanffy growth — a nonlinear mixed model

fish_growth.csv contains age and length measurements for multiple individual fish (Subject). Each fish is expected to follow the von Bertalanffy growth curve

where $L_\infty$ is asymptotic length, $\kappa$ the growth rate, and $t_0$ the theoretical age at length zero.

TheData <- read.csv("data/fish_growth.csv")
xx <- as.data.frame(cbind(Subject=TheData$Subject, Age=TheData$Age, Length=TheData$Length))
AgeLen <- groupedData(Length ~ Age | Subject, data=xx)

m_fixed <- nls(Length ~ Linf*(1-exp(-Kappa*(Age - Tzero))),
 data=AgeLen, start=c(Linf=100, Kappa=0.2, Tzero=0))

m_nlme <- nlme(Length ~ Linf*(1-exp(-Kappa*(Age - Tzero))),
 fixed = Linf + Kappa + Tzero ~ 1,
 random = Linf ~ 1 | Subject,
 data = AgeLen,
 start = c(Linf=100, Kappa=0.2, Tzero=0))

1. Why should we expect $L_\infty$ to vary by individual fish in a biological sense?
2. Plot boxplot(split(residuals(m_fixed), AgeLen$Subject)). What evidence do you see that the fixed-effects residuals are not independent within subject?
3. Compare the residual SE of m_fixed (from summary(m_fixed)$sigma) with the within-subject residual SE of m_nlme. Which is smaller, and why does that have to be the case?
4. Explain what m_nlme$coefficients$fixed[1] + m_nlme$coefficients$random$Subject returns and why it is a better per-fish estimate of $L_\infty$ than fitting a separate nls to each fish.

Fish trophic position — crossed random effects, slopes & intercepts

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: does trophic position increase with body length, and does that relationship differ by species and by lake?

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.
3. Report the AICc table. Which model has the lowest AICc? For that model, read off the variance components with summary() / VarCorr() and say whether most slope variation lives across species, across lakes, or within.
4. Only 3 species and 6 lakes is on the edge of what we should treat as random. What is the practical risk (name the warning you are likely to see), and when would you switch species to a fixed effect?

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.

What does $\Delta\text{AIC} = 1.5$ actually mean?

A common rule of thumb calls any model within 2 AIC units of the best "competitive." Test what that rule really buys you when the extra variable is pure noise.

reps <- 1000
within2 <- replicate(reps, {
 n <- 1000; x <- runif(n); x2 <- rnorm(n)   # x2 is pure noise
 y <- 2 + 3*x + rnorm(n, 0, 0.2)              # only x matters
 d <- AIC(lm(y ~ x + x2)) - AIC(lm(y ~ x))
 d < 2
})
mean(within2)

1. Before running: predict roughly what fraction of reps will have the junk model within 2 AIC units of the correct one. Explain your prediction from the size of AIC's per-parameter penalty.
2. Run it. Report mean(within2). Was your prediction right?
3. Explain, in terms of the penalty $2k$, why the $\Delta\text{AIC} \le 2$ rule is structurally biased toward keeping noise variables. What should you report about the extra variable's coefficient before claiming it matters?

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. Predict before you run: which criterion will pick "at least one noise variable" most often, and which least often?
2. Run the simulation. Report the four proportions and compare to your prediction.
3. Why does BIC avoid noise variables so much better than AIC in this $n = 50$ setting? Describe a setting in which BIC's parsimony would backfire.

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

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?