key: cord-355935-psnqrdo2
authors: Paez, Antonio; Lopez, Fernando A.; Menezes, Tatiane; Cavalcanti, Renata; Pitta, Maira Galdino da Rocha
title: A Spatio‐Temporal Analysis of the Environmental Correlates of COVID‐19 Incidence in Spain
date: 2020-06-08
journal: Geogr Anal
DOI: 10.1111/gean.12241
sha: 
doc_id: 355935
cord_uid: psnqrdo2

The novel SARS‐CoV2 has disrupted health systems and the economy, and public health interventions to slow its spread have been costly. How and when to ease restrictions to movement hinges in part on whether SARS‐CoV2 will display seasonality due to variations in temperature, humidity, and hours of sunshine. Here, we address this question by means of a spatio‐temporal analysis in Spain of the incidence of COVID‐19, the disease caused by the virus. Use of spatial Seemingly Unrelated Regressions (SUR) allows us to model the incidence of reported cases of the disease per 100,000 population as an interregional contagion process, in addition to a function of temperature, humidity, and sunshine. In the analysis we also control for GDP per capita, percentage of older adults in the population, population density, and presence of mass transit systems. The results support the hypothesis that incidence of the disease is lower at higher temperatures and higher levels of humidity. Sunshine, in contrast, displays a positive association with incidence of the disease. Our control variables also yield interesting insights. Higher incidence is associated with higher GDP per capita and presence of mass transit systems in the province; in contrast, population density and percentage of older adults display negative associations with incidence of COVID‐19.

From a small outbreak linked to a live animal market at the end of 2019 to a global pandemic in a matter of weeks, the SARS-CoV2 virus and COVID-19, the disease it causes, have threatened to overrun health systems around the world. In efforts to contain the spread of the disease, numerous governments in many regions and nations have either recommended or mandated social distancing measures, and as of this writing, millions of people in five continents shelter in place. There are encouraging signs that these measures have mitigated the spread of the virus (e.g., Lancastle 2020; Lewnard and Lo 2020; Wilder-Smith and Freedman 2020) . Even so, this has come at a high cost, and the consequences for all spheres of economic, social, and cultural life have been dire (e.g., Fernandes 2020; Luo and Tsang 2020) . As a result, there is a sense of urgency to anticipate the progression of the pandemic, in order to plan for progressive lifting of restrictions to movement and social contact (e.g., Kissler et al. 2020) . Needless to say, this has become a delicate, and politically charged, balancing act between public health and the economy (Gong et al. 2020) .

A salient question in the debate about how and when to ease social distancing measures is whether the virus will display seasonal variations. Existing research on similar pathogens suggests that the virus could be more stable and potentially easier to transmit in conditions of low temperature and low humidity. While this is encouraging, it is important to keep in mind that "not all seasonal respiratory viruses experience the same spatiotemporal patterns" (de Ángel Solá et al. 2020, sec. 4) . This urges caution when extrapolating from known viruses. The evidence in this respect is as yet inconclusive, and although easing restrictions as the weather warms may appear tempting, doing so prematurely could well undo weeks or possibly months of costly measures.

It is not surprising, given the stakes involved, that this issue has already triggered a lively debate. The current state of knowledge was well-summarized by the National Academy of Sciences, Engineering, and Medicine in the U.S. in a recent report (see National Academies of Sciences, Engineering and Medicine 2020). Engaged by the Office of the Executive for guidance on this matter, this organization concluded that: "[some] limited data support a potential waning of cases in warmer and more humid seasons, yet none are without major limitations… Additional studies as the SARS-CoV2 pandemic unfolds could shed more light on the effects of climate on transmission" (p. 6). To further complicate matters, much of the relevant work has yet to be peer-reviewed (see for instance the challenge of Harbert, Cunningham, and Tessler 2020;  to Araujo and Naimi 2020) .

With the above considerations in mind, our objective with this paper is to investigate the influence of environmental factors, concretely temperature, humidity, and sunshine, on the progression of the pandemic. We adopt a population health approach, and report results from a spatio-temporal model of the incidence of COVID-19 in the coterminous provinces in Spain, one of the countries hardest hit by the pandemic. We combine data on reported cases of the disease with metereological information, to create a spatio-temporal dataset covering a period of 30 days. We then join this dataset with provincial-level economic and demographic information to act as controls to our key environmental variables. These data are analyzed using a spatial Seemingly Unrelated Regressions (SUR) approach, which allows us to model incidence of COVID-19 as a contagion process.

The results of this research provide evidence of the effect of temperature, humidity, and sunshine on the incidence of COVID-19. The clearest result with respect to these variables is a lower incidence of COVID-19 at higher temperatures and levels of humidity, while the opposite happens with respect to hours of sunshine. Our control variables also provide some intriguing insights. Higher incidence is associated with higher GDP per capita and presence of mass transit systems in the province; in contrast, population density and percentage of older adults display negative associations with incidence of COVID-19. The results of this analysis provide support to the hypothesis of seasonality of the novel SARS-CoV2, and should be of interest to public health officials and policy makers grappling with the question of the trajectory of the pandemic.

Please note that this paper is prepared as a reproducible research document. The source R markdown document, as well as all data and code needed to reproduce/review/extend the analysis can be obtained from a public repository. 1

The global emergence of infectious diseases is mostly driven by environmental, ecological, and socio-economic factors (Jones et al. 2008 ). In the case of SARS-CoV2, the ecological factors include the interaction between humans and wildlife. Once transmission of a disease begins to happen between humans, socio-economic and environmental factors become increasingly important. As noted in the introduction, the focus of the paper is on environmental variables, concretely three related to meteorological conditions: temperature, humidity, and sunshine.

