The doubling time analysis for modified infectious disease Richards model with applications to COVID-19 pandemic

1.   Introduction

The modified Richards model describes epidemic dynamics in two phases of fast and slow infection spread with a transition point (also referred to as turning or inflection point), at which the maximum incidence occurs. In the slow phase of infection spread (after the turning point), the epidemic peaks and subsequently declines, while the cumulative number of cases eventually saturates. At this point, little can be done in terms of control and prevention. On the other hand, efficient mitigation measures implemented at the early ascending stage of the outbreak have the potential of significantly reducing humanitarian cost and economic impacts of the emerging disease.

In this study, we employ the modified Richards model to analyze the growth of epidemic in terms of 1) the number of times cumulative cases double until the epidemic peaks and 2) the rate at which the intervals between consecutive doubling times increase during the emerging stage of the disease. It is our expectation that the doubling-time research will become a major tool in the study of infectious disease epidemics, and potentially in understanding the underlying processes by which emerging outbreaks quantify their spread in distinct ways. Furthermore, we believe that the doubling-time metrics for the modified Richards model could be used to systematically characterize recurrent epidemics (e.g., multiple influenza seasons) or spatial heterogeneity for epidemics that can be disaggregated at finer spatial resolutions (e.g., regions, districts).

In what follows, we combine theoretical evaluation of doubling times with rigorous numerical experiments using synthetic and real epidemiological observations stemming from the well-documented COVID-19 pandemic. These observations are crucial to test, guide and improve epidemic modeling. This data driven enhancement of the modified Richards model leads to substantial improvements in our ability to capture epidemic trajectories and to a better characterization of the transmission potential of infectious diseases. Importantly, our results underscore how various underlying assumptions in epidemic models influence the doubling-time mathematical analysis and its numerical estimation. More generally, our study highlights the importance of confronting crucial modeling assumptions with real epidemic data. These assumptions undoubtedly play a key role in projecting possible outcomes of an outbreak, and estimating the intensity of intervention steps required for minimizing its impact on the population. Understanding the complete effect of parameters quantifying these assumptions and their variations on epidemic doubling times will require further methodological advancements.

2.   Derivation of doubling times for modified Richards model in terms of cumulative number of cases

We analyze the growth of epidemics as a function of the number of times the cumulative incidence doubles before the epidemic peaks, and the rate at which the doubling time increases (hence the epidemic slows down), owing to interventions, behavior changes and depletion of susceptible individuals. The most difficult to control are epidemic threats that double in size many times by the time the peak has been reached, while doubling times increase at a relatively low rate. Conversely, the least difficult to control are epidemic threats that only double in size very few times, while the sequence of doubling times rapidly increases ^[1]. When $C'(t)$ represents the number of new infected cases at time $t$, the progression of the epidemic is often modeled by an autonomous differential equation tracking the cumulative number of infections, $C(t)$:

Given this equation, the doubling time, $\Delta t$, can be expressed in terms of the cumulative number of cases as follows

Epidemics that span over a short period of time relative to the average lifetime of the host population can be routinely described by S-shaped logistic curves ^[2,3,4,5,6]. These differential equations rely on a small number of parameters and assume initial exponential or sub-exponential growth phase that saturates as the number of cases accumulates, a pattern that implicitly captures a gradual decrease in at-risk susceptible population. In particular, the modified Richards model with four parameters has been used to fit a range of logistic-type epidemic curves ^[6,7,8,9,10]. For this model,

where $r > 0$ represents the intrinsic growth rate in the absence of any limitation to disease spread, $p$, $0\le p\le 1$, is a deceleration of growth parameter, $K > 0$ is the size of the epidemic, and $a > 0$ is a parameter that measures the extent of deviation from the S-shaped dynamics of the classical logistic growth model. At the early stages of the epidemic, this model enables us to capture different growth profiles ranging from constant incidence ($p = 0$), polynomial growth ($0 < p < 1$), to exponential growth ($p = 1$). As shown in ^[11], the solution of the indefinite integral in (2.2) for $F(s)$ defined in (2.3) can be expressed in terms of Gaussian Hypergeometric Function, which indicates that solution to (2.2)-(2.3) will converge under mild conditions. To derive the solution, let $u: = (\frac{s}{K})^{a}$ be a new variable. Then the doubling time, $\Delta t$, for the modified Richards model (2.1), (2.3) is given by

