The epidist package enables users to estimate delay distributions while accounting for common reporting biases.
This vignette first introduces the background (Park et al. 2024) required to understand the models implemented in epidist.
It then goes on to explain each model in turn.
1 Background
Estimating a delay distribution may appear to be simple: one could imagine fitting a probability distribution to a set of observed delays. However, observed delays are often biased during an ongoing outbreak. We begin by presenting a formalism for characterizing delay distributions as well as two main biases (truncation and censoring) affecting them. We then present statistical models for correcting for these biases.
Any epidemiological delay requires a primary (starting) and a secondary (ending) event. For example, the incubation period measures the time between infection (primary event) and symptom onset (secondary event). Here, we use \(p\) to denote the time of the primary event, and \(s\) to denote the time of the secondary event.
1.1 Forwards and backwards distributions
We can measure any delay distribution in two different ways. First, we can measure the forward distribution \(f_p(\tau)\), starting from a cohort of individuals who experienced the primary event at the same time \(p\) and looking at when they experienced their secondary event. Second, we can measure the backward distribution \(b_s(\tau)\), starting from a cohort of individuals who experienced the secondary event at the same time \(s\) and looking at when they experienced their primary event. While the length of each individual delay \(\tau = p-s\) remains constant whether we look at it forward or backwards, the shape of the distribution is affected by the differences in perspectives due to the differences in cohort composition.
To illustrate their differences, let’s assume that primary and secondary events occur at rates, or equivalently incidences, \(\mathcal{P}(p)\) and \(\mathcal{S}(s)\), respectively. Then, the total density \(\mathcal{T}(p, s)\) of individuals with a primary event at time \(p\) and secondary event at time \(s\) can be expressed equivalently in terms of both forward and backward distributions: \[ \mathcal{T}(p, s) = \mathcal{P}(p) f_p(s - p) = \mathcal{S}(s) b_s(s - p). \tag{1.1} \] Rearranging Equation (1.1) gives \[ b_s(\tau) = \frac{\mathcal{P}(s - \tau) f_{s - \tau}(\tau)}{\mathcal{S}(s)}. \tag{1.2} \] The denominator of Equation (1.2), which corresponds to the incidence of secondary events, may be expressed as the integral over all possible delays \[ \mathcal{S}(s) = \int_{-\infty}^\infty \mathcal{P}(s - \tau) f_{s - \tau}(\tau) \text{d} \tau, \] such that \[ b_s(\tau) = \frac{\mathcal{P}(s - \tau) f_{s - \tau}(\tau)}{\int_{-\infty}^\infty \mathcal{P}(s - \tau) f_{s - \tau}(\tau) \text{d} \tau}. \] Here, we see that \(b_s(\tau)\) depends not only on \(f_{s - \tau}(\tau)\) but also on \(\mathcal{P}(s - \tau)\), meaning that past changes in the incidence pattern will affect the shape of the distribution. Particularly, when an epidemic is growing, we are more likely to observe shorter delays, causing an underestimation of the mean delay. Therefore, we always want to characterize epidemiological delays from the forward perspective and estimate the forward distribution. For this reason, our current methodology focuses on biases that affect the estimation of the forward distribution.
1.2 Right truncation
One key bias that affects the forward distribution is right truncation. Right truncation refers to the bias arising from the inability to observe future events and occurs when we observe data based on the secondary events. For example, assume the data are right truncated and we don’t observe secondary events past time \(T\). Then, we will only observe delays whose secondary events occurred before time \(T\), causing us to underestimate the mean of the distribution as these delays will on average be shorter.
Bias from right truncation is greater when events are more likely to be more recent. A common example of severely right truncated data is data collected during outbreaks when growth in incidence is exponential (so you are much more likely to have a recent event). On the other hand, if data collection is continued until the end of an outbreak then many fewer events are likely to be more recent and so there will be little right truncation in general.
Mathematically right truncation can be described as follows. Let \(P\) and \(S\) be random variables. Let \(F_p\) be the forward cumulative distribution. Then, the probability of observing a delay of length \(\tau\) given that the primary event occurred at time \(p\) and a truncation at time \(T\) can be written as: \[ \begin{aligned} \mathbb{P}(S = P + \tau \, | \, P = p, S < T) &= \frac{\mathbb{P}(S = P + \tau, P = p, S < T)}{\mathbb{P}(P = p, S < T)} \\ &= \frac{\mathbb{P}(S = P + \tau < T \, | \, P = p)\mathbb{P}(P = p)}{\mathbb{P}(S < T \, | \, P = p)\mathbb{P}(P = p)} \\ &= \frac{\mathbb{P}(S = P + \tau < T \, | \, P = p)}{\mathbb{P}(S < T \, | \, P = p)} \\ &= \frac{f_p(\tau)}{\int_0^{T - p} f_p(x) \text{d}x} = \frac{f_p(\tau)}{F_p(T - p)}, \quad p + \tau < T. \end{aligned} \]
Examining this equation illustrates that (right) truncation renormalises the density by the values which are possible. For example, if the distribution \(x \sim \text{Unif}(0, 1)\) were right truncated by \(T = 0.5\) then \(x \sim \text{Unif}(0, 0.5)\).
1.3 Interval censoring
The exact timing of epidemiological events is often unknown. Instead, we may only know that the event happened within a certain interval. We refer to this as interval censoring. A very common example of interval censoring in epidemiology is date censoring, where we only know, or are using, data to the day of an event rather than the precise time. Other forms of interval censoring, like weekly or monthly interval censoring, are also common. When both primary and secondary events are interval censored, this is referred to as double censoring.
Mathematically single interval censoring is defined as follows. Assume the secondary event \(S\) is censored and so we don’t know when the event exactly happened. Instead, we only know that the secondary event happened between \(S_L\) and \(S_R\). Then, \[ \begin{aligned} \mathbb{P}(S_L < S < S_R \, | \, P = p) &= \int_{S_L}^{S_R} f_p(y-p) dy\\ &= F_p(S_R-p) - F_p(S_L-p) \end{aligned} \]
Similarly, double interval censoring is defined as follows. Now, assume that both the primary \(P\) and secondary \(S\) events are truncated. We only know that the primary event happened between \(P_L\) and \(P_R\) and the secondary event happened between \(S_L\) and \(S_R\). We now write \(g_P\) to denote the unconditional distribution of primary events. Then, \[ \begin{aligned} \mathbb{P}(S_L < S < S_R \, | \, P_L < P < P_R) &= \mathbb{P}(P_L < P < P_R, S_L < S < S_R \, | \, P_L < P < P_R)\\ &= \frac{\mathbb{P}(P_L < P < P_R, S_L < S < S_R)}{\mathbb{P}(P_L < P < P_R)}\\ &= \frac{\int_{P_L}^{P_R} \int_{S_L}^{S_R} g_P(x) f_x(y-x) \,dy\, dx}{\int_{P_L}^{P_R} g_P(z)\, dz }\\ &= \int_{P_L}^{P_R} \int_{S_L}^{S_R} g_P(x\,|\,P_L,P_R) f_x(y-x) \,dy\, dx \end{aligned} \] where \(g_P(x\,|\,P_L,P_R)\) represents the conditional distribution of primary event given lower \(P_L\) and upper \(P_R\) bounds.
2 The naive model
The simplest approach to modelling epidemiological delay distributions is ignoring truncation and censoring biases and simply treating the delays as continuous fully observed data. Then, the likelihood of observing a delay \(\mathbf{Y}_i\) given parameter \(\mathbf{\theta}\) is straightforward: \[ \mathcal{L}(\mathbf{Y}_i \, | \, \mathbf{\theta}) = f(y_i - x_i). \] where \(y_i\) and \(x_i\) are the observed primary and secondary event times.
As shown in (Park et al. 2024) when the data is double censored this modelling approach biases the mean by approximately a day as well as the standard deviation. Where right truncation is also present biases can be more severe with plausible simulated scenarios leading to biased means that were >30% shorter than the true distributions.
3 The latent model
This approach aims to account for the right truncation and double censoring using a generative modelling approach.
For each event, a latent variable is used to represent the exact time of the event.
This then allows the modelling of the continuous distribution, adjusted for the right truncation.
Whilst this is an approximation (Park et al. 2024) showed good recovery of simulated distributions in a range of settings.
However, the use of two latent variables per observed delay means that this approach may scale poorly with larger datasets.
That being said this approach has been used successfully in multiple real-world outbreak settings ((Ward et al. 2022)).
If using the latent model, please cite Park et al. (2024) in addition to epidist.
Mathematically this model is described as follows. We look at the conditional probability that the secondary event \(S\) falls between \(S_L\) and \(S_R\), given that the primary event \(P\) falls between \(P_L\) and \(P_R\) and that the secondary event \(S\) occurs before the truncation time \(T\): \[ \begin{aligned} &\mathbb{P}(S_L < S < S_R \, | \, P_L < P < P_R, S<T)\\ &= \mathbb{P}(P_L < P < P_R, S_L < S < S_R, S<T \, | \, P_L < P < P_R, S<T)\\ &= \frac{\mathbb{P}(P_L < P < P_R, S_L < S < S_R, S<T)}{\mathbb{P}(P_L < P < P_R, S<T)}\\ &= \frac{\int_{P_L}^{P_R} \int_{S_L}^{S_R} g_P(x) f_x(y-x) \,dy\, dx}{\int_{P_L}^{P_R} \int_{z}^T g_P(z) f_z(w-z) \, dz \,dw }\\ &= \frac{\int_{P_L}^{P_R} \int_{S_L}^{S_R} g_P(x) f_x(y-x) \,dy\, dx}{\int_{P_L}^{P_R} g_P(z) F_z(T-z) \,dw }\\ &= \frac{\int_{P_L}^{P_R} \int_{S_L}^{S_R} g_P(x|P_L, P_R) f_x(y-x) \,dy\, dx}{\int_{P_L}^{P_R} g_P(z|P_L, P_R) F_z(T-z) \,dw }\\ \end{aligned} \] Using latent variables, we can now rewrite the observation likelihood as: \[ \begin{aligned} x_i &\sim g_P(x_i \, | \, p_{L, i}, p_{R, i}) \\ y_i &\sim \text{Unif}(s_{L, i}, s_{R, i}) \\ \mathcal{L}(\mathbf{Y} \, | \, \mathbf{\theta}) &= \prod_i \left[ \frac{f_{x_i}(y_i - x_i)}{\int_{P_{L, i}}^{P_{R, i}} g_P(z \, | \, p_{L, i}, p_{R, i}) F_z(T - z) \text{d}z} \right]. \end{aligned} \] As before, \(g_P(z \, | \, p_{L, i}, p_{R, i})\) represents the conditional distribution of the primary event given lower \(P_L\) and upper \(P_R\) bounds; this is equivalent to modelling the incidence in primary events.
3.1 Bounding the latent primary event time
epidist samples both latent offsets on the unit scale.
Write \(w_{P, i} = p_{R, i} - p_{L, i}\) and \(w_{S, i} = s_{R, i} - s_{L, i}\) for the two censoring window widths.
Then
\[
\tilde{p}_i \sim \text{Unif}(0, 1), \qquad
\tilde{s}_i \sim \text{Unif}(0, 1), \qquad
s_i = w_{S, i}\, \tilde{s}_i,
\]
and the primary event offset is ordinarily
\[
p_i = w_{P, i}\, \tilde{p}_i .
\]
When the two censoring windows overlap the delay must still be non-negative. The primary event therefore has to precede the sampled secondary event. The offset is then bounded by the sampled secondary offset rather than by the window width. \[ p_i = s_i\, \tilde{p}_i . \]
That upper bound is itself a parameter. The map \(\tilde{p}_i \mapsto p_i\) therefore has Jacobian determinant \[ \left| \frac{\partial p_i}{\partial \tilde{p}_i} \right| = s_i , \] Stan does not add this term for a transformation written this way. See the Stan user’s guide on changes of variables. The log density of these observations needs it explicitly. \[ \log \mathcal{L}_i \;\longmapsto\; \log \mathcal{L}_i + \log s_i . \] The non-overlapping case needs no such term because \(w_{P, i}\) is data, so its Jacobian is a constant that cancels. The primary event density is used unnormalised here. A constant normalising it over \([0, s_i]\) would itself depend on \(s_i\).
3.2 A non-uniform primary event
Incidence is not flat within a reporting window during rapid growth or decline.
Setting primary = "expgrowth" replaces the flat primary event with an exponentially tilted one at rate \(r\).
The tilt is placed on the unit scale offset, with the rate scaled by that offset’s bound \(b_i\), which is \(w_{P, i}\) ordinarily and \(s_i\) where the windows overlap. \[ \tilde{p}_i \sim \text{ExpGrowth}(0, 1, r_i b_i). \] Unlike the flat case, this density has to be normalised. Its constant depends on \(r_i\), which is estimated, so dropping it would let the likelihood grow without bound in \(r_i\).
\(r\) is a distributional parameter, so it takes a formula and a prior like any other. This allows the rate to vary by covariate. The delays carry little information about \(r\). It is normally taken from a separate estimate of epidemic growth and given an informative prior centred on that value, as in Brand et al. (2026), rather than learned from the delays.
This follows the implementation in Brand et al. (2026).
The latent model formulation follows Ward et al. (2022). The need for this adjustment was identified in Funk and Abbott (2026) and derived in Brand et al. (2026).
4 The marginal model
The marginal model corrects for the same biases as the latent model but integrates out the exact event times numerically, or analytically where closed-form solutions exist, rather than sampling latent variables.
This approach uses the primary event censored distribution implemented in the primarycensored package (Abbott et al. 2025).
If using the marginal model, please cite primarycensored in addition to epidist.
Under the assumption that the forward distribution does not change within the censoring interval (i.e. \(f_x = f\) for \(x \in [P_L, P_R]\)), the double censoring probability from Section 1.3 simplifies to
\[
\mathbb{P}(S_L < S < S_R \mid P_L < P < P_R) = \int_{P_L}^{P_R} g_P(x \mid P_L, P_R) \left[F(S_R - x) - F(S_L - x)\right] \text{d}x.
\]
For common delay and primary event distributions, such as gamma, lognormal, Weibull or generalised gamma delays with uniform primary events, primarycensored provides closed-form analytical solutions to this integral.
For other combinations, numerical integration is used.
The generalised gamma delay is provided by the gengamma() family of epidist, since brms has no such family, and it contains the gamma and Weibull families as special cases.
Right truncation at time \(T\) is handled by normalising the likelihood as in the latent model: \[ \mathcal{L}(\mathbf{Y} \mid \mathbf{\theta}) = \prod_i \frac{\mathbb{P}(S_{L,i} < S_i < S_{R,i} \mid P_{L,i} < P_i < P_{R,i})}{\int_{P_{L,i}}^{P_{R,i}} g_P(z \mid p_{L,i}, p_{R,i}) F(T - z) \, \text{d}z}. \]
Removing the latent variables reduces the number of parameters that must be sampled, and where analytical solutions exist the likelihood can be evaluated without numerical integration. In addition, identical observations can be aggregated and the likelihood computed once per unique combination of delay, censoring windows, and covariates. Together these make the marginal model substantially more efficient than the latent model, particularly for larger datasets with daily-censored data where many observations share the same structure.
For the mathematical details of primary event censored distributions, including the survival function derivation and closed-form solutions for specific distributions, see vignette("why-it-works", package = "primarycensored") and vignette("analytic-solutions", package = "primarycensored").
5 The meta model
Published delay estimates are usually summary statistics, and the estimation procedure behind them is itself a source of bias (Charniga et al. 2024; Park et al. 2024). The meta model fits one delay distribution to a mix of individual level data and published summaries in three steps.
- The latent delay is the forward distribution of Section 1, with density \(f(\cdot \, ; \theta)\) and distribution function \(F(\cdot \, ; \theta)\).
Its parameters \(\theta\) take a
brmsformula, so study level heterogeneity, for examplemu ~ 1 + (1 | study), is specified as for the other models. - Each study gets one representation of the delay, its estimand \(\tilde{F}\), the distribution its reported summaries estimate. For a study that estimated the delay without bias \(\tilde{F}\) is \(F\). Otherwise \(\tilde{F}\) is the result of applying the study’s procedure to \(F\): how it adjusted for censoring and right truncation, its censoring windows and the smallest delay it counted (Section 5.3).
- Every reported summary is read from \(\tilde{F}\), a mean from its mean, a standard deviation from its standard deviation and a quantile from its distribution function, with a sampling likelihood set by the study’s sample size or standard error (Section 5.1). Individual level rows use the likelihood of the marginal model (Section 4).
The sections below start from step 3, what each reported summary needs from \(\tilde{F}\), and then give how each kind of study gets there.
If using the meta model, please cite primarycensored in addition to epidist.
5.1 Sampling likelihoods
Write \(m_1\), \(\sigma\), \(\mu_3\) and \(\mu_4\) for the mean, standard deviation and third and fourth central moments of a study’s estimand \(\tilde{F}\), \(\kappa = \mu_4 / \sigma^4\) for its kurtosis and \(Q_p\) for its quantile at probability \(p\). Section 5.3 later gives these, and \(\tilde{F}\) itself, for each censoring code. Each likelihood below uses only them, so it is the same whatever the estimand. A study reports a value \(y\) computed from \(n\) delays.
5.1.1 A reported mean
A reported mean is approximately normal by the central limit theorem, \[ y \sim \text{Normal}\left(m_1, \; \text{se}_{m_1} \right), \tag{5.1} \] with \(\text{se}_{m_1}\) the reported standard error where given and \(\sigma / \sqrt{n}\) otherwise.
5.1.2 A reported standard deviation
A reported standard deviation is given a normal likelihood, a large sample approximation to its sampling distribution, with the kurtosis based standard error of a sample standard deviation, \[ y \sim \text{Normal}\left(\sigma, \; \text{se}_\sigma \right), \quad \text{se}_\sigma = \sigma \sqrt{\frac{\kappa - 1}{4 n}}, \tag{5.2} \] which follows from the asymptotic variance \((\mu_4 - \sigma^4)/n\) of the sample variance by the delta method (Cramér 1946). The usual \(\sigma / \sqrt{2(n-1)}\) assumes normal data, \(\kappa = 3\), and is too narrow for right skewed delays. Equations (5.2) and (5.10) are based on this approximation. When \(\sqrt{(\kappa - 1) / (4 n)}\) exceeds about a quarter it is inaccurate.
5.1.3 A reported quantile
A reported quantile \(y\) at probability \(p\) of a continuous estimand is fitted through the number of delays at or below it, which is binomial, \[ \text{round}(n p) \sim \text{Binomial}\left(n, \; \tilde{F}(y)\right). \tag{5.3} \]
A single quantile of integer day delays, codes 0 and 3, is a discrete statistic. “The median is 5 days” says that the empirical distribution function crossed one half between 4 and 5 days, that is \(N_{\le y - w_s} < \lceil n p \rceil \le N_{\le y}\) with \(N_{\le y}\) the number of delays at or below \(y\). It is fitted as the probability of that event, \[ P(N_{\le y} \ge k) - P(N_{\le y - w_s} \ge k), \quad k = \lceil n p \rceil, \quad N_{\le y} \sim \text{Binomial}\left(n, \tilde{F}_0(y)\right), \tag{5.4} \] where \(w_s\) is the reporting resolution and \(\tilde{F}_0\) the step distribution function of the discrete estimand of Section 5.3, before continuity correction. The information this carries saturates as \(n\) grows.
Posterior predictions for a quantile row are drawn from the normal approximation of Equation (5.3), \[ p \sim \text{Normal}\left(\tilde{F}(y), \; \text{se}_p\right), \tag{5.5} \] with \(\text{se}_p = \sqrt{p(1-p)/n}\) the binomial standard error of an empirical distribution function. A quantile supplied with a standard error \(\text{se}_y\) on the delay scale is fitted on that scale instead, \[ y \sim \text{Normal}\left(Q_p, \; \text{se}_y\right), \tag{5.6} \] because moving \(\text{se}_y\) onto the probability scale multiplies it by the density \(\tilde{f}(y)\), which is close to zero when \(y\) is far from the quantile the model implies.
Several quantiles from one study at probabilities \(p_1 < \dots < p_k\) with values \(y_1 \le \dots \le y_k\) cut the delay axis into \(k + 1\) cells, and the counts in them are multinomial, \[ (c_1, \dots, c_{k+1}) \sim \text{Multinomial}\left(n, \; \left(\tilde{F}(y_1), \; \tilde{F}(y_2) - \tilde{F}(y_1), \; \dots, \; 1 - \tilde{F}(y_k)\right)\right), \tag{5.7} \] with \(c_j = \text{round}(n p_j) - \text{round}(n p_{j-1})\) and \(c_{k+1} = n - \text{round}(n p_k)\). A single quantile reduces this to Equation (5.3). Two quantiles reported at the same value are merged into one cell. A cell whose probability underflows to zero while the study saw delays in it is floored at \(10^{-300}\).
Several quantiles of integer day delays from one study are fitted as the joint probability of the crossings they stand for. The counts \(N_{e_1} \le N_{e_2} \le \dots \le N_{e_m}\) at the integer edges the reported quantiles name, the day below each and the day itself, form a Markov chain of binomial steps, \[ N_{e_1} \sim \text{Binomial}\left(n, \; \tilde{F}_0(e_1)\right), \quad N_{e_{i+1}} \mid N_{e_i} \sim N_{e_i} + \text{Binomial}\left(n - N_{e_i}, \; \frac{\tilde{F}_0(e_{i+1}) - \tilde{F}_0(e_i)}{1 - \tilde{F}_0(e_i)}\right), \tag{5.8} \] and a quantile reported at \(p\) landing on day \(y\) puts the box \(N_{\le y - w_s} \le \lceil n p \rceil - 1\) and \(N_{\le y} \ge \lceil n p \rceil\) on two of them. The likelihood is the probability that every count fell in its box, \[ P\left(l_i \le N_{e_i} \le u_i, \; i = 1, \dots, m\right), \tag{5.9} \] computed by a forward pass over the counts, one binomial step per edge. Two quantiles reported at the same value are two constraints at one edge, and a single quantile reduces Equation (5.9) to Equation (5.4). It is used in place of the multinomial of Equation (5.7), which treats each quantile of integer day delays as a continuous cut point and so overstates what a large study reports.
5.1.4 Summaries from the same study
Summaries that one study computed from the same delays are correlated, so they are fitted jointly. Two summaries are fitted together when they agree on every field other than the summary itself, so a linear predictor cannot vary within the group. A summary supplied with its own standard error is fitted on its own.
A mean and a standard deviation from one study are given the asymptotic bivariate normal of the pair, \[ \begin{pmatrix} y_{m} \\ y_{\sigma} \end{pmatrix} \sim \text{Normal}\left( \begin{pmatrix} m_1 \\ \sigma \end{pmatrix}, \; \frac{1}{n}\begin{pmatrix} \sigma^2 & \mu_3 / (2 \sigma) \\ \mu_3 / (2 \sigma) & \sigma^2 (\kappa - 1) / 4 \end{pmatrix} \right), \tag{5.10} \] whose off diagonal is \(\text{Cov}(\bar{x}, s^2) = \mu_3 / n\) carried onto the standard deviation scale by the delta method (Cramér 1946). The correlation is \(\gamma_1 / \sqrt{\kappa - 1}\) with \(\gamma_1 = \mu_3 / \sigma^3\) the skewness, which every distribution keeps inside \([-1, 1]\). Moments taken from a grid or quadrature can sit just outside, so the correlation is clipped.
A study with a continuous estimand, codes 1, 2 and 4 of Section 5.3, that reports a mean or a standard deviation alongside quantiles has all of them fitted as one multivariate normal, \[ \begin{pmatrix} y_m \\ y_\sigma \\ y_{p_1} \\ \vdots \\ y_{p_k} \end{pmatrix} \sim \text{Normal}\left( \begin{pmatrix} m_1 \\ \sigma \\ Q_{p_1} \\ \vdots \\ Q_{p_k} \end{pmatrix}, \; \frac{1}{n} \Sigma \right), \tag{5.11} \] whose mean and standard deviation block is that of Equation (5.10). The rest follows from the Bahadur representation of a sample quantile, \(\hat{Q}_p - Q_p \approx -(\hat{F}_n(Q_p) - p) / \tilde{f}(Q_p)\) with \(\hat{F}_n\) the empirical distribution function of the study’s delays and \(\tilde{f}\) the density of \(\tilde{F}\) (Bahadur 1966), \[ \Sigma_{q_i q_j} = \frac{p_i (1 - p_j)}{\tilde{f}(Q_{p_i}) \, \tilde{f}(Q_{p_j})}, \quad p_i \le p_j, \qquad \Sigma_{m q_i} = -\frac{\int_L^{Q_{p_i}} (x - m_1) \, \text{d}\tilde{F}(x)}{\tilde{f}(Q_{p_i})}, \qquad \Sigma_{\sigma q_i} = -\frac{\int_L^{Q_{p_i}} \left((x - m_1)^2 - \sigma^2\right) \text{d}\tilde{F}(x)}{2 \sigma \, \tilde{f}(Q_{p_i})}, \tag{5.12} \] with \(L\) the smallest delay the study counted, and the last entry carried from the sample variance to the standard deviation by the delta method. The partial moments are integrated by parts over the nodes on which \(\tilde{F}\) is evaluated, and \(\tilde{f}(Q_p)\) is the closed form density of the estimand where its quantile is refined and the slope of \(\tilde{F}\) between the bracketing nodes where it is not, see Section 5.3.1. Fitting the two kinds separately would count the information they share twice. A continuous estimand reporting only a mean and a standard deviation, or only quantiles, keeps Equation (5.10) or (5.7), which the joint normal reduces to and which is exact for a single quantile. On the grid of codes 0 and 3 the quantiles are discrete statistics fitted by Equations (5.4) and (5.9), so the mean and standard deviation of such a study are fitted separately from its quantiles, and a study reporting both kinds is over weighted in the same way. Keep its mean and standard deviation and drop its quantiles.
5.1.5 A vector of summaries with a covariance
A study that fitted a distribution to its delays can publish draws of the fitted parameters, which as_epidist_multivariate() summarises by their mean and covariance and as_epidist_estimates_data() pushes through to the summaries the fitted distribution implies.
This is the reporting format we recommend, because it keeps the correlation between the reported quantities.
With \(y\) the reported vector and \(\Sigma\) the covariance over it,
\[
y \sim \text{Normal}\left(m(\theta), \; \Sigma\right), \tag{5.13}
\]
where \(m(\theta)\) holds \(m_1\) for a mean, \(\sigma\) for a standard deviation and \(Q_p\) for a quantile.
Summaries of a \(k\) parameter fit are functions of \(k\) numbers, so at most \(k\) may be reported with a covariance or standard errors, and asking for more is an error.
5.1.6 Reported distribution parameters
A study that published the parameters of a fitted distribution \(\hat{F}\) has them converted to summaries by epidist_estimates_parameters(), so the family it fitted need not match the family fitted to it.
The summaries are taken over the range of delays the study could have seen, conditioning \(\hat{F}\) on \((L, D]\),
\[
\hat{F}_{L,D}(y) = \frac{\hat{F}(y) - \hat{F}(L)}{\hat{F}(D) - \hat{F}(L)}, \quad L < y \le D, \tag{5.14}
\]
with \(L\) the smallest delay it counted and \(D\) its observation time, or \(D = \infty\) for a study that adjusted for right truncation.
Without this a study that did not correct for right truncation is charged with tail spread its data never had.
Reported parameter standard errors are carried onto the summaries by the delta method, as \(J V J^\top\) with \(V\) the diagonal matrix of squared standard errors and \(J\) the Jacobian of the map from parameters to summaries, and fitted jointly through Equation (5.13).
For \(k\) summaries of a \(k\) parameter fit this gives back the curvature \(V^{-1}\), where fitting each with its own standard error would overstate their uncertainty.
A full parameter covariance is better passed as draws to as_epidist_multivariate().
A study with no parameter uncertainty falls back to the sample size likelihoods above, and any number of its summaries may be reported.
This route assumes the reported distribution has the shape of the study’s estimand, which holds only approximately where a continuous family was fitted to integer date differences.
There, quantiles in the body of the distribution are more reliable than a standard deviation, which depends on a tail the study never saw.
5.1.7 Individual level records
Individual level rows use the marginal model likelihood of Section 4 unchanged, with \(\theta\) shared with the summary rows.
Their primary event distribution is set with primary, so the tilted primary event of Section 3.2 is available to them.
A summary row whose growth rate is unknown, or reported with uncertainty, uses the same pgrowth parameter, see Section 5.3.3.
The joint likelihood is the product of the sampling likelihoods above over the summary rows and the marginal model likelihood over the individual level rows.
5.2 What we need from each study
as_epidist_estimates_data() takes, for each summary:
- how the study adjusted for interval censoring, a code from 0 to 4 defined in Section 5.3,
- whether it adjusted for right truncation, and if not its observation time, how collection stopped and the growth rate of primary events over the study period, or
NAwhere that is to be estimated, with a standard deviation where the study reported one, - its censoring windows \(w_p\) and \(w_s\), the widths of the intervals its primary and secondary events were observed in,
- its sample size, or a standard error or covariance in its place,
- the smallest delay it counted.
Both truncation designs assume the study sampled a cohort of primary events, so a study that sampled on the secondary event reports a backward distribution of Section 1.1, which is not represented here.
Systematic reviews rarely record this metadata, so it is usually the analyst’s judgement and should be stated and varied in a sensitivity analysis.
The checklist of Charniga et al. (2024) covers it.
Where it is missing entirely, a covariate for the phase of the outbreak makes the brms formula a meta-regression that estimates the residual bias instead.
5.3 The biased estimands
Let \(f(\tau; \theta)\) and \(F(\tau; \theta)\) be the forward density and distribution function of Section 1, and \(g_P\) the distribution of a primary event within its window of width \(w_p\), uniform under constant incidence and exponentially tilted towards more recent times otherwise.
Write
\[
F_{pc}(\tau; \theta) = \mathbb{P}(\tau^\star + U \le \tau), \quad \tau^\star \sim f(\cdot \, ; \theta), \; U \sim g_P,
\]
for the primary event censored distribution function of primarycensored (Abbott et al. 2025).
Each censoring code below defines a distribution, right truncated at the observation time \(D\), with \(D = \infty\) for a study that adjusted for right truncation itself.
\(D\) is on the delay scale, the truncation time \(T\) of Section 1.2 less the time of the primary event.
This distribution is the study’s estimand \(\tilde{F}\), and its moments, distribution function and quantiles are what the likelihoods of Section 5.1 read.
Four moments are needed, because the kurtosis sets the sampling error of a reported standard deviation and the skewness its correlation with a reported mean.
Codes 0 and 3 give a discrete estimand on the study’s reporting grid and codes 1, 2 and 4 a continuous one, which is why the likelihoods for quantiles and for joint summaries differ between them.
5.3.1 Censoring adjustment
cens_adjusted records how a study handled interval censoring, one of five codes.
| Code | What the study did | Estimand |
|---|---|---|
| 0 | summarised integer date differences directly, that is summary statistics of the raw data | discrete, on the reporting grid |
| 1 | adjusted for both intervals, for example with a double interval censored likelihood | continuous, the delay itself |
| 2 | adjusted the secondary interval only, assuming a uniform delay within it | continuous, the delay plus the primary offset |
| 3 | assigned each delay to the centre of its interval | discrete, the code 0 grid moved up |
| 4 | placed the primary event at the midpoint of its window and integrated the secondary interval | continuous, code 2 moved down |
Each code below gives the estimand’s moments, \(\mathbb{E}[\tau^k]\) for \(k = 1, \dots, 4\), from which the mean, variance, skewness and kurtosis follow, its distribution function \(\tilde{F}\) and its quantiles \(Q_p\). \(\tilde{F}\) is computed at a set of points \(t_1 < \dots < t_m\), the grid cells for codes 0 and 3 and quadrature nodes otherwise. Every quantile starts from linear interpolation between the two points that bracket \(p\), \[ Q_p^{(0)} = t_i + \frac{p - \tilde{F}(t_i)}{\tilde{F}(t_{i+1}) - \tilde{F}(t_i)} \left(t_{i+1} - t_i\right), \quad \tilde{F}(t_i) \le p < \tilde{F}(t_{i+1}), \tag{5.15} \] which is exact on a grid and otherwise only as accurate as the spacing of the points. Where \(\tilde{F}\) and its density \(\tilde{f}\) have a closed form it is refined by two Newton steps, \[ Q_p^{(k + 1)} = Q_p^{(k)} - \frac{\tilde{F}(Q_p^{(k)}) - p}{\tilde{f}(Q_p^{(k)})}. \tag{5.16} \] A fixed number of steps from a fixed start keeps \(Q_p\) differentiable in \(\theta\), which a root search would not.
5.3.1.1 Code 0: summary statistics of the raw data
Code 0 is the code for summary statistics computed from the raw data, because recorded delays are always censored to the day or coarser. Its estimand is discrete, with probability \(q_j\) on the \(j\)th secondary window, \[ q_j = \frac{F_{pc}(j w_s; \theta) - F_{pc}((j-1) w_s; \theta)}{F_{pc}(w_s \lfloor D / w_s \rfloor; \theta)}, \quad j = 1, \dots, \lfloor D / w_s \rfloor, \tag{5.17} \] where bin \(j\) carries the delay \((j-1) w_s\) and the truncation point is discretised to the last full grid boundary. This is the doubly interval censored, right truncated probability mass function the marginal model of Section 4 uses for individual observations.
Moments. Sums over the grid.
Distribution function. The cumulative sum of the grid gives \(\tilde{F}_0\). For quantiles \(\tilde{F}\) is continuity corrected by interpolating \(\tilde{F}_0\) linearly through the mid points of its cells, because a reported quantile of day resolution data otherwise lands on a jump.
Quantiles. Equation (5.15) on the continuity corrected \(\tilde{F}\) is exact. The reported value is itself rounded to the grid, so a bias remains that does not shrink with \(n\).
5.3.1.2 Code 1: full adjustment
A fully adjusted study estimated \(f(\cdot \, ; \theta)\) itself, right truncated at \(D\).
Moments.
\[
\mathbb{E}[\tau^k \mid \tau \le D] = \frac{\int_0^D k t^{k-1} \left(F(D; \theta) - F(t; \theta)\right) \text{d}t}{F(D; \theta)}, \quad k = 1, \dots, 4, \tag{5.18}
\]
evaluated by Simpson’s rule, which reduces to the family moments when \(D = \infty\).
Simpson’s rule integrates numerically over an even number of equal intervals between \(L\) and \(D\), fitting a parabola to each pair.
Each study gets enough intervals to resolve the spread it reported, at least 100 or the value of options(epidist.meta_n_quad).
Distribution function. \(F(y; \theta) / F(D; \theta)\).
Quantiles. For a lognormal or Weibull delay without an accrual design the quantile has a closed form, \[ Q_p = F^{-1}\left(F(L) + p \left(F(D) - F(L)\right)\right), \] with \(L\) the smallest delay the study counted, see Section 5.3.4. Other families take the Newton steps of Equation (5.16).
5.3.1.3 Code 2: the uniform single interval approximation
A study that left the primary interval uncorrected observed \(\tau^\star + U\), right truncated at \(D\).
Moments. \[ \mathbb{E}[(\tau^\star + U)^k \mid \tau^\star + U \le D] = \frac{\int_0^D k t^{k-1} \left(F_{pc}(D; \theta) - F_{pc}(t; \theta)\right) \text{d}t}{F_{pc}(D; \theta)}, \quad k = 1, \dots, 4, \tag{5.19} \] evaluated by Simpson’s rule. Where \(D = \infty\) and the primary event is uniform the convolution is exact, \[ \mu_{pc} = \mu + \frac{w_p}{2}, \quad \sigma_{pc}^2 = \sigma^2 + \frac{w_p^2}{12}, \quad \mu_{4,pc} = \mu_4 + 6 \sigma^2 \frac{w_p^2}{12} + \frac{w_p^4}{80}, \tag{5.20} \] with \(\mu\), \(\sigma^2\) and \(\mu_4\) the mean, variance and fourth central moment of \(f(\cdot \, ; \theta)\).
Distribution function. \(F_{pc}(y; \theta) / F_{pc}(D; \theta)\).
Quantiles. With a uniform primary event \(\tilde{F}\) and \(\tilde{f}\) have a closed form, so the quantile takes the Newton steps of Equation (5.16). With a growing primary event they do not, and the quantile keeps Equation (5.15).
5.3.1.4 Codes 3 and 4: midpoint imputation
The estimand of code 3 is the grid of Equation (5.17) moved up by \(w_s / 2\). The estimand of code 4 is that of code 2 moved down by \(w_p / 2\).
Moments. A shift changes the mean alone, by \(+ w_s / 2\) for code 3 and \(- w_p / 2\) for code 4, and leaves the central moments as they are.
Distribution function. Shifted with the estimand.
Quantiles. Shifted with the estimand, from those of code 0 for code 3 and those of code 2 for code 4.
Under a uniform primary event code 4 therefore has the mean of code 1 and the variance of code 2 before truncation, because midpointing removes the mean of \(U\) but not its spread. The mirror reading, a midpointed secondary event and an integrated primary interval, has variance \(\sigma^2 + w_s^2 / 12\) and is not used, because the literature midpoints the wide exposure window of the primary event. A study that midpointed the secondary interval and left the primary alone is code 3. Each code integrates the interval it did not midpoint rather than drawing a random position in it, which would add \(w^2 / 6\) to the variance.
5.3.2 Right truncation
The truncation above conditions on the delay falling below one cutoff, which is what a cohort followed for a common observation time gives, trunc_design = "cohort".
This is truncation rather than right censoring, where a case is known to exist and contributes a survival term rather than being absent from the study.
Right censoring is not yet supported.
A study that accrued primary events over a window of length \(A\) and stopped at a calendar date, trunc_design = "accrual", saw a delay \(d\) only for primary events at least \(d\) before the stop.
With primary events arriving at a rate proportional to \(\exp(r t)\) the follow up available is
\[
w(d) = \int_0^{A - d} \exp(r t) \, \text{d}t = \frac{\exp(r (A - d)) - 1}{r}, \quad 0 \le d \le A, \tag{5.21}
\]
which tends to \(A - d\) as \(r\) tends to zero and, for a long window and a growing epidemic, to an exponential tilt by \(\exp(-r d)\).
This is the dynamical bias of Park et al. (2024).
The estimand is \(f(d; \theta) w(d)\) renormalised over \([0, A]\), so \(A\) replaces the cohort cutoff \(D\).
The weight multiplies the quadrature for Equations (5.18) and (5.19) at each node and renormalises.
For Equation (5.19) it is evaluated at \(d - w_p / 2\), to average the primary offset back out.
This is an approximation, because the study knows each primary event only to within its window, and it worsens as \(r\) grows and as \(w_p\) approaches \(A\).
The correction is only as good as the growth rate supplied.
Its quantiles keep Equation (5.15), because the estimand is defined by the interpolation between its points.
On the grid of Equation (5.17) the weight cannot be applied delay by delay, because a study that recorded dates knows each primary event only to its window. A delay is seen only if the whole primary window it came from was open long enough, so the weight is a step function of the delay that drops by one window’s growth weighted mass at each point \(A - j w_p\). Each grid cell is split at those points and each piece weighted by the mass left. When \(A\) is not a whole number of primary windows the last one is partial, of length \(l = A - w_p \lfloor A / w_p \rfloor\), and its cases are added separately, following \(F_{pc}\) with a window of \(l\). This keeps the grid exact for any \(A\), \(w_p\) and \(w_s\).
5.3.3 An estimated growth rate
The rate \(r_j\) of study \(j\) sets the tilt of \(g_P\), which every code but code 1 uses.
It also sets the weight of Equation (5.21).
Where the study’s growth_rate is NA, or is given with a growth_rate_sd, it is instead the distributional parameter pgrowth evaluated on the row, the parameter of Section 3.2, so it takes a brms formula and a prior and can be shared with individual level rows from the same outbreak.
Unless a formula is given for pgrowth the model uses pgrowth ~ 0 + study, one rate per study.
A study that reported a rate \(\hat{r}_j\) with a standard deviation \(s_j\) then gets
\[
r_j \sim \text{Normal}(\hat{r}_j, s_j^2), \tag{5.22}
\]
so the reported rate is a prior rather than a constant and its uncertainty reaches the delay.
Every other coefficient of pgrowth gets \(\text{Normal}(0, 0.25^2)\), which is weakly informative for a delay measured in days.
5.3.4 Left truncation
A study that only counted delays of at least \(L\) reported summaries conditioned on \(\tau > L\), the left truncation of survival analysis (Klein and Moeschberger 2003).
Each expression above reduces to its earlier form when \(L = 0\).
On the grid of Equation (5.17) the cells below \(L\) are dropped and the rest renormalised by their mass, which is \(F_{pc}(D) - F_{pc}(L)\) when \(L\) falls on a grid boundary.
For codes 3 and 4, \(L\) is on the reported scale, so the estimand they move is truncated at \(L\) moved back by the same shift, \(L - w_s / 2\) or \(L + w_p / 2\).
The truncated moments of Equation (5.18) pick up a boundary term,
\[
\mathbb{E}[\tau^k \mid L < \tau \le D] =
\frac{L^k \left(F(D; \theta) - F(L; \theta)\right) +
\int_L^D k t^{k-1} \left(F(D; \theta) - F(t; \theta)\right) \text{d}t}
{F(D; \theta) - F(L; \theta)}, \quad k = 1, \dots, 4, \tag{5.23}
\]
and Equation (5.19) likewise with \(F_{pc}\) in place of \(F\).
The distribution function becomes
\[
\tilde{F}(y) = \frac{F(y; \theta) - F(L; \theta)}{F(D; \theta) - F(L; \theta)}, \quad L < y \le D, \tag{5.24}
\]
zero at or below \(L\) and one above \(D\).
The accrual weight of Equation (5.21) is unchanged, since the cells and nodes it multiplies now start at \(L\).
Individual level rows pass \(L\) to primarycensored as their left truncation point, as the marginal model does.