Estimating epidemiological delay distributions for infectious diseases

Conceptual overview

For any epidemiological process, there is an underlying intrinsic distribution that describes the time difference between events in the process. In the context of epidemiological delays, the intrinsic distribution is a theoretical distribution that describes the probability of waiting a certain amount of time between a primary and a secondary event (e.g. between infection and symptom onset) under constant conditions in the population (Champredon and Dushoff, 2015). This distribution can be calculated by averaging across individuals: for example, the intrinsic generation-interval distribution only depends on the average infectiousness of infected individuals at a given time and does not depend on other factors, such as intervention measures or the proportion susceptible (Park et al., 2021). As the intrinsic distribution characterises the underlying infection characteristic, it is used for modeling, and is therefore often the distribution that we want to estimate. However, the intrinsic distribution is generally not directly observable, as it may differ from the realised distributions that are measured from the actual primary and secondary events observed during an epidemic. In particular, realised distributions can change over the course of an epidemic: for example, changes in realised generation intervals can reflect changes in transmission dynamics, including susceptible depletion. For practitioners, characterizing both the intrinsic and realised distribution is important. In most cases, the primary event always occurs before the secondary event (e.g., infection followed by symptom onset). However, there are exceptions: for example, an infectee may develop symptoms before their infector (Svensson, 2007). Negative delays can also happen for single-individual events: an infected individual may test positive before or after symptom onset (Singanayagam et al., 2020). Here we focus on non-negative, single-individual delays, i.e. where the secondary event always happens after the primary event. We also focus on the case where forward distributions do not change over the course of an epidemic—later in the Ebola virus disease example, we show that this is not always the case. Many of the lessons in this paper will however carry over to more complicated cases involving negative delays and/or delays between epidemiological events from two individuals.

We take the reporting delay—defined as the time from symptom onset (the primary event) to confirmation (the secondary event)—as an example. There are two different approaches to measuring this interval. If we group individuals based on their time of symptom onset and follow them until confirmation we are measuring forward delay. If we instead group individuals based on when they tested positive and ask when they developed symptoms, we are measuring backward delays. For any individual realisation of a delay as a paired primary and secondary events, the length of forward and backward delays are identical. However, in a cohort of individuals that experienced the primary (forward) or secondary (backward) events at the same time, the resulting forward and backward distributions can systematically differ when incidence is changing over time. Taking the cohort approach also allows us to ask how the forward and backward distributions change over the course of an epidemic.

As an example, consider a population where the incidence of infection is growing exponentially (Fig. 1A). If we take a cohort of individuals who developed symptoms on day 25, then we observe Fig. 1B, which corresponds to a forward distribution. In this case, we see that the forward distribution in Fig. 1B matches the intrinsic distribution that we used to simulate Fig. 1A. In general, we expect the forward distribution to approximate the intrinsic distribution reasonably well, in particular under relatively stable external conditions (Park et al., 2021). This means that the forward distribution (rather than the backward distribution) is most useful for downstream analysis, modeling, and decision support.

• Download figure
• Open in new tab

Figure 1: Schematic diagrams showing forward and backward distributions with the impact of truncation and censoring biases.

(A) The density plot of exponentially increasing cumulative cases of primary events (orange) and the corresponding secondary events (purple). Case IDs (y-axis) are ordered by the primary event time. The dashed vertical line indicates day 25. Panels B–E show four different delay distributions using day 25 as a reference day. The left panels show a sample data set with circles representing the timing of primary events and triangles representing the timing of secondary events. The right panels show the corresponding distributions (filled areas) against the intrinsic distribution (transparent areas, solid lines). (B) The empirical forward distribution measured from a cohort of individuals who experienced the primary event on day 25. (C) The empirical backward distribution measured from a cohort of individuals who experienced the secondary event on day 25. (D) The truncated distribution measured from all individuals who experienced both the primary and secondary events before day 25 (i.e., those with solid lines connecting their events). (E) The observed distribution under daily reporting from all individuals who experienced the primary event before day 25 (without truncation). Empty circles and triangles represent the event dates under daily reporting.

Some delay distributions are observed when the secondary event is reported, and the timing of the primary event is then recalled or inferred (e.g. in a reporting delay, typically the confirmation is first observed and the date of symptom onset recalled by the patient). Since at this point, both primary and secondary events have already happened, backward-looking cohorts, and hence backward distributions, typically do not have truncation. For example, if we take a cohort of individuals who tested positive on day 25, we observe Fig. 1C which corresponds to a backward distribution. In this case, we see that the backward distribution has a shorter mean than both the intrinsic and forward distributions. This is because the backward-looking cohort contains a mixture of individuals who developed symptoms recently (shorter reporting delays) and individuals who developed symptoms further in the past (longer reporting delays). In the case of a growing epidemic, individuals who developed symptoms recently are more abundant, and so the backward distribution is shifted toward shorter intervals. Similarly, when the epidemic is declining, individuals who developed symptoms recently are rare, and so the backward distribution is shifted toward longer intervals (and has a longer mean than the intrinsic distribution) (Xin et al., 2021). More generally, the backward distribution always depends on the timevarying incidence of primary events. We refer to this dependence as dynamical bias because it biases the backward distribution in comparison to the forward distribution (Park et al., 2021).

Additionally, the realized forward and backward distributions are both susceptible to a range of biases. Both can depend on events that are censored and both may be truncated. In the case of forward-looking cohorts, the most common bias is right truncation as the data set may be finalized for analysis before all individuals are observed. Similarly, the backward distribution suffers from left truncation when data on primary events are not available before some time point (Cain et al., 2011).

In particular, we note that, in the exponential growth case, there is a theoretical equivalence between backward and truncated forward cohorts. This can be seen in Fig. 1D, which can be calculated by taking all individuals who were infected and developed symptoms before day 25 (therefore truncating all events that happened after day 25). This results from an equivalence of the backward and truncated cohorts in exponential growth settings. In the following section, we describe the effects of truncation (including its similarities and differences with the dynamical bias) and censoring in detail.

Mathematical definition

To sharpen our discussion of forward and backward distributions, we define them mathematically. As in the conceptual overview and Fig. 1, we often define cohorts based on time intervals (e.g., a group of individuals who were infected on the same day, same week, or since the start of an epidemic). Mathematically, we can instead define incidence, and therefore cohorts, at any given, infinitely precise, time point : for example, in the standard continuoustime Susceptible-Infected-Recovered (SIR) model, the incidence of infection at time t corresponds to βS(t)I(t), where β represents the transmission rate, S represents the number of susceptible individuals, and I represents the number of infected individuals.

Here, we write f_p(τ) to represent the conditional probability density of observing a forward delay of length τ given that the primary event occurred at time p. Whilst the forward distribution can vary, for example over time or by risk or age groups, we primarily focus on cases where the forward distribution remains stable over the course of an epidemic (i.e., f_p(τ) = f (τ) for all p in this paper), and is thus identical to the intrinsic distribution. We also write 𝒫(p) to denote the incidence of primary events (i.e., the rate at which individuals experience the primary event at time p). Likewise, we write b_s(τ) to represent the conditional probability density of observing a backward delay of length τ given that the secondary event occurred at time s. We also write 𝒮(s) to denote the incidence of secondary events (i.e., the rate at which individuals experience the secondary event at time s). Then, the total density of individuals 𝒯(p, s) who experienced the primary event at time p and secondary event at time s can be equivalently expressed in terms of the forward distribution: and the backward distribution: We can then write down a relationship between the forward and backward distributions using the above relationship: By substituting p = s− τ and rearranging, we see that the backward distribution is given by the incidence of primary events normalized by the size of the relevant cohort: Here, the denominator represents the normalization constant such that the probability distribution integrates to 1. In particular, when the incidence of primary events is growing exponentially at rate r and the forward distribution is static (i.e. f_p(τ) = f (τ)), the backward distribution is also static, and given by: This simplified relationship has been used extensively in the literature (Britton and Scalia Tomba, 2019; Verity et al., 2020b; Park et al., 2022). When r > 0, the backward distribution will have a shorter mean than the corresponding forward distribution because we are more likely to observe individuals who experienced the primary event recently. This dependence on the incidence of the primary event leads to dynamical bias when the backward distribution is used as an estimate of the intrinsic distribution (Park et al., 2022). As can be seen from Eq. (5), the forward and backward distributions will be equivalent when incidence is stable (r = 0).

