How Many Parameters Did Your Data Pay For?
A fitted model tells you how many parameters it has. A second number tells you how many of them the data determines, and it comes from the same fit for the price of one extra compiled Hessian. Put two coefficients on the same regressor and the first number stays at two while the second comes back 1.000.
Curvature Has Two Owners
Everything a point fit knows about its own precision lives in the curvature of the log posterior at the mode. Steep along a direction means a narrow posterior. Flat along one, and the fit reports wherever the optimizer happened to stop. A Laplace approximation makes this literal by inverting that curvature and calling the result a covariance matrix, which is why a failing Hessian is a weak-identification flag delivered in seconds.
The question underneath that one is who paid for the curvature.
PyMC answers it almost for free. The log density decomposes as model.logp() == model.varlogp + model.datalogp, prior plus data, and second derivatives are linear, so the posterior curvature is exactly the sum of a prior part and a likelihood part:
The split
$$ H_{\text{post}} \;=\; H_{\text{prior}} \;+\; H_{\text{lik}} $$Each is the Hessian of the negative log density in PyMC's unconstrained space. One extra compiled Hessian of datalogp gives the second term, and the first follows by subtraction.
That identity is not an approximation, and the repo gates it. I ran four analytic models, three of them from tests/test_identification_diagnostics.py and a fourth with the prior narrowed to Normal(0, 1e-3), and took all three Hessians by hand: the largest absolute discrepancy between \( H_{\text{post}} \) and \( H_{\text{prior}} + H_{\text{lik}} \) was exactly 0.0 on every one of them. The test allows 1e-8 relative, which is the right thing to allow for a quantity that goes through two separate PyTensor compilations.
With the split in hand, the effective number of parameters is a trace:
$$ p_{\text{eff}} \;=\; \operatorname{tr}\!\left( H_{\text{post}}^{-1} H_{\text{lik}} \right) \;=\; \sum_j \frac{v_j^{\top} H_{\text{lik}} \, v_j}{\lambda_j} $$where \( (\lambda_j, v_j) \) are the eigenpairs of \( H_{\text{post}} \). The right-hand form is how diagnostics/identification.py computes it, and the reason is worth a sentence: summing per-direction shares needs no matrix inverse, and \( H_{\text{post}} \) is singular in exactly the cases you built the diagnostic for.
Here is what it returns on four models small enough to reason about completely. Each row is a real run of weak_identification_report(model, point) at the MAP point.
| Model | Parameters | Data paid for | Verdict |
|---|---|---|---|
| Two coefficients, two regressors, tight noise | 2 | 2.000 | ok |
| Two coefficients on the same regressor | 2 | 1.000 | weak |
| A \( [0,1] \) parameter entering the mean multiplied by zero | 2 | 0.9996 | weak |
One coefficient under a Normal(0, 1e-3) prior | 1 | 0.00021 | weak |
The second row is the one to sit with. Two coefficients, one regressor, and only their sum appears in the likelihood. Both parameters have perfectly good marginal posteriors. The model will print two numbers with two intervals, and a reader who takes them at face value will attribute revenue to each of them. The data bought one.
Eigenvalues in Whatever Units You Happened to Pick
Before any of that means anything, the curvature has to be put on a scale that does not depend on your data dictionary. Raw Hessian eigenvalues are not invariant to rescaling a parameter, so on a real MMM they rank parameters by whether spend was recorded in dollars or in thousands of dollars.
The demonstration takes four lines. Take the two-coefficient model above, multiply the regressor by 1000, and tighten that coefficient's prior by the same factor so the model is the same model in different units. The raw spectrum at the mode moves from \( \{18{,}423,\ 20{,}059\} \), a ratio of 1.089, to \( \{2.0 \times 10^{4},\ 1.8 \times 10^{10}\} \), a ratio of 924,000. Read raw eigenvalues and you have just discovered a catastrophic flat direction in a model that did not change.
Scale by the diagonal first and the problem disappears. The correlation-scaled curvature \( D^{-1/2} H D^{-1/2} \) with \( D = |\operatorname{diag}(H)| \) has eigenvalues \( \{0.9841,\ 1.0159\} \) in both parameterizations, a Belsley condition index of 1.0160 in both, and an effective parameter count of 1.999996 in both.
Why the saddle test survives the scaling
\( D \) is a positive diagonal matrix, so \( D^{-1/2} H D^{-1/2} \) is a congruence transformation of \( H \). By Sylvester's law of inertia, congruence preserves the number of positive, negative and zero eigenvalues. A negative eigenvalue after scaling is therefore still an unambiguous negative eigenvalue before it, which is what lets the same decomposition carry both the ridge test (small positive eigenvalue) and the saddle test (any negative one) without a second matrix.
Belsley's conventional threshold on the condition index \( \sqrt{\lambda_{\max} / \lambda_i} \) is 30. The shared-regressor model scores 18,548. Nothing about that is a borderline call. Its smallest scaled eigenvalue is 5.8e-9, and the direction carrying it loads on the two coefficients at \( (-0.707, +0.707) \): the difference of the two coefficients, exactly as constructed.
A Parameter Determined by Nothing in the Data
Two of the four things this decomposition separates are already familiar. A ridge is a parameter combination that trades off. It shows up as a small scaled eigenvalue, and the loadings on that eigenvector name which parameters are trading. An uninformed parameter has near-zero curvature of its own, nothing trades off against it, and it is what sends an optimizer walking out to the edge of a bounded parameter's support until a compiled Beta density divides by exactly zero.
The third case is the one that gets published.
Row three of the table is a two-parameter model where a Beta(1, 1) parameter called alpha enters the mean multiplied by zero. The likelihood is completely silent about it. Its posterior is not flat, because the Beta prior curves it. It is nowhere near a boundary, because the prior holds it in the middle. The ridge test finds nothing, the uninformed test finds nothing, and the report still comes back weak, because the attribution sees that the likelihood contributes 0% of alpha's curvature.
⚠️ The Laplace covariance is fine, and one of its intervals is the prior
On that model report.laplace_usable comes back True, and correctly: the smallest scaled eigenvalue is 1.0, there is no ridge, no uninformed parameter and no boundary hit. Laplace will invert a well-conditioned matrix and return two tidy intervals. One of them is the Beta(1, 1) prior, rendered as a finding. laplace_usable answers whether the covariance means anything as a covariance, which is a narrower question than whether every number in it came from the data.
The docstring for PriorDetermined puts the stakes plainly: the parameter is well determined, the optimizer is nowhere near trouble, “but every bit of that determination came from the prior. Reporting it as ‘fine’ would be reporting the prior back to the user as a finding.”
The trap in the per-parameter version
There is an obvious way to make this per-parameter instead of per-direction: take the ratio of diagonal entries, \( H_{\text{lik}}[i,i] / H_{\text{post}}[i,i] \), and call it the share of parameter \( i \)'s curvature the data supplied. That ratio is real and the module reports it. It can only be read in one direction.
On the shared-regressor model, both coefficients score 0.99999999 on the diagonal. The likelihood does curve in each coordinate individually, so each parameter looks like a triumph of measurement. Compute it along the flat direction instead and the eigenvector version returns an informed fraction of 0.0000. Only the sum is measured, and the diagonal has no way to see that, because seeing it requires looking at a combination of two coordinates at once.
So a small diagonal fraction is trustworthy and a large one is worthless. The asymmetry has a reason. Likelihood curvature is positive semi-definite, and a zero on the diagonal of a PSD matrix forces its entire row to zero, so a near-zero share really does mean the likelihood is flat in that coordinate. A near-one share means only that it curves there, which is compatible with the parameter being unidentified in company. tests/test_identification_diagnostics.py pins both halves, and the name of the test that pins the second is test_the_per_parameter_fraction_is_not_used_to_claim_identification.
The Same Question, Asked From the Draws
Curvature at a point is one view. The complementary one works from samples and asks how far the posterior moved from the prior it started at. diagnostics/learning.py computes it as a variance ratio:
One case in this framework makes the point concretely, and it ships. MultivariateMMM gives its cannibalization cross-effect the prior \( \psi = -\mathrm{HalfNormal}(\sigma) \), structurally non-positive, because a channel cannibalizing a sibling category cannot help it. With \( \sigma = 0.3 \) the prior puts 100.00% of its mass below zero and its own central 95% interval runs \( [-0.671,\ -0.010] \). That interval excludes zero, and the model has not seen any data yet.
Hand that prior a likelihood so weak it barely tilts anything, and the posterior comes back at \( [-0.659,\ -0.009] \) with 100.00% of its mass below zero. A report that says “the cannibalization effect is negative with 95% credibility” is repeating the prior. Contraction says so: 0.038, with a prior-posterior overlap coefficient of 0.992, and the verdict is prior-dominated.
What makes the module worth reading is that it refuses to hand you one number. Five come back per parameter, and the argument in the docstring is that any one of them alone will mislead you.
| Column | What it measures |
|---|---|
contraction | \( 1 - \operatorname{Var}_{\text{post}} / \operatorname{Var}_{\text{prior}} \), deliberately unclipped so a widening posterior stays visible |
overlap | \( \sum_i \min(p_i, q_i) \) over shared histogram bins, so a hard prior edge at zero is not smeared the way a Gaussian kernel would smear it |
shift_z | \( (\bar\theta_{\text{post}} - \bar\theta_{\text{prior}}) / \mathrm{sd}_{\text{prior}} \), location learning that a width-only reading misses |
contraction_robust | the same ratio on interquartile ranges, so heavy posterior tails cannot drag the headline number negative on their own |
post_ess_bulk | whether the width estimate is itself a sampling artifact |
The verdict rule checks shift_z before it checks width, and that ordering carries the module's sharpest idea. A parameter whose posterior moved at least one prior standard deviation without narrowing is labelled relocated rather than weak, because the evidence dominated the location. Its prior was tight in the wrong place.
The two verdicts are easy to confuse and easy to separate. Both of these use the same Normal(0, 0.10) prior on a standardized coefficient, and differ only in what the data says.
| Likelihood | Posterior | contraction | shift_z | overlap | Verdict |
|---|---|---|---|---|---|
| centre 0.05, standard error 1.00 | 0.000 ± 0.100 | 0.007 | 0.002 | 0.995 | prior-dominated |
| centre 1.50, standard error 0.35 | 0.113 ± 0.096 | 0.079 | 1.132 | 0.564 | relocated |
Both rows have a contraction under 0.1, and a width-only reading calls both of them “the data taught us nothing.” The second row is a prior and a likelihood in open disagreement. Its posterior sits at 0.113 while the data points at 1.50, the prior is winning by a wide margin, and the reported number is a compromise that neither party would have written down. Overlap tells the same story from another angle, 0.995 against 0.564.
Informativeness, not importance
A high contraction means the data was informative about the parameter, which includes confidently pinning it near zero. It says nothing about effect size or sign. Contraction tells you the data spoke; the posterior location tells you what it said. Both belong in the same sentence of a readout. reporting/evidence.py puts them there, and anything under a contraction of 0.10 gets the prior-dominated chip on the client-facing report.
This is also the check that closes a hole left by simulation-based calibration. As the SBC post covers, a broken computation whose posterior simply stays at the prior sails through parameter-by-parameter rank uniformity, because prior draws ranked against prior draws are uniform by construction. Contraction near zero is precisely that pathology viewed from the other side. The two checks are complementary, and running SBC without one leaves a known blind spot open. Both are workflow instruments in the sense Gelman, Vehtari, McElreath and their co-authors give the term, things you run around a fit rather than results you read out of one.
Before Any Fitting, Where Does the Prior Put the Elbow?
Saturation is not identified from the sales likelihood, which this blog has argued at length and which Jin et al. established for the Hill parameters a decade ago. When a parameter is unidentified, the prior is the answer, and the useful follow-up question is whether that prior is defensible.
It is answerable before a single draw. Media reaches the saturation transform already normalized by the channel's training maximum, so the elbow of the curve is directly a fraction of the largest spend you have ever recorded. Prior mass beyond 1.0 is mass on a region no observational data can reach, and it comes back as posterior essentially unchanged. saturation_prior_report draws the saturation parameter from its prior, converts it to an elbow, and reports mean(elbow > 1.0) per channel.
The conversion is closed form per family. Logistic saturation gives \( \ln 2 / \lambda \); Hill and Michaelis-Menten read the half-saturation parameter directly; tanh gives \( \operatorname{arctanh}(0.5) \cdot \kappa \). The root family has no asymptote and therefore no elbow, so it is absent from the table by design rather than by oversight.
What the number is for shows up most clearly in this framework's own history. Before 1.3 the default put sat_lam on an Exponential(0.5), which had never been reparameterized after the normalization landed.
Elbow prior on sat_lam | 90% interval (fraction of max spend) | Mass beyond max spend | Verdict |
|---|---|---|---|
Exponential(0.5) (pre-1.3) | 0.116 to 6.757 | 29.3% | unanchored |
LogNormal(log(ln2 / 0.5), 0.4) (current) | 0.259 to 0.965 | 4.2% | anchored |
Nearly a third of the old prior's mass sat past the last dollar anyone had spent, on a parameter the likelihood cannot move. The thresholds the module warns at are 5% for diffuse and 20% for unanchored, and the current default lands at 4.2% with a median elbow at half of maximum spend. Run the same report over the bundled seven-channel example, where every channel uses Hill saturation, and all seven come back anchored with a 90% elbow interval of 0.138 to 0.866 of maximum spend and no mass at all beyond it.
The check is honest about what it is worth. Its docstring says it reports whether the prior is “defensible” and warns in the same sentence that a prior can be perfectly anchored while the parameter stays unidentified. What makes a saturation parameter identifiable is dose spread, a rank condition on the spend design rather than a property of any estimator, and planning/identification refuses to claim identification below three in-support spend levels.
Running both
# illustrative
from mmm_framework.diagnostics.saturation import saturation_prior_report
before = saturation_prior_report(model) # pre-fit: is the prior defensible?
results = model.fit(method="map", identification_check=True)
report = results.diagnostics["identification"] # how many parameters did the data pay for?
The report also rides out automatically whenever a guarded MAP or Laplace fit has to retry inside its float64-safe box.
What the Number Cannot See
All of the curvature work is local. It is the shape of the log density at one point, and it answers what this fit resolves here. It cannot see a second mode across the parameter space, a better one further out, or a likelihood that is flat somewhere you did not land. For global questions the module points at NUTS, SBC and coverage. That is the right referral.
Three narrower caveats matter in practice, and the code states all of them.
datalogp is observedlogp + potentiallogp, so every pm.Potential in the graph lands on the likelihood side of the split. For a potential that carries an observation, such as the log-sum-exp mixture likelihood in the garden's latent-class model, that is exactly right. For a potential encoding a soft prior constraint it is exactly wrong, and such a model will read as better-informed than it is. If you write potentials, know which kind you wrote.
find_MAP locates its mode with jacobian=False while the Hessian is taken with jacobian=True, so the point being differentiated is not the exact mode of the density being differentiated. Flatness survives the distinction. The condition index does not survive it as a precise quantity, and should be read as an order of magnitude.
And the verdict thresholds are conventions. Belsley's 30, a contraction of 0.5 for “strong” and 0.1 for “weak,” an overlap of 0.85, a shift of one prior standard deviation: each is a defensible line drawn through a continuum, and the underlying numbers are in the frame for anyone who wants to draw a different one.
One framing worth borrowing comes from outside this repo. Kitagawa and Kuang open their 2025 paper on identification-aware MCMC by putting non-identifiability on the computation side of the ledger: so long as the posterior is proper, the inference is well posed, and the difficulty is in sampling the thing. MMM discourse usually runs those two together. Keeping them apart is what tells you that an effective parameter count of 1.0 out of 2 will survive any amount of target_accept, because there is nothing there for the sampler to fix. Their own remedy needs the observationally equivalent set written down in closed form, which an approximate trade-off between adstock and saturation does not have, so the method stays where it is. The distinction travels.
What a Check Owes You When It Cannot Run
Every diagnostic in this post is built on machinery that can fail. A Hessian may not compile. A graph may have no second derivative. What the module does then is return verdict='unknown' with a note naming the exception, on the stated principle that a diagnostic must never be the reason a fit fails. Knock out the likelihood split specifically and the report keeps the ridge it already found and sets effective_parameters to None, which is a test in its own right. Elsewhere, diagnostics/convergence.py returns None rather than True for a fit with no chains, on the same reasoning.
One module in the repo does not yet hold that line. validation/channel_diagnostics.py wraps its ArviZ summary call in a bare except, and on any failure returns every channel as rhat=1.0, ess_bulk=1000.0, converged=True. A check that reports a pass when it could not run is worse than one that raises, because the pass is indistinguishable from a real one downstream. The same class matches parameters to channels by substring, so a channel named tv claims the R-hat of tv_brand, and a channel matching no parameters at all inherits the optimistic seeds and reports as converged.
The fix is a small one and the pattern already exists twice in the same repo: a three-valued verdict, with None surfaced as “could not assess” rather than folded into “fine.” Everything else in that file is the useful part and deserves the upgrade. Its detect_collinear_clusters runs union-find over the channel correlation graph and returns connected components, the right unit of identification in a collinear MMM. The honest statement about two channels that always move together is a statement about the pair.
Give that pair one row. The clustering works on the design matrix and weak_identification_report works on the posterior curvature. The design side records np.linalg.cond of the correlation matrix and names 30 only in a docstring as a rule of thumb. Belsley's 30, the threshold identification.py actually applies, sits on \( \sqrt{\lambda_{\max} / \lambda_i} \), the square root of that quantity, so the two thirties are a square apart. A channel group the clustering flags and a flat direction that trips the report are one finding arriving twice whenever the collinearity in the design is what flattened the posterior. One row naming the group and reporting its combined effect beats two warnings and a table of individual numbers nobody should read separately.
Where This Goes Next
Report \( p_{\text{eff}} \) beside the parameter count in every readout that already carries a convergence chip. The number exists on any MAP or Laplace fit run with identification_check=True and it costs one compiled Hessian. Two counts side by side, with a plain sentence saying the second is what the data determines, is something a marketing director can act on in a way that a condition index is not.
Close the loop with the prior-widening refit that the confounding post asked for. That post wanted a measure of the prior's effect on the point estimate rather than on the width, and noted that contraction only covers the width. Half of it is already here: shift_z measures location movement from the draws, and relocated is the verdict when location moved and width did not. The other half is the sensitivity curve, refitting with the prior progressively loosened and watching where the coefficient walks.
Takeaways
model.logp() == model.varlogp + model.datalogpand Hessians add, so one extra compiled Hessian splits posterior curvature into the part the likelihood paid for and the part the prior did. \( \operatorname{tr}(H_{\text{post}}^{-1} H_{\text{lik}}) \) is 2.000 for two well-measured coefficients, 1.000 when two share one regressor, and 0.00021 under a prior tight enough to swamp the likelihood.- Correlation-scale the Hessian before you eigendecompose it. Rescaling one regressor by 1000 moved the raw eigenvalue ratio from 1.089 to 924,000 and left the scaled spectrum and the condition index unchanged to four decimals.
- A parameter can have ample posterior curvature, none of it the data's. That case leaves no ridge, no flat coordinate and no boundary hit, and a Laplace fit will report its prior as an interval.
- Read the per-parameter data share downward only. Two coefficients on one regressor both score 0.99999999 on the diagonal while the data determines one parameter between them.
- A posterior interval excluding zero is uninformative when the prior's interval also excludes zero. Under \( \psi = -\mathrm{HalfNormal}(0.3) \) that is true by construction, and contraction of 0.038 with overlap of 0.992 is what says so.
- 29.3% of the pre-1.3 saturation prior's elbow mass sat beyond maximum observed spend, against 4.2% for the current default. That is measurable before fitting, and the warning threshold is 20%.
References
- Gelman, A., Vehtari, A., McElreath, R., et al. (2026). Bayesian Workflow. CRC Press. ISBN 978-0-367-49014-0.
- Schad, D. J., Betancourt, M., & Vasishth, S. (2021). Toward a Principled Bayesian Workflow in Cognitive Science. Psychological Methods, 26(1).
- Jin, Y., Wang, Y., Sun, Y., Chan, D., & Koehler, J. (2017). Bayesian Methods for Media Mix Modeling with Carryover and Shape Effects. Google Research.
- Dew, R., Padilla, N., & Shchetkina, A. (2024). Your MMM is Broken: Identification of Nonlinear and Time-varying Effects in Marketing Mix Models. arXiv:2408.07678 [econ.EM].
- Kitagawa, T., & Kuang, Y. (2025). Identification-aware Markov chain Monte Carlo. arXiv:2511.12847.
- Source:
src/mmm_framework/diagnostics/identification.py,diagnostics/learning.py,diagnostics/saturation.py,validation/channel_diagnostics.py,tests/test_identification_diagnostics.py.