where $\delta: = 1 - \frac{1-p}{a}$ and $b: = \frac{K^{1-p}}{ar}$. Since the cumulative number of cases cannot exceed the threshold level, $K$, the new variable, $u$, is less than $1$. Thus, one can view $\frac{1}{1-u}$ as the sum of a geometric series and, therefore,

From (2.5) and the definitions of $b$ and $\delta$, one concludes that

3.   The sequence of doubling times

Consider the sequence of doubling times, $\{\Delta t_{j}\}$, for the modified Richards model (2.1), (2.3) such that

where $N$ is the finite number of times cumulative cases double until they reach half of the threshold level, $K/2$. Let $C(0) = C_0$. Then clearly,

and the $(j+1)$th element of the sequence, $\Delta t_{j+1}$, can be evaluated as

For $F(s)$ defined in (2.3), one concludes that

Identities (3.4) yield the following breakdown for $\Delta t_{j+1}$:

4.   Important special cases

In what follows, we consider some important special cases, for which the expression for the doubling time, $\Delta t$, admits a closed form representation. First, we explore the classical Richards model where

that is, a special case of (2.3) that corresponds to $p = 1$. According to (2.4), for $p = 1$,

or, alternatively, $\Delta t = \frac{1}{ar}\ln\left[\frac{(2K)^a-(2C)^a}{K^a-(2C)^a}\right]$. In (4.2), the term $\frac{\ln 2}{r}$ is the reflection of the exponential growth at the early ascending stage of the outbreak. The second term, $\frac{1}{ar}\ln\left\{\frac{1-\left(\frac{C}{K}\right)^a}{1-\left (\frac{2C}{K}\right)^a}\right\}$ is due to the contribution of the factor $1-\left(\frac{C}{K}\right)^a$ in the model, which, as opposed to a simple exponential growth, takes into account the carrying capacity of the environment. As expected, the expression for $\Delta t$ implies that the doubling of cases occurs until $C < \frac{K}{2}$. Taking into account (3.3) and (4.1), one arrives at the following equality for the sequence of doubling times, $\Delta t_{j+1}$, in the case of the classical Richards model:

Applying elementary formula $\ln(x) = \sum_{n = 1}^{\infty}\frac{(-1)^{n-1}(x-1)^n}{n}$, $0 < x\le 2$, with $x = 1-\left(\frac{2^jC_0}{K}\right)^a$ and then with $x = 1-\left (\frac{2^{j+1}C_0}{K}\right)^a$, one gets

Thus, for the classical Richards model (2.1), (4.1), one obtains

Comparing identities (3.5) and (4.5), one can easily notice that the second term in (4.5) coincides with the second term in (3.5) when $p = 1$. At the same time, the first term in (4.5), $\frac{\ln 2}{r}$, is equal to $\lim_{p\to 1}\frac{(2^jC_0)^{1-p}\left[2^{1-p}-1\right]}{r(1-p)}$. In other words, it is the limit (as $p$ approaches $1$) of the first term in (3.5).

Another special case that also gives rise to a Bernoulli differential equation is the one where

For $F(s)$ defined in (4.6) and $0 < p < 1$, one obtains a wide spectrum of early growth rates. At the same time, just like in the case of the classical Richards model, the doubling time here is a simple elementary function. From (2.4) one has

and, therefore,

We now consider the generalized-growth model, that is, $a \rightarrow \infty$ in (2.3), which results in

Similar to (2.3) and (4.6), the generalized-growth model enables us to capture the diversity of epidemic growth profiles during the first few generations of disease transmission, including constant ($p = 0$), polynomial ($0 < p < 1$), and exponential ($p = 1$) growth rates. However, as opposed to (2.3), the doubling time for (4.9) admits the closed form expression for any $0\le p\le 1$. Indeed, according to (2.2), in the case of sub-exponential growth, i.e., when $0 < p < 1$, for $F(s)$ defined in (4.9) one concludes

Thus, for an emerging outbreak, in the absence of reliable information about transmission mechanisms, the generalized-growth model could provide a starting point for the assessment of the potential outcomes, particularly when the epidemiological characteristics of the disease are poorly understood or subject to substantial uncertainty. This implies that the sequence of doubling times, $\{ \Delta t_{j}\},$ increases according to

