Diagnostics

Convergence and calibration diagnostics: Simulation-Based Calibration, recovery-coverage checks, and prior-to-posterior learning.

Model diagnostics for the MMM framework.

Currently exposes learning diagnostics – how much the data updated each parameter relative to its prior (prior-to-posterior contraction, overlap, and location shift), used to flag posteriors that are over-informed by the prior rather than the data.

exception mmm_framework.diagnostics.ConvergenceWarning

Bases: UserWarning

Warning category emitted when an MCMC fit fails standard convergence checks.

Filter or escalate it explicitly, e.g.:

import warnings
from mmm_framework.diagnostics.convergence import ConvergenceWarning
warnings.simplefilter("error", ConvergenceWarning)  # treat as failure

class mmm_framework.diagnostics.RecoveryCoverageResult(targets, n_sims_requested, n_sims_effective, n_failed_fits, truth_source, sampler, L, seed, levels, elapsed_s=0.0, caveats=<factory>)

Bases: object

Full recovery-coverage run: per-target stats + run metadata.

targets: list[RecoveryTargetStat]

n_sims_requested: int

n_sims_effective: int

n_failed_fits: int

truth_source: str

sampler: str

L: int

seed: int

levels: tuple[float, ...]

elapsed_s: float = 0.0

caveats: list[str]

property headline_level: float

property all_nominal: bool

True when no target’s headline-level interval demonstrably under-covers.

worst()

Return type:

RecoveryTargetStat | None

summary()

Return type:

str

to_dashboard()

Return type:

dict[str, Any]

__init__(targets, n_sims_requested, n_sims_effective, n_failed_fits, truth_source, sampler, L, seed, levels, elapsed_s=0.0, caveats=<factory>)

class mmm_framework.diagnostics.SBCResult(params, n_sims_requested, n_sims_effective, L, n_bins, sampler, seed, alpha, elapsed_s=0.0, n_failed_fits=0, caveats=<factory>)

Bases: object

Full SBC run: per-parameter stats + run metadata.

params: list[SBCParamStat]

n_sims_requested: int

n_sims_effective: int

L: int

n_bins: int

sampler: str

seed: int

alpha: float

elapsed_s: float = 0.0

n_failed_fits: int = 0

caveats: list[str]

property all_calibrated: bool

worst()

Return type:

SBCParamStat | None

summary()

Return type:

str

to_dashboard()

Return type:

dict[str, Any]

__init__(params, n_sims_requested, n_sims_effective, L, n_bins, sampler, seed, alpha, elapsed_s=0.0, n_failed_fits=0, caveats=<factory>)

mmm_framework.diagnostics.annotate_convergence(diagnostics)

Return diagnostics with flags and converged filled in.

Mutates and returns the same dict so it round-trips through serialization.

For a frequentist fit the sampler metrics are also nulled: az.ess happily reports ≈``n_boot`` for iid bootstrap replicates and az.rhat reports NaN, and a downstream table that formats whatever it finds would render both as if they meant something.

Return type:

dict[str, Any]

mmm_framework.diagnostics.compute_convergence(trace)

Compute divergences / rhat_max / ess_bulk_min from an ArviZ trace.

Best-effort: any piece that cannot be computed comes back as None rather than raising, so a diagnostics computation never fails a fit.

Return type:

dict[str, Any]

mmm_framework.diagnostics.compute_fit_diagnostics(mmm, results=None, *, max_parameters=40)

JSON-safe model-health snapshot (schema v1) for a fitted model.

results is the MMMResults whose diagnostics dict already carries divergences/R-hat/ESS; pass None to recompute from the trace. Each layer is independently best-effort, so a learning failure still yields the convergence block (and vice versa).

Return type:

dict[str, Any]

mmm_framework.diagnostics.convergence_flags(diagnostics)

Which checks FAILED: subset of {"divergences", "rhat", "ess"}.

Empty for a frequentist fit — not because it passed, but because none of these checks apply to it. Read alongside is_converged(), which returns None for the same fit; an empty flag list on its own must never be taken as a pass (that conflation is exactly what made a bootstrap trace read as converged).

Return type:

list[str]

mmm_framework.diagnostics.coverage_from_ranks(int_ranks, L, levels=(0.5, 0.8, 0.9, 0.95))

Empirical central-interval coverage read directly off SBC ranks.

The truth lies inside the central level posterior interval exactly when its normalized rank u = (r + 0.5)/(L + 1) falls in [(1−level)/2, 1−(1−level)/2] — so an SBC run already contains every coverage number; this just states them in user language.

Return type:

list[CoverageLevelStat]

mmm_framework.diagnostics.failure_mode_guide()

Markdown table of the ways nominal coverage fails and what to do.

Kept next to the tool so every surface (chat, Validation tab, docs) tells one story. The long-form version is technical-docs/coverage-diagnostics.md.

Return type:

str

mmm_framework.diagnostics.is_converged(diagnostics)

Convergence verdict from a diagnostics dict.

