Bayesian data analysis, end to end: from exploratory plots to a report someone can act on

E01 taught the loop prior predictive -> fit -> diagnose -> posterior predictive. That loop is the easy part. This notebook is about the analysis around it:

Part 1 needs no MCMC at all. Parts 2-4 fit one small regression several times (seconds each).

#Part 1 · EDA that chooses the model

Exploratory analysis for a Bayesian model is not a tour of every column. Each plot below is there because it decides something about the model: the likelihood family, whether a hierarchy will matter, whether coefficients will be identifiable, whether observations can be treated as exchangeable. The habit to form: finish every EDA plot with the sentence "so the model needs ...".

#1.1 · The outcome's support and shape decide the likelihood

Six real outcomes, six different answers.

None of these decisions needed a model fit, and each of them matters more to the final answer than any prior you will agonise over later.

#1.2 · The mean-variance relationship picks the family

For counts and positive outcomes the support does not finish the job: Poisson, NegativeBinomial, Gamma and LogNormal all live on the positive numbers. They differ in how the variance grows with the mean, and that is something you can plot. Split the data into groups that differ in their mean, compute each group's mean and variance, and look at the slope on log-log axes:

Three datasets, three slopes:

• EPL goals: the teams scatter around the dashed Poisson line (variance/mean close to 1). With 38 matches per team the scatter is wide, so the fitted slope itself means little - what matters is that nothing sits systematically above the line. Decision: Poisson is a fine starting point (C02).
• Hourly bike rentals: the quiet night-time cells sit just above the Poisson line, the busy ones two orders of magnitude above it, and the points follow the slope-2 line - the signature of var = mu + mu^2/alpha. A Poisson here would produce absurdly narrow predictive intervals. Decision: NegativeBinomial (its alpha is roughly the number printed above) or a LogNormal/Gamma on the positive scale; the choice between those is about zeros and discreteness, not about fit.
• Penguin body mass (not a count, so no Poisson line): the group means differ by 60% but the standard deviations barely move, and the slope is far below 2. Decision: additive Normal errors on the gram scale; taking logs would impose a constant CV that these groups do not show. We will still try the log model in Part 4 - as a sensitivity check, not as a belief - and see whether this forecast holds.

Caveat: the groups must differ in their mean for reasons inside the model. If a group mixes regimes (say, sunny and rainy days), its variance includes between-regime variation and the plot overstates dispersion. Group by the predictors you intend to use.

#1.3 · Group sizes and between/within variation forecast how much pooling will matter

Before fitting a hierarchical model you can predict what it will do. With within-group variance $\sigma^2$, between-group variance $\tau^2$ and group size $n_j$, the partial-pooling estimate of group $j$ moves a fraction $\lambda_j = \dfrac{\sigma^2/n_j}{\sigma^2/n_j + \tau^2}$ of the way from its raw mean towards the grand mean. A one-way ANOVA gives moment estimates of both variances in three lines.

The dashed funnel is what pure sampling noise around the state mean would produce. The most extreme raw means all belong to counties with fewer than ten homes, and the spread narrows as counties grow, roughly as the funnel says it should; the points outside it are the evidence that counties differ at all. The intraclass correlation is small, the typical county has five homes, and so the forecast is that the median county gets pulled more than half-way to the state mean.

Decision: a no-pooling model would mostly report noise for small counties and a complete- pooling model would erase real differences for the big ones - this is a dataset where the hierarchy changes the answers, and where the prior on $\tau$ deserves care. E02 fits it. Had the pooling factors all been below 0.05 (large groups, big between-group spread), a hierarchical model would have returned the raw means and you could have saved the effort.

#1.4 · Collinearity forecasts what the data cannot tell apart

Now the dataset for the rest of the notebook. The file is the raw field record; we do a minimal clean here (rename, split the species label, drop incomplete rows) - D01 treats the cleaning of this file properly.

Two forecasts from this picture:

1. Sign flips. Bill depth is negatively related to mass in the pooled data and positively within every species (Gentoos are heavy with shallow bills). Any model without species answers a different question from a model with it. Species goes in.
2. A ridge in the posterior. Flipper length is very nearly a species label: most of its variance lies between the species-by-sex cells. Once the cells have their own intercepts, only the remaining within-cell variation informs the flipper slope, and the slope will be negatively correlated with whichever intercepts belong to birds whose flippers are far from the centring point. Decision: centre flipper length inside the data range, expect a wider slope posterior than the pooled scatter suggests, and do not interpret the species intercepts and the slope separately. We check this forecast against the posterior in Part 3.

#1.5 · Structure left in the residuals of a quick fit

Exchangeability is a modelling assumption you can test before you rely on it. Fit something crude - least squares is fine, this is EDA - and look at the residuals along every axis the model ignores: time, space, the order of data collection.

Two findings. The 30-day mean wanders (the quick model over-predicts the first winter and under-predicts autumn 2011) and the autocorrelation is well outside the grey band for a week, not just at lag 1. Decision: a regression with independent errors would be over-confident about every coefficient - by the crude AR(1) correction the data hold fewer than half as many independent days as rows. The model needs a smooth time component (random walk, spline or GP; see E05 and C03) or an autoregressive error. And one day sits 4.5 log units below the fit: 29 October 2012, when Hurricane Sandy shut Washington down. Decision: either an explicit indicator for known events or a likelihood with tails heavy enough not to let one day move the weather coefficients. The same plot against map coordinates does the same job for spatial data (E07).

Part 1 in one line: support -> family; mean-variance slope -> dispersion; group sizes and ICC -> whether to pool; collinearity -> what is identifiable; residual structure -> what is exchangeable. All of it before any MCMC.

#Part 2 · Turning domain knowledge into priors with PreliZ

#The brief

The Palmer Station field team weighs birds to monitor body condition. To say a bird is light they need a reference: the expected mass of a healthy bird of the same species and the same structural size (flipper length). Their question for us:

Does the reference also have to depend on sex? Sexing a penguin needs a blood sample and a lab (it could not be done for nine of the birds measured here). If a male and a female of the same species and flipper length differ by less than about 200 g - roughly 5% of a bird, the size of change we would act on - one reference line per species is enough.

That fixes the estimand before any modelling: the male-minus-female difference in expected body mass at equal species and flipper length, in grams, per species, judged against 200 g. Note what it is not: it is not "how much heavier are males" (the raw difference, which includes males simply being bigger birds - 675 g vs 412 g vs 805 g in the table above). Those are different questions with different answers, and only the brief can say which one is wanted.

$$ \begin{aligned} \text{mass}_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= a_{s[i]} + d_{s[i]}\,\text{male}_i + b\,(\text{flipper}_i - 200) \end{aligned} $$

$a_s$ is the expected mass of a female of species $s$ with a 200 mm flipper, $d_s$ the sex difference we were asked about, $b$ grams per millimetre of flipper. Normal errors and an additive structure are what section 1.2 suggested; the per-species $d_s$ is there because the question is per species.

#2.1 · Priors from constraints

Nobody knows a prior standard deviation. People know ranges. pz.maxent finds the maximum-entropy member of a family with a stated mass between two bounds - the least informative distribution that still honours what you said. pz.quartile does the same from three quartiles. Here is the knowledge, stated before looking at the body-mass column:

• Adult Pygoscelis penguins weigh between about 2.7 kg (a small Adelie) and 6.3 kg (a big Gentoo): 95% of the mass for each $a_s$ in that range.
• Males are heavier in most penguins, by up to 10-20%, but we do not want the prior to decide the sign: $d_s$ centred on 0 with 95% of its mass within +-1000 g.
• Geometry bounds the slope. If a bird scaled up perfectly (mass proportional to length cubed), one more millimetre on a 200 mm, 4.2 kg bird would add 3 x 4200/200 = 63 g. Within a species, a longer flipper is a much weaker signal of overall size than that, and it will not be negative: 95% of $b$ between 0 and 65 g/mm.
• Birds of the same kind differ "by a few hundred grams": quartiles of $\sigma$ at 200, 300 and 450 g.

The dots mark the numbers we stated; the curve is what PreliZ made of them. Things to know about these helpers (PreliZ 0.28): they return a tuple (distribution, axes) - but only the distribution when called with plot=False; a parameter given in the constructor (pz.Normal(mu=0)) is held fixed and the rest are solved for; quartile warns when the family cannot hit all three numbers exactly - the Gamma above is a near miss (26/49/76% instead of 25/50/75%), which is fine for a prior. dist.summary(), dist.plot_pdf() and dist.to_pymc("name") do what they say.

#2.2 · Priors stated on the outcome scale - and a tool to triple-check

Experts find it easier still to describe observable quantities ("what would a bird weigh?") than parameters. PreliZ has an experimental helper for that direction, pz.ppe (prior predictive elicitation): give it a PyMC model and a target distribution for the outcome, and it returns model code whose priors make the prior predictive match the target. Try it on the simplest possible model, with two different numbers of observations:

The prior predictive does match the target - centre 4500 g, spread about 900 g, all of it assigned to s. But look at the prior on a: its width shrinks with the number of observations in the model (about $918/\sqrt{n}$). ppe draws datasets from the target, fits each by maximum likelihood and fits a prior to those estimates - so the prior on the mean encodes sampling noise of an n-bird average, not what we know about the mean. With 1000 rows it claims we know the species mean to +-30 g before seeing data. The function's own warning says "triple-check the results"; this is what that looks like. A statement about individual birds does not pin down a prior on a mean - you need a second statement ("the typical bird of a species could be anywhere from 2.7 to 6.3 kg"), which is exactly what we gave maxent above. We keep the 2.1 priors.

