Skip to content

Ascertainment Models for Related Observation Processes

Code 
import jax.numpy as jnp
import numpyro.distributions as dist

from pyrenew.ascertainment import JointAscertainment, RatioLinkedAscertainment
from pyrenew.deterministic import DeterministicPMF, DeterministicVariable
from pyrenew.model import PyrenewBuilder
from pyrenew.observation import NegativeBinomialNoise, PopulationCounts
from pyrenew.randomvariable import DistributionalVariable
from pyrenew.time import MMWR_WEEK

Overview

An ascertainment rate is the rate at which an infection results in an observed event. For hospital admissions, this is often called the infection-hospitalization rate (IHR). For emergency department visits, it can be called the infection-to-ED-visit rate (IEDR).

For a count observation process, the ascertainment rate scales delayed latent infections to obtain the expected number of observed events:

\[ \mu_s(t) = \alpha_s \sum_{d=0}^{D_s} I(t-d)\,\pi_{s,d}, \]

where \(s\) identifies the signal, \(\alpha_s\) is its ascertainment rate, and \(\pi_{s,d}\) is its delay distribution.

If a model treats the ascertainment rates for different signals as unrelated, each observation process can receive its own RandomVariable directly. An ascertainment component is needed when the model represents a relationship between the ascertainment rates of multiple observation processes.

PyRenew provides two such components:

  • JointAscertainment specifies a joint prior for two or more rates.
  • RatioLinkedAscertainment specifies one rate as a multiple of another rate.

The choice between them is part of the statistical model, not only an API choice. Each component expresses a different assumption about how information is shared across signals.

Independent ascertainment rates

An ascertainment component is not required when the signal rates are modeled independently. Pass each rate directly to its observation process:

hospital_obs = PopulationCounts(
    name="hospital",
    ascertainment_rate_rv=DistributionalVariable("ihr", dist.Beta(2, 198)),
    delay_distribution_rv=hosp_delay_rv,
    noise=hosp_noise,
)

ed_obs = PopulationCounts(
    name="ed_visits",
    ascertainment_rate_rv=DistributionalVariable("iedr", dist.Beta(2, 98)),
    delay_distribution_rv=ed_delay_rv,
    noise=ed_noise,
)

This specification does not encode prior dependence between IHR and IEDR; each observation process has its own ascertainment-rate prior.

Ascertainment components and the builder

An ascertainment component is a model-level component shared by multiple observation processes. Define it and register it with PyrenewBuilder’s add_ascertainment() method before constructing those observation processes. Each observation process then sets its ascertainment_rate_rv using the component’s for_signal() method.

builder = PyrenewBuilder()
builder.configure_latent(...)

ascertainment = JointAscertainment(...)
builder.add_ascertainment(ascertainment)

hospital_obs = PopulationCounts(
    name="hospital",
    ascertainment_rate_rv=ascertainment.for_signal("hospital"),
    ...
)
builder.add_observation(hospital_obs)

ed_obs = PopulationCounts(
    name="ed_visits",
    ascertainment_rate_rv=ascertainment.for_signal("ed_visits"),
    ...
)
builder.add_observation(ed_obs)

During model execution, the registered component supplies the signal-specific rates to the observation processes that reference it. The signal name passed to for_signal() must match a name declared by the ascertainment component. It does not have to match the observation-process name, but matching names usually makes the model specification and output easier to interpret.

Choosing an ascertainment model

JointAscertainment and RatioLinkedAscertainment both relate signal rates, but they represent the relationship differently.

Model Parameters with priors Relationship represented Probability bounds
JointAscertainment A vector of logit-scale rates Prior covariance controls how the rates vary together Every rate is transformed to (0, 1)
RatioLinkedAscertainment A base rate and a relative multiplier The linked rate is the base rate multiplied by the ratio The derived linked rate is not automatically constrained to (0, 1)

Use JointAscertainment when prior dependence between the rates is most naturally described through a covariance structure. Use RatioLinkedAscertainment when the relative quantity, such as IHR / IEDR, is itself meaningful and can be assigned an interpretable prior.

Joint ascertainment

JointAscertainment places a multivariate normal prior on the logit-transformed ascertainment rates:

\[ \boldsymbol{\eta} \sim \operatorname{MultivariateNormal} \left( \operatorname{logit}(\boldsymbol{\alpha}_0), \boldsymbol{\Sigma} \right), \qquad \boldsymbol{\alpha} = \operatorname{logit}^{-1}(\boldsymbol{\eta}). \]

The baseline rates \(\boldsymbol{\alpha}_0\) set the prior location. The covariance matrix \(\boldsymbol{\Sigma}\) controls the scale of variation in each logit-transformed rate and how strongly the rates vary together under the prior. It describes prior dependence, not a causal relationship between the signals.