which is precisely the first term in (3.5). As mentioned before, $\lim_{p\to 1}\frac{(2^jC_0)^{1-p}\left[2^{1-p}-1\right]}{r(1-p)} = \frac{\ln 2}{r}$, the doubling time for the simple exponential growth, $C' = rC.$

5.   Doubling times in terms of time variable

In this section, we continue to explore the progression of doubling times for the classical Richards model

The special S-shaped logistic curve, given by (5.1), is described by the following equation:

Substituting this expression for $C(t)$ into (4.2), one arrives at

This identity reinforces the fact that cumulative cases would double unless $C_0\ge \frac{K}{2}$ (and this is unlikely for most outbreaks). The cumulative cases will continue to double until the moment

that is, until the moment $C(t)$ reaches the value $K/2$. At that same moment, the argument of the logarithm in (5.3) becomes negative. As $t$ approaches $\hat t = \frac{1}{ar}\ln \mu$ (and $C(t)$ approaches $K/2$), the term $\frac{1}{ar}\ln \left\{1-\frac{(2^a-1)C_0^ae^{art}}{K^a-C_0^a}\right\}$ in (5.3) goes to negative infinity, which means the doubling time approaches infinity.

For comparison, the doubling time for the generalized-growth model, $C' = rC^p,$ (as a function of time passed from the start of the outbreak) takes the form

Obviously, $\lim_{p\to 1}\left(2^{1-p}-1\right)\left\{t+\frac{C_0^{1-p}}{(1-p)r}\right\} = \frac{\ln 2}{r}$, which is the doubling time for $C' = rC.$

6.   Numerical simulations

We now investigate how the doubling time depends on various values of the parameters $r$, $K$, and $a$ for classical Richards model (4.1). In the first numerical example, we set $C_0 = 2$, $r = 1.2$, and $a = 0.5$. Our goal is to visualize how various values of $K$ alter the doubling times and the sequence of doubling times. From a theoretical standpoint (see (4.2)), the doubling time must be decreasing with respect to $K$. The numerical experiment confirms this theoretical finding as shown in Figures 1 and 4. As expected, the number of times cumulative cases double increases as the value of $K$ goes up.

In the second numerical example, we set $C_0 = 2$, $r = 1.2$, and $K = 100$. Here the goal is to observe how various values of $a$ alter the doubling times and the sequence of doubling times. Again, from a theoretical standpoint (see (4.2)), the doubling time decreases with respect to $a$. Thus, as $a$ increases, the doubling of cases occurs more quickly and it takes less time for the cumulative cases to approach the carrying capacity as shown in Figures 2 and 5.

In the third example, we fix $C_0 = 2$, $a = 0.5$, and $K = 100$ and we study the dependence of doubling times on the parameter $r$. Numerical simulations illustrate that the dependence is very similar to the one of $a$, that is, the doubling time decreases as $r$ goes up and, as $r$ increases, the cumulative cases double more quickly and it takes less time to reach the threshold (Figures 3 and 6).

7.   Comparison of doubling times

Here our goal is to compare the sequences of doubling times across various models. In the first experiment, we examine how the sequences of doubling times compare for Richards model and the modified Richards model. In the case of both models, we set $a = 0.5$, $r = 1.2$, $C_0 = 1$, and $K = 100$, and our goal is to see how the sequence of doubling times for the modified Richards model behaves as $p$ approaches $1$. To that end, we look at the sequence of doubling times for $p = 0.5, 0.7, 0.8, 0.9$, and $0.95$. In line with our theoretical study, as $p$ approaches $1$, the sequence of doubling times for the modified Richards model approaches that of the Richards model as illustrated in Figure 7.

In the second experiment, we compare the generalized growth model to the modified Richards model. Clearly, the modified Richards model becomes the generalized growth model as $a \rightarrow \infty$. In order to simulate this theoretical result, for both models we set $K = 100$, $r = 1.2$, $p = 0.5$, and $C_0 = 1$, and the modified Richards model is considered for $a = 0.5, 0.9, 1, 1.5,$ and $2.5$. From Figure 8, one can see that as $a$ increases, the sequences of doubling times for the modified Richards model approaches that of the generalized growth model.

In the third numerical example, we compare the sequence of doubling times for the generalized growth model and the simple exponential model. For this numerical example, we set $r = 1.2$ and $C_0 = 1$, and then the generalized growth model is studied for $p = 0.5, 0.7, 0.8, 0.9$, and $0.95$. Figure 9 confirms that as $p$ approaches $1$, the sequence of doubling times of the generalized growth model approaches that of the exponential model.

8.   Synthetic data analysis

Our goal here is to estimate unknown parameters $r$, $a$, $p$, and $K$ for the modified Richards model using least squares minimization, specifically through the lsqcurvefit algorithm in MATLAB. The parameters are extracted from early wave of the outbreak. Then the doubling times are calculated from these extracted parameters.

We start by using synthetic data in order to assess if three methods for uncertainty quantification, Normal Approximation, Empirical Likelihood, and Jackknife Empirical Likelihood (see Appendix for details), employed to calculate $95\%$ confidence intervals (CIs) for the mean doubling time, will give rise to the CIs that contain true doubling times' values. Note that the Normal Approximation (NA) method provides a formulaic representation for calculating confidence intervals, and the calculation of confidence intervals is fast. On the flip side, NA method relies on the use of the normal family of distributions which can be applied when the sample size is large (typically when $n > 30$). In the case of limited emerging data, the method likely won't be applicable for accurate uncertainty quantification since the normal distribution may not be appropriate especially if the emerging incidence data is heavily skewed. Empirical Likelihood (EL) method doesn't require the use of this assumption and is more robust against non-normality and skewed data sets. Inferences are done using only the empirical data. The EL method can be easily implemented in most statistical software. However, if the sample size is small, the EL algorithm may diverge or fail to produce accurate confidence intervals depending on the variability of the data. And as the sample size increases, the computational burden increases over the NA method.

The Jackknife Empirical Likelihood (JEL) method makes no use of the underlying distribution of data, and inferences are done using only the empirical data. It applies to both normal and skewed distributions. But, unlike the EL method, the JEL performs well with both small and large samples. Unfortunately, the computational cost of JEL is greater than the cost of EL and NA methods. The results take longer to obtain. So, having consistent results with all three algorithms, allows to confirm that the CIs are reliable regardless of whether the incidence data is normal or heavily skewed and whether the sample size is large or limited.

To quantify uncertainty in the reconstructed doubling times' values, we conduct two numerical experiments. In the first experiment, we set

These parameters are used to solve the forward ODE problem and to generate cumulative data on a given time interval $[t_1, t_n]$. Afterwards, random Gaussian noise (with $0$ mean and standard deviation equal to $1$) is added to the generated cumulative curve to get "real" epidemic data. Then MATLAB lsqcurvefit is employed to fit the model and to extract estimated parameters $a$, $p$, $K$, and $r$. Since incidence data is known to be positive in a real-life setting, we use uniform noise in the case where the incidence becomes negative.

We perform uncertainty quantification for the extracted parameters by assuming Poisson error structure. A total of $400$ bootstrap iterations are carried out and the extracted parameters are used to calculate the sequence of doubling times for all $400$ bootstrap iterations. Displayed in Figure 10 are the graphs containing the reconstructed cumulative incidence curves with noisy data and the extracted parameters $a$, $p$, $K$, and $r$ (see also Table 1). Table 4 (see Appendix) shows the true moments when the doubling of time occurs given the original values of $a$, $p$, $K$, and $r$, and the 95$\%$ confidence intervals obtained by Normal Approximation, Empirical Likelihood, and Jackknife Empirical Likelihood methods (see Appendix for details). In Table 5 (see Appendix), the true intervals of time between consecutive doublings of cases and the corresponding CIs (estimated by using the Normal Approximation, Empirical Likelihood, and Jackknife Empirical Likelihood methods) are presented. Figure 11 shows the comparison of true time when the doubling occurs with reconstructed time (left) and the comparison of true intervals of time between doubling of cases with reconstructed intervals of time (right).

One can clearly see that in the first numerical simulation, the $95\%$ confidence intervals, obtained from the parameter estimation inverse problem, correctly capture the true parameters of the modified Richards model (see Table 1) as well as the true doubling times (Tables 4 and 5). Looking at the constructed intervals, the Normal Approximation method tends to produce the widest intervals.

The Empirical Likelihood and Jackknife Empirical Likelihood methods produce narrower intervals. Among the last two, the Jackknife Empirical likelihood produces narrower intervals which is to be expected. In the second experiment, we set the following original values

The simulations are repeated using the three uncertainty quantification methods described above, and the results are displayed in Figures 12 and 13 and in Tables 2, 6, and 7 (see Appendix for table results).

In the second numerical simulation, the $95\%$ confidence intervals accurately capture the true parameters and the true doubling times for the modified Richards model. The Empirical Likelihood and Jackknife Empirical Likelihood methods both produced narrower intervals over the NA method but from a practical standpoint, they're both roughly the same up to a certain decimal point as was seen in our simulation studies. Since narrower intervals tend to be more precise, the JEL method overall outperformed the NA and EL methods.

9.   Real data analysis

In the following examples, we look to estimate the doubling times using real data for the states of New Hampshire and Hawaii.

Assume that a wave of an outbreak of COVID-19 originates on day $t_1$ and ends on day $b$ (which is unknown at the early stage of disease transmission). Let cumulative data, $d_{\delta}$, be reported at $t = t_1 < t_2 < ... < t_n$ with $t_n$ being smaller than $b$. In that case, given limited data, $d_{\delta}$, at the start of a new wave one has to solve the following constrained minimization problem

To solve this minimization problem in a stable manner, we use the MATLAB built-in function lsqcurvefit. This function solves the nonlinear curve fitting problem using Levenberg - Marquardt optimization algorithm, which allows to extract the system parameters of $r$, $a$, $K$, and $p$ that best model the early wave of the outbreak as illustrated in Table 7. After $p$, $K$, $a$, and $r$ have been recovered from the early epidemic data, $400$ additional cumulative bootstrap curves are generated by adding Poisson error structure to the series of reported cases in order to quantify uncertainty in the reconstructed parameters.

In Figures 15 and 18, we show the reconstructed cumulative incidence for the two states with real cumulative incidence data for comparison (along with the histograms of the obtained values of those parameters along with their $95\%$ confidence intervals). We also calculate the doubling times for the extracted values of the parameters $r$, $a$, $K$, and $p$. Figures 16 and 19 include the moments of time the doubling occurs and the intervals of time between consecutive doubling of cases. Expression (3.5) and the extracted parameter values are used to calculate the sequence of doubling times until cumulative cases stop doubling.

In the first example, we show the results obtained for the state of New Hampshire between March 2 to August 25, 2020 (See Figures 14–16), and in the second example, we present the results obtained for Hawaii from March 5 to May 28, 2020 (Figures 17–19) ^[12].

To quantify uncertainty in the extracted doubling times, we refit the model to $M = 400$ additional data sets for cumulative cases (assuming Poisson error structure). This yields $400$ sequences of doubling times for each state. Tables 8–11 display the Normal Approximation, Empirical Likelihood, and Jackknife Empirical Likelihood methods used to construct the $95\%$ confidence intervals for the estimated doubling times (see Appendix).

The constructed 95$\%$ confidence intervals for the sequences of doubling times for the states of New Hampshire and Hawaii are displayed in Tables 8–11 (see Appendix). These tables show the resulting confidence intervals for the doubling times of cumulative cases during the ascending phase of the first wave of the outbreak as well as the times between consecutive doubling of cases during the same phase. By looking at the CIs, we can see that the intervals for doubling times are very narrow which is not surprising because the intervals constructed for our recovered parameters are also quite narrow. Thus, the recovered doubling times will be less variable as the iterative process is carried out. In addition, the NA confidence intervals are generally wider than the non-parametric EL and JEL methods which is consistent with the simulation studies conducted above for synthetic data. The non-parametric tests make no assumptions about the distribution of the doubling times at each time point so they are constructed from the data only which is important when we cannot make any claims about the distribution. Therefore they are considered more robust due to their reliance on fewer assumptions. As such, the JEL confidence intervals may be the most applicable to our purposes since they do not violate any major assumptions about the population distribution and they are generally the narrowest among all the intervals constructed.

10.   Conclusions

In this paper, we derived theoretical formulas for calculating the doubling times for the modified Richards model and important special cases. Using these theoretical calculations, we analyzed the growth of epidemics in terms of the number of times cumulative cases double until the epidemic peaks and the rate at which intervals between consecutive doubling times change. We explored the relationships between the doubling times of the modified Richards model, the exponential growth model, the generalized growth model, and the classical Richards model.

