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 assert?

Practice Exercise M1 had you draw from a vague log-scale prior and look at the multiplicative effect it implies across the latitude range; Practice Exercise M2 showed what that same prior does to the chains when the predictors are unscaled. This problem picks the M1 draw back up and takes it two steps further, onto a per-unit effect you can state in a sentence, and into an argument you would actually have with a reviewer. On a log link a prior on $\beta$ is a prior on a multiplicative effect $e^\beta$, and that translation is where "uninformative" priors stop being uninformative.

1. Draw $10^5$ values from $\mathcal{N}(0,\ \sigma^2 = 10000)$, the classic "vague" choice. Report the fraction for which $e^{\beta}$ exceeds a million-fold, and the fraction for which the effect exceeds a million-fold in either direction. State the second number in a sentence beginning "before seeing any data, this prior says…".
2. Repeat with $\mathcal{N}(0,\ \sigma^2 = 4)$. Report the same fraction, and also the central 95% range of multiplicative effects this prior allows (i.e. $e^{\pm 1.96 \times 2}$). Is that range wide enough for any ecological effect you have ever measured?
3. A reviewer objects that $\mathcal{N}(0, 4)$ is "informative" and demands $\mathcal{N}(0, 10000)$ for objectivity. Using your two numbers, write the two-sentence reply you would actually send.

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. Now let the N-input slope vary by fertilizer type, $\beta_k \sim \mathcal{N}(\mu_\beta, \sigma_\beta^2)$. Report the posterior of $\sigma_\beta$ and rank the ten types by posterior median slope. Then test the ranking before you believe it: across your draws, count how often the top-ranked type actually holds the largest $\beta_k$ of the ten, and report that probability together with the type's 95% CI and the number of measurements it rests on. On this evidence, would you write “fertilizer type X produces the steepest emissions response” in a paper? One or two sentences.
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, with the reference points inside the model. Monitor $a,b,\sigma,N_\text{eq},S_\text{msy},H_\text{msy}$. Before you look at the output, work out the three true reference points by hand, substituting $a=5$ and $b=5000$ into the formulas above. Doing that arithmetic yourself is how you sanity-check any assessment. Report posterior medians and 95% CIs for all six quantities next to your hand-computed truths. Does every interval cover?
2. Read the formulas before you read the output. $H_\text{msy} = 1-\sqrt{1/a}$ contains no $b$ at all, while $N_\text{eq}$ and $S_\text{msy}$ both scale with it. Say what that means for a real assessment. Which of the three could you still estimate for a stock whose absolute abundance you have no handle on, and which ones collapse without it?
3. A harvest rate with a risk statement attached. The posterior median of $H_\text{msy}$ is, by construction, the harvest rate with a 50% chance of exceeding the true one. Report the 5th, 10th, 25th and 50th percentiles of the posterior. Which would you advise, and write the one-sentence risk statement that has to travel with it.
4. An escapement target. Report $S_\text{msy}$'s posterior median and 95% CI, and compute $\Pr(S_\text{msy} > 1000)$, $\Pr(S_\text{msy} > 1500)$ and $\Pr(S_\text{msy} > 2000)$. A manager asks you for “the” escapement target. Answer them in two sentences that use these numbers rather than hiding them.
5. Which quantity is best determined? Compute the coefficient of variation (posterior SD ÷ posterior median) for $a$, $b$, $N_\text{eq}$, $S_\text{msy}$ and $H_\text{msy}$, and rank them. You should find that two of the derived quantities are more precisely determined than the parameters they are built from, which sounds impossible if you think of uncertainty as something that only accumulates. Explain how it happens. One reason is visible in the shape of $1-\sqrt{1/a}$; the other you can see by computing cor(a, b) across your draws.

Why you cannot build a derived quantity from summary numbers

You have posterior draws of $a$ and $b$ from Problem 7. Management quantities are functions of those parameters, and there are two ways to produce an interval for them, one right and one that looks reasonable and is not.

1. Compute $S_\text{msy} = b\sqrt{1/a} - b/a$ once per draw, then take the 5th, 50th and 95th percentiles of the result. Now do it the other way. Take the 5th, 50th and 95th percentiles of $a$ and of $b$ separately, and substitute those numbers into the formula. Report both sets of three numbers and the percentage difference at the 5th percentile.
2. Explain why the two disagree. What does cor(a, b) across your draws have to do with it, and why does summarising each parameter first throw away exactly the information the derived quantity needs?
3. Now repeat the comparison for $H_\text{msy} = 1 - \sqrt{1/a}$. The two approaches agree here, to the digit. What is different about this quantity, and what does that tell you about when the shortcut is safe? State the general rule in one sentence.
4. In one sentence, contrast how a frequentist MLE fit would produce a precautionary interval for $S_\text{msy}$ (bootstrap or the 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, using the z-init trick (zst <- matrix(1L, nS, nY)) that the walkthrough works through. Report $\hat R$ for $p$ and for the $\gamma_t$.
2. Plot the posterior medians of $\gamma_t$ against year with their 95% intervals. Does the rising colonization signal come through? Report the mean width of those five intervals, and say what that width is telling you. Think about how many sites actually contribute information to a single year's transition.
3. Two ways to tighten them. Refit twice more, changing only how $\gamma_t$ is built:
   1. a hierarchical prior, $\text{logit}(\gamma_t) \sim \mathcal{N}(\mu_\gamma,\ \sigma_\gamma^2)$, so the years borrow strength from each other;
   2. a trend, $\text{logit}(\gamma_t) = g_0 + g_1 \cdot \text{year}_t$, which is a reasonable thing to try given how you simulated the data.
   For all three fits report the mean 95% interval width for $\gamma_t$ and the RMSE of the posterior medians against the gam.t you simulated. Which structure narrowed the intervals most? Which one estimated $\gamma_t$ most accurately? They are not the same model.
4. Now the trap. Refit once more with a covariate that is pure noise in place of the year index:

   set.seed(99); noise <- as.numeric(scale(rnorm(nY - 1)))

   Report the same two numbers and compare them against (c)(ii). Then reach for the obvious defence. Is the noise covariate's slope $g_1$ credibly different from zero? Report its posterior median and 95% interval, and the same for the real trend's $g_1$. Now answer in four or five sentences. If a narrower interval was what you were hoping to achieve, what have you learned by achieving it? Given what the two slopes just did, is there anything inside this fit that would have told you which covariate was real, and if not, what would? Finish by saying which of the three structures you would reach for by default when you cannot be sure your covariate is real, and why it is the safe choice.