Estimating Incidence Rates from Cross-Sectional Serosurveys
2026-08-10
1 Infectious disease epidemiology
1.1 Salmonella enterica (“typhoid fever”)
- ~21.7 million symptomatic cases/year
- ~217,000 deaths/year (~1% case fatality with treatment; up to 20% untreated)
- Highest burden: ages 5–19
- Symptoms: sustained high fever, severe headache, abdominal pain, malaise; rose-spot rash (~30% of patients); constipation (early) or diarrhea (later)
1.2 Shigella (“dysentery”)
- ~100 million infections annually
- ~100,000 deaths annually (~1% case fatality; higher in malnourished children)
- Most cases and deaths: children under 5 in low- and middle-income countries (LMICs)
- Symptoms: bloody diarrhea (dysentery), abdominal cramps, fever, tenesmus (painful urge to defecate)
1.3 Orientia tsutsugamushi (“scrub typhus”)
- ~1,000,000 infections/year
- ~10,000 deaths/year (~1% case fatality with treatment; up to 30% untreated)
- Symptoms: high fever, severe headache, myalgia; eschar (pathognomonic black scab at mite-bite site); macular rash; cough; gastrointestinal symptoms; hemorrhage in severe cases
2 Why estimate incidence from serosurveys?
2.1 The burden-data gap
Typhoid conjugate vaccines (TCVs) are 80-90% effective and WHO-recommended since 2018, but few countries have adopted them into routine immunization.
A major barrier is the lack of burden data: most low- and middle-income countries lack robust typhoid surveillance and have little or no incidence data.
These gaps block applications for vaccine funding, widening equity gaps in access to effective vaccines.
Three large randomized trials showed that a Vi tetanus-toxoid typhoid conjugate vaccine is 80-90% effective at preventing symptomatic typhoid, and the WHO has recommended these vaccines since 2018. Yet as of 2023 only a handful of countries had introduced them into routine immunization. The binding constraint is data: without incidence estimates, countries cannot make the case for vaccine funding.
“Data inequality is our biggest challenge moving forward” — Kathy Neuzil, 2023
2.2 Seroepidemiology can fill the gap
The Seroepidemiology and Environmental Surveillance for Enteric Fever (SEES) study collected cross-sectional serosurveys across multiple countries to estimate typhoid incidence from antibody data (Aiemjoy et al. 2022).
This lecture describes the methodology behind that kind of estimate: how a single cross-sectional serosurvey, combined with a model of antibody dynamics, can recover an incidence rate.
2.3 Defining incidence
2.4 Mathematical definition of incidence
More precisely, the incidence rate at time \(t\) is the rate of new infections per person at risk:
\[\lambda_t = \frac{1}{N(t)}\,\frac{d}{dt}\,\mathbb{E}[C(t)]\]
where \(C(t)\) is the cumulative number of infections and \(N(t)\) the number of individuals at risk at time \(t\).
2.5 Scale of incidence rates
In both definitions, the units for an incidence rate are “# new infections per # persons at risk per time duration”; for example, “new infections per 1000 persons per year”.
For convenience, we can rescale the incidence rate to make it easier to understand; for example, we might express incidence as “# new infections per 1000 persons per year” or “# new infections per 100,000 persons per day”, etc.
2.6 Incidence from an individual’s perspective
From the perspective of an individual in the population:
the incidence rate (at a given time point \(t\)) is the instantaneous probability density of becoming infected at that time point, given that they are at risk at that time point.
That is, the incidence rate is a hazard rate.
Notation: let’s use \(\lambda_{t}\) to denote the incidence rate at time \(t\).
3 Study designs for estimating incidence rates
3.1 Longitudinal cohort studies
Incidence rates can be estimated from longitudinal cohort studies, but cohort studies are:
- costly to conduct
- slow to produce results,
- vulnerable to selection and censoring (drop-out) biases
3.2 Clinical case data
Incidence rates can also be estimated from clinical case rates, but clinical case rates undercount:
- asymptomatic cases
- symptomatic cases who don’t receive clinical care.
3.3 Cross-sectional serosurveys
Alternatively, incidence rates can be estimated from biomarker data collected through a single cross-sectional survey, combined with a longitudinal model of how those biomarkers respond to infection over time. Compared with cohort studies and clinical case rates, cross-sectional incidence estimation can produce estimates that are:
- accurate
- timely
- cost-efficient
See Hay et al. (2024) for an overview.
4 Estimating incidence from cross-sectional serosurveys
4.1 Goal
Easily and reproducibly translate quantitative antibody responses at the population level into meaningful and accurate epidemiological measures of infection burden.
4.2 Antibody responses are complicated
Using antibody levels to recover infection times is hard, because antibody responses:
- decay over time after infection
- vary from individual to individual (age, immune function, prior infections, vaccination)
- vary from measurement to measurement (assay noise)
- can cross-react with antibodies from other exposures
Each of these complications has to be handled somewhere in the model: waning, between-person heterogeneity, measurement error, and cross-reactivity. The rest of this section builds that model up piece by piece.
4.3 Cross-sectional antibody surveys
Typically, it is difficult to estimate changes from a single time point. However, we can sometimes make assumptions that allow us to do so. In particular, if we assume that the incidence rate is constant over time, then we can estimate incidence from a single cross-sectional survey.
We will need two pieces of notation to formalize this process.
We recruit participants from the population of interest.
For each survey participant, we measure antibody levels \((Y)\) for the disease of interest
Each participant was most recently infected at some time \((T)\) prior to when we measured their antibodies.
- If a participant has never been infected since birth, then \(T\) is undefined.
- \(T\) is a latent, unobserved variable.
- We don’t directly observe \(T\); we only observe \(Y\), which we hope tells us something about \(T\) and \(\lambda\).
4.4 Modeling assumptions
(For diseases like typhoid, the no-immunity assumption may not hold exactly, but hopefully approximately; modeling the effects of re-exposure during an active infection is on our to-do list).
4.5 Time since infection and incidence
This is a time-to-event model, looking backwards in time from the survey date (when the blood sample was collected).
The probability that an individual was last infected \(t\) days ago, \(p(T=t)\), is equal to the probability of being infected at time \(t\) (i.e., the incidence rate at time \(t\), \(\lambda\)) times the probability of not being infected after time \(t\), which turns out to be \(\operatorname{exp}\mathopen{}\left\{-\lambda t\right\}\mathclose{}\).
The distribution of \(T\) is truncated by the patient’s birth date; the probability that they have never been infected is \(\operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}\), where \(a\) is the patient’s age at the time of the survey.
4.6 Maximum likelihood estimation
Several approaches could estimate \(\lambda\) from the data. serocalculator uses maximum likelihood: it reports the value of \(\lambda\) that makes the observed data most probable. Stating that precisely takes three definitions.
The log-likelihood is maximized at the same \(\lambda\) as the likelihood, and is easier to work with because the logarithm turns products into sums (the logarithm-of-a-product theorem in Regression Models for Epidemiology).
4.7 Likelihood of latent infection times
If the infection times \(T\) were observed directly, these would all take a simple closed form:
If we could observe \(T\), then we could estimate \(\lambda\) using a typical maximum likelihood approach.
. . .
Starting with the likelihood (Definition 3), which multiplies participants’ contributions because we treat them as independent:
\[ \begin{aligned} \mathscr{L}^*(\lambda) &= \prod_{i=1}^n \operatorname{p}(T=t_i \mid \lambda) \\&= \prod_{i=1}^n \lambda \operatorname{exp}\mathopen{}\left\{-\lambda t_i\right\}\mathclose{} \end{aligned} \]
. . .
Taking the logarithm of the likelihood:
\[ \begin{aligned} \ell^*(\lambda) &= \operatorname{log}\mathopen{}\left\{\mathscr{L}^*(\lambda)\right\}\mathclose{} \\&= \sum_{i=1}^n \mathopen{}\left(\operatorname{log}\mathopen{}\left\{\lambda\right\}\mathclose{} -\lambda t_i\right)\mathclose{} \end{aligned} \]
. . .
Taking the derivative of that log-likelihood to find the score function (Definition 5):
\[\ell^{*'}(\lambda) = \sum_{i=1}^n \mathopen{}\left(\mathopen{}\left(\lambda\right)^{-1}\mathclose{} - t_i\right)\mathclose{}\]
. . .
Setting the score function equal to 0 to find the score equation, and solving the score equation for \(\lambda\) to find the maximum likelihood estimate:
4.8 Example log-likelihood curves
Here’s what that would look like:
Code
library(serodynamics)
library(serocalculator)
library(dplyr)
antibodies <- c("HlyE_IgA", "HlyE_IgG")
set.seed(1)
sim_case_data <-
serocalculator::typhoid_curves_nostrat_100 |>
sim_case_data(n = 5,
antigen_isos = antibodies,
max_n_obs = 20, followup_interval = 14)
t1 <- sim_case_data$timeindays
loglik0 <- function(lambda) {
sum(
dexp(t1, rate = lambda, log = TRUE)
)
}
loglik1 <- Vectorize(loglik0, vectorize.args = "lambda")
library(ggplot2)
ggplot() +
geom_function(fun = loglik1) +
xlim(0, .1) +
theme_bw() +
xlab("lambda") +
ylab("log-likelihood")
4.9 Standard errors
How precisely the data pin down \(\hat\lambda\) depends on how sharply the log-likelihood peaks at its maximum, which takes two more definitions to make precise.
The standard error is the square root of that variance, so:
more curvature \(\rightarrow\) a sharper likelihood peak \(\rightarrow\) more observed information \(\rightarrow\) smaller standard errors.
4.10 Hessian for the exponential model
4.11 Likelihood of observed data
Unfortunately, we don’t observe infection times \(T\); we only observe antibody levels \({Y}\). So things get a little more complicated.
In short, we are hoping that we can estimate \(T\) (time since last infection) from \(Y\) (current antibody levels). If we could do that, then we could plug in our estimates \(\hat t_i\) into that likelihood above, and estimate \(\lambda\) as previously.
We’re actually going to do something a little more nuanced; instead of just using one value for \(\hat t\), we are going to consider all possible values of \(t\) for each individual.
We need to link the data we actually observed to the incidence rate.
The likelihood of an individual’s observed data, \(\operatorname{p}(Y=y)\), can be expressed as an integral over the joint likelihood of \(Y\) and \(T\) (using the Law of Total Probability):
- \[\operatorname{p}(Y=y) = \int_t \operatorname{p}(Y=y,T=t)dt\]
Further, we can express the joint probability \(p(Y=y,T=t)\) as the product of \(p(T=t)\) and \(p(Y=y|T=t)\) the “antibody response curve after infection”. That is:
- \[\operatorname{p}(Y=y,T=t) = \operatorname{p}(Y=y|T=t) \operatorname{p}(T=t)\]
4.12 Antibody response curves
Code
sim_case_data |>
autoplot(alpha = .5)
Code
case_data <-
serodynamics_example(
"SEES_Case_Nepal_ForSeroKinetics_02-13-2025.csv"
) |>
readr::read_csv() |>
dplyr::mutate(
.by = person_id,
visit_num = dplyr::row_number()
) |>
as_case_data(
id_var = "person_id",
biomarker_var = "antigen_iso",
value_var = "result",
time_in_days = "dayssincefeveronset"
)
most_obs <-
case_data |>
count(id) |>
arrange(desc(n)) |>
head(10)
case_data |>
semi_join(most_obs, by = "id") |>
autoplot(alpha = .5, log_x = FALSE)
4.13 The per-person likelihood
Two pieces here were not yet defined:
- \(p(T=\text{NA}) = \operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}\) is already known — the same never-infected probability from Proposition 1.
- \(p(Y=y\mid T=\text{NA})\), the antibody response of a never-infected subject, needs its own definition. Without any noise, a never-infected subject’s true antibody level is exactly zero, so this is a point mass at \(y=0\) (a Dirac delta): \(p(Y=y\mid T=\text{NA}) = \delta(y)\). Once measurement noise and biological (cross-reactivity) noise are introduced below, this point mass is handled differently by each — see “Noise and never-infected subjects”.
4.14 The full-sample likelihood
Independence is what licenses the product; the latent-time likelihood \(\mathscr{L}^*\) above multiplies participants’ contributions for the same reason. It is an assumption about the sampling design, separate from the assumptions about disease dynamics in Definition 2.
If we know \(\operatorname{p}(Y=y|T=t)\), then we can maximize \(\mathscr{L}(\lambda)\) over \(\lambda\) to find the “maximum likelihood estimate” (MLE) of \(\lambda\), denoted \(\hat\lambda\).
4.15 Finding the MLE numerically
Taking logs turns that product into a sum:
\[ \begin{aligned} \ell(\lambda) &= \operatorname{log}\mathopen{}\left\{\prod_{i=1}^n \mathscr{L}_i(\lambda)\right\}\mathclose{} \\&= \sum_{i=1}^n \ell_i(\lambda) \end{aligned} \]
Each \(\ell_i\) is the log of an integral (Proposition 3), so the derivative of this sum doesn’t come out cleanly, and we use a numerical method (specifically, a Newton-type algorithm, implemented by stats::nlm()) to find the MLE. That algorithm also returns a numerical Hessian (Definition 6) at the maximum, and Theorem 1 turns that Hessian into a standard error.
4.16 Cluster-robust standard errors for clustered sampling designs
In many survey designs, observations are clustered (e.g., multiple individuals from the same household, school, or geographic area). Observations within the same cluster are often more similar to each other than to observations from different clusters, which violates the independence assumption of Proposition 4. The estimate \(\hat\lambda\) itself remains usable, but the inverse-information standard error of Theorem 1 no longer applies to it.
4.16.1 Why clustering matters
When observations are clustered:
- Individuals within the same cluster share common exposures or characteristics
- Standard errors that ignore clustering are usually too small (anti-conservative)
- Confidence intervals are then too narrow
- p-values are then too optimistic
Usually rather than always: the direction follows the sign of the within-cluster correlation, which is positive in most survey designs but need not be.
4.16.2 Cluster-robust variance estimation
To account for within-cluster correlation, serocalculator implements the sandwich estimator — also known as the Huber-White robust variance estimator (Huber 1967; White 1980), generalized to clustered/grouped data by Liang and Zeger (1986).
The sign convention does not matter here: \(H\) appears twice, so substituting the observed information \(I(\hat\lambda) = -H\) (Definition 7) leaves \(V_{\text{robust}}\) unchanged.
The sandwich generalizes Theorem 1 rather than replacing it. When the model is correctly specified and observations really are independent, \(B\) estimates the same information as \(-H\) does, and \(H^{-1} B H^{-1}\) collapses to \(\mathopen{}\left[I(\hat\lambda)\right]\mathclose{}^{-1}\) — the ordinary standard error. Clustering is what drives the two apart, and the gap between them is what the correction recovers.
4.16.3 Multi-way clustering
When cluster_var names more than one variable, serocalculator does not collapse them into a single interaction. Instead it applies the inclusion-exclusion identity of Cameron, Gelbach, and Miller (2011) (Cameron et al. 2011) for multi-way cluster-robust variance estimation.
Computing every subset costs \(2^p - 1\) one-way variance evaluations, so this is intended for the small \(p\) (2-3) typical of nested survey designs.
4.16.4 Small-sample correction
Two adjustments refine the one-way variance terms \(V_S\) that feed Definition 9.
4.16.5 Safeguards against invalid variance
The inclusion-exclusion sum’s occasional negativity is a documented property of the Cameron-Gelbach-Miller estimator, not a bug.
4.16.6 Implementation in serocalculator
Users can specify clustering using the cluster_var parameter:
Code
# Single-level clustering (e.g., by household)
est <- est_seroincidence(
pop_data = data,
cluster_var = "household_id",
...
)
# Multi-level clustering (e.g., schools within districts)
est <- est_seroincidence(
pop_data = data,
cluster_var = c("district_id", "school_id"),
...
)When cluster-robust standard errors are used, the summary() output indicates this with se_type = "cluster-robust". The multi-way behavior above is itself controlled through summary():
Code
summary(
est,
small_sample = "CR1", # the default; "none" disables the CR1 correction
floor_to_standard = FALSE, # the default; TRUE additionally floors at V_standard
debug_cluster = TRUE # prints the V_S/V_raw/V_final decomposition
)Internally the optimization is parameterized by \(\log\lambda\) rather than \(\lambda\), which keeps the estimate positive, so the sandwich above is evaluated on that scale and summary() converts the result back with the delta method (\(\mathop{\widehat{\operatorname{SE}}}\nolimits\mathopen{}\left(\hat\lambda\right)\mathclose{} = \hat\lambda \cdot \mathop{\widehat{\operatorname{SE}}}\nolimits\mathopen{}\left(\log\hat\lambda\right)\mathclose{}\)).
4.16.7 Effect on results
-
Point estimates (incidence rates) are unchanged:
cluster_varaffects only the variance, not the likelihood being maximized - Standard errors usually increase, reflecting within-cluster correlation
- Confidence intervals widen correspondingly, reflecting the reduced effective sample size
5 Modeling the seroresponse kinetics curve
Now, we need a model for the antibody response to infection, \(\operatorname{p}(Y=y|T=t)\). The current version of the serocalculator package uses a two-phase model for the shape of the seroresponse (Teunis et al. 2016).
5.1 Model for active infection period
The first phase of the model represents the active infection period, and uses a simplified Lotka-Volterra predator-prey model (Volterra 1928) where the pathogen is the prey and the antibodies are the predator:
Notation:
- \(x(t)\): pathogen concentration at time \(t\)
- \(y(t)\): antibody concentration at time \(t\)
- \(\mu_0\): pathogen growth rate; \(\beta\): pathogen-inactivation strength; \(\mu\): antibody growth rate
Model:
- \[x'(t) = \mu_0\, x(t) - \beta y(t)\]
- \[y'(t) = \mu\, y(t)\]
With baseline antibody concentration \(y(0) = y_{0}\) and initial pathogen concentration \(x(0) = x_{0}\).
Compared to the standard LV model:
the predation term with the \(\beta\) coefficient is missing the prey concentration \(x(t)\) factor; we assume that the efficiency of predation doesn’t depend on pathogen concentration.
the differential equation for predator density is missing the predator death rate term \(-\gamma y(t)\); we assume that as long as there are any pathogens present, the antibody decay rate is negligible compared to the growth rate.
the predator growth rate term \(\mu\, y(t)\) is missing the prey density factor \(x(t)\); we assume that as long as there are any pathogens present, the antibody concentration grows at the same exponential rate.
These omissions were made to simplify the estimation process, under the assumption that they are negligible compared to the other terms in the model.
5.2 Model for post-infection antibody decay
Once the immune response clears the infection, the pathogen concentration reaches zero, so antibody production is no longer stimulated and the antibody concentration decays.
Once the pathogen is cleared (\(x(t) = 0\)), the antibody concentration decays:
- \[y^{\prime}(t) = -\alpha y(t)^r\]
Antibody decay is different from exponential (log–linear) decay. When the shape parameter \(r > 1\), log concentrations decrease rapidly after infection has terminated, but decay then slows down and low antibody concentrations are maintained for a long period. If \(r = 1\), this model reduces to exponential decay with decay rate \(\alpha\).
5.3 Assembling the full response curve
The serum antibody response \(y(t)\) can be written as
\[ y(t) = y_{+}(t) + y_{-}(t) \]
where
\[ \begin{align} y_{+}(t) & = y_{0}\text{e}^{\mu t}[0\leq t <t_{1}]\\ y_{-}(t) & = y_{1}\left(1+(r-1)y_{1}^{r-1}\alpha(t-t_{1})\right)^{-\frac{1}{r-1}}[t_{1}\le t < \infty] \end{align} \]
5.4 Growth rate from peak and baseline
Since the peak level is \(y_{1} = y_{0}\text{e}^{\mu t_{1}}\) the growth rate \(\mu\) can be written as \[\mu = \frac{1}{t_{1}}\log\left(\frac{y_{1}}{y_{0}}\right)\]
5.5 From the full model to the observed curve
The antibody curve \(y(t)\) that we observe and plot falls out of the two-phase system as follows:
During infection (\(t < t_1\)): the antibody equation \(y'(t) = \mu\, y(t)\) gives exponential growth, \(y(t) = y_0\, e^{\mu t}\).
The pathogen equation sets the peak time \(t_1\): circulating antibodies inactivate the pathogen (the \(-\beta y(t)\) term), driving \(x(t)\) to zero at \(t_1\). Once the pathogen is cleared, stimulation stops and decay begins, so the peak antibody level is \(y_1 = y_0\, e^{\mu t_1}\).
After the peak (\(t \ge t_1\)): power-law decay, \(y'(t) = -\alpha y(t)^r\).
So the pathogen trajectory \(x(t)\) never appears in the final antibody curve — it enters only through the peak time \(t_1\). From antibody data alone the pathogen-side quantities (\(\mu_0\), \(\beta\), and the inoculum \(x_0\)) are not separately identifiable (Teunis et al. 2016), so the curve is summarized by five parameters: \(y_0,\ y_1,\ t_1,\ \alpha,\ r\).
Code
library(serocalculator)
library(dplyr)
library(ggplot2)
cur_ai <- "HlyE_IgG"
# `typhoid_curves_nostrat_100` bundles antibody response-curve parameters
# for typhoid, in long format with one MCMC draw per `iter` (`iter` runs
# 1:100; see `?typhoid_curves_nostrat_100`); keep the first 49 draws.
curves_all <-
typhoid_curves_nostrat_100 |>
filter(iter < 50)
curve1 <-
curves_all |>
filter(
iter == 5,
antigen_iso == cur_ai
)
curve1 |>
# Relabel as the (only) curve to plot.
mutate(iter = 1) |>
serocalculator:::plot_curve_params_one_ab(
log_y = FALSE
) +
xlim(0, 100) +
theme_minimal() +
geom_vline(
aes(
xintercept = curve1$t1,
col = "t1"
)
) +
geom_hline(
aes(
yintercept = curve1$y0,
col = "y0"
)
) +
geom_hline(
aes(
yintercept = curve1$y1,
col = "y1"
)
) +
geom_point(
data = curve1,
aes(
x = t1,
y = y1,
col = "(t1,y1)"
)
) +
theme(legend.position = "bottom") +
labs(col = "")
The antibody level at \(t=0\) is \(y_{0}\); the rising branch ends at \(t = t_{1}\) where the peak antibody level \(y_{1}\) is reached. Any antibody level \(y(t) \in (y_{0}, y_{1})\) eventually occurs twice.
An interactive Shiny app lets you manipulate the model parameters:
Interactive Shiny app: https://ucdserg.shinyapps.io/antibody-kinetics-model-2/
5.6 Biological noise
When we measure antibody concentrations in a blood sample, we are essentially counting molecules (using biochemistry).
We might miss some of the antibodies (undercount, false negatives) and we also might incorrectly count some other molecules that aren’t actually the ones we are looking for (overcount, false positives, cross-reactivity).
We are more concerned about overcount (cross-reactivity) than undercount. For a given antibody, we can do some analytical work beforehand to estimate the distribution of overcounts, and add that to our model \(p(Y=y|T=t)\).
Notation:
- \(y_\text{obs}\): measured serum antibody concentration
- \(y_\text{true}\): “true” serum antibody concentration
- \(\epsilon_b\): noise due to probe cross-reactivity
Model:
- \(y_\text{obs} = y_\text{true} + \epsilon_b\)
- \(\epsilon_b \sim \text{Unif}(0, \nu)\)
\(\nu\) needs to be pre-estimated using negative controls, typically using the 95th percentile of the distribution of antibody responses to the antigen-isotype in a population with no exposure.
Why the 95th percentile of negative controls, rather than fitting the full shape of the noise distribution?
(Teunis and Eijkeren 2020) show that a uniform noise model only needs to get the noise width right, not its exact shape: comparing \(\lambda\) estimates under the true (lognormal) noise distribution against estimates under a uniform distribution with the same one-sided 95% range, the two are similar. Getting the width right matters far more than getting the shape right, so a single width parameter – rather than a fully specified distribution – is sufficient.
That still leaves the question of how to estimate the width. The general cross-sectional sample cannot be used directly: it is a mixture of never-infected subjects (pure noise) and subjects at various stages of a real antibody response (signal plus noise), and separating the two requires already knowing who is infected – exactly what the model is trying to estimate. (Teunis and Eijkeren 2020) note this explicitly: “against a background of ongoing seroresponses, even reliable verification of the (95%) width of the distribution may be difficult.”
A negative-control panel (subjects confirmed never infected) avoids this circularity: every observed value in that panel is noise, with no signal mixed in, so its empirical 95th percentile is a direct estimate of the width the model needs.
Why the 95th percentile specifically, rather than the 99th or the 100th (the sample maximum)? (Teunis and Eijkeren 2020) motivates a width-based approach, first framing it as “the one-sided 95%” and later as the “95% range” / “(95%) width” of the noise distribution, but never compares 95% against other levels – so the choice of 95% itself is an adopted convention (it matches how reference intervals are commonly defined in clinical laboratory practice), not a result derived or optimized in the paper. The 100th percentile (the sample maximum) would be a poor choice regardless: it is driven entirely by the single most extreme observation in the panel, so one contaminated or mislabeled negative control inflates \(\nu\) directly, and for any finite panel it underestimates the true theoretical bound. A percentile further into the tail than 95% (e.g., 99%) needs a substantially larger negative-control panel to estimate with comparable precision, for the same reason – order statistics closer to the extreme have higher sampling variance. 95% is the conventional middle ground: far enough into the tail to capture the width that matters, without the fragility of chasing a more extreme quantile.
5.7 Measurement noise
There are also some other sources of noise in our bioassays; user differences in pipetting technique, random ELISA plate effects, etc. This noise can cause both overcount and undercount. We can also estimate the magnitude of this noise source and include it in \(p(Y=y|T=t)\).
Measurement noise, \(\varepsilon\) (“epsilon”), represents measurement error from the laboratory testing process.
Unlike biological noise, measurement noise is multiplicative: the error scales with the true concentration (equivalently, it is additive on a log scale), and it has mean zero, so on average it neither inflates nor deflates the measured concentration.
Notation:
- \(y_\text{obs}\): measured serum antibody concentration
- \(y_\text{true}\): “true” serum antibody concentration
- \(\xi\): relative measurement error
- \(\varepsilon\): bound on the relative error
Model:
- \(y_\text{obs} = y_\text{true} \cdot (1 + \xi)\)
- \(\xi \sim \text{Unif}(-\varepsilon, \varepsilon), \quad 0 \le \varepsilon < 1\)
Here \(\varepsilon\) is the bound on the relative error, not a coefficient of variation (CV). Assay precision is often reported as a CV, the ratio of the standard deviation to the mean for replicates, ideally measured across plates rather than within the same plate. Under this uniform model the CV equals \(\varepsilon/\sqrt{3}\), so a measured CV corresponds to \(\varepsilon = \sqrt{3}\,\text{CV}\).
This follows directly from the variance of a uniform distribution. For \(\xi \sim \text{Unif}(-\varepsilon, \varepsilon)\),
\[ \text{Var}(\xi) = \frac{(2\varepsilon)^2}{12} = \frac{\varepsilon^2}{3} \quad\Longrightarrow\quad \text{SD}(\xi) = \frac{\varepsilon}{\sqrt{3}}. \]
Since \(E[\xi] = 0\) by symmetry, for a fixed true concentration \(y_\text{true}\),
\[ E[y_\text{obs}] = y_\text{true} \cdot E[1 + \xi] = y_\text{true}, \qquad \text{SD}(y_\text{obs}) = y_\text{true} \cdot \text{SD}(\xi) = y_\text{true} \cdot \frac{\varepsilon}{\sqrt{3}}, \]
so
\[ \text{CV} = \frac{\text{SD}(y_\text{obs})}{E[y_\text{obs}]} = \frac{\varepsilon}{\sqrt{3}}. \]
The true concentration \(y_\text{true}\) cancels, so the CV depends only on \(\varepsilon\), not on the concentration being measured – as expected for a relative (multiplicative) error model.
5.8 Combined biological and measurement noise
In practice both noise sources are usually present at once. The two are applied in sequence (Teunis and Eijkeren 2020): biological noise is added to the true concentration, and measurement noise then scales the result.
Model:
- \(y_\text{obs} = (y_\text{true} + \epsilon_b)(1 + \xi)\)
- \(\epsilon_b \sim \text{Unif}(0, \nu)\)
- \(\xi \sim \text{Unif}(-\varepsilon, \varepsilon)\)
This is the model serocalculator uses whenever a noise_params row has both \(\nu > 0\) and \(\varepsilon > 0\); setting either width to zero recovers the corresponding single-source model above.
The observed concentration is a product of two independent random variables, \(y_\text{obs} = S \cdot X\), where \(S = y_\text{true} + \epsilon_b\) and \(X = 1 + \xi\). Conditional on \(y_\text{true}\) (holding it fixed, since only \(\epsilon_b\) and \(\xi\) are random),
\[ E[S] = y_\text{true} + \frac{\nu}{2}, \qquad \text{Var}(S) = \text{Var}(\epsilon_b) = \frac{\nu^2}{12}, \qquad E[X] = 1, \qquad \text{Var}(X) = \text{Var}(\xi) = \frac{\varepsilon^2}{3}, \]
the first pair from \(\epsilon_b \sim \text{Unif}(0,\nu)\), the second from \(\xi \sim \text{Unif}(-\varepsilon,\varepsilon)\) (as derived above). For independent \(S\) and \(X\),
\[ \text{Var}(SX) = \text{Var}(S)\text{Var}(X) + \text{Var}(S)\,E[X]^2 + E[S]^2\,\text{Var}(X). \]
Substituting \(E[X] = 1\) and the variances above,
\[ \text{Var}(y_\text{obs} \mid y_\text{true}) = \frac{\nu^2}{12}\cdot\frac{\varepsilon^2}{3} + \frac{\nu^2}{12}\cdot 1 + \left(y_\text{true} + \frac{\nu}{2}\right)^2 \frac{\varepsilon^2}{3}, \]
which simplifies to
\[ \text{Var}(y_\text{obs} \mid y_\text{true}) = \frac{\nu^2}{12}\left(1 + \frac{\varepsilon^2}{3}\right) + \left(y_\text{true} + \frac{\nu}{2}\right)^2 \frac{\varepsilon^2}{3}. \]
As a check, this reduces to the two single-source cases above: setting \(\varepsilon = 0\) gives \(\nu^2/12\), matching \(\text{Var}(\epsilon_b)\) for biological noise alone; setting \(\nu = 0\) gives \(y_\text{true}^2\,\varepsilon^2/3\), matching the measurement-noise-only CV derivation (\(\text{SD}(y_\text{obs}) = y_\text{true}\cdot\varepsilon/\sqrt{3}\)).
\(\text{Var}(y_\text{obs} \mid y_\text{true})\) conditions on a specific individual’s true concentration. At the population level, even at a fixed time since infection \(T=t\), different individuals have different true antibody concentrations: serodynamics returns a posterior sample of curve-parameter sets rather than a single fixed curve, capturing both estimation uncertainty and between-person heterogeneity in the seroresponse (see “Propagating uncertainty and heterogeneity” in the methodology vignette). Write \(\mu(t) = E[y_\text{true} \mid T=t]\) and \(\sigma^2(t) = \text{Var}(y_\text{true} \mid T=t)\) for this heterogeneity. Unlike \(\epsilon_b\) and \(\xi\), \(y_\text{true}\) is not a named distribution in this framework – its population distribution is the empirical Monte Carlo sample, not a stated closed form – so \(\sigma^2(t)\) has no closed form here and is carried through the derivation symbolically.
Assuming \(y_\text{true}\), \(\epsilon_b\), and \(\xi\) are mutually independent given \(T=t\), the law of total variance gives
\[ \text{Var}(y_\text{obs} \mid T=t) = E\big[\text{Var}(y_\text{obs} \mid y_\text{true}) \mid T=t\big] + \text{Var}\big(E[y_\text{obs} \mid y_\text{true}] \mid T=t\big). \]
For the second term, \(E[y_\text{obs}\mid y_\text{true}] = y_\text{true} + \nu/2\) (a constant shift of \(y_\text{true}\)), so
\[ \text{Var}\big(E[y_\text{obs} \mid y_\text{true}] \mid T=t\big) = \text{Var}(y_\text{true} \mid T=t) = \sigma^2(t). \]
For the first term, take the expectation of the \(\text{Var}(y_\text{obs}\mid y_\text{true})\) formula derived above over \(y_\text{true}\):
\[ E\big[\text{Var}(y_\text{obs}\mid y_\text{true}) \mid T=t\big] = \frac{\nu^2}{12}\left(1+\frac{\varepsilon^2}{3}\right) + \frac{\varepsilon^2}{3}\, E\!\left[\left(y_\text{true}+\frac{\nu}{2}\right)^2 \middle| T=t\right]. \]
Writing \(W = y_\text{true} + \nu/2\) (a constant shift), \(E[W^2] = \text{Var}(W) + E[W]^2 = \sigma^2(t) + (\mu(t)+\nu/2)^2\), so
\[ E\!\left[\left(y_\text{true}+\frac{\nu}{2}\right)^2 \middle| T=t\right] = \sigma^2(t) + \left(\mu(t) + \frac{\nu}{2}\right)^2. \]
Adding the two law-of-total-variance terms together,
\[ \text{Var}(y_\text{obs} \mid T=t) = \left(1+\frac{\varepsilon^2}{3}\right)\left[\frac{\nu^2}{12}+\sigma^2(t)\right] + \frac{\varepsilon^2}{3}\left(\mu(t)+\frac{\nu}{2}\right)^2. \]
This has the same functional form as \(\text{Var}(y_\text{obs}\mid y_\text{true})\) above, with \(y_\text{true}\) replaced by its conditional mean \(\mu(t)\) and an extra between-person heterogeneity term \(\sigma^2(t)\) added alongside the biological-noise variance \(\nu^2/12\). As a check, setting \(\sigma^2(t) = 0\) (no heterogeneity, a single fixed curve so \(y_\text{true} = \mu(t)\) with certainty) recovers \(\text{Var}(y_\text{obs}\mid y_\text{true}=\mu(t))\) exactly, as expected.
5.9 Noise and never-infected subjects
A subject who has never been infected has true antibody concentration \(y_\text{true} = 0\); the probability of never having been infected by age \(a\) is \(\operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}\). The two noise sources treat such a subject very differently (Teunis and Eijkeren 2020):
Biological noise is additive, so a never-infected subject is measured as \(y_\text{obs} = 0 + \epsilon_b \sim \text{Unif}(0, \nu)\). Cross-reactivity can register a positive signal even when there was no true response, so never-infected subjects still contribute a spread-out distribution of small positive values, not a single spike at zero.
Measurement noise is multiplicative, so a never-infected subject is measured as \(y_\text{obs} = 0 \cdot (1 + \xi) = 0\). A relative error cannot move a true zero, so measurement noise alone leaves these subjects exactly at zero.
So the probability mass at \(y = 0\) from never-infected subjects (with weight \(\operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}\)) survives measurement noise but is smeared into small positive values by biological noise. This also explains why the biological-noise width \(\nu\) can be estimated from a known unexposed population: those subjects are essentially all never-infected, so their measured antibody levels are, to good approximation, draws from the biological-noise distribution itself.
Under both noise sources together, the never-infected contribution is neither of these two single-source shapes. Substituting \(y_\text{true} = 0\) into the combined model gives
\[ y_\text{obs} = \epsilon_b (1 + \xi), \qquad \epsilon_b \sim \text{Unif}(0, \nu), \quad \xi \sim \text{Unif}(-\varepsilon, \varepsilon), \]
a product of two independent uniform random variables, not itself uniform. Following the same conditioning argument Teunis and van Eijkeren use to derive their combined-noise density (Teunis and Eijkeren 2020), its density (the never-infected contribution to the overall observed density, including the never-infected weight \(\operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}\)) is
\[ \rho_{bm}(y \mid \text{never-infected}) = \begin{cases} \dfrac{\operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}}{2\varepsilon\nu} \log\!\left(\dfrac{1+\varepsilon}{1-\varepsilon}\right) & 0 < y \le \nu(1-\varepsilon) \\[8pt] \dfrac{\operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}}{2\varepsilon\nu} \log\!\left(\dfrac{\nu(1+\varepsilon)}{y}\right) & \nu(1-\varepsilon) < y < \nu(1+\varepsilon) \\[8pt] 0 & y \ge \nu(1+\varepsilon) \end{cases} \]
The first piece matches Teunis and van Eijkeren’s Equation 19: it is flat (constant in \(y\)) over most of its support, then tapers to zero by \(y = \nu(1+\varepsilon)\). Integrating both pieces over their full range recovers exactly \(\operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}\), the never-infected probability, confirming the density is correctly normalized. As \(\varepsilon \to 0\) the flat piece reduces to \(\operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}/\nu\), the biological-noise-only density on \([0, \nu]\), matching the single-source case above.
5.10 Variation across the posterior
The figure above draws one curve to name the parameters. The fitted model is really a posterior sample of such curves, and autoplot() on a curve_params object shows their spread:
curves holds 49 draws per antigen-isotype here, since it was thinned to iter < 50 when it was loaded.
n_curves caps how many draws are drawn rather than setting a target, so asking for more than the object holds simply draws them all. Passing quantiles = NULL omits the overlay, and log_x = FALSE puts time on a linear scale.
This spread is what the incidence likelihood averages over, as described under “Propagating uncertainty and heterogeneity” above: each draw is one plausible antibody-response curve, and together they carry both the estimation uncertainty and the between-person heterogeneity.
5.11 Multiple biomarkers
A serosurvey usually measures more than one antigen-isotype per participant. The current model treats biomarkers as conditionally independent given the time since infection, so their response densities multiply inside each person’s integral over \(t\):
\[ \operatorname{p}(Y_1=y_1, Y_2=y_2) = \int_t \operatorname{p}(Y_1=y_1|T=t)\,\operatorname{p}(Y_2=y_2|T=t)\,\operatorname{p}_\lambda(T=t)\, dt \]
Note that this is not the product of the two single-biomarker likelihoods: the shared latent infection time \(t\) couples them.
log_likelihood() evaluates the log-likelihood at a given \(\lambda\), which lets us check that claim numerically:
Code
ll_iga <- log_likelihood(
pop_data = xs_data, curve_params = curves, noise_params = noise,
antigen_isos = "HlyE_IgA", lambda = 0.1
)
ll_igg <- log_likelihood(
pop_data = xs_data, curve_params = curves, noise_params = noise,
antigen_isos = "HlyE_IgG", lambda = 0.1
)
ll_both <- log_likelihood(
pop_data = xs_data, curve_params = curves, noise_params = noise,
antigen_isos = c("HlyE_IgG", "HlyE_IgA"), lambda = 0.1
)
c(IgA = ll_iga, IgG = ll_igg, sum = ll_iga + ll_igg, joint = ll_both) IgA IgG sum joint
-3801.688 -4918.504 -8720.192 -8720.192 The two-biomarker value is exactly the sum of the separate ones. That is worth pausing on, because the formula above says it should not be: one integral over a shared \(t\) does not factor into two separate integrals.
What it tells us is that log_likelihood() currently combines biomarkers by summing their individual log-likelihoods, which is the stronger assumption that the biomarkers are independent unconditionally, each with its own latent infection time. est_seroincidence() optimizes through the same code path, so this is the estimation model rather than a quirk of this one function.
Whether the formula or the implementation should change is an open question. Until it is settled, read a multi-biomarker standard error as assuming independence across biomarkers.
R/log_likelihood.R loops over antigen_isos, calls f_dev() once per biomarker, and adds the results, with a source comment naming the independence assumption. Summing log-likelihoods is the log of the product of the marginals.
(Section 9.2 describes work in progress to relax conditional independence given \(t\); the question above is whether the current code implements even that much.)
Using both biomarkers concentrates the log-likelihood around its maximum, and so shrinks the standard error, relative to either biomarker alone (Figure 6):
Code
lik_HlyE_IgA <- graph_loglik( # nolint: object_name_linter.
pop_data = xs_data,
curve_params = curves,
noise_params = noise,
antigen_isos = "HlyE_IgA",
log_x = TRUE
)
lik_HlyE_IgG <- graph_loglik( # nolint: object_name_linter.
previous_plot = lik_HlyE_IgA,
pop_data = xs_data,
curve_params = curves,
noise_params = noise,
antigen_isos = "HlyE_IgG",
log_x = TRUE
)
lik_both <- graph_loglik(
previous_plot = lik_HlyE_IgG,
pop_data = xs_data,
curve_params = curves,
noise_params = noise,
antigen_isos = c("HlyE_IgG", "HlyE_IgA"),
log_x = TRUE
)
print(lik_both)
5.12 Variation in antibody kinetics: by country
Antibody-decay kinetics are not identical across populations. Estimated seroresponse parameters vary by country, age, and serotype (Typhi vs Paratyphi A).
5.13 Variation in antibody kinetics: by age and serotype
This heterogeneity is why we stratify the analysis (for example by country and age group) and why extending the model to handle covariates directly is on the roadmap.
6 Estimating the curves with serodynamics
6.1 Where do the curve parameters come from?
The decay curves above are themselves estimated from longitudinal antibody measurements on confirmed cases — individuals with a known infection date who are sampled repeatedly afterward. Fitting the two-phase model to those data yields, for each antigen-isotype, the parameters the incidence model needs:
- baseline antibody level (\(y_0\))
- peak concentration (\(y_1\))
- time to peak (\(t_1\))
- decay rate (\(\alpha\))
- decay shape (\(r\))
6.2 The serodynamics package
serodynamics is an open-source R package that fits this two-phase within-host kinetics model with a Bayesian hierarchical model, sampled by MCMC (via JAGS). The hierarchical structure stabilizes individual-level estimates by borrowing strength across participants — valuable when longitudinal data are sparse.
The same modeling framework has been applied to pertussis, typhoid, scrub typhus, and Shigella. Previously each application re-implemented the JAGS model specification, data formatting, and post-processing by hand; serodynamics packages that into a single reusable, validated workflow.
6.3 A typical serodynamics workflow
Code
library(serodynamics)
# Longitudinal confirmed-case data (example data ships with the package)
data("nepal_sees")
# Fit the two-phase model by MCMC (JAGS); slow, so not run here
fit <- run_mod(data = nepal_sees, with_post = TRUE)
# Check convergence, then export serocalculator-ready curve parameters
plot_jags_Rhat(fit)
curve_params <- postprocess_jags_output(fit)run_mod() runs several MCMC chains for tens of thousands of iterations, so a real fit takes minutes to hours; the package ships ready-made example output (nepal_sees_jags_output) for experimentation.
6.4 What a fitted model looks like
The package ships a cached example fit; plot_jags_dens() shows the posterior densities of the five kinetic parameters (here HlyE IgG, Typhi), overlaid by MCMC chain:
The overlapping chains indicate good convergence.
Code
data("nepal_sees_jags_output")
plot_jags_dens(nepal_sees_jags_output, iso = "HlyE_IgG", strat = "typhi")nepal_sees_jags_output fit (HlyE IgG, Typhi), by MCMC chain.
6.5 Two packages, one pipeline
-
serodynamics(upstream): longitudinal confirmed-case data \(\rightarrow\) estimated antibody-decay curve parameters -
serocalculator(downstream): a cross-sectional serosurvey plus those curve parameters \(\rightarrow\) a seroincidence estimate
The output of serodynamics feeds directly into serocalculator’s est_seroincidence_by().
6.6 Propagating uncertainty and heterogeneity
The seroresponse model \(\operatorname{p}(Y=y|T=t)\) is not a single fixed curve. serodynamics returns a posterior sample of curve-parameter sets \(\theta^{(1)},\dots,\theta^{(M)}\) (each \(\theta^{(k)} = (y_0, y_1, t_1, \alpha, r)^{(k)}\)), which together capture both
- estimation uncertainty in the kinetics parameters, and
- between-person / between-case heterogeneity in the seroresponse.
Each draw is one plausible antibody-response curve; together they describe the distribution of responses across cases and our uncertainty about it.
The incidence likelihood averages over these draws (Monte Carlo integration), so each person’s contribution becomes
\[ \operatorname{p}(Y=y) \approx \frac{1}{M}\sum_{k=1}^{M} \int_t \operatorname{p}(Y=y|T=t, \theta^{(k)})\; \operatorname{p}_\lambda(T=t)\, dt \]
Because the curve-parameter distribution is marginalized into the likelihood this way, the resulting \(\hat\lambda\) and its standard error reflect the seroresponse heterogeneity and the parameter uncertainty — they are not conditional on a single point-estimated curve. In serocalculator this is the average over the Monte Carlo parameter sets (the iter draws) carried in the curve-parameter object.
6.7 Random effects and variance over time
The longitudinal model distinguishes two kinds of variation:
- Conditional observation noise: given a person’s kinetic parameters \(\theta_i\), residual noise has constant variance on the log-antibody scale.
- Marginal population variation: because \(\theta_i\) includes random effects — notably the waning rate \(\alpha_i\) — individual trajectories can spread apart with increasing time since infection.
Thus, even without an explicit time-varying residual variance \(\sigma^2(t)\), the variance around the population mean trajectory can vary with \(t\):
\[ \operatorname{Var}(\log Y(t)\mid T=t) = \sigma_{\log Y}^2 + \operatorname{Var}_{\theta}\!\left[\log y(t;\theta)\right]. \]
This is the distinction between conditional and marginal heteroscedasticity. The serodynamics observation model uses one residual precision per antigen-isotype, shared across subjects and observation times. It therefore does not fit an additional residual-variance function of time.
Nevertheless, its subject-level kinetic parameters are random effects with a joint distribution. For illustration, if post-peak decay were exponential (\(r=1\)), with the other parameters held fixed,
\[ \log y_i(t) = \log y_1 - \alpha_i(t-t_1), \]
so variation in \(\alpha_i\) alone contributes \((t-t_1)^2\operatorname{Var}(\alpha_i)\) to the population variance on the log scale. The full nonlinear model also allows the other kinetic parameters and their correlations to affect how the variance changes over time.
serocalculator accounts for this trajectory-induced variation by averaging the incidence likelihood over the Monte Carlo sample of kinetic-parameter sets. It does not reclassify that variation as observation noise or replace it with a fitted \(\sigma^2(t)\). Its separate biological- and measurement-noise models have constant additive and relative widths, respectively; neither width is an explicit function of time since infection.
7 Using serocalculator
7.1 An open-source R package
The methods in this lecture are implemented in the open-source serocalculator R package.
The rest of this section walks through the whole estimation pipeline on simulated data, where the true incidence rate is known and the estimate can be checked against it. Section 8 then repeats the exercise across many simulated populations at once.
7.2 Simulating a cross-sectional survey
sim_pop_data() simulates infections as a Poisson process with rate lambda, then generates each person’s antibody levels from the seroresponse model:
Code
sim_lambda <- 0.2 # true incidence rate, per person-year
sim_lifespan <- c(0, 10) # age range of the simulated survey
sim_n <- 100 # number of participants
# Antibody-decay parameters: `typhoid_curves_nostrat_100`, the same bundled
# sample used for the antibody-response curves earlier in this article.
sim_curve_params <-
typhoid_curves_nostrat_100 |>
dplyr::filter(iter < 50, antigen_iso %in% antibodies)
# Biological-noise limits used when *generating* the data
sim_noise_limits <- rbind(
"HlyE_IgA" = c(min = 0, max = 0.5),
"HlyE_IgG" = c(min = 0, max = 0.5)
)
set.seed(54321)
sim_data <- sim_pop_data(
curve_params = sim_curve_params,
lambda = sim_lambda,
n_samples = sim_n,
age_range = sim_lifespan,
antigen_isos = antibodies,
n_mcmc_samples = 0,
renew_params = FALSE,
add_noise = TRUE,
noise_limits = sim_noise_limits,
format = "long"
)Each individual is simulated separately, but their antibodies are modelled jointly, so the biomarkers share one latent infection history.
Three arguments are worth knowing:
-
renew_paramsdecides whether each infection draws a fresh parameter set from the posterior.FALSE(used here) picks one set at birth and keeps it, updating only the baseline levely0before each new infection.TRUEis more realistic but is not what the estimation method assumes — see the bias check at the end of Section 8. -
age_fixedfixes the age at which longitudinal parameters are evaluated; setting it toNAuses the age at infection instead. -
n_mcmc_samples = 0draws a random posterior sample for each infection. Fixing it to a value in1:4000reuses that one sample throughout, which is useful for diagnostics.
7.3 Noise parameters
Estimation needs its own noise parameters, matching Section 5.6. Here they are declared directly, consistent with the noise used to generate the data above:
Code
sim_noise_params <- tibble::tibble(
antigen_iso = antibodies,
nu = 0.5, # biological noise
eps = 0, # measurement noise
y.low = 1, # lower limit of detection
y.high = 5e6 # upper limit of detection
)7.4 The simulated antibody distribution
Code
library(ggbeeswarm)
sim_data |>
ggplot() +
aes(x = as.factor(antigen_iso), y = value) +
geom_beeswarm(
size = .5, alpha = .5,
aes(color = antigen_iso), show.legend = FALSE
) +
geom_boxplot(outlier.colour = NA, fill = NA) +
scale_y_log10() +
theme_linedraw() +
labs(x = "antigen - isotype")
7.5 Estimating incidence from the simulated survey
est_seroincidence() maximizes the log-likelihood \(\log\mathcal{L}(\lambda)\) derived above:
Code
sim_est <- est_seroincidence(
pop_data = sim_data,
sr_params = sim_curve_params,
noise_params = sim_noise_params,
lambda_start = .1,
build_graph = TRUE,
print_graph = FALSE,
antigen_isos = antibodies
)
summary(sim_est)# A tibble: 1 × 11
est.start incidence.rate SE CI.lwr CI.upr se_type coverage log.lik
<dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
1 0.1 0.194 0.0213 0.157 0.241 standard 0.95 -469.
# ℹ 3 more variables: iterations <int>, antigen.isos <chr>,
# nlm.convergence.code <ord>The true rate used to generate these data was 0.2 infections per person-year.
autoplot() draws the log-likelihood curve, with the maximum marked:
Code
autoplot(sim_est, log_x = TRUE)
7.6 Estimating seroincidence in real data
The same call works on a real serosurvey. est_seroincidence_by() estimates separately within each stratum, and summary() reports each stratum’s \(\hat\lambda\) with a standard error and confidence interval (cluster-robust, if cluster_var was supplied):
Code
# Load antibody-decay curve parameters and cross-sectional population data
curves <- "https://osf.io/download/rtw5k/" |> load_sr_params()
xs_data <- "https://osf.io/download/n6cp3/" |> load_pop_data()
noise <- url("https://osf.io/download/hqy4v/") |> readRDS()
# Visualize the cross-sectional antibody distribution
xs_data |> autoplot(strata = "Country", type = "density")
# Estimate incidence, stratified by country and age group
est <- est_seroincidence_by(
pop_data = xs_data,
sr_params = curves,
noise_params = noise,
strata = c("Country", "ageCat"),
antigen_isos = c("HlyE_IgG", "HlyE_IgA")
)
summary(est)Stratifying this way also makes the constant-and-homogeneous-incidence assumption more plausible within each stratum, and lets each stratum use the seroresponse parameters appropriate to it.
7.7 Interactive Shiny app
A point-and-click interface is available at https://ucdserg.shinyapps.io/shiny_serocalculator/.
8 Validation: recovering known incidence rates
8.1 Simulating clusters with known incidence
The single survey above recovered a rate close to the one that generated it. To check that this was not luck, we simulate many clusters across a range of known incidence rates and ask whether the estimates track them.
Code
sim_lambdas <- c(.05, .1, .15, .2, .5, .8) # true incidence rates
sim_nclus <- 20 # clusters per rate
sim_df <- sim_pop_data_multi_cached(
curve_params = sim_curve_params,
lambdas = sim_lambdas,
nclus = sim_nclus,
sample_sizes = sim_n,
age_range = sim_lifespan,
antigen_isos = antibodies,
renew_params = FALSE,
add_noise = TRUE,
noise_limits = sim_noise_limits,
format = "long",
cache_path = sim_cache,
cache_id = "sim_df"
)That is 6 rates times 20 clusters times 100 participants — 1.2^{4} simulated individuals.
Simulating and estimating on that many people is by far the most expensive work in this document, and it renders to three formats, so the work would otherwise be repeated three times per build. sim_pop_data_multi_cached() and est_seroincidence_by_cached() write their results to a committed cache directory keyed on their arguments, so a build reuses them and only recomputes when the simulation settings actually change.
Code
sim_df |>
ggplot() +
aes(x = as.factor(cluster), y = value) +
geom_beeswarm(size = .2, alpha = .3, aes(color = antigen_iso)) +
geom_boxplot(outlier.colour = NA, fill = NA) +
scale_y_log10() +
facet_wrap(~ antigen_iso + lambda.sim, nrow = 2) +
theme_linedraw() +
theme(legend.position = "bottom") +
labs(x = "cluster")
8.2 Estimating incidence in each cluster
Code
ests <- est_seroincidence_by_cached(
pop_data = sim_df,
sr_params = sim_curve_params,
noise_params = sim_noise_params,
strata = c("sample_size", "lambda.sim", "cluster"),
curve_strata_varnames = NULL,
noise_strata_varnames = NULL,
antigen_isos = antibodies,
# build_graph = TRUE attaches a log-likelihood plot to every stratum's
# fit, which grows large across many strata; re-fit individual curves on
# demand instead (see below).
build_graph = FALSE,
cache_path = sim_cache,
cache_id = "ests"
)
ests_summary <- ests |> summary()Code
if (knitr::is_html_output()) {
ests_summary |>
DT::datatable(options = list(scrollX = TRUE)) |>
DT::formatRound(
columns = c("incidence.rate", "SE", "CI.lwr", "CI.upr", "log.lik")
)
} else {
ests_summary
}The interactive table is HTML-only; other formats print the plain summary, since converting the widget to a static image would need a headless browser.
8.3 Inspecting individual clusters
The cached fit above carries no log-likelihood curves, so the few clusters we want to look at are re-fit with build_graph = TRUE. Strata are selected by the columns that define them, because stratum labels are assigned positionally and so are not stable across a subset:
Code
stratum_ids <- ests_summary |>
tibble::as_tibble() |>
dplyr::select(sample_size, lambda.sim, cluster)
ests_first_five <- refit_strata(
pop_data = sim_df,
strata_ids = stratum_ids |> dplyr::slice(1:5),
sr_params = sim_curve_params,
noise_params = sim_noise_params,
curve_strata_varnames = NULL,
noise_strata_varnames = NULL,
antigen_isos = antibodies,
num_cores = 1
)Code
autoplot(ests_first_five, log_x = TRUE)
8.4 Checking nlm() convergence
est_seroincidence() maximizes the log-likelihood numerically with stats::nlm(). Exit codes 3 to 5 flag a possible non-convergence, so they are worth checking before trusting any of the estimates:
Code
[1] 13Remedies, all passed through to nlm():
- code 3: decrease
stepmin - code 4: increase
iterlim - code 5: increase
stepmax
Pass them with num_cores = 1. est_seroincidence_by() currently drops these arguments in its parallel branch (#629), so under parallelism they are silently ignored.
Code
if (nrow(problem_ids) > 0) {
refit_strata(
pop_data = sim_df,
strata_ids = problem_ids,
sr_params = sim_curve_params,
noise_params = sim_noise_params,
curve_strata_varnames = NULL,
noise_strata_varnames = NULL,
antigen_isos = antibodies,
num_cores = 1
) |>
autoplot(log_x = TRUE)
}Warning: Some strata are completely missing one or more biomarkers.
, , cluster = 1, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 100 0 100 0
, , cluster = 2, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 0 100
, , cluster = 5, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 0 100
, , cluster = 7, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 100 0 0 0
, , cluster = 8, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 100 0 0
, , cluster = 9, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 100 0 0
, , cluster = 12, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 100 100 0
, , cluster = 14, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 100 0
, , cluster = 15, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 100 0 0 0
, , cluster = 16, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 100 0
, , cluster = 20, antigen_iso = HlyE_IgA
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 100 0
, , cluster = 1, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 100 0 100 0
, , cluster = 2, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 0 100
, , cluster = 5, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 0 100
, , cluster = 7, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 100 0 0 0
, , cluster = 8, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 100 0 0
, , cluster = 9, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 100 0 0
, , cluster = 12, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 100 100 0
, , cluster = 14, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 100 0
, , cluster = 15, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 100 0 0 0
, , cluster = 16, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 100 0
, , cluster = 20, antigen_iso = HlyE_IgG
lambda.sim
sample_size 0.1 0.15 0.2 0.8
100 0 0 100 0Warning: `nlm()` may not have reached the maximum likelihood estimate.
`nlm()` completed with the following convergence code:
3: Last global step failed to locate a point lower than x. Either x is an
approximate local minimum of the function, the function is too non-linear for
this algorithm, or `stepmin` in `est_seroincidence()` (a.k.a. `steptol` in
`nlm()`) is too large.
`nlm()` may not have reached the maximum likelihood estimate.
`nlm()` completed with the following convergence code:
3: Last global step failed to locate a point lower than x. Either x is an
approximate local minimum of the function, the function is too non-linear for
this algorithm, or `stepmin` in `est_seroincidence()` (a.k.a. `steptol` in
`nlm()`) is too large.
8.5 Estimates recover the simulated rates
Code
ests_summary |>
autoplot(
type = "scatter",
xvar = "lambda.sim",
CI = TRUE,
dodge_width = .05
) +
ggplot2::geom_function(
fun = function(x) x,
col = "red",
aes(linetype = "data-generating incidence rate")
) +
labs(linetype = "") +
scale_x_log10()
Each point is one simulated cluster. The estimates scatter around the identity line, and the confidence intervals widen at higher incidence.
analyze_sims() summarizes the simulation results across clusters:
Code
ests_summary |> analyze_sims()# A tibble: 6 × 8
lambda.sim sample_size Bias Mean_Est_SE Empirical_SE RMSE Mean_CI_Width
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.05 100 0.00292 0.00905 0.0102 0.0103 0.0362
2 0.1 100 0.00894 0.0144 0.0204 0.0218 0.0571
3 0.15 100 0.0222 0.0199 0.0183 0.0285 0.0786
4 0.2 100 0.0269 0.0242 0.0319 0.0411 0.0957
5 0.5 100 0.0857 0.0568 0.0758 0.113 0.224
6 0.8 100 0.172 0.0989 0.115 0.205 0.390
# ℹ 1 more variable: CI_Coverage <dbl>Code
ests_summary |>
analyze_sims() |>
autoplot(statistic = "Empirical_SE")
8.6 A caveat: renewing parameters at each infection
Everything above sets renew_params = FALSE, matching what the estimation method assumes. Setting it to TRUE — drawing fresh seroresponse parameters at every infection — is more realistic, and is not accounted for by the current method, so estimates can be biased in high-incidence populations:
Code
sim_df_renew <- sim_pop_data_multi_cached(
curve_params = sim_curve_params,
lambdas = sim_lambdas,
nclus = sim_nclus,
sample_sizes = sim_n,
age_range = sim_lifespan,
antigen_isos = antibodies,
renew_params = TRUE,
add_noise = TRUE,
noise_limits = sim_noise_limits,
format = "long",
cache_path = sim_cache,
cache_id = "sim_df_renew"
)
ests_renew_summary <-
est_seroincidence_by_cached(
pop_data = sim_df_renew,
sr_params = sim_curve_params,
noise_params = sim_noise_params,
strata = c("sample_size", "lambda.sim", "cluster"),
curve_strata_varnames = NULL,
noise_strata_varnames = NULL,
antigen_isos = antibodies,
build_graph = FALSE,
cache_path = sim_cache,
cache_id = "ests_renew"
) |>
summary()Code
ests_renew_summary |>
autoplot(
type = "scatter",
xvar = "lambda.sim",
CI = TRUE,
dodge_width = .05
) +
ggplot2::geom_function(
fun = function(x) x,
col = "red",
aes(linetype = "data-generating incidence rate")
) +
labs(linetype = "") +
scale_x_log10()
9 Ongoing and future work
9.1 In-progress work
Extending and improving existing Shiny apps for these methods (e.g., the serocalculator Shiny app, source)
Multivariate modeling of biomarkers (relaxing conditional independence)
A graphical (Shiny) app for
serodynamicsModeling time-varying incidence rates
Accounting for re-exposure
Accounting for latent immunocompromised subpopulations
Calibrating to population demographics
9.2 Multiple biomarkers: beyond conditional independence
With several biomarkers (e.g. HlyE IgA and IgG), the current method treats them as conditionally independent given the time since infection, so the joint likelihood factors into a product of per-biomarker terms:
\[ \operatorname{p}(Y_1 = y_1, Y_2 = y_2 \mid T=t) = \operatorname{p}(Y_1 = y_1 \mid T=t)\,\operatorname{p}(Y_2 = y_2 \mid T=t) \]
This is the same conditional-independence simplification that naive Bayes makes (features independent given the class): convenient, but it ignores any within-person correlation between biomarkers.
. . .
Kwan Ho Lee (UCD-SERG) is relaxing this to allow covariance among biomarkers — a multivariate seroresponse model with a Kronecker-structured covariance, fit in Stan (UCD-SERG/shigella#13).
References
Aiemjoy, K., Seidman J. C., Saha S., Munira S. J., Islam Sajib M. S., and Sarkar Sium S. M. al. 2022. “Estimating Typhoid Incidence from Community-Based Serosurveys: A Multicohort Study.” The Lancet Microbe 3 (8): e578–87. https://doi.org/10.1016/S2666-5247(22)00114-8.
Cameron, A. Colin, Jonah B. Gelbach, and Douglas L. Miller. 2011. “Robust Inference with Multiway Clustering.” Journal of Business & Economic Statistics 29 (2): 238–49. https://doi.org/10.1198/jbes.2010.07136.
Hay, James A., Isobel Routledge, and Saki Takahashi. 2024. “Serodynamics: A Primer and Synthetic Review of Methods for Epidemiological Inference Using Serological Data.” Epidemics 49: 100806. https://doi.org/10.1016/j.epidem.2024.100806.
Huber, Peter J. 1967. “The Behavior of Maximum Likelihood Estimates Under Nonstandard Conditions.” In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Statistics. University of California Press.
Liang, Kung-Yee, and Scott L. Zeger. 1986. “Longitudinal Data Analysis Using Generalized Linear Models.” Biometrika 73 (1): 13–22. https://doi.org/10.1093/biomet/73.1.13.
Noordzij, Marlies, Friedo W. Dekker, Carmine Zoccali, and Kitty J. Jager. 2010. “Measures of Disease Frequency: Prevalence and Incidence.” Nephron Clinical Practice 115 (1): c17–20. https://doi.org/10.1159/000286345.
Teunis, P. F. M., and J. C. H. van Eijkeren. 2020. “Estimation of Seroconversion Rates for Infectious Diseases: Effects of Age and Noise.” Statistics in Medicine 39 (21): 2799–814. https://doi.org/10.1002/sim.8578.
Teunis, P. F. M., J. C. H. van Eijkeren, W. F. de Graaf, A. Bonačić Marinović, and M. E. E. Kretzschmar. 2016. “Linking the Seroresponse to Infection to Within-Host Heterogeneity in Antibody Production.” Epidemics 16 (September): 33–39. https://doi.org/10.1016/j.epidem.2016.04.001.
Volterra, Vito. 1928. “Variations and Fluctuations of the Number of Individuals in Animal Species Living Together.” ICES Journal of Marine Science 3 (1): 3–51.
White, Halbert. 1980. “A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity.” Econometrica 48 (4): 817–38.