Worked solutions — Hierarchical & latent-state Bayes

Exercise 1 — "Vague" priors on a log link

The point: a $\mathcal{N}(0, \sigma=100)$ prior on a regression coefficient is wildly informative once you exponentiate it.

set.seed(1)
beta <- rnorm(1e5, 0, 100)
eff_raw <- exp(beta * 3) # 3-degree latitudinal change
mean(abs(log10(eff_raw)) > 6) # 0.963
eff_sc <- exp(beta * 1) # 1 SD on scaled x
mean(abs(log10(eff_sc)) > 6) # 0.890

Interpretation. With this seed, 96.3% of the prior on the raw-scale model puts the implied effect outside a million-fold range (the exact analytic value is $2[1-\Phi(\log(10^6)/300)] = 0.9633$) — far beyond anything physically possible. Scaling reduces this only to 89.0% (analytic 0.8901), which is nowhere near enough; you also need a tighter prior on $\beta$ itself. The fix that works in practice is both: scale predictors AND use $\mathcal{N}(0, \sigma \approx 2)$ priors.

Common mistake: Students who use a tighter prior but keep raw predictors often still see chains stall, because $\beta \cdot x$ can still get large when $|x|$ is large. The point is the product, not either piece alone.

Exercise 2 — Fit the ant model two ways

library(nimble); library(MCMCvis); library(coda)
ants <- read.csv("data/ants.csv")

antCode <- nimbleCode({
 for (i in 1:n) {
 y[i] ~ dpois(lam[i])
 log(lam[i]) <- b0 + b1*forest[i] + b2*lat[i] + b3*elev[i]
 }
 b0 ~ dnorm(0, tau); b1 ~ dnorm(0, tau)
 b2 ~ dnorm(0, tau); b3 ~ dnorm(0, tau)
})

## (a) Raw predictors with very vague prior -- chains crawl
## constants= holds only things that never appear on the LHS of a `~`;
## the observed counts y are the one and only data= entry.
constsA <- list(n=nrow(ants), forest=ants$forest, lat=ants$latitude,
 elev=ants$elevation, tau=1/10000)
dataA <- list(y=ants$richness)
out_a <- nimbleMCMC(antCode, constants=constsA, data=dataA,
 monitors=c("b0","b1","b2","b3"),
 nchains=3, niter=10000, nburnin=500,
 samplesAsCodaMCMC=TRUE)
# multivariate=FALSE: the default multivariate psrf needs a positive-definite
# covariance and errors out when these chains are this badly stuck.
# With these default inits the raw chains stall outright: NaN log densities,
# R-hat returns NaN. (See the instructor note for what happens with good inits.)
gelman.diag(out_a, multivariate=FALSE)

## (b) Scaled predictors, same prior -- chains converge instantly
ants$forest <- as.numeric(scale(ants$forest))
ants$latitude <- as.numeric(scale(ants$latitude))
ants$elevation <- as.numeric(scale(ants$elevation))

constsB <- list(n=nrow(ants), forest=ants$forest, lat=ants$latitude,
 elev=ants$elevation, tau=1/10000)
dataB <- list(y=ants$richness)
out_b <- nimbleMCMC(antCode, constants=constsB, data=dataB,
 monitors=c("b0","b1","b2","b3"),
 nchains=3, niter=10000, nburnin=500,
 samplesAsCodaMCMC=TRUE)
MCMCsummary(out_b)
glm(richness ~ ., data=ants, family=poisson)$coef # compare

## (c) Tighter prior
constsC <- list(n=nrow(ants), forest=ants$forest, lat=ants$latitude,
 elev=ants$elevation, tau=1/4)
out_c <- nimbleMCMC(antCode, constants=constsC, data=dataB,
 monitors=c("b0","b1","b2","b3"),
 nchains=3, niter=10000, nburnin=500,
 samplesAsCodaMCMC=TRUE)
MCMCsummary(out_c) # nearly identical to (b)

Posterior means from a representative run (scaled predictors, 3 chains × 10000 iter, 500 burn-in): $b_1$ (forest) ≈ 0.32 [0.21, 0.44], $b_2$ (latitude) ≈ −0.25 [−0.39, −0.13], $b_3$ (elevation) ≈ −0.18 [−0.30, −0.07], all with R-hat = 1.00 and n.eff ≈ 4700–5500. These match the scaled glm() coefficients (0.32, −0.25, −0.18) to two decimals. The intercept is $b_0$ ≈ 1.84 [1.71, 1.96] (glm 1.85). Because the scaled model converges cleanly, these summaries are essentially seed-independent.

Instructor note. With RAW predictors and the same vague prior the intercept and latitude coefficient fail to converge, and — this is the point worth showing — even starting all three chains at the frequentist MLE does not rescue them. Using the dispersed glm()-based inits from ant.priors.R (list(as.list(freq_glm), as.list(freq_glm + runif(4,-1,1)), ...)), the raw fit still returns R-hat in the tens for b0 and b2 (7–30 depending on the seed) with n.eff ≈ 12–17, while forest (b1) mixes fine. With the default inits instead you get NaN log densities and a stalled chain outright. The cause is geometric, not a bad starting point: un-centred latitude ($\approx 42$) is nearly collinear with the intercept, so the posterior is a long narrow ridge the sampler cannot traverse. Scaling the predictors (part b), with the same inits and the same vague prior, gives R-hat = 1.00 across the board — the one thing that changed was the parameterisation. Scale first; show it live.

Exercise 3 — Lizards: convergence, prediction, and prior sensitivity

library(nimble); library(MCMCvis); library(coda)
liz <- read.csv("data/IslandsLizards.csv")
x.scaled <- as.numeric(scale(liz$perimeterAreaRatio))

lizCode <- nimbleCode({
 a ~ dnorm(0, 1/2.71); b ~ dnorm(0, 1/2.71)
 for (i in 1:n) {
 logit(p[i]) <- a + b*x[i]
 y[i] ~ dbern(p[i])
 }
 x10 <- (10 - mx)/sx; x20 <- (20 - mx)/sx
 y10 <- ilogit(a + b*x10)
 y20 <- ilogit(a + b*x20)
 diff <- y20 - y10
})
consts <- list(n=nrow(liz), x=x.scaled,
 mx=mean(liz$perimeterAreaRatio),
 sx=sd(liz$perimeterAreaRatio))
samples <- nimbleMCMC(lizCode, constants=consts,
 data=list(y=as.integer(liz$presence)),
 monitors=c("a","b","y10","y20","diff"),
 nchains=3, niter=6000, nburnin=1000,
 samplesAsCodaMCMC=TRUE)
gelman.diag(samples[, c("a","b")])
heidel.diag(samples[, c("a","b")])
b.post <- as.matrix(samples)[, "b"]
mean(b.post < 0) # ~ 0.99+
MCMCsummary(samples, params=c("y10","y20","diff"))

Result (actual run; set.seed(2025), 3 chains, niter=6000, nburnin=1000, n=19 islands). All R-hat = 1.00 (gelman.diag upper C.I. ≤ 1.01) and heidel.diag stationarity passes. $P(b < 0)$ = 0.9995, confirming higher PA reduces occupancy. Predicted occupancy falls from $y_{10}$ = 0.73 [0.47, 0.91] to $y_{20}$ = 0.41 [0.17, 0.68], so diff has posterior median −0.30, 95% CI [−0.54, −0.10] — credibly negative. (Refitting with the vague dnorm(0, 1e-6) prior instead pushes the slope to $b$ mean −4.97 [−10.2, −1.5], versus the tight-prior −2.44 [−4.5, −0.8]: with only 19 near-separated islands the prior visibly regularizes $b$, though both agree $P(b<0)\approx1$.)

Biology. Higher perimeter:area ratio means more edge habitat per unit interior, which on small or irregular islands corresponds to more predation pressure and less stable microhabitat — consistent with the lizards' lower occurrence probability there.

Exercise 4 — All five N2O models

Use Multilevel_Nimble.R top to bottom. Key reports:

Actual run: set.seed(10), 3 chains, niter=10000, nburnin=1000; 563 measurements, 107 sites, 10 fertilizer types; all $\beta$ R-hat = 1.00.

The site-level covariate is $w_j = \log(c_j/(100-c_j))$, the logit of site $j$'s mean soil organic carbon percentage $c_j$ (the carbon column).

(b) Shrinkage plot. Plotting a.no (model 2) against a.part (model 3): data-rich sites sit on the 1:1 line, single-measurement sites are pulled hard toward $\mu_\alpha$. Concretely, mean $|\Delta\alpha_j|$ = 0.46 for $n_j=1$ sites vs 0.06 for $n_j\ge10$ sites (correlation of $n_j$ with $|\Delta\alpha_j|$ = $-0.49$).

