COVID-19 attributed mortality and ambient temperature: a global ecological study using a two-stage regression model.

A negative correlation between ambient temperature and COVID-19 mortality has been observed. However, the World Meteorological Organization (WMO) has reinforced the importance of government interventions and warned countries against relaxing control measures due to warmer temperatures. Further understanding of this relationship is needed to help plan vaccination campaigns opportunely. Using a two-stage regression model, we conducted cross-sectional and longitudinal analyses to evaluate the association between monthly ambient temperature lagged by one month with the COVID-19 number of deaths and the probability of high-level of COVID-19 mortality in 150 countries during time t = 60, 90, and 120 days since the onset. First, we computed a log-linear regression to predict the pre-COVID-19 respiratory disease mortality to homogenize the baseline disease burden within countries. Second, we employed negative binomial and logistic regressions to analyze the linkage between the ambient temperature and our outcomes, adjusting by pre-COVID-19 respiratory disease mortality rate, among other factors. The increase of one Celsius degree in ambient temperature decreases the incidence of COVID-19 deaths (IRR = 0.93; SE: 0.026, p-value<0.001) and the probability of high-level COVID-19 mortality (OR = 0.96; SE: 0.019; p-value<0.001) over time. High-income countries from the northern hemisphere had lower temperatures and were most affected by pre-COVID respiratory disease mortality and COVID-19 mortality. This study provides a global perspective corroborating the negative association between COVID-19 mortality and ambient temperature. Our longitudinal findings support the statement made by the WMO. Effective, opportune, and sustained reaction from countries can help capitalize on higher temperatures' protective role including the timely rollout of vaccination campaigns.

Introduction

The COVID-19 pandemic caused by the rapid global spread of novel coronavirus (SARS-CoV-2) continues to harm population health at an unprecedented rate [1]. Early and effective containment measures significantly reduce the virus spread, protect the public, and prevent capacity-constrained healthcare systems from becoming overloaded [2-4]. Yet, with no effective treatment or vaccine available at the onset of the pandemic, many countries resorted to different public health measures to reduce the disease’s spread, including closing schools and public places to complete lockdowns [5]. Notwithstanding, several countries have experienced a higher disease burden with disproportionate numbers of cases and deaths without further consensus on the effect of temperature and seasonality on transmission and consecutive mortality [6-9].

Some studies have found a negative relationship between temperature with both -COVID-19 infectivity [10-16], -and the risk of death due to COVID-19 [15-19], while a cross-country analysis using early released data found a modest relationship between average temperature and COVID-19 reproduction rate [20]. Due to the limited evidence currently available, the World Meteorological Organization (WMO) has reiterated the importance of government interventions and warned against relaxing control measures because of higher temperatures [21]. Furthermore, without sufficient knowledge of infection, community engagement, public health system capacity and adequate border control measures in place [22], many countries relaxed their measures going into the warmer periods in 2020 despite warnings from experts [6,22]. Nowadays, with several countries with ongoing vaccination campaigns, only a few countries have entirely controlled the spread of infections and disease severity [23]. Further understanding of temperature’s role in countries with a high mortality attributable to COVID-19 may help plan vaccination campaigns, especially before cold seasons begin.

Cross-country comparisons using COVID-19 mortality rates can be challenging and there is a notable lack of global studies aiming to address these factors using a broader and comparative perspective [18,24-26]. First, countries have employed different testing methods, standards to report numbers, diagnostics definitions, and criteria for COVID-19 related death declaration [18,26]. Second, the cases onset dates are dissimilar between countries [7,9,26-28]. Third, countries have different individual compositions, epidemiological profiles, and public health resources that may determine the severity of COVID-19 [24,27]. Despite these shortcomings, a study carried out by Sornette et al. (2020) analyzed mortality by classifying countries according to their geography. Apart from government measures taken and demographic and cultural factors, they found that climatic features such as temperature may explain the variation in mortality rates, especially in western countries [19].

