Estimating Lower Bounds for COVID-19 Mortality from Northern Italian Towns

1 Summary of Results

Towns in the four northern Italian provinces of Milano, Bergamo, Brescia, and Lodi were heavily infected with the coronavirus in the first few months of 2020. Some towns had substantial excess mortality:

• Nembro, 147 extra deaths, population 11,526, excess mortality 1.275 percent

• Alzano Lombardo, 102 extra deaths, population 13,655, excess mortality 0.75 percent

If the towns had been fully infected during this period the excess mortality would provide an estimate of mortality due to COVID-19. Since the towns may have been less than fully infected these provide lower bounds.

The excess mortality provides a valuable bound for COVID-19 mortality and it would

Generalizing the idea to multiple towns provides better precision in the estimates, and opens the way for estimating mortality by age and sex. A Poisson count model provides the statistical framework to pool data across towns, estimating age- and sex-specific mortality rates, while controlling for differences across towns in the degree of excess mortality. There are 612 towns with deaths by age and sex available from ISTAT for January 1 through April 15 for the years 2017-2020 in the northern provinces of Milano, Bergamo, Brescia, and Lodi – a population of 5,753,296 or 98 percent of the regional population.

Town-specific random effects control for differences across towns in the level of excess mortality which presumably results from differences in infection rates, and also provide a method for estimating mortality for the most-infected towns. The upper tail (99th percentile) of the town distribution corresponds to the highest-mortality towns – presumably those most infected by the coronavirus. The 99th percentile out of 612 towns is roughly the 606th town, so at the upper tail but not simply the most-infected town.

Table 1 shows overall excess mortality of 1.09 percent (95% CI 1.06-1.14),¹ an IFR higher than many published estimates (e.g. Ioannidis [2020] (surveying published studies) reports 0.02 percent to 0.40 percent; CDC [2020] uses an “upper bound” for disease severity of case fatality rate 1 percent and asymptomatic case percentage of 50 percent, implying an IFR of 0.5 percent; Grewelle and Leo [2020] does estimate a higher global IFR (1.04 percent), but this is equal to our lower bound). The IFR shows a dramatic age and sex profile: men age 65-74 have excess mortality 12-times higher than younger men and 2.7-times higher than women of the same age. Applying the age- and sex-specific excess mortality to the U.S. age profile produces a bound for the U.S. overall IFR of 0.87 percent – Table 2.

These estimates of all-cause excess mortality provide particularly firm and useful bounds for COVID-19 mortality because they do not depend on the assumptions often necessary with other measures. They do not depend on ascribing cause of death or determining exposure or proportion asymptomatic infections.

These excess mortality estimates provide lower bounds for COVID-19 IFR in realistic circumstances, but we need to discuss their interpretation. The Infection Fatality Rate (IFR) is a key variable in epidemiological modeling and is usually considered a constant, but in reality IFR can vary with circumstances. There are three channels through which a COVID-19 epidemic will affect mortality. The first is death caused directly by COVID-19 but under “good medical treatment.” We might call this the “clinical” or “best-case” IFR and this is what is usually meant when discussing the IFR. The second is higher mortality still caused directly by COVID-19 but due to ineffective or lesser medical care that may exist during the initial outbreak or disorganized medical conditions. The third is excess mortality from other causes due to disorganized medical treatment, for example cardiac cases not seeking or receiving adequate care – not technically part of the IFR but still mortality resulting from the disease.² The three channels combined we might label the “initial response,” “real-world,” or “worst-case” IFR – due not only to the clinical conditions of infection but also to real-world treatment challenges and the “fog of war” in fighting a new disease.

Reports on the response of hospitals and ICU beds in Lombardy (e.g. Manca [2020], Grasselli et al. [2020]) indicate that the medical system was stressed but not overwhelmed and hospitals responded with large increases in ICU capacity. From anecdotal reports it appears that increased mortality due to the third channel (non-treatment of non-COVID conditions) was small, and due to the second channel (ineffective COVID care) was probably not substantial. Nonetheless, as time, knowledge, and medical practice progresses this initial response IFR will likely decrease.

The estimates in Table 1 provide benchmarks against which to compare parameters used in epidemiological modeling. In this respect it is instructive to examine recent CDC simulations (CDC [2020]). Re-weighting our age-specific estimates for US population (as shown in Table 2) we find 0.87 percent overall and 3.95 percent for age 65+. The CDC scenario 5 (“current best estimates”) assumes 35 percent asymptomatic cases and case fatality rates of 0.4 percent overall, 0.2 percent for 55-64, and 1.3 percent for age 65+, implying IFR of 0.26 percent overall, 0.13 percent for 55-64, and 0.845 percent age 65+. These are substantially lower than the observations from northern Italy and raise the question of whether U.S. mortality will in fact be improved relative to Italy by a factor of three.³

By estimating excess mortality using fixed effects for demographic and town effects, we can examine town-by-town mortality while still pooling data across demographics and towns. Table 3 shows six towns between the 98th and 99th percentile with the highest predicted overall mortality.⁴ There are a number of towns in the tail of the overall mortality distribution (for exampleGazzaniga, Castiglione d’Adda, and Nembro) that have moderate populations and statistically significant estimated fixed effects. This indicates that the estimates from the tail of the random effects estimator shown in Table 1 are not simply a statistical artifact:

Estimates for the upper tail of the random effects estimator and for individual towns tell the same story of high excess mortality for highly-infected towns.⁵

The “Actual” excess mortality for males aged 55-64 and 65-74 also show why pooling across the towns is valuable. Towns with low population will have substantial random variation, particularly for groups with low mortality. Even if the overall town mortality is enough to provide reliable overall mortality estimates, the age-specific mortality measures may not be reliable. The town of Valbondione provides an example. The town has a population of 606 with only 57 men aged 55-64. With one death in 2019 versus none in 2020 the town has a large negative excess mortality, when in fact the variation in number of deaths is likely due to random chance given the small size of the relevant population.

2 Background and Summary Statistics

There are 612 towns (Comunes) having data provided by ISTAT in the four provinces of Milano, Bergamo, Brescia, and Lodi. Table 4 shows summary statistics. The median town population is 4,108. Median observed mortality for 2020 was 0.58 percent, with the median excess mortality 0.25 percent.⁶

There is substantial variation in 2020 mortality across towns, and substantially more than in 2019. The interquartile range is 0.39 percentage points, versus 0.19 percentage points in 2019. This variation in observed mortality, which carries over to the excess mortality, is presumably due to the variation across towns in overall infection rates. This variation, and particularly the upper tail of the distribution of excess mortality across towns, is what we can exploit to calculate lower bounds for IFR.

Using the raw mortality for estimating IFR bounds by specific age and sex groupingsx presents problems, however. Some towns have low population and for groups with low underlying mortality rates, purely random variation will sometimes produce a small number of excess deaths. In towns with small populations such excess deaths translate into unusually high mortality rates, and these high rates will bias upwards the upper tail of the across-town distribution. Table 4 shows the problem. For the median town Males 15-54 have low overall mortality and zero excess mortality. In contrast, the upper tail (the 99th percentile) shows a high observed excess mortality of 0.33 percent. This is not a reliable estimate. The Poisson mixed effects model provides more reliable estimates.

Figures 1 and 2 shows the location of the towns. Figure 1 is a simple quantile map, showing towns by their predicted 2020 excess mortality. Figure 2 clusters towns by the excess mortality and contiguity (Euclidean distance). The maps indicate that mortality was high around Bergamo and particularly in the valleys north-east of Bergamo. Milan and surroundings appear to have been more lightly infected.

• Download figure
• Open in new tab

Figure 1: Quantile Map for Predicted Excess Mortality by Town in the Italian provinces of Milano, Bergamo, Brescia, and Lodi

• Download figure
• Open in new tab

Figure 2: Clustered KMeans Map for Predicted Excess Mortality by Town in the Italian provinces of Milano, Bergamo, Brescia, and Lodi

Predicted excess mortality rate (from fixed effect estimation) clustered with GeoDa using KMeans clustering with 5 clusters, geometric centroids weighted at 0.33. The cluster center excess mortality (after standardizing) is: 1: 0.0038; 2: 0.0008, 3: 0.0008; 4: 0.0070; 5: 0.0126

3 Poisson Count Mixed Effects Model

I use a mixed fixed and random effects Poisson count model (a form of general linear mixed model – see Raudenbush and Bryk [2002] Chapter 10 or Stroup [2013]). Data are grouped by:

• i – 612 towns

• j – ten demographic groups (Male 0-14, Male 15-54, Male 55-64, Male 65-74, Male 75+, Female 0-14, …)

• t – year (2017, 2018, 2019, 2020)

We have observations on:

• Y_ijt – deaths for each town, demographic group, and year for January 1 through April 15th for each of the years 2017-2020

• POP_ijt – population for each town, demographic group, for the years 2017-2019 (the 2019 population is used for 2020)

• C_ijt – COVID treatment (indicator variable for pre-2020 vs 2020)

This is a mixed effects model, with grouping by town (i) and demographic group (j), with the data being “treatment” of COVID versus pre-2020: where the intercepts and slopes are specified as

β_0i = α₀₀ + α_0i Town pre-2020 effect (intercept)

β_1i = α₁₀ + α_1i Town 2020 effect (slope)

β_0j = γ₀₀ + γ_0j Demographic pre-2020 effect (intercept)

β_1j = γ₁₀ + γ_1j Demographic 2020 effect (slope)

I am going to treat all the demographic factors as fixed and (potentially) allow the town effects to be random. This seems reasonable because there are a large number of towns, each with a relatively small number of observations (eight demographic groups pre- and post-2020) but a small number of demographic groups each with a large number of observations (307 towns for each of pre- and post-2020). I can of course also examine this by treating all factors as fixed.) In this case the combined model is

4 Details on Estimation

Table 5 shows the coefficients for estimating Equation 1 as a mixed effect Poisson count model, with the town effects treated as random and the demographic effects as fixed. Using these estimates I can calculate predicted counts and thus predicted rates for the eight demographic groups both pre-COVID and for 2020. If we set the random effects at zero we will estimate excess mortality for the middle of the random effects distribution – i.e. the median town in terms of both pre-COVID and 2020 mortality. These results are shown in Table 8 below.⁸

We want the predicted excess mortality for the upper tail of the distribution – for towns that had high excess mortality and were thus, presumably, heavily infected. The estimated distribution of the 2020 random effect (α_1i) measures this right tail (more infected towns) and will provide a lower bound for age-specific IFR. We use the 99th percentile of the α_1i distribution. If the estimated 2020 demographic fixed effect for a group j is and the variance of the α_1i distribution is then the 99th percentile (2.5σ) of the 2020 mortality will be larger than the pre-2020 mortality by the factor . The excess mortality for group j will be

The coefficients in Table 5 show the age and sex profile of mortality in ratio terms (while the excess mortality in Table 1 show it as differences in rates – Equation 2). Mortality for demographic group j for the 99th-percentile town increased by the factor . These increases, as a ratio to the pre-2020 mortality, are shown in Table 6; they are large and larger for older men.

Male mortality is higher than female pre-COVID at every age, with the ratio given by exp (γ_{0F j} − γ_0Mj). The ratio exp (γ_{1F j} − γ_1Mj) measures how much the ratio increases in 2020.

Table 7 shows summary statistics and test for the random and fixed effects models. I estimate three nested models for each the random and fixed effect specification. The first includes no estimates for 2020 effects – setting α₁₀ + γ₁₀ + γ_1j and α_1i to zero. The second includes town-specific effects, estimating α_1i. For both random and fixed effects models a likelihood ratio test shows that town-specific 2020 effects are important. The third model includes 2020 age-specific excess mortality effects – α₁₀ + γ₁₀ + γ_1j. Again, a likelihood ratio test shows that these effects are highly significant.

The residual deviance in Table 7 provides a test for whether the Poisson regression adequately accounts for the observed variation – whether there is “overdispersion” relative to the Poisson assumption. The residual deviance is asymptotically χ²-distributed: a large value implies that the model does not account for the data very well. For both the random and the fixed effects model, the residual deviance is small relative to a χ² random variable – the probability of observing a value as large or larger than reported is 1.0.

For reference, Table 8 shows the estimated excess mortality for the median town, towns with presumably modest overall infection rate.

Calculating standard errors for the excess mortality is difficult. For a specific age group this would be the 99th percentile 2020 rate less the base-case pre-2020 rate shown in Equation 2. Standard errors for that expression is not easy because it is a non-linear expression of the coefficients. One might try doing so by the delta-method but that is complicated by two factors: first the standard error for is not easy to obtain, and second we want to aggregate age- and sex-specific rates to overall rates using regional population weights, leading to weighted averages of terms like Equation 2. The most feasible alternative is bootstrapping. Non-parametric bootstrapping where towns are sampled from the original distribution and Equation 1 re-estimated for each sample. Tables 1 and 8 show the results for bootstrapping with 400 samples (sampling with replacement from the town distribution).

Figure 3 shows a histogram of estimated 2020 excess mortality across the 612 towns. This shows how towns varied in their 2020 excess mortality, presumably resulting from different coronavirus infection rates.

• Download figure
• Open in new tab

Figure 3: Histogram of Predicted Excess Mortality by Town