(c) Group-level effect. $\eta$ posterior median = 0.30, 95% CI [0.01, 0.59], $\Pr(\eta>0)=0.98$ — credibly positive: sites with higher soil organic carbon emit more N₂O at a given N input, consistent with more labile carbon fuelling denitrification. The effect is modest and the interval only just clears 0, because $\eta$ is informed by the number of sites (107), not the 563 observations. $\eta$ and $\kappa$ are correlated and mix slowly, so the default niter=10000 gives R-hat $\approx 1.04$ for $\eta$; the numbers above are from a longer run (3 chains, 60k iter, thin 5, R-hat = 1.00), or equivalently from centering $w$.

(d) Slopes by fertilizer. Steepest is UM (fert.index 5) at $\beta$ = 1.41 [0.76, 2.33], a urea-containing fertilizer; the other urea/mixed types (U 0.94, AM 0.98, UA 0.90, unknown 1.02) are also high, while the mineral N type is shallowest (0.60). Urea-type fertilizers do have the steepest emission-vs-input slope, as expected.

Instructor note. In this run the random-slope model converged fine at niter=10000 (all R-hat $\le$ 1.01, including $\sigma_\beta$). If a student's hyperparameters $\sigma_\alpha, \sigma_\beta$ show R-hat > 1.05, lengthening niter and thinning fixes it.

Exercise 5 — Beverton-Holt: fit, predict, derive

Run SpawnerRecruit_Nimble.R top to bottom.

(a) Posterior summaries (set.seed(2025), 3 chains, niter=10000, nburnin=1000; Nyears=20). Every 95% CI contains the truth (R-hat = 1.00–1.01):

(b) Conservative management value. 5th percentile of $H_\text{msy}$ = 0.46 (median is 0.57). That is the "if I want to be 95% confident I'm not over-fishing" harvest rate — noticeably below the point estimate.

(c) More data, tighter posteriors. Refitting with $N=100$ years shrinks the $H_\text{msy}$ 95% CI from width 0.228 (20 yr: [0.44, 0.67]) to 0.069 (100 yr: [0.53, 0.60]), a factor of 3.3×. That is a bit more than the $\sqrt{5}\approx2.24$ the pure $\sqrt{N}$ rule predicts: at 20 years the BH curvature is weakly identified ($a$'s posterior runs to 9.25, near its prior ceiling), inflating the small-sample interval; 100 years both add data and pull the posterior into its well-identified regime. The lesson still holds — more data buy tighter management quantities, on the order of $\sqrt{N}$.

Instructor note. The log() calls in this model cannot produce NaN: with $a \in [1,10]$, $b \in [1000,10000]$ and $S_t > 0$, the argument $1 + (a/b)S_t$ is always strictly greater than 1, so log(1 + a/b * S[t]) is always well-defined and positive. (And the script supplies fixed inits — $a$ = 3, 5, 7 across the three chains — so there is no "random initial $a$" to worry about either.) The failure modes that do show up: a student who deletes inits entirely and lets NIMBLE draw from the priors will still run, just with a slower burn-in; and a student who changes the simulation truth to $\sigma > 1$ will silently hit the sigma ~ dunif(0, 1) ceiling and wonder why $\sigma$'s posterior stacks up at 1 — the same prior-bound lesson as the callout above.

Exercise 6 — Dynamic occupancy

Run DynamicOccupancy_Nimble.R.

(a-c) From an actual run (set.seed(2025); 3 chains, 20000 iter, 10000 burn-in, thin 5; of the 45 monitored parameters, R-hat = 1.00 for all but the derived $\lambda_3$, which is 1.01): the quantities informed by all the data are recovered well — detection $p$ = 0.63 [0.55, 0.70] (truth 0.60) and baseline occupancy $\psi$ = 0.33 [0.18, 0.51] (truth 0.37). Posterior medians for $\phi_t, \gamma_t, \psi_t$ plotted vs truth scatter around the 1:1 line but with large per-parameter uncertainty (e.g. $\phi_3$ = 0.52 [0.16, 0.89], $\gamma_1$ = 0.47 [0.26, 0.70]) due to the small sample size (30 sites, 3 reps, 7 transitions). This is honest small-sample uncertainty, not a convergence failure.

(d) Without z-init trick. Try:

bad.z <- matrix(rbinom(Nsites*Nyears, 1, 0.3), Nsites, Nyears)
inits.bad <- list(list(z=bad.z, phi=runif(Nyears-1,.1,.9),
 gamma=runif(Nyears-1,.1,.9), p=.5, psi=.4))
nimbleMCMC(DynOccCode, constants=constants, data=data, inits=inits.bad, ...)

NIMBLE will report:

Error: Dynamic Index out of bounds, or initial values produce a density of 0

or similar — because for any (i,t) where the random $z[i,t] = 0$ but the data has $y[i,j,t] = 1$, the Bernoulli likelihood is 0 and $\log 0 = -\infty$. Always initialise $z$ at 1.

(e) Why CIs are wide. Each transition $t \to t+1$ has only 30 sites of evidence. A standard fix is to give $\phi, \gamma$ themselves a hierarchical prior (e.g. logit$(\phi_t) \sim \mathcal{N}(\mu_\phi, \sigma_\phi^2)$), which lets years borrow strength from each other and shrinks year-specific estimates. You can also add covariates like temperature or land-use change to explain year-to-year variation.

Exercise 7 — N-mixture build-up

Actual run: set.seed(2025), 3 chains, niter=10000, nburnin=1000, thin 5; all R-hat = 1.00–1.02.

(a) Basic. Posterior medians $\lambda$ = 6.24 [5.06, 7.78], $p$ = 0.80 [0.67, 0.87] (truth $\lambda=6, p=0.8$) — both cover truth, intervals wide with only 20 sites.

(b) Covariate on abundance. $\beta_1$ = 1.00 [0.94, 1.07], credibly positive and matching the simulation truth ($\beta_1 = 1$); $\beta_0$ = 2.99 [2.91, 3.07].

(c) Covariate on both. $\beta_1$ barely moves — from 1.00 to 1.04 [0.97, 1.10] — and the detection coefficients are recovered ($a_0$ = 0.92 [0.70, 1.11], $a_1$ = 1.00 [0.91, 1.09]). These data are clean, so modeling detection here neither biases nor much tightens $\beta_1$; in messier data, ignoring detection heterogeneity inflates apparent abundance heterogeneity.

(d) Overdispersion + Bayesian p-value. Data simulated from a negative binomial, fit with Poisson: Bayesian p-value = 0.008 ($T(\text{data})$ median 367 vs $T(y^\text{rep})$ median 332), far from 0.5 — clear evidence the Poisson model can't capture the variance. The discrepancy test catches a model that still returns a confident-looking $\lambda$ = 5.81 [5.46, 6.17].

Exercise 8 — Add a Poisson-lognormal random effect to fix overdispersion

Nmix5Code <- nimbleCode({
 for (i in 1:Nsites) {
 eps[i] ~ dnorm(0, tau.eps)
 log(lambda[i]) <- b0 + eps[i]
 N[i] ~ dpois(lambda[i])
 for (j in 1:Nreps) {
 y[i,j] ~ dbinom(prob=p, size=N[i])
 y.pred[i,j] ~ dbinom(prob=p, size=N[i])
 e[i,j] <- p*lambda[i]
 resid[i,j] <- pow(pow(y[i,j], 0.5) - pow(e[i,j],0.5), 2)
 resid.pred[i,j] <- pow(pow(y.pred[i,j],0.5) - pow(e[i,j],0.5), 2)
 }
 }
 fit.data <- sum(resid[1:Nsites, 1:Nreps])
 fit.pred <- sum(resid.pred[1:Nsites, 1:Nreps])

 b0 ~ dunif(-5, 5)
 p ~ dbeta(1, 1)
 sig.eps ~ dunif(0, 5)
 tau.eps <- 1/sig.eps^2
})

Expected. Bayesian p-value moves toward 0.5 (usually 0.3-0.7). $\sigma_\epsilon$ posterior median is around 0.5-0.7 with 95% CI well away from 0, capturing the overdispersion the Poisson alone missed.

Conceptual. The Poisson-lognormal is the workhorse fix for overdispersed N-mixture data. Some authors prefer a negative-binomial latent abundance (N[i] ~ dnegbin(...)), which is equivalent in the limit but parameterised differently.