Much of what is known about the potential seasonality of SARS-CoV2 is a result of research on other pathogens. Earlier, diverse studies have shown the effect of temperature and humidity on the incidence of influenza (e.g., Mäkainen et al. 2009; Jaakkola et al. 2014; Kudo et al. 2019) . Jaakkola et al. (2014) , for example, found that a decrease of temperature and absolute humidity precedes the onset of symptoms of influenza A and B viruses by 3 days in places where the temperature is low. After the 2002-2004 outbreak of SARS, researchers also began to investigate the relationship between these factors and SARS-CoV (Casanova et al. 2010; Chan et al. 2011 ). Casanova et al. (2010 , for instance, used surrogates to find that virus inactivation was likely more rapid at higher temperatures; in terms of humidity, these researchers reported that survival of the virus was lower at moderate relative humidity levels. Chan et al. (2011) also found that viability of the virus that causes SARS is also lost at higher temperatures (> 38 • C) and relative humidity superior to 95%.

Whether results from laboratory experiments will hold when the virus circulates in the community remains uncertain. At a global scale, de Ángel Solá et al. (2020) see less risk from SARS-CoV2 in the Caribean Basin; on the other hand, Coelho et al. (2020) warn that at least during the exponential phase, expansion of the virus is not driven by climate. Similarly, whereas Araujo and Naimi (2020) argue that spread of SARS-CoV2 will likely be constrained by climate, Harbert, Cunningham, and Tessler (2020) remain unconvinced that spatial modelling can currently discriminate the distribution of the disease on the basis of climate, at least in the United States. Yao et al. (2020) , examined data from China and came to the conclusion that neither temperature nor ultraviolet indices had an association with transmission of COVID-19. This is despite previous research that has linked less exposure to UVB radiation to higher prevalence and severity of acute respiratory tract infections (Zittermann et al. 2016; Dabrowska-Leonik et al. 2018; Dinlen et al. 2016; Mathyssen et al. 2017; Esposito and Lelii 2015; Jat 2017; Moriyama, Hugentobler, and Iwasaki 2020) .

In addition to the environmental variables above, from a population health perspective it is also important to account for potential socio-economic and demographic confounders.

To account for population-level factors, the first variable that we consider is GDP per capita. Much has been written about globalization and the spread of infectious disease 2 . The growth in global connections has presented a challenge to spatial approaches in the initial stages of disease management, when the cause of a disease may still be unclear but the plane has already departed (Zhou and Coleman 2016) . In reference to the earlier outbreak of SARS, Van Wagner (2008) chronicles how Toronto's status as a global city turned out to be a vulnerability in this respect. In our case, we think of GDP per capita as a marker of a region's relative position in a network of global cities, and its potential to be further ahead in the trajectory of the pandemic. Furthermore, wealthier regions also tend to concentrate more activities that produce non-traded goods, including building and construction (Hallet 2002) . Therefore, it is possible that wealthier regions remain relatively more active even during a lockdown. On the other hand, we cannot discount the possibility that less wealthy regions have a higher proportion of workers in manual occupations who cannot telework, and therefore have more difficulties complying with shelterin-place orders.

Secondly, we consider percentage of older adults (over 65) in a region. Early evidence regarding COVID-19 suggests that the case rate mortality is higher at older ages (e.g., The Novel Coronavirus Pneumonia Emergency Response Epidemiology Team 2020). However, it is not clear that a relatively large population of older adults necessarily translates into higher transmission rates of the infection. The tool of choice in containing the spread of the disease has been social distancing. In this respect, the evidence from the field of transportation is that older adults tend to travel less frequently, for shorter distances, and have higher rates of immobility than most everyone, except the youngest members of the public (e.g., Roorda et al. 2010; Morency et al. 2011; Sikder and Pinjari 2012) . In other words, many older adults are, whether by preference or otherwise, already in a form of social isolation. Social distancing during the pandemic may actually reinforce that condition for them, as suggested by the analysis of age-structured social contact in India, China, and Italy of Singh and Adhikari (2020). Since the age-structured matrix of social contact in Spain is similar to Italy (see Prem, Cook, and Jit 2017) , our expectation is that populations with higher percentages of older adults will tend to have lower levels of social contact and hence of incidence.

Population density is also relevant since it directly affects the contact patterns and contact rates between individuals in a population (Hu, Nigmatulina, and Eckhoff 2013) . The evidence available suggests a positive relationship between the transmission of COVID-19 and population density (e.g., cumulative incidence in urban areas like NYC). For this reason, we anticipate a positive relationship between population density and the incidence of the disease.

The last variable that we consider as a control is the presence of mass transit systems in a province. Every province in Spain offers some form of public transportation, but only five provinces have higher order systems of mass mobility (e.g., metro or subway), namely Barcelona, Madrid, Sevilla, Valencia, and Bizkaia. Public transportation has been hypothesized to relate to the spread of contagious disease by some researchers using agent-based approaches and simulation (e.g., Perez and Dragicevic 2009; Wang et al. 2010) , and while we find scant evidence of a link in the literature, the idea is intuitively appealing. After all, unlike the isolation that a car offers to travellers, most mass transit system are cauldrons of social contact.

The first reported case of COVID-19 in Spain was on January 31st, 2020, when a German tourist in the Canary Islands tested positive for the virus. After this case, it was still a few weeks before the first domestic case was reported, on February 27th in Sevilla province (Andalusia). In a short period of time, as testing started to ramp up, it became clear that an outbreak was flaring. By March 11th the World Health Organization (WHO) declared COVID-19 officially a pandemic. This declaration marked a turning point for the public in Spain too. As of March 13th, the number of cases of COVID-19 reported in Spain was 4,473, with a majority of cases (1,990) concentrated in Madrid: these numbers were at the time the worst outbreak in Europe after Italy. In response to the situation, on March 13th the Spanish National Government declared a state of emergency, to go into effect on Saturday March 14th. As part of the state of emergency restrictions to most activities were imposed, with the exception of essential services (e.g., food, health) and some economic subsectors of industry and construction. A few days later, on March 17th, Spain closed its lands borders to allow entry only to returnee nationals and permanent residents. The lockdown was further hardened between March 30th and April 12th (including the Easter weekend of April 10th-12th) and during this period only essential activities were allowed. During this period, there was a dramatic reduction in overall mobility, both within provinces as between. 3

Our dataset includes information about the daily number of cases of COVID-19 reported in Spain at the provincial level (NUTS3 in Eurostat terminology) for the period between March 13th and April 11th, inclusive. For our purposes, we consider positive cases reported, but exclude symptomatic cases diagnosed by a doctor without a Polymerase Chain Reaction (PCR) test. The Spanish National Government publishes periodic updates at the regional level (NUTS2) and the information is also released at the provincial level as part of a collaborative project by geovoluntarios.com, 4 ProvidencialData19, 5 and ESRI España (see DATADISTA, 2020) . This information is compiled from several sources, mainly the official web pages of the Spanish regional goverments, as documented in Centro de Datos COVID-19. 6 We consider two sets of explanatory variables. The first one, and the focus of this research, are the three environmental variables, collected from official sources (i.g., AEMET, the state meteorology agency, and MAPA, the ministry of agriculture, fisheries, and food). The second set provides some relevant controls as discussed above, and are also collected from official sources, namely INE, the national statistics institute. Table 1 shows the descriptive statistics and the provenance of the data used in this research.

The spatial and temporal coverage of the data is as follows. Our dataset begins on March 13, which is the first date when every province had reported at least one case of COVID-19, and continues until April 11, for a period of 30 days. The unit of analysis is the province. Provinces in Spain are the equivalent of states, and are embedded in Autonomous Communities. As an example, Cataluña is an Autonomous Community and consists of four provinces, namely Barcelona, Gerona, Lerida, and Tarragona. The size of the provinces is relatively large, as seen in Table 1 . The smallest province is 1,978.12 km 2 (this is smaller than Rhode Island in the U.S.) and the largest province is 21,767.93 km 2 (slightly smaller than New Jersey in the U.S.). While this is a relatively large degree of spatial aggregation, reporting on COVID-19 is not yet consistent at smaller geographies, or cases are not reported at that level at all.

An important aspect of working with environmental data such as temperature, humidity, and hours of sunshine, is the incubation period of the disease. Lauer et al. (2020) report for the case of COVID-19 a median incubation period of 5.7 days (with a confidence interval between 4.9 to 7.8 days). The vast majority of cases (95%) develop symptoms between 2.6 days (CI, 2.1 to 3.7 days) and 12.5 days (CI, 8.2 to 17.7 days). For this reason, we judge it best to use lagged values of the environmental variables. We test different time lags as follows. We consider a lagged 8-day average, from date-minus-12 to date-minus-5 days (hereafter lag8). Secondly, we consider a lagged 11-day average, from date-minus-12 to date-minus-2 days (hereafter lag11). Finally, to account for the likely duration of incubation, we consider a lagged 11-day weighted average, from date-minus-12 to date-minus-2 days (hereafter lag11w). The weights for this variable are calculated using the parameters of the log-normal distribution reported by Lauer et al. (2020) , that is, a log-mean of 1.621 and a log-standard deviation of 0.418. With these weights, the environmental variables at date-minus-2 and date-minus-12 days are weighted as 0.041 and 0.009 respectively, whereas the environmental variables at date-minus-5 days are weighted as 0.194. These weights have the effect of changing the contribution of daily values to the lagged moving average. For instance, the temperature at date-minus-4-days is weighted more heavily than the temperature at date-minus-10-days, as a closer approximation of the conditions at the most likely time of contagion before the disease became manifest.

The Seemingly Unrelated Regression equations model (SUR hereafter) is a multivariate econometric model used in different fields when the structure of the data consists of cross-sections for different time periods. The basis of this approach is well-known since the initial works of Zellner (1962) , and has become a popular methodology included in several econometrics textbook (e.g., Greene 2000) . To our knowledge, Anselin (1988) was the first author to discuss SUR from a spatial perspective, in the context of spatio-temporal analysis. In his landmark text, Anselin discussed a model made of "an equation for each time period, which is estimated for a cross section of spatial units" (p. 141). From this milestone, a large body of research has developed to extend the classical SUR into a spatial framework (e.g., Rey and Montouri 1999; Lauridsen et al. 2010; Le Gallo and Dall'Erba 2006; López, Martínez-Ortiz, and Cegarra-Navarro 2017) . The classical SUR model without spatial effects (from here, SUR-SIM) is a stack of equations as follows:

where y t = (y 1t , … , y Nt ) is a N × 1 vector, and in our case y st is the incidence ratio in the province s (s = 1, …, N) the day t (t = 1, …, T); X t is a N × k t matrix of the k t independent variables, with associated vector of coefficients β t ,; β t = (β 1t , … , β Nt ) is a vector of coefficients and ϵ t = (ϵ 1t , … , ϵ Nt ) is the vector of residuals.

A key feature of the SUR model is the temporal dependence structure among the vectors of residuals, namely:

Note that this specification is very flexible, in that it allows changes in the coefficients β it in order to modulate the effect of the independent variables on y t . This flexibility can be reduced and it is posible to impose restrictions considering some β coefficients as being constant over time.

In this case, we can reformulate the coefficients expression β t = (β 1 , … , β r−1 , β r , β r+1 , … , β Nt ) to restrict the first r coefficients to be constant across equations. This is equivalent to specifying some effects to be invariant over time.

Equation (1) can be rewriten in compact form:

where Y is now a vector of dimension NT × 1, X is a block-diagonal matrix NT × K (with K = ∑ t k t ) and ε is an NT × 1 vector. Using the Kronecker product notation (⊗), the error matrix structure is expressed concisely as:

As is the case with cross-sectional data, it is possible to test the residuals of Model (3) for spatial autocorrelation, and several tests have been developed to test the null hypothesis of spatial independence (see López, Mur, and Angulo 2014) . When the null hypothesis is rejected, several alternative specifications have been proposed to include spatial effects (Anselin 1988) . In this paper we consider a spatial SUR model that incorporates a spatial lag of the dependent variable as an explanatory factor. Spatial analytical approaches have been used to understand contagion-difussion processes in the case of infectious disease in general (e.g., Cliff, Haggett, and Smallman-Raynor 1998) and the 2003-2004 SARS outbreak in particular (e.g., Meng et al. 2005; Cao et al. 2010 ). While we are mindful of the same caveat that the novel SARS-CoV2 may not follow the patterns of previous diseases, there is still evidence from the United States that COVID-19 displays spatial patterns that are consistent with a diffusion process (Desjardins, Hohl, and Delmelle, 2020) . For this reason, the spatial lag model is appropriate to model incidence of COVID-19 geographically, since it accounts for potential spatial patterns that result from a process of contagion, as explained next.

The stack expresion for the SUR model with a spatially lagged dependent variable (SUR-SLM) is as follows:

This specification assumes that outcome (y st ) at location s and time t is partially determined by the weighted average (Wy st ) of the outcome in neighboring provinces, with neighborhood defined by matrix W of spatial weights. In other words, the spatially lagged dependent variable represents a process of contagion, where the disease in neighboring provinces can spillover in a spatial way. The coefficients of the spatially lagged variable are estimated for each time period ρ t and identify the intensity and the sign of the contagion effect. It is possible test the null hypothesis of identical levels of spatial dependence (ρ i = ρ j , ∀i, j). The correspond Wald test is available in the R package spsur.

The SUR-SLM model can be estimated using maximum likelihood (López, Mur, and Angulo 2014) or instrumental variables (Mínguez, López, and Mur 2019) .

A note regarding the interpretation of the model is in order. It is well-known that coefficients in linear regression models are partial derivatives of the dependent variable with respect to the independent variables, and therefore directly give the marginal effects or rates of change. Spatially lagged models, however, are no longer linear. The intuition behind the non-linearity is that the spatial lag expands the information set to include information from neighbouring regions: in other words, the value of an explanatory variable in a spatial unit can have influence in other spatial units via the spatial lag. This makes interpretation of the coefficients less straightforward but also richer (Golgher and Voss 2016) . The results of LeSage and Pace (2009) for cross-sectional spatial lag models can be extended to the spatial SUR framework. Note that, according to Elhorst (2014) , the partial derivatives have the following interpretation: if an explanatory variable (X k ) in a particular province changes, not only the incidence rate in that province changes, also incidence rates in other provinces change via the contagion effect. Therefore, a change in X k in a particular province has a direct effect on that province, but also an indirect effect on neighbouring provinces. In this way, the ith diagonal element of the matrix of partial derivatives represents the (5)

direct effect on the ith province, whereas the non-diagonal elements of ith row of the matrix of partial derivatives represent the indirect effects on other provinces. In order to obtain a global indicator, the direct effect is calculated as the mean of the diagonal elements and captures the average change in incidence ratio caused by the change in X k . Likewise, a global indicator of the indirect effects is calculated as the mean of the non-diagonal elements. The total effect is the sum of direct and indirect effects. Finally, note that if ρ k = 0, the indirect effects are 0 and the direct effects are equal to β k t. Fig. 1 shows the geographical variation in the incidence of COVID-19 in Spain, as well as the temporal progression of the disease in weekly averages. Our analysis begins on March 13. Albeit still low, the highest incidence at this early date was in the provice of Álava, in the North of Spain. Álava is not the most populous province, with a population of only 331,549, but it has the highest GDP per capita of all provinces. Vitoria, its main city, is the capital of the Basque Country and has been the focus of efforts, along with San Sebastian and Bilbao, to develop a "Global Basque City" (Meijers, Hoekstra, and Aguado 2008) . The other early focus of the pandemic in Spain was Madrid, which is the most populous province in the country and has the second highest GDP per capita after Álava. The early outbreaks in these two provinces can be traced throughout the progression of the pandemic over time, although by the end of the period under study, other provinces had matched and even surpased their incidence rates, including Segovia and Soria to the north of Madrid, and Ciudad Real and Albacete to the south. Fig. 2 shows the distribution of the environmental variables in Spain. For ease of visualization we aggregate the provinces by Autonomous Community. Each box-and-whisker in the figure represents the distribution of the variable for a community over the 30-day period. In the plot, the communities have been sorted by latitude, so that Principado de Asturias is the northernmost community, and Andalucia the southernmost. As seen in the figure, there is a relatively wide range of values both within and between provinces over the 30-day period examined. The top panel of the figure shows the distribution of mean temperatures. The lowest mean temperature for a community on any given day was approximately 3 • C, and the highest about 20 • C, for a range of approximately 17 • . Likewise, there is a great deal of variability in humidity, as seen in the middle panel of the figure, where the lowest mean humidity for any community is approximately 48% and the highest is close to 100%. Finally, the bottom panel displays mean daily hours of sunshine in the community. This variable displays the most variability within communities over time, but remains relatively constant across communities. It is important to note that these are summaries by community, and the actual values of these variables for the provinces display somewhat more variability. Fig. 3 includes three maps that display the spatial variation of our control variables, namely GDP per capita, percentage of older adults in province, population density, and presence of mass transit systems. As seen there, GDP per capita is higher in Madrid and the northeast part of the country, mainly in Pais Vasco and Cataluña. Percentage of older adults is somewhat more checkered, with high values in Madrid and other provinces in the center-west part of the country, but also in some provinces in the north. Outside of provinces with major cities (e.g., Madrid; Bizkaia and Gipuzkoa in Pais Vasco; Pontevedra in Galicia), population density tends to be higher in provinces along the Mediterranean coast. The final panel in the figure shows the five provinces in the country that have higher order mass transit systems (e.g., metro). Fig. 4 shows the distribution of daily simple correlations of incidence of COVID-19 with the independent variables (with the exception of Transit, which is a categorical variable). These correlations are calculated after log-transforming all variables. As previously discussed, the environmental variables have a temporal lag and were calculated using different time windows.

It can be seen in the figure that temperature (in its three forms) has the highest simple correlation with incidence of COVID-19. After temperature, GDP per capita has the highest positive correlation with the dependent variable. The distribution of these correlations is also quite tight over the 30-day period of the study. Hours of sunshine tends to have a moderately high correlation with incidence of COVID-19, but the distribution of these correlations is more spread, and in some cases strays into negative values. A similar thing happens with humidity, which also tends to display more day to day variation in the correlation with the dependent variable. The percentage of older adults shows a relatively tight distribution of day-to-day correlations, and Autonomous Communities in Spain (sorted from north to south) Value is negative. Population density, in contrast, tends to be negative, but is relatively spread, and on some days, the simple correlation between density and incidence of COVID-19 is weakly positive. Overall for the period under examination, the pairwise correlations between these variables and incidence are significant at P < 0.05, with the exception of the three Sunshine_hours variables.

Correlation analysis in the preceding section provides some insights about the potential associations between incidence of COVID-19 and the various environmental and control variables.

In this section we estimate three spatial SUR models to test the differences between the various temporal lags and weighting schemes for the environmental variables. Accordingly, we define three models: Model 1, which is estimated using the lagged 8-day averages of the environmental variables (lag8); Model 2, which is estimated using the lagged 11-day averages of the environmental variables (lag11); and finally, Model 3, which is estimated using the lagged 11-day weighted averages of the environmental variables (lag11w). 

To implement the SUR approach, we must define a matrix of spatial weights W. In this case, we consider neighborhoods based on adjacency, based on the well-known queen criterion (two provinces are adjacent if they share a boundary or touch at a vertex). The analysis is of the coterminous provinces. 7 For estimation, we log-transform the dependent and quantitative independent variables. The only variable that is not transformed is the categorical variable for transit systems. A log-log transformation is appropriate to capture non-linear relationships between variables and provides a straightforward interpretation of the coefficients as percentage change. Furthermore, we introduce restrictions so that the coefficients of two of our control variables are constant over time, namely GDP per capita and percentage of older adults. We do not see an a priori reason to let those two variables vary across equations, and the correlation analysis in Fig. 4 also suggest little temporal variation. In contrast, we let the spatial autocorrelation parameter, as well as the parameters of the rest of the independent variables (including the constant) to vary over time. 8 This will be useful to detect whether there are behavioral adaptations at the population level over the course of the period examined. As an example of behavioral adaptations, the effect of density might weaken over time, in the measure that the effects of the lockdown are felt: at full compliance with the lockdown, with people practicing social avoidance, density might matter less than other factors.

After estimation, we compare the goodness of fit of the three SUR models. Fig. 5 shows the equation-level coefficient of determination R 2 , one for each time period/day. As well, the overall coefficient of determination for the system is reported for each model pooled-R 2 . The general trend is identical for the three models, with the goodness-of-fit improving over time and plateuing around a value of R 2 of 0.55. Model 1 (lag8) performs somewhat better in the first few equations/days, when the goodness-of-fit is relatively poor, and then again in the last few equations/ Variable Correlation days. Model 3 (lag11w), in contrast, does not perform well towards the end of the study period. The most balanced model in terms of equation-level goodness-of-fit appears to be Model 1 (lag8), and this impression is further supported by a slightly higher value of the pooled-R 2 . The analysis using a lagged moving average of the environmental variables is generally in line with the incubation period reported by Lauer et al. (2020) , although the results do not support the use of a weighted average. For the remainder of the paper, we will adopt Model 1 (lag8) as our best model. In the following section we discuss the results of the analysis in more depth. Table 2 presents a summary of the results of Model 1 (lag8). Recall that two coefficients were constrained and are estimated only for the first equation of the system, and thus are assumed to be constant over time. These are the coefficients corresponding to GDP per capita (P ≤ 0.10) and percentage of older adults (P ≤ 0.05). The sign of the coefficient for GDP per capita is positive, which indicates that wealthier regions tend to have a higher incidence of COVID-19. This is in line with the idea that the epidemic started earlier in wealthier places due to their connections to a globalized world. The sign of the coefficient for percentage of older adults, on the other hand, is negative. As previously discussed, the level of social contact of older adults even under normal circumstances tends to be lower than for younger people. As a consequence, places with larger populations of older adults appear to have a natural level of social distancing in place. It is important to note that this does not detract from evidence that older adults are more vulnerable individually and in institutional settings, where their case mortality rates are perhaps the highest of all age groups. Instead, this result indicates that their presence in the community at large tends to depress transmission of the virus.

Of the two other control variables, the coefficient of population density is significant at P ≤ 0.05 in 12, at P ≤ 0.10 in 3 equations, and not significant in 15. The coefficient for transit 

Notes: Significance: This is the number of coefficients with p-values as indicated. Non-constrained: coefficient was allowed to vary across equations.

Constrained: coefficient as constant across equations.

is significant at P ≤ 0.10 in 20 equations, and of those, significant at P ≤ 0.05 in 18 equations. The next four variables are environmental factors. The coefficient for humidity is significant at P ≤ 0.19 in 20 equations, and of those, significant at P ≤ 0.05 in 15 equations. Of the environmental variables, temperature is the only variable that has significant coefficients in every equation at P ≤ 0.05. Finally, sunshine has significant coefficients at P ≤ 0.05 in 24 equations. To better understand the results, we proceed to plot the coefficients in their temporal sequence. At this point it is worth recalling that the state of emergency went into effect on March 14. In the following figures, the periods of time indicated in yellow starting on March 14 correspond to the state of emergency, with only essential travel and selected industrial activities allowed; the period of time in orange was the stricter lockdown when only essential travel was allowed.

We begin our discussion with the evolution of the spatial autocorrelation coefficient (ρ) in Fig. 6 (left panel) . We notice that the magnitude of the spatial autocorrelation coefficient ρ declines over the period under analysis, and is not significant for some days. This is an interesting result: immediately prior to the declaration of the state of emergency, there appears to have been a strong inter-provincial contagion effect. Keeping in mind that the incubation period ranges between 2 and 12 days with a median of 5, it is reasonable to expect that the effect of the lockdown will be observed with some delay. Indeed, as seen in the figure, the autocorrelation coefficient remains high for several days, then begins to decline around March 23, and continues to weaken over time. At the end of the period under examination, the strength of this effect is much diminished and we would expect that under full compliance with strict lockdown conditions (meaning no inter-provincial mobility) the spatial contagion effect would be zero-as seems to be the case.

The intercept (right panel in Fig. 6) is indicative of the variation of the incidence of COVID-19, other things being equal. Here we see that at the incidence declines somewhat immediately after the state of emergency, only to begin increasing again over time. Then, the incidence declines again after the stricter lockdown and rebounds to a higher level by April 11. Figure 6 . Temporal evolution of the spatial autocorrelation coefficient (rho) and the intercept of the model; dots are the point estimates and vertical lines are 95% confidence intervals. In yellow is the period after the declaration of the state of emergency, and in orange is the period when only essential activities were allowed. Fig. 7 shows the temporal evolution of the coefficients for the two control variables that were not fixed over time, that is, log (Density) (left panel) and Transit (right panel).

In Section "Background" we had anticipated a positive sign for the coefficient of density, and indeed, at the beginning of the period the coefficient is positive, albeit not significant, and then remains mostly non-significant for the earlier part of the lockdown. We are somewhat surprised by the way this coefficient turns significant and negative in the later part of the lockdown, after April 1st. This effect, we surmise, is the result of risk compensation, a situation where people adapt their behavior according to the perceived level of risk, becoming more careful when the perceived risk is higher and viceversa (e.g., Noland 1995; Phillips, Fyhri, and Sagberg 2011; Richens, Imrie, and Copas 2000) . Consequently, residents in high density regions may perceive the risk of infection as being higher, and adapt their behavior accordingly-while the opposite may be true of residents in low density regions. The coefficient for Transit is positive, as expected, but with very wide confidence intervals, and in fact not significant in the earlier part of the period.

The evolution of the coefficients for the three environmental variables is shown in Fig. 8 . Despite a mostly positive simple correlation with incidence (see Fig. 4 ) once that we control for other factors, humidity has a negative association with incidence of . This is in line with the literature that describes the lower viability and transmission of different viruses at higher levels of humidity. The coefficients for temperature (topright panel) are consistently negative and this variable is, besides the intercept, the only one with significant coefficients in all equations. The range of variation of this coefficient during the period examined is approximately between −1 and −2, although it is important to recall that these values should not be interpreted directly as effects; more on this below. Finally the plot for the coefficients associated with hours of sunshine (bottom panel) is more ambiguous: prior to the lockdown, the coefficient was negative, but not significant. However, five days into the lockdown, the coefficient becomes significant and positive. This result stands in contrast to previous findings regarding influenza, where more hours of sunlight reduced the strength and duration of epidemic durations (Yu et al. 2013) . A difference with previous studies is the temporal scale of the analysis: where Yu et al. (2013) use monthly averages, we use daily data for a much shorter period of time. The positive sign of sunshine may well be another instance of behavioral adaptations, whereby compliance with lockdown orders weakens on sunny days. The preceding discussion helps to establish the inferential contributions of the analysis, and indicate which variables display significant statistical associations with incidence of COVID-19. The remaining question is, what are the implications.

As discussed in Section "Methods: The spatial SUR model" the effect of a variable is not clear from its coefficient alone, since a change to a variable in a province influences, via the contagion effect, its neighbors. For this reason, the appropriate way to estimate the effects is to calculate both the own effect and the effect due to contagion, or in other words the direct and indirect effects, respectively. The total effect is the sum of the two. A summary of the effects appears in Table 3 . All effects in the table are interpreted as percentage change in the incidence of COVID-19 as a consequence of a one percent change in the variable. The exception to this is Figure 8 . Temporal evolution of coefficient for the environmental variables; dots are the point estimates and vertical lines are 95% confidence intervals. In yellow is the period after the declaration of the state of emergency, and in orange is the period when only essential activities were allowed. Transit (which was not log-transformed). This variable instead represents the percentage change in incidence between provinces without and with mass transit systems. Two variables had temporally constrained coefficients. The estimated effect of GDP per capita is to increase the incidence of COVID-19 by 0.449% for each percent increase of this variable (in €1,000 s). In our view, this is a measure of inertia, as provinces with higher GDP per capita where among the first to see exponential growth in the pandemic. Percentage of older adults has a negative effect, and each percent increase in this variable is associated with a relatively small reduction of the incidence of approximately 0.67%.

The temporal variation of the effects for the rest of the variables is shown in Fig. 9 . The largest positive direct effect is Transit, and the largest direct negative effects are temperature and humidity. The direct effects of these variables are as follows: for each percent point increase in temperature, there is between a 1% and 2% reduction in the incidence of the disease. This effect is compounded via contagion, as seen in the central panel in the figure, and the indirect effect can further reduce the incidence by up to 0.75%. The effect of humidity is also to reduce the incidence: each percent point of increase in humidity is associated with a reduction of up to 2% in incidence. With the addition of the indirect effect, the total effect of a 1% increase in humidity is to reduce incidence by up to 3%. As seen in the figure, the indirect (i.e., contagion) effects are stronger at the before and at the beginning of the lockdown period. Nonetheless, by the end of the period under study, the indirect effects have weakened considerably. What do these effects mean? Under a situation of lockdown, inter-regional contagion is reduced, as expected, and the total effects of the variables tend towards their direct effects. In the first few days covered by our analysis the total effect of all variables is greater due to the spatial contagion effect. Contagion makes analysis and intervention more complex: the contagion effect essentially acts like a multiplier, whereby developments in each province spill over to their neighbors. Once the contagion effect has been tamed, each province can be "treated" independently from its neighbors.

It is important, before concluding our discussion, to highlight some limitations of this study. First, the analysis was conducted mostly under a situation of lockdown, and therefore, besides first days of the period examined, one must exercise caution when trying to extrapolate the findings to a situations without lockdown. Secondly, it is well known that there is in many countries substantial underreporting of cases of COVID-19 due to limited testing. In this case we do not suspect systematic provincial bias in reporting, and as long as the underreporting is consistent across units of analysis, we do not expect biased results; it is still important, however, to remain aware that the number of true cases is likely higher. Finally, we defined neighborhoods based on adjacency. It would be interesting to compare adjacency to other connectivity criteria, for instance based on domestic transportation instrastructure and services. We flag this as a matter for future research.

In this paper we presented a spatio-temporal analysis of incidence of COVID-19 in Spain. The analysis covers a 30-day period that begins immediately before the state of emergency was declared in the country. The focus of the research has been on the environmental correlates of incidence of the disease. There is a need for more empirical evidence, as policy makers, public health practitioners, and the public begin planning for the months ahead at this early stage of the pandemic.

Our results offer strong support for the hypothesis that incidence of COVID-19 at the population level is lower at higher temperatures and levels of humidity: the estimated effect is a reduction in the neighborhood of 3% in incidence per each 1% increase in temperature, and a 3% reduction in incidence per 1% increase in humidity under conditions of contagion. These reductions are lower when contagion has ceased (i.e., due to lockdown conditions). The question here seems to be whether these environmental variables can yield a bigger reduction of more cases, or a smaller reduction of fewer cases.

Our control variables also offer some interesting insights. In particular, there is evidence of behavioral adaptations during lockdown in the form of risk compensation (density) and compliance with the lockdown (sunshine). These results offer a cautionary tale with regards to the effectiveness of the lockdown in more dense areas, and also the implications for compliance with stay-at-home orders as the northern hemisphere moves towards summer and more hours of sunshine during the day. If risk compensation is a factor, then efforts should be made to reduce or eliminate risk compensation in less densely populated regions.

A key aspect of the analysis using spatial SUR models is that we were able to model incidence of COVID-19 as an interregional contagion process. Here, we find that the strength of the contagion effect was dramatically reduced by the lockdown.

Needless to say, the analysis presented here is at the level of population health. For this reason, the analysis does not make any claims with respect to the effect of ultraviolet light on the virus, but rather about transmission of the virus in the population. For example, the analysis does not imply that the virus moves less effectively in places where more people live in close proximity to each other, but rather that humans are more contagious when they feel safe in less dense regions. Similarly, more sunshine does not mean that the virus thrives, but rather that humans are more contagious to each other when their behavior adapts to this environmental condition.

Some directions for future research include investigating other modelling frameworks, such as geographically and temporally weighted regression and/or space-time conditional autoregressive models. In addition, the environmental variables examined here relate to meteorological conditions only, and did not include other environmental factors that may incide in the transmission of the virus, such as pollution. These other factors should be incorporated in future studies. (Bivand, Pebesma, and Gomez-Rubio 2013) , spsur (Angulo et al. 2020 ) tidyverse (Wickham et al. 2019) , units (Pebesma, Mailund, and Hiebert 2016) . Notes 1 https://github.com/paezh a/covid 19-envir onmen tal-corre lates. 2 As the Globe and Mail, Canada's paper of record, put it in relation to the SARS outbreak in 2003:

"Globalization means that if someone in China sneezes, someone in Toronto may one day catch a cold" (Editorial, March 28, 2003, p. A18) 3 https://www.mitma.gob.es/minis terio/ covid-19/evolu cion-movil idad-big-data/movil idad-provi ncial. 4 https://www.geovo lunta rios.org/. 5 https://www.datos covid.es/pages/ provi denci aldata19. 6 https://www.datos covid.es/pages/ sobre-la-inici ativa. 7 As a check for robustness, we also tested the rook criterion (two provinces are adjacent if they share a boundary, but not if they only touch at a vertex), and included the three islands in the sample. In this case we made an allowance for adjacency between the two islands in the Autonomous Community of Canarias in the Pacific (Las Palmas and Santa Cruz de Tenerife), and assumed that Islas Baleares in the Mediterranean are adjacent to three provinces in Pais Catalans (i.e., Barcelona, Tarragona, and Castello) The results (which can be consulted in the source R markdown document) are robust to the specification of W. 8 We conducted sensitivity analysis letting all parameters vary over time, and while the results are qualitatively similar, the resulting models are less parsimonious. These results are available in the source R markdown document.

Rticles: Article Formats for R Markdown. Software Manual

Spsur: Spatial seemingly unrelated regression models. Software Manual

Spatial Econometrics: Methods and Models, Studies in Operational Regional Science

Spread of sars-cov-2 coronavirus likely to be constrained by climate

Ggthemes: Extra Themes, Scales and Geoms for 'ggplot2'. Software Manual

GridExtra: Miscellaneous Functions for "grid

Weathering the pandemic: How the Caribbean Basin can use viral and environmental patterns to predict, prepare and respond to COVID-19

Applied Spatial Data Analysis with R

Spatio-Temporal Evolution of Beijing 2003 Sars Epidemic

Effects of Air Temperature and Relative Humidity on Coronavirus Survival on Surfaces

The Effects of Temperature and Relative Humidity on the Viability of the SARS Coronavirus

Detecting Space-Time Patterns in Geocoded Disease Data. Cholera in London, 1854 Measles in the United States, 1962-95

Exponential Phase of Covid19 Expansion is not Driven by Climate at Global Scale

Vitamin D deficiency in children with recurrent respiratory infections, with or without immunoglobulin deficiency

Rapid surveillance of COVID-19 in the United States using a prospective space-time scan statistic: Detecting and evaluating emerging clusters

Association of vitamin D deficiency with acute lower respiratory tract infections in newborns

Spatial Econometrics: From Cross-Sectional Data to Spatial Panels

Vitamin D and respiratory tract infections in childhood

Economic Effects of Coronavirus Outbreak (COVID-19) on the World Economy. Available at SSRN 3557504

How to Interpret the Coefficients of Spatial Models: Spillovers, Direct and Indirect Effects

A balance act: minimizing economic loss while controlling novel coronavirus pneumonia

Econometric Analysis

Dates and Times Made Easy with Lubridate

Regional Specialisation and Concentration in the EU

Spatial Modeling Cannot Currently Differentiate SARS-CoV-2 Coronavirus and Human Distributions on the Basis of Climate in the United States

The Scaling of Contact Rates with Population Density for the Infectious Disease Models

Decline in Temperature and Humidity Increases the Occurrence of Influenza in Cold Climate

Vitamin D deficiency and lower respiratory tract infections in children: a systematic review and meta-analysis of observational studies

Global Trends in Emerging Infectious Diseases

Projecting the Transmission Dynamics of SARS-CoV-2 Through the Postpandemic Period

Low Ambient Humidity Impairs Barrier Function and Innate Resistance Against Influenza Infection

Is the Impact of Social Distancing on Coronavirus Growth Rates Effective Across Different Settings? A Non-Parametric and Local Regression Approach to Test and Compare the Growth Rate

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

A Spatiotemporal Analysis of Public Pharmaceutical Expenditure

Evaluating the Temporal and Spatial Heterogeneity of the European Convergence Process, 1980-1999

Introduction to Spatial Econometrics

Scientific and Ethical Basis for Social-Distancing Interventions Against COVID-19

Spatial Spillovers in Public Expenditure on a Municipal Level in Spain

Spatial Model Selection Strategies in a SUR Framework. The Case of Regional Productivity in EU

How Much of China and World GDP has the Coronavirus Reduced? Available at SSRN 3543760

Janssens Wim (2017) Vitamin D supplementation in respiratory diseases -evidence from RCT. Polish Archives of Internal Medicine

Cold Temperature and Low Humidity are Associated with Increased Occurrence of Respiratory Tract Infections

Strategic Planning for City Networks: The Emergence of a Basque Global City?

Understanding the Spatial Diffusion Process of Severe Acute Respiratory Syndrome in Beijing

Robust Standard Error Estimators for Panel Models: A Unifying Approach

ML Versus IV Estimates of Spatial SUR Models: Evidence from the Case of Airbnb in Madrid Urban Area

Distance Traveled in Three Canadian Cities: Spatial Analysis from the Perspective of Vulnerable Population Segments

Rapid Expert Consultation on SARS-CoV-2 Survival in Relation to Temperature and Humidity and Potential for Seasonality for the COVID-19 Pandemic

Perceived Risk and Modal Choice-Risk Compensation in Transportation System

Simple Features for R: Standardized Support for Spatial Vector Data

An Agent-Based Approach for Modeling Dynamics of Contagious Disease Spread

Risk Compensation and Bicycle Helmets

Projecting Social Contact Matrices in 152 Countries Using Contact Surveys and Demographic Data

US Regional Income Convergence: A Spatial Econometric Perspective

Condoms and Seat Belts: The Parallels and the Lessons

Aemet: Obtain Climatic and Meteorological Data from Spanish Meteorological Agency (Aemet). Software Manual

Trip Generation of Vulnerable Populations in Three Canadian Cities: A Spatial Ordered Probit Approach

Immobility Levels and Mobility Preferences of the Elderly in the United States Evidence from 2009 National Household Travel Survey

The Epidemiological Characteristics of an Outbreak of 2019 Novel Coronavirus Diseases (COVID-19)-China, 2020

Toward a Dialectical Understanding of Networked Disease in the Global City: Vulnerability, Connectivity, Topologies

Use of GIS and Agent-Based Modeling to Simulate the Spread of Influenza

Welcome to the Tidyverse

Isolation, quarantine, social distancing and community containment: pivotal role for old-style public health measures in the novel coronavirus (2019-nCoV) outbreak

Knitr: A Comprehensive Tool for Reproducible Research in R

Dynamic Documents with R and Knitr

No association of COVID-19 transmission with temperature or UV radiation in Chinese cities

Characterization of Regional Influenza Seasonality Patterns in China and Implications for Vaccination Strategies: Spatio-Temporal Modeling of Surveillance Data

An Efficient Method of Estimating Seemingly Unrelated Regressions and Tests for Aggregation Bias

Accelerated Contagion and Response: Understanding the Relationships Among Globalization, Time, and Disease

KableExtra: Construct Complex Table with 'Kable' and Pipe Syntax. Software Manual

März Winfried (2016) Vitamin D and airway infections: a European perspective