Statistical techniques to estimate the SARS-CoV-2 infection fatality rate

1 Introduction

A critical task in the context of the SARS-CoV-2 pandemic is to determine the infection fatality rate (IFR), defined as the proportion of deaths among all infected individuals. The IFR is one of several characteristic measures that form the basis for epidemiological models such as [1], which are subsequently used to shape and justify government policy on public health interventions. The basic mean reproduction number R₀ and its full distribution extensions characterize the multiplicative process rate on the production side of the epidemic and the IFR is defined on the decay side of the individual infections. In this paper we concentrate our efforts on the IFR and do not consider the estimation of R₀, even if the total risk and harm caused by the virus is driven by both factors.

A more widely reported metric is the case fatality rate (CFR), defined as the ratio of the number of deaths attributable to SARS-CoV-2 and the number of documented infections. The reported CFR can vary significantly between regions and countries, due in a large part to the varying ability of local authorities to comprehensively screen all suspected cases of infection. Further, there is significant evidence for a substantial cohort of asymptomatic carriers of SARS-CoV-2, an important subpopulation that is only partially (if at all) represented by the CFR. Studies to determine the IFR are typically supported by cross-sectional seroprevalence surveys in population samples to also account for asymptomatic (and mildly or atypically symptomatic) infections. Hence, the IFR is considered to be a more reliable variable than the CFR for the purposes of epidemiological modelling.

There are numerous serological surveys underway, or recently concluded, that aim to estimate the IFR for SARS-CoV-2. It is crucial that these studies consider all relevant sources of uncertainty, of both statistical and systematic origin, to establish the confidence intervals on individual estimates of the IFR. This in turn allows for meaningful comparisons between (and potential combinations of) independent results. The reported confidence interval in Ref. [2] appears to neglect the dominant uncertainty in the study, namely the variance in the underlying statistical distribution used to model the number of fatalities. This omission may have implications for policy decisions made on the basis of estimates from these types of study. An accurate and unbiased estimate of IFR can be also difficult to obtain during an evolving epidemic due to various time delays that must be correctly accounted for: from exposure to the virus to the onset of symptoms following the incubation period, to the reporting of a positive case following a PCR test (polymerase chain reaction), to the development of sufficient antibodies to be identified by a test (seroconversion), to the recording of a fatality.

The body of research on SARS-CoV-2, and the resulting COVID-19 disease, grows at an astonishing rate. The number of studies from which an estimate of the IFR can be drawn is now plentiful and meta-analyses are now being performed that aggregate information from several sources. For example, Ref. [3] considers 25 estimates of the IFR that are derived from studies of various types, including serological surveys and epidemiological modelling. It reports a point-estimate of 0.68 [0.53, 0.82]% for the IFR, with the interval stated at a 95% confidence level.¹ However, the study acknowledges a high degree of heterogeneity in the individual estimates. Indeed, there are many factors that may influence the results of the individual studies. Ref. [3] comments on the variability in the quality and rigour of the surveys, and the lack of peer review for some studies. Also noted is the challenge of determining the IFR from serological surveys for which there is an absence of an associated, reliable fatality count. There are many local factors, such as population density, demographics, and health, and the ability of healthcare services and government policy to protect the local population. Perhaps one of the most important factors is the age demographic of a population, as the IFR appears to be highly dependent on age. Accurate estimates of IFR stratified by age are highly desirable in the near future.

The meta-analysis reported in Ref. [3] relies on a common ‘method of moments’ method by DerSimonian and Lard [4] to provide a point-estimate of the IFR. We review several approaches to combine results from individual studies of the IFR into a single estimate. These include the method of moments and a normal likelihood based classic meta-analysis, a joint likelihood combination, and methods to combine full probability densities such as optimal transport and the product or mean of Bayesian posteriors.

In Section 2, we introduce the Gangelt field study [2], which we use as an example in the first part of the paper. Section 3 reviews several methods that can be used to determine a confidence interval for the IFR, based on a single binomial proportion, a ratio of binomial proportions, a profile likelihood ratio, Monte Carlo simulation, and Bayesian constructs. Section 4 compares the confidence intervals from the various methods. Section 5 presents a time-dependent deconvolution calculus that accounts for intrinsic time delays through the determination and application of a correction factor (and associated uncertainty) to the IFR estimate. Section 6 provides first a general perspective on how multiple data points can be combined, before a set of seroprevalence studies from around the globe are introduced, which are then used as concrete examples for the aforementioned combination techniques. Finally, Section 7 summarises this work and concludes by providing “best practice” recommendations in the context of confidence interval calculations for the IFR of the SARS-CoV-2 coronavirus.

2 The Gangelt field study

A sero-epidemiological study of a small German community exposed to a super-spreading event was undertaken to determine the IFR [2]. The municipality of Gangelt is located in the district of Heinsberg in the German state of North Rhine-Westphalia. The municipality reported unusually high levels of SARS-CoV-2 infections following local Carnival festivities held on 15^th February 2020. Strict social distancing measures, which included an advisory curfew, were introduced on 26^th February to suppress further infections. The Carnival festivities held in Gangelt are attended predominantly by people living locally and the community is relatively closed, with low levels of tourism and travel.

An estimate of the IFR is obtained from the double ratio where the raw fatality rate r_F is the ratio of the number of fatalities n_F and individuals n_P in a given population, and the raw infection rate r_I is the ratio of the number of infections n_I∧T identified in a representative cross-sectional sample of tested individuals n_T. Assuming the test sample of n_T individuals is representative of the population under study, in terms of both demographic and seroprevalence measures, the product n_Pr_I provides an estimate of the total number of infected individuals in the population n_P. The authors of Ref. [2] identified the Gangelt municipality as a unique opportunity to accurately estimate the IFR, because of its stable, relatively isolated population and an appreciable number of fatalities arising from the high level of infections present in its community.

The Gangelt municipality has a population of n_P = 12597. Inhabitants were randomly polled to participate in testing for both SARS-CoV-2 virus RNA (PCR) and antibodies to identify the number of both present and past infections, respectively. The number of inhabitants that were tested and satisfied all survey requirements is n_T = 919, and the number of identified infections (n_I∧T) in this sample is n_I∧T = 138. Within the duration of the study, the number of deaths recorded in the Gangelt municipality that were associated with a SARS-CoV-2 virus infection is n_F = 7. Hence the study yields a raw infection rate of 15.0% and an IFR of 0.37%.

Following the application of corrections for various identifiable biases and the associated statistical and systematic uncertainties, Ref. [2] reports r_I = 15.5 [12.3, 19.0]%. Hence, the total number of infected individuals in the Gangelt municipality is estimated to be 1950 [1550, 2390]. The reported interval for the IFR of 0.37% is [0.29, 0.45]%, which can be used to estimate the number of infected individuals in other places with population characteristics similar to that used in the Gangelt study. Several other key findings are reported in Ref. [2], which are beyond the scope of this note. We restrict the discussion in this note to the reported confidence interval for the IFR estimate.

A total of 6575 deaths associated with SARS-CoV-2 were recorded in Germany by 2^nd May 2020. Assuming the Gangelt sample population is representative of the German population as a whole, the IFR can be used to estimate the total number of infections – at some point earlier in the pandemic – that led to the reported death toll on 2^nd May 2020. The study yields an estimate of 1.8 [1.5, 2.3] million infected individuals in the German population. Using the methods described below, we determine a confidence interval (CI95) of [0.9, 3.7] million on the estimate of 1.8 million infected individuals, based on a simple method (the Wilson score) applicable to a single binomial proportion.