The derived formulas have been investigated numerically using, first, synthetic data sets with predetermined values of the parameters governing the modified Richards model to assess model accuracy and validate our use of various uncertainty quantification methods. We employed three methods to quantify uncertainty: Normal Approximation, Empirical Likelihood, and Jackknife Empirical Likelihood. Normal Approximation tended to always produce the widest intervals due to the assumed structure of the probability distribution underlying the doubling times at each moment cumulative cases double. Empirical Likelihood and Jackknife Empirical Likelihood produced narrower intervals and are generally considered more precise due to this property. The Normal approximation method can be calculated quickly using common statistical functions built into R. Empirical Likelihood and Jackknife Empirical Likelihood computationally take a bit longer to construct as there are no built in functions to calculate them quickly. Nevertheless, due to their non-parametric nature, whenever the calculations are entirely based on the data, we recommend using Jackknife Empirical Likelihood since the constructed intervals are the most precise, and in practice have better small sample performance.

In the future, we would like to try other forms of Empirical Likelihood methods to assess their accuracy and to see how they perform in determining the true doubling time for each time point. Another important topic that needs to be studied next is the extension of the modified Richards model to the case of elaborate epidemic trajectories aggregating multiple asynchronous sub-epidemics, since a large number of COVID-19 incidence data curves do not have a simple bell-shape behavior for each phase of the outbreak ^[13,14].

Acknowledgments

This study is supported by NSF awards 1818886, 2011622, 1610429, and 1633381, by NIH grant R01 GM 130900, and by Simons Foundation grant 638679.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix

A.   Generalized logistic growth model

All epidemic models considered above are nested in the generalized logistic model ^[5] defined as

where $r > 0$ represents the intrinsic growth rate in the absence of any limitation to disease spread, $0\le p\le 1$, is a "deceleration of growth" parameter, $K > 0$ is the size of the epidemic, $a > 0$ is a parameter that measures the extent of deviation from the S-shaped dynamics of the classical logistic growth model, and $\gamma > 0$ controls the flexibility of the inflection point. Just like the modified Richards model, the generalized logistic model (A.1) enables us to capture different growth profiles ranging from constant incidence ($p = 0$), polynomial growth ($0 < p < 1$), to exponential growth ($p = 1$) at the early accenting stage of the outbreak. While it is unclear if the addition of extra parameter, $\gamma$, gives any considerable advantage from biological standpoint, mathematically the study of doubling times can easily be extended to the generalized logistic model (A.1). Indeed, using the same substitute, $u: = (\frac{s}{K})^{a}$, for the doubling time, $\Delta t$, one obtains

where $\delta: = 1 - \frac{1-p}{a}$ and $b: = \frac{K^{1-p}}{ar}$. Given that

one arrives at the following expression for the doubling times

In the above, $\Gamma(\gamma+1): = \int_0^\infty \mbox{e}^{-t}\, t^\gamma\, dt$ is the gamma function. From (A.3) and the definitions of $b$ and $\delta$, one concludes that

This brings about the doubling times sequence in the form:

B.   Uncertainty quantification

In practice, in order to estimate the doubling times of an outbreak, one has to extract the unknown parameters by fitting model predictions to the epidemic data. Using these extracted parameters, one employs (2.6) to calculate the sequence of doubling times until cumulative cases cannot double anymore.

In order to quantify uncertainty in the extracted doubling times, we refit the model to $M = 400$ additional data sets for cumulative cases (assuming Poisson error structure), which results in $M$ approximate sequences of doubling times. Our goal is to construct 95$\%$ confidence intervals for the mean doubling time at every moment where cases double, i.e., such confidence intervals where the true doubling time falls within the interval 95$\%$ of the time. Define the true doubling time at each point where cases double as $\mu_{double}$. We calculate and compare three different confidence intervals for estimating $\mu_{double}$: Normal Approximation, Empirical Likelihood, and Jackknife Empirical Likelihood.

For the Normal Approximation method ^[15,16], which is based on the Central Limit Theorem, with $M = 400$, we have that the distribution of the mean doubling time at each point is approximately normal. Let $X_1, ..., X_M$ be an i.i.d random sample of doubling times at each moment cumulative cases double. Let $\bar{x}_{double}$ be the sample mean of $X_1, ..., X_M$ at each moment cumulative cases double. Therefore, we calculate 95$\%$ confidence intervals for the mean doubling time, $\mu_{double}$, by adding and subtracting $1.96$ standard errors from the mean. In other words, our 95$\%$ confidence intervals for the mean doubling time is calculated by:

Empirical Likelihood is a non-parametric method for constructing confidence intervals without having to assume the form of the underlying distribution introduced by Owen ^[17,18]. Empirical likelihood enables us to successfully incorporate the advantages of the likelihood methods. Unnecessary assumption of a family of distributions can be avoided in empirical likelihood inference in the sense that the object is driven by data.

Empirical Likelihood has become very popular in recent years. Huang and Zhao ^[19] used empirical likelihood inference on a bivariate survival function under univariate censoring. Cheng et al. ^[20] introduced empirical likelihood inference for semi-parametric additive isotonic regression. In addition, Zhang et. al. ^[21] used empirical likelihood for testing the symmetry of statistical distributions. Let $X_1, ..., X_M$ be an i.i.d. random sample of doubling times from an unknown distribution function $F(x) = P(X \leq x)$. For notation, define $F(x-) = P(X < x)$. Thus, $P(X = x) = F(x)-F(x-)$. For simplicity, we denote $P(X = X_i) = p_i$, $i = 1, ..., M$, where $p_i \geq 0$ and $\sum_{i = 1}^{M} p_i = 1$. The empirical cumulative distribution function (ECDF) of $X_1, ..., X_M$ is

The notation $1_{X_i \leq x}$ outputs a value of 1 if $X_i \leq x$ is true, and 0 otherwise. In addition, given $X_1, ..., X_M$, which are assumed independent with common CDF $F_0$, the non-parametric likelihood of the CDF $F$ is

The empirical maximum likelihood estimator (EMLE) based on the empirical likelihood, puts equal probability mass $1/M$ on the $M$ observed values $X_1, X_2, ..., X_M$. This allows us to calculate the empirical likelihood ratio. Using Lagrange multipliers, we can show that the $p_i$'s which maximize $L(F) = \prod_{i = 1}^{M} p_i$ subject to $p_i \geq 0$, $\sum_{i = 1}^{M} p_iX_i = \mu_{double}$, and $\sum_{i = 1}^{M} p_i = 1$ are given by:

The empirical likelihood ratio is therefore:

where $\lambda$ is the unique solution of

Let $E(X_1^2) < \infty$. Then, the empirical loglikelihood ratio for the true mean doubling time is

Lastly, let $\mu_{0}$ be the true doubling time at a certain moment in time. We have as $n \rightarrow \infty$

According to Owen ^[18], the $100(1-\alpha)\%$ confidence interval for $\mu_{double}$ is given by

where $\chi^2_{1, 1-\alpha}$ is the $1-\alpha$ quantile of the $\chi^2_1$ distribution.

There are some computational challenges for the application of empirical likelihood methods in practice. To avoid these, Jing et al. ^[22] introduced the "jackknife", as a kind of re-sample method, and combined it with empirical likelihood and called it the Jackknife Empirical Likelihood (JEL). This method can greatly simplify the optimization procedures of traditional empirical likelihood methods. Various research articles have been published using jackknife empirical likelihood. Zhao et al. ^[23] used JEL inference for the mean absolute deviation. Sang et al. ^[24] proposed JEL method for estimating Gini correlations. Lin et al. ^[25] published the method for the error variance using JEL in linear regression models. In addition, the JEL method can also be applied to Bayesian inference, which was published by Cheng and Zhao ^[26].

Let $\hat{x}_{i, double}$ denote the mean doubling time calculated with the $i^{th}$ observation deleted where $i = 1, .., M$. We define $\bar{x}_{double}$ as the mean doubling time of the entire i.i.d random sample $X_1, ..., X_M$. Suppose that $\hat{V}_i$ denotes the jackknife pseudo-value, which is obtained by

Thus, the jackknife empirical likelihood evaluated at $\mu_{double}$ can be defined as

where $\lambda$ is the solution of the following equation

With $\lambda$, we can now calculate the $-2\log$ of the empirical likelihood ratio as

Let $E(X_1^2) < \infty$. Then, by Jing et al. by ^[22] the empirical loglikelihood ratio for the true mean doubling time is

Lastly, let $\mu_{0}$ be the true doubling time at a certain moment in time. We have as $n \rightarrow \infty$

Thus, the Jackknife Empirical Likelihood confidence interval for $\mu_{double}$ is obtained by

where $\chi^2_{1, 1-\alpha}$ is the $1-\alpha$ quantile of $\chi^2_{1}$.