From "I fit one Bayesian model" to hierarchical, latent-state, management-relevant inference

set.seed(1)
beta <- rnorm(1e5, 0, 100)
eff <- exp(beta * 3) # raw scale: 3 deg latitude
hist(log10(eff), breaks=200, xlab="log10(effect)")
mean(abs(log10(eff)) > 6) # 0.963

eff2 <- exp(beta * 1) # scaled scale: 1 SD
mean(abs(log10(eff2)) > 6) # still 0.890 (!) -- prior is too wide

With this seed, 96.3% of the prior implies more than a million-fold effect over the raw 3° latitudinal range, and standardising only drops that to 89.0% — still almost the whole prior.

Where the analytic 0.9633 and 0.8901 come from — and the general rule. Don't take the simulation's word for it; the fraction is a one-line normal-tail calculation, and the same line handles any prior width and any predictor range. Start from what the prior literally says: $\beta_2 \sim \mathcal{N}(0,\ \sigma_\beta = 100)$. Over a predictor change of $\Delta x$, a log-scale coefficient produces a multiplicative effect $e^{\beta_2 \Delta x}$ (the same power-law move as Day 2's back-transform), so

A "million-fold" change means $|\log_{10}(\text{effect})| > 6$. Solve that inequality for the coefficient and the whole question becomes a statement about $\beta_2$ alone:

Now standardise the coefficient. Because $\beta_2/\sigma_\beta \sim \mathcal{N}(0,1)$, divide that threshold by $\sigma_\beta$ to read it as a $z$-score, and the prior mass is just a normal two-tail:

That is the reusable result: the prior mass beyond a million-fold effect depends only on the single product $\sigma_\beta\,\Delta x$ — prior SD times predictor range. Widen the prior or stretch the range and the same formula answers it. Only now do the numbers enter. Raw latitude spans $\Delta x = 3$ degrees, so $z = 6\ln 10/(100\cdot 3) = 4.605/300 = 0.0461$ and $2\,\Phi(-0.0461) = 0.9633$. Standardising collapses the range to $\Delta x = 1$ SD, tripling the threshold on $|\beta_2|$: $z = 6\ln 10/100 = 0.138$, and the fraction still only falls to $2\,\Phi(-0.138) = 0.8901$. Standardising helped a little, but nowhere near enough — and it is worth being clear about why, because this is the whole point of the exercise. Both moves you could make act through the same product $\sigma_\beta\,\Delta x$: standardising the predictor shrinks the range $\Delta x$, and tightening the prior shrinks the SD $\sigma_\beta$. What differs is the leverage. Standardising latitude shrinks $\Delta x$ from 3° to 1 SD — a factor of 3 — which only triples the $z$-score threshold, from $0.046$ to $0.138$. On a $\mathcal{N}(0,1)$ both of those numbers sit almost on top of the mean, where the bell curve is flat, so tripling the threshold barely moves the tail mass ($0.963 \to 0.890$). To buy the protection that $\sigma_\beta = 2$ gives, you would have to shrink $\Delta x$ by a factor of about 50 — which no amount of rescaling a real predictor can do. That is the sense in which the cure is a sane prior, not standardising: at $\sigma_\beta = 2$ the scaled $z$ jumps to $6\ln 10/2 = 6.9$, and $2\,\Phi(-6.9) \approx 5\times10^{-12}$ — effectively no prior mass on million-fold effects at all.

One more point, because it is a common source of confusion when you reach Exercise 2: standardising and priors solve two different problems, and you usually want both. Standardising centres the predictor and fixes the sampler — in Exercise 2 the raw-latitude chains stall because latitude near 42° is nearly fused to the intercept, and centring breaks that. A sane prior fixes the prior-predictive — the absurd tails you just computed. They are not substitutes for each other. And a too-wide prior only actually distorts the posterior when the data are weak: with the 44 bogs in Exercise 2 the likelihood swamps even the $\sigma = 100$ prior, so the scaled fit lands in the right place regardless (that is exactly what Exercise 2(b) and 2(c) show). Tight priors earn their keep when data are scarce — which, in ecology, is often.

The takeaway: a $\mathcal{N}(0,\ \sigma = 100)$ prior on a log-scale coefficient is almost always far too wide. Use $\mathcal{N}(0,\ \sigma \approx 2)$ or $\mathcal{N}(0,\ \sigma \approx 3)$ for standardised predictors. Tight, weakly informative priors are what most of the modern Bayesian regression literature now recommends.

Numbers below are from an actual run (set.seed(10), 3 chains, niter=10000, nburnin=500): posterior mean [95% CI] with R-hat and effective sample size from MCMCsummary. Coefficients are b1=forest, b2=latitude, b3=elevation.

(a) Raw predictors, dnorm(0, tau=1/10000) (i.e. $\sigma=100$). The intercept and latitude are badly non-identified because raw latitude ($\approx$ 42°) is far from 0 and strongly correlated with the intercept, so the posterior is a long narrow ridge the sampler cannot traverse. Even with the dispersed glm()-based inits from ant.priors.R, b0 and b2 come back with R-hat in the tens (roughly 7–30, and it jumps around from run to run — itself a sign of non-convergence) and only n.eff ≈ 12–17 effective draws; their posterior means are unstable and disagree with glm() (which gives intercept 11.94, latitude −0.24). Forest (b1 ≈ 0.64, R-hat 1.00) happens to mix, but the fit is broken. If you instead leave the inits at their defaults the raw chains stall outright, with NaN log densities. Longer chains help only slowly; the real fix is scaling — and note that scaling, not the choice of inits, is what does it.

(b) Scaled predictors, same vague prior. Everything converges (R-hat = 1.00, n.eff ≈ 5000–6000) and the posterior means match glm() to two decimals:

(c) Scaled predictors, tight prior dnorm(0, 1/4) (i.e. $\sigma=2$). The posterior is essentially unchanged: b0 1.84 [1.71, 1.96], b1 0.32 [0.20, 0.44], b2 −0.25 [−0.39, −0.13], b3 −0.19 [−0.31, −0.07]. With 44 bogs the data dominate both priors, so tightening $\sigma$ from 100 to 2 moves nothing. Tight priors matter most when the data are weak.

Numbers below are from an actual run. Lizards_Nimble.R sets no seed, so we fixed set.seed(2025); both the vague and the tight model use 3 chains, niter=6000, nburnin=1000 on the $n=19$ islands.

(a) Convergence. For the vague-prior model gelman.diag returns R-hat = 1.00 for both $a$ and $b$ (upper C.I. 1.01) and heidel.diag passes the stationarity test for both ($p=0.62$ for $a$, $0.84$ for $b$). The tight-prior model behaves the same (R-hat = 1.00, stationarity passes); its only blemish is that the halfwidth precision test on $a$ flags, purely because $a$'s posterior mean sits almost exactly at 0. All three chains agree — the model converges.

(b) Sign of the slope. The slope is negative with near-certainty: $P(b<0)$ = 1.00 under the vague prior and 0.9995 under the tight prior. Higher perimeter:area ⇒ lower occupancy.

(c) Management difference. Using the tight-prior model, predicted occupancy falls from $y_{10}$ = 0.73 [0.47, 0.91] at PA = 10 to $y_{20}$ = 0.41 [0.17, 0.68] at PA = 20, so diff $= P(\text{PA}=20)-P(\text{PA}=10)$ has posterior median −0.30, 95% CI [−0.54, −0.10] — credibly negative (the interval excludes 0).

(d) Does the prior matter? Here it does. With only 19 islands and near-separation (small-PA islands are almost all occupied, large-PA almost all empty), the vague dnorm(0, 1e-6) lets the logit slope run large and diffuse: $b$ mean −4.97, 95% CI [−10.2, −1.5]. The weakly-informative dnorm(0, 1/2.71) regularises it to $b$ mean −2.44, 95% CI [−4.5, −0.8] — roughly half the magnitude and half the width. Both priors agree on the sign and the biological conclusion, but the marginal posterior of $b$ itself is much tighter and less extreme under the weakly-informative prior. This is exactly the small-$n$ regime, foreshadowed by the U-shaped-prior figure above, where the prior does real work rather than washing out.

From an actual run of Multilevel_Nimble.R (set.seed(10); 563 measurements across 107 sites and 10 fertilizer types; 3 chains, niter=10000, nburnin=1000). Every $\beta$ has R-hat = 1.00.

Model 4's site-level covariate is $w_j = \log\!\big(c_j / (100 - c_j)\big)$, the logit of site $j$'s mean soil organic carbon percentage $c_j$ (the carbon column of the data), entered on top of the within-site N-input slope $\beta$.

(a) The N-input slope $\beta$ is positive and stable across the whole pooling spectrum (posterior median [95% CI]):

Pooling barely moves the average slope; it mostly tightens the per-site intercepts. (Model 5 reports the mean slope $\mu_\beta$; individual fertilizers spread around it — see (d).)

(b) Shrinkage is largest for the data-poor sites. Comparing per-site $\alpha_j$ medians from Model 2 (no pooling) to Model 3 (partial pooling, grand mean $\mu_\alpha=0.18$), the sites that move most toward the mean are all single-measurement sites: e.g. site 51 ($n_j=1$) shifts 0.72, site 39 ($n_j=1$) 0.63, site 42 ($n_j=1$) 0.61. Averaged over all sites, $|\Delta\alpha_j|$ = 0.46 for $n_j=1$ sites versus only 0.06 for sites with $n_j\ge10$, and the correlation between $n_j$ and $|\Delta\alpha_j|$ is $-0.49$: more data ⇒ less shrinkage, exactly as the hierarchical prior intends.

(c) $\eta$ is credibly positive — higher-carbon sites emit more N₂O. Posterior median $\eta = 0.30$, 95% CI [0.01, 0.59], with $\Pr(\eta > 0) = 0.98$. The interval clears zero, though only just: a site's soil organic carbon carries real, if modest, explanatory power for its baseline emission on top of the within-site N-input slope $\beta$ — consistent with more labile carbon fuelling denitrification. Notice $\eta$ is informed by the number of sites (107), not the 563 observations, which is why it stays fairly uncertain even with plenty of data. Note also that $\eta$ and the group-level intercept $\kappa$ are correlated and mix slowly, so the default niter=10000 leaves R-hat $\approx 1.04$ for $\eta$; the numbers above come from a longer run (3 chains, 60k iterations, thin 5, R-hat = 1.00). Centering $w$ achieves the same mixing for far less compute.

(d) Fertilizer with the steepest N-input response (Model 5). Per-fertilizer slopes range from 0.60 (fertilizer N) up to 1.41 for UM (fert.index 5; 95% CI [0.76, 2.33]) — a urea-containing type. The other urea/mixed types (U 0.94, unknown 1.02, AM 0.98, UA 0.90) are also among the steepest, while the plain mineral-N type is shallowest. So urea-based fertilizers do show the steepest emission response, consistent with the agronomic literature. (At these settings $\sigma_\beta$ and the per-fertilizer slopes all reached R-hat $\le$ 1.01, so the default 10k iterations sufficed.)

From an actual run of SpawnerRecruit_Nimble.R (set.seed(2025); 3 chains, niter=10000, nburnin=1000).

Where the truth values $N_\text{eq}=4000$, $S_\text{msy}=1236$, $H_\text{msy}=0.553$ come from. These are not three separate facts to look up — each is one of the reference-point formulas above with the simulation truth $(a,b)=(5,5000)$ dropped in, and doing that arithmetic by hand is exactly how you sanity-check any assessment's output. Equilibrium is where recruits exactly replace spawners:

The harvest rate depends on the productivity $a$ alone — the stock's scale $b$ never enters, because $\sqrt{1/a}$ is the fraction of recruits that must escape to replace themselves:

and $S_\text{msy}$ carries that same $\sqrt{1/a}$ but rescaled by the capacity $b$:

The general reading is worth keeping: $a$ (productivity) alone fixes the harvest rate, while $b$ (capacity) sets the scale — so $H_\text{msy}$ is a pure function of $a$, and both $N_\text{eq}$ and $S_\text{msy}$ ride linearly on $b$. That is why the $H_\text{msy}$ posterior below inherits its shape straight from the posterior of $a$, and why the prior ceiling $a<10$ translates directly into the $H_\text{msy}$ ceiling of $1-\sqrt{1/10}=0.684$ flagged earlier. The MCMC simply runs this same arithmetic at every draw, so the correlated uncertainty in $(a,b)$ flows automatically into a full joint posterior for all three management quantities.

(a) Posterior medians and 95% CIs (Nyears = 20). Every CI contains the truth (R-hat = 1.00–1.01):

(b) Conservative harvest rate. The 5th percentile of the $H_\text{msy}$ posterior is 0.46 — the harvest rate you could allow while staying 95% confident you are not exceeding the true $H_\text{msy}$. Note how much lower this risk-based number sits than the median 0.57.

(c) More data, tighter posterior. Regenerating from the same truth with Nyears = 100 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 actually a bit more than the $\sqrt{100/20}=\sqrt5\approx2.24$ the pure $\sqrt N$ rule predicts: at only 20 years the Beverton-Holt curvature is weakly identified (note $a$'s posterior runs to 9.25, up against its prior ceiling of 10), so the 20-year interval is inflated; 100 years both add information and pull the posterior into its well-identified, roughly-Gaussian regime, so it tightens faster than $\sqrt N$ alone. The qualitative lesson holds — more data buy tighter management quantities, on the order of $\sqrt N$.

Numbers below are from actually running the script (set.seed(2025): 30 sites, 8 years, 3 reps, true $p=0.6$, baseline occupancy $\approx 0.37$; 3 chains, 20000 iter, 10000 burn-in, thin 5). Across all 45 monitored parameters, R-hat is 1.00 everywhere except $\lambda_3$, which is 1.01 — convergence is not the problem here.

(a–c) The two 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). Year-specific persistence and colonization scatter around the 1:1 line but with very wide 95% CIs — e.g. $\phi_3$ = 0.52 [0.16, 0.89], $\phi_5$ = 0.08 [0.00, 0.36], $\gamma_1$ = 0.47 [0.26, 0.70] (its truth 0.47, recovered almost exactly; others land anywhere in their wide interval). The derived yearly occupancy $\psi_t$ tracks the true trajectory, including the dip in year 3 (true 0.13, posterior 0.18 [0.08, 0.34]).

(d) Replacing the all-ones zst with random 0/1 inits triggers an "initial values produce a density/likelihood of 0" error (a $-\infty$ log-prob), because any site-year where a randomly-drawn $z=0$ collides with an observed $y=1$ is logically impossible.

(e) Each year-transition $\phi_t,\gamma_t$ is informed by only ~30 sites across a single two-year step, so the CIs are wide even though R-hat is essentially 1.00 throughout — Bayes is being honest about a genuinely small sample, not failing to converge. The standard fixes: put a hierarchical prior on the transitions, e.g. $\text{logit}(\phi_t) \sim \mathcal{N}(\mu_\phi, \sigma_\phi^2)$ so years borrow strength, and/or add a covariate model for $\phi$ and $\gamma$.

From an actual run of NMixture_Nimble.R (set.seed(2025); 3 chains, niter=10000, nburnin=1000, thin 5). All R-hat = 1.00–1.02.

(a) Basic model (20 sites, 3 reps; truth $\lambda=6,\ p=0.8$). Recovered cleanly: $\lambda$ = 6.24 [5.06, 7.78] and $p$ = 0.80 [0.67, 0.87]. With only 20 sites the intervals are wide, but both comfortably cover truth.

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

(c) Covariate on both. Adding the detection covariate ($a_0$ = 0.92 [0.70, 1.11], $a_1$ = 1.00 [0.91, 1.09], both recovering truth 1) leaves the abundance-covariate inference essentially unchanged: $\beta_1$ moves only from 1.00 to 1.04 [0.97, 1.10]. Because these data are simulated clean, modeling detection heterogeneity here neither biases nor much tightens $\beta_1$; in messier data, ignoring detection structure is what would distort it.

(d) Overdispersion + Bayesian p-value. The 200-site data were generated from a negative binomial but fit with a Poisson N-mixture. The discrepancy sits almost entirely on one side — observed $T(\text{data})$ median 367 vs replicated $T(y^\text{rep})$ median 332 — giving a Bayesian p-value of 0.008, far from 0.5. The Poisson model cannot absorb the extra variance, so it is decisively rejected; the fix is a Poisson-lognormal or negative-binomial N-mixture (Exercise 8). Note $\lambda$ = 5.81 [5.46, 6.17] still looks confident and reasonable — the GOF check is what exposes the misfit the point estimate alone would hide.

The random effect is the new part, but keep the whole goodness-of-fit apparatus from Nmix4Code — the replicate y.pred, the expected value e, the two Freeman-Tukey residual arrays, and the two fit.* sums — or part (c) has nothing to compute. Note that e[i,j] now uses the site-specific lambda[i]:

Nmix5Code <- nimbleCode({
 for (i in 1:Nsites) {
 eps[i] ~ dnorm(0, 1/sig.eps^2)
 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]) # posterior replicate
 e[i,j] <- p * lambda[i] # expected count under the model
 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]) # T(y, theta)
 fit.pred <- sum(resid.pred[1:Nsites, 1:Nreps]) # T(y.rep, theta)

 b0 ~ dunif(-5, 5)
 p ~ dbeta(1, 1)
 sig.eps ~ dunif(0, 5)
})

Monitor c("b0", "p", "sig.eps", "fit.data", "fit.pred"), then read the p-value off exactly as in Model 4:

post <- as.matrix(samples5)
mean(post[, "fit.pred"] > post[, "fit.data"])

After fitting, the Bayesian p-value moves much closer to 0.5 and the posterior of $\sigma_\epsilon$ is bounded away from 0.