Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has rapidly spread across the globe, traversing diverse climatic and environmental conditions. Sustained local transmission has occurred in most countries, leading to political, social and economic challenges and devastating loss of life. From the early phase of the pandemic, there has been speculation that weather conditions could modulate SARS-CoV-2 transmission patterns. The debate has been driven by analogy with existing seasonal endemic respiratory viral infections, such as influenza and other human coronaviruses, which tend to peak in the drier and colder winter months in temperate climates1. However, specific mechanisms behind this seasonality, in terms of host immunity and susceptibility, viral stability or weather-sensitive human behaviour are poorly understood2. Dynamic transmission modelling has shown that meteorological variables, such as temperature and humidity, are unlikely to have been a dominant transmission risk factor in the early stages of the COVID-19 pandemic, given high population susceptibility3,4. As SARS-CoV-2 is a new virus to humans, with <1 year of data available at the time of writing, ascertaining the potential for weather modulated transmission is challenging. Several studies have attempted such analyses. However, many such studies had methodological weaknesses and the results were at times conflicting5,6. Study findings for temperature resulted in either a positive7,8, negative9,10, non-linear11,12 or non-significant association13,14 with the COVID-19 response variable. For example, most studies did not control for key modulating factors, such as varying government restrictions, socio-economic indicators, population density or age structure15,16,17.