So far, all the studies have focused on the association between temperature and cumulative deaths or cumulative mortality [15,19,29]; however, little is known on the role of temperature in the probability of countries facing high levels of COVID-19 mortality. Identifying potential factors for the likelihood of COVID-19 attributable mortality may help understand the characteristics of those countries that have been more at risk and had more burden over time.

To address the knowledge gap in the current literature on the cross-country association between ambient temperature and COVID-19 attributable mortality, we developed the present study from a global ecological perspective using a two-stage modeling approach to balance countries’ differences in resources and compositions.

Methods

Study design and sample

We employed cross-sectional and longitudinal analyses to explore the association between ambient temperature and COVID-19 number of deaths; and the link between ambient temperature and the likelihood of high-level of COVID-19 mortality in a sample of 150 countries.

We harmonized data through a deterministic data linkage process by combining records from different sources with the same three-letter ISO 3166–1 code. We extracted data from four sources. For further information, see supplementary materials (section A).

We initially included 152 countries in the analysis for those with recorded available data for a minimum of 90 days since the first confirmed case (missing rate = 17.6%; N = 152). Countries with complete information on the independent variables were therefore kept in the analyses, resulting in an analytical sample of 150 countries (see sample definition protocol in Figure A1, supplementary material). We also employed Cook’s distance test to analyze the most influential data points within the sample (see supplementary materials, section D) [30].

Dependent variable

The primary outcomes were COVID-19 number of deaths and the probability of high-level of COVID-19 mortality obtained from the Jhon Hopkins University (JHU) data repository (date of access: October 31st, 2020). We counted deaths attributable to COVID-19 at different points in time since the first reported case was confirmed in each country. First, we computed it at the 60th day and continued calculating period mortality using continuous 30-day intervals (i.e. t = 60, 90 and 120). We calculated mortality rate by the time after the first case was confirmed to make the countries comparable while accounting for the time-lapse in COVID-19 number of deaths [31]. Information on the number of daily reported cases is publicly available at Our World Data [7]. Afterward, we defined mortality by dividing the number of deaths by the 2019 population reported for each country. Also, we categorized COVID-19 mortality rate into two groups – ‘low to moderate level of COVID-19 mortality’ (coded as 0) and ‘high-level of COVID-19 mortality’ (coded as 1). We used different cutoff points for categorization, starting from the median in each timepoint (t = 60, 90 and 120) until the median plus two standard deviations. We analyzed the distribution of the unadjusted and adjusted ambient temperature Odds Ratios (ORs) when using the different definitions of a high-level of COVID-19 mortality (HLCM). We finally selected the cutoff fulfilling the statistical convention of at least 10 events-per-variable (EPV) at each time point. Additionally, we examined a longitudinal model using at least 2 EPV [32,33] (see more details at supplementary materials, section H).

Independent variables

1.- The monthly temperature lag: we calculated a monthly ambient temperature lag expressed in Celsius degrees (°C) using one month lagged temperature (tcurrent – 30 days) according to t = 60, 90, and 120 days since the first confirmed case and by country (e.g. Afghanistan reached t = 60 in June; therefore, we used the average temperature of May). Temperature data from November 2019 to October 2020 were extracted from the available data on ERA5 and analyzed in Copernicus. The data was cropped by country using GADM version 2.8 shapefiles [34] (date of access: 1 December 2020).

2. The probability of mortality due to respiratory diseases: Countries were balanced using the pre-COVID-19 estimated probability of respiratory disease mortality. The variable was constructed using the number of deaths attributed to respiratory diseases reported in 2017 from the Global Burden of Disease Study (GBD) and country’s local population size [35].

