key: cord-1045749-fnk92wpp
authors: Fong, Fang Chyi; Smith, Daniel Robert
title: Exposure-lag response of air temperature on COVID-19 incidence in twelve Italian cities: A meta-analysis()()
date: 2022-03-16
journal: Environ Res
DOI: 10.1016/j.envres.2022.113099
sha: f45c483051d2eaf3682441381baf4bb207b18108
doc_id: 1045749
cord_uid: fnk92wpp

The exposure-lag response of air temperature on daily COVID-19 incidence is unclear and there have been concerns regarding the robustness of previous studies. Here we present an analysis of high spatial and temporal resolution using the distributed lag non-linear modelling (DLNM) framework. Utilising nearly two years’ worth of data, we fit statistical models to twelve Italian cities to quantify the delayed effect of air temperature on daily COVID-19 incidence, accounting for several categories of potential confounders (meteorological, air quality and non-pharmaceutical interventions). Coefficients and covariance matrices for the temperature term were then synthesised using random effects meta-analysis to yield pooled estimates of the exposure-lag response with effects presented as the relative risk (RR) and cumulative RR (RR(cum)). The cumulative exposure response curve was non-linear, with peak risk at 15.1 °C and declining risk at progressively lower and higher temperatures. The lowest RR(cum) at 0.2 °C is 0.72 [0.56,0.91] times that of the highest risk. Due to this non-linearity, the shape of the lag response curve necessarily varied by temperature. This work suggests that on a given day, air temperature approximately 15 °C maximises the incidence of COVID-19, with the effects distributed in the subsequent ten days or more.

As of October 2021, more than 240 million cases of COVID-19 infections have been confirmed worldwide, and the global death toll has now exceeded 4 million people (World Health Organization, 2021) . Aside from the direct impact of COVID-19, healthcare in other fields unrelated to infectious disease have also suffered from its knock-on effects. Woolf et al. (Woolf et al., 2020) reported a higher number of deaths in the United States than expected when comparing data to previous years and estimated that 35 % of the excess mortality is due to conditions not directly attributed to . Though difficulty in seeking medical treatment due to movement restriction measures and reluctance to risk exposure in healthcare settings are speculated reasons for the excess mortality, an overwhelmed healthcare system likely played a significant role, especially in the formative stages of the outbreak.

Italy was among the worst affected countries near the start of the pandemic, with daily cases peaking at the end of March 2020. Italy remained within the top ten countries with the highest cumulative number of confirmed cases by June 2020, at more than 240,000 cases, with a recorded case fatality rate exceeding 10% (Worldometer, 2020) . Over a year since the start of the pandemic, widespread lockdown was reimposed due to transmission rates rising once more (Day, 2021) .