In this study, we overcome methodological issues of previous approaches by using a two-stage ecological modelling approach to examine the impact of meteorological variables on SARS-CoV-2 transmission by comparing cities located across the globe, while accounting for confounding of non-pharmaceutical interventions (NPIs) and city-level covariates. The study is based on an extensive dataset, collected by the Multi-Country Multi-City MCC Collaborative Research Network (https://mccstudy.lshtm.ac.uk/), consisting of time series of daily COVID-19 cases registered between 11 January and 28 April 2020 in 409 locations (cities or small regions) in 26 countries. In the first stage, we estimated the effective reproduction number (Re), in each city, over a city-specific time window early in the epidemic. We use a renewal equation-based approach that estimates latent infections and then map these infections to observed notifications via an incubation period, a report delay and a negative binomial observation model with a day of the week effect18. Focusing on the early phase of the pandemic allows us to minimise possible biases coming from factors impacting Re (in particular non-pharmaceutical interventions (NPIs)), which developed as the pandemic progressed. These include change of ascertainment methods and strategies, the implementation of strong NPIs (e.g. travel bans, school closures and lockdown), the appearance of new variants and ultimately vaccination campaigns. Also, in the first stage we define our exposure variables as mean values of meteorological variables (including daily mean temperature, relative and absolute humidity, solar radiation, wind speed and precipitation), for each city, over the early-phase time window, using the ERA5 fifth-generation European Centre for Medium-Range Weather Forecast atmospheric reanalysis of the global climate19. In a second ‘cross-sectional’ stage, we estimate the association of city-level Re, calculated for the city-specific window (allowing for standard errors), with each meteorological variable, controlling for confounding by total population, population density, gross domestic product (GDP) per capita, percentage of population >65 years, pollution levels (i.e. particulate matter, PM2.5) and the lagged Oxford COVID-19 Government Response Tracker (OxCGRT) Government Response Index at the endpoint of the selected time window (lagged by 10 days), allowing for the two-level (cities and countries) structure of the data using a multilevel meta-regression model20 (see ‘Methods’ for further details). We believe the data used and the analysis performed in this study improves upon previous approaches. Specifically, the fine spatial scale of the city-level data and the methodological design, accounting for confounding of NPIs and city-level covariates, allows us to accurately quantify the relationship between meteorological variables and Re.

Results

Descriptive analysis of meteorological variables and R e

The bivariate distribution of mean temperature and the effective reproduction number (Re) across the 409 study cities is shown in Fig. 1, and the characteristics of the 26 countries are reported in Table 1. The mean effective reproduction number (Re) across all cities was 1.4, ranging from 0.7 to 2.1, with all but ten cities experiencing an epidemic curve with a reproduction number >1. Mean temperatures over the observation period (between January and April 2020) reflect the late winter/early spring in 381 cities situated in the northern hemisphere and the summer/early autumn seasons in 28 cities in the southern hemisphere. Of the 136 cities classified as having high Re values, 35 cities experienced low temperatures, 64 medium temperatures and 34 high temperatures (Fig. 1). When visualising the unadjusted association of Re with mean temperature, relative humidity (RH), absolute humidity (AH), solar radiation at the surface and stratified by climate zone, we found no clear pattern (Fig. 2).

Fig. 1: Effective reproduction number and mean temperature in the observation window for 409 cities.
figure 1

Bivariate plot of effective reproduction number (Re) and mean temperature (Ta) (°C) in the observation window for each of the 409 study cities. Dark purple circles represent cities with both high Re and high Ta, while pale purple circles show areas with both low Re and low Ta. Red circles represent cities with low Re and high Ta and blue circles depict areas with high Re and low Ta. The bar chart (bottom right) represents the number of cities in each category defined in the bivariate legend (bottom left).

Full size image
Table 1 Characteristics of the 26 countries included in the study.
Full size table
Fig. 2: Effective reproduction number vs key weather variables by climate zone.
figure 2

a Mean temperature (°C), b relative humidity (%), c absolute humidity (g/m3) and d solar surface radiation (J/m2) vs effective reproduction number (Re) by climate zone (409 cities). The area of the circles is proportional to the precision (inverse of the variance) of Re estimates.

Full size image

Associations between meteorological variables and R e

Using a two-stage meta-regression model, we quantified the influence of meteorological variables, including mean temperature, on Re between cities, while controlling for confounding factors including government interventions. After adjusting for the city-level characteristics (e.g. socio-economic and demographic factors) and the country’s OxCGRT Government Response Index, we found a modest, non-linear association of mean temperature and AH with Re (Table 2). Less strong evidence of association was found for RH, with no evidence of association for solar radiation, wind speed and precipitation (Table 2). The association between mean temperature and Re is non-linear, with Re initially rising to a peak at 10.2 °C, then falling to a trough at 20 °C, 0.087 (95% confidence interval (CI): 0.025; 0.148) lower than the peak, and finally rising again (Fig. 3). AH has a similar non-linear shape with a maximum difference of 0.061 (95% CI: 0.011; 0.111) between the peak at 6.6 g/m3 and the trough at 11 g/m3.

Table 2 Association between weather variables and Re.
Full size table
Fig. 3: Associations between weather variables, non-pharmaceutical interventions and the effective reproduction number.
figure 3

Non-linear associations between (a) mean temperature (°C), (b) relative humidity (%), (c) absolute humidity (g/m3) and (d) OxCGRT Government Response Index and predicted Re difference. Curves and their 95% confidence intervals show the predicted difference in Re with respect to a reference value set to the value at the trough of the curve for meteorological variables (a–c), or for the OxCGRT Government Response Index = 50 (d). Two-sided Wald test p values and adjusted curves with 95% confidence intervals were obtained from multivariable meta-regression multilevel models adjusted by population (log scale), population density (log scale), GDP (log scale), % population >65 years of age, PM2.5 (μg/m3, log scale) and OxCGRT Government Response Index, with cities nested within countries. The marginal distribution along the x-axis represents the observed data for that covariate.

Full size image

The effect of NPIs on R e

Although we calculated Re over a time window in which the OxCGRT Government Response Index, lagged by 10 days, had not yet reached 70, we included the value of the lagged OxCGRT Government Response Index at the end of the city-specific window in the model, to control for residual confounding. Despite being capped at 70, the OxCGRT Government Response Index had a strong association with the reproduction number (p < 0.0001) (Supplementary Table 4), explaining 13.8% of its variability (Fig. 3 and Supplementary Table 4, Models D1–D7) with an estimated reduction of Re equal to 0.285 (95% CI: 0.223; 0.347) when levels of the Government Response Index increase from 21 (5th percentile) to 66 (95th percentile). Mean temperature explained 2.4% and AH 2.0% of the variation in Re, and the five city-level characteristics explained 1.4% of the variability of the reproduction number (Supplementary Table 4, Models D1–D8).

Sensitivity analyses

We performed several sensitivity analyses to evaluate the robustness of the results considering alternative analytic or selection choices (see Supplementary Table 5). The main results are stable when including a country-level fixed effect in the meta-regression model, i.e. considering the only within-country variation of covariates and outcome. Restricting the analysis to cities with weaker interventions (OxCGRT Government Response Index <60) also gives similar results to the main analysis, apart from wind speed and precipitation also showing an association with Re. The association between mean temperature and the effective reproduction number holds across all the sensitivity analyses, apart from in tropical and southern hemisphere cities, when stratifying by tropical and non-tropical or northern and southern hemisphere regions. However, this may be explained by the small number of cities and the resulting low power in the tropical and southern hemisphere sub-group. The association between AH and the effective reproduction number is somewhat less robust with no association observed when excluding tropical or southern hemisphere cities, when excluding China and Brazil (countries with earlier and later observation periods) and when considering meteorological variables lagged by 10 days. Excluding the ten cities with Re < 1 shows a tendency of an increased Re for cities with low RH (p = 0.009) and a lower Re in cities with higher solar radiation at the surface (p = 0.047) (Supplementary Figure 5). We observed similar overall tendencies to our main results when we did not control for the OxCGRT Government Response Index in our model, although the effect of temperature and AH was enhanced (Supplementary Figure 6), and when considering meteorological variables lagged by 10 days (Supplementary Table 5). We found no evidence of an interaction between mean temperature and RH categorised in two levels (≤65% and >65%) using the median value of 65% as the category threshold (p = 0.428).

Discussion

We combined datasets of COVID-19 transmission with meteorological, demographic, socio-economic and intervention data for 409 cities in 26 countries across the world to estimate the association between meteorological factors and Re in the early phase of the COVID-19 pandemic. We found evidence of a modest non-monotonic association of outdoor mean temperature and AH with early-phaseRe, after controlling for potential confounders, including NPIs. Temperature explained 2.4% and AH 2.0% of the variation in Re, compared to 13.8% explained by the OxCGRT Government Response Index in the adjusted analysis. The associations of temperature and AH with Re were not independent; the high correlation between them precluded control of one for the other. Overall, there was little evidence for any change in the Re of COVID-19 associated with RH and no evidence for precipitation and wind speed.

Associations between temperature, humidity and SARS-CoV-2 transmission might be explained by three mechanisms. First, like other viruses with a lipid envelope, SARS-CoV-2 has been found to be sensitive to temperature, humidity and solar radiation under laboratory conditions21,22,23,24,25, which affects its ability to survive on surfaces and in aerosols. The droplet behaviour in aerosols changes with different temperature and humidity levels. Low RH promotes the accumulation of aerosol particles (since evaporation leaves behind floating droplet nuclei) increasing the likelihood to be inhaled26,27. Second, innate and adaptive immune response mechanisms have been shown to be modulated by seasonal changes. Lower levels of vitamin D, mediated by decreased ultraviolet B radiation exposure during winter might lead to impaired antiviral innate immune defences28,29,30. Breathing dry air can impair mucociliary clearance, reducing the ability of cilia cells to secrete mucus and remove viral particles (innate immune response)27,31. Interferon-stimulated genes, usually inducing an antiviral state as part of the innate immune response have been found to be impaired at low RH32. High temperatures have been shown to hinder virus-specific CD8+ T cell responses and antibody production (adaptive immune system)33. Third, human mobility, contact patterns and time spent indoors are affected by weather conditions34. Very hot and very cold conditions can lead to more time spent in enclosed spaces, which might increase the likelihood of SARS-CoV-2 transmission.

Findings from this study are only partly consistent with findings from other global studies using statistical approaches to investigate meteorological effects on COVID-19 transmission. Meyer et al.9 found that mean temperature had a modest negative association with COVID-19 incidence for temperatures above −15 °C based on a dataset of 100 countries, after controlling for surveillance capacity, time since first reported case, population density and median population age, whereas RH had a negative non-significant association with case incidence. Jüni et al.13 covering 144 geopolitical areas showed that temperature and humidity measures were not significantly associated with epidemic growth while significant associations were found for restrictions of mass gatherings, school closures and measures of social distancing, which are consistent with our findings of a stronger impact of the OxCGRT Government Response Index compared to climatic conditions. Wu et al.35 incorporating data from 166 countries found that a 1 °C increase in temperature and RH was associated with a 3% and 0.85% decrease in daily new cases, respectively, after controlling for wind speed, median population age, Global Health Index, Human Development Index and population density. Interestingly, non-linear associations between mean daily temperature and the instantaneous reproduction number (Rt) in the United States of America were found in a study by Rubin et al.12 with Rt decreasing to a minimum as temperatures rose to 11 °C, increasing between 11 and 20 °C and then declining again at temperatures >20 °C. The shape of the association may be influenced by the indirect effect of weather in varying the likelihood of social interactions, e.g. at higher temperatures people may congregate in public cities, such as beaches and festivals12, while colder temperatures could limit social activities, such as sporting events34. Runkle et al.11 concluded from varying longitudinal associations in four cities that specific humidity in the range of 6–9 g/kg (i.e. AH range of 7.6–11.4 g/m3) was a significant predictor of the COVID-19 growth rate, in line with our findings.

Unclear and inconsistent findings related to temperature and humidity may be due to methodological challenges and data limitations. Similar methodological challenges were highlighted when evaluating the association between air pollution and COVID-19 outbreaks36,37. The novelty of the virus, with less than a full annual cycle of data available in most places, makes it difficult to disentangle a seasonal signal or inter-annual trends from meteorological factors using time-series models38. Moreover, different interventions (e.g. restrictions of mass gatherings, international travel and school closures) adopted by countries at different times after the onset of local outbreaks potentially confound the association between weather variables and COVID-19 spread. These challenges have led us to consider an ecological approach where we compared the outbreak curve early in the epidemic, minimising the confounding effect of NPIs. Despite this, we found a strong association of the OxCGRT Government Response Index with Re, confirming the importance of interventions implemented early on in the epidemic in controlling COVID-19 dynamics39.

Comparing the early-phase outbreak curves in different countries is challenging given that countries have different case definitions, and early COVID-19 data only captured a small portion of cases, mainly hospitalised patients or individuals with severe symptoms. The estimated high proportion of asymptomatic cases compromises the use of COVID-19 case counts to estimate transmission dynamics40. We used an estimated response variable, i.e. the effective reproduction number, calculated accounting for reporting delays and other sources of uncertainty. The 20-day duration was chosen as a compromise between needing enough days for a more precise Re estimation, while, at the same time, limiting the window to provide more constant weather, case ascertainment and Re estimates within the window. A larger window would bias estimates in ways that cannot be readily adjusted for. Our meta-analysis approach accounts for the uncertainty in Re estimates, which in turn reduces the level of certainty in the results. Further, 20 days is ~4 generations of infections, which, under most reporting scenarios, is sufficient to be confident about estimates in the level of transmission. We assume that within the 20-day time window, the case definition is constant within a city or country and Re is not affected by differences in case definition between countries.

A clear strength of this study is the use of an extensive dataset of 409 cities, representing 44.8% of all cumulative reported COVID-19 cases registered by 31 May 2020 in the John Hopkins University Coronavirus Resource Center. Our analysis covers all major climate zones across the globe, ranging from temperate, continental to tropical and dry climate settings. Another strength is our flexible methodological and statistical approach. We used multilevel meta-analytic models that take into account uncertainty of the response variable, i.e. the effective reproduction number. The model allowed for possible non-linearity of the exposures, and we adjusted for a selection of key socio-economic and demographic factors, as well as using a random effect to account for the country- and city-level differences. We chose covariates based on their potential impact on viral transmission that might confound the examined association of weather and COVID-19 dynamics. Indeed, population density leads to higher contact rates, potentially increasing the likelihood of transmission41. The age structure of a population is relevant given that elderly people were found to be more susceptible to infection and more likely to experience clinical symptoms of COVID-19 compared to younger age groups, increasing the likelihood of seeking medical care and getting tested42. Moreover, differences in contact patterns among different age groups can further affect the number of COVID-19 cases in each age group42. Socio-economic indicators, such as GDP per capita, are important to consider as more deprived populations might be at higher risk of infection due to potential conditions of overcrowded accommodation, inability to work from home or limited access to medical care43. Also, among air quality factors, a positive association between PM2.5 and COVID-19 incidence and mortality has been reported44,45.

We investigated model uncertainties with several sensitivity analyses, e.g. excluding cities with R < 1, excluding China and Brazil, cities in the southern hemisphere, cities with a latitude lower than 45°, and cities with an OxCGRT Government Response Index of more than 60. Previous studies compared cities within a country or considered large geographical units13,35,46, which could lead to a limited exposure range with narrow temperature and humidity variability reported during winter seasons, or high measurement errors for meteorological variables defined over large geographical units. We considered small area/city units distributed among 26 countries worldwide, allowing a good exposure range and minimising the measurement error of the exposures.

Our study has several important limitations in addition to those already discussed. Cities in the northern hemisphere were overrepresented compared to southern hemisphere cities, which indicates that the findings might be more representative for cities in the global north. Our results need to be put into the context of complex uncertainties surrounding characteristics of the novel virus, such as incomplete knowledge on possible underlying mechanisms between weather conditions and the virus itself, the role of host immunity and the potential influence of weather-sensitive human behaviour, such as indoor crowding47. However, AH was found to demonstrate the strongest indoor-to-outdoor correlation, indicating that outdoor AH measures could reflect indoor conditions48,49. Data limitations regarding the novel virus, including varying accuracy of COVID-19 case numbers, limited data availability across cities and temporal constraints of an incomplete seasonal cycle of SARS-CoV-2 contribute to the limitations of this analysis.

Despite these limitations, the associations of weather with Re in this study suggests that such effects are likely to be small compared to other drivers of transmission. NPIs had a stronger impact on variation in transmission between cities than meteorological variables. We found no weather conditions in which transmission is impeded if precautions (social distancing, mask use, etc.) are not taken. These results support the statement that, to date, COVID-19 interventions are critical regardless of meteorological conditions.

Methods

Data

Data in this study were obtained from a well-established MCC Collaborative Research Network50. The current MCC network covers 750 locations (cities or regions) in 43 countries/regions. For this study, 26 countries provided a daily time series of COVID-19 cases for a total of 502 locations (cities or small regions). COVID-19 data were downloaded from a publicly available repository or obtained from health agencies (Supplementary Table 1) and data management was performed using Microsoft Excel 2019. The time series from 1 January 2020 to 31 May 2020 comprises 2,771,137 COVID-19 cases, representing 44.8% of the cumulative cases registered by 31 May 2020 in the Johns Hopkins database (https://coronavirus.jhu.edu/). Supplementary Table 1 shows the sources used for each country along with the definition of COVID-19 cases.

To limit potential confounding by NPIs and temporal variation in case ascertainment, we selected a 10–20-day window early in the epidemic, starting after at least ten confirmed cases had occurred in a 10-day period, to reduce noise introduced by imported cases. We excluded days for which the OxCGRT Government Response Index exceeded 70, accepting reduced windows down to 10 days in length. The OxCGRT collates publicly available information on 18 indicators about governments’ policy responses to the COVID-19 pandemic (https://www.bsg.ox.ac.uk/research/research-projects/coronavirus-government-response-tracker). These indicators are categorised as containment or closure policies (e.g. school and workplace closures, restrictions on gatherings and movement), economic policies (e.g. income support) or health policies (e.g. COVID-19 testing programmes). The OxCGRT Government Response Index aggregates these indicators into a single score between 0 and 100 and provides a measure of how many policies a government has enacted, and to what degree. We chose 70 as the maximum value of OxCGRT Government Response Index allowable as a compromise between limiting confounding by government interventions and including enough cities to enable estimation of the associations studied (see Supplementary text 1). Applying these conditions/restrictions reduced our dataset to 409 cities or small regions in 26 countries with an observation period between 11 January 2020 and 28 April 2020.

Most of the 409 cities are situated in the northern hemisphere (n = 381), and in temperate (n = 292) or continental (n = 65) climatic zones, with few cities located in tropical (n = 23) and dry (n = 29) climatic zones. The COVID-19 cases were observed in the early phase of the epidemics, ranging from the first week of February 2020 in China to mid-April 2020 in Brazil (Supplementary Figure 1). This early epidemic phase is characterised in many countries (except Uruguay and Singapore) by a reproduction number >1 (Table 1).

We estimated Re for infections in the time window of interest using EpiNow2 1.3.218. This R package implements a Bayesian latent variable method for estimating Re, where infections at time t are estimated using the sum of previous infections, weighted by an uncertain, gamma-distributed, generation time probability mass function, and multiplied by an estimate of Re51,52. Initial infections (prior to the first reported case) were estimated using a log-linear model with priors based on the observed growth in cases. Complete infection trajectories were then mapped to reported case counts by first convolving over the incubation period distribution and an estimated distribution representing the delay between symptom onset and case report (both assumed to be log-normal). Reporting noise was then added using a negative binomial observation model combined with a multiplicative day of the week effect (modelled using a simplex). Re was considered to be piecewise constant with a breakpoint 3 days into the time window. The Re estimate from the first 3 days of the window was discarded with the Re estimate from the remainder of the window used in all analyses. Each region was fitted independently using Markov chain Monte Carlo. Four chains were used with a warmup of 1000 samples and 4000 samples post warmup. Convergence was assessed using the R hat diagnostic53.

We used a gamma-distributed generation time with a mean of 3.6 days (standard deviation (SD) 0.7) and a SD of 3.1 days (SD 0.8)54,55. This generation time was slightly shorter than the consensus estimate reported by Ferretti et al.56, leading to our Re estimates and subsequent effect sizes being conservative. We used a log-normally distributed prior for the incubation period with a mean of 5.2 days (SD 1.1) and SD of 1.52 days (SD 1.1)57. The log-normal prior for the delay from symptom onset to case report was estimated globally using a subsampled Bayesian bootstrapping approach (with 100 subsamples each using 250 samples) using data from an international line list of cases. The resulting distribution had a mean of 6.4 days and a standard deviation (SD) of 17 days (or a log mean of 0.83 (SD 0.15) and a log SD 1.44 (SD 0.12). The subsampled bootstrap approach was used to incorporate both the temporal and spatial uncertainty in the reporting distribution as data specific to each setting and time point was not available.

To define our exposures, we considered the following time series from the ERA5 dataset: 2 m temperature, 2 m dewpoint temperature, surface solar radiation downwards, precipitation, and 10 m eastward (u) and northward (v) components of wind. These are published by the Copernicus Climate Change service on a regular latitude/longitude grid of 0.25° (~25 km × 25 km) in NetCDF format19. The hourly 2 m temperature, 2 m dewpoint temperature and surface solar radiation downwards were averaged for each day to derive daily mean temperature, dewpoint temperature and surface solar radiation. From dewpoint temperature and the corresponding temperature (T; °C) we obtained RH (%) using the R ‘humidity’ 0.1.5 package58. The following formula was used to calculate AH (g/m3), which represents the mass of water vapour in the air mixture59:

$${{{{{\mathrm{AH}}}}}}=(6.112\times {e}^{(17.67\times T)/(T+243.5)}\times 2.1674\times {{{{{\mathrm{RH}}}}}})/(273.15+T).$$

The hourly 10 m u and v components of wind were averaged for each day, and the daily average u and v components were used to compute the wind speed using the formula wind speed = sqrt(u2 + v2). Hourly precipitation data were summed to derive daily totals. The daily variables were calculated for each 25 km2 grid cell and assigned to a city if the city centroid fell within the grid cell.

Mean temperature (and other meteorological variables, Supplementary Table 3) observed during the city-specific time window reflect the late winter/early spring observation period in cities situated in the northern hemisphere and in temperate or continental climatic zones. We found a high correlation between mean temperature and AH (Supplementary Figure 2). Socio-economic and demographic characteristics were extracted from the OECD Regional and Metropolitan database60 and Worldcities database61 (Supplementary Table 2). We selected, a priori, the following set of confounders: total population, population density, % elderly population (>65 years) and GDP (per capita). Pollution data (PM2.5) for the observation period (10–20 days) was obtained from the Copernicus Atmosphere Monitoring Service global near-real-time service62,63,64. This product provides hourly modelled values of surface PM2.5 (μg/m3) at a 0.4 × 0.4 arc degrees grid cell resolution. The hourly time series were averaged over the observation period and linked to the city using the city centroid coordinates. Cities vary in terms of socio-demographic characteristics (Supplementary Table 3). The correlation between socio-demographic characteristics is shown in Supplementary Figure 3 and the correlation between meteorological variables, OxCGRT Government Response Index, day of the year and Re in Supplementary Figure 4. To account for differences in NPIs we used the OxCGRT Government Response Index65. In this study, we considered the 10 days lagged value of the OxCGRT Government Response Index, and for each city, we assigned the index on the last day of the specified window for each city39. Note, in our analysis, meteorological variables and socio-demographic covariates were collated and summarised at the city level, while the COVID-19 time series were defined at the smallest administrative level containing the city. We only included cities for which the COVID-19 time series were available for an area in which most of the population resided in that city. We, therefore, refer to our unit of analysis as a city.

Statistical analysis

For descriptive purposes, the following statistics (mean, standard deviation and range) were calculated for meteorological variables (mean temperature, AH, RH, surface solar radiation, wind speed, total precipitation) and covariates considered in this study (total population, population density, % elderly population (>65 years), GDP (per capita), PM2.5, OxCGRT Government Response Index). We also calculated the correlation (Pearson coefficient) among meteorological variables and among covariates.

The association between city-level covariates and climatic variables with the effective reproduction number were evaluated using multilevel meta-regression models with two levels (cities nested within countries)20 using the R ‘mixmeta’ 1.1.0 package. The inclusion of country as a random effect allowed the model to account for country differences (e.g. data reporting) with efficient use of the within- and between-country information. Moreover, the meta-regression models allowed us to consider the precision of the Re estimates, as estimated by its variance, giving less weight to more imprecise estimates for shorter time windows.

Firstly, we used two-level meta-regression models to evaluate the possible non-linear association between each meteorological variable and the reproduction number Re. We considered possible non-linearity in the association with Re using a natural spline parameterisation of the meteorological variables with a variable number of internal knots from 0 (linear term) to 5, placed at respective percentiles of the variable. We compare the models with different non-linear parameterisations of the meteorological variable using the Akaike Information Criteria (AIC), choosing models with the lowest AIC.

We fitted the following two-level random-effects meta-regression models with cities nested within countries and an increasing number of predictors; Model A with two random effects (cities and countries) and the intercepts, Model B including the OxCGRT Government Response Index, Model C considering also total population, density, GDP, % population older than 65 years, and PM2.5 (total population, density, GDP and PM2.5 were log-transformed due to the skewness of their distribution).

Then for each meteorological variable, we fitted two-level meta-regression models (D1–D6) with the meteorological variable as exposure and total population, density, GDP, % population older than 65 years, PM2.5 and the OxCGRT Government Response Index as covariates. We considered non-linearity in the association with Re using a natural spline parameterisation of the climatic variables with the number of internal knots as determined in the univariate analysis. The coefficients of the natural spline parameterisation of the meteorological variable were used to derive the plot of the association between the meteorological variable and Re in the 5–95th percentile of the meteorological variable distribution (Fig. 3 and Supplementary Figures 5 and 6). The coefficients of the natural spline parameterisation of the meteorological variable were also used to test the association between the meteorological variable and Re using the multivariate Wald test. All the tests were two-sided. Given the small number of pre-defined exposures variables, no adjustment was made for multiple comparisons.

We quantified heterogeneity between cities with standard measures of I2 66. These measures are estimated once from a meta-regression model without meta-predictors (Model A) and once from the meta-regression models (Models B, C and D1–D6) to assess the reduction in residual heterogeneity provided by the different set of predictors. For each model, we calculated the likelihood ratio test (RLR2) statistic67. RLR2 is calculated as 1 − exp(−2/409 × (log Likm − log Lik0), where log Likm is the log-likelihood of the model of interest and log Lik0 is the log-likelihood from a null model including only city and country random effect (i.e. Model A). For each meteorological variable, we calculated the difference in the likelihood ratio test R2 (RLR2) with respect to Model C (including random effects, OxCGRT Government Response Index and city-level covariates). For OxCGRT and city-level covariates, the RLR2 represents the reduction in RLR2 when dropping OxCGRT or city-level covariates from Model D1 with temperature and all other terms (i.e. random effects, OxCGRT and city-level covariates).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.