Estimating the COVID-19 infection fatality ratio accounting for seroreversion using statistical modelling.

Background: The infection fatality ratio (IFR) is a key statistic for estimating the burden of coronavirus disease 2019 (COVID-19) and has been continuously debated throughout the COVID-19 pandemic. The age-specific IFR can be quantified using antibody surveys to estimate total infections, but requires consideration of delay-distributions from time from infection to seroconversion, time to death, and time to seroreversion (i.e. antibody waning) alongside serologic test sensitivity and specificity. Previous IFR estimates have not fully propagated uncertainty or accounted for these potential biases, particularly seroreversion.
Methods: We built a Bayesian statistical model that incorporates these factors and applied this model to simulated data and 10 serologic studies from different countries.
Results: We demonstrate that seroreversion becomes a crucial factor as time accrues but is less important during first-wave, short-term dynamics. We additionally show that disaggregating surveys by regions with higher versus lower disease burden can inform serologic test specificity estimates. The overall IFR in each setting was estimated at 0.49-2.53%.
Conclusion: We developed a robust statistical framework to account for full uncertainties in the parameters determining IFR. We provide code for others to apply these methods to further datasets and future epidemics.

Introduction

One of the most contested statistics during the coronavirus disease 2019 (COVID-19) pandemic has been the infection fatality ratio (IFR): the proportion of those infected who will go on to die from that infection. In the first general wave of the pandemic, estimates of the overall COVID-19 IFR ranged from <0.01 to 2.3%, with a review combining estimates across studies reporting an overall estimate of 0.68% (0.53–0.82%)[1-3]. In addition, an analysis using pooled data from national serologic surveys to estimate age-specific IFRs found that the IFR rose steeply with age, ranging from <0.01% in those aged under 30 to 7.3% in the 80 and older age group[2], broadly consistent with previous estimates[4-6]. IFRs are expected to vary across populations due to: the age distribution of the population, the distribution of infection across age groups, access to healthcare resources, the prevalence of underlying health conditions in the population, biological sex, and other factors. In addition, the overall population IFR may differ depending on the magnitude of outbreaks in care-home settings, where mortality has often been high[7]. As a result, heterogeneity is expected between locations and reflecting this variation is paramount for an accurate representation of the global COVID-19 IFR.

Estimating the IFR requires two key pieces of information: data on deaths and data on the number of infections in the population. Although there are challenges with quantifying and defining COVID-19 deaths, these data are widely reported and one of the more reliable indicators of COVID-19 burden in countries with good testing and reporting systems. However, determining the cumulative number of people infected in a population has proved to be far more challenging. Testing capacity has often been limited and many infections are asymptomatic[8], which makes laboratory confirmed symptomatic case numbers a poor estimate of infection attack rates. As a result, serologic tests (detecting antibodies) have been used to estimate cumulative infections among populations. These tests have several limitations: (1) tests rely on a humoral immune response and will miss infections that do not mount a detectable antibody response or recent infections where antibodies have not yet developed; (2) antibodies naturally wane over time, which can lead to seroreversion (defined in this context as an individual with a confirmed infection and positive serologic test later testing negative); (3) tests will produce imperfect results (i.e. sensitivity and specificity are <100%). Many published studies reporting IFRs did not account for uncertainty in serologic test sensitivity and specificity, nor delays from onset to death and onset to seroconversion (although there are exceptions[2,9-11]) and the possibility of seroreversion has not usually been considered (again with exceptions[11]). Failing to account simultaneously for these factors could potentially lead to biased estimates of the IFR in directions that are hard to predict.

Here, we develop a novel flexible Bayesian statistical framework for estimating the IFR that accounts simultaneously for all the factors listed above. We show that accounting for these factors is critical in accurately estimating the IFR, and that seroreversion starts to significantly affect IFR estimates some months after the start of the pandemic. Similar to previous studies, we find that although overall IFR estimates vary substantially, with age-specific IFRs demonstrating a nearly log-linear pattern. From these updated calculations, we also show that early IFR estimates were relatively accurate despite not incorporating seroreversion. Our method and open-access code provide a tool for analysing IFR using further serologic datasets in the future.

Methods

Crude and test-adjusted IFR estimates

The crude IFR was calculated by dividing the number of observed cumulative deaths at the serologic study midpoint by the cumulative number of infections at the same time point. The number of infections was estimated as the observed seroprevalence multiplied by the population size, plus COVID-19 deaths occurring up to the midpoint of the serosurvey to avoid survival bias. The 95% confidence intervals on the crude IFR were calculated using a Monte Carlo sampling approach, where the uncertainty in the seroprevalence was propagated by drawing 100,000 values of the expected seroprevalence based on the binomial distribution (i.e. the number of test-positives given the total tested). For Denmark, Italy, and Sweden where only the seroprevalence and confidence intervals were reported (i.e. counts of test-positives and total tested were not available) intervals were logit-transformed and used to calculate variances directly. Test-adjusted simple IFR estimates were calculated in the same way, but first adjusting the seroprevalence for the sensitivity and specificity of the serologic test used in the study[12].

Statistical model for estimating IFR

Daily and age-stratified deaths

For individuals who die following infection, we assume that the time from infection to death follows a gamma distribution with shape and rate . If an individual is infected at time then the probability that they die at time is:

We make the simplifying assumption that time is discrete and measured in days, defining to be the probability of death on day given infection at the start of day where :

(the term in the above comes about because we assume infections occur at the start of the day, but deaths can be registered until the end of the day, hence returns a positive value).

Our population is split into different age strata, each with their own probabilities of infection and death. Let there be age groups in total, and let be the proportion of the total population in age group . In the simplest model we would expect infections to occur in a given age group in proportion to the number of people in that group. To allow for variation in age-specific attack rates, and in order to fit to age-specific seroprevalence data, we include a multiplicative attack rate scalar within each group, allowing the final attack rate to be higher or lower than expected from proportions alone. Hence the overall probability of infection in age group , which will be written , is given by:

Once infected, the probability of death in age group (i.e. the IFR in this age group) is defined as . Hence, the overall probability of an individual in age group dying on day given infection on day can be written .

Our raw data do not consist of individual-level outcomes, but rather aggregate counts. Specifically, two marginal distributions were available for each study: (1) daily counts of the number of COVID-19 deaths, summed over all age groups, and (2) the cumulative number of COVID-19 deaths at a single point in time, but broken down by age. Both marginal distributions were fit within a single statistical framework.

Let be the number of new SARS-CoV-2 infections in the population on day . The true infections curve is unknown, and was modelled using an exponentiated natural cubic spline, subject to the constraint that the total number infected (i.e. the area under the curve) could not exceed the total population size . It follows from the definitions above that the number of infections in age group on day is given by , and the number of ultimately fatal infections is given by . The expected total deaths on day , denoted , is obtained by summing over all age groups and all possible times of infection as follows:

The observed number of COVID-19 deaths on day , denoted , is assumed to be Poisson distributed around this expectation:

The likelihood for this part of the model is simply the product of Poisson probabilities over all days in our time series:

Moving on to the second marginal distribution, the expected cumulative deaths in age group up until time can be written:

These expected values are converted into expected proportions of deaths in each age group as follows:

Finally, the observed cumulative COVID-19 deaths up until day , denoted by the vector with elements for , are assumed to be multinomially distributed with these proportions:

This is the second component of the likelihood:

Incorporating serology data

The third data type used in fitting comes from serological studies. For a given individual infected on day we model the probability of having seroconverted by day using the following formula:where is a binary variable that equals 1 if the individual has seroconverted and 0 otherwise. This is equivalent to assuming seroconversion with a constant hazard . Translating to the population level, the expected number of people to have seroconverted by time in age group , denoted , is given by:

This can be translated to an expected proportion via the expression , where is the total population size in age group , such that .

The observed prevalence of seropositive individuals (the seroprevalence) is expected to deviate from this proportion due to both sampling effects and imperfect test characteristics. If is the sensitivity of the test directly after seroconversion, before antibody waning, and is the specificity then the test-adjusted expected seroprevalence, , can be calculated using the classic Rogan-Gladen correction[12]:

Let the total number of people tested on day in age group be denoted , and let the observed number of seropositives be denoted . We model the observed counts as binomially distributed around the Rogan-Gladen-corrected proportion:

Finally, the likelihood for this component of the model is the product of the binomial probability over all age groups, and over all serology study dates :

The full likelihood is the product of the individual likelihood components listed above.

Extension for seroreversion

As part of a sensitivity analysis, we allowed for individuals to serorevert over time under an assumption of natural waning antibodies. We assumed that individuals experience a constant hazard of seroconverting, followed by a probability of seroreverting characterised by a Weibull distribution with shape and scale . Under these conditions, the probability of being seropositive by the end of day following infection on day is given by:

All subsequent steps are identical to those described above in Eqs. (12–15), resulting in an alternative version of the likelihood component .

Model fitting

We used informative priors for key parameters where they were well characterised, such as the delay from symptom onset to death. We fit the model using Metropolis-Coupled Markov Chain Monte Carlo (MC3) using the drjacoby R package (version 1.2.0)[13]. Full details of priors and model fitting are provided in Supplementary Table 1 and Supplementary Methods.

We re-estimated test specificity for serologic studies where regional data were available, by fitting a simplified version of the main model described above to seroprevalence and cumulative regional deaths at the midpoint of the most recent serosurvey, adjusting for age demographic differences within regions using RStan[14] (Supplementary Methods). These estimates were then used as informative priors for the subsequent IFR analyses of each survey.

Convergence of models was assessed by visualising the posterior distributions as well as requiring the Gelman-Rubin’s convergence diagnostic to be lower than 1.1[15]. For the IFR model using MC3, the metropolis coupling acceptance rate between rungs was also examined.

Application to first-wave data

To estimate the time of seroreversion after symptom onset from longitudinal serology data (see above), we fit a Weibull survival model using interval censoring to account for the uncertainty in the observed time of seroreversion. As a comparison to our parametric fit, we also fit a Kaplan–Meier survival curve with interval censoring. Models were fit using the ‘survival’ R-package[16,17]. The ‘survminer’ R-package was used in plotting the Kaplan–Meier survival curve[18] (Supplementary Methods).

Serologic studies were selected from an existing, continuously updated systematic review: the ‘SeroTracker’ dashboard[19]. Estimates of the sensitivity and specificity of the serologic assay were obtained preferentially from validation conducted as part of each serosurvey, rather than external validation (e.g. by manufacturers). We preferentially obtained data on COVID-19 deaths by age and date of death from Ministries of Health and national public health agencies (Supplementary Table 2), and when otherwise not available, used data from the COVID-19 Data Repository by the Centre for Systems Science and Engineering at Johns Hopkins University (JHU CSSE COVID-19 Data) up to August 17, 2020 (accessed September 14, 2020)[20,21]. Similarly, demographic information was extracted from both governmental and non-governmental websites. Ethical approval was not required because the data were publicly available. Datasets are archived on Github[22].

We calculated pooled-IFR estimates using a weighted log-linear regression on the age-specific IFR posterior estimates. Weights were incorporated as the precision from the age-specific 95% credible intervals. Prediction intervals were calculated from the log-normal density function using the mean from the model fit and model variance. Overall pooled-IFR estimates were calculated by standardising to the demographics of representative countries within the low-income country (LIC), low-middle income country (LMIC), upper-middle income country (UMIC), and high-income country (HIC) bracket, respectively[23] (Supplementary Note 4).

Table 1

Overall infection fatality ratio estimates among the included studies.

	Data					Model estimates