#2.3 · Many "weak" priors are one strong prior

The most common prior mistake in applied regression is invisible on any single coefficient. Standardise everything, put an innocent-looking Normal(0, 1) (or a "vague" Normal(0, 10)) on each of $p$ coefficients, and ask what that implies for a quantity you have an opinion about: the share of variance explained, or a predicted probability. We use the real penguin design matrix - species, sex, their interaction, three body measurements, island and year: 12 columns - once for a regression on standardised body mass and once for a logistic regression of sex.

With twelve independent Normal(0, 1) coefficients the prior's median $R^2$ is 0.95; with Normal(0, 10) it is certain the model explains everything, and as a classifier it believes that for nearly nine birds in ten the sex is known in advance with 99% confidence (note the log axis on the right: the end bins hold almost all the mass). No one holds those beliefs. They are an accident of stating priors one coefficient at a time: the variance of the linear predictor is the sum over coefficients, so it grows with $p$ unless the scales shrink like $1/\sqrt{p}$ (third curve), or unless you put the prior on $R^2$ itself and divide it among the coefficients (the R2D2 family of priors does that). Whether this matters depends on how much data you have - with 333 birds and 12 columns it would wash out; with 40 birds it would not - and prior predictive simulation is the only way to see it, because every marginal prior looks harmless.

#2.4 · The prior predictive check is the test

Now the running model. It is a function because Part 4 will refit it under other defensible choices; the arguments are those choices. expected is the expected mass on the gram scale, whatever the likelihood - which is what makes the estimand code in Part 3 reusable.

The elicited priors generate penguin-sized penguins: datasets spread between roughly 2 and 8 kg, practically no simulated bird outside 1-10 kg, and the observed histogram is one unremarkable member of the family. The "x20" version produces birds of minus 50 kg. It would give nearly the same posterior here (333 birds swamp it), but it is not a statement of ignorance - it is a statement that four penguins in five weigh less than 1 kg or more than 10 kg, and it would hurt as soon as a cell had three birds in it.

The implied $R^2$ deserves an honest look, given section 2.3: it leans towards 1 (median near 0.9) with a long left tail. Here that is not an accident of counting coefficients but a belief we can state and defend: the three intercepts are allowed to differ by a kilogram or more while birds of one kind differ by a few hundred grams, so of course species, sex and size are expected to explain most of the variance in a mixed sample. The "x20" prior is certain of $R^2 = 1$, which nobody believes. If the elicited lean were not defensible, the fix would be a hierarchical prior on the intercepts, not wider marginals.

#Part 3 · Analysing a posterior properly

#3.1 · Fit, and one paragraph of diagnostics

E01 covers convergence diagnostics; in a report they deserve one sentence, so we compute the ingredients of that sentence. We also store the pointwise log-likelihood and the log-prior - Part 4's sensitivity analysis needs both (pm.compute_log_likelihood, pm.stats.compute_log_prior).

Clean. (One thing to get right in the write-up: with nutpie, pm.sample() in PyMC 6.3 warms up for 400 steps, not the 1000 of older versions - the sampler records what it did in idata.posterior.attrs, so quote that rather than your memory of the defaults.)

Before interpreting anything, check the forecast from section 1.4 - that flipper length and the group intercepts would be hard to separate:

As forecast: the slope is strongly correlated with the intercepts of the birds whose flippers are far from 200 mm (negatively with Gentoo, positively with Adelie) and negatively with every sex difference - some of the male-female gap can be re-labelled as "males have longer flippers" and the data only partly resolve which. The within-group slope is less than half the pooled one. This is why the table above is not the result: the coefficients are entangled, their meaning depends on the centring point, and none of them is on its own the answer to the brief.

#3.2 · Report estimands, not coefficients

The recipe that works for any model: set the inputs to the scenarios you were asked about, push the posterior draws through, and summarise differences per draw. Here the scenario pair is "every bird as a male" versus "the same bird as a female", keeping species and flipper length as observed - an average predictive comparison.

In this linear model the comparison reproduces d exactly - a useful sanity check of the function, which will earn its keep in Part 4 where one of the models (LogNormal) has no coefficient in grams at all. The last column is the raw male-female difference from section 1.4: 100-200 g larger in every species, because it also counts males being bigger birds. Reporting it would have answered a question the team did not ask.

More quantities on the outcome scale, each one a line of arithmetic on draws:

Gentoo and Adelie differences are indistinguishable (the contrast straddles zero); the Chinstrap difference is smaller than the Adelie one with near certainty. The last line is worth dwelling on. A sex difference of several hundred grams with a tight interval sounds decisive, yet individuals overlap - about three Chinstrap pairs in ten have the female heavier: a statement about means (what d is) and a statement about individuals (what the team will meet in the field) are different estimands, separated by $\sigma$. Stakeholders nearly always hear the second when you say the first.

Predictions for stated scenarios. The team wants reference values, so give them the table they will use: expected mass (uncertainty about the mean), the range for an individual healthy bird, and the 5th percentile as a candidate "light for its kind" flag.

Two interval types, side by side, labelled. The expected mass is known to within +-60 to 90 g; an individual bird ranges over more than a kilogram. Mixing those up is the single most common way a correct analysis turns into a wrong decision.

#3.3 · Probability statements that answer the question, and the ROPE

The brief contained a number: 200 g. So the posterior should be read against it. The direct answer is $P(d_s > 200)$, already in the table above. The region of practical equivalence (ROPE) is the same idea made symmetric: declare in advance the range of values that are "as good as zero for our purposes" (here +-200 g) and ask how much of the credible interval lies inside it. az.ci_in_rope returns that percentage.

For Adelie and Gentoo the whole interval is far outside the ROPE: sex belongs in the reference. Chinstrap is the instructive case. Its interval excludes zero comfortably - a significance test would report "a significant sex difference" and stop - but it straddles 200 g, so part of the interval is inside the ROPE and the honest answer to the question that was asked is undecided: the data do not yet say whether the Chinstrap difference is big enough to matter.

Why the ROPE is not a significance test in disguise:

• It has three outcomes, not two. It can accept practical equivalence (the whole interval inside), something a p-value can never do, and it can say "undecided" without that being read as "no effect".
• The reference is a range on a meaningful scale, stated before the analysis by someone who owns the decision - not the point value zero, which no one believes exactly and which any large dataset will reject.
• It is not a calculation on a sampling distribution: more data narrows the interval and moves you out of "undecided" in either direction.

Its weak point is the same as its strength: the answer depends on the ROPE, so report the limits, who chose them, and (cheaply, from the same draws) how the conclusion changes for other limits. ci_in_rope measures a share of the interval, which depends on ci_prob; the last column, a plain posterior probability, does not.

#3.4 · Prior versus posterior: what did the data teach us?

A posterior that looks like its prior is a finding about the data (it was uninformative for that quantity) and must be reported as such. az.plot_prior_posterior overlays the two; the contraction $1 - \text{var}_\text{post}/\text{var}_\text{prior}$ puts a number on it.

Every parameter's variance shrank by at least 97%: the conclusions are the data's, not the prior's (Part 4 tests that claim from a second direction). The smallest contractions are for the slope b - the parameter section 1.4 predicted would be hardest to learn - and for the Chinstrap difference, which rests on 68 birds against 146 and 119. Note also where the posterior for b sits: in the lower half of a prior whose upper end was "perfect geometric scaling". The data agree with the reasoning that a flipper is a weak proxy for size.

#3.5 · Variance explained

$R^2$ has a Bayesian version with a posterior: per draw, the variance of the fitted means over that plus the residual variance (Gelman, Goodrich, Gabry & Vehtari 2019). az.bayesian_r2 wants the names of the mean and scale variables; az.loo_r2 is the leave-one-out version, which is the honest one when a model has many parameters.

The two agree, as they should for 8 parameters and 333 birds (a gap between them is a symptom of overfitting). Report it with its context: most of that $R^2$ is "Gentoos are big" - only a third of the variance in body mass is within species to begin with. A high $R^2$ here says nothing about how well the model resolves the sex difference.

#3.6 · Calibration - and what to write when it is not perfect

A density overlay (E01) cannot tell you whether the predictive intervals can be trusted. The probability integral transform can: if the model is calibrated, the position of each observation within its own predictive distribution is uniform. az.plot_ppc_pit uses the posterior predictive; az.plot_loo_pit the leave-one-out predictive, so that each bird is judged against a model that has not seen it. In ArviZ 1 both draw the difference between the ECDF of the PIT values and that of a uniform (a flat line at zero is perfect calibration), print the p-value of a uniformity test that remains valid for dependent PIT values, and highlight the points responsible for any departure. (The older simultaneous envelope, method="envelope", raises a TypeError in ArviZ 1.3.)

The curve never strays more than about 0.02 from zero, no points are flagged and the test finds nothing: taken over all birds, the predictive distribution is calibrated. A global check can hide compensating errors in subgroups, and section 1.4's table hinted at one - males looked more variable than females. So check the statistic the team's 5th-percentile flag depends on, the spread within each sex: