Estimating Incidence Rates from Cross-Sectional Serosurveys

UC Davis Seroepidemiology Research Group (UCD-SERG)

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)

Salmonella Typhi, flagellar stain. Photo: CDC/PHIL (public domain; PHIL ID 2115).

Global distribution of typhoid fever. Source: Wikimedia Commons, CC BY-SA 3.0.

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)

Shigella in a stool sample. Photo: CDC/PHIL (public domain), via Wikimedia Commons.

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

Trombicula mite larva with stylostome. Photo: Alan R Walker / Wikimedia Commons, CC BY-SA 3.0.

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.

Countries with routine TCV introduction as of 2023. Source: CDC MMWR 72(7):182–188 (public domain).

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

Definition (Population incidence rate)

Definition 1 (Population incidence rate) The incidence rate of a disease over a specific time period is the rate at which individuals in a population are acquiring the disease during that time period (Noordzij et al. 2010).

Example (Population incidence rate)

Example 1 (Population incidence rate) If there are 10 new cases of typhoid in a population of 1000 persons during a one month time period, then the incidence rate for that time period is 10 new cases per 1000 persons per month.

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

Definition (Constant-incidence model)

Definition 2 (Constant-incidence model) In the constant-incidence model:

• incidence is constant over time and homogeneous across the population, so a single number describes the whole sample, \[\lambda_{i,t} = \lambda, \quad \forall i,t\] (subpopulations can be analyzed separately to make this more plausible); and

• there is no lasting immunity: participants are always at risk of a new infection, regardless of how recently they were last infected.

(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

Proposition (Distribution of time since last infection)

Proposition 1 (Distribution of time since last infection) Under the constant-incidence model (Definition 2), a participant’s time since last infection, \(T\), is exponentially distributed with rate \(\lambda\), truncated at their age \(a\): no one can have been infected before they were born, so the density is supported on \([0,a]\) and the remaining probability, \(\operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{}\), falls on “never infected”:

\[ \operatorname{p}(T=t|A=a) = 1_{t \in[0,a]} \textcolor{red}{\lambda \operatorname{exp}\mathopen{}\left\{-\lambda t\right\}\mathclose{}} + 1_{t = \text{NA}} \operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{} \]

Its rate parameter is exactly the incidence rate of Definition 1, and that identity is what makes \(T\) worth recovering.

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.

Definition (Likelihood)

Definition 3 (Likelihood) The likelihood is the probability (or density) of the observed data, read as a function of the parameter rather than of the data:

\[\mathscr{L}(\lambda) \stackrel{\text{def}}{=}\operatorname{p}(\text{data} \mid \lambda)\]

Definition (Log-likelihood)

Definition 4 (Log-likelihood) The log-likelihood is the log of the likelihood (Definition 3):

\[\ell(\lambda) \stackrel{\text{def}}{=}\operatorname{log}\mathopen{}\left\{\mathscr{L}(\lambda)\right\}\mathclose{}\]

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).

Definition (Score function)

Definition 5 (Score function) The score function is the derivative of the log-likelihood with respect to the parameter (in general, its gradient):

\[\ell'(\lambda) \stackrel{\text{def}}{=}\frac{\partial}{\partial \lambda}\ell(\lambda)\]

A maximum likelihood estimate \(\hat\lambda\) solves the score equation \(\ell'(\hat\lambda) = 0\).

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:

Proposition (MLE when infection times are fully observed)

Proposition 2 (MLE when infection times are fully observed) If the latent infection times \(t_1, \ldots, t_n\) were fully observed, solving the score equation (Definition 5) would estimate \(\lambda\) as the reciprocal of their mean:

\[ \begin{aligned} \hat{\lambda}_{\text{ML}} &= \frac{n}{\sum_{i=1}^n t_i} \\&= \frac{1}{\bar{t}} \end{aligned} \]

4.8 Example log-likelihood curves

Here’s what that would look like:

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")

Figure 1: log-likelihood curve for latent data

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.

Definition (Hessian)

Definition 6 (Hessian) The Hessian is the second derivative of the log-likelihood with respect to the parameter (in general, the matrix of its second partial derivatives):

\[\ell''(\lambda) \stackrel{\text{def}}{=}\frac{\partial^2}{\partial \lambda^2}\ell(\lambda)\]

Definition (Observed information)

Definition 7 (Observed information) The observed information is the negative Hessian (Definition 6):

\[I(\lambda) \stackrel{\text{def}}{=}-\ell''(\lambda)\]

A sharper peak at the maximum means a larger observed information.

Theorem (Asymptotic variance of the MLE)

Theorem 1 (Asymptotic variance of the MLE) Under standard regularity conditions, the variance of a maximum likelihood estimate is approximately the inverse of the observed information (Definition 7) at that estimate:

\[\mathop{\widehat{\operatorname{Var}}}\nolimits\mathopen{}\left(\hat\lambda\right)\mathclose{} \approx \mathopen{}\left[I(\hat\lambda)\right]\mathclose{}^{-1}\]

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

Example (Latent-time exponential model)

Example 2 (Latent-time exponential model) For the latent-time model of Proposition 2, the log-likelihood is \[\ell^*(\lambda) = n\operatorname{log}\mathopen{}\left\{\lambda\right\}\mathclose{} - \lambda\sum_{i=1}^n t_i\]

so its score (Definition 5) and Hessian (Definition 6) are

\[ \begin{aligned} \ell^{*\prime}(\lambda) &= \frac{n}{\lambda} - \sum_{i=1}^n t_i \\ \ell^{*\prime\prime}(\lambda) &= -\frac{n}{\lambda^2} \end{aligned} \]

and its observed information (Definition 7) is \(-\ell^{*\prime\prime}(\lambda) = n/\lambda^2\).

At the MLE \(\hat\lambda = 1/\bar t\), Theorem 1 then gives

\[ \begin{aligned} \mathop{\widehat{\operatorname{Var}}}\nolimits\mathopen{}\left(\hat\lambda\right)\mathclose{} &\approx \left[-\ell^{*\prime\prime}(\hat\lambda)\right]^{-1} \\&= \frac{\hat\lambda^2}{n} \\ \mathop{\widehat{\operatorname{SE}}}\nolimits\mathopen{}\left(\hat\lambda\right)\mathclose{} &\approx \frac{\hat\lambda}{\sqrt{n}} \end{aligned} \]

so the standard error shrinks like \(1/\sqrt{n}\): quadrupling the sample size halves it.

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

sim_case_data |>
  autoplot(alpha = .5)

Figure 2: Antibody response curves, \(p(Y=y|T=t)\), for typhoid

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)

Figure 3: Observed antibody measurements over time since fever onset, for typhoid

4.13 The per-person likelihood

Proposition (Marginal likelihood of one participant’s data)

Proposition 3 (Marginal likelihood of one participant’s data) Substituting \(\operatorname{p}(Y=y,T=t) = \operatorname{p}(Y=y|T=t)\,\operatorname{p}(T=t)\) into the expression for \(\operatorname{p}(Y=y)\) above gives one participant’s contribution:

\[\operatorname{p}(Y=y) = \int_t \operatorname{p}(Y=y|T=t)\,\operatorname{p}(T=t)\, dt\]

That integral is shorthand, because \(T\) is not purely continuous: it has the mixed distribution of Proposition 1, continuous on \([0,a]\) plus an atom at \(T=\text{NA}\). Written out over both parts, the integral is

\[ \operatorname{p}(Y=y) = \int_0^a \operatorname{p}(Y=y\mid T=t)\,\lambda\operatorname{exp}\mathopen{}\left\{-\lambda t\right\}\mathclose{}\,dt \;+\; \operatorname{p}(Y=y\mid T=\text{NA})\cdot \operatorname{exp}\mathopen{}\left\{-\lambda a\right\}\mathclose{} \]

The first term is the contribution from subjects infected at least once; the second is the contribution from subjects never infected.

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

Proposition (Full-sample likelihood, assuming independent participants)

Proposition 4 (Full-sample likelihood, assuming independent participants) Write \(\mathscr{L}_i\) and \(\ell_i\) for participant \(i\)’s own likelihood and log-likelihood, which is Proposition 3 evaluated at \(y = y_i\):

\[ \begin{aligned} \mathscr{L}_i(\lambda) &= \operatorname{p}(Y = y_i) \\ \ell_i(\lambda) &= \operatorname{log}\mathopen{}\left\{\mathscr{L}_i(\lambda)\right\}\mathclose{} \end{aligned} \]

If participants are independent of one another, the likelihood (Definition 3) of the whole sample \(\tilde{y} = (y_1, y_2, \ldots, y_n)\) is their product:

\[\mathscr{L}(\lambda) = \prod_{i=1}^n \mathscr{L}_i(\lambda)\]

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.

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.

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).

