MCMC Diagnostics in a Hierarchical State Frequency Model
Your team is producing state-level frequency indications for a 15-state homeowners program using a hierarchical Poisson–lognormal model: claim counts ~ Poisson(exposure × θₛ), with log θₛ ~ Normal(μ, σ). Sampling used four chains of 2,000 warm-up and 2,000 sampling iterations each. The diagnostics for the run are: R̂ = 1.09 for σ and 1.01 for μ; effective sample size 45 for σ and 1,850 for μ; 312 divergent transitions after warm-up. The prior on σ was Inverse-Gamma(0.001, 0.001), chosen because “it is the standard non-informative choice.” The posterior mean of σ is 0.62 with a 90% interval of (0.19, 1.35). Three states have fewer than 400 house-years of exposure each. The filing needs indications in three weeks, and the model owner says the results are fine because “R̂ is close to 1.”
1. Are these draws usable? Explain what each diagnostic is telling you and what you would change.
Read the worked answer
No - the run is not usable, and R̂ is the least informative of the three diagnostics.
- R̂ = 1.09 on σ. Between-chain and within-chain variance disagree for the scale parameter; threshold is 1.01. μ mixing at 1.01 while σ does not is the classic signature of a funnel: the chains agree about the location and disagree about the spread.
- ESS = 45 on σ. The 90% interval (0.19, 1.35) is effectively built from about 45 independent draws. That interval is unreliable and almost certainly too narrow - the opposite of what a committee assumes when it sees a wide range.
- 312 divergences out of 8,000 post-warm-up draws (3.9%). This is the strongest signal. Divergent transitions occur when the leapfrog integrator cannot follow the posterior's curvature, which in hierarchical models happens in the funnel where σ → 0 and the state effects become unidentified. Estimates from runs with divergences are biased, not merely noisy.
Fixes, in order of importance:
- Non-centred parameterisation: θₛ = μ + σ·zₛ with zₛ ~ Normal(0,1). This removes the funnel geometrically and is usually decisive.
- Replace Inverse-Gamma(0.001, 0.001) with a weakly informative half-Normal(0,1) or half-Cauchy(0,1) on σ. The inverse-gamma is not non-informative: it puts substantial mass on absurdly large σ and its shape near zero distorts the posterior for the shrinkage.
- Raise adapt_delta to 0.99 and increase iterations.
Accept only when R̂ ≤ 1.01, ESS ≥ 400 for every reported quantity, and divergences are zero.
2. What would you compare the Bayesian output against, and what result would make you abandon it for this filing?
Read the worked answer
The natural benchmark is empirical Bayes / Bühlmann–Straub credibility. For a Poisson-lognormal structure, EB gives θ̂ₛ = Zₛ·xₛ + (1−Zₛ)·μ with Zₛ = mₛ/(mₛ + v/a), where v is the within-state process variance and a the between-state variance of hypothetical means, both estimated from the marginal data. This is the closed-form analogue of what the MCMC is doing, and it should agree closely. If the posterior means differ materially from the EB estimates, or if the ordering of states changes, the Bayesian answer is being driven by the hyperprior or by a poorly mixed sampler, not by the data.
Second comparison: σ. The posterior mean of 0.62 implies a between-state standard deviation on the log scale of about 62%, which is very large for homeowners frequency across 15 states. Compare it with the EB moment estimate of √a. If EB gives, say, 0.25 and the posterior says 0.62, the inverse-gamma prior is almost certainly placing mass where the data do not support it. A heterogeneity of that size would also imply almost no shrinkage - check whether the state indications are in fact close to raw experience, which would contradict the whole purpose of the model.
Third: the three small states. Their posterior means should sit near μ with visibly wider posteriors. If they do not, the shrinkage mechanism is not functioning.
Abandon or defer if: after the reparameterisation σ still has ESS < 100, or divergences persist, or the indications move by more than a few percent between the two specifications. A filed rate should not rest on a hyperprior the data cannot discipline.
3. What do you tell the committee, and what is your fallback under the three-week deadline?
Read the worked answer
Message to the committee, stated plainly: the numbers from this run exist but their uncertainty is understated and the estimates may be biased; they cannot go into a filing as they stand. The specific failure is not “the model is wrong” but “the computation has not converged” - a fixable, well-understood problem.
Plan and timeline:
- Day 1–2: re-code with the non-centred parameterisation, half-Normal(0,1) prior on σ, adapt_delta = 0.99, four chains of 4,000 iterations.
- Day 3–4: re-check R̂, ESS and divergences, plus posterior predictive checks against observed state frequencies.
- If clean, compare against the EB/Bühlmann–Straub result and document the comparison. Small differences are expected and reassuring; large differences need explanation before filing.
Fallback. If the re-run does not converge cleanly, file on the empirical Bayes / Bühlmann–Straub calculation. It is closed-form, reproducible by a reviewer, requires no sampler, and is defensible as the standard credibility approach. Document in the memo that a full hierarchical model was built, that its diagnostics did not support use for the filed indications, and that the EB result was used as the basis, with the Bayesian result shown as a supporting cross-check where the two agree.
Regardless of method: present the three small states with explicit credibility or interval ranges rather than point indications, and do not reproduce the (0.19, 1.35) interval - an interval from ESS 45 is not an interval a regulator should rely on.