Right truncation

Right truncation arises from the inability to observe secondary events that happen after the observation time T (see Fig. 1D). For example, during an ongoing epidemic, we tend to underestimate the infection fatality risk because we cannot observe all of the outcomes for infected individuals in real time due to the delay between infection and death (or recovery) and because we typically base observation on individuals having died (or recovered) (Ghani et al., 2005; Lipsitch et al., 2015). In general, we don’t know how much of the data is truncated until all outcomes have been fully observed (Gelman et al., 2013, p. 227). Since right truncation limits our ability to observe long intervals (i.e. those in the right tail of a distribution), it will generally bias us towards observing shorter delays (and therefore to underestimate the mean).

To define this mathematically, we introduce random variables P and S representing primary and secondary events. Then, the conditional probability of observing a forward delay of length τ (i.e., S = P + τ) given primary event time P = p and observation time T (i.e., S < T) can be written as: where f_p is the forward distribution and F_p is the forward cumulative probability distribution. Conceptually, Eq. (9) is saying that the probability of observing a given delay from a truncated distribution is equal to the probability of observing that delay from the untruncated forward distribution normalised by the probability of observing any delay before the end of observation. From a cohort perspective, when we measure forward delays from individuals who experienced primary events at time p < T, we can only observe delays that are shorter than T− p. This means that the amount of truncation is higher in a cohort of individuals who experienced their primary events more recently.

The population-level effect of right truncation will be exaggerated when an epidemic is growing because there will be more individuals who experienced primary events recently, meaning their secondary events may not have happened yet. When growth is exponential, the effect of right truncation is quantitatively equivalent to dynamical bias, as shown in Fig. 1C–D. In contrast, when the epidemic is declining quickly, fewer individuals experienced primary events recently, and the effect of right truncation can become negligible as few delays are unobserved. On the other hand, dynamical bias will still be an issue because a cohort of individuals who experienced secondary events during the decay phase will include a relatively large proportion of those who had their primary events further in the past, causing the backward distribution to have a longer mean than the forward distribution. As dynamical biases specifically refer to the dependence of backward distribution on epidemic dynamics, they need not be taken into account when analyzing the forward distribution as long as biases due to truncation are considered.

Interval censoring

Censoring occurs when the exact event time is unknown. In contrast to truncation, we know that the focal event happened in the case of interval censoring; we just do not know precisely when. For example, we know that a symptomatic individual was infected sometime before symptom onset but we usually don’t know when. Even if the exact date of the event is known, implicit censoring remains because the event could have theoretically occurred at any time within a 24-hour time window. This implicit censoring discretises the observed delay distribution (Fig. 1E). To define censoring of both primary and secondary events, which is known as double censoring (Reich et al., 2009), we start by first defining the simple single censoring case. Here the secondary event is censored between time S_L and S_R, and the primary event (p) is known. Then, we have the following relationship, Integrating this leaves This relationship is commonly used in practice to discretise delay distributions, particularly into daily intervals (Flaxman et al., 2020; Abbott et al., 2020). However, as we have conditioned on the primary event being known, this approach is not appropriate for the common case where P is also censored, for example, if only a date is known. The degree of bias this will introduce depends on event probability within the primary censoring interval.

For the more common setting where the primary event is also censored between time P_L and P_R the conditional probability that the secondary event occurs between time S_L and S_R can be written as: Here, g_P (x) is the probability distribution of the primary event time, and g_P (x| P_L, P_R) is the conditional probability distribution of the primary event time given its lower and upper bounds, As g_P (x) depends on the incidence of primary events, and therefore epidemic dynamics, this can be rewritten as, where P (x) is a continuous-time incidence, as explained earlier; in practice, accounting for this factor will require estimating continuous-time incidence from discrete-time incidence. Making the simplifying assumption of a fixed epidemic growth rate r, the conditional probability can then be rewritten as, The shape of the conditional probability distribution g_P (x|P_L, P_R) is shown in Fig. 2 for different epidemic growth rates.

• Download figure
• Open in new tab

Figure 2: Schematic diagrams showing different assumptions about interval censoring.

