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 96% and 89% come from, one small step at a time. You can get both numbers without simulating. It takes one idea about logarithms and one about the bell curve.

Step 0 · What the two functions are doing. Two functions appear, and each has a plain meaning.

• $e^{x}$ (in R, exp(x)) turns a number on the log scale into a multiplier. With a log link, a slope $\beta$ means "multiply the expected count by $e^{\beta}$ for each unit of $x$". So $e^{0.5 \times 3} = 4.5$ means richness is 4.5 times higher across the 3 degrees.
• $\log_{10}$ just counts zeros. $\log_{10}(1{,}000{,}000) = 6$ because a million is 1 followed by 6 zeros, and $\log_{10}(0.000001) = -6$. So "$|\log_{10}(\text{effect})| > 6$" is shorthand for "the effect is more than a million times bigger, or more than a million times smaller".

Step 1 · Turn "a million-fold effect" into a question about $\beta$. Over the 3-degree range the effect is $e^{3\beta}$. When is that bigger than a million? The natural log, log() in R, is the function that undoes exp(), so ask what power of $e$ gives a million:

Divide both sides by 3: the effect passes a million when $\beta > 4.6$. Check it: exp(3 * 4.6) is about 985,000, just short of a million. The same argument on the other side says the effect shrinks by more than a million when $\beta < -4.6$. So the whole question becomes: how often does the prior draw a slope more than 4.6 away from zero?

Step 2 · Compare 4.6 with the width of the prior. The prior is a bell curve centred at 0 with standard deviation 100. Most of its draws land somewhere between about $-200$ and $+200$. The "safe" band, $-4.6$ to $+4.6$, is a thin sliver in the middle of that: 4.6 is less than one-twentieth of one standard deviation (4.6 ÷ 100 = 0.046). Only about 3.7% of the prior falls inside the sliver, so about 96.3% falls outside it. In R, the share of a Normal curve beyond a cut-off is 2 * pnorm(-cutoff / SD):

log(1e6) / 3                   # 4.6, the slope that gives a million-fold effect over 3 degrees
2 * pnorm(-4.6 / 100)         # 0.963, the share of the prior beyond it

Step 3 · Standardise latitude and repeat. Now one unit of $x$ is one standard deviation of latitude, so the effect of a one-unit step is $e^{\beta}$, with no factor of 3. The cut-off becomes $\beta > 13.8$ instead of 4.6. That is three times further out, but it is still small next to a prior SD of 100 (13.8 ÷ 100 = 0.14), so the sliver is still thin: about 11% inside, 89% outside.

2 * pnorm(-13.8 / 100)        # 0.890

Step 4 · Now shrink the prior instead. Keep latitude standardised, so the cut-off is still 13.8, but use a prior SD of 2. Now 13.8 is almost 7 standard deviations from the centre (13.8 ÷ 2 = 6.9), far out in the tail where the bell curve is essentially zero:

2 * pnorm(-13.8 / 2)          # about 5 in a trillion

Under this prior half of all draws imply less than a 3.9-fold change per standard deviation of latitude, and 95% imply less than about 52-fold. Those are effects a biologist can believe.

What the table says. Standardising the predictor moved the cut-off further out, but only by a factor of 3, and 13.8 is still tiny next to 100, so almost nothing changed (96% to 89%). Shrinking the prior from 100 to 2 moved the cut-off from a fraction of one SD to seven SDs, and the absurd effects vanished. The cure for absurd prior effects is a sensible prior, not standardising. The general rule behind every row: compare the slope that would give an absurd effect with the prior's standard deviation. If that slope is well inside one SD, the prior is claiming absurd effects are likely.

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 samples 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.

library(nimble); library(MCMCvis); library(coda)
set.seed(10)
ants <- read.csv(here::here("Day5_BayesII", "data", "ants.csv"))

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

# Fit the model to a data frame with a given prior precision tau.
# Inits start chain 1 at the glm() estimates and jitter the other two.
fit_ants <- function(dat, tau) {
  start <- coef(glm(richness ~ forest + latitude + elevation, data = dat, family = poisson))
  names(start) <- c("b0", "b1", "b2", "b3")
  inits <- list(as.list(start), as.list(start + runif(4, -1, 1)), as.list(start + runif(4, -1, 1)))
  nimbleMCMC(antCode,
             constants = list(n = nrow(dat), forest = dat$forest, latitude = dat$latitude,
                              elevation = dat$elevation, tau = tau),
             data = list(y = dat$richness), inits = inits,
             monitors = c("b0", "b1", "b2", "b3"),
             nchains = 3, niter = 10000, nburnin = 500, samplesAsCodaMCMC = TRUE)
}

# (a) raw predictors, sigma = 100 (tau = 1/10000)
raw <- fit_ants(ants, tau = 1/10000)
MCMCsummary(raw, round = 2)          # look at the Rhat and n.eff columns
MCMCtrace(raw, pdf = FALSE)          # b0 and b2 wander instead of overlapping
gelman.diag(raw)

