Estimating parameters for the COVID-19 epidemic

Summary

COVID19 is an emerging pathogen and much of its epidemiology remains poorly understood. Identifying the key processes that shape transmission and estimating the relevant model parameters is therefore an important task. This document presents arguments and analysis to support the estimation of a number of key quantities:

• Incubation period (\(1/\sigma\))
• Recovery/isolation rate (\(1/\gamma\))
• Interval from hospitalization to notification (\(\tau_1\)) and notification rate (\(\eta\))
• Case fatality rate (\(m\))
• Basic reproduction number (\(R_0\))
• Case detection rate (\(q\))
• Transmissibility (\(\beta\))
• Additional parameters

To estimate these parameters, we use the following resources:

1. A record of case reports at the province level in China maintained by the CEID at the University of Georgia
2. A global linelist maintained by Be Outbreak Prepared (BOP) This dataset was preceded by the “Moritz Kraemer” linelist and has since been published in the Lancet

The Be Outbreak Prepared Linelist contains individual-level information on COVID19 cases. The BOP dataset is not comprehensive, there are currently ~1.2 million cases in it, but many cases have limited information and many cases are likely not represented in the dataset. Russia, Italy, Germay, the US, France, Belgium, and China have the most cases in the linelist. New linelist cases in the US cut off at the beginning of April. The number of new cases over time in 5 countries is shown below.

Findings are preliminary and subject to change, pending changes in the underlying data. Results have not been peer-reviewed, but have been prepared to a professional standard with the intention of providing useful information about a rapidly developing event.

Incubation period

Work on distribution of the incubation period (time from exposure to symptom onset) has estimated it to vary from around 4-12 days (MIDAS network parameter list). Some peer-reviewed estimates of the incubation period include:

• Gamma (5.8, 0.94); (Lauer et al. 2020 in Annals of Internal Medicine). They also fit Weibull and Erlang distributions.
• Mean 5.2 (95% CI 0-14 days); (Zhang et al. 2020 in Lancet Infectious Diseases).
• Median 4 days (95% CI 3.5-5.1 days); (Sanche et al. 2020 in EID)

We use the Be Prepared linelist to estimate incubation period. From this dataset, we assume incubation period was the delay from day of travel from Wuhan and the symptom onset. Only cases outside Wuhan were used in this calculation since we used travel to Wuhan as the time of exposure. We fit gamma and Erlang distributions to these intervals. The Erlang distribution is a special case of a gamma distribution with an integer shape parameter. The PDF of the Erlang distribution is

\[\begin{equation} f(x; k, \lambda) = \frac{\lambda^{k}x^{k-1}e^{-\lambda x}}{(k-1)!} \end{equation}\]

for rate parameter \(\lambda \geq 0\). We estimate the Erlang distribution of isolation intervals numerically and analytically by optimizing the negative log likelihood. The MLE of this distribution can be solved by finding the derivative of the log likelihood function:

\[\begin{equation} L = n k log(\lambda) + (k-1) \sum_{i=1}^{n}{x_i} - \lambda \sum_{i=1}^{n}{x_i} - n log(k-1)! \end{equation}\]

which has solution \(\hat{\lambda}=\frac{k}{\sum_{i=1}^{n}{x_i}}\).

The below histograms show incubation period from the linelist outside of Wuhan, with periods longer than 0 and shorter than 50 days. We filtered the linelist to include only cases that had known exposures in: Wuhan, high risk country, contact with a confirmed case. Hubei, Hubei province, Italy, and USA.

Incubation period for US cases in the Be Outbreak Prepared linelist.

Incubation period for US cases in the Be Outbreak Prepared linelist.

The above plot shows the negative log likelihood curve of Erlang shape parameter (k) given optimal rate parameter for the incubation period based on hospitalization date (right). Blue line shows the MLE of the rate parameter (lambda) and the dashed grey line shows the estimates within 2 log units of the MLE.

In total the incubation period is based on records from 547 case histories. The incubation period estimated here is consistent with published data.

Isolation rate

Isolation rate is defined as the reciprocal of the interval between onset of symptoms and hospital admission.

Work on distribution of the isolation rate (1/time from onset to isolation) has estimated it to vary from around 2.7-4.6 days (MIDAS network parameter list). Isolation has also been defined as confirmation or hospitalization. The relevance of which measure to use varies depending on the population of interest. Whether you use hospitalization, isolation, or confirmation could contribute to confusion when comparing estimates of this parameter.

• Between December 29 and January 23, Liu et al. 2020 (pre-print on Biorxiv) reported cases in China had an average period from onset of symptoms to isolation of 2.9±3.0 days, and it decreased from 6.7 days in cases before January 9, to 0.7 days in cases after January 19.
• Between January 14 and February 12, Bi et al. 2020 (pre-print on MedRxiv) reported cases in Shenzen detected through symptom-based surveillance were confirmed on average 5.5 days (95% CI 5.0, 5.9) after symptom onset (Figure 3, Table S2); compared to 3.2 days (95% CI 2.6,3.7) in those detected by contact-based surveillance.

We investigate the isolation interval using the BOP dataset as the time from symptom onset to hospitalization and confirmation. Only cases that had positive number of days from symptom onset to hospitalization or confirmation were included in this analysis (a significant number of cases showed they were hospitalized prior to having symptom onset - many in Germany).

In total the isolation delay is based on records from 1208 case histories.

Negative log likelihood curve of Erlang shape parameter (k) given optimal rate parameter for the isolation interval based on hospitalization date (right). Blue line shows the MLE of the rate parameter (lambda) and the dashed grey line shows the estimates within 2 log units of the MLE.

Next plots show isolation interval (onset to hospitalization) over time with a loess regression line to detect trends. In China, isolation interval decreased from around 2 weeks in early January to a few days in February. Data is limited after that point. Outside of China, it looks like isolation interval was flat around one week January to March and then decreases linearly after that.

To obtain an estimate of the rate of change over time, we use a stepwise linear regression. For China breakpoints, we used two values. On Jan 19, China increased testing in Wuhan suggesting a natural breakpoint when the interval to admission starts to decline. After plotting the data, we noticed another potential breakpoint that seemed to show a starker contrast between interval times: Jan 15th. We have shown results for both breakpoints.

## 
## Call:
## lm(formula = rate ~ days, data = filter(d, country == "China", 
##     date_onset_symptoms <= as.Date("2020-01-15")))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12412 -0.06437 -0.02914  0.01508  0.81868 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.026541   0.084171   0.315    0.753
## days        0.003599   0.002168   1.660    0.100
## 
## Residual standard error: 0.1371 on 94 degrees of freedom
## Multiple R-squared:  0.02848,    Adjusted R-squared:  0.01815 
## F-statistic: 2.756 on 1 and 94 DF,  p-value: 0.1002

## 
## Call:
## lm(formula = rate ~ days, data = filter(d, country == "China", 
##     date_onset_symptoms > as.Date("2020-01-15")))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4407 -0.2226 -0.1096  0.1054  0.6376 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.115537   0.118876   0.972   0.3315  
## days        0.005367   0.002145   2.502   0.0126 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3097 on 535 degrees of freedom
## Multiple R-squared:  0.01157,    Adjusted R-squared:  0.009719 
## F-statistic: 6.261 on 1 and 535 DF,  p-value: 0.01264

In China, isolation rate was not significantly related to time before Jan 15, but significantly increased AFTER to Jan 15.

Outside of China, we saw a natural breakpoint around Feb 15, when the rate began increasing.

## 
## Call:
## lm(formula = rate ~ days, data = filter(d, country != "China", 
##     date_onset_symptoms < as.Date("2020-02-15")))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28771 -0.20492 -0.11923  0.01578  0.68336 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.3093538  0.1259327   2.457   0.0147 *
## days        0.0001519  0.0021064   0.072   0.9426  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3017 on 258 degrees of freedom
## Multiple R-squared:  2.015e-05,  Adjusted R-squared:  -0.003856 
## F-statistic: 0.005198 on 1 and 258 DF,  p-value: 0.9426

## 
## Call:
## lm(formula = rate ~ days, data = filter(d, country != "China", 
##     date_onset_symptoms >= as.Date("2020-02-15")))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4639 -0.2366 -0.1134  0.1899  0.7159 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.382982   0.131594  -2.910  0.00386 ** 
## days         0.008664   0.001362   6.361 6.75e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.317 on 328 degrees of freedom
## Multiple R-squared:  0.1098, Adjusted R-squared:  0.1071 
## F-statistic: 40.46 on 1 and 328 DF,  p-value: 6.748e-10

Similarly outside of China, isolation rate was not significantly related to time before Feb 15, but significantly increased AFTER to Feb 15.

For models, the isolation rate IN China is:

\[\begin{equation} \gamma(t) = \begin{cases} 0.165 & \text{if } t \leq Jan 15\\ 0.165 + .005 t, & \text{otherwise} \end{cases} \end{equation}\]

For models, the isolation rate OUTSIDE China is:

\[\begin{equation} \gamma(t) = \begin{cases} 0.318 & \text{if } t \leq Feb 15\\ 0.318 + .0008 t, & \text{otherwise} \end{cases} \end{equation}\]

Notification rate

Notification rate is the reciprocal of the average interval between hospital admission and case notification. Here we investigate the interval between hospital admission and case notification (at which time the case appears in the data set.) We denote the interval between admission and notification by \(\tau_1\) and define the case notification rate as \(\eta=1/\tau_1\), which is presumed to be a function of time. Again, we estimated this parameter using the the BOP dataset.

## 
## Call:
## lm(formula = rate ~ days, data = filter(d, country == "China"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50385 -0.21436 -0.07536  0.08018  0.59131 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.148030   0.073425  -2.016    0.044 *  
## days         0.011134   0.001283   8.679   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2926 on 1059 degrees of freedom
## Multiple R-squared:  0.0664, Adjusted R-squared:  0.06552 
## F-statistic: 75.32 on 1 and 1059 DF,  p-value: < 2.2e-16

## 
## Call:
## lm(formula = rate ~ days, data = filter(d, country != "China"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67690 -0.27170 -0.08921  0.36965  0.53423 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.8582325  0.0477031  17.991  < 2e-16 ***
## days        -0.0031650  0.0005245  -6.034 2.44e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3333 on 808 degrees of freedom
## Multiple R-squared:  0.04311,    Adjusted R-squared:  0.04193 
## F-statistic: 36.41 on 1 and 808 DF,  p-value: 2.438e-09

In China, the notification rate increased significantly in early January to mid February. Outside China, the notification rate (significantly) increased from February on.

In China, we used the linear regression to obtain an estimate of the rate of change over time. The inverse of this is a hyperbolic function for the increase in case notification rate, which diverges around Day 73.

For notification rate, we assume that the maximum case notification rate is one (i.e. one day between admission and notification) since this is the reporting interval of almost all case reports. Thus, our model of the notification rate in China is

\[\begin{equation} \eta(t) = \begin{cases} (-0.1335t+10.8)^{-1}, & \text{if } t <73\\ 1, & \text{otherwise} \end{cases} \end{equation}\]

Time to death

We look at the interval between hospitalization or confirmation and date of death in BOP data. Other estimates of this interval:

• 11.2 days (95% CI 8.7-14.9 days); (Linton et al. 2020 in medrxiv)
• 8.3 days (95% CI 6.4-10.5 days); (Sanche et al. 2020 in EID) Estimated median with Weibull distribution. This paper also fit lognormal and gamma distributions to hospitalization to death interval.

There is not that much information in the BOP dataset that we can use to calculate this interval. We remove negative values and intervals above 50 days as these seemed unreasonable.

Hospitalization or confirmation to death.

## 
##        Canada         China      Ethiopia        France        Gambia 
##             2            44             1             2             1 
##         Italy         Japan   Philippines    San Marino United States 
##             3             4            43             1            11 
##      Zimbabwe 
##             1

Hospitalization or confirmation to death.

Unfortunately, there is not that much information about hospitalization to death in this dataset (56 rows of data). By using the minimum date between date of hospitalization and date of confirmation, we have 115 rows of data. Much of this data is from China and the Philipines.

Estimating \(R_0\)

## 
## Call:
## lm(formula = log(new_cases + 1) ~ as.numeric(day), data = filter(d, 
##     day > 40, day < 65))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7943 -0.7053  0.1976  0.5649  2.2970 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -17.26523    1.44704  -11.93 4.44e-11 ***
## as.numeric(day)   0.41434    0.02733   15.16 3.95e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9267 on 22 degrees of freedom
## Multiple R-squared:  0.9127, Adjusted R-squared:  0.9087 
## F-statistic: 229.9 on 1 and 22 DF,  p-value: 3.953e-13

## [1] "R0=2.28"

Case fatality rate

Here we estimate case fatality rate using the number of reported deaths (\(D\)) and recoveries (\(R\)) in each province.

\[\begin{equation} CFR = \frac{D}{D+R} \end{equation}\]

For the purposes of this analysis, a “case” is a clinically relevant, detected infection. To calculate the infection fatality rate, these estimates would need to be adjusted for undetected infections (which may not have the same fatality rate).

The following is a plot of CFR by province or region (with Tibet removed because we had limited data - 0 deaths and 1 recovery).

Hubei has the highest estimated CFR with 6.6% (CI: [6.4, 6.8]).

The standard of care for COVID has evolved over the course of the epidemic and the case fatality rate outside Hubei is considerably lower. We estimate the case fatality rate outside of Hubei to be .9% (CI: [.7%, 1%]). The grand mean is 5.5% (confidence interval: [0.054, 0.057]).

Case detection rate

Here we present an argument for how to estimate the case detection rate of SARS-CoV-2 infections in Wuhan during January 2020. The idea is to use outcome data on a cohort of patients from early in the epidemic to estimate the size of the population from which they were sampled by comparing with estimates of case fatality rate from later in the epidemic.

The logic is as follows.

As time \(s\) there are \(I_s\) infections in the population, of which \(D_s\) will result in death. The case fatality rate is therefore defined as

\[\begin{equation} \label{eq:cfr} CFR = \frac{D_s}{I_s}. \end{equation}\]

It is assumed that all deaths are known cases, but that there are an unknown number of infections that have gone undetected. The detection rate is \(q\) such that

\[\begin{equation} \label{eq:cases} C_s = qI_s \rightarrow I_s = C_s/q. \end{equation}\]

A fraction \(p\) of the population were are enrolled in cohort \(m\) at time \(s\). Since the outcome (death or recovery) isn’t known at the time of enrollment, we assume that patients with death and recovery outcomes are included in the cohort in proportion to their fequency in the population, i.e., \(D_m = pD_s\) is the number of patients in the cohort whose infections will result in death and \(C_m = pC_s\) is the number of patients in the cohort whose infections will result in recovery. Rearranging these equations, we have \(D_s = D_m/p\) and \(C_s=C_m/p\). Just making substitutions, we have

\[\begin{equation} CFR = \frac{D_s}{I_s} = q\frac{D_s}{C_s} = q \frac{pD_m}{pC_m} = q\frac{D_m}{C_m} = q (CFR_m). \end{equation}\]

where the ratio \(CFR_m = \frac{D_m}{C_m}\) is just the case fatality rate in cohort \(m\). Rearranging, we have a formula for \(q\).

\[\begin{equation} q = \frac{CFR}{CFR_m}. \end{equation}\]

Ghani et al have shown that, during the initial phase of an outbreak when the outcomes of all cases are not yet known, the simple case fatality estimator

\[\begin{equation} \label{eq:e1} e_1=D_m/C_m \end{equation}\]

typically provides an underestimate due to right censoring. A better estimator is provided by the formula

\[\begin{equation} \label{eq:e2} e_2(s) = \frac{D_m}{D_m+R_m}, \end{equation}\]

where \(D_m\) is understood to be the number of patients in cohort \(m\) that have died at the time of estimation and \(R_m\) is the number that have recovered, while the outcomes of the remaining \(X_m = C_m-D_m-R_m\) are unknown.

Chen et al. report on \(C_m=99\) patients treated in Wuhan Jinyintan Hospital from Jan 1 to Jan 20, 2020. All 2019-nCov patients were eligible for treatment in the hospital and all treated patients were enrolled in the study. The majority (possibly all) patients would have been infected prior to widespread testing which began on January 19 and before high case fatality rate was widely known. Thus, it is reasonable to assume that, at least for the purpose of statistical analysis, this cohort represents a random sample of the case detection process during this time.

Chen et al. report that, at the time of their writing, \(D_m=11\) of 99 patients died, \(R_m=31\) patients had been discharged, and \(C_m-D_m-R_m=57\) patients remained in the hospital.

Using equation , the cohort case fatality rate on January 25 (the final day in the study) may be estimated to be \(11/(11+31)=0.262\). Several news sources (https://www.nytimes.com/interactive/2020/world/asia/china-coronavirus-contain.html#virulence, https://www3.nhk.or.jp/news/html/20200124/k10012257631000.html) report that the actual case fatality rate may be around 3% (attributing this information to the WHO). It is unclear if this estimate is based estimator or an analysis that either excludes or accounts for right censoring. If the estimate of 3% is close to the truth, then we have a case detection rate for January 1 - January 20 of \(q=0.03/0.262 \approx 0.11\).

We can estimate the uncertainty in this number by randomly sampling from among the \(11+31=42\) cases for which outcome is known and computing a distribution.

##       2.5%      97.5% 
## 0.07411765 0.21000000

Transmissibility (\(\beta\))

Transmissibility (\(\beta\)) is defined in the usual way by rearranging the equation \(R_0 = \frac{\beta}{\gamma}\) to give \(\beta = R_0 \times \gamma\). Taking our estimate from the early stages of the outbreak and fixed infectious period (\(\gamma = 1/7\)) we conclude that \(\beta\) ranges from \(\beta = R_0 \times \gamma = 1.4/7 = 0.2\) to \(\beta = R_0 \times \gamma = 3.5/7 = 0.5\). For comparison, if we use the Imperial College estimate of \(R_0=2.6\), we obtain \(\beta = 2.6/ 7 = 0.371\).

Additional parameters

December 1, 2019 is taken to be the notional time of first infection for the purposes of this analysis.

Population size of Wuhan: \(N \approx 11,081,000\) (Wikipedia)

Population size of Hubei: \(N \approx 59,002,000\) (Wikipedia)