Hierarchical & latent-state Bayesian inference on real data

Overdispersed ant richness with a bog-level random effect

Load the ant-richness data and fit a Poisson regression of richness on the three scaled covariates, but add a bog-level (observation-level) random effect $\epsilon_i$ on the log mean to absorb overdispersion:

ants <- read.csv("data/ants.csv")
ants$forest    <- as.numeric(scale(ants$forest))
ants$latitude  <- as.numeric(scale(ants$latitude))
ants$elevation <- as.numeric(scale(ants$elevation))

1. Fit the model in NIMBLE with weakly informative priors ($\beta \sim \mathcal{N}(0, \sigma^2{=}4)$, $\sigma_\epsilon \sim \text{Uniform}(0,5)$). Report posterior medians and 95% CIs for $\beta_1,\beta_2,\beta_3$ and for $\sigma_\epsilon$.
2. Is $\sigma_\epsilon$ credibly bounded away from 0? What does that tell you about overdispersion relative to a plain Poisson GLM?
3. Compare your $\beta$ posterior medians to the frequentist glm(richness ~ forest + latitude + elevation, family = poisson) fit. Where do they agree or differ, and why?

What does a "vague" log-scale prior actually do?

Using the same ants.csv Poisson regression (you may drop the random effect and use the plain log-link GLM form here), investigate how the prior on the coefficients interacts with predictor scaling.

1. Fit the model with RAW (unscaled) predictors and a $\mathcal{N}(0, \sigma^2{=}10000)$ prior on each $\beta$. Report what happens to chain mixing and R-hat. (You may hit Inf log-densities or a stuck chain — describe what you see.)
2. Fit again with SCALED predictors and the same vague prior. Report R-hat and effective sample size, and compare the posterior medians of $(\beta_1,\beta_2,\beta_3)$ to the frequentist glm().
3. Keep the scaled predictors but tighten the prior to $\mathcal{N}(0, \sigma^2{=}4)$ (i.e. dnorm(0, 1/4)). Does the posterior move meaningfully? Explain in one sentence why (or why not).
4. Draw $10^5$ values of a slope from $\mathcal{N}(0,10000)$ and report the fraction for which the implied multiplicative effect $e^{\beta}$ exceeds a million-fold. Use this to justify a weakly informative default.

The pooling spectrum for soil N₂O emissions

Load the N₂O data and model $\log(\text{emission})$ as a function of centred $\log(\text{n.input})$, with sites indexed by group.index.

N2O   <- read.csv("data/N2OEmission.csv")
y     <- log(N2O$emission)
x     <- log(N2O$n.input) - mean(log(N2O$n.input))  # centred predictor
group <- N2O$group.index
nsites <- length(unique(group))

1. Fit the complete-pooling model $\mu_i = \alpha + \beta x_i$ and the partial-pooling random-intercept model $\alpha_j \sim \mathcal{N}(\mu_\alpha, \sigma_\alpha^2)$, $\mu_i = \alpha_{j[i]} + \beta x_i$. Report the posterior of $\beta$ from each.
2. Also fit the no-pooling model (independent $\alpha_j$). Plot the per-site $\alpha_j$ posterior medians from no-pooling (x-axis) vs. partial-pooling (y-axis) with the 1:1 line. Which sites shrink the most toward $\mu_\alpha$, and why?
3. Report the posterior of $\sigma_\alpha$. In plain language, what does its magnitude say about how different the sites are once N input is accounted for?

Modeling the random intercept with a site-level covariate

Extend the partial-pooling model from Problem 3 so the per-site mean intercept is itself regressed on a site-level covariate $w_j$ (a per-site summary derived from the data), and let the N-input slope vary by fertilizer type (fert.index):

Build $w_j$ from each site's soil organic carbon (the carbon column, a percentage) — e.g. its logit $\log\!\big(c_j/(100-c_j)\big)$ — and state how you scaled it.

1. Fit the group-level-regression model in NIMBLE. Is the posterior of $\eta$ (the effect of the site covariate on baseline emissions) credibly different from 0?
2. Add the fertilizer-varying slope $\beta_k$. Which fertilizer type has the largest absolute N-input effect? Report the posterior of $\sigma_\beta$.
3. Show that writing the covariate as a regression on the intercept, $\alpha_j \sim \mathcal{N}(\kappa+\eta\, w_j,\sigma_\alpha^2)$, is the same model as adding $w_j$ as an ordinary term in the linear predictor with a mean-zero site random effect (substitute $\alpha_j=\kappa+\eta\, w_j+u_j$). Which parameter carries the between-site effect of $w_j$, and why does shrinkage still happen?

Island lizard occupancy vs. perimeter:area ratio

Load the island-lizard data and fit logistic regression of presence on the scaled perimeter:area ratio:

liz      <- read.csv("data/IslandsLizards.csv")
x.scaled <- as.numeric(scale(liz$perimeterAreaRatio))
y        <- as.integer(liz$presence)

1. Fit in NIMBLE with the weakly informative prior $a, b \sim \mathcal{N}(0, \sigma^2{=}2.71)$ (i.e. dnorm(0, 1/2.71), chosen because it implies a near-uniform prior on $p$). Report R-hat and heidel.diag for $a$ and $b$. Use starting values near 0 (why is that necessary here?).
2. What is the posterior probability that the slope $b$ is negative (i.e., more edge per unit area $\Rightarrow$ lower occupancy)?
3. Re-fit with a genuinely vague prior $a,b \sim \mathcal{N}(0, 10^6)$. With only 19 islands, does the prior visibly change the posterior of $b$? Overlay the two posterior densities.

A management-relevant probability difference

Building on your Problem 5 fit, compute derived quantities inside the model so their full posterior is delivered by MCMC. Add the deterministic nodes

where $x_{10}$ and $x_{20}$ are the scaled versions of raw perimeter:area ratios 10 and 20 (scale them with the same mean/SD you used for the data).

1. Report the posterior median and 95% CI of $\text{diff}$. Does the interval exclude 0?
2. State in one sentence why computing $\text{diff}$ as a node inside nimbleCode gives correct uncertainty, whereas plugging point estimates of $a,b$ into the formula would not.

Beverton-Holt fit with management quantities

Simulate a spawner-recruit time series from a known truth, then fit the Beverton-Holt model in NIMBLE and define the management quantities as deterministic nodes so their posteriors come for free.

set.seed(7)
Nyears <- 25; a.t <- 5; b.t <- 5000; sig.t <- 0.3
S <- rlnorm(Nyears, 7.5, 0.8)
R <- numeric(Nyears)
for (t in 1:Nyears)
  R[t] <- (a.t*S[t])/(1 + a.t/b.t*S[t]) * exp(rnorm(1, 0, sig.t))

1. Fit the model (monitor $a,b,\sigma,N_\text{eq},S_\text{msy},H_\text{msy}$). Report posterior medians and 95% CIs. Do the CIs cover the simulated truth?
2. Plot the posterior of $H_\text{msy}$ as a histogram and mark its median.
3. Compare $S_\text{msy}$ from Bayes to what you'd get by plugging the posterior-median $a,b$ into the formula. Are they the same? Explain any difference.

The 5th-percentile ("precautionary") harvest rate

Using the posterior draws of $H_\text{msy}$ from Problem 7, quantify risk-based management targets.

1. Report the 5th percentile of the $H_\text{msy}$ posterior. Interpret it in management language: the harvest rate you'd allow if you wanted only a 5% chance of exceeding the true MSY harvest rate.
2. Refit Problem 7 with Nyears = 100 (regenerate the data). By roughly how much does the 95% CI width of $H_\text{msy}$ shrink relative to the 25-year fit? Is the gain consistent with a $\sqrt{N}$ rule?
3. In one sentence, contrast how a frequentist MLE fit would produce this precautionary number (bootstrap / delta method) with how Bayes produces it directly.

Repeat-count abundance and a goodness-of-fit check

An N-mixture model estimates the abundance of unmarked animals when you cannot count every individual: repeated visits to the same closed site let it separate true abundance ($\lambda$) from detection probability ($p$), which a raw count confounds. (It is the abundance analog of the occupancy model; see the practice lab for a fuller explanation.) Here you will simulate repeated counts, fit the Royle (2004) N-mixture model in NIMBLE, and then test the Poisson-abundance assumption with a posterior-predictive Bayesian p-value.

library(MASS)
set.seed(9)
sites <- 150; reps <- 3; lambda <- 6; p.true <- 0.5
N <- rnegbin(sites, mu = lambda, theta = 2)   # NB, not Poisson
y <- matrix(NA, sites, reps)
for (i in 1:sites) y[i, ] <- rbinom(reps, N[i], p.true)

1. Fit the basic (constant $\lambda$, constant $p$) N-mixture. Remember: dbinom(prob, size) in NIMBLE and initialise N = apply(y, 1, max). Report the posteriors of $\lambda$ and $p$.
2. Add posterior-predictive replicates $y^\text{rep}$ and a discrepancy statistic (e.g. a Freeman-Tukey or variance/mean statistic). Compute the Bayesian p-value $p_B$. What does a value far from 0.5 indicate here?
3. Propose (in words, or as a modified nimbleCode block) a fix — e.g. a Poisson-lognormal random effect on $\log\lambda_i$ — and explain how it would move $p_B$ toward 0.5.

Multi-season occupancy: recover persistence and colonization

Simulate a dynamic (multi-season) occupancy dataset and recover the year-specific persistence $\phi_t$ and colonization $\gamma_t$ from replicated detection/non-detection data with imperfect detection.

set.seed(10)
nS <- 50; nY <- 6; nR <- 3; p.t <- 0.5
phi.t <- c(0.70, 0.70, 0.80, 0.85, 0.85)
gam.t <- c(0.05, 0.15, 0.25, 0.30, 0.30)  # rising colonization
z <- matrix(NA, nS, nY); z[, 1] <- rbinom(nS, 1, 0.2)
for (i in 1:nS) for (t in 2:nY)
  z[i,t] <- rbinom(1, 1, z[i,t-1]*phi.t[t-1] + (1-z[i,t-1])*gam.t[t-1])
y <- array(NA, c(nS, nR, nY))
for (t in 1:nY) for (i in 1:nS) for (j in 1:nR)
  y[i,j,t] <- rbinom(1, 1, p.t) * z[i,t]

1. Fit the model in NIMBLE. Use the z-init trick (zst <- matrix(1L, nS, nY)) — explain in one sentence why initialising $z$ at all-1 is required.
2. Plot the posterior medians of $\gamma_t$ against year. Does the rising-colonization signal come through?
3. The individual $\phi_t$ posteriors are wide. Why? Name one change to the model that would tighten them in a real analysis.