Returns True/False for NUTS fits, and None when convergence is not assessable — a frequentist fit (no chain exists), an approximate fit (diagnostics["approximate"] truthy), or one with no usable R-hat/ESS/divergence signal. None is NOT “converged”; callers should surface it as “N/A”.

The frequentist branch is load-bearing rather than defensive. A (chain=1, draw=B) bootstrap trace passes every check as True: az.rhat on one chain is NaN → filtered to None, a None metric does not raise a flag, and az.ess returns ≈B because bootstrap replicates are iid. So without this the estimator is silently green everywhere the verdict is consumed — measured and recorded in technical-docs/frequentist-estimation.md §8.

Return type:

bool | None

mmm_framework.diagnostics.parameter_learning(prior, posterior, var_names=None, *, bins=60, c_strong=0.5, c_weak=0.1, ovl_dominated=0.85, z_relocated=1.0)

How much did the data teach us about each parameter, beyond the prior?

Return type:

DataFrame

Parameters

prior, posterior:

Sample containers. Each may be an arviz InferenceData (the "prior" / "posterior" group is read), an xarray Dataset, or a dict mapping a parameter name to a sample array (sample axis first). Both must describe the same parameters; vector/matrix parameters are compared element-wise.

var_names:

Optional restriction. Names match either a flattened element ("beta[0]") or a base parameter ("beta" keeps every "beta[...]" element). None compares every parameter present in both containers.

bins:

Histogram resolution for the overlap coefficient.

c_strong, c_weak, ovl_dominated, z_relocated:

Heuristic verdict thresholds (conventions, not law). Checked in order: |shift_z| >= z_relocated and contraction < c_weak -> "relocated" (the posterior moved at least one prior-sd without narrowing – evidence dominated the location); then contraction < c_weak and overlap > ovl_dominated -> "prior-dominated"; contraction >= c_strong -> "strong"; contraction >= c_weak -> "moderate"; otherwise "weak".

Returns

pandas.DataFrame

One row per parameter, sorted by contraction ascending (so the least-learned, most prior-dominated parameters sort to the top). Columns: parameter, prior_mean, prior_sd, post_mean, post_sd, contraction, contraction_robust, overlap, shift_z, post_ess_bulk, verdict.

Notes

contraction is intentionally not clipped: a negative value (posterior wider than prior) is a genuine warning sign, but it is ambiguous on its own – it can mean prior-data conflict, a likelihood that is genuinely flat in the newly-favored region, or sd inflation from heavy tails left by chains that mixed slowly into that region. Disambiguate with the companions: a large |shift_z| says the evidence dominated the location (verdict "relocated", not a failure to learn); contraction_robust >= 0 alongside contraction < 0 says the widening is tail-driven; a low post_ess_bulk says the width estimate itself may be a sampling artifact – re-run with more tune/draws before interpreting.

Informativeness, not importance. High contraction / low overlap means the data was informative about the parameter – which includes confidently pinning it near zero. It is not a statement about effect size or sign. Always read the posterior post_mean / interval alongside the diagnostic: contraction tells you the data spoke; the posterior location tells you what it said. (A sign-constrained prior such as psi = -HalfNormal makes “the posterior is below zero” near-automatic; a high contraction to a value of ~0 means the data confidently found a negligible effect, not a confirmed one.)

mmm_framework.diagnostics.plot_parameter_learning(learning, ax=None, *, top=None, metric='contraction', threshold=None)

Horizontal bar chart of per-parameter learning, colored by verdict.

metric is the column to plot ("contraction" by default; "overlap" also works). top keeps only the top least-learned parameters (the head of the sorted frame). threshold draws a reference line.

mmm_framework.diagnostics.plot_prior_posterior_overlay(prior, posterior, parameter, ax=None, *, bins=60, transform=None)

Overlay the prior and posterior sample histograms for a single parameter.

The clearest way to see whether the data moved/narrowed a parameter relative to its prior. transform (callable) is applied to both sample sets before plotting – e.g. lambda x: -x to view the signed cannibalization effect from its psi_..._raw HalfNormal parameter.

mmm_framework.diagnostics.run_mmm_sbc(model, *, n_sims=64, L=100, n_bins=20, sampler='numpyro', params=None, tune=200, chains=2, seed=0, alpha=0.05, progress=None)

Run SBC for a (built, unfitted) BayesianMMM.

Draws n_sims paired (θ*, y_sim) from the model’s prior, refits the posterior on each y_sim by swapping the observed data on the SAME model graph (pm.observe — keeps θ* and the posterior on one fixed scale), and computes per-parameter ranks.

EXPENSIVE: one refit per simulation. Use a background job; defaults are tuned for a tractable national MMM (numpyro NUTS, n_sims=64, L=100).

Return type:

SBCResult

mmm_framework.diagnostics.run_recovery_coverage(model, *, truth='posterior_mean', n_sims=24, levels=(0.5, 0.8, 0.9, 0.95), sampler='numpyro', L=200, tune=200, chains=2, seed=0, params=None, include_contributions=True, progress=None, refit=None, extra_caveats=())

Fixed-truth recovery coverage for a BayesianMMM.

