Exposure-lag Response of Air Temperature on COVID-19 Incidence in Italy

Abstract


Introduction
As of May 2021, more than 150 million cases of COVID-19 infections have been confirmed worldwide, and the global death toll has now exceeded 3 million people (World Health Organization, 2020).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. (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 COVID-19.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 recently 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 (Dong et al 2021).Studies in Italy are no exception with reports of both positive (Passerini et al., 2020;Zoran et al., 2020) and negative (Pirouz et al., 2020;De Angelis et al 2021) associations.Several studies have also attempted to allow for non-linear and delayed effects of climatic predictors, though until very recently (Nottemyer and Sera 2021) these have been limited to select regions (Runkle et al., 2020), relatively narrow temperature ranges (Passerini et al., 2020;Bashir 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 significant methodological flaws in the majority of papers studied (Dong et al 2021).This is perhaps unsurprising given the urgency for published works.In addition to the limitations highlighted above, concerns highlighted by Dong et al (2021) included: low spatial resolution, failure to account for autocorrelation, failure to account for various categories of confounding variables (e.g.meteorological, policy, mobility, seasonal trends) and failure to account for lag effects with regards to confounders.
Here we attempt to address these concerns by performing an analysis of high spatial and temporal resolution using over a full year's worth of data.We first fit statistical models to select Italian cities, accounting for several categories of potential confounders (policy, mobility, meteorological and pollution variables).Estimates from these models are then synthesised using meta-analysis to yield pooled estimates of the exposure-lag response of air temperature on COVID-19 incidence.

Data sources
Nine Italian cities were included in this analysis, namely: Bologna, Brescia, Livorno, Milan, Modena, Parma, Prato, Rome and Trieste, selected based on the availability of daily city level data, between 01 Jan 2020 and 24 Apr 2021 inclusive.Daily confirmed COVID-19 cases for each city were obtained using the COVID-19 R package (Guidotti and Ardia, 2020).Daily country-level public policy measures were also obtained from the above source, including: (3) parks, (4) transit stations, (5) workplaces; and (6) residential.

Study design
We used 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 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 exposurelag response is not city-specific and has been deployed in a recent similar study (Nottmeyer and Sera 2021).

First stage analysis
Due to the large number of potential confounding variables, we performed a redundancy analysis using R 2 =0.8 as the cut-off, to remove variables that could be strongly predicted from those remaining (Harrell 2015).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) of the form: eq. 1 where   denotes the COVID-19 cases on day t and  0 the model intercept.  ,  denotes a crossbasis for temperature on day t considering lag l, constructed using natural splines with 3 prespecified degrees of freedom for both the exposure and lag bases (Nottmeyer and Sera 2021 found that larger values for the degrees of freedom incurred a loss of precision for the estimates) with  1 the corresponding vector of coefficients.We considered lag range 0 through 10 days inclusive as in recent studies (Lauer et al., 2020;Nottmeyer and Sera 2021).dow denotes the day of week fitted as dummy variables, with   the corresponding vector of coefficients.
denotes additional predictors fitted as penalised smoothing splines.Of these smooth terms, time, t was fitted as a cubic regression spline with 9 degrees of freedom to control for seasonal and long term trends, whilst the remaining terms were fitted as penalized cross-basis splines (Gasparrini et al. 2017) in order to account for lag effects in these predictors.Specifically, grocery, work, park, humidity, PM10, pressure, wind-speed and wind-gust were fitted using penalized cubic regression splines with 3 degrees of freedom for both the exposure and lag dimensions (total = 9 degrees of freedom for each term).Transport closing and internal movement restrictions were fitted using dummy coding for the exposure dimension and natural splines for the lag dimension, with 3 pre-specified degrees of freedom.To avoid overfitting and potential multicollinearity, additional penalties were added to the null space of the bases for all smooth terms, allowing them to shrink to zero (Marra and Wood 2011).However, we did not shrink the cross-basis for air temperature since this was the predominant term of interest in our model (Harrell 2015).Heteroscedasticity and Autocorrelation Consistent (HAC) estimation of the covariance matrix was performed using the sandwich estimator (Zeileis 2004(Zeileis , 2006)).Parameters were estimated using quasi-restricted maximum likelihood (REML) (McCullagh and Nelder, 1989).We performed a sensitivity analysis by altering the lag period to 0-5 and 0-15 days respectively as in Nottmeyer and Sera (2021).

Second stage analysis
For each of the city-specific GAM's, the cross-basis   ,  was reduced to simpler sets of onedimensional coefficients and corresponding covariances 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 computed the Relative Risk (RR) (±95 % confidence intervals) and cumulative RR (RRcum) to quantify the direction and magnitude of the exposure-and lag-response respectively.For the exposure response, predictions were centred at the air temperature which yielded the minimum RRcum (Nottmeyer and Sera 2021).

Results
Across all city-specific time series', daily median air temperature ranged from -2.2 to 31.2 °C.
The air temperature yielding minimum RRcum was 4.7 °C; all exposure response functions were therefore centred at this point.
The exposure response curves varied according to the lag period (fig 1).At early lags, the RR increased from 1 at 4.7°C (reference level) to peak at 19.8 °C, before declining at higher temperatures.At lag 6 days, the RR increased linearly from the reference level as temperatures increased.At late lags, the RR largely remained close to 1 with increasing temperatures, though there was a decline below 1 at high temperatures by lag 10 days.The RR remained close to 1 for decreases in temperature compared below the reference level, though there is the possibility of a slight increase by the end of the lag period.
The lag-specific effects were summed to yield the cumulative exposure response curve (fig 2).
For temperatures above the reference level, the cumulative exposure response function resembled a bell-shaped curve; with RRcum increasing to peak at 19.8 °C before declining to RRcum ~ 1 at 30 °C.For temperatures below the reference level, RRcum increased slightly.
The exposure-specific lag-response curves are shown in fig 3. Small increments either below or above the reference temperature produced little effect and the RR's were constant across the lag period.For larger increases in temperature, the RR increased to some peak before declining at the end of the lag period, with the suggestion that the peak occurs later for larger increments in temperature compared to the reference level.Examining the cumulative effects plot (fig 4), the RRcum appears to monotonically increase over the lag period for 14 and 21 °C, though there is the suggestion of a plateauing effect in the latter.For air temperature = 28 °C, the RRcum peaks at around lag=8 days before declining to lag=10 days.
Shortening the modelled lag period from 0-10 to 0-5 days resulted in a sigmoidal shaped cumulative exposure response, with air temperature corresponding to maximal risk somewhat higher (approx.25 °C) (Supplementary information fig 1).However, the confidence intervals do not rule out a declining RRcum at higher temperatures (as in the main analysis) and there appeared to be increasing uncertainty at high temperatures with widening confidence intervals.
However, the maximal RRcum is similar in magnitude to the main analysis.Lengthening the modelled lag period from 0-10 to 0-15 days did not change the shape of the cumulative exposure response curve (Supplementary information fig 2).However, it appeared somewhat flatter, with a slightly lower maximal RRcum = ~2.However, this maximal risk occurred at approximately the air temperature as in the main analysis.There also appeared to be a loss of precision incurred by extending the lag period, with wider confidence intervals than in the corresponding main analysis.

Discussion
In this contribution, we quantified the exposure-lag response of air temperature on COVID-19 incidence in Italy.We fitted a series of city-specific statistical models to extended time series accounting for numerous confounders, before pooling the effects using meta-analysis.We found a non-linear exposure-lag response with minimal and maximal risk at 4.9 and 19.8 °C respectively.
Our study supports a growing body of evidence of a non-linear association between herein is in agreement this study suggesting "lower" and "higher" temperatures might reduce COVID-19 transmission.However, it appeared as though the peak risk was approximately 6 °C higher in our study (19.8 °C) compared to theirs (11.9 °C).Moreover, the maximal RRcum reported herein (2.39 [95% CI: 1.13; 2.94]) is also somewhat larger than that of the above study 1.62 [95% CI: 1.44; 1.81).The maximum RR for our exposure-specific lag response curves falls between approximately 3 and 6 °C and is consistent with Nottmyer and Sera (2021).
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 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 (Zhao et al., 2020;Lowen et al., 2007), whilst high temperature and low humidity promote the accumulation of aerosol particles, increasing the likelihood of inhalation (Nottemeyer and Sera 2021).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, so caution is warranted in drawing definite conclusions (Kratzel et al., 2020).
Our study has several strengths which we believe make it a valuable contribution to the literature, particularly in the light of criticism raised by Dong et al (2021).First, we employ the DLNM framework to simultaneously assess non-linear and delayed dependencies in the exposure and lag dimensions (Gasparini et al 2010).Second, we used city-level daily data for COVID-19 incidence and a number of confounders.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 a large number of confounders including meteorological, pollution, policy and mobility and allow for delayed dependencies in these effects.Fifth, the precision of our estimates are robust to heteroskedasticity and autocorrelation.Sixth, we controlled for overfitting using shrinkage estimators.In sum, we have addressed all the main concerns highlighted in a recent review of the relationship between air temperature and COVID-19 incidence (Dong et al 2021).
However, our work is not without limitations.The first concerns the accuracy of the estimated daily counts of COVID-19.Notwithstanding the rate of testing, it is postulated that estimates are biased low due to undocumented infections (Li et al., 2020) or other reporting errors (De Natale et al., 2020).Secondly, our study is limited to the lag effects of median daily temperature.A recent study (Babin, 2020) reported that extremes of daily temperature might actually be a better predictor of COVID-19 incidence; this would be a worthy avenue for future research.Third, our analysis was based on a smaller number of cities than similar work (Shi et 2020;Nottemyer and Sera 2021); this was a direct consequence of the availability of city level data for confounders which we considered a key feature of our analysis.Although metaanalysing 9 cities in this framework is valid (the example in Gasparini and Armstrong 2013 used 10 regions), further work might consider a larger sample size of cities to yield more precise effect estimates.Forth, although meteorological and pollution data were city level, policy and mobility data were country level which may incur some aggregation bias.
In conclusion, our results support previous work of a non-linear exposure-lag response of air temperature on COVID-19 incidence, with "low" and "high" temperatures less favourable to outdoor transmission.However, our finding of maximal risk at 19.8 °C is somewhat higher than similar recent work.If our findings are correct, this implies that only at higher temperatures (>20 °C) might the risk of COVID-119 infection start to reduce.Due to this uncertainty, this work underscores the need for facemasks and social distancing even in warm temperatures.

Figure 1 .
Figure 1.Lag-specific effects of the exposure-response association.The effects are centred at the reference level (solid diamond = 4.7 °C).Grey shaded area represents 95% confidence intervals.

Figure 2 .
Figure 2. Overall cumulative exposure-response association.Cumulative RR is computed by summing the log RR over the lag period (0 through 10 days inclusive).The effects are centred at the reference level (solid diamond = 4.7 °C).Grey shaded area represents 95% confidence intervals.

Figure 3 .
Figure 3. Exposure-specific effects of the lag-response association.The effects are centred at lag=0 and temperature = 4.7 °C.Grey shaded area represents 95% confidence intervals.

Figure 4 .
Figure 4. Overall exposure-specific lag-response association.The effects are centred at lag=0 days and temperature = 4.7 °C.Cumulative RR is computed by summing the log RR over the lag period (0 through 10 days inclusive).Dashed lines represent 95% confidence intervals.
temperature and COVID-19 incidence (Runkle et al 2020; Shi et al 2020; Gao et al 2021; Nottmyer and Sera 2021; Yuan et al 2021).Nottmyer and Sera (2021) performed a similar study on cities in England.The shape of the non-linear exposure-response relationship reported

Figures Figure 1
Figures

Figure 3 Exposure
Figure 3