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 M2. Standardising and priors solve two different problems, and you usually want both. Standardising centres the predictor and fixes the sampler. In Exercise M2 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 M2 the likelihood swamps even the $\sigma = 100$ prior, so the scaled fit lands in the right place regardless (that is exactly what Exercise M2(b) and 2(c) show). Tight priors earn their keep when data are scarce, which, in ecology, is often.

The takeaway is that 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 badly unconverged, and how badly is itself unstable, which is the point. Over ten different seeds at these settings, R-hat for b0 and b2 ranged from about 1.8 to 17 (median ≈ 4.6) and the effective sample size from 6 to 19 draws (median ≈ 10) out of 28,500 saved iterations. Do not memorise a single number here. The honest summary is "R-hat between roughly 2 and 20, n.eff in the single or low double digits, different every run." Their posterior means are correspondingly unstable. The intercept came back anywhere from 10.0 to 27.0 across those ten seeds, and disagree with glm() (which gives intercept 11.94, latitude −0.24). b3 (elevation) is the milder casualty, R-hat 1.01–2.3 depending on the seed. 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, with 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 fits 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 all three chains pass the stationarity test for both parameters, with p-values of 0.62, 0.23, 0.18 for $a$ and 0.84, 0.40, 0.48 for $b$. So yes, the chains agree and the model has converged.

(b) Three tables, not one. Chain 1 alone would have you reporting 0.62 and 0.84, comfortably the most reassuring pair of the nine numbers, and the honest summary is the whole spread, 0.18 to 0.62 for $a$. Nothing here is failing; the point is that a stationarity p-value is not an estimate of anything. It means only "no evidence against stationarity" whenever it exceeds 0.05, so it carries real sampling variation from chain to chain, and quoting the best of three makes the agreement look tighter than it is. The habit generalises. Any diagnostic reported per chain should be read across all of them, because the one that fails is the one you need to see.

(c) A flag that is not a failure. The tight-prior model converges just as well (R-hat = 1.00; stationarity p-values 0.83/0.19/0.89 for $a$ and 0.93/0.77/0.69 for $b$). The halfwidth test flags anyway, and the reason is arithmetic rather than sampling. That test asks whether the Monte Carlo error of the posterior mean is small relative to the mean itself. Roughly, is the half-width of the mean's confidence interval under some fraction of $|\hat\mu|$. Under dnorm(0, 1/2.71) the posterior mean of $a$ sits almost exactly at zero, so that denominator is near zero and the criterion is close to unsatisfiable no matter how long you run the chain. Stationarity is the diagnostic that answers "has it converged?"; halfwidth answers the different question "have I run it long enough for the precision I want?", and it is unreliable whenever the parameter is centred on zero. The right response is to ignore it here, not to add iterations.

The general lesson. A diagnostic flag is a prompt to work out what the test measures, not an instruction. Two of the three results on this page look alarming at first glance and neither one indicates a problem.

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). R-hat is 1.00 for every $\beta$ in models 1–3.

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

This is worth pausing on, because it is the usual result and it is not the one students expect. Pooling is not a knob that changes your headline slope; it changes the per-site intercepts, and it changes how much you should believe each one. The average within-site response to nitrogen input is about the same however you handle the sites.

(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 biggest movers are dominated by the data-poorest sites, though not exclusively. The top five are site 51 ($n_j=1$) shifting 0.72, site 76 ($n_j=2$) 0.64, site 39 ($n_j=1$) 0.63, site 21 ($n_j=3$) 0.62, and site 42 ($n_j=1$) 0.61. Shrinkage depends on how far a site's own estimate sits from the mean as well as on how little data it has, so a site with two or three measurements that lands far out (site 76's no-pooling $\alpha$ is 2.64, against a grand mean of 0.18) can be pulled harder than a one-measurement site that was already near the middle. The tendency is what matters, and in the aggregate it is unmistakable. 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.

Where this goes next. Partial pooling pulls a data-poor site toward the grand mean because the grand mean is the only other information the model has about it. Problem 4 asks what happens when you give the model something better to pull it toward.

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

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

(b) $a$ is barely identified by 20 years. Its interval, [3.22, 9.25], covers almost two-thirds of the prior's entire range and its upper end runs very nearly into the prior ceiling of 10. That is the tell. When a posterior presses against a bound you chose, the bound is part of the answer. Productivity is the hardest thing to estimate from a short stock-recruit series, because $a$ is the slope of the curve at the origin and 20 years rarely include many observations at low spawner abundance, the very region that identifies it. Everything built on $a$ inherits that, which is worth remembering when Problem 7 turns $a$ into a harvest rate.

(c) More data, tighter posterior, faster than $\sqrt{N}$. Regenerating from the same truth with Nyears = 100 shrinks $a$'s 95% interval from [3.17, 9.19] (width 6.02) to [4.63, 6.34] (width 1.71), a factor of 3.5×. That beats the $\sqrt{100/20}\approx 2.24$ the $\sqrt N$ rule predicts, and (b) explains why. At 20 years the posterior is inflated by that collision with the prior ceiling, so it is wider than a well-identified posterior would be. The extra data do two jobs at once. They add information and pull the parameter into the regime where it is properly identified and roughly Gaussian, so the interval tightens faster than counting observations alone would suggest. $b$ behaves the same way (width 3269 → 878, a factor of 3.7), while $\sigma$, which was well identified from the start, tightens by 3.3×. The qualitative lesson is the safe one: more data buy tighter management quantities, on the order of $\sqrt N$ once you are out of the badly-identified regime.

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 does not stop NIMBLE 1.4.2. Every site-year where a randomly-drawn $z=0$ collides with an observed $y=1$ is logically impossible, so each such observation starts at a $-\infty$ log-density, and NIMBLE says so, printing one warning: logProb of data node y[i, j, t]: logProb is -Inf. line per offending observation before it runs. It then runs to completion, and because $z$ is binary the samplers move those cells to 1 on essentially the first sweep. With a random 0/1 start the posterior comes back the same as with zst (detection $p$ median 0.63 either way). Run the identical model in JAGS and it is fatal: jags.model() aborts with Error in node y[1,1,1]: Node inconsistent with parents. The lesson to take from (d) is therefore not "it crashes" but "read the warnings". A wall of $-\infty$ log-probabilities means your initial values contradict your data; you cannot rely on the sampler rescuing you in a model with continuous latent states, and the all-ones zst makes the whole class of problem disappear.

Notice what you are left with. The intervals on the individual $\phi_t$ and $\gamma_t$ are very wide, and R-hat says that is not a convergence problem. The model is being honest about a genuinely small sample. What to do about it is Problem 10, where you measure two different fixes and find that one of them can fool you.

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 M8). 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, $\sigma_\epsilon$ comes back clearly bounded away from 0 (median 0.47, 95% CI [0.39, 0.57]) and the Bayesian p-value climbs from 0.008 to about 0.17: a big improvement, and inside the unremarkable 0.1–0.9 band, but not all the way to 0.5. It should not be, because the lognormal random effect is not the negative binomial the data actually came from, so it soaks up most of the extra variance and leaves a little behind. Two practical notes. Give this run more iterations than the script's default (the numbers above are 3 chains, 60000 iterations, 20000 burn-in, thin 10: sig.eps mixes slowly), and expect a wall of problem initializing stochastic node y.pred[i, j]: logProb is -Inf warnings at startup, which are about the posterior-predictive replicates being initialised at values impossible under the current N[i] and clear themselves within the first sweeps.