(A) The distribution of primary event times within a one-day censoring interval across different growth rates. Note that a sudden decrease in probability densities around 0 and 1 are numerical artifacts from plotting. (B) The corresponding distribution of weights for the interval-reduced censoring, which is a convolution between the distribution of event times of primary and secondary events within each censoring interval. (C) The distribution of weights for the interval-reduced censoring for different approximations.

Censoring and truncation

Finally, to understand the joint effects of right truncation and interval censoring on the forward distribution, we can quantify the conditional probability of the secondary event occurring between time S_L and S_R given the primary event occurring between P_L and P_R and the observation time T : By dividing both the numerator and denominator by , we have: Here, expressing the probability distribution of primary events g_P (x) in terms of the conditional probability distribution g_P (x|P_L, P_R) allows for a more flexible inferentialframework later on. We also see that the truncation problem also depends on the censoring problem: uncertainties in the primary event time further affect the amount of truncation. This framework generalises the work of Seaman et al. (2022) who reviewed various methods for dealing with right truncation bias but did not account for interval censoring.

2.3.1 Simplifying assumptions for the following methods

We only consider single-individual delays, where we have independent observations of event times across different individuals. This means that in this work, we do not consider methods that account for non-independence in generation or serial intervals. These pair delays also introduce additional biases due to uncertainties about who infected whom—we do not address these issues here. Nonetheless, the methods we present here can still be applied to adjust for truncation, censoring, and dynamical biases in estimating generation-and serial-interval distributions. Finally, we also assume that the forward distribution stays constant over time (i.e., f_p(τ) = f_forward(τ)). All of the methods we consider can be generalised to delays that vary in discrete time using our implementations.

2.3.2 Exact censoring and truncation method

Under daily reporting, the exact timing of the primary and secondary event times are unknown. Instead, the corresponding observed event times are P_L = ⌊P ⌋ and S_L = ⌊S⌋, respectively. Likewise, the corresponding upper bounds for the censoring interval are P_R = P_L + 1 and S_R = S_L + 1, respectively. Following Eq. (21), the joint likelihood can be written as: where θ represents the parameter vector, and Y represents the data vector.

From a practical standpoint, it is difficult to use this likelihood for inference, especially in real-time epidemic monitoring applications where estimation must be fast and robust. Solving double integrals analytically may be impossible, and calculating them numerically can be computationally costly and numerically unstable. Instead, we can implement this model using an equivalent Bayesian latent variable approach. Here we can treat the primary and secondary event times as latent variables (denoted x_i and y_i, respectively), conditional on the known upper and lower bound of these events for each individual, and then integrate across the uncertainty. Then, by taking away both integrals in the numerator, Eq. (22) can be re-written as: Here, we use the uniform distribution for the secondary event time y_i, so that y_i does not affect the likelihood, which allows Eq. (25) to be equivalent to Eq. (22). Using any other distribution for y_i will cause Eq. (25) to deviate from Eq. (22) and therefore will be incorrect.

In practice, modeling the conditional distribution of primary event time g_P (x| p_L,i, p_R,i) is expected to be a difficult problem as it depends on the changes in the incidence of primary events (Eq. (16)). However, if we assume that the incidence of the primary event is changing at a fixed growth rate r within the censoring interval we can further simplify the problem: Generalizing the problem beyond stable exponential growth within the interval requires a way of modeling realistic changes in incidence. Even with this simplifying assumption in most real-world settings this would require joint estimation of the growth rate of the primary events in order to properly propagate uncertainty; however, such a model is likely to be computationally more costly.

The exponential scenario provides insights for understanding the method introduced by Linton et al. (2020) during the beginning of the SARS-CoV-2 pandemic. Under Eq. (27), our likelihood is similar in form to the likelihood presented in Linton et al. (2020) with one major distinction (among several other minor differences). The conditional probability g_P (x p_L,i, p_R,i) we present here converges to a uniform distribution as r→ 0 as expected when incidence is stable (Fig. 2A). On the other hand, the corresponding component in Linton et al. (2020) is modeled as r exp(−ru)/(1 − exp(−ru)), where u is a variable of integration; this term does not converge as r → 0.

2.3.3 Approximate latent variable censoring and truncation method

