Community factors and excess mortality in the COVID-19 pandemic in England, Italy and Sweden

Abstract Background Analyses of coronavirus disease 19 suggest specific risk factors make communities more or less vulnerable to pandemic-related deaths within countries. What is unclear is whether the characteristics affecting vulnerability of small communities within countries produce similar patterns of excess mortality across countries with different demographics and public health responses to the pandemic. Our aim is to quantify community-level variations in excess mortality within England, Italy and Sweden and identify how such spatial variability was driven by community-level characteristics. Methods We applied a two-stage Bayesian model to quantify inequalities in excess mortality in people aged 40 years and older at the community level in England, Italy and Sweden during the first year of the pandemic (March 2020–February 2021). We used community characteristics measuring deprivation, air pollution, living conditions, population density and movement of people as covariates to quantify their associations with excess mortality. Results We found just under half of communities in England (48.1%) and Italy (45.8%) had an excess mortality of over 300 per 100 000 males over the age of 40, while for Sweden that covered 23.1% of communities. We showed that deprivation is a strong predictor of excess mortality across the three countries, and communities with high levels of overcrowding were associated with higher excess mortality in England and Sweden. Conclusion These results highlight some international similarities in factors affecting mortality that will help policy makers target public health measures to increase resilience to the mortality impacts of this and future pandemics.


Introduction
T he coronavirus disease 19 (COVID- 19) pandemic has been spreading throughout Europe since early 2020 causing an estimated 580 000 deaths by the end of 2020. 1 The impact of the pandemic has varied widely across Europe resulting in speculation as to the possible reasons for higher mortality in some countries. The multitude of factors affecting mortality and the heterogeneity of data collected in different countries make accurate international comparisons challenging. Existing comparison studies have examined COVID-19 mortality between countries down to the NUTS3 regional level [2][3][4] but not for smaller areas that would allow the identification of the community characteristics which affect the risk of excess mortality in local populations. 5 Here, we study England, Italy and Sweden, three European countries differently affected by the pandemic and characterized by very different public health responses and policies. Italy was the first European country to be hit by large numbers of SARS-CoV-2 infections, soon followed by England and the rest of Europe. By 28 February 2021, after a year of COVID-19 in Europe, cumulative mortality attributed to COVID-19 in England, Italy and Sweden had reached 1930, 1618 and 1262 deaths per million people respectively. 6,7 Here, we aim to quantify community-level variations in excess mortality within the three countries and to identify how such spatial variability was driven by community-level characteristics. We analysed data on all-cause mortality at ages 40 years and over for 6791, 7900 and 5985 community-level areas in England (middle-layer super output areas, MSOAs), Italy (municipalities) and Sweden (Demografiska statistikområden, DeSOs) respectively, from 1 March 2020 to 28 February 2021.

Excess mortality
Excess mortality, typically calculated as a counterfactual comparison using data from years prior to the pandemic, is an overall measure of the mortality impact of COVID-19. It includes not just deaths known to be caused by SARS-CoV-2 infection, but also deaths where SARS-CoV-2 infection was undiagnosed and deaths resulting from behavioural, social and healthcare changes that accompanied interventions to reduce SARS-CoV-2 infections.

Study period and units of analysis
We considered mortality data over a 1-year study period (1 March 2020-28 February 2021) and a 5-year comparison period (1 March 2015-29 February 2020). The study period and comparison period were chosen to reflect the different dynamics of the pandemic in England, Italy and Sweden and to minimize the effect of the vaccination rollout which commenced on 8 December 2020 in England and 27 December in Italy and Sweden. [8][9][10] Due to the low numbers of deaths in people under the age of 40 years, 2,11 we limited the analysis to people over that age.
To investigate the association of community characteristics with excess mortality, we included data at small area level with a view to maintaining consistency across England, Italy and Sweden. For England, we used MSOAs (mean population: 8288, range: 2224-26 513). In Sweden, we used DeSOs (mean population: 1726, range 663-5291). 12 Both MSOAs and DeSOs have been defined based on population size. For Italy, the only comparable geographical units for which data on mortality and covariates were available was the municipality which has a much broader population range (mean: 7548, range: 30-2 808 293). 13

Characteristics of communities
We considered community characteristics measuring deprivation, air pollution, living conditions, population density and movement of people. These were used as covariates in our model to quantify their associations with excess mortality. A summary of community characteristics is presented in table 1. Covariates were provided as five categories (overcrowding and deprivation for Italy) or continuous data which we divided into quintiles, giving, in each category, $1358 MSOAs (England), $1580 municipalities (Italy) and $1197 DeSOs (Sweden).

Data sources
For England, mortality data were obtained from the UK Small Area Health Statistics Unit's (SAHSU) national mortality registrations supplied by the Office for National Statistics (ONS). We used mortality data for 2015-20 and provisional data for 2021 supplied by the ONS in July 2021 which included date of death and place of residence of the deceased. Annual population was from ONS mid-year (June) population estimates for MSOAs in England, 2015-19. For 2020, the population was a combination of ONS population estimates and linear regressions designed to minimize the effect of the deaths in the first wave of the pandemic in England (March-May 2020) (Supplementary section S1).
For Italy, mortality and population data were retrieved from the National Institute of Statistics (ISTAT) website. 25 Population data assessed on 31 December 2015 was used for the period 1 March 2015-28 February 2016, and so on.
For Sweden, Statistics Sweden supplied geocoded mortality and population data to each DeSO, based on the mortality registrations and residential address data in the Total Population Register. 18 Population data assessed on 31 December 2015 was used for the period 1 March 2015-29 February 2016, and so on.

Statistical methods
This study extends previous work by SAHSU examining excess mortality in England during the first wave of the COVID-19 pandemic. 5 We carried out all analyses for males and females and for each country separately. We split the data into four age groups: 40-59 years; 60-69 years; 70-79 years; and 80þ years.
Given the small population at MSOA and DeSO levels, as well for many of the Italian municipalities, we built a two-stage Bayesian spatial model, which overcomes potential issues of instability of mortality estimates to obtain estimates of excess death rates based on data for each community, its neighbouring communities and the national average. For the first stage model, we assumed that the number of deaths y itk for the ith area (i ¼ 1,. . .6791 for England; 1,. . .7900 for Italy; 1,. . .5985 for Sweden), the tth year (t ¼ 2015,. . .,2019) of the comparison period and kth age group (k ¼ 40-59, 60-69, 70-79, 80þ) is modelled as a Poisson distribution: where the log-transformed death rates are modelled as fixed and random effects: In (1), a 0 is the common intercept which represents the overall mortality rate, while b 0k is the age fixed effect for the kth age group. To account for the multiple years in the comparison group, c t provides a year-specific Gaussian random effect. Area-level effect b 0 are modelled as a combination of a conditional autoregressive structure (u 0i ), to capture local dependency and independent and identically distributed Gaussian random effects (v 0i ) to account for the similarities of all the areas in the study region (i.e. each country): The model has a precision parameter (s b ), while u represents the weight of the spatially structured component. 26 This model provides both local and global smoothing on the underlying death rate k ikt1 .
We obtained the posterior distributions of the death rates in each area, age group and year, k ikt1 ; we then averaged over the 5 years of the comparison period to obtain k ik1 , the expected death rate for the ith area and kth age group during the study period under the alternative scenario that the pandemic did not take place.
In the second-stage model, we estimated the ratio between death rates in the study period and the death rates we would have expected had there been no pandemic, using data for the comparison period. For the number of deaths in the ith area and kth age group in the study period, we specified the following model: Where the rate ratio, q ik represents the age-specific ratio between death rates in the study period and the comparison period k ik1 ð Þ. The product of the rates and the 2020 population (k ik1 .Pop i.2020.k ) can be considered the age-specific expected number of deaths for each spatial unit.
We modelled q ik in a similar way to stage one using terms to account for both space and age: Covariates were incorporated into this second stage log-linear model to evaluate their effect on the mortality rate ratio. For onevariable-as-a-time effects in (2) we added the term dX i where X i is the category of the variable in the ith area and d is the associated effect. Similarly for the full multivariable model evaluating the joint effect of all variables, we added P 6 n ¼ 1 . .6 covariates and j ¼ 2,. . .5 categories, as the first category is the reference one.
Uncertainty was considered in each stage by drawing 100 samples from the posterior distribution of each parameter, giving a total of 10 000 samples from which we defined 95% credible intervals (CI, 2.5th-97.5th percentiles). The number of samples ensures accuracy of the credible interval estimates, while limiting computational burden as discussed in Ref. [5]. We fitted the models using the Integrated Nested Laplace Approximation (INLA), through the R-INLA software package (http://www.r-inla.org/). 27,28 To show the change in death rates associated with the covariate categories, we report posterior mean and 95% CIs of mortality rate ratios. We obtained the posterior distribution of the estimated number of deaths across ages for the study period, calculated aŝ y i:2020 : ¼ P k k ik1 Â q ik Â Pop i:2020:k , and subtracted the corresponding number for 2020 using the rates for 2015-19, calculated aŝ y i:2015À19 : ¼ P k k ik1:2015À19 Â Pop i:2020:k where the samples for k ik1:2015À19 are drawn in equal proportions from each of the 5 years of the comparison period allowing the uncertainty to be fully represented. This difference was then divided by the 2020 population over 40 years old in that area and multiplied by 100 000 to give the excess deaths per 100 000 people, which we map along with the posterior probability of an excess above zero. We further calculated the % of variation explained by the covariates from the ratio of the variance of the fixed effect for the covariate (from all 10 000 samples) to the total variance.
As sensitivity analyses, we ran models: (i) without accounting for local spatial correlation; (ii) using continuous variables instead of quintiles for the covariates which were directly comparable across the three countries (NO 2 , PM 2.5 and population density).  Table 2 gives a summary of the distribution of excess deaths across the community-level areas. England and Italy showed positive excess mortality in more than 90% of areas, while for Sweden, values were lower at 71% and 55% for males and females, respectively. When considering >300 excess deaths per 100 000 population over 40, almost half of areas matched this criterion for males in England and Italy, but only 23% of areas did in Sweden.

Spatial patterns of excess mortality
The maps (figure 1) show that for England and Sweden, the high excess mortality was spread across the countries. In England, the lowest excesses are in rural areas with lower population densities. The maps suggest the largest excesses mortality in England were Where not specified as categorical data, the variables are continuous which we then divided into quintiles. Categories or quintiles were ordered based on the expected relationship between the covariate and excess mortality, with the highest category/quintile associated with higher excess mortality (e.g. high air pollution associated with higher excess deaths, but lowest change in workplace visits associated with higher excess deaths).
concentrated in the two largest cities, London and Birmingham especially for males. In Italy, the largest excesses in mortality were mostly in the north of the country in both urban and rural municipalities where the pandemic first hit and the impact was already significant before the lockdown in March 2020. Other areas of high excess include parts of Puglia and the islands of Sicily and Sardinia as well as large urban areas. In Sweden, the maps show little or no regional clustering of excess mortality. The patterns visible for England and Italy are more pronounced for males than females.

Associations of community characteristics with excess mortality
When considering the association between the covariates and excess mortality, there were striking similarities between the results for England and Sweden that were not always shared with Italy (figure 2). Considered individually (i.e. in univariate analysis), for England and Sweden, all covariates, except the workplace visits data showed positive associations with excess mortality across the categories for both sexes, although the association was generally weaker for females (Supplementary tables S13, S15 and S17). For Italy, only PM 2.5 and deprivation showed a positive association with excess mortality for both sexes. There were positive correlations between some of the covariates (Supplementary tables S10-S12); for example, Kendall's Tau for England was 0.56 and 0.61 between overcrowding and NO 2 and population density respectively, and for Sweden the equivalent correlations were 0.50 and 0.54. These contrast markedly with the same correlations for Italy (0.02 and 0.12). When the covariates were considered together (multivariate analysis), associations between excess deaths and NO 2 , and PM 2.5 mostly or completely disappeared for England and Sweden. The association with population density and overcrowding/living area was attenuated, but remained, especially for males. For Italy, the association with PM 2.5 and deprivation remained in the multivariable analysis.
Notably, the relationship with deprivation was strong under both univariate and multivariate analyses for all countries, both sexes. In particular, for Sweden there was a $15% (95% CI 10-20%) higher excess mortality for the most deprived category (compared to least deprived) for males and it was $13% (95% CI 8-19%) higher for females.
Aside from the positive association for deprivation, the results for Italy were distinct from the other two countries. Specifically, a weak negative association for both sexes with population density was present in the multivariable analysis for Italy whereas both England and Sweden showed a weak positive association; for both sexes, PM 2.5 showed a strong positive association for Italy. This was much weaker for England, and non-existent for Sweden.
For all countries, the covariates accounted for a higher percentage of the variation in excess mortality across the community-level areas for males than females (Supplementary tables S13-S18).

Results of sensitivity analyses
In sensitivity analysis, removing the random effect for local spatial correlation (Supplementary tables S26-S28), while altering the total estimated excess deaths very little (< 1% for all counties), indicated a positive association, particularly for males for both air pollution covariates in the multivariable model and there also seemed to be an association for workplace visits for England. The results for Italy revealed a stronger association with NO 2 , and negative associations for population density and workplace visits for both sexes. The results for Sweden showed no material differences to the main model for this analysis. The sensitivity analyses using continuous variables for NO 2 , PM 2.5 and population density (Supplementary tables S20-S25) showed no material differences to the main model, except a stronger association between PM 2.5 and excess mortality for males in England (Supplementary table S20).

Discussion
Our study provides insights into the comparison of spatial distribution of excess mortality and the role of community-level characteristics in three European countries with dissimilar social structures which applied different public health strategies. As illustrated in figure 1, our maps revealed strong geographical patterns in England and Italy, but less so in Sweden. It supports our choices of geographical areas across the three countries to identify such patterns.
The similar patterns of results between England and Sweden suggest that the characteristics of communities affecting excess mortality in the first year of the COVID-19 pandemic are shared between some countries, in particular deprivation, even if the overall rates of excess mortality varied between countries. The contrasting results from Italy could be explained in at least two ways. First, the way the disease first spread within the country, specifically its rapid spread in the north in the initial phases before the national lockdown and clear response measures were in place, and a consequent high excess in mortality as shown in figure 1 and reported previously. 29 This was not the case for England or Sweden where the disease spread later and more evenly. 30 Second, because spatial units in Italy are based solely on administrative area, resulting in more diverse population sizes, with large metropolitan areas (such as Rome or Milan) aside Our results suggest a stronger deprivation gradient in excess mortality for Sweden than the other countries. Although we used different definitions of deprivation for the three countries and our data did not provide any direct explanations, previous research has shown a positive correlation between income and compliance with nonpharmaceutical interventions, such as handwashing, social distancing, reducing mobility and remote work. 31,32 Sweden's reliance on recommendations rather than a strict lockdown 33 could explain the stronger deprivation gradient, although more research is needed to confirm this hypothesis.
Sweden is a smaller country (2019 population 10.3 million) than both England (2019 population 56.3 million) and Italy (2019 population 59.7 million) with lower excess mortality per 100 000 population and DeSOs are smaller than both MSOAs and Italian municipalities. These differences implied wider 95% CIs around the estimates for Sweden than around those for England and Italy, for both the area-specific excess mortality estimates and the relationships between excess mortality and the covariates. The lower February 2021 compared to the same period for the preceding 5 years. Posterior probability that excess deaths > 0 for males (c) and females (d). The posterior probability measures the extent to which an estimate of excess/fewer deaths is likely to be a true increase/decrease. Where the entire posterior distribution of estimated excess deaths for an area is greater than zero, there is a posterior probability of $1 of a true increase, and conversely where the entire posterior distribution is less than zero, there is a posterior probability of $0 of a true increase. This posterior probability would be $0.5 in an area in which an increase is statistically indistinguishable from a decrease. Contains OS data ß Crown copyright (2020). Data available under the UK Open Government Licence v3 percentage of variance explained by local clustering (Supplementary table S18) and lack of regional patterns of excess mortality (figure 1) for Sweden could be explained by greater homogeneity across communities in the country. This study not only provides insights into the community level factors affecting excess mortality, in particular the consistent association between deprivation and higher excess mortality shown for all three countries with different levels of overall deprivation, but also showed significant overlap in other important effects, for example, for both Sweden and England the multivariable results (figure 2) suggested that population density itself was a much less important driver of excess mortality than the one-variable-at-atime results at first suggested.

Strengths and limitations
Our analysis aimed to use similar community-level areas and a range of covariates previously linked to excess mortality in COVID-19 studies across three European countries which experienced different pandemic dynamics and implemented different time-varying, public health interventions. 34 The methodology applied combines a rigorous two-stage modelling approach, a large number areas, and expertise from each of the three countries.
The large range of population sizes for Italian municipalities is likely to affect the relationship between excess mortality and the covariates because the larger municipalities are unlikely to be homogenous communities with similar characteristics. The influence of the highly populous municipalities also affected the categories of the community characteristics. For example, 59.6% of the population over 40 years was in the top category of the population density covariate (Supplementary table S19).
Some community characteristics previously reported to show associations with COVID-19 related excess mortality 5 (proportion of population that is non-white, the number of care homes) could not be used here because the data (or close equivalents) were not available for all three countries. Our model did not include a temporal dimension, which made the inclusion of time varying Figure 2 The relationship between community covariates of small areas in England, Italy and Sweden and excess mortality from 1 March 2020 to 28 February 2021 compared to the same period for the preceding 5 years. Proportional increase in death rates shown as rate ratios (data are presented as posterior mean with 95% credible intervals) for categories of the distributions relative to lowest category. Multivariable relationship between covariates and excess mortality after adjustment for the other covariates, numerical values reported in Supplementary tables S14, S16 and S18 covariates known to affect mortality, such as temperature, 35 unnecessary. The inclusion of a year-specific random effect in our model should at least partially account for other confounders not directly included in our model. Although covariates were selected for consistency, it was not possible to use the same measures of deprivation or overcrowding for the three countries, meaning that they will be measuring different things. Furthermore, for the covariates which were the same or very similar (NO 2 , PM 2.5 , population density and workplace visits), the quintiles were defined per country, so cut-off points between quintiles were different for each country with Sweden having lower NO 2 and PM 2.5 quintile boundaries than England and Italy (Supplementary tables S4-S6) which should be considered when making direct comparisons between the results for the three countries for example in figure 2.
The process of estimating the counterfactual mortality rates for the community-level areas involves inherent uncertainty which might be reflected in some communities showing excess mortality even in the counterfactual case, so it was important that we accounted for all uncertainties in our model and reported results with 95% CIs. Consequently, our two-stage model fully propagated uncertainty from the first stage (comparison period) to the second stage (study period).
We used mortality data for the comparison period to estimate the expected number of deaths in the study period. Prior to the start of the study period, 29 COVID-19-related deaths had been recorded in Italy, and none in England and Sweden. 6 Additional factors beyond the scope of this study, like temperature and healthcare provision are likely to have had an impact on excess mortality, both between the countries and between areas within countries. In addition, we used data for England from the last national census in 2011 to obtain information on sociodemographic characteristics of communities. To the extent that there have been demographic changes in the decade since then, this may have led to misclassification of areas with respect to their community characteristics.
The study period was chosen partly to avoid the effect of the vaccination programmes on COVID-19 mortality. By the end of the study period, approximately 2.1%, 4.0% and 6.5% of the populations over 40 years of England, Italy and Sweden had been given two vaccine doses, but the time delay between full vaccination and reaching maximum immunity combined with the likely time delay between infection and death meant that the vaccine programme will likely not have had a measurable effect on excess mortality by 28 February 2021.
We did not account for the spread or extent of infection in our model because seroprevalence data were not available at the high spatial resolution that we used. It is reasonable to assume that over the course of the year of the study period, the disease had spread to all regions included in our study, although the higher excess mortality concentrated in the north of Italy strongly implies heterogeneity of infection rates in this country during the first wave of the pandemic.
Improving the comparability of community-level geographies and definitions of key covariates for epidemiological studies would greatly facilitate future international comparisons within Europe, allowing further research to build on our findings for future larger studies.

Conclusion
In what we believe to be the first international study examining the role of community-level factors on excess mortality during the COVID-19 pandemic, we have shown striking similarities between England and Sweden in spite of the different government responses to the pandemic in the two countries as well as the very different population sizes and demographics and behaviours. Although the pandemic in Italy was unique, both in the timing and localized spread in the initial phases, some common aspects remain such as the role of deprivation. This strengthens the message to address deprivation-related health inequalities before the next pandemic. In spite of the challenges inherent in making international comparisons with heterogeneous data, these results present a compelling reason to implement a cross country comparison in real time in future pandemics to help identify successful interventions. We believe the results presented here provide invaluable guidance for policy makers to target interventions such as localized public health messaging and provision of personal protective equipment that assist communities most at risk from mortality from this or future pandemics and that more interventions such as localized surge testing could have been implemented in high risk communities from an earlier stage.

Supplementary data
Supplementary data are available at EURPUB online.
personal information in Italy. The Italian mortality data, made available online by the national Institute of Statistics (ISTAT) in aggregated form, thus is GDPR compliant and therefore is exempt for the requirement of ethics approval.

Results
• Community level results (excess deaths, CIs, posterior probabilities) used in figure 1 are available on request to the authors.

Code availability
The computer code written in R 36 for the two stages of Bayesian models used in this work is available on request to the authors.