Abstract

Background

Several nonpharmaceutical interventions, such as masking, were mandated or recommended during the Coronavirus disease 2019 (COVID-19) pandemic. This study鈥檚 primary objective is to investigate the relationship between population-level mask usage and excess mortality across Europe.

Methods

We collected data on mask usage and other relevant variables from 24 European countries during 2020-2021, a period in which mask policies varied widely across nations, providing an ideal basis for a natural experiment. To assess the association between mask usage and relevant medical and socioeconomic data at the country level, we conducted both bivariate and multivariate regression analyses. Confounding factors were accounted for in the regression models, and numerous sensitivity tests were performed to ensure robustness.

Results

Statistically significant correlations were found between mask usage rate and age-adjusted excess mortality in both bivariate (Spearman coefficient = 0.477, p = 0.018) and multivariate (Standardized coefficient = 0.52, p = 0.0012) regressions. Likewise, vaccination rates showed negative and significant bivariate (Spearman coefficient = -0.659, p < .001) and multivariate (Standardized coefficient = -0.48, p = 0.0016) correlations with age-adjusted excess mortality.

Conclusions

No correlation was observed between mask usage rates and COVID-19 morbidity. However, significant associations were identified between mask usage rates, COVID-19 mortality, and excess deaths. Various hypotheses have been proposed to explain these associations, with thorough consideration given to potential confounders, such as socioeconomic factors and the severity of COVID-19 waves.

Background

The COVID-19 pandemic, which lasted nearly three years and officially ended in March 2023, provided an opportunity to assess the impact of nonpharmaceutical interventions (NPIs). Despite extensive research on the effectiveness of NPIs, findings have often been conflicting [1,2,3,4]. At the population level, the most definitive way to assess the impact of an NPI in the context of the pandemic is by evaluating its effect on COVID-19 mortality. However, official death counts remain a subject of ongoing revision by national health authorities and independent researchers [5,6,7].

Mask recommendations and mandates were among the most widely implemented yet highly debated measures during the COVID-19 pandemic. Research on the effectiveness of masking and mask mandates has yielded mixed conclusions. Randomized controlled trials (RCTs), widely regarded as the gold standard for evaluating drugs or therapies, have found little to no effect of masking on viral transmission [8,9,10]. Conversely, several observational studies 鈥� naturally more prone to confounding factors and biases, have claimed that mask-wearing is associated with reduced viral transmission [11,12,13,14]. Excess mortality has gained attention as a proxy for pandemic impact, with studies highlighting its advantages over reported COVID-19 deaths [7, 15, 16]. Nevertheless, no study on mask effectiveness has yet used excess mortality as an endpoint.

A previous analysis of European countries during the 2020鈥�2021 winter suggested that mask usage did not correlate with COVID-19 cases but showed a weak positive correlation with COVID-19 deaths. To further investigate this finding, the present study expands the analysis to a two-year period (2020鈥�2021) and incorporates multiple excess mortality estimation approaches [17,18,19,20,21,22], focusing on 24 European countries for which reliable excess death data are available. Using bivariate and multivariate statistical analyses, we systematically examine the relationship between mask usage and mortality, while carefully accounting for potential confounders through multiple statistical methods.

Methods

Study design

We examined the correlation between mask adherence and excess mortality in 24 European countries with populations over one million - Austria, Belgium, Croatia, Czechia, Denmark, Estonia, Finland, France, Germany, Greece, Hungary, Italy, Latvia, Lithuania, Netherlands, Norway, Poland, Portugal, Slovakia, Slovenia, Spain, Sweden, Switzerland and the United Kingdom, during the years 2020鈥�2021, a period during which different European countries had widely varying mask policies.

Data sources

All data was obtained from public sources. Data on mask usage was downloaded on 25th March 2023 from the Institute for Health Metrics and Evaluation (IHME) at the University of Washington [23]. The prevalence of mask wearing in a population is based on self-declaration represented by the 鈥減ercent of population reporting always wearing a mask when leaving home鈥� [24]. Data on cumulative percentage of excess mortality in 2020鈥�2021 (Per Levitt age adjusted (PLAA), Per Levitt not age adjusted (PLNAA), eLife, Economist, Lancet and WHO) were taken from Levitt et al. [17]. Multiverse excess mortality was from Levitt et al. [18]. Except for the Gini index, obesity rate and population density all data shown in Table听1 and Table S1 (Supplementary file 1) were downloaded from the Our world in data (OWID) website [25] on 31th March 2023. The Gini index data was downloaded from the WorldBank Open Data [26] on 7th May 2023; the obesity rate was downloaded from the Central Intelligence Agency (CIA) Fact Book [27] on 25th April 2023 and the population urban density was from the Global Human Settlement Layer [28], downloaded on 18th May 2023. The entire dataset used in this study can be found in the fulldataset.csv file deposited in the Zenodo repository [29].

Table 1 Correlation between excess deaths (PLAA) and relevant variables

Statistical analyses

All statistical analyses - bivariate and multivariate regressions, were performed in R (version 4.1.0 (2021-05-18)) or in the R interface software JASP [30]. All R codes are available in the AllCode.R file in the Zenodo database [29]. More specifically, for the bivariate regressions, Spearman鈥檚 coefficients were obtained for all variables. Shapiro-Wilk analyses showed that most variables were not normally distributed. Scatter plots of each variable against % mask usage (shown in Supplementary file 1, Figure S2) used the lm function of R to plot the regression line. In the multivariate regression analyses, 鈥渇itted values鈥� versus residuals and 鈥渇itted values鈥� versus standardized residuals plots were used to check for the presence of nonlinear effects. The basis for the selection of the independent variables in our main regression model was a statistically significant bivariate correlation with PLAA at the 0.05-significance level. The variables selected for the model were mask usage, vaccination rate (people fully vaccinated/hundred), HDI and CardLife. A scale location plot and a studentized Breusch-Pagan test were used to check for the presence of heterocedasticity. Normality of residuals was checked by computing skewness and excess kurtosis (with confidence intervals) by a Shapiro-Wilk test and a QQ plot. All data was evaluated with plots obtained using simulated data in which linear model assumptions were satisfied exactly. Presence of multicollinearity was checked by regressing each independent variable against all others and computing the corresponding R2. Cook鈥檚 distance was used to detect possible overly influential data points. For the sensitivity analyses of our main regression, we used both a bootstrap and Huber-White heterocedastic-consistent robust sandwich estimators for standard errors (HC3), which works well with violations of heterocedasticity even in small samples (see Long and Ervin [31]). Confidence intervals for standardized regression coefficients were computed using the new betaDelta and betaSandwich R packages [32]. A total of 63 alternative regressions with different choices of covariates were performed to investigate whether inclusion of other potential confounding variables could change our main conclusions. To determine the maximum number of covariates per regression we conducted a power estimate and established a minimum power of 80%.

Results

The inclusion criteria for this study was any European country with a population larger than one million and for whom reliable excess death data was available [17], a total of 24 countries: Austria, Belgium, Croatia, Czechia, Denmark, Estonia, Finland, France, Germany, Greece, Hungary, Italy, Latvia, Lithuania, Netherlands, Norway, Poland, Portugal, Slovakia, Slovenia, Spain, Sweden, Switzerland and the United Kingdom. Many countries adopted mask mandates or mask recommendations during the COVID-19 pandemic. Mask usage rate, defined as the 鈥減ercent of the population reporting always wearing a mask when leaving home鈥� [23], varied country-wise in the two-year period analyzed in this study (2020鈥�2021). While Italy, Spain and Portugal implemented strict mask mandates already in the spring of 2020, Norway and the Netherlands intervened only in the winter of 2020鈥�2021, and a third group of countries, notably Denmark and Sweden displayed very low levels of mask usage during the entire two-year period (Figure S1 鈥� Supplementary file 1).

Bivariate correlations

Our first step was to analyze the associations between excess mortality, medical and socioeconomic variables (Table听1). The preferred method for assessing excess mortality was the 鈥橮er Levitt Age Adjusted鈥� (PLAA) calculation, as described by Levitt et al. [17]. This approach stood out as the sole excess deaths calculation technique that factored in the age structure of a country鈥檚 population. Since COVID-19 has a clear age-related risk mortality an approach that takes age stratification into account is likely to express more faithfully the rate of excess deaths in a population [17, 33]. Notwithstanding, all excess mortality approaches (Levitt Age-Adjusted [17], Levitt Not Age-Adjusted [17], Economist [22], eLife [20], Lancet [21], Multiverse [18] and WHO [19]) strongly correlated with each other with correlation coefficients around 0.9 or higher and p < 0.001 (Table S1 鈥� Supplementary file 1).

Table听1 shows that the PLAA excess death index strongly correlated with COVID-19 deaths (per million) ($\uprho$ = 0.878), indicating that an important share of 2020鈥�2021 excess deaths was due to COVID-19 mortality. The cardiovascular death rate (CVD) also correlated with PLAA ($\uprho$ = 0.667). Health and socio-economic parameters, such as Vaccination rate, life expectancy at birth (in 2019), GDP per capita and the Human development index, but not the Gini coefficient, that measures income inequality, were negatively associated with PLAA. Mask usage rate showed a moderate positive association with PLAA ($\uprho$ = 0.477).

Several other variables were considered for their potential relevance to excess mortality. These included obesity rates, urban population density, and the stringency index of COVID-19 containment measures (a composite measure based on 9 response indicators: school closures; workplace closures; cancellation of public events; restrictions on public gatherings; closures of public transport; stay-at-home requirements; public information campaigns; restrictions on internal movements; and international travel controls (). However, upon analysis, none of these factors demonstrated a statistically significant association with excess mortality rates. Additionally, our analysis revealed no statistically significant correlation between a country鈥檚 proportion of seniors (individuals aged over 65 years) and its excess mortality rates (Table听1). Similarly, no significant correlation was found between this demographic factor and COVID-19 deaths, as shown in Table S2 of Supplementary file 1. This is true for PLAA, that is age-adjusted, but also for all other excess death metrics that are not age-adjusted (Table S4 鈥� Supplementary file 1). A plausible explanation for this finding is that European populations exhibit similar levels of seniority, making it unlikely that small variations in age structure would lead to significant differences in mortality rates across these countries.

It has been previously shown that the average rate of mask usage in Europe during the 6-month period encompassing the winter of 2020鈥�2021 was not associated with the rate of COVID-19 cases, but was positively correlated with COVID-19 mortality [34]. Our current analysis, spanning a two-year period, indicates that mask usage exhibited no significant association with COVID-19 case rates ($\uprho$ = 鈭�0.011). This finding suggests that, when examined at a population level, the widespread use of masks did not appear to have a substantial impact on the transmission of COVID-19 (Table听2). Conversely, the rate of mask usage in Europe showed a moderate positive correlation with PLAA ($\uprho$ = 0.477) (Fig.听1) as well as with all other excess death metrics (PLNAA, eLife, Economist, Lancet, WHO and Multiverse) (Table听2). Mask usage also showed a significant correlation with COVID-19 mortality ($\uprho$ = 0.415).

Table 2 Bivariate correlations between mask usage (predictor/explanatory variable) and potential response variables

Fig. 1

Scatterplot of mask usage and excess mortality. The average percentage of mask usage in the two-year period 2020鈥�2021 in each of 24 European countries was confronted with the cumulative percentage of excess death (PLAA) during the same period. $\uprho$, Spearman鈥檚 coefficient correlation

Of the other potentially relevant response variables masks significantly correlated with the Stringency Index ($\uprho$ = 0.518), suggesting that countries that masked more also mandated stricter NPIs. A correlation between mask usage rate and socioeconomic factors such as 鈥楪DP per capita鈥�, HDI and the Gini coefficient was also observed. Mask usage was not significantly associated with vaccination rate, age (% of >65), cardiovascular death rate, population (urban) density, life expectancy and obesity rate. Some of these variables are associated with COVID-19 mortality and/or excess deaths (see Table听1 and Table S2 鈥� Supplementary file 1). Scatterplots of the bivariate correlations between mask usage and each response variable are shown in the supplement (Figure S2 鈥� Supplementary file 1). Overall, the most notable finding from the bivariate analyses is that mask usage was positively associated with excess mortality and showed no correlation with COVID-19 case rates.

Multivariate regression

Multivariate linear regression was used to estimate the independent effect of each variable on the cumulative percentage of excess deaths (PLAA). We initially considered the following variables as candidates for independent variables in our main regression model: mask usage, vaccination rate (people fully vaccinated/hundred), HDI, 鈥橤DP per capita鈥�, stringency, total tests per thousand inhabitants, population (urban) density, 鈥機VD rate鈥�, life expectancy, obesity rate, diabetes prevalence and Gini index. It is important to note that since PLAA is age-adjusted, variables reflecting the age pyramid can be omitted from the regression (although, as a matter of caution, they were considered in the sensitivity analyses). Since we only have 24 datapoints, inclusion of too many independent variables should be avoided, otherwise the tests for the regression coefficients will be underpowered. For this reason and in order to avoid multicolinearity, we performed a principal component analysis to aggregate highly (Pearson)-correlated pairs of variables, as follows: (I) cardiovascular death rate and life expectancy (r = 鈭�0.96) were aggregated into a new variable denoted CardLife. We arbitrarily chose the sign of CardLife to make it positively correlated with cardiovascular death rate. (II) HDI was used as a proxy for 鈥橤DP per capita鈥�. In addition to the fact that these two variables showed a strong correlation (r = 0.89), 鈥橤DP per capita鈥� is very similar to another index, 鈥橤NI per capita鈥�, which has been originally incorporated in the HDI. The CVD rate and HDI were not combined, despite the fact that this pair shows the third-highest correlation (r = 鈭�0.69), indicating a moderate relationship between the two.

The results of the multivariate regression are shown in Table听3. Vaccination was associated with a decrease in PLAA (p = 0.0016) and mask usage was associated with an increase in PLAA (p = 0.0012). HDI and CardLife were not significant at the 0.05 level (p = 0.42 and p = 0.21, respectively), but the signs of the coefficients align with what would be expected, i.e., HDI has a negative coefficient (suggesting that larger HDI is associated with less deaths) and Cardlife has a positive coefficient (suggesting that high 鈥機VD rate鈥� and lower life expectancy are associated with more deaths). To facilitate the analysis and to enable a meaningful comparison of coefficient sizes, Table听3 displays standardized regression coefficients expressed in standard deviation units for both independent and dependent variables. Figure听2 provides the point estimates and 95% confidence intervals for the standardized coefficients of the multivariate linear regression of PLAA on the other variables.

Table 3 Multivariate regression between independent variables

Fig. 2

Point estimates and 95% confidence intervals for the standardized coefficients of the multivariate linear regression of PLAA on vaccination, HDI, Cardlife and mask usage. Positive values for the coefficients indicate that an increase of the value of the given variable is associated with an increase in age-adjusted excess mortality

It is important to take into consideration that variables not included in our main model due to lack of statistical significance may still have an influence on PLAA. To account for this possibility we performed several sensitivity analyses with additional variables (see Section 4 in Supplementary file 2). We established a maximum of 7 independent variables per regression in order to satisfy a power estimate of at least 80% (Section 3 in Supplementary file 2). We also performed standard diagnostic tests to evaluate the main linear multivariate regression (Section 1 in Supplementary file 2) and sensitivity analyses in which linear model assumptions were relaxed (Section 1 in Supplementary file 2). In all sensitivity analyses considered, vaccination was statistically significantly associated with a decrease in PLAA and mask usage was statistically significantly associated with an increase in PLAA.

Discussion

This is the first study to examine the impact of mask usage on excess mortality across European countries, encompassing a population of about 500 million people. By applying bivariate and multivariate regression analyses we retrospectively examined the impact of population-level mask usage on excess mortality across Europe. The main conclusions of this study are twofold: at the population level (1) masks did not reduce COVID-19 transmission, and (2) mask usage is significantly associated with excess mortality.

Why Europe? The European continent is particularly well-suited for this analysis due to several key factors. Firstly, it comprises numerous densely populated countries within a relatively compact geographic area, each exhibiting distinct approaches to the COVID-19 pandemic, especially regarding mask mandates and recommendations. Additionally, European democracies offer reliable data on various public health topics, including mortality rates and health-related metrics. Furthermore, European populations share similar age distributions and generally have good access to healthcare services.

On the effectiveness of masks in reducing viral transmission, our retrospective observational study aligns with the comprehensive and authoritative systematic Cochrane review and meta-analysis by Jefferson et al. [8], that concluded that high-quality randomized controlled trials did not demonstrate a clear reduction in respiratory viral infection with the use of surgical masks. Accordingly, it has recently been shown an inverse correlation between the quality of studies included in meta-analyses on masks and the reported effect size of masks on viral transmission [35].

Bivariate regression analyses showed that all seven excess mortality approaches (Economist, eLife, Lancet, multiverse, PLAA, PLNAA and WHO) displayed moderate positive correlations with mask usage rate and were also associated with COVID-19 death rate. The multivariate regression analyses and sensitivity tests confirmed the association between masks and excess mortality.

What could explain the association between masks and excess deaths? One hypothesis is that mask use may lead to increased mortality in COVID-19 patients by promoting the re-inhalation of virions, potentially worsening patient prognosis [36]. However, there is to date no experimental data that confirms this mechanism. Additionally, several prospective interventional studies reported on mask adverse effects (see, for instance, [37,38,39,40]), but none of them was as severe as to cause death. However, the prolonged use of masks may adversely affect the health of particularly vulnerable individuals - an effect that could not have been clinically tested, potentially contributing to their death.听It is important to emphasize that, while a retro- spective observational study such as this cannot establish causality, it may contribute to questioning the effective- ness and safety of widespread mask use in the population.

Another possibility is that mask usage is actually a proxy for other NPIs, such as lockdowns and mobility restrictions or for socio-economic discrepancies. However, these variables were taken into account in the multivariate regression, so the association between mask usage and excess mortality cannot be explained by them.听The bivariate regressions showed a听negative correlation听between socioeconomic factors (GDP per capita, HDI) and excess mortality.听Why are rich European countries less likely to experience excess deaths? Leaving aside the tautological aspects of this question, one possibility is that richer societies generally enjoy better healthcare systems, increased access to medical resources, and healthier lifestyles. Although this may be true on a global scale, where first-world nations typically offer better healthcare than developing countries, all 24 nations examined in this study are categorized as high-income, meaning their economic disparities are relatively small. Furthermore, the ultimate endpoint for determining the health of a population is the life expectancy at birth. In the Nordic countries of Denmark, Norway, and Sweden (DNS), where masking rate was, on average, 7.7 times lower than in Italy, Portugal, and Spain (IPE), the average life expectancy (in 2019) was 82 years, while in the IPE countries the average life expectancy was 83 years, suggesting that the IPE countries are at least as healthy as the richer DNS countries. All of this suggests that economic disparities cannot explain the low excess mortality of the DNS countries in the pandemic years.

Another potential confounder is the proportion of seniors in a population. In this study, however, the excess death approach of choice, PLAA, had already been adjusted for age. Furthermore, the percentage of seniors did not correlate with any excess mortality metric (Table S4 鈥� Supplementary file 1) or COVID-19 mortality (Table S2 鈥� Supplementary file 1). This lack of correlation can be attributed to the fact that European countries have similar patterns of age stratification. Consequently, it is unlikely that small differences in the proportion of seniors between these countries could explain the differences in mortality rates.

The multivariate regression analysis revealed a statistically significant positive association between excess deaths and masks, as well as a statistically significant negative association between excess deaths and vaccination rates. These results held consistently across all sensitivity analyses, including those that relaxed linear model assumptions and those that incorporated additional variables into the regression. Standard diagnostic tests did not reveal any evidence of violations of linear model assumptions. The scatter plot of residuals in particular does not indicate any influence from nonlinear functions of the independent variables.

It could be argued that high-masking countries had a relatively high excess mortality because they experienced stronger COVID-19 waves than low masking countries, and that, therefore, their high mask compliance reflects their epidemiological circumstances. However, Figure S1 (in Supplementary file 1) shows that there were almost no cases in which masking rates followed an increase in deaths, which would characterizes a clear response to the death peaks. Furthermore, mask mandates and high mask compliance were already in place before the strong second COVID-19 wave (winter of 2020鈥�2021), yet countries with high mask compliance did not fare better than those with low masking rates [34]. This indicates that masks did not come as a response to pandemic peaks, but rather as a preventive measure.

Our analysis ended on 31^st December 2021. We did not continue through 2022 because by then many mask mandates have already been abolished and the general rate of mask usage dropped considerably (Table S3 鈥� Supplementary file 1). On the other hand, excess mortality trends听kept improving for the low masking countries (DNS) during 2022. According to 鈥淭he Economist鈥� [22], Sweden, with an average masking rate of just 3.3% (Fig.听1), was ranked 20^th out of 24 in cumulative excess mortality by the end of 2021. By the end of 2022, it had dropped to the 22^nd position. Other excess mortality metrics hailed Sweden as having the lowest excess mortality country of all Europe [41] or of all OECD countries [42].

A primary limitation of this study is its observational retrospective design, which, as mentioned earlier, cannot infer causality. Although efforts were made, particularly in the multivariate regression analysis, to rule out potential confounding factors, these cannot be totally eliminated in a non-randomized study. Thus, unknown confounders may be present and were not taken into account. Additionally, we cannot offer a definitive explanation for the correlation between mask usage and excess mortality; at this stage we can only propose hypotheses and present some arguments regarding their feasibility. Other limitations were that no in-country comparisons with subsets of the populations have been performed and that the linear regression model that we used did not include non-linear terms and interaction terms between the variables. In addition, although the sensitivity analyses and the scatter plot of residuals in particular does not indicate any influence from nonlinear functions of the independent variables this possibility cannot be completely excluded due to the small number of data points (24 countries) in the analysis.

Conclusions

In conclusion, our results suggest that, at a population level, masks not only fell short of preventing COVID-19 transmission in Europe but may have also contributed to unforeseen adverse effects. Nonetheless, there remains the possibility that an unmeasured variable is driving the observed correlation between mask use and excess mortality 鈥� a limitation inherent to all observational/ecological studies and a key reason this type of study is unable to establish causality. Still, this study can and should instigate further inquiry into the effectiveness and safety of masks. High-quality randomized trials would be invaluable to conclusively determine whether masks effectively prevent respiratory viral transmission, have no significant impact, or potentially cause adverse effects.