3 Methods

Estimation of confidence intervals of parameters with statistical tests goes as follows. The test statistic for the null hypothesis H₀ is known to asymptotically follow an analytic distribution, which is taken as an approximate proxy to the problem. The exact solution is available only in special cases. The analytic approximation can be relaxed with more modern numerical bootstrap and Monte Carlo techniques which provide the most appropriate tools when data needs to be also re-weighted, propagated through a chain of analysis algorithms or manipulated in more complex ways.

3.2 Generating non-asymptotic test statistics

Confidence intervals without relying on the asymptotic χ²-statistic can be obtained using Monte Carlo. A common choice with optimal properties in this context is the likelihood ratio based construction or ‘ordering principle’ of acceptance sets, which has been studied theoretically and computationally since the Neyman-Pearson lemma, see e.g. Refs. [15–19]. It relies on the (exact) duality between statistical tests and confidence intervals. Now, the well known brute-force algorithm to construct the so-called (Neyman) confidence belt is as follows.

1. The parameter θ₀ space is discretized or randomized uniformly.

2. For each point value of θ₀, a sample of toy MC values is generated by drawing from the corresponding sampling model density, for example . Drawing from the underlying density is typically strictly necessary (only) if no analytic or parametric pdf is available, or when their evaluation is difficult.

3. Ordering principle (acceptance region) in the sample space: A. Generated MC points can be used without any intermediate test to construct the exact empirical CDF quantiles, similar to Clopper-Pearson central intervals. B. Targeting optimal properties induced by Neyman-Pearson lemma, the (profile)-likelihood ratio such as Eq. 3.8 can be calculated as an intermediate step to provide an ordering, using for different . The exact empirical CDF quantile points are then taken using the generated test statistics. Return to Step 2.

4. Finally, the parameter uncertainty region is obtained by intersecting the Neyman belt in ‘orthogonal direction’ at the observed value of k, i.e. in the parameter space. This is an essential part of the construction and implements the (test) inversion, by taking a union of acceptance sets. Special care must be taken at this point while looping over the acceptance sets, because their union may yield sometimes discontinuous topologies depending on the specific ordering or acceptance principle of Step 3.

For more information on this, see Appendix L. We illustrate this algorithm in Figure 1 for a single binomial uncertainty using the exact (MC based) log-likelihood ratio test statistic [16, 17, 19] as described above, where we see that the asymptotic χ²-approximation is quite good in this case, but the relative discrepancy grows reasonably large with small numbers. This MC based construction was introduced to high energy physics by Feldman and Cousins [19]. Clopper-Pearson procedure is also compared, which gives slightly different intervals than the exact LLR based. The basic property of the asymptotic LLR approximation is that the acceptance region threshold given by the χ²-distribution quantile, is independent of the local θ₀ value, unlike the exact Monte Carlo driven LLR test statistic and the CP construction. These both are exact procedures in a sense that their coverage is always larger than or equal to 1 − α. The lattice structure shows why < and ≤ operations make a difference with discrete numbers but not with continuous parameters. When constructing Figure 1 for LLR based variants, in Step 3., we used t ≤ t_C, where t is the log-likelihood ratio test statistic and t_C the cut value (χ²-quantile or Monte Carlo based). Similar care is required with the ‘vertical direction’ in Step 4.

• Download figure
• Open in new tab

Figure 1.

(Gangelt setup) CI95 confidence belts using a Monte Carlo based (exact) likelihood ratio test, an asymptotic χ²-based likelihood ratio test and a central binomial pdf based construction (Clopper-Pearson). Each horizontal gray slice is computed with MC, the rest can be obtained without random numbers.

The approach described here is in principle generic, however, the construction in higher parameter dimensions can be technically challenging. This depends on the construction of the likelihood functions (parametric vs non-parametric) and computational complexity of the Monte Carlo procedures. With several nuisance parameters, the profile likelihood approximation is typically used in Step 3. to keep the whole approach feasible. For K nuisance parameters, the profiled likelihood ratio complexity may scale (naively speaking) linearly 𝒪(K) using e.g. simultaneous stochastic gradient search at each likelihood evaluation point, whereas a full hyperbelt scan is exponentially hard 𝒪(N ^K), where N is the number of discretization points in each dimension. The multiple minima need to be in principle handled by the profiling procedure (an example in Section 3.4). Analogous computational challenges arise in Bayesian solutions with high dimensional integration, typically implemented with various Markov Chain MC and variational approximations.

3.4 Profile likelihood ratio

The profile likelihood method splits the parameters into two groups: true parameters of interest θ and nuisance parameters ξ, and maximizes the full likelihood over the nuisance parameters where sup denotes the supremum, the least upper bound, which is almost the same as the maximum but takes into account the possibility that the likelihood cannot be evaluated exactly at that point ξ. The main idea behind profiling is the dimensional reduction over the nuisance parameters, which then allows one to infer the uncertainty on θ by formally proceeding as with a usual likelihood, for example by using the score test or the likelihood ratio test which are asymptotically equivalent. Solutions based on the score test for the two binomial case have been proposed in [31, 32]. We shall now derive the likelihood ratio test based solution.

Let us parametrize r ≡ p₁/p₂ and write down the joint likelihood function for two independent binomials

This re-parametrization does not involve the change of variables formula (Jacobian), because the likelihood as a sampling function and its volume normalization is over k₁ and k₂, which we left intact. We are interested in the parameter r and treat the parameter p₁ as a nuisance parameter, which we profile out by finding a value for p₁ which maximizes the joint likelihood for every single value of r. This procedure is in principle readily generalized to arbitrary number of nuisance and true parameters of interest, which however is only a formal statement. Possible singularities depend on the exact type of nuisance parameters and models. A practical problem in higher dimensional cases is also the parameter optimization problem itself.

The profile log-likelihood ratio test statistic follows here asymptotically

The fact that χ²-asymptotics holds also for the profile likelihood inference is a non-trivial result, but propagates from Wilks’ theorem under certain assumptions. After maximizing Eq. 3.27, the profile log-likelihood ratio is where the maximum likelihood estimates are and the local profile extremum solution conditioned at point r₀ has two roots where the negative (−) branch gives the right solution in our problem. It is not guaranteed that every profile likelihood problem is differentiable and has a closed form solution, but this turned out to be the case here. The profile log-likelihood ratio is illustrated in Figure 2, which has asymmetric 95% confidence interval endpoints around the maximum likelihood value.

• Download figure
• Open in new tab

Figure 2.

(Gangelt setup) On the left, a profile likelihood ratio based IFR confidence interval.

3.5 3-dimensional Monte Carlo simulation

Notation

Q68, Q95 are used for pure density quantiles and CI68, CI95 for a parameter estimate confidence (frequentist) or credible (Bayesian) intervals.

This elementary simulation approach starts with three Bernoulli random numbers: tests ∼ B_T, infections ∼ B_I and deaths ∼ B_F, which are together per person modelled as a 3-dimensional Bernoulli distribution. To remind, the Bernoulli distribution is the underlying distribution behind the binomial distribution, which turns into a Poisson distribution when p is small and n is large. Thus, this approach is ab initio in this hierarchy of distributions. Now in general, a D-dimensional Bernoulli requires 2^D − 1 free parameters. However, we do not have enough measurements here to constrain all the parameters. To simplify this problem, we factorize B_T to be independent of B_I and B_F