3. The stringency of government measures in response to the COVID-19 outbreak: The Oxford COVID-19 Government Response Tracker (OxCGRT) [36] stores data on eleven indicators representing the stringency level in the government response against COVID-19. Data were recorded daily by each country (date of access: October 31st, 2020). The OxCGRT index ranges from 0 (no government stringency) to 100 (very strict government). We set the government stringency index at t = 30, 60, 90, and 120 since the first case was reported and used the variation between periods (e.g. t = 90, then government measure = Δt90, t60). For further details, see supplementary materials section A.

4. The countries’ hemisphere: Countries’ geographical location wasclassified as Northern or Southern hemisphere based on their winter and summer seasons. Countries facing winter between December and March were categorized as ‘Northern’ countries and countriesfacing summer as ‘Southern’ countries. Rainy or dry seasons were not considered.

Auxiliary variables

We used seven variables for the first stage of our analysis (explained in the statistical analysis section). These variables included the percentage of women and the disability-adjusted life years (DALYs) attributed to asthma extracted from the GBD; the prevalence of obesity, and the level of air pollution (pm 2.5) obtained from the World Health Organization (WHO); and the human development index (HDI), population density, and the proportion of people aged 65 and older extracted from the World Bank (WB). See further details on the sources in the supplementary materials section A.

Statistical analysis

Firstly, we used heatmaps to describe the cross-country variation in temperature and the crude COVID-19 attributed mortality (logged and population-adjusted). Secondly, to avoid potential biases in the analysis, we employed a two-step regression model to study the link between ambient temperature and COVID-19 number of deaths, and the relationship between ambient temperature and the probability of HLCM. This method is based on Heckman’s approach [37-39], which has been widely used in previous studies to correct nonrandomly selected observations and to avoid potential biases in the analysis (i.e. available countries in this study).

In the first stage, we balanced our sample by estimating the pre-COVID-19 probability of respiratory diseases mortality to have a homogeneous sample given countries’ baseline epidemiological characteristics, avoiding multicollinearity while maintaining relevant factors previously seen as risk factors toward COVID-19. This step permits us to correct our second-stage models to account for specification errors. In the first stage, a log-linear transformation to compute the respiratory disease mortality using robust standard errors presented the best goodness-of-fit according to the R2 and Akaike Information Criterion (AIC) (see supplementary materials, section B, C, and D). We explored different models by adding characteristics related to respiratory disease mortality from the GBD study [35]. Other variables related to respiratory disease mortality were tested and discarded as predictors due to the lack of fit and multicollinearity issues. Finally, we predicted the mortality rate attributed to pre-COVID-19 respiratory diseases as a result of the first-stage process (Equation 1).

This model uses a fixed time point.

In the second stage, we corrected the estimates for selection-bias by adding the predicted probabilities from the preceding step (i.e. pre-COVID respiratory mortality) as an additional independent variable. Therefore, based on Equation 2, we ran cross-sectional and longitudinal negative binomial regression models for the period incidence risk of COVID-19 deaths and logistic regression models for the likelihood of HLCM.

The timestep was fixed to the specific days (t = 60, 90, and 120) for cross-sectional models. Longitudinal models included the same covariates throughout the three different time points. No collinearity was present amongst our predictors and the dependent variables (see supplementary materials, section G). We used Bootstraps errors with 1,000 iterations to account for sampling biases.

‘t’ stands for time = 60, 90, and 120 days since the onset

Y refers to ‘number of deaths’ in negative binomial regression and ‘high-level of COVID-19 mortality’ in the logistic regression analyses, respectively. Cross-section and panel data models were used. ∆ stands for variation between period t and t-1. ‘Region’ stands for countries’ hemisphere.

All analyses were performed using STATA 16.1 [40], QGIS 3.6 (QGIS Geographic Information System) [41], and R software 4.0.2 [42]. An online repository for data management and consolidation is available at https://bit.ly/36IKhhJ and https://bit.ly/2UXAz8B for data visualization and examination.

Results

