Skip to content

A negative count is scored as log p = 0 (a perfect fit) and counted in n for AIC/BIC/LOO #523

Description

@wshlavacek

Summary

The count families guard against a negative observation by having it contribute nothing:

def data_fit(self, prediction, observation, noise, extra=None):
    if observation < 0:
        return 0

That guard is right — a negative value has no probability under a negative binomial, and real
surveillance data contains them (a state revising its cumulative total downward produces a
negative daily increment). But "contribute nothing to the cost" becomes something different
one layer up. NoiseModel.log_density is -(data_fit + normalizers + constant), and the count
family is self-normalizing, so:

>>> nb = NegBinomial(location=MEAN); r = 9.06252
>>> nb.log_density(3000.0, -5948.0, r)      # negative count
-0.0
>>> nb.log_density(3000.0,     0.0, r)
-52.61
>>> nb.log_density(3000.0,  3000.0, r)      # a perfect prediction
-7.83

A negative observation is assigned log p = 0, i.e. probability 1 — a better per-point
density than any real observation can achieve, including one the model predicts exactly.

_pointwise_suffix skips a point only on np.isnan(observation), so these points are also
counted in n:

for col_name in sorted(compare_cols):
    observation = exp_data.data[rownum, exp_data.cols[col_name]]
    if np.isnan(observation):
        continue

So they enter likelihood_information_criteria (n = len(values)) and the ids/values arrays
that back LOO/WAIC (ADR-0056).

Impact

For AIC/BIC the effect is small in the case that surfaced it. Fitting NYT MSA daily case counts
(Mallela et al. 2024) with noise_model = neg_bin, ..., location = mean:

MSA rows negative counts
nyc 649 1 (t=460, −5948)
houston 649 4
phoenix 649 2
dallas 649 0

The nyc fit reports AIC=10475.6 BIC=10484.5 AICc=10475.61733 (k=2, n=649, lnL=-5235.8).
Correcting n to 648 moves BIC by k·ln(649/648) ≈ 0.003 — negligible.

Two reasons I think it is still worth fixing:

  1. It scales badly and silently. Each negative point contributes the maximum possible
    per-point log-density. A series with many downward revisions — not unusual in surveillance
    data — looks increasingly well-fit for a reason that has nothing to do with the model, and
    nothing in the output says so.
  2. LOO/WAIC are more exposed than AIC. They consume the per-point densities directly, so a
    handful of points pinned at log p = 0 distorts the pointwise distribution (and the Pareto-k
    diagnostics computed from it) more than it distorts a summed criterion.

Suggested fix

Treat an out-of-domain observation the way NaN is already treated — as not a scored point
rather than as a perfectly-scored one:

  • skip it in _pointwise_suffix (and the parallel walk in aligned_prediction_data), so it is
    excluded from n, from LOO/WAIC, and from the ids axis;
  • warn once per observable with the count excluded, e.g.
    neg_bin: excluded 4 measurement(s) of 'fDCs_Cum' with a negative count (out of the count domain) — today a user has no way to learn this happened;
  • leave data_fit / the derivative guards exactly as they are: not contributing to the cost
    is the correct behavior for the fit itself, and this issue is only about what happens
    downstream of that.

One wrinkle worth deciding deliberately: evaluate_pointwise documents that "the emitted
observation set is fixed by the experimental data (a point is skipped only on a NaN
observation, never on anything draw-dependent), so every draw yields the same ids in the same
order". A negative-observation skip preserves that property — the observation is data, not
draw-dependent — but the docstring should be updated to name the second condition.

Environment

PyBNF 7913d037 (main), Python 3.12.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions