Does climate help modeling COVID-19 risk and to what extent?

RESEARCH ARTICLE Does climate help modeling COVID-19 risk and to what extent? Giovanni ScabbiaID1, Antonio SanfilippoID1*, Annamaria Mazzoni1, Dunia Bachour1, Daniel Perez-Astudillo1, Veronica Bermudez1, Etienne Wey2, Mathilde Marchand- Lasserre2, Laurent Saboret2 1 Qatar Environment and Energy Research Institute, Hamad Bin Khalifa University – Qatar Foundation, Doha, Qatar, 2 Transvalor S.A, Sophia Antipolis, France * asanfilippo@hbku.edu.qa Abstract A growing number of studies suggest that climate may impact the spread of COVID-19. This hypothesis is supported by data from similar viral contagions, such as SARS and the 1918 Flu Pandemic, and corroborated by US influenza data. However, the extent to which climate may affect COVID-19 transmission rates and help modeling COVID-19 risk is still not well understood. This study demonstrates that such an understanding is attainable through the development of regression models that verify how climate contributes to modeling COVID- 19 transmission, and the use of feature importance techniques that assess the relative weight of meteorological variables compared to epidemiological, socioeconomic, environ- mental, and global health factors. The ensuing results show that meteorological factors play a key role in regression models of COVID-19 risk, with ultraviolet radiation (UV) as the main driver. These results are corroborated by statistical correlation analyses and a panel data fixed-effect model confirming that UV radiation coefficients are significantly negatively corre- lated with COVID-19 transmission rates. Introduction

The COVID-19 pandemic has triggered an extensive amount of research across different fields, from epidemiology to the social sciences. One of the issues under investigation pertains to the heterogeneous character of COVID-19 diffusion across diverse geographic areas. The debate on this matter has given rise to two main approaches. The first is primarily based on the epidemiological explanation of contagion dynamics in terms of factors such as frequency and intensity of physical contact among people and their viral charge. According to this approach, the only variables capturing the viral spread are contagion dynamics factors and the ability of policymakers to reduce the contagion rate by restricting population mobility. The second approach evaluates the possibility that contagion dynamics result from multiple fac- tors, including meteorological and environmental variables. Since the early stages of the COVID-19 pandemic, climate has provided an important refer- ence point to explain the spread of the virus. Just three months after the first outbreak in PLOS ONE PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 1 / 27 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Scabbia G, Sanfilippo A, Mazzoni A, Bachour D, Perez-Astudillo D, Bermudez V, et al. (2022) Does climate help modeling COVID-19 risk and to what extent? PLoS ONE 17(9): e0273078. https://doi.org/10.1371/journal.pone.0273078 Editor: Daniel Dunea, Valahia University of Targoviste, ROMANIA Received: January 22, 2022 Accepted: August 1, 2022

Published: September 7, 2022 Copyright: © 2022 Scabbia et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: Code for the statistical and machine learning analyses has been implemented in Python, while the econometric analysis was implemented in STATA. Source data and code to ensure the reproducibility of the analysis is available at the GitHub public repository: https://github.com/QEERI/DEFCOV. Funding: This study reported in this paper was funded by grant RCC-2-044 from the Qatar National Research Fund (QNRF, https://www.qnrf. org), awarded to Dr. Antonio Sanfilippo. QNRF RCC-2 call was specifically focused on COVID

Wuhan, China, Bukhari and Jameel [1] reported that 90% of COVID-19 cases were recorded through 22 March 2020 in world areas with colder and less humid late winters and early springs (temperatures of 3–17˚C, and absolute humidity of 4–9 g/m3). As recognized by the authors, these early data were likely to be biased by minimal testing per capita in tropical coun- tries. After over two years into the pandemic, 226 countries across the globe have reported a total of over 500 million confirmed COVID-19 cases [2]. We now have a clearer picture of the global distribution of this disease. There are numerous examples indicating that the cooler sea- son in the northern hemisphere may have favored the spread of the disease, while warmer and more humid weather in late spring and summer has seen a substantial and rapid decline in transmission numbers, once the different containment strategies adopted worldwide are taken into consideration [3, 4]. Prior work on the betacoronavirus genus shows that viruses similar to SARS-CoV-2, such as HCoV-HKU1 and HCoV-OC43, also display seasonal peak trends [5]. The spread of previ- ous coronavirus epidemics, such as SARS-related and the Middle East respiratory syndrome (MERS)-related coronaviruses, have also displayed correlations with seasonal weather changes [6–9]. More specifically, outbreaks of respiratory virus infections are commonly associated with seasonality, with peaks during the winter months. Data from the 1918–19 Flu Pandemic support this correlation. Peak infection/mortality during the 1918–19 Flu Pandemic occurred

in the winter months and waned as solar radiation and absolute humidity increased from late March onward. The resurgence of mortality in the winter of 1919 and its decline with the arrival of spring confirm this trend. These observations are corroborated by US influenza data relative to the last decade where the percentage of patient visits for influenza-like illness consis- tently grows in the winter months. Regression analyses correlating meteorological factors with the spread of SARS-CoV-2 have reported contrasting results. A number of studies have found a negative correlation between temperature and the spread of COVID-19 in China and the US [10–16], Italy [17, 18], Spain [19], Mexico [20], Brazil [21], Latin America and the Caribbean [22], and worldwide [23–27]. By contrast, Xie and Zhu [28] suggest a positive relationship between temperature and the spread of COVID-19 (up to 3˚C), and city-level data from Oslo, Jakarta and, five Brazilian cit- ies show a positive correlation between COVID-19 transmission and higher temperatures, and negative correlation with precipitations [29–31] which typically lead to lower temperatures. Studies on the impact of humidity have also given conflicting results with reports of positive correlation by Jiang et al. [10] for China, and the opposite in studies by Ward et al. [32] for Australia, Qi et al. [11] and Wang et al. [12] for China, Yin et al. [21] in Brazil, and Ju¨ni et al. [33] and Aboura [34] worldwide. Several global analyses and reviews on the impact of climate

on COVID-19 spread [35–39] have also reported conflicting results. Other meteorological factors investigated with reference to COVID-19 transmission include wind and solar radiation. Most studies focusing on the impact of wind speed on the incidence of COVID-19 cases have found a negative correlation [40–42]. Studies focusing on solar radiation in Brazil [43], Australia [44], and worldwide [27, 45, 46] have found that high exposure to solar radiation including UV is negatively correlated with the occurrence of COVID-19 cases. S1 Table provides a summary of the peer-reviewed research on the interrelationship between COVID-19 and meteorological/climatic factors discussed in this section. Overall, there is no conclusive evidence that climate plays a role in the spread of COVID-19 [47–50]. This is probably due to the fact that studies on the impact of climate on COVID-19 transmis- sion have been piecemeal (e.g., limited to country or administrative division-level data), have only taken into account a few climatic parameters, and have not considered the impact of socioeconomic factor, as remarked in Mecenas et al. [4]. The present study addresses these PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 2 / 27 research, and QNRF encouraged wide dissemination of the study’s results. The study’s results and findings are solely the responsibility of the authors. Competing interests: The authors have declared that no competing interests exist.

impediments by investigating the contribution of meteorological factors to modeling COVID- 19 transmission at the global scale with reference to reported COVID-19 cases in 196 countries over a 14-month period, using socioeconomic, environmental, and global health factors as control variables. The study employs three complementary approaches to measure the correla- tion between reported rates of COVID-19 transmission and selected meteorological, socioeco- nomic, environmental, and global health factors. The first approach is based on the statistical analysis of the rank correlation of every factor with the number of daily confirmed COVID-19 cases. Following is a feature importance analysis that uses Shapley Additive Explanations (SHAP) to elucidate the COVID-19 rate predictions of a regression model by computing the contribution of each feature in the model to the prediction task. The third and last approach is based on econometric analysis with panel data fixed-effect regression models. The three approaches supply independent analytic evidence of the relationship between confirmed COVID-19 cases and meteorological factors. Of particular interest is the complementarity between the machine learning analysis, which is intent on prediction, and the econometric analysis, which focuses on explanation. As discussed in the literature, the use of a hybrid approach where machine learning modeling is paired with econometric analysis can help address relative weaknesses in the two methods by leveraging relative strengths. Further details

are provided in the Methodology subsection below. Materials and methods Data sources and pre-processing The data used in this study include epidemiological, socioeconomic, environmental, global health indicators, and meteorological variables–see S2 Table. All population-related variables are converted to percentages of the total population per country. Parameters that behave as time-invariant variables during the period of focus for this study, e.g., socioeconomic variables, are used as control variables. Epidemiological data on the cumulative number of confirmed COVID-19 cases were retrieved for the period from 23 January 2020 to 21 March 2021 from two main sources: the data repository from the Johns Hopkins Center for Systems Science and Engineering [51], and the Corona Data Scraper online data service that retrieves COVID-19 Coronavirus cases data from verified sources worldwide and adds population data on a daily basis (coronadatascraper. com). Data from these two sources were merged to create the first dataset. In addition to coun- try-level records, data at the regional or state level were included when available. We derived the number of daily registered COVID-19 cases by differencing entries in the initial dataset. Data points characterized by a modified Z-score > 3.5 and with values above the third quartile or below the first quartile and 1.5 times greater than the interquartile range were discarded as outliers together with other inconsistent data points (e.g., negative values) [52]. From the daily

values, we retrieved 3-day and 7-day moving averages as features to incorporate an auto- regressive component in the learning model. Only records reporting more than 10 daily cases were included in the analysis. The time of exposure to the pandemic for each country was also calculated as the cumulative temporal distance from the first registered case in the country. We selected data through March 2021, when most countries began their vaccination cam- paigns, to avoid including vaccination as an additional impact on COVID-19 transmission rates. Such inclusion would have created discrepancies with earlier data and introduced incon- sistencies emerging from the adoption of diverse vaccination strategies worldwide. We used several socioeconomic time-invariant data sources including demographic infor- mation, technology adoption rates, and Gross Domestic Product per-capita (GDPP) as control factors in the cross-sectional fixed effect analysis. Demographic, population density, and PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 3 / 27

population age data were derived from the 2019 Population Division dataset compiled by the Department of Economic and Social Affairs of the United Nation (UNDESA) [53] and the World Bank indicators database [54]. Information for geographical locations not included in the UNDESA dataset was retrieved online from national official sources. Rates of internet users, subscribers to mobile telephony services, and the number of secure Internet servers were retrieved from the World Bank indicators database. These technology adoption variables are used as proxies for the capacity of different countries to provide smart-work environments under lockdown, create awareness and keep the population updated about the development of the pandemic, and support effective contact tracing applications. GDPP data at constant price purchasing-power-parity were sourced from the International Monetary Fund’s World Eco- nomic Outlook Database [55]. Environmental indicators retrieved from the World Bank database include population- weighted exposure to ambient PM2.5 pollution, carbon dioxide, methane, nitrous oxide emis- sions, and greenhouse gas emissions. These time-invariant variables were used as control indi- cators of the degree of pollution in each country, on the assumption that long-term exposure to pollutants may increase the risk of contracting COVID-19 given that outdoor air pollution has been positively correlated with respiratory diseases. Health indicators included time-invariant variables such as the general Global Health Secu-

rity (GHS) index, the GHS detect and prevent scores, diabetes prevalence, and the number of hospital beds for both acute and chronic care. GHS provides a country-level score of health security and was used as a proxy variable for a country’s capability to prevent and mitigate infectious diseases. For the purpose of this study, only the “detect” and “prevent” GHS catego- ries were used, which focus on a country’s readiness to promptly identify, report, and antici- pate disease outbreaks of potential international concern [56]. Health indicators relative to diabetes prevalence and the number of hospital beds for both acute and chronic care were retrieved from the World Bank database. These variables serve as proxies for population health status and public health preparedness. Differences in intervention responses by governments to mitigate the pandemic were accounted through a variety of indicators from the Oxford COVID-19 Government Response Tracker (OxCGRT) project [57, 58]. The OxCGRT variables include policy information on school closures (C1), workplace closures (C2), cancellation of public events (C3), restriction on gatherings (C4), closure of public transports (C5), lockdowns (C6), and restriction on inter- nal (C7) and international (C8) movements and travels. In addition, the OxCGRT database provides health system policy data on the presence of public information campaigns (H1), test- ing policy (H2), contact tracing (H3), and facial covering policies (H6). Finally, we used two

more variables that are calculated as a weighted aggregation of the single C and H indices: the stringency and containment & health indices. The first variable reflects the strictness of “lock- down style” policies that primarily restrict people’s behavior. The second combines “lock- down” restrictions with measures such as testing policy and contact tracing, short-term investment in healthcare, and investments in vaccines. OxCGRT data are provided at a coun- try-level, with subnational data available for Canadian provinces, US states, and UK regions (New England, Northern Ireland, Scotland, and Wales). Meteorological variables were obtained from two main sources: the Modern-Era Retrospec- tive analysis for Research and Applications, Version 2 (MERRA-2), a gridded reanalysis data- set produced by NASA’s Global Modeling and Assimilation Office (GMAO), and the Copernicus Atmosphere Monitoring Service (CAMS) with specific reference to the McClear Clear-Sky Irradiation model. The MERRA-2 data include daily averages of the original hourly data at a spatial resolution of approximately 50 km, temperature in ˚C, relative humidity at 2 m above ground in %, short-wave solar irradiation (total of the day) in Wh/m2, pressure at 2 PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 4 / 27

m from ground level (station pressure) in hPa, wind speed at 10 m above ground in m/s, and rainfall in mm. Temperature is combined with relative humidity to derive measures of absolute humidity in g/m3 [59]. Data obtained from CAMS include downward UV Radiation at the Surface (UV) in J/m2, and Particulate Matter (PM) concentrations (PM2.5 and PM10) for 3-hour periods at a spatial resolution of 40 km. UV exposure can have a sterilizing effect [60] and ultraviolet B light (UVB), which is present in small amounts in natural sunlight, is known to rapidly inactivate SARS-CoV-2 on surfaces [61]. Data on particulate matter [62], originally in kg/m3, is converted to micrograms/m3 and it can provide preliminary evidence on the rela- tion between air quality and the chronicity of exposure to the viral infection. Coccia [63] and Bloise and Tancioni [64] both suggest that air pollution may have accelerated the transmission rate of COVID-19 in northern Italy, even though the viability of infectious viruses embedded in suspended aerosol particles is still under debate [65]. Streaming access to MERRA-2 and CAMS was provided by Transvalor S.A. For each loca- tion considered in the study, we derived the geographical centroids of the country’s most pop- ulous cities [66]. For country-level locations, we considered the ten most populated cities, while for admin-level locations (sub-country, i.e., states, territories, provinces, cantons, etc. as appropriate per country) we considered the five most populous cities. We used the corre-

sponding latitude and longitude coordinates of each city to query the climatic information from MERRA-2 and CAMS through Transvalor’s SoDa data service (http://www.soda-pro. com/) and we finally derived the resulting time-variant meteorological data by averaging across the different cities for each location. After merging all the various sources, the resulting dataset includes data on 196 countries covering 96% of the world population and 97.6% of worldwide confirmed COVID-19 cases (123,491,126 –at the period of the study). Data for 28 of these countries are detailed in the data- set at a state or regional level (i.e., admin-level) for the available periods (see S1 File for a detailed list). We only consider country-level epidemiological data for the remaining 152 countries, even though admin-level data are available from coronadatascraper.com, in order to keep a certain level of minimum comparison between the locations under study in terms of overall population size. S2 and S3 Tables provide detailed information on the variables used for the study and their descriptive statistics. The COVID-19 mean incubation period, defined as the time period ranging between expo- sure to the virus and the onset of the illness, is estimated by WHO at 5–6 days (median 5.1 days, 95% Confidence Interval (CI): 4.5 to 5.8 days) [67]. According to Lauer et al. [67], 97.5% of those who develop symptoms will do so within 11.5 days (CI: 8.2 to 15.6 days) from the day of infection. Moreover, the results of a COVID-19 PCR test have been known to take up to an

average of 3–4 days (3.6 days according to Cereda et al. [68]), particularly during the initial months of the pandemic. For these reasons, the number of new COVID-19 cases that are offi- cially announced each day corresponds to a time window of infection that spans from a few days up to potentially two weeks earlier. To account for this timeframe uncertainty and test the robustness of our results, the analysis is carried out with moving averages for both the time- variant meteorological and the policy variables of different duration: 5, 7, 10, 12, and 14 days (minimum length of 3 days). In reporting results, the number of days determining the moving average of time-variant variables is encoded as either a number suffix (e.g., Temperature_7), or a suffix variable indexed to a specific numeric value (e.g., Temperature_T, . . . T = 7). Limitations and assumptions As for other data-driven studies on COVID-19 transmission, the present analysis relies on rec- ords whose quality varies across sources, due to heterogeneous collection and reporting PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 5 / 27

practices worldwide. Data quality, extensiveness, and uniformity are therefore subject to a cer- tain degree of uncertainty. Moreover, reports of confirmed COVID-19 cases tend to underesti- mate the actual number of infections because of asymptomatic patients and undetected COVID-19 deaths. For ease of purpose, we will assume that the number of confirmed COVID-19 cases is monotonically related to the true number of infections, recognizing that this is a simplification that may limit the significance of this study’s results. An additional but less impactful limitation emerges from the normalization of meteorologi- cal variables across different points within each area of interest where exposure to the virus has occurred. As previously described, weather conditions have been averaged across the ten most populous cities for each country-level location and the five most populated cities for the admin-level locations. This approach should provide a good approximation of the overall weather conditions under which the viral spread has occurred. We reason that population den- sity may have played a significant role in the spread of COVID-19. Furthermore, it is also likely that the majority of testing was carried out in major urban areas. Finally, selecting up to ten country-level and five admin-level cities can be expected to evenly spread the sampling across the most populated areas, as such an approach provides an average coverage of 75% or more of the entire population for most of the locations considered (S1 Fig). Methodology

The method used in this study combines three approaches to capture dependencies between confirmed COVID-19 cases and climate factors (see S2 Fig for a visual diagram summarizing the approach): • A statistical analysis based on Spearman’s and Kendall’s correlation coefficients • A machine learning model based on the Gradient Boosted Regression Tree algorithm paired with the Tree SHAP algorithm to perform feature importance analysis • An econometric model based on fixed effect panel regression analysis. The use of three distinct approaches is intended to provide independent analytic evidence. Of particular interest is the complementarity between machine learning and econometric approaches, where the first is intent on prediction while the latter focuses on explanation. As discussed in the literature, the use of a hybrid approach where machine learning modeling is paired with econometric analysis can help address relative weaknesses in the two methods by leveraging relative strengths [69–73]. For example, machine learning is better equipped to take advantage of structural heterogeneity in training data to make short-term predictions, whereas econometric methods are better at capturing long-term trends [73]. We can therefore expect that the results of machine learning and econometric analysis are not always going to coincide [73]. This lack of overlap points to the areas of relative improvement that can be obtained through a functional integration of the two methods. While appealing, such integration is chal-

lenging and largely remains a goal to be achieved, for which a better understanding of the dif- ferences and relative strengths/weaknesses is required [71]. In this regard, the present study contributes to advancing our understanding of the specific complementarities between machine learning and econometrics in a new domain of inquiry. The statistical analysis in turn provides the baseline due to its more basic analytic capacity in dealing with non-station- ary processes compared to machine learning and accounting for long terms trends compared to econometric analysis. PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 6 / 27

Finally, the validity of the approach adopted in this study is corroborated by the research framework for linking environmental and weather factors to the incidence of COVID-19 pro- posed in a recent study published by Zaitchik et al. [74]. Statistical analysis The statistical dependency of confirmed COVID-19 cases from environmental and meteoro- logical regressors is first performed by calculating and comparing Spearman and Kendall rank-order correlation coefficients. These coefficients provide a nonparametric measure of the monotonicity (i.e. strength and direction) of the relation between the number of confirmed COVID-19 cases (the output/dependent variable) and the input environmental and meteoro- logical variables. Unlike the Pearson correlation, the Spearman and Kendall rank-order corre- lations do not carry any assumptions about the normal distribution of the data and the linearity between the variables. The statistical significance of the association between input and output variables is determined using the two-sided p-value in order to measure both decreasing and increasing departures from the null hypothesis. Spearman’s and Kendall’s cor- relation coefficient values (ρ and τ) can range from +1 to -1. The sign of the coefficient indi- cates the direction of the association of ranks (+ positive,—negative), while its absolute value expresses magnitude. The closer the coefficient is to zero, the weaker the association between the ranks: an absolute value between 0.5 and 1 is considered to provide a strong correlation,

0.3 to 0.5 a moderate correlation, 0.1 to 0.3 a weak correlation, and <0.1 no correlation. Machine learning modeling The aim of the machine learning analysis is to assess the relative feature impact of factors con- tributing to COVID-19. Feature impact is computed by applying the Tree SHAP algorithm to a Gradient Boosted Regression Tree (GBRT) model. GBRT is an additive stochastic model that combines multiple sequentially connected weak learners (regression trees) in a way that each new learner fits the residuals from the previous step to optimize the overall predictive perfor- mance [75]. The resulting model can describe multiple nonlinear interactions and partial dependency with sufficient flexibility, remarkably high predictive accuracy, and robustness to missing data and outliers. The study uses the open source xgboost Python library which offers a highly efficient, flexi- ble, and portable implementation of GBRT [76]. The xgboost algorithm provides several ways to control overfitting, i.e., when the model fits the training data so closely that it fails to provide useful predictions when applied to new data. The first is to constrain the maximum depth of individual trees used in the boosting process to modulate the degree of feature interactions that the model can fit. The second is to control the number of samples that each tree leaf can contain to avoid forming imbalanced leaves that have a single or too few data points. The third and most important way is to control the learning rate. Overfitting is also reduced through the

use of randomization into the tree building process by subsampling the training set before deriving each tree, and subsampling features before searching for the best node split. Finally, xgboost provides a parameter that enables model regularization using "across trees" information. As a first step, we optimized the model hyperparameters using a grid-search method com- bined with a cross-validation approach specifically designed for this study. We first selected only the data records of geographical locations that presented at least 90 data points (about 3 months’ worth of data) filtering out about 12,000 observations from a total of 77,300. This strategy is intended to select a time window size that presents sufficient seasonal variation for each location. We then randomly selected 24 locations to be designated as a test set to measure PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 7 / 27

the final unbiased performance of the optimized model. Each location in the filtered dataset has on average 273 data points, thus resulting in an overall test set size of ~6,500 observations (about 10% of the starting dataset) and a training set of 60,000 observations. We used a grid- search approach to find the best model hyperparameters validated on a further 10% split share of the training set. This validation set was derived by randomly selecting 22 locations for a total of ~6,000 observations, thus leaving 54,000 data points for training the model with the specific parameters under validation. For each parameter, we repeated this procedure 5 times to assure that the resulting validation error score (in terms of symmetric mean absolute per- centage error—SMAPE) would converge. During each reiteration, we re-selected a new ran- dom set of 22 locations from the overall training set as the new validation set. This step led to a measurable improvement in the prediction accuracy over the same algorithm initialized with the default hyperparameters values (base model). Fig 1 shows the performance gain in terms of lower SMAPE error of the optimized model compared to the base model computed on the first randomized test set, which was never seen by both models. On average this optimization procedure results in a ~5–10% lower SMAPE. Table 1 summarizes the set of hyperparameters leading to the best evaluation results. Table 2 reports the training, validation, and test performance of the GBRT model (default and

optimized) compared to two deterministic baselines: the prediction obtained by using a 7-day moving average, and a persistence model, where the value of the predicted dependent variable is assumed to be the same as the previous day. For comparison, we also report the performance of other regression models such as Lasso, Elastic Net, and Random Forest after a cross-vali- dated tuning of their hyperparameters. The optimized GBRT model outperforms all other models on the test set with a mean SMAPE error of 7.5% lower than the base model, 2.4% smaller than the Random Forest model, 11% better than the Lasso model, and 24% lower than the Elastic Net model. High error values for Lasso and Elastic Net are likely related to their lower model complexity that prevents the proper learning of data interrelationships. Random Forest produces a comparable accuracy to xgboost, but it shows overfitting on the train set. Finally, when compared to the two deterministic baselines the tuned xgboost model produces predictions on the test set that are 20.5% and 4% more accurate than the 1-week moving aver- age and persistence models respectively. Fig 2 shows model accuracy distribution as a function of the daily cases grouped in different intervals for all the locations considered in this study. The width of each boxplot is propor- tional to the number of observations included in the specific range (n), which is also reported below the label of each interval. The cross-sectional median error of the model decreases with Table 1. Xgboost hyperparameter tuning result.

Hyperparameter Tuning range Best value Learning rate (eta) 0.001–0.3 0.1 Maximum depth 3–10 10 Minimum sum of instance weight (hessian) 1–10 7–8� Gamma 0–0.4 0.2� Subsample ratio of the training instances 0.5–1 1� Subsample ratio of columns when constructing each tree 0.3–1 1� Lambda 1 1 Alpha 0 0 Number of boosting rounds (validated) ~30 Learning objective function (regression) Squared error Custom evaluation metric (for training and validation) SMAPE Early stopping rounds 10 �Irrelevant to the model performance. https://doi.org/10.1371/journal.pone.0273078.t001 Table 2. Regression modeling performance comparison. Model Mean SMAPE Train Validation Test GBRT (with hyperparameter optimization) 28.3% 30.6% 29.0% GBRT (no optimization) 33.5% - 36.5% Lasso� 45.7% 41.9% 40.1% Elastic Net� 60.2% 60.5% 53.3% Random Forest� 11.9% 31.2% 31.4% 1-week moving average - - 49.5% Persistence (previous day) - - 34.9% �With hyperparameter optimization https://doi.org/10.1371/journal.pone.0273078.t002 PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 9 / 27

the increase of the values of the dependent variable, going from 35.1% when daily cases are < 100 down to 17.5% with daily cases >500. To obtain a more complete estimate of model accuracy we tested the performance of the optimized GBRT model for each distinct location. We selected areas from different climatic zones that presented high numbers of daily cases during one or several contagion waves throughout the year. The results of this evaluation reveal error rates ranging from 12.9% to 26.8%, as shown in Fig 3. Once the GBRT model is trained, the relative ranking of model parameters is obtained through the SHAP method. Tree SHAP is an algorithm that computes SHAP values for Deci- sion Trees models such as GBRT. SHAP (SHapley Additive exPlanation) [77, 78] uses a game theoretic approach to explain the prediction for each instance as a sum of contributions from its individual feature values. This type of analysis does not identify causal correlation, but it is still a useful metric to capture relative feature importance. Econometric analysis The econometric analysis of the association between confirmed COVID-19 infections and cli- matic factors is carried out using the multivariate equation in (1) which implements a panel data approach based on a fixed-effects model [79]. In this equation, the dependent variable ln_daily_casesi,t expresses the number of daily cases of COVID-19 cases on a logarithmic scale for location i and time index t. We use the log-transformed version of the dependent variable

on the assumption that by doing so the variable becomes log-normal conditional on all the covariates and thus allows us to limit the heteroscedasticity of the estimated residuals. β0 is the Fig 2. Boxplot of the SMAPE distribution as a function of intervals of number of COVID-19 daily cases. https://doi.org/10.1371/journal.pone.0273078.g002 PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 10 / 27

regression intercept, while βn represents the regression slope coefficient of each respective regressor. We include a vector of cross-sectional unit fixed effects ci to account for all time- invariant factors across a location that affect the local growth rate of infections, such as differ- ences in demographics, socioeconomic status, culture, and health systems. This is an impor- tant feature since it allows to partial out heterogeneous omitted factors that might be correlated with the dynamics of contagion and the daily cases. We also include a vector of (daily) time fixed effects λt to absorb the autoregressive component specific to the COVID-19 spread growth and to account for the presence of any potential seasonal bias. Since the effect of variables that behave as time-invariant factors for the period of focus (e.g., socioeconomic, environmental, and some global health indicator variables) would be absorbed in the intercept for collinearity due to the use of time fixed effect regression, these variables were omitted from the analysis. Finally, we cluster the standard errors ui,t at the entity-level to account for error correlation within each location. ln daily casesi;t ¼ b0 þ X N n¼1 bn Explanatory variablesn;i;t þ ci þ lt þ ui;t ð1Þ Fig 3. Examples of modeling performance for the optimized GBRT model. https://doi.org/10.1371/journal.pone.0273078.g003 PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 11 / 27

Results We designed a cross-national data-oriented study using global records of confirmed daily cases of COVID-19 to examine the association between the pandemic growth and climatic conditions, using several socioeconomic, environmental, and global health factors as control variables. We first apply Spearman’s and Kendall’s rank-order correlation analysis to the selected data to derive a first estimate of the statistical relationship between each explanatory variable and COVID-19 transmission at each geographic location under study. We then use machine learning techniques to model and understand the relative importance that climatic and control variables have on the spread of COVID-19. We finally use panel data econometric analysis to estimate the impact of climatic conditions on COVID-19 daily rates and test the efficacy of different containment policies. Statistical correlation analysis Fig 4 provides the distribution of Spearman’s rank correlation coefficients (ρ) that model the dependency of COVID-19 daily rates on the environmental and meteorological explanatory Fig 4. Spearman’s coefficient. https://doi.org/10.1371/journal.pone.0273078.g004 PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 12 / 27

variables used in this study. Correlation coefficients are calculated for each location and time series. Results are clustered based on the specific location and consider only geographical areas with at least 90 data records (at least 3 months of data). Correlation coefficients greater than 1.5 IQR (interquartile range) that are below the first quartile or above the third quartile are considered outliers and reported in Fig 4 as scatter points. Only locations that showed statisti- cal significance (P < 0.01) were considered and displayed in the descriptive analysis reported in Fig 4 (see S2 and S3 Tables for variable description and statistics). The Spearman’s rank correlation results suggest that solar irradiation and UV emission are strongly negatively correlated with COVID-19 spread (median ρ equal to -0.55 and -0.56, respectively). Temperature and absolute humidity also show a negative correlation, but with a weaker amplitude (at median values of -0.42 and 0.39, respectively; moderately correlated). The other meteorological and air-quality factors do not show a significant association with COVID-19 transmission (low |ρ| and large IQR: high standard deviation). PM2.5 concentra- tions and pressure register a positively weak correlation. All other variables present a weakly negative correlation. These results are corroborated by Kendall’s rank correlation analysis shown in Fig 5. Fig 5. Kendall’s coefficient. https://doi.org/10.1371/journal.pone.0273078.g005 PLOS ONE Does climate help modeling COVID-19 risk and to what extent?

PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 13 / 27

Feature importance analysis Relative feature predictive value is computed by applying the Tree SHAP (SHapley Additive explanation for tree-based machine learning models) algorithm to a Gradient Boosted Regres- sion Tree (GBRT) model (see Methodology section for details). The GBRT model was trained on the whole dataset including all geographical locations independently from the number of its data records. Epidemiological, meteorological, socioeconomic, environmental, and global health indicator indices were used as explanatory variables, while the number of Covid-19 daily cases served as the dependent variable. See the Methodology section for a detailed defini- tion of the variables. S2 and S3 Tables present a brief description of the variables included in the study and their summary statistics. The results of feature importance analysis suggest that climate plays a meaningful role in modulating the dynamics of the COVID-19 pandemic, as shown in Fig 6 where feature impor- tance is ranked in terms of logarithmic mean absolute SHAP values. SHAP values relative to the average number of previous COVID-19 cases as a predictor for the current number of COVID-19 cases were computed but have been omitted in Figs 6 and 7 due to their obvious relevance in order to focus on the other variables. All meteorological and air quality factors score at similar levels of importance, showing that there is no dominant predictor. UV radia- tion is the meteorological factor with the greatest SHAP value, confirming the results of the

statistical analysis where UV radiation was the most highly correlated factor with COVID-19 cases. Socioeconomic, environmental, and global health indicator variables all show minor impacts other than population density, and the annual carbon dioxide emissions (a time- invariant proxy for the country’s overall air quality). Intervention and health system policies, described by the different OxCGRT indices, all score similar to or slightly lower than meteoro- logical factors. Although not included in Fig 6, the average number of previous COVID-19 cases results having the greatest impact (i.e. highest SHAP value) in line with results reported in the current literature on COVID-19 and other coronaviruses [80]. Fig 7 shows the intensity and polarity of specific SHAP variable values for each data point (a dot in the plot) with reference to the dependent variable (daily COVID-19 cases). The red- to-blue color scale indicates magnitude (high/red vs. low/blue). Position on the horizontal axis signals polarity (negative vs. positive). The resulting analysis suggests that higher UV radiation is significantly correlated with lower occurrence of COVID-19 cases (the dependent variable), while the other meteorological factors show a much weaker contribution. Temperature shows a weakly positive correlation, but the results are not consistent. Absolute humidity also displays a weakly positive correla- tion, despite the expected impact according to earlier studies. Rainfall appears to be negatively

correlated with the dependent variable, but there are too few observations for high rainfall val- ues in the dataset to properly confirm this result. All the other meteorological factors do not exhibit a distinct correlation directionality or significant impact. Among the intervention policies, both the OxCGRT stringency and containment & health indices display a weakly negative correlation suggesting that more stringent prevention mea- sures have been beneficial in mitigating the spread of COVID-19, at least for a significant number of countries. Policies on school closures, stay-at-home requirements, and testing reveal a somewhat lower weakly negative correlation. In order to provide a more detailed analysis, we have also compared feature importance analysis results of locations in the northern and southern hemispheres. As the SHAP values in S3 Fig show, the results for each hemisphere are in line with those of the global analysis (Fig 6). Meteorological conditions are still the more crucial factors when compared to intervention policies, with UV still among the most prominent features for both datasets. As shown in S4 PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 14 / 27

Fig 6. Feature importance summary plot. Mean absolute SHAP value (in log scale) of each variable showing the average impact on the model output magnitude. https://doi.org/10.1371/journal.pone.0273078.g006 PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 15 / 27

Fig 7. Feature impact scatter plot. SHAP value of each variable for all the single observations as a function of their relative value. The color transition on the vertical axis indicates value strength (red/high to blue/low). https://doi.org/10.1371/journal.pone.0273078.g007 PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 16 / 27

Fig, the intensity and polarity of SHAP values for the explanatory variables in separate hemi- sphere are also in line with those for the entire globe (Fig 7), although for some variables in the southern hemisphere polarity this similarity is not as explicit as for the northern hemisphere. For example, while UV still shows a negative association in both hemispheres, the distribution of SHAP score for the southern hemisphere is less marked than in the northern hemisphere. By contrast, rainfall displays a clear negative correlation in the southern hemisphere that is not clearly visible in the northern hemisphere. For completeness, we include the feature importance results for the Lasso, Elastic Net, and random forest models at a global scale (S5 Fig). All point to UV as one of the most important parameters. For Lasso and Elastic Net we report regression coefficients. For the random forest tree, we use the Gini importance score. Econometric analysis The econometric analysis was carried out using a panel data fixed-effects model. Confirmed daily cases of COVID-19 in log-scale were regressed against climate and air-quality factors, with reference to cross-sectional and time fixed effects. We could not add in the regression the moving averages of the dependent variable as we did for the machine learning analysis because this inclusion would violate the assumptions underlying the fixed effects estimator. If the inde- pendent variable is correlated with the error term in a regression model (endogeneity), then

the estimate of the regression coefficient would be inconsistent. Moreover, adding one or more autocorrelated terms to the regression would remove most of the model variance, mak- ing the effects of the other independent variables less significant (leading to smaller βn and larger standard errors). Note that the nominal magnitude of the regression coefficient of every single explanatory variable is likely to be biased due to the undetermined confounding effects. For this reason, we mainly focus our discussion of the results on the significance and polarity of the coefficients. Before proceeding with the econometric analysis, we test our data for stationarity, since non-stationary data may lead to spurious regression results thus falsely indicating the existence of a relationship between two variables [81]. For each time series variable considered in the econometric study, we run the Fisher-type unit-root test based on augmented Dickey-Fuller tests at 0 and 1 lag. Table 3 provides a summary of the results for the variables with moving averages at 7 days. We omit the test results for other window sizes for which we obtain the same outcome. The four tests all strongly reject the null hypothesis that all the panels contain unit roots for each variable under consideration and we can therefore proceed with the analysis. The statistical significance of the regression coefficients is computed by clustering the regression standard error at the country/admin level, to account for error correlation within

the geographical areas where our unit of observation was collected. For each location, we select days with a minimum of 10 confirmed COVID-19 cases, and we limit our analysis to the loca- tions with at least 90 data records (3 months’ worth of data), in line with the other analyses. The results are reported in Table 4 for all the T-day moving averages and time-variant regres- sors (see Methodology section). The regression model has a R2 value of 0.73 over 65,369 obser- vations, which shows that the independent and dependent variables are significantly correlated. UV radiation shows strong negative correlation with COVID-19 spread, while tem- perature has a positive association, in line with the statistical correlation and feature impor- tance analyses discussed in the previous two sections. Both results are statistically significant (P < 0.01). For other climatic factors, the econometric analysis is congruent with the feature importance analysis, but either the coefficients are not statistically significant, or the PLOS ONE Does climate help modeling COVID-19 risk and to what extent? PLOS ONE | https://doi.org/10.1371/journal.pone.0273078 September 7, 2022 17 / 27

magnitude of their regression coefficients is comparable to the standard error (e.g., PM2.5 and PM10). Rainfall seems to have significance only at shorter moving averages (5 days). For the remaining meteorological factors, different moving averages (5, 7, 10, 12, and 14 days), which relate climatic variables to incubation periods of diverse duration, do not seem to influence the overall result of the econometric analysis. Testing for the added effect of intervention policies requires a more in-depth analysis. The enactment of restrictions and the pandemic peak tend to vary from country to country due to the diversity, severity, and enforcement ability of the containment policies implemented. For this reason, we limit our analysis to the OxCGRT stringency and containment & health indices which allow us to capture the general level of restriction without focusing on the intervention policies of each government. Moreover, most countries have applied distinct levels of restric- tions at a much deeper level of granularity than what the available data would allow us to test for (city, province, or regional level). Therefore, we only present an analysis of the north Amer- ican region including all the US states and Canadian territories, for which OxCGRT provides the highest level of detail. Table 5 provides the results of panel data analysis on the impact of stringency and contain- ment levels for COVID-19 and climatic factors. To evaluate the impact of the stringency factor, we created the two variables high_stringency and high_containment that take values equal to 1