joint_ascertainment = JointAscertainment(
    name="he_ascertainment",
    signals=("hospital", "ed_visits"),
    baseline_rates=jnp.array([0.01, 0.02]),
    scale_tril=jnp.array(
        [
            [0.7, 0.0],
            [0.35, 0.606],
        ]
    ),
)

The order of signals determines how entries in baseline_rates and the covariance parameter are interpreted. The inverse-logit transformation guarantees that every sampled rate lies between zero and one.

During model execution, this component creates one vector sample site, he_ascertainment_eta, and signal-specific deterministic sites such as he_ascertainment_hospital and he_ascertainment_ed_visits.

Ratio-linked ascertainment

RatioLinkedAscertainment samples a base rate and a ratio of the linked rate to the base rate:

\[ \alpha_{\mathrm{linked}} = \alpha_{\mathrm{base}}\rho, \qquad \rho = \frac{\alpha_{\mathrm{linked}}}{\alpha_{\mathrm{base}}}. \]

For a hospital and ED model, the IEDR can be the base rate and the ratio can represent IHR / IEDR. A ratio below one means that hospitalization is less likely than an ED visit for a given infection; a ratio above one means that it is more likely.

ratio_linked_ascertainment = RatioLinkedAscertainment(
    name="he_ascertainment",
    base_signal="ed_visits",
    linked_signal="hospital",
    base_rate_rv=DistributionalVariable("iedr", dist.Beta(1, 100)),
    ratio_rv=DistributionalVariable(
        "ihr_rel_iedr",
        dist.LogNormal(jnp.log(0.5), 0.35),
    ),
)

Note: These prior parameters are illustrative; they should be chosen using domain knowledge appropriate to the pathogen, population, observation definitions, and study period.

This model does not specify a covariance between the two signal rates directly. However, the signal rates are not independent: the linked rate contains the sampled base rate. The priors are assigned to the base rate and the ratio, and together they induce a distribution for the linked rate.

The ratio may be any positive value. Unlike JointAscertainment, this component does not automatically constrain the derived linked rate to the unit interval. If both rates should be \(<1\) (for instance because both are probabilities), make sure to choose priors that make \(\alpha_{\mathrm{base}}\rho > 1\) negligibly likely. The component does not clip an invalid linked rate, because clipping would change the specified statistical model.

Building a model an ascertainment component

The complete builder pattern is the same for either ascertainment component - here we use a RatioLinkedAscertainment but the sequence of definitions is the same for JointAscertainment.

builder = PyrenewBuilder()
builder.configure_latent(
    PopulationInfections,
    gen_int_rv=gen_int_rv,
    I0_rv=I0_rv,
    log_rt_time_0_rv=log_rt_time_0_rv,
    single_rt_process=rt_process,
)

ascertainment = RatioLinkedAscertainment(
    name="he_ascertainment",
    base_signal="ed_visits",
    linked_signal="hospital",
    base_rate_rv=DistributionalVariable("iedr", dist.Beta(1, 100)),
    ratio_rv=DistributionalVariable(
        "ihr_rel_iedr",
        dist.LogNormal(jnp.log(0.5), 0.35),
    ),
)
builder.add_ascertainment(ascertainment)

hospital_obs = PopulationCounts(
    name="hospital",
    ascertainment_rate_rv=ascertainment.for_signal("hospital"),
    delay_distribution_rv=hosp_delay_rv,
    noise=NegativeBinomialNoise(hosp_concentration_rv),
    aggregation="weekly",
    reporting_schedule="regular",
    start_dow=MMWR_WEEK,
)
builder.add_observation(hospital_obs)

ed_obs = PopulationCounts(
    name="ed_visits",
    ascertainment_rate_rv=ascertainment.for_signal("ed_visits"),
    delay_distribution_rv=ed_delay_rv,
    noise=NegativeBinomialNoise(ed_concentration_rv),
    day_of_week_rv=ed_day_of_week_rv,
)
builder.add_observation(ed_obs)

model = builder.build()

Shared assumptions and limitations

Both ascertainment components produce one scalar ascertainment rate for each signal over the entire model period. Consequently, neither component allows the rates or their implied relationship to change over time. Whether this assumption is reasonable should be evaluated in the context of the observation processes and study period.

The choice of ascertainment model should consider:

  • whether the rates should be related at all
  • whether a covariance structure or a relative multiplier is easier to justify
  • whether the priors keep all rates, including derived linked rates, within scientifically plausible ranges
  • whether treating ascertainment rates as constant over the model period is scientifically reasonable

Ascertainment rates are generally difficult to distinguish from the scale of latent infections using surveillance counts alone. Allowing the rates to vary over time introduces additional confounding unless external information or strong structural assumptions identify that variation. The built-in ascertainment components therefore use scalar rates over the model period. When this assumption is questionable, shorter modeling windows and sensitivity analyses may be more appropriate than introducing a weakly identified time-varying process.