The effect of cold waves on mortality in urban and rural areas of Madrid

While many studies analyze the effect of extreme thermal events on health, little has been written about the effects of extreme cold on mortality. This scarcity of papers is particularly relevant when we search studies about extreme cold on the health of rural population. Therefore, we tried to analyze the effect of cold waves on urban areas and rural areas from Madrid and to test whether differentiated effects exist between both population classes. For this purpose, we analyzed data from the municipalities with over 10,000 inhabitants for the period from January 1, 2000 through December 31, 2013. Municipalities were classified as urban or rural (Eurostat), and they were grouped into similar climatological zones: Urban Metropolitan Centre (UMC), Rural Northern Mountains (RNM), Rural Centre (RC) and Southern Rural (SR). The dependent variable was the daily mortality rate due to natural causes per million inhabitants (CIE-X: A00-R99) that occurred between the months of November and March for the period. The independent variable was minimum daily temperature (ºC) (Tmin). Social and demographic contextual variables were used, including: population > age 64 (%), deprivation index and housing indicators. The analysis was carried out in three phases: (1) determination of the threshold temperature (Tthreshold) which defines the cold waves; (2) determination of the relative risk (RR) for cold waves using Poisson linear regression (GLM); and (3) using GLM of the binomial family, Odds Ratios (OR) were calculated to analyze the relationship between the frequency of the appearance of cold waves and the socioeconomic variables. The UMC zone experienced 585 extreme cold events related to attributable increases in the mortality rate. The average number of cold waves in the rural zones was 319. The primary risk factor was the percentage of population over age 64, and the primary protective factor was housing rehabilitation. As a whole, the period experienced more cold waves (1542) than heat waves (1130). The UMC was more vulnerable than the rural areas. Furthermore, the results support the development of prevention policies, especially considering the fact that cold wave events were more frequent than heat waves.


Introduction
Climate change is related to a global increase in temperatures. However, this will not necessarily translate into the disappearance of cold waves and their impact on mortality in the future [1,2]. On the contrary, there is ample scientific evidence that confirms greater risks to health related to cold than to heat at the global level [3][4][5][6][7][8]. In addition, the risks associated with cold waves remain constant in time, in contrast with what happens with heat waves [9].
On one hand, this phenomenon has climatological origins. Modifications in atmospheric and oceanic currents related to global warming result in the intensification of cold waves in middle latitudes [10,11]. On the other hand, in general, populations are adapting to higher and higher temperatures. This phenomenon is reflected in the increasing trend shown in Minimum Mortality Temperatures (MMT) observed in different populations [1,[12][13][14], which simultaneously could lead to poor adaptation to low temperatures. Thus, the most likely scenario is that cold waves will affect health at lower extreme temperatures in the future [1].
Given that different social factors play a key role in the evolution of this adaptation, it is important to consider these factors in the search for mitigation and adaptation strategies for climate extremes [2,15]. However, there are few studies that address this question in relation to cold waves. Those that do exist tend to focus on urban populations, possibly due to the absence of a universal definition of rurality and data limitation [16][17][18].
There are differences between urban and rural populations that justify differential behavior in the risks of temperature extremes [15]. One of these is the distribution of the population in age groups. In Europe, the rural population is older, and there are fewer people in working ages [19], which is precisely the group that is more vulnerable to low temperatures [20]. In terms of infrastructure and services, rural areas receive less investment than do urban areas [21]. There is also evidence of a lack of access to health services and digital health services. In this way, rural areas tend to have deficits in transportation infrastructure that make it difficult for the population to travel to health centers in other municipalities. These deficits include lack of public transportation and the high cost of gasoline involved in traveling long distances [22]. These add to the effects of exposure implied by the length of the workday and physical labor carried out in rural areas [18,23]. It is for these reasons that rural areas have been found to be more vulnerable to cold than urban areas [15,[24][25][26].
On the other hand, low economic status is related to difficulty covering energy costs, which prevents adequate insulation from exterior temperatures [27]. In addition, the housing of people with low economic status tends to present worse thermal properties [28], such that both factors impact the vulnerability of the population to low temperatures. These social inequalities are cause of concern in the different urban and rural zones of Madrid, along with the predicted negative health impact of the cold [27,29].
The aim of this study was to determine whether extremely low temperatures have a differential impact on population mortality in urban and rural zones in the province of Madrid. The secondary aim was to analyze whether the impact of extreme cold on health follows a pattern related to environmental and socioeconomic contextual variables such as housing conditions, the deprivation index (IP), the urban character (rurality) of the zone and the percentage of population over age 64.

Classification of the municipalities included in the study
The analysis included, as a criterion, all those province municipalities in Madrid with over 10,000 inhabitants. These were then classified into two levels, first, by the rurality of the municipality, and second, by its climatic patterns.
The rurality of municipalities was determined using the DEBURGA criterion as defined by Eurostat [30] (which is explained in greater detail in the bibliography [17,31]). The criterion establishes as "urban" all those municipalities that comply with two conditions. First, the population of the municipality must be concentrated primarily in the center of the municipality, and second, the densely populated central zone must contain at least 50,000 inhabitants. Using this definition of "urbanness", Eurostat provides a map of urban municipalities in Europe [32]. This map was used to classify the municipalities in this study. In contrast, those municipalities that are not defined as "urban" are "non-urban" or "rural", considered equivalent concepts in this study.
Municipalities were also grouped that had similar climatological behavior. "Isoclimatic" zones, were previously established by AEMET, which established those municipalities that share a common climate. This allowed us to assign exterior temperatures to the population analyzed that were representative of true exposure. In the case of Madrid, the province is divided into the zones of "Madrid Mountains", "Metropolitan and Henares" and "South, Vegas and West" [33]. Four groups were created by bringing these classifications together: "Rural Northern Mountains" (RNM), "Urban Metropolitan Center" (UMC), "Rural Center" (RC) and "Southern Rural" (SR). These groups made up the sample units of the analysis.
Once municipalities were classified, we developed an ecological, retrospective longitudinal time series analysis of the daily mortality rate (TM). The time series period concerned was 14 years, beginning on January 1, 2000 and ending on December 31, 2013. The statistical strategy was carried out in three phases, described in "Social context variables" section and beyond.

Variables used Main dependent variable
Mortality rate due to natural causes (TM): the variable of interest was daily TM per million inhabitants. This variable was determined for each municipality based on daily registries of deaths and population census data. Of all deaths registered, only those due to natural causes were used (CIE-10:A00-R99). The TM were calculated at the group level by aggregating the average of municipal data. Both the death count as well as the population census data were provided by the National Statistics Institute (INE).

Independent variable
Minimum daily temperature (T min ): minimum daily temperature was used (expressed in degrees Celsius), because it is the variable with the strongest association with winter mortality [3]. Values of T min were calculated as the daily average of the observations registered at all of the observatories found in the same meteorological prediction zone.
In this case, we used temperatures from 21 observatories. Six were situated in the "Madrid Mountains" isoclimate zone that includes only the RNM group. There are 11 observatories in the "Metropolitan and Henares" group, which includes UMC and RC. Finally, the isoclimate zone "South, Vegas and West", which includes the SR group, includes the four remaining observatories. The geographical distribution of the zones is shown in Fig. 1.
These data were provided by the State Meteorological Agency (AEMET).

Social context variables
The socioeconomic and demographic context variables used included: Deprivation index (IP): based on the transformation of the deprivation index carried out by the multi-center MEDEA project (IP2011) [34]. The index was originally calculated by census tracts (administrative sections below the municipal level). Based on these tracts, a value was calculated for the groups analyzed by aggregating the data for the census areas with the average. Later, the variable was transformed such that the zone with the lowest IP was equal to 1 (IP = 1). In this way, we established a reference framework in which precariousness could be directly interpreted based on the least socially vulnerable group. Hypothetically, an IP = 3 would indicate that for this group the deprivation is three times greater than in the reference group.   [35].
Housing rehabilitation: we had access to the number of housing rehabilitation licenses by municipality and by year included in the study. Based on these data, the value per group was calculated by aggregating the data for the whole period at the municipal level with the average. These data were provided by the registries published by the open access database of the regional statistics institute of the Community of Madrid [36].
Housing in good conditions: in the same way as what was described in the prior paragraph for "housing rehabilitation", housing in good condition was calculated for each group. These data were also provided by the Community of Madrid [36].
The following was included as a demographic variable: People aged > 64: the percentage of population over age 64 was available, which is the group most vulnerable to cold waves [20]. This variable was calculated for each group by aggregating with the average the annual population census of this age group at the municipal level. These data come from the continuous population census in the years 2000-2013 published by the INE [37].

Control variables
Variables for seasonality, trend and the autoregressive component of mortality were generated as control variables for the time series data.
Autoregressive component of mortality (torg1): a variable was generated of the first order of mortality to control for the autoregressive component [38]. Seasonality (Seas): oscillatory variables were calculated for the annual, biannual, cuatrimestral, trimestral, bimonthly and monthly periodicities, using sine and cosine functions. This group of variables provided for controlling the seasonal behavior of the mortality rate.
Trend: this variable acquires a value of 1 on January 1, 2000 and grows by a unit in each new observation. In this case, this variable allowed us to control for the trend of the time series.
Time: the month and year were included as two factor variables. Both variables offered an additional control for the seasonality of the time series of TM data.
Flu: flu was included as one of the variables that could explain winter mortality [24,39]. Epidemiological data on the flu in the province of Madrid were included. In this case, the variable was given as the weekly rate of the number of cases per 100,000 inhabitants. These data were provided by the Sentinel Surveillance System of Flu in Spain (SVGE)-National Network of Epidemiological Surveillance (RENAVE), which depends on the National Center for Epidemiology (CNE).
Finally, variables that were transformed based on those variables described here were necessary in some of the study phases, and they are mentioned hereafter. These new variables are described in their respective sections.

Phase 1: calculation of the threshold temperatures (T threshold ) of the definition of a cold wave
This study uses an epidemiological definition for a cold wave. This means that cold waves are characterized based on the relationship between mortality and temperature. Specifically, a threshold temperature (T threshold ) is established for which the TM is statistically greater than the season average. Once T threshold is known, it is considered that a cold wave exists for all of the days for which T min is lower than T threshold .
The T threshold was determined for each group based on the methodology described by Carmona et al. [3]. Therefore, ARIMA models were adjusted, using the complete time series (January-December) of the mortality rate as a dependent variable. Thereafter, our interest concerned the model errors (ε) given that they are free of periodicity, trend and the autoregressive component of the times series of TM.
Based on the literature, it was determined that the temperature's effect on mortality is not instantaneous. On the contrary, there would most likely be a first mortality peak attributable to a cold wave between two and six (including both) days after the date the population is exposed to low temperatures [40]. Therefore, to register this lagged effect, each observation of T min was assigned the value of the moving average of ε (εmovil) for the above-mentioned interval (2-6 observations later). After this, the observations from November to March were retained (winter period) and εmovil were grouped into intervals of two and two degrees. The global averages of εmovil, together with their 95% probability confidence intervals, are shown in Fig. 2, both at the global level as well as in temperature intervals. With the help of these graphs, it was impossible to determine for which values of T threshold TM was statistically greater than the seasonal average of the period.

Phase 2: determination of the relative risks (RR) and attributable risks (AR) associated with low temperatures
In phase 2, the relative risks (RR) and attributable risks (AR), that associated cold waves with mortality, were calculated. These risks were calculated based on the estimators from Generalized Linear Models (GLM) of the Poisson family (Link = log), which were calculated with the winter observations (November to March). In these models, the mortality rate was used as a dependent variable. In terms of the principal independent variable, we used the new variable represented by the cold wave (T cold ). This variable is calculated using the following expression: If T min ≥ T threshold then T cold = 0. If T min < T threshold then T cold = T threshold − T min . As previously mentioned, the effect of a cold wave on mortality is not immediate. To register the lagged effect of the impact of a cold wave in time, lag variables for T cold of orders 1 through 13 were included in the model [3]. Given that mortality describes markedly seasonal behavior across time, we incorporated variables for seasonality, trend and the autoregressive component of the 1st order for mortality [38]. Finally, flu incidence was included as a control variable. Given that the effect of flu on mortality is also long term, lag variables of orders 1 through 15 were incorporated [41]. 0 1 2 3 4 5 6 7 8 9 10 11 12 Temperatures (Cº)  0 1 2 3 4 5 6 7 8 9 10 11 12 Temperature (Cº) Residuals RNM Fig. 2 ARIMA model residual graph for urban metropolitan center. Axis X represents minimum daily temperature and axis Y represents the mean of the residual. The central band shows the overall mean residuals confidence interval with a 95% probability. This graph is plotted for each study group: Rural Northern Mountains (RNM), Rural Center (RC), Urban Metropolitan Center (UMC) and Southern Rural (SR) These models followed the general formula: where TM represents the mortality rate, a is the intercept, β k represents the coefficient of each of the k variables, lag(TM,1) is the autoregressive component of the first order of mortality in observation i, n1 is the value of the trend in observation i, seas represents seasonality j in observation i, time represents the time variable k in observation i, Flu represents the epidemiological data on flu in observation I, and lag(Flu, p) represents the lagged components of the variable Flu of order 0-15, and finally, T cold represents the degrees by which T threshold exceeds T min for observation I, and lag (T cold , g) represents the lagged components of T cold of order g, from 0 to 13, for observation i.
An independent model was calculated for each group analyzed. We began by including all of the variables in the model and sequentially discarding those that failed to reach statistical significance (p value < 0.005) by descending p value order. When found, overdispersion of the data was corrected using binomial negative regression. Based on the coefficients of the models described, RR and RA were calculated with usual methodologies [42,43]. Once the risks were determined, the same models were again calculated without the flu and its lag variables to determine the intensity of confounding.

Phase 3: quantification of the social and economic risk factors
In the final phase of the analysis, we analyzed to what degree the contextual and sociodemographic variables explained the results we found. In this case, our event of interest was the frequency of episodes of extreme cold associated with increases in mortality. Therefore, a new dichotomous variable was generated, called ColdW. ColdW was defined by the following expression: If T min < T threshold then ColdW = 1.
If T min > T threshold then ColdW = 0. We proceeded to analyze the association between the frequency of these events and the socioeconomic variables using GLM with logit link (linear regression models of the binomial family). Again, only the winter observations were retained (November through March). In these models, the dependent variable was ColdW, and the independent variables were the socioeconomic and demographic indicators described above.
First, the functioning of each of the variables with the ColdW variable was analyzed using univariate models. Later, new variables were sequentially introduced. In the multivariate models, variables were either retained log(TM) = a + β 1 lag(TM, 1) i + β 2 n1 i + β 3j seas ij + β 4k time ik + β 5p lag(Flu, p) i + β 6g lag T cold , g i or discarded based on biological meaning and statistical significance (p value < 0.05). These models were described using the following formula: where ColdW is the event of interest, a represents the intercept, β ij represents the index of context variable i in group j; ctxt j represents the context variables described above, in group j.
In this case, only the model including the observations of all of the groups was adjusted. Finally, based on the coefficients of these models, Odds Ratios (OR) were calculated associated with the variables considered, using the following expression: where β is the coefficient of context variable i.
The database was managed using the free software R version 5.3.1. To calculate the ARIMA models, we used the time series module of the IBM SPSS Statistics 25 software program. In the case of the Poisson and Binomial family models, we used STATA 14 software. The maps that show the results were generated using the free software QGIS3.10.2. Table 1 provides a descriptive picture of the groups analyzed. The table shows that each of the groups includes at least 15 municipalities, which also include a significant portion of the population: around 300,000 rural inhabitants and 5 million urban inhabitants. A high number of observatories was also included (21), dispersed widely across the whole province (Fig. 1). Figure 2 represents the graphs used to determine T threshold for each of the groups. Figure 3 shows a geographic representation of the values of T threshold , together with the percentiles of minimum daily temperatures for the corresponding winter months. In this way, higher temperatures imply effects on mortality and less extreme temperatures. Thus, the map shows that the most vulnerable area is the Metropolitan Center (P27). In relation to the rural groups, the most vulnerable is the Southern Rural (P18) which is closely followed by the Rural Northern Mountains (P16). Finally, the place with the greatest tolerance for extreme cold is the Rural Center, although the percentile is close to those of the other rural zones (P12). Table 2 shows the basic descriptive statistics of the variables used in the models of the Poisson regression.

Results
The greatest rates of winter mortality (deaths/10 6 inhabitants) occur in the urban area (TM = 20.0), followed by the Rural Northern Mountains (19.7). In terms of the climatology, the northern mountain zone tends to be the zone with the lowest T min . The climate of the rest of the province is less cold, especially in the central area. Finally, in quantifying the cold wave in terms of degrees (T cold ), it can be observed that in the urban center the cold wave average was 0.6 ºC. This same average oscillated between 0.2 and 0.4 in the rural periphery. It would seem, at the descriptive level, that temperatures were more intense in the urban zone. Table 3 shows that the risks related to extreme cold in the Rural Northern Mountains are greater than (RA = 5.2 (3.3, 7.1)) (and statistically significant) those in the other rural zones. In terms of the lagged effect of cold waves, it can be seen to act in terms of peaks in the short (lag 3) and long term, as shown in the urban zone and rural  northern mountain zone. There is a long-term effect (lag12) in the rural southern zone, and a medium term effect in the rural center (lag 5). In terms of the confounding effect of the epidemiological incidence of flu, after removing this variable from the models, we found variations in RR of less than 0.00% in all and in each one of the groups analyzed. Table 4 shows the variables included in the binomial models. In this case, the events of interest were those days in which temperatures were low enough to be associated with increases in TM (ColdW). The number of these events in the urban zone (585) was statistically higher than the average of the same in the rural areas (319) (p value < 0.05). Globally, in all of the province, 1,542 episodes of extreme cold were found, which is greater (and statistically significant) than the number of heat waves (1,130) previously found for the same groups in the same period [44]. Table 4 also shows the Southern Rural zone as the most socioeconomically vulnerable, with an IP almost 3.5 times larger than that observed in the Rural Center. It is also where there were fewer housing rehabilitations and fewer housing units in good conditions. Finally, the Metropolitan Center was the zone with the greatest population over age 64. The zone with the least population in this age group was the Rural Center. In contrast, the population over age 64 was similar in the northern and southern rural areas. Table 5 shows the results of the binomial regression models. The univariate models show that the deprivation index and the percentage of the population over age 64 are risk factors. In addition, the quality of housing and housing rehabilitation are protective factors against the effects of cold waves. In this case, the odds of episodes analyzed increased by 20% with the IP and decreased by 12% in the case of rehabilitated housing. The model with economic variables indicates that the most determinant variable was housing rehabilitation, which displaced IP   and the percentage of housing in good conditions. However, of all of the socioeconomic variables analyzed, the most determinant was population over age 64, which displaced all of the other variables in the model.

Discussion
The results shown above reveal that cold waves do not have a homogenous influence on mortality in the province of Madrid and that the population living in the metropolitan center area is more vulnerable. Table 1 and Fig. 1 show the representativeness of these results and are supported by a large volume of population included in the study as well as the wide geographic dispersion of the meteorological observatories that registered the temperature exposures. Figure 2 shows that rural zones experienced statistically significant increases in the TM below the threshold temperature of − 2 ºC. This threshold coincides with what has been shown previously for the province in work by Carmona et al. [3]. Thus, this threshold was used as a reference in the prevention plans to address cold waves in the Community of Madrid. In terms of prevention plans, alert level 1 against cold is activated for the whole province when predicted temperatures reach or go below − 2 ºC for three consecutive days [45]. However, this threshold does not coincide with findings from the urban area (0 ºC). Thus, these plans would not be activated between the degrees of 0 and − 2 ºC, despite the fact that they are within the range of temperatures that pose a risk to mortality. Therefore, prevention plans should be reevaluated to make geographical adjustments over time [33,46,47].
On the other hand, situating these threshold temperatures in their respective temperature percentiles allowed for a comparison of relative vulnerability between the analyzed areas. Figure 2 presents a descriptive representation that shows that the urban area is most vulnerable to cold waves (P27). This percentile was higher than those of similar studies that are restricted to the city of Madrid [1,3,46]. However, these studies monitored temperatures in the city of Madrid using as a reference only those data registered by the meteorological observatory of Madrid Retiro. In contrast, this study used the average temperature of all of the isoclimate area based on data from 11 observatories.
In translating these percentiles into events, the urban zone experienced a higher number of episodes of extreme cold associated with increases in mortality (585). These same events occurred with lower frequency in the rural periphery, with an average of 319. In comparing the differences between both groups, the Chi-square test indicated that these are statistically significant (p value < 0.05). Therefore, these results show that the urban zone of Madrid is more vulnerable to cold than the nonurban areas they are compared to.
Although there are few publications on this topic, these results contradict part of what has been suggested by the literature. In general, it has been assumed that rural areas are more vulnerable to meteorological disasters, including temperature extremes [15]. A study by Jegasothy [25] in New South Wales (Australia) found a more pronounced effect of cold waves on mortality due to all causes in remote areas compared to large cities. On the other hand, Urban et al. [48] found greater mortality attributable to cold waves for deaths due to heart attacks (a severe pathology whose prognosis depends in large part on having access to a healthcare center nearby) in South Bohemia (a rural zone) than in Prague (an urban zone). However, it should be noted that their climate and social context is different from that of the population analyzed. This tends to be interpreted in terms of the demographic weight of those over age 64 in rural areas and problems of access to healthcare centers [15,24]. Although this might be the general case for Europe [19], it is not the case for the province of Madrid (Table 4). This could be due to the effect of "bedroom communities", or traditionally rural or non-urban municipalities that provide housing and are close to cities that provide work and services [49].
In the case of the creation of the metropolitan area of Madrid, the decade of the 1960s was marked by a process in which immigrant groups settled in rural areas of the southern and eastern zones in their search for work [50]. This trend has been perpetuated in time given that housing prices tend to decline with greater distance to the city center of Madrid [51]. Therefore, this phenomenon could explain why there is a younger population in the non-urban periphery than in the metropolitan area (Table 4). On the other hand, a comparison of percentiles shows that, among rural zones, the Southern Rural is the most vulnerable (P18). A greater precariousness (IP = 3.42), a greater percentage of housing in poor conditions (95.55%), and a lower number of housing rehabilitations are all present in this zone (9.79 licenses/million inhabitants). On the opposite extreme, the region least vulnerable to cold waves is the Rural Center, with a lower percentage of population over age 64 (7.88%). Furthermore, it is the zone with the lowest levels of precariousness (IP = 1) and where more licenses for housing rehabilitation were provided (14.93). Therefore, it would seem that these variables explain the results we obtained.
In terms of the social and demographic variables, Table 4 confirms the observations at the descriptive level. The univariate models show that the percentage of the population over age 64 and the IP are statistically significant risk factors. These results coincide with what has been observed in the neighborhoods of Madrid and are in line with scientific consensus [2,29].
On the other hand, the univariate models show that the percentage of housing in good condition and housing rehabilitation licenses are protective factors. In addition, in the multivariate model of economic variables, housing rehabilitation replaced the other variables in the regression model. Thus, housing rehabilitation seems to be the most determinant variable in explaining the results found. This, therefore, supports cold mitigation policies supported by the identification of "vulnerable buildings" to which rehabilitation, insulation and energy efficiency programs can be applied. In addition, inasmuch as these factors are associated with inefficient energy consumption, this would contribute to reducing energy poverty and the emission of greenhouse gases [24,26,27].
In the global model (Table 4), the percentage of population over age 64 displaced the rest of the socioeconomic context variables considered. This result agrees with the etiology of cold, for which the most vulnerable population is those over age 64 [20,24,29,40,52]. Table 2 shows RR and that greater percentiles are associated with more reduced RR. In this way, the urban zone experienced a greater number of cold waves, but the risks associated with each degree tended to be lesser than in the rural periphery-although the differences are statistically significant only in the case of the Rural Northern Mountains. In terms of the lagged effect of cold waves, the short-term effect observed in the Rural Northern Mountains and the metropolitan area (lag3), as well as in the Rural Center (lag5), tends to be associated with worsening of cardiovascular illnesses. On the other hand, a long-term effect was detected in the Northern Rural Mountains (lag13), the urban metropolitan area (lag9) and the Southern Rural zones (lag12). This longer term effect over time did not manifest in the Rural Center, where the demography is of a younger population. The longer term effect mentioned also tends to be related to negative prognosis for respiratory illnesses [3,24,40,53]. It is important to note that the confounding effect of flu on cold waves was negligible, with modifications lower than 0.00% in terms of the RR associated with T cold among all of the groups. Therefore, it would seem that controlling for this variable is not necessary.
Finally, a similar study determined that there were 1,130 heat waves in the same areas analyzed in this study [44]. This number is lower than the 1,542 cold waves we detected (Table 4). In comparing these data, the Chisquare test indicated that these differences are statistically significant (p value < 0.005). This differs from what was found in another study [3], where it was determined that there were more heat waves (211) than cold waves (30) in the capital of the province of Madrid during the 2000-2009 period. This discrepancy shows the limitations of extrapolating the results of a capital city to an entire province, which is composed of heterogeneous subpopulations. Thus, it seems important to carry out further studies to determine the effects of temperature extremes at the local level [54].
On the other hand, the greater vulnerability to cold compared to heat, shown by the results described, agrees with the finding that the risks associated with cold waves remain constant across time, or are even greater than those related to heat [1, 4-6, 8, 9]. This contradicts the lower level of development of action plans for cold in comparison with what has been developed for heat. Although many countries in the European Union do not have these plans in place at the national level (for example, in Spain), evaluations of Heat Health Action Plans (HHAP) have provided conclusions that can be extrapolated to the development of action plans to address cold [55]. First, it is important to note that these plans have been effective [56,57], although they could be improved by making them more specific and directing them more to vulnerable populations according to income, gender or occupational health risks, among other characteristics [29,[58][59][60][61][62]. In addition to the implementation of these plans at the national level, it is also important to adapt them to the local level to make them more precise and efficient [63,64]. Finally, to achieve these objectives, further scientific study is needed. Such study should cover the still unknown aspects related to the impact of temperature extremes on health, the adaptation process and mitigation strategies [54].
As Cold health Action Plans does not exist at country level, their effectiveness cannot be established [3,20]. However, there is a "Plan for Surveillance and Control of the Cold effects on Health" specifically in Madrid province [65]. Thus, this protocol provides actions like to inform and vaccinate against flu the vulnerable population, or social services attention to homeless people. Nevertheless, its effectiveness has not yet been evaluated.
In relation with these plans, the first conclusion from the results obtained by us is that they should be enhanced and implemented in both Madrid and the rest of the country. One area for improvement is the promoting a "Cold culture" which have gotten good results with the heat waves [1,66,67]. Furthermore, they must imply vulnerable people's dwellings rehabilitations plans and ensure the right to access to heater systems in coldest months [29]. An taking account of the recent heat health action plan assessment [68], which is also suitable for cold waves, these plans must be implemented at the local scale and periodically updated.
It is important to note that there are limitations to this study, one of which relates to its ecological nature [69]. These types of studies present a disadvantage in that it is not possible to make inferences about individual patients. That is to say, the results obtained are only valid at the population level. Another problem related to these studies is that the exposure variables are measured in a way that is delocalized from the exposed population. Finally, the interpretation of the effect of the social indicators used was limited by the effect of the ecological fallacy. Despite these considerations, these limitations are common among all ecological studies of this type.
On the other hand, this study did not include confounding environmental variables, because of the absence of suitable data. This could be relevant in urban results because of pollution and noise, which are typical urban environmental issues due to traffic [70,71]. This limitation, which is common in these studies [72][73][74][75], can be partially minimized by the use of control variables such as seasonality and trend in the models [76]. There are other similar studies, in which these ambient variables had a reduced effect on the lagged combined coefficient associated with extreme temperatures [73][74][75]77]. Therefore, other environmental variables of interest, such as air pressure and relative humidity, were not included in the analysis owing to their negligible relevance in the temperature-mortality relationship [75,77]. Insofar as particulate air pollution is concerned, the lack of quality data at a subregional level rendered the use of such data unadvisable, owing to the risk of introducing instability into the complete series, in the sense that introducing air pollutants into the GLM models could slightly modify the impact of extreme temperatures on mortality [74].
In relation to the meteorological stations, its number is unbalanced between the isoclimatic areas. This can led biased measures of temperature among the analyzed groups so that the higher the stations number, the higher the accuracy of the daily minimum temperatures values. Moreover, although this study included information of more meteorological stations than other similar ones, there are few observatories in relation to the study area. Therefore, this study suffers from Berkson-type measurement error, which decreases statistical power. Nevertheless, it is another common issue in ecological studies. Unfortunately, there are no better measurements since we have included all existing stations in the area.
In addition, there is some limitation in relation to the Privation Index which have been used in the analysis. On the one hand, this index was calculated with data corresponding to the 2011 census, while the study period runs from 2000 to 2014. Therefore, one must note that the economic crisis, which had an important impact in the indicators represented by this index, broke out between 2007 and 2008 in Spain. But the relative ranks of the study areas probably did not change in a determinant way even so all them got worse social indicators. On the other hand, the index was validated only for largest cities, but not for small cities neither rural areas. Nevertheless, a sensitive analysis was made for small and rural areas [34].
In relation with the rehabilitation, licenses may be a selection bias. The number of new-dwelling licenses increased much faster than rehabilitation licenses between 2000 and 2007 [78]. Thus, the economic crisis in 2008 led the fall out of both type of licenses [78,79]. However, the new housing market went into crisis while the number of rehabilitation licenses decreased slightly and kept relatively stable. Therefore, this variable may captures dwellers badly affected by the crisis who kept enough wealth for improving their homes. This type of dwellers probably had houses with better thermal properties than poorer dwellers even before doing the reforms. In other words, if government's campaigns focuses on social vulnerable population the potential protective effect could be higher.

Conclusions
The principal conclusion of this study is that the metropolitan area of Madrid is more vulnerable to cold waves than the areas in the rural periphery. The two factors detected that explain the pattern of cold waves were housing rehabilitations and the demographic population over age 64. The results of this study suggest the need for greater development and implementation of action plans to address cold, both at the national and local levels.