When the primary event time is assumed to be uniformly distributed and censoring interval is narrow, the integral in the denominator of Eq. (25) can be approximated by F_forward(T x_i) for some p_L,i < x_i < p_R,i. Under these conditions, the exact method can be rewritten as (Ward et al., 2022): The likelihood for this method corresponds to the conditional probability under right truncation (Eq. (9)). This is a convenient approximation as it does not require additional integration. While it would be possible to extend this method to account for more complex prior distributions that capture epidemic growth or daily activities (e.g., circadian rhythms), the exact method (Eq. (25)) should be considered to avoid potential biases arising from the approximations.

2.3.4 Interval-reduced censoring and truncation methods

Without using latent variables, we can equivalently write down the approximate latent variable model as follows: where g_S(y_i|s_L,i, s_R,i) corresponds to the uniform distribution between s_L,i and s_R,i. When we rewrite in terms of the delay, d_i = y_i − x_i, we know that d_i should range between s_L,i − p_R,i and s_R,i − p_L,i. Further assuming F_forward(T − x_i) ≈ F_forward(T− p_R,i), we can write the above likelihood using a single integral: Here, w(d_i) is a convolution between g_P (x_i|p_L,i, p_R,i) and g_S(y_i||s_L,i, s_R,i): When both g_P (x_i|p_L,i, p_R,i) and g_S(d_i − x_i|s_L,i, s_R,i) are assumed to follow uniform distributions, w(d_i) will have a trapezoid shape (Fig. 2C), as discussed in Reich et al. (2009).

In Fig. 2B–C, we show how the shape of the weight function w(d_i) changes fordifferent epidemic growth rates (B) and for different assumptions (C). In particular, assuming a uniform distribution for both events censoring intervals approximates the interval-reduced censoring window well when the growth rate is small compared to the width of the censoring window and less well when this is not the case. In Fig. 8, we present the consequences of different assumptions about w(d_i) on the resulting continuous-time distributions.

We can further simplify the likelihood by assuming a uniform distribution for w(d_i) (Fig. 2C): This method is simple to implement as well as being computationally convenient and supported by many software packages, such as the brms package (Bürkner, 2018) in the R programming language (R Core Team, 2019).

However, it does not fully account for the censoring bias as the uniform assumption across the combined interval is not a good match for the prior that results from combining two uniforms or in more complex settings a growth-mediated prior and a uniform (Fig. 2B–C) or the impact of censoring on the amount of truncation due to the simplified denominator. This means that, depending on the amount of censoring and the shape of the delay distribution, this method will have a residual bias in both the mean and standard deviation of the estimated distribution (Fig. 8). In the special case when the primary event time is exactly known (therefore, p_L,i = p_i = p_R,i), this method would be accurate, and also equivalent to the latent variable method.

In the daily reporting setting, we consider in the rest of this paper, the likelihoodcan be reformulated in terms of the observed, discrete-time delay (d = S_L,i − P_L,i): Rephrasing the likelihood in this way highlights that the correct censoring interval for this approach is [d_i − 1, d_i + 1]. Instead, many researchers tend to rely on the [d_i, d_i + 1] interval (Flaxman et al., 2020; Abbott et al., 2020), which in effect only accounts for the censoring of a single event and hence induces additional bias beyond that introduced by assuming a uniform distribution. Other potential approximations are [d_i − 0.5, d_i + 0.5], which assumes the primary event is known to have taken place at the midpoint of the interval, and that all primary and secondary event times are known (He et al., 2020). The implications of these approximations are also explored in Fig. 8.

2.3.5 Interval-reduced censoring method

If we ignore the presence of right truncation, then Eq. (35) reduces to: which, again, assumes a uniform distribution for w(d_i). Reich et al. (2009) tested the consequences of this approximation and showed that assuming w(d_i) is uniform results in low coverage probabilities (defined as the proportion of confidence intervals that contain the true value across multiple simulations). However, the uniform assumption remains in regular use in the applied literature.

2.3.6 Truncation method

If we ignore the presence of censoring, Eq. (31) reduces to which can also be derived directly from Eq. (9) (Sun, 1995). Here, we use the lower bound for the primary event time to estimate the amount of truncation T− p_L,i as it represents the maximum degree of truncation. Other choices include using the upper bound T− p_R,i or the midpoint T − (p_L,i + p_R,i)/2. The use of an upper bound may be numerically less stable because the upper bound of primary event time can be equal to the observation time, in which case the denominator becomes zero.