Table 1 summarizes the descriptive characteristics of our sample. The COVID-19 number of deaths drastically increased through the time points. COVID-19 mortality increased over time, but it decelerated between t = 90 and t = 120. Specifically, it increased twofold in t = 90 compared to t = 60 (Meant:90 = 2.56; 95% CI: 1.37–3.74) and decreased at t = 120 compared to t = 90 (Meant=120 = 2.05; 95% CI: 0.89–2.13). Ambient temperature increased by 1.5°C per defined time (Meant:60 = 20.83; 95% CI: 17.13–24.54; Meant:90 = 22.33; 95% CI: 18.69–25.97; Meant:120 = 23.78; 95% CI: 20.16–27.39), whereas the index of government measures increased drastically between t = 0 and t = 60; however, it shrank after that period.

Table 1.

Descriptive statistics of the sample (N = 150)

Country-level characteristics	MEAN	SD	IQR
First-stage variables
Women (%)	49.98	2.91	1.12
People aged 65 and older (%)	5.65	4.42	5.40
Obesity (%)	18.10	9.47	17.9
DALYs Asthma (standardized)	0.21	0.16	0.17
Low HDIa	53.24	7.57	13
Medium HDIa	74.18	3.99	5.9
High HDIa	88.34	4.36	7.6
Population density (population/km2)	136.78	216.48	109.20
Low air pollutionb	11.72	3.36	5.34
Medium air pollutionb	21.99	3.26	5.63
High air pollutionb	50.16	17.83	19.75
Respiratory disease mortality	0.03	0.02	0.023
Second-stage variables
COVID-19 Deaths at t = 60	366.5467	1310.381	144
COVID-19 Deaths at t = 90	1279.553	4967.145	235
COVID-19 Deaths at t = 120	1214.267	5293.42	346
COVID-19 Mortality at t = 60	1.1.51	3.86	0.85
COVID-19 Mortality at t = 90	2.56	7.37	1.16
COVID-19 Mortality at t = 120	2.05	4.13	1.65
Ambient T°C at t = 30	17.58	10.65	17.15
Ambient T°C at t = 60	19.06	9.62	14.83
Ambient T°C at t = 90	20.57	8.43	12.90
Δ Government measures t = 60/t = 30	7.77	2.18	12.04
Δ Government measures t = 90/t = 60	−6.99	1.07	12.96
Δ Government measures t = 120/t = 90	−8.76	0.99	14.35

Notes: stands for variation between two periods. SD is standard deviation, while IQR is for the Interquartile range. a Countries level of Human Development Index [HDI] was divided using terciles. b Countries level of air pollution were classified using terciles.

Figure 1 shows the average ambient temperature by country, whereas Figure 2 depicts the death rate (adjusted to the population size per 100,000 habitants) since the onset. The Figures indicate that countries close to the Equator and the southern hemisphere had the highest temperatures but low or medium levels in deaths on average (e.g. Australia, Democratic Republic of the Congo, Ghana, Singapore). On the contrary, northerly countries faced the highest number of deaths attributed to COVID-19 but the lowest average temperatures over the timespan (e.g. France, Denmark, Italy, Spain, Sweden, Switzerland, the UK, and the US).

Figure 1.

Average ambient temperature in °C per country and time (t = 60, 90, 120), (N = 150 countries).

Figure 2.

Average COVID-19 attributed mortality per country and time (t = 60, 90, 120), (N = 150 countries).

Prediction of the pre-COVID-19 respiratory disease mortality

Table 2 displays the results of the first-stage modeling from Equation 1 (see supplementary material, section B for modeling diagnostics and predictors eligibility). The percentage of women and people aged 65 and older, the DALYs attributed to asthma, obesity prevalence, population density, HDI, and air pollution (pm 25), accounted for 59% of the variation of the pre-COVID-19 respiratory mortality. We predicted the adjusted pre-COVID-19 respiratory disease mortality based on Table 2 results.

Table 2.

First-stage analysis: log-linear regression results (N = 150 countries)

