Using Early Data to Estimate the Actual Infection Fatality Ratio from COVID-19 in France

Abstract

The number of screening tests carried out in France and the methodology used to target the patients tested do not allow for a direct computation of the actual number of cases and the infection fatality ratio (IFR). The main objective of this work is to estimate the actual number of people infected with COVID-19 and to deduce the IFR during the observation window in France. We develop a ‘mechanistic-statistical’ approach coupling a SIR epidemiological model describing the unobserved epidemiological dynamics, a probabilistic model describing the data acquisition process and a statistical inference method. The actual number of infected cases in France is probably higher than the observations: we find here a factor ×8 (95%-CI: 5–12) which leads to an IFR in France of 0.5% (95%-CI: 0.3–0.8) based on hospital death counting data. Adjusting for the number of deaths in nursing homes, we obtain an IFR of 0.8% (95%-CI: 0.45–1.25). This IFR is consistent with previous findings in China (0.66%) and in the UK (0.9%) and lower than the value previously computed on the Diamond Princess cruse ship data (1.3%).

1. Introduction

The COVID-19 epidemic started in December 2019 in Hubei province, China. Since then, the disease has spread around the world reaching the pandemic stage, according to the WHO [1], on 11 March. The first cases were detected in France on 24 January. The infection fatality ratio (IFR), defined as the number of deaths divided by the number of infected cases, is an important quantity that informs us on the expected number of casualties at the end of an epidemic, when a given proportion of the population has been infected. Although the data on the number of deaths from COVID-19 are probably accurate, the actual number of infected people in the population is not known. Thus, due to the relatively low number of screening tests that have been carried out in France (about five in 10,000 people in France to be compared with 50 in 10,000 in South Korea up to 15 March 2020; sources: Santé Publique France and Korean Center for Disease Control) the direct computation of the IFR is not possible. Based on the PCR-confirmed cases in international residents repatriated from China on January 2020, Verity et al. [2] obtained an estimate of the infection fatality ratio (IFR) of 0.66% in China, and, adjusting for non-uniform attack rates by age, an IFR of 0.9% was obtained in the UK [3]. Using data from the quarantined Diamond Princess cruise ship in Japan and correcting for delays between confirmation and death, Russel et al. [4] obtained an IFR of 1.3%.

Using the early data (up to 17 March) available in France, our objectives are: (1) to compute the IFR in France, (2) to estimate the number of people infected with COVID-19 in France, and (3) to compute a basic reproduction number $R_0$.

2. Materials and Methods

2.1. Data

We obtained the number of positive cases and deaths in France, day by day from Johns Hopkins University Center for Systems Science and Engineering [5]. The data on the number of tests carried out was obtained from Santé Publique France [6]. As some data (positive cases, deaths, number of tests) are not fully reliable (example: 0 new cases detected in France on 12 March 2020), we smoothed the data with a moving average over 5 days. Official data on the number of deaths by COVID-19 in France only take into account hospitalised people. About 728,000 people in France live in nursing homes (EHPAD, source: DREES [7]). Recent data in the Grand Est region (source: Agence Régionale de Santé Grand Est [8]), report a total of 570 deaths in these nursing homes, which have to be added to the official count (1015 deaths on 31 March).

2.2. Mechanistic-Statistical Model

The mechanistic-statistical formalism, which is becoming standard in ecology [9,10,11] allows the analyst to couple a mechanistic model that describes a latent variable, here an ordinary differential equation model (ODE) of the SIR type, and uncertain, non-exhaustive data. To bridge the gap between the mechanistic model and the data, the approach uses a probabilistic model describing the data collection process. A statistical method is then used for the estimation of the parameters of the mechanistic model.

Mechanistic model. The dynamics of the epidemic are described by the following SIR compartmental model:

$$\left\{\begin{array}{c}{\displaystyle S^{\prime }\left(t\right)=-\frac{\alpha }{N}S\left(t\right)I\left(t\right),}\hfill \\ {\displaystyle I^{\prime }\left(t\right)=\frac{\alpha }{N}S\left(t\right)I\left(t\right)-\beta I\left(t\right),}\hfill \\ {\displaystyle R^{\prime }\left(t\right)=\beta I\left(t\right),}\hfill \end{array}\right.$$

(1)

with S the susceptible population, I the infected population, R the recovered population (immune individuals) and $N=S+I+R$ the total population, supposed to be constant. The parameter $\alpha$ is the infection rate (to be estimated) and $1/\beta$ is the mean time until an infected becomes recovered. Based on the results in [12], the median period of viral shedding is 20 days, but the infectiousness tends to decay before the end of this period: the results in [13] show that infectiousness starts from 2.5 days before symptom onset and declines within 7 days of illness onset. Based on these observations we assume here that $1/\beta =10$ days.

The initial conditions are $S\left(t_0\right)=N-1$, $I\left(t_0\right)=1$ and $R\left(t_0\right)=0$, where $N=67·{10}^6$ corresponds to the population size. The SIR model is started at some time $t=t_0$, which will be estimated and should approach the date of introduction of the virus in France (this point is shortly discussed at the end of this paper). The ODE system (1) is solved thanks to a standard numerical algorithm, using Matlab^® ode45 solver.

Next we denote by $D\left(t\right)$ the number of deaths due to the epidemic. Note that the impact of the compartment $D\left(t\right)$ on the dynamics of the SIR system and on the total population is neglected here. The dynamics of $D\left(t\right)$ depends on $I\left(t\right)$ trough the differential equation:

$$D^{\prime }\left(t\right)=\gamma \left(t\right)I\left(t\right),$$

(2)

with $\gamma \left(t\right)$ the mortality rate of the infecteds.

Observation model. We suppose that the number of cases tested positive on day t, denoted by ${\widehat{\delta }}_t$, follow independent binomial laws, conditionally on the number of tests $n_t$ carried out on day t, and on $p_t$ the probability of being tested positive in this sample:

$${\widehat{\delta }}_t\sim Bi(n_t,p_t).$$

(3)

The tested population consists of a fraction of the infected cases and a fraction of the susceptibles: $n_t={\tau }_1\left(t\right)I\left(t\right)+{\tau }_2\left(t\right)S\left(t\right)$. Thus,

$$p_t=\frac{\sigma {\tau }_1\left(t\right)I\left(t\right)}{{\tau }_1\left(t\right)I\left(t\right)+{\tau }_2\left(t\right)S\left(t\right)}=\frac{\sigma I\left(t\right)}{I\left(t\right)+{\kappa }_{t}S\left(t\right)},$$

(4)

with ${\kappa }_t:={\tau }_2\left(t\right)/{\tau }_1\left(t\right)$, the relative probability of undergoing a screening test for an individual of type S vs. an individual of type I (probability of being tested conditionally on being S/probability of being tested conditionally on being I). We assume that the ratio $\kappa$ does not depend on t at the beginning of the epidemic (i.e., over the period that we use to estimate the parameters of the model). The coefficient $\sigma$ corresponds to the sensitivity of the test. In most cases, RT-PCR tests have been used and existing data indicate that the sensitivity of this test using pharyngeal and nasal swabs is about 63–72% [14]. We take here $\sigma =0.7$ (70% sensitivity).

2.3. Statistical Inference

The unknown parameters are $\alpha$, $t_0$ and $\kappa$. The parameter $\gamma \left(t\right)$ is computed indirectly, using the estimated value of $I\left(t\right)$, the data on $D\left(t\right)$ (assumed to be exact) and the relationship (2). The likelihood $\mathcal{L}$ is defined as the probability of the observations (here, the increments $\left\{{\widehat{\delta }}_t\right\}$) conditionally on the parameters. Using the observation model (3), and assuming that the increments ${\widehat{\delta }}_t$ are independent conditionally on the underlying SIR process and that the number of tests $n_t$ is known, we get:

$$\mathcal{L}(\alpha ,t_0,\kappa ):=P\left(\left\{{\widehat{\delta }}_t\right\}\right|\alpha ,t_0,\kappa )=\prod _{t=t_i}^{t_f}\frac{n_t!}{\left({\widehat{\delta }}_t\right)!(n_t-{\widehat{\delta }}_t)!}{p}_t^{{\widehat{\delta }}_t}{(1-p_t)}^{n_t-{\widehat{\delta }}_t},$$

(5)

with $t_i$ the date of the first observation and $t_f$ the date of the last observation. In this expression $\mathcal{L}(\alpha ,t_0,\kappa )$ depends on $\alpha ,t_0,\kappa$ through $p_t$.

The maximum likelihood estimator (MLE, i.e., the parameters that maximise $\mathcal{L}$), is computed using the BFGS constrained minimisation algorithm, applied to $-ln\left(\mathcal{L}\right)$, via the Matlab ^® function fmincon. In order to find a global maximum of $\mathcal{L}$, we apply this method starting from random initial values for $\alpha ,t_0,\kappa$ drawn uniformly in the following intervals: $\alpha \in (0,1)$, $t_0\in (1,31)$, (1–31 January) and $\kappa \in (0,1)$. The minimisation algorithm is applied to 10,000 random initial values of the parameters.

The posterior distribution of the parameters $(\alpha ,t_0,\kappa )$ is computed with a Bayesian method, using uniform prior distributions in the intervals given above (a more informative prior distribution has also been tested, see Appendix A). This posterior distribution corresponds to the distribution of the parameters conditionally on the observations:

$$P(\alpha ,t_0,\kappa |\left\{{\widehat{\delta }}_t\right\})=\frac{\mathcal{L}(\alpha ,t_0,\kappa )\pi (\alpha ,t_0,\kappa )}{C},$$

(6)

where $\pi (\alpha ,t_0,\kappa )$ corresponds to the prior distribution of the parameters (therefore uniform) and C is a normalisation constant independent of the parameters. The numerical computation of the posterior distribution is performed with a Metropolis–Hastings (MCMC) algorithm, using four independent chains, each of which with ${10}^6$ iterations, starting from random values close to the MLE (only the second half of the iterations are used to generate the posterior). The Matlab^® codes are available as Supplementary Materials.

The data ${\widehat{\delta }}_t$ used to compute the MLE and the posterior distribution are those corresponding to the period from 29 February to 17 March.

3. Results

Model fit. To assess model fit, we compared the observations, i.e., the cumulated number of cases ${\Sigma }_t:={\Sigma }_{s=1,\dots ,t}{\widehat{\delta }}_s$, with the expectation of the observation model associated with the MLE (expectation of a binomial). Namely, we compared ${\Sigma }_t$ and ${\Sigma }_{s=1,\dots ,t}{n}_s{p}_s^{\ast }$ with

$$p_s^{\ast }=\frac{\sigma I^{\ast }\left(s\right)}{I^{\ast }\left(s\right)+{\kappa }^{\ast }{S}^{\ast }\left(s\right)},$$

(7)

and $I^{\ast }\left(s\right)$, $S^{\ast }\left(s\right)$ the solutions of the system (1) (at time s) associated with the MLE. The results are presented in Figure 1. We observe a good match with the data.

Infection fatality ratio and actual number of infected cases. Using the posterior distribution of the model parameters (the pairwise distributions are presented in Appendix A, see Figure A1), we computed the daily distribution of the actual number of infected peoples. Using the relation (2) together with the data on $D\left(t\right)={\Sigma }_t$, we deduce the distribution of the parameter $\gamma \left(t\right)$, at each date. The IFR corresponds to the fraction of the infected who die, that is:

$$IF{R}_t:=\gamma \left(t\right)/(\gamma \left(t\right)+\beta ).$$

(8)

We thus obtain, on 17 March an IFR of 0.5% (95%-CI: 0.3–0.8), and the distribution of the IFR is relatively stable over time (see Figure A3 in the Appendix A).

Additionally, the distribution of the cumulated number of infected cases ($I\left(t\right)+R\left(t\right)$) across time is presented in Figure 2. We observe that it is much higher than the total number of observed cases (compare with Figure 1). The average estimated ratio between the actual number of individuals that have been infected and observed cases $(I\left(t\right)+R\left(t\right))/{\Sigma }_t$ is eight (95%-CI: 5–12) over the considered period.

Taking into account the data in the nursing homes. The above computation of the IFR is based on the official counting of deaths by COVID-19 in France, which does not take into account the number of deaths in nursing homes. Based on the local data in Grand Est region, we infer that the IFR that we computed has been underestimated by a factor about $(1015+570)/1015\approx 1.6$, leading to an adjusted IFR of 0.8% (95%-CI: 0.45–1.25).

Basic reproduction number. With SIR systems of the form (1), the basic reproduction number $R_0$ can be computed directly, based on the formula $R_0=\alpha /\beta$ [15]. When $R_0<1$, the epidemic cannot spread in the population. When $R_0>1$, the infected compartment I increases as long as $R_{0}S>N=S+I+R$. We computed the marginal posterior distribution of the basic reproduction number $R_0$. This leads to a mean value of $R_0$ of 3.2 (95%-CI: 3.1–3.3). The full distribution is available in the Appendix A (Figure A4).

Sensitivity of the results with respect to the fixed model parameters. We computed the MLE with a larger infectious period ($1/\beta$) of 20 days estimated by [12]. This leads to a much larger basic reproductive number $R_0=4.8$ and a factor $\times 15$ between the reported cases and the actual number of cases. However, the value of the IFR remains unchanged (0.5%). We also checked if the window width of the smoothing (moving average over 5 days) had an impact on our results. Computations of the MLE with a window width of 3 days (and $\beta =1/10$) led to the same results as those presented above, namely a $R_0=3.2$ and an IFR of 0.5%.

4. Discussion

On the IFR and the number of infected cases. The actual number of infected individuals in France is probably much higher than the observations (we find here a factor ×8), which leads at a lower mortality rate than that calculated on the basis of the observed cases: we found here an IFR of 0.5% based on hospital death counting data, to be compared with a case fatality rate (CFR, number of deaths over number of diagnosed cases) of 2% on 17 March. Adjusting for the number of deaths in the nursing homes, we obtained an IFR of 0.8%. These values for the IFR are consistent with the findings in [2] (0.66% in China) and [3] (0.9% in the UK). The value of 1.3% estimated on the Diamond Princess cruise ship [4] falls above the top end of our 95% CI. This reflects the age distribution on the ship, which was skewed towards older individuals (mean age: 58 years), among whom the IFR is higher [3,4].

The objective of our study was to estimate the IFR based on early data, before large scale surveys become available. By late April, new data and preliminary studies are available and can be compared to our results. An antibody study in New York released on 24 April 2020 shows about 14 percent tested positive, corresponding to 2.7 million cases, to be compared with the 271,000 confirmed cases and a statewide total of 15,500 deaths. This corresponds to an IFR of 0.6%. In France, another preliminary study conducted by Pasteur Institute [16], and based on a joint estimate from French data up to 14 April (hospital death counting data) and from the Diamond Princess cruise ship data finds an IFR of 0.5% thus confirming our result. By 28 April, the number of deaths from COVID-19 in Lombardy (Italy) is 13,575 (source: Ministero della Salute), for a population of 10 million people, showing that the IFR is at least 0.14%. On the other hand, in South Korea where the number of detected cases rapidly reached a plateau, suggesting a small proportion of undetected cases, the ratio between the number of deaths and the number of positive cases is 244/10,752 ≈ 2.3% (source: Johns Hopkins University Center for Systems Science and Engineering [5]), which can be considered as an upper bound for the IFR, though overestimated.

If the virus led to contaminate 80% of the French population [3], the total number of deaths to deplore in the absence of variation in the mortality rate (increase induced for example by the saturation of hospital structures, or decrease linked to better patient care) would be 336,000 (95%-CI: 192,000–537,000), excluding the number of deaths in the nursing homes. This estimate could be corroborated or invalidated when 80% of the population will be infected, eventually over several years, assuming that an infected individual is definitively immunised. It has to be noted that measures of confinement or social distancing can decrease both the percentage of infected individuals in the population and the degree of saturation of hospital structures.

On the value of $R_0$. The estimated distribution in France is high compared to recent estimates (2.0–2.6, see [3]) but consistent with the findings in [17] (2.24–3.58). A direct estimate, by a non-mechanistic method, of the parameters $(\rho ,t_0)$ of a model of the form ${\widehat{\delta }}_t=e^{\rho (t-t_0)}$ gives $t_0=$ 36 (5 February) and $\rho =$ 0.22. With the SIR model, $I^{\prime }\left(t\right)\approx I(\alpha -\beta )$ for small times ($S\approx N$), which leads to a growth rate equal to $\rho \approx \alpha -\beta$, and a value of $\alpha \approx 0.32$, that is to say $R_0=3.2$, which is consistent with our distribution of $R_0$. Note that we have assumed here a infectiousness period of 10 days. A shorter period would lead to a lower value of $R_0$.

On the uncertainty linked to the data. The uncertainty on the actual number of infected and therefore the IFR are very high. We must therefore interpret with caution the inferences that can be made based on the data we currently have in France. In addition, we do not draw forecasts here: the future dynamics will be strongly influenced by the containment measures that will be taken and should be modelled accordingly.

On the sensitivity of the results with respect to the fixed model parameters. We deliberately chose a parsimonious model with a few parameters to avoid identifiability issues. However, we needed to fix some parameter values. In particular, we assumed a mean duration of the infectious period ($1/\beta$) of 10 days. A much larger infectious period of 20 days (corresponding to the median period of viral shedding found in [12]) would lead to a much larger basic reproduction number $R_0=4.8$ (but still within the range 1.4–6.49 described in [18]) and a factor $\times 15$ between the reported cases and the actual number of cases. However, our main result on the IFR would remain unchanged (0.5%). We also assessed the sensitivity of the inference with respect to the prior knowledge, by proposing a set of more informative uniform prior distributions than the set specified in the main text. Overall, this prior modification does not significantly influence the posterior distributions; see Appendix A.

On the hypotheses underlying the model. The data used here contain a limited amount of information, especially since the observation period considered is short and corresponds to the initial phase of the epidemic dynamics, which can be strongly influenced by discrete events. This limit led us to use a particularly parsimonious model in order to avoid problems of identifiability for the parameters. The assumptions underlying the model are therefore relatively simple and the results must be interpreted with regard to these assumptions. For instance, the date of the introduction $t_0$ must be seen as an efficient date of introduction for a dynamics where a single introduction would be decisive for the outbreak and the other (anterior and posterior) introductions would have an insignificant effect on the dynamics.

A more complex epidemiological model of the COVID-19 epidemic in China has been proposed in [19], with an infectious class divided into several compartments (asymptomatic individuals, unobserved symptomatic infectious and observed symptomatic infectious). The authors use this model in [20] to make forecasts on the cumulative number of cases in China, while taking into account management strategies. In these two studies the authors emphasise the importance of being able to estimate the fraction of infectious cases that are not observed in order to forecast the dynamics of the epidemic. Our study, though based on a simpler SIR model, shows that this fraction can be estimated based on early data.