# (b) scaled predictors, same prior
ants_s <- ants
for (v in c("forest", "latitude", "elevation")) ants_s[[v]] <- as.numeric(scale(ants[[v]]))
scaled <- fit_ants(ants_s, tau = 1/10000)
MCMCsummary(scaled, round = 2)
round(coef(glm(richness ~ forest + latitude + elevation, data = ants_s, family = poisson)), 2)

# (c) scaled predictors, tight prior sigma = 2 (tau = 1/4)
tight <- fit_ants(ants_s, tau = 1/4)
MCMCsummary(tight, round = 2)

(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 about 20 draws (median ≈ 10) out of 28,500 saved iterations. The code above, at set.seed(10), gives R-hat 3.4 and n.eff 21 for both. 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.00–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 samples 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.

set.seed(2025)
source(here::here("Day5_BayesII", "Lizards_Nimble.R"))
library(coda)
# (a) the vague-prior fit
gelman.diag(samples1[, c("a", "b")])
heidel.diag(samples1[, c("a", "b")])

# (b) heidel.diag returns one table per chain; line up the stationarity p-values
sapply(heidel.diag(samples1[, c("a", "b")]), function(h) h[, "pvalue"])

# (c) the tight-prior fit: stationarity passes, halfwidth flags on a
heidel.diag(samples2[, c("a", "b")])
MCMCsummary(samples2, params = c("a", "b"), round = 3)   # a's posterior mean is near 0

(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.41/0.78/0.35 for $a$ and 0.36/0.24/0.22 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$ is close to zero, about −0.20 against a posterior SD of 0.58, so that denominator is small. Each chain's halfwidth is about 0.036, roughly 18% of $|\hat\mu|$, over the default 10% threshold, and because the target itself is near zero the criterion stays hard to satisfy 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.

# Run models 1-3 of Multilevel_Nimble.R first (everything above the MODEL 4 header),
# so samples1, samples2, samples3, group and nsites exist.

# (a) the N-input slope beta under each pooling choice: median and 95% CI
cols <- c("50%", "2.5%", "97.5%")
rbind(complete = MCMCsummary(samples1, params = "beta")[, cols],
      none     = MCMCsummary(samples2, params = "beta")[, cols],
      partial  = MCMCsummary(samples3, params = "beta")[, cols])

# (b) shrinkage: no-pooling vs partial-pooling site intercepts
a.no   <- MCMCpstr(samples2, params = "alpha", func = median)$alpha
a.part <- MCMCpstr(samples3, params = "alpha", func = median)$alpha
plot(a.no, a.part, pch = 19, xlab = "alpha[j], no pooling", ylab = "alpha[j], partial pooling")
abline(0, 1, col = "grey")

n_j   <- as.vector(table(group))      # measurements per site (group.index 1..107)
shift <- abs(a.no - a.part)           # how far each site moved
shrink <- data.frame(site = 1:nsites, n_j, no_pool = a.no, partial = a.part, shift)
head(shrink[order(-shrink$shift), ], 5)                 # the biggest movers
tapply(shift, cut(n_j, c(0, 1, 9, Inf), labels = c("1", "2-9", "10+")), mean)
cor(n_j, shift)
MCMCsummary(samples3, params = "mu_alpha", round = 2)

(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).

source(here::here("Day5_BayesII", "SpawnerRecruit_Nimble.R"))
# (a) medians and 95% CIs against the truth (a = 5, b = 5000, sigma = 0.2)
MCMCsummary(samples, params = c("a", "b", "sigma"), round = 3)

# (b) how much of the prior's range, U(1, 10), does a's 95% interval cover?
ci_a <- quantile(as.matrix(samples)[, "a"], c(0.025, 0.975))
ci_a
diff(ci_a) / (10 - 1)

# (c) simulate and fit a fresh dataset of any length, same truth and priors
fit_bh <- function(Nyears) {
  S <- rlnorm(Nyears, meanlog = 7.5, sdlog = 0.8)
  R <- (a.true * S) / (1 + (a.true / b.true) * S) * exp(rnorm(Nyears, 0, sigma.true))
  nimbleMCMC(BHcode, constants = list(Nyears = Nyears, S = S), data = list(log.R = log(R)),
             inits = inits, monitors = c("a", "b", "sigma"),
             nchains = 3, niter = 10000, nburnin = 1000, samplesAsCodaMCMC = TRUE)
}
width <- function(s, par) diff(quantile(as.matrix(s)[, par], c(0.025, 0.975)))

set.seed(2025); fit20  <- fit_bh(20)
set.seed(2025); fit100 <- fit_bh(100)
MCMCsummary(fit20,  params = c("a", "b", "sigma"), round = 3)
MCMCsummary(fit100, params = c("a", "b", "sigma"), round = 3)
sapply(c("a", "b", "sigma"), function(par) unname(width(fit20, par) / width(fit100, par)))
sqrt(100 / 20)

(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 6 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.22, 9.25] (width 6.03) to [4.58, 6.28] (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 3224 → 893, a factor of 3.6), while $\sigma$, which was well identified from the start, tightens by 3.0×. 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.

source(here::here("Day5_BayesII", "DynamicOccupancy_Nimble.R"))
cols <- c("2.5%", "50%", "97.5%")
# (a)-(c) posterior medians and 95% CIs beside the simulated truth
MCMCsummary(samples, params = c("p", "psi"), round = 2)
cbind(truth = round(phi.true, 2),   MCMCsummary(samples, params = "phi",   round = 2)[, cols])
cbind(truth = round(gamma.true, 2), MCMCsummary(samples, params = "gamma", round = 2)[, cols])
cbind(truth = round(colMeans(z), 2), MCMCsummary(samples, params = "psi.t", round = 2)[, cols])

# (d) refit with random 0/1 starting values for z instead of the all-ones zst
bad.init <- function() {
  i <- init.fun()
  i$z <- matrix(rbinom(Nsites * Nyears, 1, 0.3), Nsites, Nyears)
  i
}
samples.bad <- nimbleMCMC(DynOccCode, constants = constants, data = data,
                          inits = list(bad.init(), bad.init(), bad.init()),
                          monitors = c("p", "psi", "phi", "gamma"),
                          nchains = 3, niter = 20000, nburnin = 10000, thin = 5,
                          samplesAsCodaMCMC = TRUE)
MCMCsummary(samples.bad, params = c("p", "psi"), round = 2)

(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 9, 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.

source(here::here("Day5_BayesII", "NMixture_Nimble.R"))
# (a)-(c) the three fits, against the simulated truth
MCMCsummary(samples1, round = 2)   # truth: lambda = 6, p = 0.8
MCMCsummary(samples2, round = 2)   # truth: b0 = 3, b1 = 1, p = 0.8
MCMCsummary(samples3, round = 2)   # truth: b0 = 3, b1 = 1, a0 = 1, a1 = 1

# (d) the Bayesian p-value, and the picture behind it
MCMCsummary(samples4, params = c("lambda", "p"), round = 2)
post <- as.matrix(samples4)
mean(post[, "fit.pred"] > post[, "fit.data"])
apply(post[, c("fit.data", "fit.pred")], 2, median)
plot(post[, "fit.data"], post[, "fit.pred"], pch = ".",
     xlab = "T(observed data)", ylab = "T(replicate data)")
abline(0, 1, col = "red")          # a good fit scatters evenly around this line

(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.06].

(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. With the chi-square-type sum as the discrepancy, the real data sit consistently further from the model's expected counts than the replicates do (observed $T(\text{data})$ median 1247 vs replicated $T(y^\text{rep})$ median 1163), giving a Bayesian p-value of 0.025, far from 0.5. Both numbers are worth reading against the yardstick from the box above: with 600 counts, a model that fitted would give $T \approx 600$. The replicates come in at about twice that, because the fitted $\lambda$ and $p$ are themselves distorted by the misfit, and the real data are higher still. The direction is the one that matters: the observed counts are more variable than the model can produce, never less. The Poisson model cannot absorb the extra variance, so it is rejected; the fix is a Poisson-lognormal or negative-binomial N-mixture (Exercise M8). Note $\lambda$ = 5.81 [5.46, 6.16] 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 chi-square 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(y[i,j] - e[i,j], 2) / e[i,j] # chi-square contribution
 resid.pred[i,j] <- pow(y.pred[i,j] - e[i,j], 2) / e[i,j]
 }
 }
 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.025 to about 0.36 (real-data $T$ median 608, replicates 597): comfortably inside the unremarkable 0.1–0.9 band. It need not land exactly on 0.5, because the lognormal random effect is not the negative binomial the data actually came from, but it absorbs the extra variance well enough that the real counts no longer look unusual. Note too that the real-data $T$ has dropped from about 1247 to 608, right on the yardstick of 600 counts. For part (d), the fitted $\beta_0 = 1.68$ and $\sigma_\epsilon = 0.47$ give a mean abundance of $m = e^{\beta_0 + \sigma_\epsilon^2/2} \approx 6.0$ and a variance of $6.0 + 6.0^2\,(e^{0.47^2} - 1) \approx 6.0 + 8.9 = 14.9$, about two and a half times the 6.0 a Poisson with the same mean allows. The negative binomial that generated the data ($\theta = 3$) has variance $6 + 6^2/3 = 18$, so the lognormal random effect recovers most, though not all, of the real extra variation, which is also why the p-value settles at 0.36 rather than 0.5. 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.