2.3.7 Filtering method

A simple solution to mitigate truncation bias is to drop (or filter out) the most recent data (as these data will be most impacted by truncation) and then fit a probability distribution on the remaining observations. Given a filtering time t_filter, the likelihood can be written as: This approach can be readily extended to include double censoring as the filtering and distribution fitting steps are independent.

Unfortunately, whilst the filtering method may be convenient, it is sensitive to our choices of filtering time t_filter. If t_filter is too early, we might not have enough data to estimate the focal delay distribution with certainty. If t_filter is too late, truncation bias might be introduced. More broadly, dropping existing data is often not a good practice and at the very least is likely to increase the uncertainty of any estimated distribution. This can be a particular issue during epidemics when data can be sparse and as accurate as possible estimates are needed to inform decision-making. For simplicity, we assume t_filter = 10 days throughout the paper.

2.3.8 Naive method

A simple approach to estimating a delay distribution is to directly fit a probability distribution f_naive to the observed discrete delays. In this case, the likelihood is given by: where θ represents the parameter vector and Y represents the data vector. This method does not account for truncation or censoring biases but was historically used in practice due to convenience and lack of awareness of its limitations. An even more basic form of this approach is to simply report the empirically observed summary statistics such as the mean, standard deviation, or range —this approach is common in practice (Nolen et al., 2016; Gostic et al., 2020; Miura et al., 2022). The no interval method in Fig. 2 highlights that this approach does not approximate the underlying generative process well in the presence of censoring.

2.3.9 Discrete time methods for accounting for censoring and right truncation

So far, we have focused on fitting probability distributions to continuous forward delays directly by accounting for right truncation and censoring, presenting methods in order from more exact to most approximate. In this section, we outline discrete-time approaches to estimating epidemiological delay distributions whilst accounting for right truncation and censoring. These methods make the further simplifying assumption that censoring is daily for both primary and secondary events. To counterbalance this they both have a range of advantages to other methods in certain settings. In Fig. 9, we also present how different assumptions about the prior distributions within censoring intervals affect the discretisation of continuous-time distributions.

2.3.10 Overview table

Simulating infections

We explored both fixed growth rate and epidemic scenarios. For all scenarios, we explored 3 lognormal distributions representing “short” (log mean of 1.2 and log standard deviation of 0.4), “medium” (log mean of 1.6 and log standard deviation of 0.6), and “long” (log mean of 1.8 and log standard deviation of 0.8) time delays. These assumptions correspond to distributions with the following means (and standard deviations) on the natural scale: 3.6 (1.5) days, 5.9 (3.9) days, and 8.3 (7.9) days (Fig. 4B).

To simulate settings with a fixed growth rate, we used continuous-time stochastic simulations of exponential growth with a range of daily growth rate assumptions (assuming 10,000 individuals): “fast decay” (−0.2), “decay” (−0.1), “stable” (0), “growth” (0.1), and “fast growth” (0.2).

To simulate epidemic scenarios, we used a stochastic Susceptible-Infected-Recovered (SIR) model implemented using a Gillespie algorithm assuming an early exponential growth rate of 0.2 per day and a recovery rate of 1/7 per day (corresponding to a mean of 7 days) with 50 initial cases and a total population of 10,000.

For all scenarios, we took the primary event time to be the time of infection, and we simulated secondary events by drawing delays from the assumed continuous lognormal distribution for that scenario for each individual and then adding this delay to their primary event time. This simulation approach assumes no observation error for the primary event or secondary events beyond daily censoring and thus represents an idealised system.

Scenarios investigated

For the simulations with fixed growth rates, we assumed a sample size of 200 randomly drawn pairs of primary and secondary events with observation cut-off 30 days after the start of the simulation. We repeated each sampling step 20 times independently so that we ended up with 20 replicates of each simulated growth rate scenario.

For the epidemic wave simulation, we explored a range of sample sizes (10, 100, 200, and 400 samples) and observation time scenarios (including all observations up to 15 days, 30 days, 45 days, and 60 days from the start of the simulation).

Data

