1. Setup of the example

  • Population: 804 municipalities in NY State (1960), total population known: about 13.8 million.
  • Two simple random samples of size 100 were drawn.
  • Goal: estimate the population total from just a sample, pretending the full population is unknown.
  • This is a deliberately tricky setting because municipality sizes are highly skewed: many small towns, a few very big places.

The key thing: the total is dominated by the largest units. But a simple random sample of 100 might or might not grab those large ones.

2. First try: normal model on raw data

Assume the 804 values are drawn from

$y_i \mid \theta \sim \text{Normal}(\mu, \sigma^2), \quad \text{flat prior on } (\mu, \log \sigma).$

Then use standard finite-population/Bayesian machinery to get an interval for the total.

Result for sample 1:

a 95% posterior interval like $[-5.4 \times 10^6, \; 3.7 \times 10^7]$

which is absurd because the lower end is below the observed sample total. You could “fix” it by replacing the lower end with the sample total, but that’s patching, not modeling.

This happens because the normal model is a bad description of those skewed municipality sizes.

3. Second try: lognormal model

So, transform:

$\log y_i \mid \theta \sim \text{Normal}(\mu, \sigma^2)$

with a flat prior on $(\mu, \log \sigma)$. This makes sense: populations are positive and skewed.

Now the interval for the total from sample 1 becomes something like $[5.4 \times 10^6, \; 9.9 \times 10^6]$

Much nicer and tighter. Looks good… until you check.

Posterior predictive check

Use the sample total as the test statistic:

  • simulate replicated samples from the fitted lognormal model,
  • compute their totals,
  • compare to the real sample total.

Result: every single simulated sample total is smaller than the actual sample total. That means: even though the model looked reasonable, it systematically underpredicts big totals like the one you saw. So the model is misfitting in the exact way that matters for estimating the total.

So: good-looking interval, bad model for this estimand.

4. Third try: power-transformed normal family

To get more flexibility in the tail, widen the family. Instead of just normal (φ=1) or lognormal (φ=0), use a Box–Cox–style power transform:

  • transform data by $y^{(\phi)}$,
  • model the transformed data as normal,
  • estimate φ from the data.

This family includes:

  • φ = 1 → raw normal
  • φ = 0 → lognormal
  • other φ → in-between or heavier/lighter tails

For sample 1, the posterior for φ is sharply peaked around $-1/8$. Under this extended model:

  • the 95% interval for the total becomes $[5.8 \times 10^6, \; 3.18 \times 10^7]$,
  • and the posterior predictive check now looks fine: about 15% of replicated sample totals exceed the observed total → plausible.

So for sample 1, this looks like a happy ending.

5. But then sample 2 ruins it

Repeat the exact same logic on sample 2:

  • Normal model → extremely wide interval.
  • Lognormal model → tight interval again, but posterior predictive check fails again (replicated sample totals never as large as the observed one).
  • So go to the power family again. For sample 2, φ is estimated around $-1/4$. Posterior predictive check looks fine (about half of replications bigger than observed). So it seems okay.

But: when actually simulating total population values from this fitted power family, the draws blow up:

  • median of simulated totals: $5.7 \times 10^8$ (way too big),
  • upper simulated totals up to $10^{15}$ or $10^{17}$ — completely unrealistic for NY municipalities in 1960.

So the model fit the data and passed the posterior predictive check on the sample total, but still gave nonsense for the actual population total.

Why? Because the model was allowed to have very heavy right tails, and the data (just 100 values) do not pin down the extreme tail. Yet the total is very sensitive to that tail. So tiny differences in that tail, driven by the main-body fit, produce enormous differences in the total.

6. What actually fixes it: add real prior information

The problem wasn’t Bayesianism or transformations — it was missing substantive constraints.

If one adds a simple, realistic piece of information:

“No municipality has population > $5 \times 10^6$.”

and truncate simulated values at that bound, suddenly the inferences become reasonable:

  • Sample 1 → 95% interval: $[6 \times 10^6, \; 20 \times 10^6]$
  • Sample 2 → 95% interval: $[10 \times 10^6, \; 34 \times 10^6]$

The true total (13.7 million) is inside both. So the bounded power-transformed model works.

This shows: to predict a tail-dependent quantity, you must tell the model what tails are plausible.

7. Why the plain normal model didn’t fail as catastrophically

Oddly, the original untransformed normal model gave “not terrible” totals. That’s mostly because it was implicitly assuming a sort of boundedness: for the central limit theorem to give a decent sampling distribution for the sample mean of 100 units, the underlying yi can’t be arbitrarily huge. So even though the data aren’t normal, the normal model was sneaking in a tail assumption that kept the total from exploding.

But that same normal model would do badly for other estimands (like the population median). So it wasn’t the better model — it just had a tail that happened to keep totals sensible.

8. Big lessons from this example

  1. A model can fit the observed data and still be useless for a tail-sensitive estimand. Posterior predictive checks must use test quantities relevant to the estimand (here: the sample total).
  2. A more flexible model can make things worse if it frees up the part of the distribution (the far right tail) that the data don’t inform but the estimand depends on.
  3. For tail-driven estimands (like totals, means of heavy-tailed populations, total cost, total pollution), data alone are not enough. You have to add substantive prior constraints (like an upper bound, or a realistic tail form).
  4. Sometimes the right move is not to “improve the fit” but to “state the real-world constraint.” That’s model expansion in the good sense: not just more parameters, but more realism.
  5. Some questions are robust, some aren’t.
    • Median of municipality size → similar answers under several models, no strong need for tail constraints.
    • Total population → wildly different answers unless tail is constrained.
  6. So robustness is a property of (data + model + question), not just the model. If the real-world question is inherently non-robust (total cost, total pollution, rare but huge claims), then a good Bayesian analysis should show that sensitivity, not hide it.