Pre-COVID-19 respiratory disease mortality per 100,000 people	β	SE
Women	0.049***	0.018
People aged 65 and above	0.098***	0.020
DALYs Asthma	1.089***	0.261
Obesity prevalence	−0.014*	0.007
Population density (population/km2)	0.000***	0.000
HDIa
Medium	−0.020	0.138
High	−0.138	0.191
Air Pollutiona
Medium	−0.126	0.085
High	−0.208*	0.122
Constant	0.270	0.913
R2	0.593
AIC	172.015

Notes. * 0.1 ** 0.05 *** 0.01. Robust standard errors were used. IL stands for inferior limit while SE standard error. a Terciles, using low groups as the reference category.

Cross-sectional and longitudinal analysis of temperature and COVID-19 number of deaths

Table 3 (section A) shows the main results of the cross-sectional multivariate analysis using negative binomial regression. Countries with higher ambient temperature had a significantly lower incidence risk ratio of COVID-19 death at t = 60, t = 90, and t = 120. An additional 1°C from the previous timestep decreased the incidence risk of COVID-19 death by 10% at t = 60 (IRR = 0.90; SE: 0.036) and 8% at t = 90 and t = 120 (IRRt=90 = 0.92; SEt=90: 0.04; IRRt=120 = 0.92; SEt=120:0.04). Our longitudinal analysis (section B) showed that after adjusting the model by pre-COVID-19 respiratory mortality, the variation in government’s stringency measures, countries’ hemisphere, and ambient temperature remained as a protective factor for the incidence risk of death attributable to COVID-19 over time (models 6 and 7). The approach derivation and tests for the main assumptions of the model are found in supplementary material, section E.

Table 3.

Second stage analysis: negative binomial regression results for the incidence of COVID-19 deaths

Section A. Cross-sectional negative binomial regression models (N = 150)
	Model 1 (t = 60)			Model 2 (t = 90)			Model 3 (t = 120)
	IIR	SE		IIR	SE		IIR	SE
Ambient temperature	0.902**	0.036		0.919*	0.041		0.917*	0.043
PRDM (%)	1.018	0.032		1.029	0.029		0.988	0.028
Δ Government measuresb	1.014	0.011		1.090***	0.033		1.016	0.023
Regionc	2.892	2.252		1.570	1.176		2.017	2.296
Constant	214.81***	382.321		1194.97***	2146.91		3308.99***	8299.39
Ln(alpha)	1.417	0.104		1.482	0.100		1.649
Pseudo R2:	0.0206			0.0320			0.010
AIC:	1633.501			1806.938			614,894.2
Section B. Longitudinal negative binomial regression models (N = 450)
	Model 4		Model 5		Model 6		Model 7
	IIR	SE	IIR	SE	IIR	SE	IIR	SE
Ambient temperature	0.930***	0.020	0.949**	0.029	0.932**	0.027	0.926***	0.026
PRDM (%)			1.026*	0.017	1.030*	0.018	1.032*	0.018
Δ Government measuresc					0.976**	0.011	0.977**	0.010
Regionc							2.207	1.094
Constant	2753.29***	1013.58	1006.26***	911.25	1322.91***	1272.83	513.005***	556.941
Chi2 (p-value):	10.44(<0.001)		30.24(<0.001)		37.55(<0.001)		44.48(<0.001)

Notes. * 0.1 ** 0.05 *** 0.01. IRR stands for incidence risk ratios. aPRDM stands for pre-COVID respiratory disease mortality adjusted. b Δ stands for the variation in the stringency government index between timepoints. cRegion stands for hemisphere of the country, ‘Southtern’ was used as reference. Sections B and D display results using GEE population-averaged model. Bootstrap standard errors calculated with 1000 iterations were used in all models.

Cross-sectional and longitudinal analysis of temperature and a high-level COVID-19 mortality