Since the beginning of the outbreak, a number of studies have investigated the relationship between air temperature and COVID-19 incidence, though the evidence is inconclusive . Studies in Italy are no exception, with reports of both positive (Passerini et al., 2020 , Zoran et al., 2020 and negative (De Angelis et al., 2021 , Pirouz et al., 2020 associations. Several studies have also attempted to allow for non-linear and delayed effects of meteorological predictors, though until very recently (Nottmeyer and Sera, 2021) these have been limited to select regions (Runkle et al., 2020) , relatively narrow temperature ranges (Bashir et al., 2020 , Passerini et al., 2020 , Tosepu et al., 2020 , relatively short time series (Shi et al., 2020) and lag periods (Tobías and Molina, 2020) , all of which may induce significant biases. Indeed, a recent review of papers which have investigated the association between air temperature and COVID-19 incidence identified critical methodological flaws in the majority of papers studied . This is perhaps unsurprising given the urgency for published works. In addition to the limitations mentioned above, concerns highlighted by Dong et al. (Dong et al., 2021) included: low spatial resolution, failure to assess J o u r n a l P r e -p r o o f and / or account for autocorrelation, and failure to account for various categories of confounding variables (e.g., meteorological, air quality, non-pharmaceutical interventions (NPI)), including their potentially delayed effects.

Here we attempt to address these concerns by performing an analysis of high spatial and temporal resolution using almost two years' worth of daily data. We fit statistical models to select Italian cities to quantify the delayed effect of air temperature on daily COVID-19 incidence, accounting for confounding variables. Coefficients and covariance matrices for the temperature term are then synthesised using random effects meta-analysis to yield pooled estimates of the exposure-lag response.

Twelve Italian cities were included in this analysis, namely: Bologna, Brescia, Firenze, Livorno, Milano, Modena, Napoli, Parma, Prato, Roma, Torino and Trieste, selected based on the availability of daily city level data, between 24 February 2020 and 30 December 2021 inclusive. Daily confirmed COVID-19 cases for each city, with respective populations at risk, were accessed using the COVID-19 R package (Guidotti and Ardia, 2020) , an interface to the COVID-19 Data Hub (https://covid19datahub.io/). For daily city-level data, confirmed cases and populations are sourced from Ministero della Salute (2021) and Istituto Nazionale di Statistica, Italia (Istituto Nazionale di Statistica, 2018) respectively. Country-level daily values of the stringency index (SI) are sourced from the Oxford Coronavirus Government Response Tracker (OxCGRT) project (Hale et al., 2021) . This is a composite metric assessing the number and strictness of government policies and is comprised of nine response indicators including travel bans, workplaces closures and school closures (https://www.bsg.ox.ac.uk/research/research-projects/covid-19-government-responsetracker).

Daily median city-level weather and air quality data were obtained from local meteorological stations (https://aqicn.org/data-platform/covid19/ (2022)) and included: (1) air temperature, (2) relative humidity, (3) air pressure, and (4) PM2.5. Daily city level residential mobility data was obtained from https://www.google.com/covid19/mobility/ (2020) . Finally, daily COVID-19 vaccination administration counts at regional level were obtained from J o u r n a l P r e -p r o o f https://raw.githubusercontent.com/italia/covid19-opendatavaccini/master/dati/somministrazioni-vaccini-latest.csv with corresponding populations from Istituto Nazionale di Statistica, Italia (Istituto Nazionale di Statistica, 2018).

We utilised the two-stage time series design (Berhane and Thomas, 2002) . In the first stage, a series of city-specific statistical models were fitted, with terms representing the exposure-lag response with respect to air temperature, as well as additional confounding variables. In the second stage, coefficients and covariance matrices from the first stage models were pooled using meta-analysis. This approach controls for city-level confounders, relaxes the strong assumption that the exposure-lag response is not city-specific, and has been deployed in a recent similar study (Nottmeyer and Sera, 2021) .

In order to model the temporal dependency between daily temperature and COVID-19

incidence, we adopted the Distributed Lag Non-linear Model (DLNM) framework (Gasparrini et al., 2010) . For each city, we fitted a Generalized Additive Model (GAM) (Wood, 2011 ) of the form:

Where denotes the number of COVID-19 cases on day t and 0 the model intercept. The logarithm of the population at risk (P) is fitted as a model offset. , denotes a cross-basis term for air temperature on day t considering lag l, constructed using natural splines with 1 the corresponding vector of coefficients. For the exposure dimension, internal knots were placed at the 27.5 th and 72.5 th percentiles, corresponding to 10.4 and 21.9 °C respectively (informed by minimising the AIC over a range of knot placement settings; see R code in the SI). For the lag basis, we pre-specified 3 degrees of freedom, with knots placed at the 33.3 th and 66.7 th percentiles (Nottmeyer and Sera, 2021) found that larger values for the degrees of freedom incurred a loss of precision for the estimates). We considered lag range 0 through 10 days inclusive as in recent studies (Lauer et al., 2020, Nottmeyer and Sera, 2021) . DOW denotes J o u r n a l P r e -p r o o f the day of week fitted as dummy variables, with the corresponding vector of coefficients.

denotes additional predictors fitted as penalised smoothing splines, with the degree of smoothing estimated from the data. Time, t, was fitted as a highly flexible penalised natural spline with 20 knots to allow for seasonal and long-term trends, as well as sudden changes such as testing policy. We fitted cross-basis terms for additional confounders using penalised natural splines with 3 degrees of freedom for the exposure and lag dimensions and lag effects up to 10 days. This included: (1) log(daily cases + 1) in order to allow for autocorrelation (Imai 2015, Nottmeyer and Sera, 2021) ; residential mobility, cumulative percentage of vaccines administered; stringency index, relative humidity, air pressure and PM2.5 (time series plots of these variables are provided in the SI). By fitting these as penalised terms, corresponding parameter estimates can shrink to simple linear effects to avoid over-parameterization (Benedetti and Abrahamowicz, 2004, Chen et al., 2016) .

Heteroscedasticity and Autocorrelation Consistent (HAC) estimation of the covariance matrix was performed using the sandwich estimator (Zeileis, 2004 , Zeileis, 2006 . Parameters were estimated using restricted maximum likelihood (REML). Simulated Quantile-Quantile ( Figure   SI .9) (Hartig, 2021) and partial autocorrelation ( Figure SI .10) plots are provided as supplementary information. We performed a sensitivity analysis by adjusting the parameterisation of the cross-basis term for temperature since this was the predominant term of interest. Accordingly, we assessed the impact of changing the lag period to 5 and 15 days, as well as increasing the degrees of freedom for the lag basis to 4 and 5 respectively.

For each of the city-specific GAM's, the cross-basis term for air temperature, , , was reduced to simpler sets of one-dimensional coefficients and corresponding covariance matrices for the exposure-and lag-dimensions respectively (Gasparrini and Armstrong, 2013) .

These were then pooled using random effects meta-analysis with REML estimation (Sera et al., 2019) . We assessed for heterogeneity using the multivariate Cochran's Q test (using prespecified α=0.05) and the I 2 . Statistic. The pooled exposure (i.e., air temperature) response summed over the lag period was presented as the cumulative relative risk (RRcum), centred on the temperature of the maximal RRcum. The pooled lag response at select temperatures (i.e., "exposure specific") was presented as the relative risk (RR) centred on the current day (i.e. lag J o u r n a l P r e -p r o o f day 0). We computed 95 % confidence intervals for both RR and RRcum to quantify uncertainty in these point estimates.

All statistical analyses were performed in R studio version 2021.09.0, using R version 4.1.2 (R Core Team, 2013) and relying heavily on the packages 'dlnm' (Gasparrini, 2011) , 'mgcv' (Wood, 2011) and 'mixmeta' (Sera et al., 2019) . R code and datasets are included in the supplementary material to fully reproduce the analysis. The lag response curves were centred at lag=0 days and temperature=15.1°C (Figure 2 ). The shape of the lag response necessarily varied by temperature due to the non-linear cumulative exposure response curve (fig 2) . The magnitude of the effect becomes greater for increasing or decreasing temperatures compared with the reference temperature. Accordingly the lag J o u r n a l P r e -p r o o f response curves are relatively flat at temperatures 10 and 20°C compared with 1 and 30° C (fig2). At 1°C (top left fig 2) the decline in RRcum is statistically significant (P<.05) from day0 to approx. lag 5.5 days with peak reduction at approx. 3 days lag. In contrast, For 30°C the reduction in RRcum is delayed with point estimates <1 only occurring beyond day lag 2 and so that significant between approx. 7 and 9.5 days lag. Aside from very high temperatures the curves retract to the null at the end of the modelled lag (10 days).

Changing the lag period of the cross-basis term for temperature from 10 to either 5 or 15 days in the sensitivity analysis incurred a loss of precision. For lag 5 days, the decline in RRcum above the reference temperature was less certain with confidence intervals that include the null and becoming wider with increasing temperature (Figures SI.11) . Upon increasing the lag to 15 days, the temperature of the maximum RRcum reduced slightly to approx. 11°C ( Figure   SI .12). Under this scenario point estimates hint of a possible upward trend at high temperatures, though again confidence intervals are too wide to conclude with any degree of certainty.

The shape and precision of the RRcum curve was relatively robust to increasing the degrees of freedom allocated to the lag basis from 3, to 4 and 5 respectively ( Figures SI.13 , SI.14).

J o u r n a l P r e -p r o o f 

In this contribution, we quantified the exposure-lag response of air temperature on daily COVID-19 incidence in Italy. We fitted city-specific statistical models to extended time series accounting for numerous confounders, before pooling the temperature effects using metaanalysis, which was justified by low heterogeneity. We show that the cumulative exposure response curve is non-linear with peak risk at 15.1°C, and declining risk at progressively lower and higher temperatures.

The decline in RRcum at temperature >12.9°C was slightly shallower and exhibited more uncertainty than the corresponding effect at temperatures <12.9°C. The air temperature yielding minimum RRcum was 0.2°C, i.e., the minimum temperature used in the pooled analysis.

Due to the non-linear cumulative exposure response curve, the shape of the lag response curves necessarily vary by temperature. RR falls below 1 within the lag period which is consistent with delayed effects of an increase or decrease in temperature (relative to reference) on a given day. Whilst the lag response for low temperatures declines and then returns to the null by day lag 10, this was not observed at high temperatures where the effects seemed further delayed with RR <1 at lag day 10 (though the effect is not statistically significant). Moreover, our sensitivity analysis suggests effects may be distributed beyond 10 days. It is possible that, at the population level, temperature changes on a given day may affect daily incidence up to and beyond 15 days hence, and further work might consider investigating longer lag periods.

Nottmeyer and Sera (Nottmeyer and Sera, 2021) performed a similar study on cities in England.

The shape of the non-linear exposure-response relationship reported herein shows a similar trend to their study (though note the different centering temperature), suggesting progressively lower and higher temperatures beyond that at the maximum RRcum, are associated with reduced COVID-19 incidence. However, it appeared as though the temperature at which peak risk occurred was approximately 3.2°C higher in our study (15.1 °C) compared to theirs (11.9 °C).

Moreover, the largest absolute RRcum reported herein (1.39 [95% CI: 1.10; 1.79]) is also slightly smaller than that of the above study 1.62 [95% CI: 1.44; 1.81]). The temperature for which risk is minimised was 0.2°C in the current study compared with 21.8°C in their study (though confidence intervals widened at higher temperatures in their study).

Experiments on SARS-CoV-2 have demonstrated a temperature-dependent effect on the stability of virus particles on smooth surfaces (Chin et al., 2020) . Since SARS-CoV-2 is shown to be capable of surviving in aerosols for at least 3 hours, airborne transmission is highly J o u r n a l P r e -p r o o f plausible (van Doremalen et al., 2020) . Aerosol dynamics, and subsequently disease transmission, appear to be affected by air temperature and humidity. Low humidity and low temperature result in smaller infected droplet nuclei that can travel further, remain airborne for a longer duration, and cause delayed virus inactivation (Chen, 2020 , Moriyama et al., 2020 .

Experimental evidence for the effects of climate on transmission dynamics using influenza virus in a guinea pig model lends credence to this hypothesis (Lowen et al., 2007) . Low temperature and high humidity favour contact transmission (Lowen et al., 2007 , whilst high temperature and low humidity promote the accumulation of aerosol particles, increasing the likelihood of inhalation (Moriyama et al., 2020) . As highlighted in the above reference, this ambiguity makes non-linear associations plausible. However, not all experimental studies agree that stability of SARS-CoV-2 is greatly affected by temperature (Kratzel et al., 2020) . More recent evidence showed that outdoor airborne transmission is unlikely in non-crowded areas (Belosi et al., 2021 , Chirizzi et al., 2021 , Pivato et al., 2021 .

Though outdoor air temperature likely has limited direct impact on transmission via virus stability in non-crowded outdoor settings, host factors, such as human behaviour influenced by meteorological variables, may underlie the observed trend. Several studies have suggested an association of human mobility with warmer temperatures (Shao et al., 2021 .

An exploratory study in Italy found that mobility correlated with increased spread of COVID-19 (Carta et al., 2021) . Although we attempt to control for this confounder by including city level residential mobility, omitted variable bias may have occurred if other mobility indices are confounded with temperature. Furthermore, the stringency index we included is only available on a national level. However, restrictions enforced by the Italian government varied across regions depending on the evolving situation of the pandemic. As such, we acknowledge this as a major limitation of our analysis, and, increased human mobility as influenced by temperature remains one plausible explanation for our results.

Our study has several strengths which we believe make it a valuable contribution to the literature, particularly in the light of criticisms raised by Dong et al. (Dong et al., 2021) . First, we employ the DLNM framework to simultaneously assess non-linear and delayed dependencies of air temperature in the exposure and lag dimensions (Gasparrini et al., 2010) .

Second, we used city-level daily data for COVID-19 incidence and several confounders providing high spatial resolution. Third, we used the powerful two-stage design; this allows for both the background risk and the exposure-lag response functions to vary by city which would otherwise incur ecologic bias. Forth, we account for many confounders, including J o u r n a l P r e -p r o o f meteorological, air quality and NPI variables, and allow for delayed dependencies in these effects. Fifth, analysing data from a single country, Italy, which has coordinated nationwide COVID-19 interventions, as compared to analysis using data from different countries, allows for better control over confounding variables. This is especially significant, because confounders, such as government interventions, likely have a greater influence over an immunologically naïve population more than meteorological factors (Baker et al., 2020) . For the later period of the study, we have also accounted for vaccination efforts by including the cumulative percentage of vaccine doses administered in each region. Sixth, we have accounted for and assessed autocorrelation and model fit in each of the city level models. Finally, we controlled for potential overfitting using shrinkage estimators for confounding variables. In sum, we have addressed the major concerns highlighted in a recent review of the relationship between air temperature and COVID-19 incidence .

Our work has several limitations. The first concerns the accuracy of the estimated daily counts of COVID-19. In addition to varying rates of testing, it is postulated that estimates are biased low due to undocumented infections or other reporting errors (De Natale et al., 2020) . Secondly, our study is limited to the lag effects of median daily temperature. It has been reported that extremes of daily temperature might be a better predictor of COVID-19

incidence (Babin, 2020) ; this would be a worthy avenue for future research. Third, our analysis was based on a smaller number of cities (albeit longer time series and more confounders considered) than similar work (Nottmeyer and Sera, 2021, Shi et al., 2020) ; this was a direct consequence of the availability of good quality city-level data for confounders which we considered a key feature of our analysis. Although meta-analysing 12 cities in this framework is valid (the example in Gasparini and Armstrong (Gasparrini and Armstrong, 2013) used 10 regions), further work might consider a larger sample size of cities, increasing statistical power and yielding more precise effect estimates. Forth, although meteorological, air quality and residential mobility data were city level, vaccine data was regional, whilst SI data was country level policy; this may incur some aggregation bias and presents a limitation for our analysis.

Fifth, we also had not considered the possible change in virus infectivity throughout the course of the study period.

In conclusion, our results support recent works of a non-linear cumulative exposureresponse of air temperature on COVID-19 incidence. Peak risk for daily incidence is 15.1°C, with the effects seen over the subsequent 10 days or more due to the delayed response. It should be stressed that risks presented herein are population level risks since the statistical units J o u r n a l P r e -p r o o f analysed are population groups. This information is particularly useful for policy makers assessing population level trends, though might not necessarily correspond to individual level risk due to the ecological fallacy. 

We declare no competing interests, or other interests that might be perceived to influence the results and/or discussion reported in this paper.

COVID-19 community mobility reports

Dati COVID-19 Italia

Air Quality Open Data Platform

Use of Weather Variables in SARS-CoV-2 Transmission Studies

Susceptible supply limits the role of climate in the early SARS-CoV-2 pandemic

Correlation between climate indicators and COVID-19 pandemic

On the concentration of SARS-CoV-2 in outdoor air and the interaction with pre-existing atmospheric particles

Using generalized additive models to reduce residual confounding

A two-stage model for multiple time series data of counts

The COVID-19 incidence in Italian regions correlates with low temperature, mobility and PM10 pollution but lethality only with low temperature

Effects of ambient temperature and humidity on droplet lifetime -A perspective of exhalation sneeze droplets with COVID-19 virus transmission

Too many covariates and too few cases? -a comparative study

Stability of SARS-CoV-2 in different environmental conditions

SARS-CoV-2 concentrations and virus-laden aerosol size distributions in outdoor air in north and south of Italy

COVID-19 incidence and mortality in Lombardy, Italy: An ecological study on the role of air pollution, meteorological factors, demographic and socioeconomic variables

The COVID-19 Infection in Italy: A Statistical Study of an Abnormally Severe Disease

Data-related and methodological obstacles to determining associations between temperature and COVID-19 transmission

Distributed Lag Linear and Non-Linear Models in R: The Package dlnm

Reducing and meta-analysing estimates from distributed lag non-linear models

Distributed lag nonlinear models

COVID-19 Data Hub

A global panel database of pandemic policies (Oxford COVID-19 Government Response Tracker)

DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models

DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models

Population and household: Data and microdata

Temperature-dependent surface stability of SARS-CoV-2

The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)

Influenza virus transmission is dependent on relative humidity and temperature

Seasonality of Respiratory Viral Infections

Influence of temperature, and of relative and absolute humidity on COVID-19 incidence in England -A multi-city time-series study

A Preliminary Investigation on the Statistical Correlations between SARS-CoV-2 Spread and Local Meteorology

Development of an Assessment Method for Investigating the Impact of Climate and Urban Parameters in Confirmed Cases of COVID-19: A New Challenge in Sustainable Development

Short-term effects of specific humidity and temperature on COVID-19 morbidity in select US cities

An extended mixed-effects framework for meta-analysis

Mediation by human mobility of the association between temperature and COVID-19 transmission rate

Is temperature reducing the transmission of COVID-19 ?

Aerosol and Surface Stability of SARS-CoV-2 as Compared with SARS-CoV-1

Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models

Excess Deaths From COVID-19 and Other Causes

WHO Coronavirus (COVID-19) Dashboard [Online

WORLDOMETER. 2020. Italy Coronavirus

Influence of outdoor thermal environment on clothing and activity of tourists and local people in a severely cold climate city

Econometric Computing with HC and HAC Covariance Matrix Estimators

COVID-19: Effects of Environmental Conditions on the Propagation of Respiratory Droplets

Assessing the relationship between surface levels of PM2.5 and PM10 particulate matter impact on COVID-19 in