Definition (Sandwich variance estimator)

Definition 8 (Sandwich variance estimator) Let \(C\) be the total number of non-overlapping clusters in the sample, and let \(H = \ell''(\hat\lambda)\) be the Hessian (Definition 6) at the MLE — serocalculator obtains it from stats::nlm()’s own hessian output. Write \(U_i\) for observation \(i\)’s score contribution (Definition 5) — the gradient of its own log-likelihood \(\ell_i\) (Proposition 4) — and \(U_c\) for the total over the observations in cluster \(c\):

\[ \begin{aligned} U_i &= \nabla_\lambda \ell_i(\lambda) \\ U_c &= \sum_{i \in c} U_i \end{aligned} \]

where \(\nabla_\lambda\) denotes the gradient with respect to \(\lambda\). These are gradients of the same \(\ell_i\) that Proposition 4 defines: the full-sample log-likelihood sums \(\ell_i\) over every observation, and the sandwich accumulates their gradients cluster by cluster. The cluster-robust variance-covariance matrix for the parameter estimates is then

\[ \begin{aligned} V_{\text{robust}} &= H^{-1} B H^{-1} \\ B &= \sum_{c=1}^C U_c U_c^T \end{aligned} \]

with \(B\) the “meat” of the sandwich.

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.

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.

Definition (Multi-way clustering-robust variance)

Definition 9 (Multi-way clustering-robust variance) Let \(\mathcal{K} = \{k_1, \dots, k_p\}\) denote the cluster variables. For a non-empty subset \(S \subseteq \mathcal{K}\), let \(V_S\) be the one-way sandwich variance \(V_{\text{robust}}\) (Definition 8), computed by clustering on the interaction of the variables in \(S\). The multi-way variance is the signed sum over every such subset:

\[V_{\text{multiway}} = \sum_{\emptyset \neq S \subseteq \mathcal{K}} (-1)^{|S| + 1} V_S\]

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.

Example (Two cluster variables)

Example 3 (Two cluster variables) For two cluster variables (\(p = 2\)), Definition 9 reduces to the familiar three-term formula:

\[V_{\text{multiway}} = V_1 + V_2 - V_{1 \cap 2}\]

where \(V_1\) and \(V_2\) cluster on each variable separately, and \(V_{1 \cap 2}\) clusters on their interaction, the most granular grouping.

Small-sample correction

Two adjustments refine the one-way variance terms \(V_S\) that feed Definition 9.

Definition (Finite-difference score approximation)

Definition 10 (Finite-difference score approximation) serocalculator approximates each cluster’s score contribution \(U_c\) (Definition 8) by a forward finite difference of the cluster’s log-likelihood, \(\ell_c(\lambda) = \sum_{i \in c} \log p(Y_i \mid \lambda)\):

\[U_c \approx \frac{\ell_c(\hat \lambda+ \epsilon) - \ell_c(\hat \lambda)}{\epsilon}, \quad \epsilon = 10^{-6}\]

Definition (CR1 small-sample correction)

Definition 11 (CR1 small-sample correction) When small_sample = "CR1" (the default), each one-way subset variance \(V_S\) (Definition 9) is inflated by the usual cluster degrees-of-freedom correction:

\[V_S^{\text{CR1}} = \frac{G_S}{G_S - 1} \, V_S\]

where \(G_S\) is the number of unique clusters in subset \(S\)’s grouping.

Safeguards against invalid variance

Definition (Flooring and fallback for the final variance)

Definition 12 (Flooring and fallback for the final variance) The inclusion-exclusion sum \(V_{\text{multiway}}\) (Definition 9) is not guaranteed to be non-negative. serocalculator therefore floors the raw estimate at 0 (emitting a warning), independent of floor_to_standard, before it can reach \(\sqrt{V}\) as a silent NaN:

\[V_{\text{raw}}^{+} = \max\mathopen{}\left(0, V_{\text{multiway}}\right)\mathclose{}\]

Setting floor_to_standard = TRUE additionally floors the result at the model-based (non-cluster-robust) variance \(V_{\text{standard}} = H^{-1}\), the inverse of the observed information at the MLE:

\[V_{\text{final}} = \max\mathopen{}\left(V_{\text{standard}}, V_{\text{raw}}^{+}\right)\mathclose{}\]

If any one-way subset term is unavailable, for example because its Hessian is non-finite or non-positive, the inclusion-exclusion sum is itself non-finite. serocalculator warns and returns a missing standard error in that case. With floor_to_standard = TRUE it falls back to \(V_{\text{standard}}\) instead, unless \(V_{\text{standard}}\) is itself unusable (the same degenerate Hessian poisons \(V_{\text{standard}} = H^{-1}\) too), in which case it also returns a missing standard error rather than silently falling back to a negative or infinite “variance”.

The inclusion-exclusion sum’s occasional negativity is a documented property of the Cameron-Gelbach-Miller estimator, not a bug.

Implementation in serocalculator

Users can specify clustering using the cluster_var parameter:

# 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():

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{}\)).

Effect on results

• Point estimates (incidence rates) are unchanged: cluster_var affects 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.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\).

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 = "")

Figure 4: An example kinetics curve for HlyE IgG

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/

QR code for the antibody-kinetics Shiny app

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 |> autoplot(quantiles = c(0.1, 0.5, 0.9), log_x = TRUE)

Figure 5: Posterior draws of the antibody-decay curve per antigen-isotype, with the 10%, 50%, and 90% point-wise quantiles overlaid.

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:

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):

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)

Figure 6: Log-likelihood curves for HlyE IgA alone, HlyE IgG alone, and both biomarkers combined.

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).

Decay curves (A), peak antibody levels (B), and decay rates (C) for HlyE and LPS (IgG/IgA), by country (Bangladesh, Pakistan, Nepal, Ghana). SEES study.

5.13 Variation in antibody kinetics: by age and serotype

Decay curves by age stratum (<5, 5–15, 16+) and serotype (Typhi vs Paratyphi A), for HlyE and LPS (IgG/IgA). SEES study.

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

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.

data("nepal_sees_jags_output")
plot_jags_dens(nepal_sees_jags_output, iso = "HlyE_IgG", strat = "typhi")

Posterior densities of the five kinetic parameters from the cached 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:

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_params decides 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 level y0 before each new infection. TRUE is more realistic but is not what the estimation method assumes — see the bias check at the end of Section 8.
• age_fixed fixes the age at which longitudinal parameters are evaluated; setting it to NA uses the age at infection instead.
• n_mcmc_samples = 0 draws a random posterior sample for each infection. Fixing it to a value in 1:4000 reuses 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:

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

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")

Figure 7: Distribution of simulated antibody responses, by antigen-isotype.

7.5 Estimating incidence from the simulated survey

est_seroincidence() maximizes the log-likelihood \(\log\mathcal{L}(\lambda)\) derived above:

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:

autoplot(sim_est, log_x = TRUE)

Figure 8: Log-likelihood curve for the simulated survey, on a log x-axis.

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):

# 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/.

The serocalculator Shiny app.

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.

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.

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")

Figure 9: Simulated antibody responses by cluster, faceted by antigen-isotype and by the true incidence rate used to generate each cluster.

8.2 Estimating incidence in each cluster

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()

summary() returns a tibble() carrying some extra metadata about the stratification.

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
}

Table 1

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:

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
)

autoplot(ests_first_five, log_x = TRUE)

Figure 10: Log-likelihood curves for the first five simulated clusters.

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:

problem_ids <- ests_summary |>
  tibble::as_tibble() |>
  dplyr::filter(nlm.convergence.code > 2) |>
  dplyr::select(sample_size, lambda.sim, cluster)

nrow(problem_ids)

[1] 13

Remedies, 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.

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   0

Warning: `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.

Figure 11: Log-likelihood curves for any clusters whose optimizer reported a possible convergence problem.

8.5 Estimates recover the simulated rates

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()

Figure 12: Estimated incidence rates (points, with 95% confidence intervals) for each simulated cluster, against the true rate used to generate it. The red line is the identity.

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:

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>

ests_summary |>
  analyze_sims() |>
  autoplot(statistic = "Empirical_SE")

Figure 13: Empirical standard error of the incidence estimates, by true rate.

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:

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()

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()

Figure 14: Recovery of the true incidence rate when seroresponse parameters are renewed at each infection. Compare with Figure 12.

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 serodynamics

• Modeling 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.