We compared different definitions for ‘high-level of COVID-19 mortality’ (HLCM) (see supplementary materials, section G, for model comparisons). We used the median + 0.4 SDs to analyse the countries with a HLCM because it fulfills the statistical criteria of 10 EPV in the model. Countries’ hemisphere was not added as independent variable in the models because of the low number of countries classified as being in the ‘Southern’ area at t = 60. Nevertheless, we ran an exploratory analysis by adjusting our main models by countries’ hemisphere, and it showed a consistent relationship between ambient temperature and HLCM over time (supplementary materials, section H).

At t = 60, t = 90, and t = 120, there were 30 (20%), 36 (2 4%) and 64 (43%) countries classified as HLCM, respectively (see supplementary materials, section G for the list of countries). For instance, Belgium, Switzerland, Spain, Luxembourg, Netherlands, Italy, the United Kingdom, and Ireland had a high level of COVID-19 mortality at the three timesteps.

Table 4 (section C) shows the main results of the cross-sectional multivariate analysis using the logistic regression approach detailed in Equation 2. Countries with higher ambient temperature had lower likelihood of HLCM regardless of the period. An additional 1°C from the previous timestep decreased the likelihood of HLCM by 7% at t = 60 (OR = 0.930; SE: 0.026) and at t = 90 (OR: 0.925; SE: 0.028), and by 8% at t = 120 (OR = 0.915; SE: 0.031). Pre-COVID-19 respiratory mortality was significantly related to a HLCM at t = 60. The variation in the government’s stringency measures was not significantly related to HLCM.

Table 4.

Second-stage analysis: logistic longitudinal regression results for high-level of COVID-19 mortality

Section C. Cross-sectional logistic regression models (N = 150)
	Model 8 (t = 60)		Model 8 (t = 60)		Model 10 (t = 120)
	OR	SE	OR	SE	OR	SE
Ambient temperature	0.930**	0.026	0.925***	0.028	0.915***	0.031
PRDM (%)a	1.050*	0.031	1.037	0.026	1.006	0.022
Δ Government measuresb	0.977*	0.012	0.995	0.018	0.098	0.016
Constant	0.224	0.231	0.426	0.43	3.495	3.55
Pseudo R2:	0.19		0.18		0.16
AIC:	130.2326		143.64		192.17
Section D. Longitudinal logistic regression models (N = 450)
	Model 11		Model 12		Model 13
	OR	SE	OR	SE	OR	SE
Ambient temperature	0.961**	0.019	0.985	0.021	0.964*	0.019
PRDM (%)			1.048***	0.010	1.040***	0.012
Δ Government measuresc					0.977***	0.006
Constant	0.827	0.299	0.158***	0.087	0.257**	0.142
Chi2 (p-value):	3.95 (0.06)		43.63 (<0.01)		49.15 (<0.01)

Notes. * 0.1 ** 0.05 *** 0.01. OR stands for odds ratios. aPRDM stands for pre-COVID respiratory disease mortality adjusted. b Δ stands for the variation in the stringency government index between timepoints. Sections B and D display results using GEE population-averaged model. Bootstrap standard errors calculated with 1000 iterations were used in all models.

Table 4 (Section D) presents the longitudinal model for COVID-19 mortality detailed in Equation 2. The unadjusted and fully adjusted model showed a relationship between ambient temperature and the probability of HLCM. Model 13 shows the relationship adjusted by pre-COVID-19-respiratory mortality, and the time point variation in government’s stringency measures. Over time, an additional 1°C from the previous month decreased the likelihood of HLCM by 4% (ORmodel14 = 0.96; SE: 0.019). Furthermore, between periods variation in the government stringency measures was a country protective factor for HLCM over time (ORmodel14 = 0.98; SE = −0.006), while pre-COVID 19-respiratory mortality was a risk factor for HLCM over time (ORmodel13 = 1.040; SE = 0.012). In sensitivity analyses, we changed the threshold for HLCM to the median + 2.0 SDs (EPV = 2); however, similar results were found (see supplementary material, section H).