That is, tests do not (hopefully) affect infections or fatalities. This leaves us with one and two dimensional sub-problems which require together 1 + 3 parameters. The two-dimensional problem can be parametrized directly with four probabilities of (B_I, B_F)-binary combinations, which sum to one. Another parametrization uses the expectation values E[B_I], 𝔼[B_F] and the correlation coefficient λ[B_I, B_F] ∈ [−1, 1]. We use both in order to be able to sample correlated Bernoulli variables in the direct (multinomial) basis, which satisfy by definition the conservation of probability, and to give an interpretation of the problem in the correlation basis. The formulas are given in Appendix A and B.

The four ‘dynamic’ parameters of the simulation are fixed in the correlation basis according to their maximum likelihood values where the count variables and their associations to the underlying sets of tested T, all infected I in the city and all fatal F are as described in Section 2. Here, one must pay attention to Eq. 3.34, where the study sample infection rate is assumed to be representative in the simulation for the whole population by assuming homogeneity between samples. A systematic uncertainty could be associated here. By choosing in Eq. 3.36 the maximum possible positive correlation coupling (see Appendix B), the simulation output in Table 2 reproduces the event count observables, which enter as the input variables. That is, its value is fixed by data. Also, a boundary condition is used which states that no fatalities happen without getting infected. This forbids the combinations 1 and 5 in Table 2 from appearing. The total number of people n_P = |P| is kept fixed for each MC run. We have also fixed the test sample size n_T = |T| to be constant in these simulations, to follow more closely the Gangelt setup, but turning on the Bernoulli fluctuations is implemented in the code as an option. However, with the given event counts the difference is not significant for the IFR. Type I and II errors of tests are not simulated here. Once calibrated, including their effect as a post processing step is trivial with two free parameters and an additional coin flipping per tested person. For more details, see Appendix G.

Using the generated Monte Carlo samples, arbitrary observables such as the IFR are computed by simply counting numbers from an (8 × N_MC)-dim matrix, here N_MC = 10⁶, and accumulating numerically the relevant point estimates such as mean values and percentiles. Note that this aggregated matrix is the fully sufficient statistic and contains all simulation information, due to binary random variables. These results are given Table 2. The simulated distributions for the infection fatality rates are shown in Fig. 3, which illustrates the non-trivial Dirac’s comb discrete characteristics of the problem, but also the finite sample smearing effect on the extrapolation estimate. The smeared IFR-distribution is calculated from the test samples and the reference IFR-distribution from the inaccessible (full) population statistics, which are both obtained simultaneously in the simulation. See Appendix A for the exact definitions.

• Download figure
• Open in new tab

Figure 3.

(Gangelt setup) The IFR distributions obtained from Bernoulli Monte Carlo simulations on the left. The blue distribution simulates the full population with a complete test coverage and the orange distribution is obtained using the finite test sample, and its extrapolation to the full population using Eq. 2.1. The same folding effect is visualized on the right in 2D, where the full population realized IFR value is on the horizontal axis and the test sample driven extrapolated estimate on the vertical axis.

To this end, we may summarize that the power of this simulation is the ‘full phase space’ modelling of partially overlapping sets of tested, infected and fatal, which is not possible with independent binomial ratios. The simulation is based on describing the most elementary classical stochastic process involved, namely correlated Bernoulli coins. The optimal confidence interval estimator can be based on the simulations as described in Sec. 3.2 where each simulation with fixed input parameters simply generates a sample for a single null hypothesis H₀. Alternatively, faster bootstrap approximations can be used.

3.6 Bayesian inference

In the exact Bayesian inference within two independent binomial distributions, we keep the number of observed events k₁, k₂ as fixed numbers, also n₁, n₂, and calculate the joint posteriori density for the binomial parameters p₁ and p₂. Their joint posteriori density is described with a product of two beta distributions given generic beta priors Beta(α_i, β_i) and binomial likelihoods. The derivation of this and priors are given in Appendix C. Given this joint posterior, the ratio r ≡ p₁/p₂ density is obtained via change of variables such that p₁ = ry and p₂ = y. Writing down the Jacobian determinant gives |∂(p₁, p₂)/∂(r, y)| = |y|. Then we substitute these new variables in Eq. 3.38, include the determinant and integrate over y

Using Mathematica, we obtain for this integral a representation where is the regularized Gauss hypergeometric ₂F₁ function and B is the Euler beta function. The regularized version is . Equation 3.40 represents the master formula, which can be evaluated numerically with high precision special function libraries and credible intervals can be obtained with standard numerical integration techniques. However, we found that numerically it is easier to use directly Eq. 3.39.

The results are given Figure 4, where the joint posteriori distribution includes the credible regions (CR), which encapsulate 68 and 95 percent of the probability mass. The shape is constructed according to the natural contour lines. The right figure shows the ratio posteriori density and the corresponding cumulative distribution by using two different prior distributions. The solid line is obtained using the non-informative Jeffreys prior Beta(1/2, 1/2), which is invariant under coordinate transformations. It is proportional to the square root of the Fisher’s information determinant , where the determinant represents abstract information volume (here in one dimension the determinant is trivial). The results with dashed lines are obtained using a unit flat prior, which is not completely non-informative and results in slightly larger values. Its maximum (mode) gives numerically the same estimate for the IFR as the simple ML estimate.

• Download figure
• Open in new tab

Figure 4.

(Gangelt setup) On the left, the Bayesian joint posteriori distribution using Jeffreys prior. On the right, the ratio posteriori density in red and the cumulative distribution in blue for Jeffreys and flat priors.

Nuisance parameters and systematic uncertainties

Bayesian framework allows one to add nuisance parameters and systematic uncertainties into the formulation. For example: the death counts k₁ may be need to be scaled with a parameter γ due to time delays. Note that scaling the parameter p₁ instead is ambiguous, which is seen using the binomial pdf and by computing the Fisher information matrix, which will turn out to be singular. Which means that the corresponding parameter estimation problem would be rank deficient. Similarly the positive test counts k₂ may be multiplied with another scale λ. Using auxiliary measurements or prior judgement, the uncertainty information on the nuisance parameter is often modelled using a Gaussian prior π_γ(γ; µ_γ, σ_γ) with fixed µ_γ, σ_γ. However, with a positive definite scale, using the gamma prior could be more suitable, although practical difference can be small. This is applied on the Bayesian inference master formula as an additional prior constraint term and by replacing k₁ → γk₁ everywhere. Computationally, each nuisance parameter requires typically an additional integral when marginalizing the posterior and in the normalization (Bayes denominator). For more details, see Appendix D.

4 Estimator comparisons

Table 3 shows the numerical results for different confidence interval estimators, using count data of the Gangelt study and Figure 5 similarly, but as a function of death counts (rather than for just the observed F = 7). Clear outliers in the group of these estimators are the normal (Wald) test based and the bootstrap percentile estimator. Their interval is shifted toward smaller IFR values, which is especially visible with CI95 intervals. The Wald test will also give negative (unphysical) values at small F. This can be expected from their mathematical construction. The impact of the bias correction and acceleration for the bootstrap is clear. The rest of the estimators yield numerically similar values for the Gangelt input data and small differences are more easily seen from Figure 5. Compared with the Monte Carlo simulations of Figure 3, the Bayesian distributions of Figure 4 are completely smooth because the observed event counts are considered fixed and the continuous binomial parameters p₁, p₂ are considered random. In contrast, the simulations are closer to a frequentist inference, because the model parameters are fixed and the discrete event counts are random. Though technically speaking, the underlying simulation can be used within both philosophies.

• Download figure
• Open in new tab

Figure 5.

(Gangelt setup) Relative interval widths compared for a subset of the estimators with running death counts, other numbers kept fixed. At the top, CI68 endpoints and at the bottom, CI95 endpoints. The likelihood ratio is based on χ²-approximation and the bootstrap is constructed using the percentile approximation.

Coverage probability simulations and interval widths are given in Appendix L for the single binomial based estimators. These illustrate significant undercoverage of the Wald test based estimator, the conservative coverage of the Clopper-Pearson based estimator and the ‘bracketing’ coverage behavior of the Wilson score and Bayesian estimators. The likelihood ratio with an asymptotic χ² approximation has behavior similar to the Wilson score, but significantly undercovers at very small values of p, similarly to the the Wald test.

5 Time evolution

The previous discussion in Sections 2 and 3 is only fully applicable under the asymptotic time t → ∞ limit or instantaneous action, i.e., no time delays. To be concrete, by time asymptotic we mean the tail of a single epidemic outbreak and neglect additional possible complications (immunology evolution, different viral strains) which result from overlapping epidemic ‘waves’. However, purely mathematical overlap is automatically handled by our description, that is, we do not assume any specific epidemic shape for the time-series input. In this section, we briefly outline how the IFR estimation and its uncertainty is implemented during an evolving epidemic with finite time delays. For the interested reader, further details of our time evolution study can made available upon request.

A combined double delay effect can be summarized with one ratio function where K_F is the delay kernel (pdf) from infections to deaths, K_S is the delay kernel from infections to seroprevalence (antibodies) and the symbol ∗ denotes a linear time-lag convolution integral. These kernels are extracted from data by fitting them typically with Weibull or log-normal distributions, see e.g. supplementary material of the Geneva study in [33] and references there, where a similar convolution calculus was used. The kernels are given in Appendix H, which shows explicitly the expected delays. The relative differences between K_S(t) and K_F (t) drive Eq. 5.1. The calculus here is general, and one may replace the antibody type tests with PCR type tests by replacing the delay kernel K_S. In that case an additional multiplicative effect to include is the ‘viral shedding’ (loading) period probability, i.e., how long an infected person gives a positive test result. That evaporation factor can be neglected with antibody based serology (seroreversion effect), if the half-life involved is large enough on the scale of epidemic. For some specific antibodies, this may not be the case. The denominator of Eq. 5.1 can be extended to incorporate it, see Appendix H for that construction.

Basic numerical integration is used here for the convolutions. Time t denotes the time of the seroprevalence determination and Δt ≥ 0 denotes how many days later the population cumulative death count is taken. Because our problem here is essentially an inverse problem, the underlying cumulative infection count Î(t) can be estimated computationally by regularized deconvolution of the reported positive cases dC(t)/dt of PCR viral tests. Daily counts need to be used in the inversion instead of cumulative counts, to conserve all information. Although even if one assumes a constant reporting rate, Î(t) can be estimated only up to an unknown scale (probability), which however fortunately cancels in Eq. 5.1. In principle, one may also use the reported daily death counts dF (t)/dt to obtain the deconvolution inverse estimate of Î(t), albeit the statistics might be too limited. For technical details about the deconvolution algorithm, see Appendix H. The algorithm is based on non-negative linear least squares with Tikhonov smoothness regularization. Regularization is needed, because basically all naive inversion procedures always amplify the counting fluctuations (noise). No fine structure recovery is needed, thus smoothness is a good functional prior in this problem.

By using Eq. 5.1, the delay corrected non-equal time IFR estimate is now where F (t) is the population cumulative death count and Î_S (t) is the population level seroprevalence (extrapolated) estimate Î_S (t) = n_P × n_I∧T/n_T. Here n_P is the population size, n_I∧T is the number of infection positive in the demographically randomized test sample and n_T is the test sample size. By non-equal time we refer here to the shift by Δt, which can be optimized after the seroprevalence test.

No delay correction is needed, if t or Δt are chosen (or happen to be) with certain lucky values. This depends on interplay between three factors: 1. the delay kernel K_F (t) of deaths, 2. the delay kernel K_S(t) of antibodies (seroprevalence) and 3. the cumulative epidemic curve Î(t). In our SARS-CoV-2 case, using kernels from [33], the kernels give a functional shape for ψ(t, Δt) which peaks above one for small t, then decreases below one, and asymptotically approaches one when t → ∞. However, eventually the antibodies will vanish from the body (seroreversion), so realistic times scales must be used, also for other obvious reasons.

The systematic uncertainty estimates should include perturbation of the kernels and the estimated Î(t) function, most easily studied via toy Monte Carlo, propagated through the deconvolution algorithm and Eqs. 5.1 and 5.2. The full procedure to estimate ψ(t, Δt) is illustrated in Figure 6, based on kernel data from [33] and Switzerland data from [34]. For the uncertainties, we used approximately 20 % Gaussian equivalent uncertainties in the Weibull kernel parameters and fluctuated the input data with Poisson uncertainties, propagated via toy Monte Carlo. A good re-projection of the deconvolved Î(t) to deaths is observed shape wise, as a ‘closure test’. We emphasize that this closure test would be trivial, if the daily death count dF (t)/dt would have been used as the algorithm input. But we used the reported daily PCR cases dC(t)/dt as the input, so the result is non-trivial. The absolute normalization for Î(t) and is matched to data for the visualization, because it is not obtained as a part of the procedure. The scale function ψ(t, Δt) has larger uncertainty and larger values earlier in the epidemic, which is natural.

• Download figure
• Open in new tab

Figure 6.

Switzerland data deconvolution results (gray dotted), forward re-projection (red dotted) and in the bottom, the double delay scale function ψ(t, Δt) contour estimates (CI95). The overall scale normalization of estimates (dashed) is matched with the measured (solid) functions of cases C(t) and deaths F (t), for visualization purposes. The top figure illustrates that our regularized deconvolution estimate of the infection rate dÎ(t)/dt is stable and that the delay kernels given in Appendix H are realistic – a good match is obtained between the measured cumulative death counts F (t) and the forward projected shape in the middle figure. Based on these observations, the obtained delay scale function ψ(t, Δt) values in the bottom figure can be considered realistic. Raw data is from [34].

6 Combination analysis

In general, the value of IFR for a given individual is dependent on a number of factors such as age, sex, viral load, diet, genetics etc. For example, estimates of IFR in SARS-CoV-2 were found to have a strong age dependence in [33] and [35], varying over two to three orders of magnitude as a function of age. Mathematically, let there be a multivariate IFR function which takes as an input a random human feature vector X and returns the expected probability to die Y if infected.

In Bayesian modelling, the random feature vector is often composed into measured X, and identified to be important but not necessarily measured variables Z, such that,

Given enough, well sampled data, the function IFR(X) ≡ P (Y |X) can be modelled using simple histograms, (conditional) logistic regression, deep learning learning techniques or other methods. However, given the naturally limited data available early on in a pandemic, the IFR function can be integrated over its dependents to obtain a local expected IFR for a particular population, where f_j(X|city) is the normalized sampling density of the population. A city is chosen here to represent a realistic system size in terms of sample statistics, which can be considered as independent from other systems.

Any significant difference in the population densities f (X|city), between cities, will yield different empirical ⟨IFR⟩_j values. If the IFR function has a strong dependence over a particular feature, a bias with respect to the other city will be present. When comparing studies implemented in different cities, this intrinsic sampling bias can be compensated in two ways:

1. physically, a priori, using carefully designed sampling (selection) of the population

2. mathematically, a posteriori, using an inverse weight or sampling function, which can be modelled using detailed demographic statistics of the test sample and the city.

Using either the strategy 1 or 2, one must also choose a reference population density or a ‘standard template’ to represent a typical demography. This provides the golden reference for the sampling procedures. Once this sampling or stratification bias is accounted for, a truthful comparison of IFR values can be obtained. The expected raw global value without any re-sampling schemes can be modelled with a mixture density model as where the weights are w_j = N_j/∑i N_i and N_j is the total number of people in the city.

6.2 Combination strategies

The datasets from the studies in Table 4 can be considered simultaneously to study the global IFR. In this section, we describe several strategies to do so. In general, we assume that the IFR values in each dataset are integrated values over all physical properties, and that the sample used in the study is chosen such that the single IFR value obtained will be representative of the whole population. Furthermore, we assume statistically independent infection rates in each city, because the systems are physically long distance isolated and do not contain common infection sources. An exact decomposition of the sources of variance in the IFR, e.g. due to different demographies both within (local) and across (global) populations, is clearly not uniquely solvable. However, meaningful statistical estimates can be obtained. For comparisons, see [46].

Fundamentally, we can write down an illustrative global-local-random-sampling additive model to study the datasets as a whole,

The following is a brief overview of the methods used in Section 6.3. In the following, r_j represents the observed IFR values for j = 1…K independent studies. For some early work on the comparison analysis of similar experiments, we refer to Cochran [47].

Method of Moments

This non-parametric estimator for the meta-analysis was in its simplest form proposed by DerSimonian and Lard [4]. This method aims to determine the mean and variance of the parent distribution of r. The mean is determined as, with the global variance or ‘heterogeneity’ given by [48] where the test statistic is and represents the estimated variance within each study. We use a two step, iterative procedure in which and are first determined using, and then both and Δ² are updated by setting

Local convergence is obtained typically after a few iterations. The asymptotic standard error on is given by which provide Wald test-like confidence intervals. This method does not yield uncertainty on . In general, a finite sample error is to be expected directly based on the sampling error in the study-specific variance estimates . Several weighting scheme variants of the DL estimator have been proposed. See Ref. [48] for a recent comparison study, where the two-step DL estimator was found to be among the best, but the original DL performed weakly in some scenarios.

Normal Likelihood model

Hardy and Thomson [49] used parametric normal-normal hierarchy with sampling densities After integrating out the latent θ_j, this gives a normal marginal distribution . The= total joint log-likelihood with K contributing studies is

It is easy to change the underlying sampling densities e.g. to a log-normal which takes effectively into account the physical boundary r > 0. The maximum likelihood solution can be obtained via standard optimization techniques or by iterating the following equations

The two-dimensional confidence region on these parameters is obtained with the log-likelihood ratio

By profiling, we obtain the individual confidence intervals for r and Δ². It is also straight-forward to extend this normal-normal model to fully Bayesian hierarchies as described in [50]. In that case Markov Chain MC sampling is typically used for obtaining the posterior density, which requires special technical care and is most easily dealt with specialized libraries.

As an alternative to the simple normal likelihood based, the so-called Restricted Maximum Like-lihood (REML) method was introduced by Patterson and Thomson [51] for unbiased estimates of variance components in linear mixed models. See Ref. [52] for a detailed derivation within the correlated and full multivariate formulations.

Wasserstein-Fréchet mean

The Wasserstein metric barycenter or the Fréchet mean [53, 54] is an optimal transport (OT) based approach, which solves the optimization problem where is the optimally combined new density under the p-Wasserstein metric W_p, which is a geodesic transport metric in the space of densities. See Appendix I for more details and Ref. [55] for statistical properties. The weights can be taken as , if the solution is taken to be (inversely) proportional to the sample variances, for example. The solution is found by discretizing the 1D-posterior densities for each dataset and constructing the barycenter as an average in the inverse CDF space of quantile functions. We tried also a Sinkhorn iteration based algorithm using an entropy regularized transport cost formulation known as Bregman projections [56], with minimum regularization set such that a numerically stable output was obtained. However, this approach resulted seemingly in an over-smoothed output which is mathematically expected due to the entropic approximation.

Arithmetic mean of posteriors

The arithmetic mean of posterior densities is which has a mixture model interpretation and the weights can be taken as in the optimal transport case. The density interpolation properties can be more limited compared to the optimal transport case, which can be crucial if the idea behind combining the posterior densities is to find one common data generating distribution. For multimodal densities, the mean estimator is an inclusive, probability mass covering estimator.

Product of posteriors

The normalized product (geometric mean) of posterior densities is where provides the re-normalization. This approach is known in machine learning as the product of experts model by Hinton [57], where several simple models are combined to ‘vote’ together. It is also called logarithmic pooling. Thus for multimodal densities, the product estimator is an exclusive, single mode seeking estimator. By using information theory, this approach can be derived using e.g. the so-called α-divergence of Chernoff, which is a generalization of the standard Kullback-Leibler (KL) divergence (relative entropy). The work by Amari [58] shows how the product is an optimal solution to a divergence risk minimization problem with α = 1 (reverse KL) and that the arithmetic mean of Eq. 6.19 is also a minimal risk solution, but under its dual α = −1 (forward KL).

Joint likelihood ratio

This approach uses a product over the likelihood for each independent dataset, where X_j = {k₁, n₁, k₂, n₂}_j is the j-th dataset. We compute the profile likelihood ratio test statistic of Eq. 3.27 independently for each dataset. The combined IFR maximum likelihood value and confidence intervals are then obtained easily by comparing the total product (sum) of Eq. 6.21 against the -distribution quantiles as described in Section 3.4.

We can summarize that the optimal choice of the combination method depends on the underlying assumptions, which are encoded by the implicit or explicit algebraic, information theoretic or probabilistic aspects of the method. Here the first class of combiners is more inclusive, the second class more exclusive, in their output decision. See Appendix J for some illustrative properties of the probabilistic risk functions, which may be used for formal motivations.

As already mentioned, the methods we considered can be roughly classified into two scenarios for combining the data from the different studies. The first scenario aims to determine a parent distribution, from which each IFR observed for a particular city is assumed to be sampled from. In the second scenario, we assume δ_j ≡ 0 for all j in Equation 6.4 which represents a model in which the the true global and local IFR values are the same. Under this assumption, the individual measurements of IFR can be combined to provide a more precise but possibly overoptimistic measurement of IFR. In the following section, we use the Method of moments, the Normal likelihood, the Wasserstein-Fréchet mean, and the Arithmetic mean of posteriors methods for the first scenario, and the Product of posteriors and Joint likelihood ratio for the second scenario.

6.3 Results

The individual observed IFR result for each dataset is given in Table 6, where the 95% confidence (credible) intervals are obtained using the Bayesian estimator described in Section 3.6. The results are given using four different choices of the fixed time delay Δt, and with an adaptive delay determined with the inverse machinery described in Section 5. In general, some datasets are strongly dependent on the chosen read-out delay. The reason for this is the underlying local epidemic evolution and its time derivatives. The adaptive delays are given in Table 8 where the confidence intervals are based on propagating Poisson fluctuations and kernel uncertainties through the whole deconvolution chain with Monte Carlo sampling as described in Section 5. This estimated uncertainty δγ on the death count scale is included as a Gaussian prior in the Bayesian estimator as described in Section 3.6 and Appendix D, whereas the fixed read-out delay estimates here are ‘bare’ and do not include time scale uncertainties. This difference is manifest in the width of the individual densities. Measurements done in the very end of the local epidemic, generate larger optimal read-out delay values and correspondingly smaller delay uncertainties on the death counts because the daily increase in deaths has reached the slow asymptotic regime. The test inversion systematic scale uncertainty δλ is included as a Gaussian prior affecting the IFR denominator, numerically larger for low prevalence datasets (see Table 8).

View this table:

• View inline
• View popup
• Download powerpoint

Table 6.

Infection Fatality Rate (IFR) [%] estimates (CR95) for each dataset, where columns denote different fixed read-out delays Δt [days] and a data adaptive (inverse solved) read-out Δt ← ψ(t, Δt) using global kernels. The credible interval takes into account the statistical counting uncertainty in the double ratio via Bayesian posterior estimator under non-informative Jeffreys prior, the test inversion related uncertainty δλ and also the delay uncertainty δγ in the case of adaptive delay estimation in the rightmost column.

View this table:

• View inline
• View popup
• Download powerpoint

Table 7.

Infection Fatality Rate (IFR) [%] combined results with different strategies. The columns are the distribution mode, mean and 68 and 95 quantile intervals. The intervals have different interpretations and formal definitions depending on the underlying method (see text for details). N.B. Joint LLR is without systematic scale uncertainties δγ and δλ of Table 8 and includes only statistical counting uncertainties. The fixed Δt delay estimates are without including the systematic scale uncertainty δγ on the death count.

View this table:

• View inline
• View popup
• Download powerpoint

Table 8.

Left column: the optimal read-out delay in days for each dataset based on the deconvolution inversion and its uncertainty. The Gangelt data includes no detailed time-series for deconvolution, thus Δt = 1 is used instead. Center column: deconvolution and Monte Carlo propagated relative systematic scale uncertainty δγ ≡ σ[n_F(Δt)]/n_F(Δt) on the cumulative death counts due to causal time-delays at time t + Δt. Right column: Type I and II test error inversion relative systematic scale uncertainty δλ constructed with an error propagated uncertainty and a pure binomial Wilson uncertainty on the corrected prevalence fraction p₂, according to Eq. G.6. Test sensitivity and specificity values are as given in Table 4, which are used in the error propagation described in Appendix G, to obtain individual δλ values given here. These systematic uncertainties are finally applied with the Bayesian priors as described in Appendix D.

The posterior distributions from each city are presented in Figures 7, 8 and 9 for the choices of adaptive, Δt = 7 and Δt = 14 days, respectively. The sensitivity to the time delay effects highlights the importance of including these effects into any complete study of the IFR. Philadelphia yields significantly larger IFR peak values than the rest, while Los Angeles, Santa Clara and Finland all yield much smaller IFR peaks. This could point to a relatively strong IFR dependence on the local population, and further highlights the importance of sampling the population within an individual study. The Gangelt study is approximately in the middle and fairly constant with choice of Δt. It is worth noting that the kernels used in adaptive delay inversion are the same (global) for each dataset, which makes it locally biased for PCR test based prevalence data (Stockholm, Iceland) and due to local reporting delays.

• Download figure
• Open in new tab

Figure 7.

On the left, Bayesian posteriori densities for different datasets under non-informative Jeffreys prior. On the right, the method of moments (MoM) and the normal likelihood (NL) model based point estimates of visualized with Gaussian uncertainties, the combined densities using unweighted and variance weighted optimal transport (OT), the mean of posteriors (SUM) and the normalized product (PROD). With the optimal Δt ← ψ(t, Δt) read-out delay solution from deconvolution.

• Download figure
• Open in new tab

Figure 8.

Same as Fig. 7, but using the read-out delay Δt = 7 days.

• Download figure
• Open in new tab

Figure 9.

Same as Fig. 8, but using the read-out delay Δt = 14 days.

The results of the different combination strategies are given in Table 7 and in Figures 8 and 9, obtained by combining the individual Bayesian posterior densities. The meta-analysis method of moments (MoM), and normal likelihood (NL) based strategies are presented assuming a Gaussian distribution, with the mean being and the standard deviation representing the uncertainty on , obtained using the methods described in Section 6.2. The global dispersion parameter Δ and its uncertainty are shown instead in Figure 11 together with the IFR parameter r, using the NL model full likelihood information. The optimal transport is fully non-parametric, thus its output distribution does not explicitly try to disentangle the IFR and global dispersion like components and their individual uncertainties. Instead, it yields distributions which look in their functional form similar to individual distributions. The inverse variance weighted variant prefers smaller values of r, which is well expected with this set of data. However, often there can be other more suitable weighting strategies, based on e.g. some auxiliary risk minimization and the available domain knowledge.

The mean of posteriors method gives a very different distribution to the others. This method is sensitive to the datasets chosen as input and in general will not give a single representative (unimodal) distribution when the individual distributions have small overlap regions. In particular this is seen in the close by double peak structure around r = 0.2 in Δt = 7 days results, which disappears in the Δt = 14 days result. This double peak structure is caused mainly by the New York distribution being narrow enough and moving significantly between these two choices. The simplest distribution peak structure is obtained with the adaptive delays, shown in Figure 7, which has only one strong peak between r = 0.1 and r = 0.2.

The product of posteriors and joint likelihood methods yield distributions which are much more narrow than the others. Figure 10 show the combination results using the joint likelihood method, where the delay Δt = 14 days gives slightly smaller joint log-likelihood ratio than the delay Δt = 7 days. However, this cannot stand on its own as a method for determining a good effective global Δt, because some datasets are invariant under the choice of Δt, thus their IFR values would have a pivotal role. Based on our time-delay calculus, individual datasets all have different delay scale functions ψ(t, Δt) and are evaluated at different t values, which give optimal local Δt values shown in Table 8.

• Download figure
• Open in new tab

Figure 10.

Same as Figs. 8 and 9, but using the joint LLR combination and without scale uncertainties δγ and δλ. The vertical lines are the maximum likelihood estimate of the IFR and its CI95 confidence intervals.

• Download figure
• Open in new tab

Figure 11.

The normal likelihood model parameters r (IFR) and Δ (global dispersion) uncertainty regions (CI68, CI95 and CI99) based on contouring the corresponding log-likelihood 2D-function against distribution quantiles. The parameter maximum likelihood values are shown with black dots.

Finally, the high similarity between the two product methods is a natural consequence of the underlying assumption of these strategies that the IFR is a global quantity and that individual studies can be combined to improve the precision of the measurement. We caution however that this is a strong assumption and in particular highlight the fact that several of the individual distributions show tension with the product like combinations, indicating that the assumption may be unjustified.

7 Conclusions

The demographic averaged infection fatality rate (IFR) of SARS-CoV-2 in Western societies has been estimated here to be approximately 0.4 [0.2, 0.8] % (Q95), when using a fixed Δt = 7 day read-out delay between the prevalence determination and public death counts and the optimal transport based posterior combination of individual datasets. In general, our methods included the Bayesian double ratio based binomial counting uncertainty, deconvolution of the underlying time-series for the optimal read-out delay determination, systematic uncertainties in the death counts and systematic uncertainties in the corrected positive test counts. This estimate is a factor of 3 − 5 larger than the IFR of a typical seasonal influenza, if one assumes its often referred value of ∼ 0.1 %. Our result is numerically similar to analyses e.g. in [59] or [3], but differs in the methodology. However, even if we constructed our methodology rigorously, one should not take our estimate as an ultimate measurement of the IFR. In addition of requiring better control of the underlying details of collected data, we identified also other factors which more extensive analysis would take into account such as demographic differences, population counts used in the normalization and regional time-delays.

Our IFR estimation is based on several random sampled seroprevalence determination datasets combined with different statistical techniques, such as Bayesian estimation of counting uncertainties and modern algorithmic optimal transport driven probability density fusion. We recommend the Wilson estimator for fast but reasonable confidence interval estimation of binomial counting experiments and the Bayesian double ratio estimator for more extensive counting uncertainty estimates in the context of the IFR, because it provides access to a full posterior distribution. Also additional effects can be more easily incorporated with the Bayesian formulation. When analysing multiple datasets, the choice of data combination tools depends on the underlying scenarios. If it is reasonable to assume that there exists one global integrated IFR value, the product of posterior distributions or the joint likelihood method are perhaps the most natural approaches to use for improving the individual estimates. However, if that is not the case, then the Wasserstein-Fréechet mean (optimal transport) of posteriors, the classic meta-analysis models and the linear mean of posteriors can be more suitable, as we also reasoned with results from the information theory.

An accurate estimate also requires careful consideration of the significant time delays observed between infection and death. The required time delay convolutions necessitate the extraction (i.e. fitting) of delay kernel functions that may need to be locally adjusted because of differences in health care and administrative procedures. Poorly estimated kernels used within (de)convolution results in unknown bias, naturally. A more transparent solution is to always show in addition results with several fixed delays such as one or two weeks, as presented here. However, we provided and analyzed the necessary inverse problem methods for advanced, close to optimal delay corrections and demonstrated this machinery with data using fixed global kernels. The best solution experimentally is a cross-sectional seroprevalence trial that minimizes time-domain effects, namely, not done too early in the epidemic and which is tightly localized (not spread) over time, if possible.

Based on studies in Stockholm [35] and Geneva [33], and global comparisons [45], binning (stratifying) over age seems to be a crucial selection variable for the SARS-CoV-2 IFR, as expected. All statistical methodology developed here can be applied also in a stratified analysis. However, age stratification is not enough; although it gives strong ordering in IFR, it does not provide a proper explanation of the underlying dynamics. Possible body response differences and uncertainties in the antibody type tests are crucial factors, as well as the crucial extrapolation to the total population level. A recent study has found a new type of genetic defect risk in type I interferon (IFN) path-ways, inducing a life-threatening COVID-19 pneumonia [60], but not being an active mechanism with influenza viruses.

From a larger perspective, a direct one-to-one comparison e.g. with seasonal flu is non-trivial, because for that the population has more natural immunity and there are seasonal vaccinations. However, an indirect comparison is possible and for understanding the total risk and harm on the society, one should understand the multiplicative reproduction differences between these viruses. A virus with a seemingly relatively small IFR can still be of high risk, if it is easily transmitted. A driving factor might not just be the mean R₀, but also independently the variance and tails of the transmission chain multiplicities. This possibly overdispersed case (compared to e.g. Poisson) is often modelled with a negative binomial distribution, which in physics describes at a phenomenological level the charged particle final state multiplicities of high energy proton-proton collisions. By using sharp analogies, tools from high energy physics can be useful on modelling and analyzing the epidemic production side of the problem.

Physically, the ultimate solution for the future is to increase massively the testing capabilities for real-time monitoring. This basically requires new type of non-invasive personal health care technology, perhaps based on promising new techniques such as CRISPR based diagnostics [61]. More measurements is not just a requirement to obtain minimally extrapolated IFR estimates and reliable values for the epidemic parameters such as R₀ and understanding the critical role of multiplicative fluctuations, but more crucially, to obtain control of the epidemic with minimal lock-down measures. All the analyzed and developed tools here were constructed essentially from the first principles, and thus these methods should stay highly relevant also in the future.

Notes

Open source Python code which reproduces all algorithms, figures and tables shown here, and beyond, is available at: github.com/mieskolainen/covidgen.

A Infection fatality rate observable

Definitions of the infection fatality rate estimators are where F, I, T, P denote the sets of fatal, infected, tested and all people in the city, respectively. The number of elements of a set is denoted with | · | and the intersect of two sets with ⋂. To show that the construction is consistent, consider the limit where all people in the city are tested by substituting T → P in Eq. A.2. Then the limited statistics IFR coincides with the full statistics IFR by construction, because also always holds that I ⋂F ≡ F and I ⋂ P ≡ I. These definitions are purely formal and assume perfect test sensitivity & specificity and no time-delays. These necessary corrections are discussed in other sections of this paper.

The implicit assumption made in the extrapolation is that the infections observed in the test sample represent truthfully the stochastic infection process in the full sample. In essence some ergodicity (time-average equals ensemble average) and sample homogeneity must be assumed.

B Sampling two dimensional Bernoulli random numbers

Using the expectation values 𝔼[X], 𝔼[Y] and the correlation coefficient ρ[X, Y] between two correlated Bernoulli random variables X and Y, the direct (hypercube) basis parametrization is

The sampling of vectors (X, Y) is now multinomial, such that [(0, 0), (0, 1), (1, 0), (1, 1)] ∼ [P₀, P₁, P₂, P₃] are four corners of the hypercube. Any multinomial distribution sampling algorithm can be used.

C Details of the Bayesian estimator

D Systematic uncertainties via Bayesian priors

Additional systematic uncertainties on counts k₁ and k₂, are applied by multiplying and integrating over the Bayesian posteriori ratio IFR formula of Eq. 3.39 with where G(x,µ, σ) is a normal density, ϵ a small positive scalar and the integrals are computed numerically. These additional Gaussian distributed parameters model multiplicative scale corrections γ and λ on the death counts k₁ and on the positive test counts k₂, respectively. The triple integral gives the posteriori ratio probability density up to the overall normalization, which is obtained numerically. The normal prior densities here can be replaced with gamma densities, for example.

The mean values are taken µ_γ = µ_λ = 1, typically, if the the counts are corrected prior this formula. The 1-sigma uncertainties on these corrections are described by δγ and δλ, which come from auxiliary procedures or calibration measurements. In our case, δγ is obtained via Monte Carlo error propagation of the deconvolution procedure and its kernel uncertainties in Section H, and δλ as described in Section G. One needs to pay attention to possible double counting of statistical uncertainties in Equation D.1, when estimating these parameters.

E Credible and confidence intervals

F Acceptance set ordering principles

The optimal frequentist confidence interval acceptance set construction, used in the inverse construction of the Neyman confidence belts, can be derived briefly as follows [16].

1. There is a one-to-one mapping between tests and confidence intervals.

2. Uniformly most accurate (UMA) confidence region minimizes the probability of false coverage.

3. By using Property 1, UMA set is found by inverting the uniformly most powerful (UMP) test.

4. According to the Neyman-Pearson lemma [63], the likelihood ratio test is the UMP when both

5. the hypothesis H₀ and alternative H_A are simple (not composite). The UMP also exists for a composite H_A, if the underlying distributional family has the so-called monotone likelihood ratio property. In the most general case, no UMP test is guaranteed to exist.

Several other acceptance set constructions or ordering principles also exist, such as the shortest expected length and various pdf based orderings, perhaps optimal under some very specific condition such as certain interval topology. Also randomized intervals can be constructed, but which are mostly used only in theoretical analysis of (discrete) problems.

G Type I and type II test error inversion

Let p = P (V₊) be the true viral prevalence of the population, let q = P (T₊) be the fraction of positive tests in the test sample. Let specificity be s ≡ P (T₋|V₋) = 1 − α = 1 − P(type I error) and let sensitivity be v ≡ P (T₊|V₊) = 1 −β = 1 − P(type II error). Using alternative terminology, α is known as False Positive Rate and 1 − β as True Positive Rate. These symbols should not be mixed with the parameters of the Beta priors, to be clear. The following derivation uses pure probability calculus, without specifying the underlying density or mass functions.

The four different conditional probabilities can be combined under the Bayes’ rule with the law of total probability expanded for the true prevalence

Using these, a well-known inverse estimator (see e.g. [64]) for the true prevalence is which has a physical solution , if and only if

False Positive Rate α ≤ Positive Test Fraction q ≤ True Positive Rate (1 − β).

Otherwise the problem is physically ill-posed. Especially the FPR lower bound is problematic when the viral prevalence is low. As a simple estimate of the related uncertainty, we can use the first order Taylor expansion (error propagation) with q, s, v taken independent. We get where are the individual 1-sigma uncertainties squared. The first one is driven by thebinomial counting, the two other by the uncertainty in the laboratory calibration of the test error rates. Instead of using dichotomic (binary) test output decisions and Eq. G.3, alternative inversion strategies can be based on a test-by-test weighted inversion, according to the conditional probabilities of Bayes’ rule and Expectation-Maximization (maximum likelihood) iteration of the prevalence fraction. This requires that the test provides a probabilistic output (e.g. multivariate analysis). Different strategies should be simulated with Monte Carlo sampling

H Regularized non-negative deconvolution

The deconvolution here is implemented as a non-negative least squares with Tikhonov regularization. We found this classic approach to be by far the most stable of standard methods in this problem, including regularized Fourier space methods and early stopping regulated maximum likelihood EM-iteration (Richardson-Lucy). The EM-iteration driven formulation assumes Poisson noise, which in principle should be more optimal, however, the explicit regularization properties of the method shown here seemed to play a bigger role.

The regularized least squares solution for the discretized infection rate x ∼ dI(t)/dt is obtained by inverting the linear convolution equation Ax = y, by minimizing where y is the measurement vector and λ_R controls the regularization strength. The measurement vector is constructed from the daily PCR infection counts y ∼ dC(t)/dt, where the vector domain is extended (padded) with zeros before the first counts, in order to be able to describe the ‘pull-back’ of deconvolution without a limiting boundary. The matrix A is the convolution operation Toeplitz matrix constructed from the corresponding discretized kernel function K(t). The auxiliary vector is x₀ = 0 in our problem formulation. The regulate the solution smoothness (curvature), we use a finite difference second order derivative matrix

Other typical options for L are the identity matrix and first order derivatives. The minimization is done through an active set method [65] which enforces the necessary Karush-Kuhn-Tucker (KKT) constrained optimization conditions. To be able to use standard optimization algorithms with the regularization term included, we use an augmented matrix formulation

Thus the regularization is fully explicit here. We use minimal parameter values for λ_R yielding smooth inversion results without oscillatory behavior and remark that the statistical uncertainties of the inverse estimate are affected by the regularization procedure, due to the bias-variance trade-off. The regularization adds a small bias term into the solution, and correspondingly suppresses the statistical fluctuations. This makes the statistical uncertainty properties of inverse estimates non-trivial.

The kernel extraction from data itself is its own problem, typically approached using e.g. Kaplan-Meier type non-parametric estimators [66] and functions generalizing the basic exponential (memory-less process) delay kernel, such as the Weibull pdf. Sequential delays are easy to model via cascaded convolutions, but the factorization and identifiability of the component kernels in terms of the underlying physically independent delay sources is not necessarily possible. The fitted kernels are shown in Figure 12, which are also causal such that they are defined only for t > 0.

• Download figure
• Open in new tab

Figure 12.

Delay convolution Weibull kernels K(t) fitted using the mean and standard deviation values given in [33], where λ (scale) and k (shape) denote the corresponding Weibull pdf parameters. The individual delays are: K_I→O is from the infection to the symptom onset (incubation period), K_O→C is from the symptom onset to the case report, K_O→S is from the symptom onset to seroconversion (antibodies) and K_C→F is from the case report to death. The combined delays are: K_C = K_I→O ∗ K_O→C is from the infection to the case report, K_S = K_I→O ∗ K_O→S is from the infection to seroconversion and K_F = K_I→O ∗ K_O→C ∗ K_C→F is from the infection to death. The combined kernels are solved by numerical convolution of the individual kernels, with tilded variables denoting the resulting Weibull parameters.

I Wasserstein optimal transport

Using standard notation, the p-Wasserstein metric [67] for p ≥ 1 is given by where d(x, y) is the basic cost between two points x and y, for example d(x, y) = lx − yl. The so-called transport map T : µ ↦ ν, maps a measure density µ to another measure density ν over the space χ. The set of all possible couplings is Γ(µ, ν), which has marginals µ and ν, with a realization γ(dx, χ) = µ(dx) and γ(dy, χ) = ν(dy). The case p = 2 and χ = RD has a unique minimum solution. In one dimension D = 1, the metric can be written as where U (z) and V (z) are the cumulative distribution functions (CDF) of µ and ν, and the comparison in Eq. I.2 is between inverse CDFs (quantile functions).

J Optimality under risk functions

It may be tempting to choose only one of the estimator results. However, optimality of this decision depends on the risk function definition. To ease out with possible interpretations, here we list shortly some typical risk functions. Let our parameter of interest be θ, the random variable of data be X, the decision function (estimator) be δ(X) and the loss function be ξ(θ, δ(X)), which encodes our cost definition. Beyond these probabilistic risks, there are related principles in information theory, such as the minimum description length (MDL) [68] and other formulations e.g. in economics.

K Overview of systematic uncertainties

This section is a general summary of possible unknowns.

L Coverage simulations

• Download figure
• Open in new tab

Figure 13.

On the left, simulations of a single binomial proportion confidence interval CI95 coverage as a function of the binomial probability p. On the right, the corresponding average confidence interval CI95 relative widths (endpoints). Each row for n number of binomial trials. The likelihood ratio is with χ²-approximation.