We downloaded line list data from Fang et al. (2016) which contained age, sex, symptom onset date, sample test date, the district of the case, and the Chiefdom of the case for Ebola virus disease (EVD) cases in Sierra Leone from May 2014 through September 2015. We then processed these data to keep only the date of symptom onset and the date of the sample test. We assumed that censoring for each of these dates was daily with a day defined as being from 12:00 AM to 11:59 PM.

Empirical context

We calculated various summary statistics using the available line list data, which we used to provide a context for comparisons of different inference methods. These statistics included changes in the mean forward and backward distributions, the empirical forward distribution at each observation cut-off date, as outlined in the ‘Scenarios Investigated’ section, in both real-time and retrospective settings; and, the proportion of secondary events that were unobserved for a cohort of individuals whose primary events took place within a rolling 60-day window.

Scenarios investigated

We estimated the delay from symptom onset to sample test using each method across four different observation windows of 60 days each (0–60 days, 60–120 days, 120–180 days, and 180–240 days after the first symptom onset). In particular, for each window, we only considered individuals who developed symptoms and also got tested within the time period to match a real-time analysis setting. Note that in this case study we do not account for additional potential delays from the sample test date to the reporting of these tests as these data were not available in Fang et al. (2016). However, our methods could naturally account for these additional delays if data were available. As a comparison, we then re-estimated the delay using each method for each observation period by including all individuals who developed symptoms within the observation period (regardless of when they got tested) and estimated the delay distribution. These retrospective estimates were then used to represent the “true” distribution when calculating the relative difference between estimates using data available in real-time and data available retrospectively for each method. Note that in this setting the “true” distribution may still be biased relative to the underlying forward distribution (e.g. methods that do not account for censoring may still be biased in both real-time and retrospective settings).

Summary

• More complex methods for estimating epidemiological delay distributions have increased computational requirements relative to simpler methods. In most settings, these computational requirements are expected to still be feasible for routine usage.

• The approximate-latent-variable method had exponential scaling in resource requirements when the sample size increased.

• The interval-reduced-censoring-and-truncation method required the least computational resources by an order of magnitude of the methods that accounted for both censoring and truncation.

• The dynamical adjustment method is unstable for short delays with larger sample sizes both with and without a truncated time series, though the instability increased when a truncated time series was used (i.e., a real-time one).

All the methods we considered are implementable with modest (i.e., laptop scale) computational hardware at the time of writing, and we are able to run our full analysis pipeline of several thousand model fits within this resource budget at a practical time scale (i.e, within several hours). However, more complex methods required greater computational resources. For example, the approximate latent variable, dynamic correction, and joint modeling approaches required an order of magnitude more resources for the same sample size than the interval-reduced-censoring-and-truncation method (Fig. 7A). The resource requirements for all models scaled with sample size, with the approximate latent variable model scaling the worst of the methods we explored. The dynamical correction approach had the highest variance in its computational requirements with some fits taking 10 to 100 times longer despite having the data used having the same sample size across fits. This was likely due to numerical instability from the integration step. It was a particular issue when a real-time incidence time series was used.

• Download figure
• Open in new tab

Figure 7:

(A) Distributions of run-times for model fitting using a 2019 16-core AMD Threadripper on the log 10 scale. Diamonds represent the overall mean run times and vertical lines represent the 5%, 35%, 65%, and 95% quantiles respectively. Whilst run times are specific to the hardware used, the relative differences between models and scenarios should be more readily generalisable to other hardware. (B) Percentage of model fits with more than 0.1% divergent transitions on the logit scale. (C) Percentage of divergent transitions for fits with more than 0.1% of samples that had divergent transitions on the logit scale.

In general, all the methods were numerically stable, excluding the dynamical correction approach in some settings, and converged. Methods that better captured the data-generating process (i.e. that accounted for known biases) were the most stable and had fewer diagnostic warnings. A notable exception to this was the dynamical model which was the only model to cause sampling in stan to fail completely in six instances and more generally had the highest proportion of fits with divergent transitions (Fig. 7B). These issues occurred most frequently in simulations with short delays and during periods of epidemic decline. In the case study, both retrospective and real-time fits failed at the 180-day observation point. Another driver of model issues was low sample sizes with the majority of instances of divergent transitions occurring in settings with a sample size of 10.