Study location	Cumulative COVID-19 deaths	Reported Seroprevalence (dates)	Serostudy Sensitivity (%) (T+/D+ or 95% CrI)	Serostudy Specificity (%) (T−/D− or 95% Crl)	Crude IFR (95 CI%)	Sensitivity (%) (95% CrI)	Specificity(%) (95% CrI)	IFR without Seroreversion (95% CrI)	IFR with Seroreversionβ (95% CrI)
Brazil*	51,179	2.42% (Jun. 04–Jun. 07)	85.14 (81.93, 87.97)	99.72 (99.55, 99.85)	0.99 (0.92, 1.06)	85.28 (82.12, 88.14)	99.76 (99.62, 99.87)	1.03 (0.93, 1.15)	0.99 (0.89, 1.12)
Denmark*	463	2.4% (Apr. 27–May 03)	82.09 (75.51, 87.58)	99.25 (98.94, 99.56)	0.33 (0.23, 0.48)	82.45 (76.11, 87.8)	99.16 (98.73, 99.46)	0.54 (0.38, 1.02)	0.51 (0.36, 0.97)
England*	48,301	5.94% (Jun. 20–Jul. 13)	78.4 (65.68, 88.15)	99.44 (99.11, 99.71)	1.42 (1.39, 1.46)	79.48 (68.74, 88.93)	99.59 (99.34, 99.78)	1.18 (1.02, 1.34)	1.07 (0.84, 1.24)
Italy*, α	34,610	2.44% (May 25–Jul. 15)	96.04 (89.84, 99.05)	99.7 (99.59, 99.79)	2.3 (1.94, 2.72)	96.42 (90.93, 99.13)	99.69 (99.57, 99.78)	2.53 (2.31, 2.78)	2.40 (2.18, 2.63)
Netherlands	5767	5.5% (May 10–May 20)	98.28 (171/174)	99.65 (281/282)	0.6 (0.58, 0.63)	98.23 (95.61, 99.52)	99.83 (99.43, 99.98)	0.62 (0.58, 0.69)	0.59 (0.55, 0.65)
Spain*	28,116	5.27% (Jun. 08–Jun. 22)	81.84 (75.67, 87.01)	98.79 (98.55, 99.02)	1.12 (1.08, 1.16)	84.72 (83.08, 88.36)	99.05 (98.86, 99.21)	1.14 (1.08, 1.22)	1.08 (1.01, 1.16)
Sweden	4992	7.1% (Jun. 08–Jun. 12)	99.36 (156/157)	98.89 (267/270)	0.68 (0.46, 1)	99.28 (97.23, 99.93)	99.17 (98.12, 99.77)	1.02 (0.87, 1.37)	0.98 (0.83, 1.35)
Geneva, Switzerland	262	10.84% (May 03–May 10)	91.16 (165/181)	100 (176/176)	0.48 (0.42, 0.56)	91.47 (87, 94.89)	99.89 (98.82, 100)	0.49 (0.42, 0.59)	0.47 (0.4, 0.57)
Zurich, Switzerland	124	1.59% (May 01–May 31)	90.74 (49/54)	99.89 (5,497/5,503)	0.51 (0.45, 0.58)	91.77 (83.39, 96.89)	99.87 (99.74, 99.95)	0.52 (0.41, 0.67)	0.50 (0.39, 0.64)
New York State, USA*	17,718	12.1% (Apr. 19–Apr. 28)	89.39 (85.57, 92.55)	98.73 (98.15, 99.27)	0.75 (0.74, 0.76)	89.66 (85.9, 92.68)	98.7 (98.05, 99.2)	0.78 (0.73, 0.84)	0.76 (0.7, 0.81)