Discussion

We analyzed the association between monthly ambient temperature and COVID-19 mortality across countries at the beginning of the pandemic accounting for epidemiological factors. We used the adjusted COVID-19 number of deaths and the high-level of COVID-19 mortality group of countries to understand how temperature has been related to them. We found that ambient temperature was associated with the adjusted COVID-19 number of deaths and a high-level of COVID-19 mortality (HLCM) using different model specifications.

The results of the relationship between ambient temperature and COVID-19 number of deaths are in line with the observational studies that have reported a negative relationship between them [15-17,20,43]. Our results showed that this negative relationship was present after 60 days since the first case confirmed. At the same time, ambient temperature was a protective factor for COVID-19 deaths over time. Our results were not altered using the variation in the stringency of government measures taken and countries’ hemisphere. Our results corroborate the WMO call that considered the previous literature as inconclusive and called for further analysis on the matter [21].

The ambient temperature might be a protective factor against the HLCM over time. Based on the dynamics of existing and previous infectious diseases, SARS-CoV-2 mortality may differ according to environmental changes because seasons have factors that determine the pathogen’s abundance, reproduction, and survival time within the environment, therefore, the community [44]. Three underlying hypotheses may drive the negative association between COVID-19 mortality and ambient temperature. Firstly, colder ambient temperatures could be linked to changes in population behavior; people spend more time indoors during colder weather. Consequently, crowded and poorly ventilated spaces, such as urban transit systems, could increase the viral load by the increased exposure to airborne and droplet transmitted pathogens from one person to another [45,46]. Secondly, the seasonal variability of the immune system’s functions affecting the host’s susceptibility to infection, such as seasonal variation of vitamin D and melatonin levels, are integral to upholding a strong immune system [47-51]. Low levels of these factors have been linked to significantly increased risks of viral upper respiratory tract infections, pneumonia, severe inflammations, and thrombosis, all of which have been frequently observed in patients with severe COVID-19 [52]. Thirdly, co-occurrence of infections may increase the severity of COVID-19 cases [53]. Colder seasons increase the morbidity and mortality of low respiratory tract infections and chronic respiratory diseases [46,54].

In the longitudinal analysis, we found that the variation of the government measures was significantly related to the incidence risk of COVID-19 deaths and the likelihood of HLCM. Government measures, including strict lockdown, may not be sufficient to stop the spread and reduce mortality especially considering that after t = 60 they were decreased by the countries observed. However, countries’ effective, opportune, and sustained reaction can help capitalize on higher temperatures’ protective role [2,3,6,19,20,22,24,26,27,35,55-57]. Other important factors must be considered to reduce the odds of high-level of COVID-19 mortality. These factors include countries’ economic resources, quality of care, healthcare coverage, demographic distribution, air quality, population-specific underlying conditions, and the prevalence of other respiratory diseases [4,16,19,26]. The relationship between obesity prevalence and respiratory disease mortality was negative in our results, and it might be possible driven by the negative relationship within countries with moderate HDI (r = −0.367). Previous literature has related obesity with increased risk of mortality [58]; however, a study analyzing mortality risk of COVID-19 patients in ICUs found a potential obesity paradox [59]. Moreover, our longitudinal findings suggest that countries should improve their efforts by implementing effective preventive measures to reduce respiratory disease mortality, which accounts for a vast disease burden in high-income countries (HICs) [35].

We also found that mainly northern HICs exhibited higher mortality rates attributed to COVID-19 during the observed period (e.g. France, Italy, Spain, and the UK). Some shared features are highlighted across these countries. They have had an earlier onset of infections, a greater proportion of older people, a higher burden of disease from chronic conditions, including cardiovascular diseases and pre-COVID-respiratory disease mortality, and the lowest ambient temperatures observed since their onset. Most of these characteristics represent high-risk factors for severe COVID-19 and attributed fatality rate [9,19,28,60-62]. Additionally, countries that reported HLCM and the greatest temperatures were mostly low- and middle-income countries (LMICs) with lower HDIs (e.g. Guatemala, Honduras, Panama) or with poor government performance in managing the pandemic (e.g. Brazil, Mexico). On the other hand, most African countries reported high or medium levels in ambient temperatures and the lowest number of deaths adjusted to their population size due to having younger populations, rapid action through the implementation of large-scale containment measures, low prevalence of chronic cardiovascular conditions, favorable climate and good community health systems, and lack of resources for epidemiological vigilance [63].

This article has limitations. First, we did not include other potential external variables which may impact transmission, and therefore, the number of deaths. For instance, population mobility might be related to local weather conditions, and to the stringency level derived from the government measures implemented. However, data on mobility was not widely available. Second, the existing missing data for LMICs may bias interpretations toward socioeconomic disparities. Third, we did not use excess deaths attributed to COVID-19 nor age-standardized mortality due to lack of data availability, especially for LMICs [64,65]. Further analyses should look at both measures combined to disentangle the links between them while trying to correct, contrast, and interpolate mortality estimates, specifically in countries with insufficient or null data published. Fourth, discretizing a continuous variable may complicate the results, so they must be interpreted cautiously [66,67]. Fifth, ambient temperature was an average measure for the entire country. Therefore, indoor temperatures may represent an unmeasured confounder, while countries with high variability of ambient temperature and wide geographical areas might be underrepresented. Sixth, given the complexity of the relationship examined, there were potential unassessed cofounders involved in the association between ambient temperature and COVID-19 mortality (E-value coefficient. = 1.36; Inferior CI = 1.23; see supplementary materials, section H) [68].

Considering these limitations, the strengths of the present study outweigh the shortcomings. We attempted to eliminate endogeneity biases accounting for pre-COVID characteristics. We have contrasted cross-sectional and longitudinal methods to test the linkage between our variables over specific time points and over time and using two different outcomes: the risk of COVID-19 deaths and the high-level of COVID-19 mortality. Using ecological data, we included a vast number of countries conducting a global analysis of the relationship between ambient temperature and COVID-19 mortality. Previous articles have only focused on incidence of deaths or continuous mortality; however, we have included high-probability risk, which has been often overlooked. Besides, we used COVID-19 attributed mortality ratio due to the limited testing capacity in some countries, especially in LMICs. Mortality measures may serve as an accurate indicator for COVID-19 spread in highly impoverished countries, but also in HICs. For instance, France, Italy, Spain, and the US have evidenced great numbers of underreported and undetected COVID-19 cases due to the great number of tests taken during patient’s hospitalization or before death occurrence [69]. Massive-scale testing to the wide population should be implemented instead.

Finally, we decided not to include any further time to analyze mortality at the early pandemic and avoid the potential indirect effects carried by the vaccination process started at the end of 2020. Since the vaccination process began, the Oxford Coronavirus Government Response Tracker has been updated. The evidence of this report should be useful to take faster and effective decisions under similar scenarios related to MERS-CoV and SARS-CoV infections [70].

The present study attempts to understand the cross-country relationship between COVID-19 mortality and temperature, accounting for government containment measures to reduce its spread. The protective role of ambient temperature on the incidence of COVID-19 deaths and the probability of a high-level of COVID-19 mortality over time remained when considering the stringency level of the governments’ measures to tackle the disease spread. We provided preliminary evidence for the relationship between the lag of monthly ambient temperature and the probability of high-level of COVID-19 mortality through a global study. Our findings support the call from the WMO to not taking government COVID-19 infectious containment decisions only derived from meteorological factors [21]. Conversely, the relaxation of COVID-19 related government measures should be based on the country’s public health capacity, community engagement, health system, and border control measures [6,19,22,24]. Moreover, a reinforcement on vaccine campaigns should be in place during warmer seasons, especially in those countries where vaccination strategies are still slow and incomplete.