Fixes every free parameter at θ* (pm.do), simulates n_sims datasets from the likelihood at θ*, refits the ORIGINAL model on each (pm.observe on the same graph — θ* and every posterior share one fixed standardized scale), and measures how often each central interval contains the truth — for scalar free parameters and per-channel total contributions.

EXPENSIVE: one refit per simulation, like SBC. Use a background job for thorough runs; n_sims < 30 only flags gross failures (the binomial CI is reported either way).

Answers: “IF the world matched my fitted model exactly, would my reported intervals cover?” Under-coverage here is a mechanical problem (approximate inference, sampler, priors); it deliberately cannot see real-world misspecification — see the module docstring.

Parameters:

• refit (Optional[Callable[[ndarray, int], Any]]) – Estimator injection. None (default) refits the PyMC model with pm.observe + NUTS, which is the only thing this used to be able to do. Supply a callable (y_sim, seed) -> idata to grade a non-PyMC estimator — the frequentist ridge + bootstrap path (#186) is why the hook exists. The simulate-from-the-model half is unchanged, so both paths are graded against the same θ* and the same central equal-tailed intervals. The callable may return anything carrying a .posterior mapping, or a plain {name: draws} mapping.

• extra_caveats (tuple[str, ...]) – Statements appended verbatim to the result’s caveats. Required for the frequentist path: this function auto-attaches an uncertainty caveat for advi/fullrank_advi only, so ridge shrinkage bias and conditional-on-selection intervals must be passed explicitly or they go unstated.

Return type:

RecoveryCoverageResult

Note

Parameters the injected estimator holds fixed (the frequentist path fixes adstock and saturation) have degenerate draws, and their coverage is not meaningful. Pass params=[...] to restrict the table to the quantities the estimator actually estimates.

mmm_framework.diagnostics.run_sbc(draw_and_fit, *, n_sims, L, n_bins=20, alpha=0.05, seed=0, sampler='generic')

Model-agnostic SBC loop.

draw_and_fit(rng) returns (theta_star, posterior_draws) where theta_star maps param→scalar truth and posterior_draws maps the same param→a 1-D array of posterior draws (length ≥ L; thinned to L).

Return type:

SBCResult

mmm_framework.diagnostics.saturation_learning(model, results=None, *, draws=2000, random_seed=0)

Did the data move the saturation parameters? (post-fit)

Thin wrapper over mmm_framework.diagnostics.learning.parameter_learning(), restricted to the saturation block and with the prior drawn for you.

Parameters:

• model (BayesianMMM) – A fitted model.

• results (Any) – Unused; accepted so the call reads like the other diagnostics.

• draws (int) – Prior draws to compare the posterior against.

• random_seed (int | None) – Seed for the prior draw.

Return type:

DataFrame

Returns:

The parameter_learning frame for the saturation parameters, sorted least-learned first. A "prior-dominated" verdict means the fitted value is the prior you chose.

Raises:

ValueError – If the model has not been fitted.

mmm_framework.diagnostics.saturation_prior_report(model, *, draws=4000, random_seed=0)

Where does the saturation prior put the curve’s elbow? (pre-fit)

Draws each channel’s saturation parameter from its prior and converts it to the half-saturation point as a fraction of the channel’s maximum observed spend — the scale saturation actually sees, since media is max-normalized before the transform.

Parameters:

• model (BayesianMMM) – A constructed BayesianMMM. It need not be fitted; the graph is built if it has not been.

• draws (int) – Prior draws per channel.

• random_seed (int | None) – Seed for reproducibility.

Return type:

DataFrame

Returns:

One row per channel with channel, saturation, parameter, elbow_q05 / elbow_median / elbow_q95 (fractions of maximum observed spend), mass_beyond_support and verdict (anchored / diffuse / unanchored).

Channels whose family has no elbow (root, none) are omitted.

Note

A prior can be perfectly anchored and the parameter still unidentified — this reports whether the prior is defensible, not whether the data will move it. Pair it with saturation_learning() after fitting.

mmm_framework.diagnostics.warn_if_not_converged(diagnostics, label='model', *, stacklevel=3)

Emit a ConvergenceWarning if the fit failed convergence checks.

Returns True if a warning was emitted.

Return type:

bool

mmm_framework.diagnostics.warn_if_saturation_prior_is_unanchored(model, *, draws=2000, random_seed=0)

Emit a UserWarning per channel whose elbow prior is unanchored.

Returns the same frame as saturation_prior_report() so a caller can act on it. Silent when every channel is anchored — this is meant to be safe to call unconditionally.

Return type:

DataFrame

Sbc

Simulation-Based Calibration (SBC) for the inference engine (Talts et al. 2018).

SBC machine-checks whether the posteriors this framework reports are calibrated — i.e. whether their credible intervals have nominal coverage — rather than merely asserting it. It belongs with the prior-predictive checks: it runs on data generated from the model’s own prior, so it validates the inference engine itself, independent of any real dataset.

Algorithm (for each of n_sims iterations):

1. draw a “true” parameter vector θ* from the prior and simulate a dataset y_sim from the likelihood at θ* (one sample_prior_predictive call gives both, jointly);

2. refit the posterior on y_sim — on the SAME model graph (so θ* and the posterior live on the identical, fixed standardized scale);

3. record the rank r = #{posterior draws ≤ θ*} of the true value within its L (thinned) posterior draws.

Under a correctly-calibrated inference procedure those ranks are Uniform{0..L}. The shape of the deviation diagnoses the failure (rank counted as draws ≤ θ*, the Talts convention):

• ∪ / U-shape (mass at the edges) → posterior too narrow / overconfident; the true value lands in the posterior tails too often → reported intervals under-cover. Fix: widen priors / loosen the noise prior; do not trust the intervals.

• ∩ / frown (mass in the centre) → posterior too wide / overdispersed; intervals over-cover (conservative). Fix: tighten priors.

• left-skew (mass at low ranks, mean rank < L/2) → posterior sits above the truth → estimates biased high. Fix: check prior centring / standardization.

• right-skew (mass at high ranks) → posterior biased low.

(NB: this is the standard Talts 2018 direction — ∪ = underdispersed/overconfident. It is the opposite of the labels in some informal write-ups; the conjugate unit test in tests/test_sbc.py pins the direction.)

SBC is expensive — one model refit per simulation — so it is an offline verification tool gated behind a background job, run once per model architecture, never per fit.

mmm_framework.diagnostics.sbc.uniformity_chisq(int_ranks, L, n_bins=20)

χ² goodness-of-fit of the rank histogram against the discrete uniform.

Returns (chi2, p_value, bin_counts). p_value large ⇒ ranks look uniform ⇒ calibrated on this parameter.

Return type:

tuple[float, float, ndarray]

mmm_framework.diagnostics.sbc.normalized_ranks(int_ranks, L)

u = (r + 0.5) / (L + 1) ∈ (0, 1); Uniform(0, 1) under calibration.

Return type:

ndarray

mmm_framework.diagnostics.sbc.classify_shape(int_ranks, L, *, chi2_p=None, alpha=0.05)

Classify the rank-histogram geometry into a calibration verdict.

Uses N-aware z-scores so thresholds scale with the number of simulations: a bias z-score on the mean normalized rank (uniform mean = 0.5) and a dispersion z-score on its variance (uniform var = 1/12). Direction follows Talts 2018 (rank = #draws ≤ θ*).

Returns {shape, mean_norm_rank, var_norm_rank, skewness, excess_kurtosis, bias_z, dispersion_z} where shape ∈ {uniform, smile(∪), frown(∩), left-skew, right-skew}.

Return type:

dict[str, Any]

mmm_framework.diagnostics.sbc.miscalibration_score(bin_counts, L, n_bins)

Total-variation distance of the rank histogram from uniform, in [0, 1).

0.5 · Σ_b |O_b/N − p_b| with p_b the exact discrete-uniform bin probability. 0 = perfectly uniform; → 1 − min_b p_b when all mass is in one bin. Monotone and comparable across parameters, so the agent can rank the worst offenders by a single number instead of “reading the picture”.

Return type:

float

mmm_framework.diagnostics.sbc.rank_hist_band(n_sims, L, n_bins=20, prob=0.95)

Simultaneous confidence band for the rank-histogram bin counts.

Returns (lower, upper) arrays of length n_bins (per-bin counts the histogram should fall inside prob of the time, jointly across bins, under calibration). For equal-width bins the band is near-constant. Built from the Multinomial null with the multiplicity-adjusted level from _simultaneous_gamma().

Return type:

tuple[ndarray, ndarray]

mmm_framework.diagnostics.sbc.ecdf_diff_band(n_sims, prob=0.95, n_points=100, n_mc=2000, seed=12345)

Säilynoja-style simultaneous band for the ECDF-difference plot.

Under calibration the normalized ranks are Uniform(0,1), so at grid point z the ECDF satisfies N·ECDF(z) ~ Binomial(N, z). Returns (z, lower, upper) for ECDF(z) − z such that the whole curve stays in the band prob of the time jointly across the grid (multiplicity-adjusted per-point level found by Monte-Carlo + binary search).

Return type:

tuple[ndarray, ndarray, ndarray]

class mmm_framework.diagnostics.sbc.SBCParamStat(name, int_ranks, L, n_sims, n_bins, bin_counts, chi2_stat, chi2_pvalue, shape, mean_norm_rank, var_norm_rank, skewness, excess_kurtosis, bias_z, dispersion_z, miscalibration, calibrated, coverage=<factory>)

Bases: object

Per-parameter SBC verdict + the integer ranks behind it.

name: str

int_ranks: ndarray

L: int

n_sims: int

n_bins: int

bin_counts: ndarray

chi2_stat: float

chi2_pvalue: float

shape: str

mean_norm_rank: float

var_norm_rank: float

skewness: float

excess_kurtosis: float

bias_z: float

dispersion_z: float

miscalibration: float

calibrated: bool

coverage: list[CoverageLevelStat]

coverage_at(level)

Return type:

CoverageLevelStat | None

to_dashboard(*, max_ranks=0)

JSON/msgpack-safe summary (numpy scalars cast to float/int).

Return type:

dict[str, Any]

__init__(name, int_ranks, L, n_sims, n_bins, bin_counts, chi2_stat, chi2_pvalue, shape, mean_norm_rank, var_norm_rank, skewness, excess_kurtosis, bias_z, dispersion_z, miscalibration, calibrated, coverage=<factory>)

class mmm_framework.diagnostics.sbc.SBCResult(params, n_sims_requested, n_sims_effective, L, n_bins, sampler, seed, alpha, elapsed_s=0.0, n_failed_fits=0, caveats=<factory>)

Bases: object

Full SBC run: per-parameter stats + run metadata.

params: list[SBCParamStat]

n_sims_requested: int

n_sims_effective: int

L: int

n_bins: int

sampler: str

seed: int

alpha: float

elapsed_s: float = 0.0

n_failed_fits: int = 0

caveats: list[str]

property all_calibrated: bool

worst()

Return type:

SBCParamStat | None

summary()

Return type:

str

to_dashboard()

Return type:

dict[str, Any]

__init__(params, n_sims_requested, n_sims_effective, L, n_bins, sampler, seed, alpha, elapsed_s=0.0, n_failed_fits=0, caveats=<factory>)

mmm_framework.diagnostics.sbc.compute_param_stat(name, int_ranks, L, *, n_bins=20, alpha=0.05, tau=0.15)

Build the full per-parameter verdict from its integer ranks.

Return type:

SBCParamStat

mmm_framework.diagnostics.sbc.build_sbc_result(int_ranks_by_param, *, L, n_sims_requested, sampler='generic', n_bins=20, alpha=0.05, tau=0.15, seed=0, elapsed_s=0.0, n_failed_fits=0)

Assemble an SBCResult from per-parameter integer ranks.

Return type:

SBCResult

mmm_framework.diagnostics.sbc.run_sbc(draw_and_fit, *, n_sims, L, n_bins=20, alpha=0.05, seed=0, sampler='generic')

Model-agnostic SBC loop.

draw_and_fit(rng) returns (theta_star, posterior_draws) where theta_star maps param→scalar truth and posterior_draws maps the same param→a 1-D array of posterior draws (length ≥ L; thinned to L).

Return type:

SBCResult

mmm_framework.diagnostics.sbc.run_mmm_sbc(model, *, n_sims=64, L=100, n_bins=20, sampler='numpyro', params=None, tune=200, chains=2, seed=0, alpha=0.05, progress=None)

Run SBC for a (built, unfitted) BayesianMMM.

Draws n_sims paired (θ*, y_sim) from the model’s prior, refits the posterior on each y_sim by swapping the observed data on the SAME model graph (pm.observe — keeps θ* and the posterior on one fixed scale), and computes per-parameter ranks.

EXPENSIVE: one refit per simulation. Use a background job; defaults are tuned for a tractable national MMM (numpyro NUTS, n_sims=64, L=100).

Return type:

SBCResult

Coverage

Interval-coverage diagnostics — does the 90% interval contain the truth 90% of the time?

Coverage is the property users implicitly assume when they read a credible interval. This module makes it checkable. Three distinct notions map to three tools:

1. Predictive coverage (observation scale): do posterior-predictive intervals contain the observed data at the nominal rate? → posterior predictive checks and the report’s calibration curve (elsewhere).

2. Engine calibration (SBC, sbc): over datasets simulated from the model’s own prior, are posterior ranks uniform? Uniform ranks ⇔ every central interval has nominal coverage on average over the prior. coverage_from_ranks() turns SBC ranks into explicit per-level coverage numbers (“your 90% interval empirically covers 78% [68–86%]”).

3. Recovery coverage at a fixed truth (run_recovery_coverage(), this module’s headline): pick ONE parameter vector θ* — the fitted posterior mean, a prior draw, or user-supplied values — simulate n_sims datasets from the model at θ*, refit each, and count how often the X% interval contains θ*, for raw parameters AND per-channel contribution estimands. Reports per-target empirical coverage with Monte-Carlo error bars plus a bias/width diagnosis: is the failure a location problem (posterior sits away from the truth) or an interval-too-narrow problem (overconfidence)?

What recovery coverage CANNOT see: it simulates data FROM the model, so it can never detect real-world misspecification (wrong adstock/saturation family, unobserved confounding, time-varying effects). Read it jointly with an external-truth check:

• under-covers here → the problem is mechanical — an approximate fit (MAP/ADVI/Pathfinder uncertainty is not calibrated), a broken sampler, or priors so tight the data cannot move the posterior to the truth;

• covers here but failed against an external answer key (e.g. the mmm_framework.synth scenarios’ synthetic_truth.json) → the gap IS structural: misspecification, confounding, or an estimand mismatch (contribution ROI and counterfactual ROI are different numbers).

The full failure-mode table lives in technical-docs/coverage-diagnostics.md and failure_mode_guide().

Interval convention: central equal-tailed percentile intervals, matching the framework’s compute_hdi_bounds (which is percentile-based despite the name), so the coverage measured here is the coverage of the intervals actually reported.

class mmm_framework.diagnostics.coverage.CoverageLevelStat(level, n, hits, coverage, ci_low, ci_high, verdict)

Bases: object

Empirical coverage of the central level interval, with MC error.

level: float

n: int

hits: int

coverage: float

ci_low: float

ci_high: float

verdict: str

to_dashboard()

Return type:

dict[str, Any]

__init__(level, n, hits, coverage, ci_low, ci_high, verdict)

class mmm_framework.diagnostics.coverage.RecoveryTargetStat(name, kind, truth, n_sims, levels, bias_z, z_spread, spread_p, rmse, mean_post_sd, verdict, flags=<factory>)

Bases: object

Per-target recovery-coverage verdict (a parameter or an estimand).

name: str

kind: str

truth: float

n_sims: int

levels: list[CoverageLevelStat]

bias_z: float

z_spread: float

spread_p: float

rmse: float

mean_post_sd: float

verdict: str

flags: list[str]

coverage_at(level)

Return type:

CoverageLevelStat | None

to_dashboard()

Return type:

dict[str, Any]

__init__(name, kind, truth, n_sims, levels, bias_z, z_spread, spread_p, rmse, mean_post_sd, verdict, flags=<factory>)

class mmm_framework.diagnostics.coverage.RecoveryCoverageResult(targets, n_sims_requested, n_sims_effective, n_failed_fits, truth_source, sampler, L, seed, levels, elapsed_s=0.0, caveats=<factory>)

Bases: object

Full recovery-coverage run: per-target stats + run metadata.

targets: list[RecoveryTargetStat]

n_sims_requested: int

n_sims_effective: int

n_failed_fits: int

truth_source: str

sampler: str

L: int

seed: int

levels: tuple[float, ...]

elapsed_s: float = 0.0

caveats: list[str]

property headline_level: float

property all_nominal: bool

True when no target’s headline-level interval demonstrably under-covers.

worst()

Return type:

RecoveryTargetStat | None

summary()

Return type:

str

to_dashboard()

Return type:

dict[str, Any]

__init__(targets, n_sims_requested, n_sims_effective, n_failed_fits, truth_source, sampler, L, seed, levels, elapsed_s=0.0, caveats=<factory>)

mmm_framework.diagnostics.coverage.jeffreys_interval(k, n, prob=0.95)

Jeffreys (Beta(1/2,1/2)) binomial interval for an empirical proportion.

The Monte-Carlo error bar on an estimated coverage: with n simulations and k hits, the coverage estimate is k/n with this uncertainty.

Return type:

tuple[float, float]

mmm_framework.diagnostics.coverage.coverage_from_ranks(int_ranks, L, levels=(0.5, 0.8, 0.9, 0.95))

Empirical central-interval coverage read directly off SBC ranks.

The truth lies inside the central level posterior interval exactly when its normalized rank u = (r + 0.5)/(L + 1) falls in [(1−level)/2, 1−(1−level)/2] — so an SBC run already contains every coverage number; this just states them in user language.

Return type:

list[CoverageLevelStat]

mmm_framework.diagnostics.coverage.coverage_from_draws(truth, draws_by_sim, levels=(0.5, 0.8, 0.9, 0.95))

Coverage of central equal-tailed percentile intervals across refits.

draws_by_sim[i] is the posterior sample for one simulated dataset; the truth is fixed. Matches the framework’s reported intervals (compute_hdi_bounds is percentile-based).

Return type:

list[CoverageLevelStat]

mmm_framework.diagnostics.coverage.recovery_diagnosis(truth, post_means, post_sds)

Decompose a coverage failure into bias vs. interval width.

With z_i = (mean_i − θ*)/sd_i across refits: under calibrated recovery the z_i are ≈ standard normal, so :rtype: dict[str, Any]

• bias_z = mean(z)·√n ≈ N(0,1) — large |bias_z| ⇒ the posterior location is systematically off (positive = estimates sit above truth);

• z_spread = sd(z) ≈ 1 — significantly > 1 ⇒ intervals too narrow (overconfident), < 1 ⇒ too wide (conservative). Significance via (n−1)·s² ~ χ²(n−1).

mmm_framework.diagnostics.coverage.build_recovery_result(truths, draws_by_target, *, kinds=None, levels=(0.5, 0.8, 0.9, 0.95), truth_source='fixed', sampler='generic', L=0, seed=0, n_sims_requested=0, n_failed_fits=0, elapsed_s=0.0, extra_caveats=())

Assemble a RecoveryCoverageResult from per-target refit draws.

Pure (no PyMC) — the fast unit tests drive this directly with analytic posteriors.

Return type:

RecoveryCoverageResult

mmm_framework.diagnostics.coverage.failure_mode_guide()

Markdown table of the ways nominal coverage fails and what to do.

Kept next to the tool so every surface (chat, Validation tab, docs) tells one story. The long-form version is technical-docs/coverage-diagnostics.md.

Return type:

str

mmm_framework.diagnostics.coverage.run_recovery_coverage(model, *, truth='posterior_mean', n_sims=24, levels=(0.5, 0.8, 0.9, 0.95), sampler='numpyro', L=200, tune=200, chains=2, seed=0, params=None, include_contributions=True, progress=None, refit=None, extra_caveats=())

Fixed-truth recovery coverage for a BayesianMMM.

Fixes every free parameter at θ* (pm.do), simulates n_sims datasets from the likelihood at θ*, refits the ORIGINAL model on each (pm.observe on the same graph — θ* and every posterior share one fixed standardized scale), and measures how often each central interval contains the truth — for scalar free parameters and per-channel total contributions.

EXPENSIVE: one refit per simulation, like SBC. Use a background job for thorough runs; n_sims < 30 only flags gross failures (the binomial CI is reported either way).

Answers: “IF the world matched my fitted model exactly, would my reported intervals cover?” Under-coverage here is a mechanical problem (approximate inference, sampler, priors); it deliberately cannot see real-world misspecification — see the module docstring.

Parameters:

• refit (Optional[Callable[[ndarray, int], Any]]) – Estimator injection. None (default) refits the PyMC model with pm.observe + NUTS, which is the only thing this used to be able to do. Supply a callable (y_sim, seed) -> idata to grade a non-PyMC estimator — the frequentist ridge + bootstrap path (#186) is why the hook exists. The simulate-from-the-model half is unchanged, so both paths are graded against the same θ* and the same central equal-tailed intervals. The callable may return anything carrying a .posterior mapping, or a plain {name: draws} mapping.

• extra_caveats (tuple[str, ...]) – Statements appended verbatim to the result’s caveats. Required for the frequentist path: this function auto-attaches an uncertainty caveat for advi/fullrank_advi only, so ridge shrinkage bias and conditional-on-selection intervals must be passed explicitly or they go unstated.

Return type:

RecoveryCoverageResult

Note

Parameters the injected estimator holds fixed (the frequentist path fixes adstock and saturation) have degenerate draws, and their coverage is not meaningful. Pass params=[...] to restrict the table to the quantities the estimator actually estimates.

Learning

Prior-vs-posterior learning diagnostics.

A posterior credible interval only tells you what you believe after seeing the data – it does not tell you how much of that belief came from the data versus the prior. When a prior is informative or sign-constrained, a posterior can look “conclusive” while encoding almost nothing the prior did not already assert.

The canonical example in this framework is the cannibalization cross-effect of MultivariateMMM, whose prior is psi = -HalfNormal(sigma) – structurally non-positive. Reporting that the posterior of psi is “entirely below zero” is then nearly vacuous: so is the prior. The honest question is whether the data moved or narrowed that parameter relative to its prior.

This module quantifies that with three complementary, sample-based diagnostics:

• contraction c = 1 - Var_post / Var_prior (Betancourt; Schad, Betancourt & Vasishth 2021). c -> 1 the posterior is far narrower than the prior (the data pinned it); c ~ 0 the posterior is as wide as the prior (data uninformative); c < 0 the posterior is wider than the prior – a red flag for prior-data conflict or poor sampling (we deliberately do not clip it). Contraction is a pure variance ratio, so it is robust even for bounded/half-normal priors.

• overlap – the prior-posterior overlap coefficient OVL = sum_i min(p_i, q_i) over shared histogram bins (probability-mass form, in [0, 1]). OVL ~ 1 the posterior is indistinguishable from the prior (nothing learned); OVL ~ 0 they barely overlap (strong learning, whether by narrowing or by shifting). Histogram binning – not a Gaussian KDE – is used on purpose so a hard prior edge (e.g. the HalfNormal boundary at 0) is not smeared.

• shift_z = (mean_post - mean_prior) / sd_prior – how far, in prior standard deviations, the posterior mean moved. Catches pure location learning that contraction alone misses (a posterior can shift without narrowing).

• contraction_robust – the same variance-ratio idea computed on the interquartile range instead of the standard deviation. The sd is tail-sensitive: heavy posterior tails (a genuine feature of the likelihood, or an artifact of chains that mixed slowly into a region the prior disfavored) inflate post_sd and drag contraction negative even when the posterior bulk narrowed. If contraction is negative but contraction_robust is not, the widening is tail-driven – corroborate with post_ess_bulk before reading it as prior-data conflict.

• post_ess_bulk – the posterior bulk effective sample size per parameter (when chain structure is available). A negative contraction with low ESS may say more about the sampler than the posterior.

Each parameter gets a heuristic verdict – "strong" / "moderate" / "weak" / "relocated" / "prior-dominated" – from contraction, overlap, and shift_z (the thresholds are tunable; treat them as conventions, not law). A "prior-dominated" verdict is the diagnostic the user is after: the posterior is essentially the prior. "relocated" is the opposite failure of a width-only reading: the posterior moved at least one prior-sd without narrowing – the evidence dominated the location (the prior was tight in the wrong place); the width reflects the likelihood’s flatness in the new region (or tail/mixing artifacts – see contraction_robust / post_ess_bulk), so do not read it as “the data taught nothing”.

The core parameter_learning() works on raw samples (arviz InferenceData, xarray Dataset, or plain dicts of arrays), so it is unit-testable without a fit. The model classes expose compute_parameter_learning(...) which draws prior samples and calls it for you.

mmm_framework.diagnostics.learning.parameter_learning(prior, posterior, var_names=None, *, bins=60, c_strong=0.5, c_weak=0.1, ovl_dominated=0.85, z_relocated=1.0)

How much did the data teach us about each parameter, beyond the prior?

Return type:

DataFrame

Parameters

prior, posterior:

Sample containers. Each may be an arviz InferenceData (the "prior" / "posterior" group is read), an xarray Dataset, or a dict mapping a parameter name to a sample array (sample axis first). Both must describe the same parameters; vector/matrix parameters are compared element-wise.

var_names:

Optional restriction. Names match either a flattened element ("beta[0]") or a base parameter ("beta" keeps every "beta[...]" element). None compares every parameter present in both containers.

bins:

Histogram resolution for the overlap coefficient.

c_strong, c_weak, ovl_dominated, z_relocated:

Heuristic verdict thresholds (conventions, not law). Checked in order: |shift_z| >= z_relocated and contraction < c_weak -> "relocated" (the posterior moved at least one prior-sd without narrowing – evidence dominated the location); then contraction < c_weak and overlap > ovl_dominated -> "prior-dominated"; contraction >= c_strong -> "strong"; contraction >= c_weak -> "moderate"; otherwise "weak".

Returns

pandas.DataFrame

One row per parameter, sorted by contraction ascending (so the least-learned, most prior-dominated parameters sort to the top). Columns: parameter, prior_mean, prior_sd, post_mean, post_sd, contraction, contraction_robust, overlap, shift_z, post_ess_bulk, verdict.

Notes

contraction is intentionally not clipped: a negative value (posterior wider than prior) is a genuine warning sign, but it is ambiguous on its own – it can mean prior-data conflict, a likelihood that is genuinely flat in the newly-favored region, or sd inflation from heavy tails left by chains that mixed slowly into that region. Disambiguate with the companions: a large |shift_z| says the evidence dominated the location (verdict "relocated", not a failure to learn); contraction_robust >= 0 alongside contraction < 0 says the widening is tail-driven; a low post_ess_bulk says the width estimate itself may be a sampling artifact – re-run with more tune/draws before interpreting.

Informativeness, not importance. High contraction / low overlap means the data was informative about the parameter – which includes confidently pinning it near zero. It is not a statement about effect size or sign. Always read the posterior post_mean / interval alongside the diagnostic: contraction tells you the data spoke; the posterior location tells you what it said. (A sign-constrained prior such as psi = -HalfNormal makes “the posterior is below zero” near-automatic; a high contraction to a value of ~0 means the data confidently found a negligible effect, not a confirmed one.)

mmm_framework.diagnostics.learning.plot_parameter_learning(learning, ax=None, *, top=None, metric='contraction', threshold=None)

Horizontal bar chart of per-parameter learning, colored by verdict.

metric is the column to plot ("contraction" by default; "overlap" also works). top keeps only the top least-learned parameters (the head of the sorted frame). threshold draws a reference line.

mmm_framework.diagnostics.learning.plot_prior_posterior_overlay(prior, posterior, parameter, ax=None, *, bins=60, transform=None)

Overlay the prior and posterior sample histograms for a single parameter.

The clearest way to see whether the data moved/narrowed a parameter relative to its prior. transform (callable) is applied to both sample sets before plotting – e.g. lambda x: -x to view the signed cannibalization effect from its psi_..._raw HalfNormal parameter.

Snapshot

Fit-time model-health snapshot: the JSON block a model_run carries so the UI can show whether to believe the fit — not just what it estimated.

Two layers, mirroring the docs’ honesty contract (green diagnostics validate computation, not causality):

• convergence — did the sampler work? Divergences, max R-hat, min bulk ESS (already computed by BayesianMMM.fit; recomputed from the trace as a fallback for models that don’t stamp them).

• learning — did the data inform the parameters, or are the posteriors re-stating the priors? Per-parameter prior→posterior contraction / overlap / shift verdicts from mmm_framework.diagnostics.parameter_learning(), reusing the prior group already attached to the trace at fit time (no extra sampling).

Like planning.history.compute_run_metrics, this runs kernel-side at fit time, must touch only the fitted model, and is best-effort: callers wrap it in try/except and a failure never fails a fit.

mmm_framework.diagnostics.snapshot.compute_fit_diagnostics(mmm, results=None, *, max_parameters=40)

JSON-safe model-health snapshot (schema v1) for a fitted model.

results is the MMMResults whose diagnostics dict already carries divergences/R-hat/ESS; pass None to recompute from the trace. Each layer is independently best-effort, so a learning failure still yields the convergence block (and vice versa).

Return type:

dict[str, Any]