The data columns (left) contain data and parameters used to calculate the crude IFR, the model columns (right) contain the posterior estimates from the full model. Citations for all the data sources are in Supplementary Table 2. The reported seroprevalences are listed along with the most recent dates for the seroprevalence survey. Cumulative deaths are summed to the mid-date of the most recent seroprevalence survey, and were usually confirmed COVID-19 test-positive patients except in England, which also reported probable COVID-19 deaths (individuals without test results but with COVID-19 on the death certificate). For the six studies with regional data, estimates of specificity and sensitivity were from analysis of regional data: posterior distributions with the median and 95% credible intervals are provided in place of the serologic test validation numbers (*). Model-estimated posterior sensitivity and specificity are indicated for the model with seroreversion, although these estimates were similar for both models (Supplementary Table 5). Overall IFR estimates were calculated by standardising the age-specific IFR estimates according to the inferred age-specific attack rate and the population demography with respect to the age groups used in the model (median, (95% Credible Intervals)). For comparison, the overall IFR estimates calculated by standardising for solely the demography and assuming the same attack rate in each age group are provided in Supplementary Table 3.

Serologic test performance is measured by the sensitivity and specificity (T+ test positive, D+ true positives, T− test negative, D− true negatives).

αSerovalidation data for the Italian serosurvey using the Abbott assay were not validated within the same study; here we used an alternative study testing the same assay.

βAssuming an extreme rate of seroreversion for sensitivity analysis based on the Abbott assay. The true seroreversion rates in these studies are unknown, but are likely less extreme, particularly if the Abbott assay was not used (only the Italy study used the Abbott assay.

Table 2

Pooled estimates of the infection fatality ratio.

Age-band (years)	IFR (%) without seroreversion (95% PI)	IFR (%) with seroreversion (95% PI)
0–4	0 (0, 0.04)	0 (0, 0.04)
5–9	0.01 (0, 0.07)	0.01 (0, 0.07)
10–14	0.01 (0, 0.12)	0.01 (0, 0.11)
15–19	0.02 (0, 0.2)	0.02 (0, 0.19)
20–24	0.03 (0, 0.32)	0.03 (0, 0.31)
25–29	0.04 (0, 0.5)	0.04 (0, 0.48)
30–34	0.07 (0.01, 0.75)	0.06 (0.01, 0.72)
35–39	0.1 (0.01, 1.09)	0.1 (0.01, 1.05)
40–44	0.16 (0.02, 1.54)	0.16 (0.02, 1.47)
45–49	0.25 (0.03, 2.11)	0.24 (0.03, 2.02)
50–54	0.4 (0.06, 2.84)	0.38 (0.05, 2.7)
55–59	0.62 (0.1, 3.75)	0.59 (0.1, 3.56)
60–64	0.96 (0.19, 4.9)	0.92 (0.18, 4.64)
65–69	1.5 (0.35, 6.38)	1.43 (0.34, 6.03)
70–74	2.34 (0.66, 8.31)	2.23 (0.63, 7.85)
75–79	3.66 (1.23, 10.9)	3.47 (1.18, 10.27)
80–84	5.71 (2.26, 14.44)	5.41 (2.16, 13.59)
85–89	8.9 (4.09, 19.37)	8.43 (3.91, 18.21)
90+	17.36 (9.73, 30.97)	16.4 (9.25, 29.08)
Overall (LIC)	0.24 (0.15, 0.43)	0.23 (0.14, 0.41)
Overall (LMIC)	0.4 (0.27, 0.68)	0.39 (0.25, 0.65)
Overall (UMIC)	0.62 (0.41, 1.01)	0.59 (0.39, 0.97)
Overall (HIC)	1.16 (0.79, 1.82)	1.1 (0.75, 1.72)

Bold and Italic values represent the overall numbers at the end.

IFR estimates were calculated by combining study- and age-specific IFR estimates in a log-linear model. The median predicted estimate and corresponding 95% prediction intervals (PIs) are shown above. Predictive intervals were used to express the plausible range of IFRs that can be expected in a new study population, rather than showing our degree of certainty of our estimates with confidence intervals. For the 90+ age group, we assumed a maximum age of 100 years. The overall IFR estimates were standardised by the population structure in a representative low-income country (LIC), low-middle income country (LMIC), upper-middle income country (UMIC), and high-income country (HIC), assuming equal attack rates across age groups.