Evidence for Urban-Rural Disparity in Temperature-Mortality Relationships in Zhejiang Province, China.

BACKGROUND: Temperature-related mortality risks have mostly been studied in urban areas, with limited evidence for urban-rural differences in the temperature impacts on health outcomes.
OBJECTIVES: We investigated whether temperature-mortality relationships vary between urban and rural counties in China.
METHODS: We collected daily data on 1 km gridded temperature and mortality in 89 counties of Zhejiang Province, China, for 2009 and 2015. We first performed a two-stage analysis to estimate the temperature effects on mortality in urban and rural counties. Second, we performed meta-regression to investigate the modifying effect of the urbanization level. Stratified analyses were performed by all-cause, nonaccidental (stratified by age and sex), cardiopulmonary, cardiovascular, and respiratory mortality. We also calculated the fraction of mortality and number of deaths attributable to nonoptimum temperatures associated with both cold and heat components. The potential sources of the urban-rural differences were explored using meta-regression with county-level characteristics.
RESULTS: Increased mortality risks were associated with low and high temperatures in both rural and urban areas, but rural counties had higher relative risks (RRs), attributable fractions of mortality, and attributable death counts than urban counties. The urban-rural disparity was apparent for cold (first percentile relative to minimum mortality temperature), with an RR of 1.47 [95% confidence interval (CI): 1.32, 1.62] associated with all-cause mortality for urban counties, and 1.98 (95% CI: 1.87, 2.10) for rural counties. Among the potential sources of the urban-rural disparity are age structure, education, GDP, health care services, air conditioners, and occupation types.
CONCLUSIONS: Rural residents are more sensitive to both cold and hot temperatures than urban residents in Zhejiang Province, China, particularly the elderly. The findings suggest past studies using exposure-response functions derived from urban areas may underestimate the mortality burden for the population as a whole. The public health agencies aimed at controlling temperature-related mortality should develop area-specific strategies, such as to reduce the urban-rural gaps in access to health care and awareness of risk prevention. Future projections on climate health impacts should consider the urban-rural disparity in mortality risks. https://doi.org/10.1289/EHP3556.

Introduction

Nonoptimum temperatures (either heat or cold) have been widely documented to be associated with increased risks of cause-specific mortality, such as cardiovascular and respiratory mortality, mostly in developed countries (Anderson and Bell 2009; Gasparrini et al. 2015; Guo et al. 2014). Most of these studies on temperature–mortality relationships have been conducted in urban areas (Analitis et al. 2008; Basu 2009; Madrigano et al. 2015b; Medina-Ramón and Schwartz 2007). In contrast, few studies have been performed in rural areas (Hashizume et al. 2009; Todd and Valleron 2015) because of the lack of sufficient meteorological and health data. In addition, most global, national, and regional temperature-related mortality projections usually use the same exposure–response associations for the whole population (Ballester et al. 2011; Huang et al. 2011; Takahashi et al. 2007), which ignores possible urban–rural differences and may result in an incorrect estimation of the temperature-related health burden.

Due to the urban heat island (UHI) effect, it is usually assumed that urban residents are at a higher risk of extreme heat than rural dwellers (Heaviside et al. 2017; Tan et al. 2010; Tomlinson et al. 2011). However, rural residents are also sensitive to nonoptimum temperatures and exhibit different patterns of vulnerability from those of urban populations (Kovach et al. 2015; Sheridan and Dolney 2003). For example, their living conditions and their outdoor occupations mean that they are more frequently exposed to extreme temperatures. Moreover, they might have limited risk awareness and limited access to health services, particularly in developing countries (Bai et al. 2016; Li et al. 2017; Williams et al. 2013).

To date, very few studies have compared temperature-related mortality risks between urban and rural areas (Henderson et al. 2013; Li et al. 2016; Urban et al. 2014), and the results are mixed. For example, Gabriel and colleagues reported a greater increase in mortality during heat waves in the city of Berlin, Germany, compared with neighboring nonurban areas (Gabriel and Endlicher 2011). In contrast, in Jiangsu and Hubei provinces in China, findings indicated a higher mortality increase associated with high temperatures in less urbanized counties compared with urban counties (Chen et al. 2016b; Zhang et al. 2017b). Of particular note is that most of these studies focused exclusively on heat-related mortality without assessment of health risks associated with cold temperatures. However, cold effects were reported to account for most of the temperature-related attributable mortality burden in a multicountry study (Gasparrini et al. 2015). Therefore, the urban–rural disparity in mortality risks associated with both cold and heat deserves to be further investigated with equal attention. A better understanding of these risks is important for effective decision support to design spatially targeted interventions and mitigation policies. This is especially relevant for developing countries, which are more sensitive to extreme temperature events and often lag behind developed countries in health risk management capacities (Laboy-Nieves et al. 2010; Mendelsohn et al. 2006; Tol et al. 2004).

Epidemiological studies have predominantly examined temperature–mortality associations in a city using temperatures from one site or the average from a network of sites (Guo et al. 2013), which induces exposure measurement error and biases the estimates. Satellite-measured land surface temperature (LST) has been used to identify the temperature variations at a high spatial resolution (Madrigano et al. 2015a; Wan 2008). However, LST cannot serve as a proper proxy for the daily mean temperature because only two images are generally available within a day, and cloud-contaminated values often lead to incomplete data (Wan 2008). Moreover, despite sometimes high correlations, LST cannot be used as a direct substitute for ambient temperature due to the complex relationship between them (Vancutsem et al. 2010). Some studies have spatially interpolated temperature data from multiple sites and then used the average temperature of the study unit (Madrigano et al. 2015a; Urban et al. 2014), but this was often limited by the sparse distribution of weather stations that cannot accurately measure temperature variations. Here we used data from a highly dense network of weather stations in Zhejiang Province, China, to investigate how temperature–mortality associations vary between urban and rural counties and whether the temperature-related mortality risks are higher in urban counties compared with their rural counterparts. The potential sources of the urban–rural differences were also explored.

Materials and Methods

Study Area and Population

This study was conducted in Zhejiang Province in China (Figure 1), which has a total area of and had a population of 54.4 million in 2010. Based on the Zhejiang Planning Bureau’s categorization of districts and counties, there are 89 counties in Zhejiang Province, which include 29 “main city zones” (termed as urban counties below) and 60 rural counties (Table 1). Urban counties are home to 22.2 million inhabitants, which account for 40.7% of the total population of Zhejiang (2010). The climate in Zhejiang is humid subtropical, with four distinct seasons, characterized by long, very hot, humid summers, and chilly, cloudy, and dry winters with occasional snow. The mean annual temperature of Zhejiang is 17.0°C, with the monthly mean temperature ranging from 4.8°C in January to 29.1°C in July.

Figure 1.

(A) Locations of study area and automatic weather stations (some outside the province not shown), and (B) average of daily mean temperature across Zhejiang Province, 2009–2015.

Table 1

Summary data of number of counties, number of residents, cause-specific death counts, daily weather conditions, and daily air pollution of urban and rural counties in Zhejiang Province, 2009–2015.

	Urban	Rural
Total number of counties	29	60
Total number of residents (10^6)	22.2	32.2
Mean temperature (°C)
1st	0.4	−1.2
25th	9.5	8.7
50th	18.8	17.5
75th	25.1	23.9
99th	33.4	31.5
Mean (SD)	17.4 (9.1)	16.3 (9.0)
Range	−2.0, 35.3	−4.5, 34.3
Relative humidity [%, mean (SD)]	74.4 (12.8)	75.7 (12.5)
Air pollution [μg/m^3, mean (SD)]
PM_10	81.2 (44.5)	69.3 (40.2)
Ozone	87.0 (35.9)	86.6 (35.8)
Total death count
All cause	671,177	1,398,340
Nonaccidental	614,323	1,254,434
Age 0–64	129,003	327,069
Age 65+	485,320	1,194,761
Males	345,759	869,567
Females	268,564	652,262
Cardiopulmonary	311,720	671,694
Cardiovascular	101,856	226,467
Respiratory	209,864	445,227
Daily death count [mean (SD)]
All cause	9.4 (6.2)	9.0 (5.8)
Nonaccidental	8.6 (5.7)	8.0 (5.3)
Age 0–64	1.8 (1.7)	1.8 (1.7)
Age 65+	6.8 (4.7)	6.5 (4.5)
Males	4.8 (3.5)	4.7 (3.4)
Females	3.8 (2.9)	3.5 (2.8)
Cardiopulmonary	4.4 (3.4)	4.3 (3.4)
Cardiovascular	2.9 (2.5)	2.9 (2.4)
Respiratory	1.4 (1.6)	1.5 (1.6)

Note: 1st, 25th, 50th, 75th, and 99th refer to the percentiles of mean temperature. SD, standard deviation.

Data Collection

Meteorological and air pollution data.

Temperature and relative humidity data, obtained from a highly dense network of 4,007 automatic weather stations (AWSs) across Zhejiang Province (2,430 AWSs) and its neighboring provinces (1,577 AWSs) (Figure 1A), were acquired from the Zhejiang Meteorological Bureau and underwent an extensive automated quality control to eliminate random errors. Systematic errors of meteorological data were removed using the quality control system integrating a) laboratory calibration, b) periodic maintenance services, c) automated routines, and d) manual inspection. The new European Centre for Medium-Range Weather Forecasts Interim Reanalysis (ERA-Interim), which is the latest global atmospheric reanalysis dataset, was used to provide gridded estimates of three-dimensional temperature and relative humidity at a resolution of (Dee et al. 2011). Hourly temperature and relative humidity data from January 2009 to December 2015 were interpolated to a 1-km resolution using a method described in a previous study (Chen et al. 2016a), based on meteorological observations, ERA-Interim, and a digital elevation model (DEM). This method divides the model prediction error into that attributable to the modeling of weather systems (ERA-Interim) and that which describes the topography and then uses the observation data to revise these two components separately. Validation was carried out by using 80% of the observations as the training data and retaining 20% of the observations as the validation data. The cross-validation for the predicted daily mean temperature shows that the interpolation method produced accurate predictions with little bias (; ).

The human settlement index (2010) at a 1-km resolution, which can be an accurate proxy for population density, was obtained from a previous study based on the Defense Meteorological Satellite Program/Operational Linescan System nighttime light imagery, enhanced vegetation index, and DEM data (Yang et al. 2013). The population-weighted averages of the daily 24-h mean temperature and the daily 24-h mean relative humidity in each county were then estimated in attempt to accurately measure human exposure (Qi et al. 2012). The daily 24-h average (particulate matter with an aerodynamic ) and the daily maximum 8-h average ozone from January 2013 to December 2015 were calculated by averaging the data within each county from the 174 air pollution monitoring stations across Zhejiang Province, obtained from the Zhejiang Meteorological Bureau.

Mortality data.

County-specific daily mortality data during the study period (2009–2015) were obtained from the Zhejiang Center for Disease Control and Prevention based on the National Death Registration Reporting Information System. According to the Tenth Revision of the International Classification of Diseases (ICD-10; WHO 2016), the mortality data for all ages were classified into the following five categories: all-cause mortality, nonaccidental mortality (codes A00–R99), cardiopulmonary mortality (codes I00–I99 and J00–J99), cardiovascular mortality (codes I00–I99), and respiratory mortality (codes J00–J99). Additionally, nonaccidental mortality was stratified by age (0–64 y, ) and sex (male, female).

Census data.

County-specific socioeconomic and demographic characteristics were obtained from 2010 Population Census of China (Population Census Office and National Bureau of Statistics of China 2012) and 2013 Hangzhou Statistical Yearbook (Hangzhou Statistical Bureau 2013) and summarized by urban and rural categories in Table 2.

Table 2

Summary of socioeconomic and demographic characteristics in 29 urban and 60 rural counties of Zhejiang Province.

Mean (SD)	Definition	Sources	Unit	Urban counties	Rural counties
Male	Percent population of males	Population Census Office and National Bureau of Statistics of China (2012)	%	51.5 (46.5)	51.2 (50.6)
Children	Percent population <15y of age		%	12.0 (2.6)	12.7 (3.3)
Elderly	Percent population >65y of age		%	9.3 (2.3)	10.6 (1.9)
Education	Average years of schooling		year	9.7 (1.2)	8.1 (0.5)
Primary industry employment	Percent population employed in primary industry (agriculture, forestry and fisheries)		%	7.3 (7.1)	26.8 (16.0)
Hospital beds	Number of hospital beds per 10^4 people		/10^4 people	88.0 (74.8)	32.4 (14.7)
Air conditioner	Number of air conditioners per household	Hangzhou Statistical Bureau (2013)	/household	1.7 (0.2)	1.2 (0.2)
GDP	Gross domestic product per capita		10^3 RMB	87.1 (43.2)	58.5 (31.5)

Statistical Analysis

Two-stage time-series analyses.

We used a two-stage analytic approach to perform the time-series analysis. In the first stage, we estimated county-specific increases in mortality risks using a standard quasi-Poisson generalized linear model (GLM) allowing for overdispersion (Bhaskaran et al. 2013). The relationship between temperature and mortality was modeled using a distributed lag nonlinear model (DLNM) (Gasparrini 2014), which provided a modeling framework to flexibly describe both nonlinear and delayed effects. A cross-basis function was defined using a quadratic B-spline with two internal knots of temperature and a natural cubic spline for the space of lag days with 4 degrees of freedom (df). We placed the temperature knots at equally spaced quantiles of the daily urban or rural average mean temperatures, as urban and rural counties may have different exposure and vulnerability characteristics. We placed the lag knots along the logarithmic scale to account for a higher variability at lower lags up to a maximum of 21 lags (Guo et al. 2014; Lokys et al. 2018). The choice of 21 lag days was due to the fact that cold effects often appear some days after exposure and persist for several days, whereas hot effects are immediate and possibly have a harvesting effect (Analitis et al. 2008; Guo et al. 2012). A previous study suggested there was no one temperature measure (e.g., daily mean, maximum or minimum temperature) that was superior to the others (Barnett et al. 2010); thus, we computed daily 24-h mean temperature to estimate population exposure. The GLM regression also included the following covariates: a) a natural cubic spline for time with 7 df per year to control for seasonal and long-term trends, and b) categorical variables for day of the week and public holidays.

In the second stage, we used multivariate meta-analysis to pool the urban estimates and the rural estimates of overall temperature–mortality associations by combining three sets of county-specific parameters obtained from the reduction of the first-stage model (Gasparrini and Armstrong 2013). The meta-analyses were fitted using a random effects model by maximum likelihood.

The minimum mortality temperature (MMT) and minimum mortality temperature percentile (MMP), corresponding to the minimum mortality during the study period (2009–2015), were derived from overall pooled temperature–mortality associations for cause-, age- and sex-specific mortality (Gasparrini and Armstrong 2013; Gasparrini and Leone 2014). To calculate the cumulative relative risks (RRs) at the 1st and the 99th percentiles of mean temperature with 95% confidence intervals (CIs), we centered DLNMs at the MMTs as references. The 1st and 99th percentiles of mean temperature for urban counties and rural counties were calculated from the combined temperature data of 29 urban counties and the 60 rural counties, respectively. RRs of rural counties vs. urban counties were estimated using meta-regression of county-specific RRs by a binary variable (i.e., , and ).

Calculation of fractions and number of deaths attributable to nonoptimum temperatures.

The attributable counts of deaths and their corresponding fractions caused by nonoptimum temperatures were calculated by the sum of the contributions from the whole study period using a forward perspective (Gasparrini and Leone 2014). For each county, the number of deaths attributable to heat or cold in city i on day t () was estimated using the method described by Gasparrini and Leone (2014), with empirical CIs (eCIs) estimated using Monte Carlo simulation (5,000 random samples): where is the cumulative risk of cause-specific mortality in the following 0–21 d associated with daily mean temperature in city i on day t, in comparison with the MMT. is the moving average of daily cause-specific death counts in city i in the following 0–21 d since day t. The total number of attributable numbers of deaths in city i () was calculated by summing corresponding to days with temperatures lower or higher than the MMT during 2009 to 2015. The attributable fraction () in city i was calculated by dividing () by the total number of cause-specific deaths during the corresponding study period.

Potential sources of the heterogeneity between counties.

In addition, a meta-regression with the percentage of the urban population (termed as the urbanization level, obtained from the 2010 Population Census of China) as an independent variable was developed to check whether the urban–rural differences in the temperature–mortality associations could be explained by the urbanization level. Heterogeneity between counties was statistically assessed using the index and Cochran Q test (Gasparrini and Armstrong 2013). The statistics were used to quantify the extent of heterogeneity between counties by measuring the percentage of variability due to the true differences across counties rather than chance. The Wald test was also performed to determine if the urbanization level describes a significant modification to the original model. These tests were employed in the RRs for heat (RRs at the 99th percentile of temperature vs. MMT) and RRs for cold (the 1st percentile vs. MMT), and the cumulative temperature–mortality associations.

In order to explore the potential sources of the heterogeneity of the RRs for heat (99th percentile vs. MMT) and RRs for cold (1st percentile vs. MMT) across 89 counties, a mixed-effects model meta-regression analysis was also performed. Single predictors in meta-regression models are nine socioeconomic, demographic, and meteorological variables including: a) male, b) children, c) elderly, d) education, e) GDP, f) primary industry employment, g) air conditioner, h) hospital beds, and i) annual mean temperature. The percentage changes of RRs for heat (99th percentile vs. MMT) and RRs for cold (1st percentile vs. MMT) on nonaccidental mortality per interquartile range (IQR) increase of the abovementioned nine county-level characteristics were calculated. A meta-regression with the abovementioned nine multiple variables was also performed. Variance inflation factors (VIFs) were calculated to explore the potential collinearity in the multiple meta-regression. Moreover, Spearman’s correlation between urbanization level and other county-level characteristics was checked.

Sensitivity analyses.

Sensitivity analyses were performed to check the robustness of our findings by changing lag days from 7 to 28 d, changing the df for time (3–10 df per year) and lag days (3–6 df), and controlling for relative humidity at 0–2 lag days. We also used the 25th and 75th percentiles of the temperatures as reference values for comparing the results with a previous study (Chen et al. 2016b). In addition, we calculated the RRs of the 95th percentiles vs. MMT and the 5th percentiles vs. MMT. Moreover, because air pollution data were not available before 2013, we controlled for and ozone between 2013 and 2015 by including and ozone at 0–2 lag days (Hu et al. 2018).

All the analyses were performed with the R software (version 3.3.2; R Project). The dlnm package was used to create the DLNM (Gasparrini 2011), and the mvmeta and metafor packages were used to conduct the meta-analysis (Gasparrini and Armstrong 2013; Viechtbauer 2010).

Results

Descriptive Statistics

Daily meteorological variables, daily air pollution, and daily cause-specific death counts for urban and rural counties in Zhejiang are summarized in Table 1. The average daily mean temperature in urban counties is 1.1°C higher than in rural counties, indicating the presence of the UHI effect (Figure 1B). This analysis included total counts of 2.1 million all-cause deaths, 1.9 million nonaccidental deaths, 1.0 million cardiopulmonary deaths, 0.7 million cardiovascular deaths, and 0.3 million respiratory deaths (Table 1). The average all-cause, nonaccidental, cardiovascular, respiratory, and cardiopulmonary mortality rates were 1.5-fold, 1.5-fold, 1.5-fold, 1.6-fold, and 1.5-fold higher in rural counties than urban counties, respectively.

Urban–Rural Disparity in Temperature–Mortality Relationships

The associations between temperature and cause-specific mortality have different characteristics. However, the pooled urban and rural temperature–mortality curves show that the RRs at both the same absolute extreme temperatures and the same extreme percentile (1st and 99th) of the temperature are generally higher in rural counties compared with urban counties (Figure 2). Table 3 further displays the RRs for heat (99th percentile vs. MMT) and RRs for cold (1st percentile vs. MMT) on cause-specific mortality. A comparison of the relative urban–rural risks shows that all-cause mortality risks are 1.30 (95% CI: 1.22, 1.39) times higher for cold (99th percentile vs. MMT) and 1.01 (95% CI: 0.82, 1.24) times higher for heat (1st vs. MMT) among rural residents compared with urban residents.

Figure 2.

Pooled temperature–mortality associations along lag 0–21 d for cause-specific mortality for urban and rural counties in Zhejiang Province, 2009–2015, with 95% confidence intervals (CIs). Note: The vertical lines represent the minimum mortality temperature (MMT, solid) and the 1st and 99th percentiles of the temperature distribution (dashed) for 29 urban counties and 60 rural counties in Zhejiang Province, 2009–2015. The histograms represent the distributions of the daily averages of mean temperatures of urban and rural counties in Zhejiang Province, 2009–2015. The shading lines represent the 95% CI areas for risk estimates. Distributed lag nonlinear models (DLNMs) were used to model the exposure–lag–response associations between temperature and mortality. A cross-basis function was defined using a quadratic B-spline with two internal knots of temperature and a natural cubic spline for the space of 21 lag days with 4 degrees of freedom. RR, relative risk.

Table 3

Cumulative relative risks for cold (1st vs. MMT) and for heat (99th vs. MMT), relative risks of rural counties vs. urban counties, p-value for urban–rural difference along lag 0–21 days for cause-specific mortality for urban and rural counties in Zhejiang Province, 2009–2015, with 95% confidence intervals.

Cause-specific mortality and subgroups	1st vs. MMT (cold)				99th vs. MMT (heat)
	Relative risks		Relative risks of rural vs. urban	p-Value for urban–rural difference	Relative risks		Relative risks of rural vs. urban	p-Value for urban–rural difference
	Urban	Rural			Urban	Rural
All cause	1.47 (1.32, 1.62)	1.98 (1.87, 2.10)	1.30 (1.22, 1.39)	0.0004	1.15 (1.07, 1.24)	1.18 (1.14, 1.23)	1.01 (0.82, 1.24)	0.2
Nonaccidental	1.50 (1.35, 1.66)	2.04 (1.92, 2.16)	1.31 (1.23, 1.39)	0.003	1.12 (1.04, 1.20)	1.16 (1.12, 1.21)	1.03 (0.87, 1.20)	0.06
Age 0–64	1.24 (1.05, 1.46)	1.56 (1.39, 1.76)	1.20 (0.81, 1.79)	0.3	1.04 (0.91, 1.20)	1.02 (0.96, 1.08)	0.98 (0.63, 1.52)	0.6
Age 65+	1.59 (1.36, 1.87)	2.18 (2.05, 2.32)	1.34 (1.14, 1.57)	0.005	1.11 (0.98, 1.26)	1.20 (1.15, 1.25)	1.05 (0.73, 1.51)	0.5
Males	1.53 (1.27, 1.84)	2.00 (1.86, 2.14)	1.26 (1.06, 1.49)	0.04	1.09 (0.97, 1.22)	1.11 (1.06, 1.16)	1.01 (0.64, 1.60)	0.6
Females	1.48 (1.28, 1.71)	2.09 (1.91, 2.28)	1.37 (1.14, 1.64)	0.003	1.12 (0.96, 1.31)	1.24 (1.17, 1.30)	1.09 (0.90, 1.33)	0.1
Cardiopulmonary	1.64 (1.45, 1.85)	2.35 (2.17, 2.54)	1.39 (1.27, 1.52)	0.009	1.16 (1.04, 1.29)	1.26 (1.19, 1.33)	1.06 (0.76, 1.47)	0.6
Respiratory	1.67 (1.47, 1.91)	2.31 (2.05, 2.61)	1.36 (1.21, 1.53)	0.01	1.24 (1.02, 1.50)	1.32 (1.20, 1.45)	1.04 (0.73, 1.48)	0.5
Cardiovascular	1.65 (1.41, 1.92)	2.36 (2.13, 2.61)	1.38 (1.23, 1.55)	0.001	1.16 (1.03, 1.29)	1.22 (1.15, 1.30)	1.03 (0.71, 1.50)	0.7

Note:1st and 99th refer to the percentiles of daily mean temperatures in 29 urban counties (combined) and 60 rural counties (combined). MMT refers to the minimum mortality temperature for 29 urban counties (combined) and 60 rural counties (combined).

We estimated the mortality fraction (%) and death counts attributable to nonoptimum temperatures for urban and rural counties (Table 4). The attributable fractions of mortality and attributable death counts are much higher for cold effects (temperature below MMT) than hot effects (temperature above MMT). An estimated 7.0% (95% eCI: 3.9, 10.1) of all-cause mortality was attributable to cold effects for urban counties, while 0.8% (95% eCI: 0.4, 1.2) of all-cause mortality was attributable to hot effects. An estimated 16.4% (95% eCI: 14.9, 17.9) of all-cause mortality was attributable to cold effects, while 0.9% (95% eCI: 0.7, 1.1) of all-cause mortality was attributable to hot effects.

Table 4

Attributable fraction of mortality and annual attributable death counts (with 95% empirical confidence intervals) of cold and hot effects (nonoptimum temperatures) for urban and rural counties in Zhejiang Province, 2009–2015.

Cause-specific mortality and subgroups	Attributable fraction of mortality (%)				Attributable death counts (cases/year)
	Cold effects (temperature below MMT)		Hot effects (temperature above MMT)		Cold effects (temperature below MMT)		Hot effects (temperature above MMT)
	Urban	Rural	Urban	Rural	Urban	Rural	Urban	Rural
All cause	7.0 (3.9, 10.1)	16.4 (14.9, 17.9)	0.8 (0.4, 1.2)	0.9 (0.7, 1.1)	6,741 (3,739; 9,684)	32,801 (29,765; 35,758)	748 (384; 1,151)	1,878 (1,398; 2,197)
Nonaccidental	7.5 (4.1, 10.7)	17.6 (15.9, 19.1)	0.7 (0.2, 1.2)	0.9 (0.7, 1.1)	6,556 (3,598; 9,390)	31,540 (28,494; 34,228)	606 (176; 1,053)	1631 (1,254; 1,971)
Age 0–64	4.7 (−2.5, 11.4)	10.6 (6.4, 14.5)	0.1 (−0.4, 0.7)	0.1 (−0.1, 0.3)	870 (−461, 2101)	4,967 (2,990; 6,775)	24 (−74, 129)	37 (−47, 140)
Age 65+	10.7 (4.3, 16.3)	19.6 (17.7, 21.1)	0.6 (−0.1, 1.2)	1.0 (0.7, 1.2)	7,398 (2,981; 11,301)	33,368 (30,210; 36,014)	416 (−69, 832)	1639 (1195, 2048)
Males	10.1 (3.5, 15.8)	16.9 (15.0, 19.0)	0.4 (−0.2, 1.0)	0.6 (0.3, 0.8)	4,979 (1,729; 7,804)	20,994 (18,634; 23,603)	212 (−99, 494)	683 (373, 994)
Females	7.6 (2.5, 12.3)	18.5 (16.0, 20.8)	0.9 (−0.2, 1.8)	1.1 (0.9, 1.4)	2,900 (959; 4,719)	17,257 (14,909; 19,381)	326 (−77, 691)	1,062 (839; 1,305)
Cardiopulmonary	9.6 (5.4, 13.1)	22.1 (19.8, 24.1)	1.0 (0.2, 1.8)	1.2 (0.9, 1.4)	4,257 (2,405; 5,834)	21,197 (18,999; 23,125)	450 (89, 802)	1103 (864; 1,343)
Respiratory	8.1 (4.5, 11.5)	22.1 (18.2, 25.8)	0.9 (−0.3, 2.2)	1.4 (1.0, 1.8)	1,174 (655; 1,673)	7,160 (5,888; 8,347)	135 (−44, 320)	453 (324, 582)
Cardiovascular	11.5 (6.7, 16.0)	20.0 (17.4, 22.4)	0.9 (0.1, 1.6)	1.7 (1.3, 2.2)	3,448 (2,009; 4,797)	12,714 (11,067; 14,247)	255 (30, 480)	1,100 (827; 1,399)

Additionally, cause-specific mortality shows different optimum temperatures between urban and rural counties (Table S1). The MMPs for cause-specific mortality are lower in urban counties (from 49th to 71st) than rural counties (from 76th to 78th).

Stratified Urban–Rural Temperature–Mortality Relationships

The stratified analysis shows the overall associations between temperature and nonaccidental mortality varied by age and sex (Figure 3). The increased risk of nonaccidental mortality during extreme cold and extreme heat among males and females are both higher in rural counties than in urban counties, and the urban–rural disparity is larger among females than males (Table 4). For the elderly (), there is a distinct urban–rural gap of RRs for cold, which is 1.59 (95% CI: 1.36, 1.87) for urban counties vs. 2.18 (95% CI: 2.05, 2.32) for rural counties (Table 3). Furthermore, RRs for heat are also larger for the rural elderly (1.20; 95% CI: 1.15, 1.25) than urban elderly (1.11; 95% CI: 0.98, 1.26) (Table 3). In contrast, it is interesting to note that the apparent increase in mortality among people aged 0–64 only occurs in low temperatures, but not in high temperatures. There were a weaker urban–rural difference of cold mortality risks among people aged 0–64 compared with the elderly () (Figure 3).

Figure 3.

Pooled temperature–mortality associations along lag 0–21 d for nonaccidental mortality stratified by age and sex for urban and rural counties in Zhejiang Province, 2009–2015, with 95% confidence intervals (CIs). Note: The vertical lines represent the minimum mortality temperature (MMT, solid) and the 1st and 99th percentiles of the temperature distribution (dashed) for 29 urban counties and 60 rural counties in Zhejiang Province, 2009–2015. The shading lines represent the 95% CI areas for risk estimates. Distributed lag nonlinear models (DLNMs) were used to model the exposure–lag–response associations between temperature and mortality. A cross-basis function was defined using a quadratic B-spline with two internal knots of temperature and a natural cubic spline for the space of 21 lag days with 4 degrees of freedom. RR, relative risk.

Heterogeneity Test

County-specific associations between temperature and cause-specific mortality are shown in Figure S1. Results of the heterogeneity analysis are illustrated in Table 5. In the meta-analytical model without meta-predictors, the estimated heterogeneity () in the overall exposure–response associations for all-cause mortality between counties was 20.9% (Cochran Q test ). Adding the urbanization level to the model (Wald test ) decreased the heterogeneity to 16.8% (Cochran Q test ), suggesting that the modification effect of the urbanization level on temperature–mortality associations is statistically significant. The results for cause-specific mortality, and for the subgroup of elderlies, males and females, also show statistically significant Wald test p-values.

Table 5

Second-stage random-effects meta-analysis and meta-regression models of 89 county-specific results: Wald test on significance of urbanization level in explaining variations in relative risks for heat (99th percentile vs. MMT) and for cold (1st percentile vs. MMT), and overall cumulative temperature–mortality curves, Cochran Q test for heterogeneity, statistics for residual heterogeneity.

Cause-specific mortality and subgroups	I^2 and Q test	RR for heat (99th vs. MMT)			RR for cold (1st vs. MMT)			Overall temperature–mortality associations
		Intercept only	Urbanization level	Wald test (p-value)	Intercept only	Urbanization level	Wald test (p-value)	Intercept only	Urbanization level	Wald test (p-value)
All cause	I^2 (%)	30.5	22.9	<0.001	20.1	15.4	<0.001	20.9	16.8	<0.001
	Q test (p-value)	<0.001	0.004		<0.001	0.014		<0.001	0.005
Nonaccidental	I^2 (%)	25.7	18.9	<0.001	20.9	16.7	<0.001	17.4	14.2	<0.001
	Q test (p-value)	<0.001	0.003		<0.001	0.008		0.004	0.018
Age 0–64	I^2 (%)	16.3	14.7	0.884	1.0	1.0	0.118	7.3	6.4	0.139
	Q test (p-value)	0.033	0.018		0.456	0.513		0.163	0.180
Age 65+	I^2 (%)	23.8	17.3	<0.001	20.2	16.8	<0.001	18.6	15.8	<0.001
	Q test (p-value)	<0.001	0.007		<0.001	0.006		0.002	0.010
Males	I^2 (%)	12.3	8.7	0.003	4.9	1.0	<0.001	5.0	3.8	0.120
	Q test (p-value)	0.051	0.100		0.237	0.499		0.221	0.276
Females	I^2 (%)	18.8	12.1	<0.001	17.2	13.5	<0.001	7.9	3.6	<0.001
	Q test (p-value)	0.002	0.038		0.004	0.024		0.133	0.312
Cardiopulmonary	I^2 (%)	29.5	24.1	<0.001	19.4	14.9	0.019	9.3	6.8	<0.001
	Q test (p-value)	<0.001	<0.001		<0.001	0.014		0.090	0.168
Cardiovascular	I^2 (%)	18.1	14.3	<0.001	13.1	8.7	<0.001	6.8	5.2	0.037
	Q test (p-value)	0.003	0.016		0.026	0.106		0.166	0.238
Respiratory	I^2 (%)	21.9	17.9	0.004	9.5	8.0	0.028	12.2	10.1	0.061
	Q test (p-value)	<0.001	0.004		0.084	0.127		0.041	0.062

Note: RR, relative risk.

Potential Sources of Urban–Rural Disparity

There exist distinct gaps in the percentage of the elderly and primary industry employment, education, GDP, air conditioner, and hospital beds between urban and rural counties (; Table 2). Table S2 shows the results of the effect modification on nonaccidental mortality risks by county-level characteristics. Counties with higher percentages of elderly and primary industry employment, lower education, less GDP, fewer hospital beds, and fewer air conditioners have higher mortality risks related to both heat and cold. For instance, an IQR increase in percentage of people old (i.e., 3.3% increase) was associated with a 7% (95% CI: 2, 13) increase in the effect of cold temperature (99th percentile vs. MMT) on nonaccidental mortality (Table S2). The modification effects of the sex structure and the percentage of children are not observed. In addition, counties with lower mean temperature have higher mortality risks to both heat and cold. However, these county-level characteristics are highly correlated to the urbanization level (Tables S3 and S4), and meta-regression with multiple variables was checked to have collinearity issue according to VIFs (Table S5). Thus, it is hard to specifically explore their role in the urban–rural disparity.

Sensitivity Analysis

Sensitivity analyses were performed by changing the df for temperature, time, and lag structure, and varying the maximum lags. Similar exposure–response curves and MMTs were obtained (see Figure S2). Similar urban–rural differences were found when using the 95th percentile vs. MMT and the fifth percentile vs. MMT to calculate RRs for cold and heat (Table S6). We also performed a sensitivity analysis by adjusting for relative humidity, , and ozone. Although the variances of the RRs become larger because of the reduced study period of 2013–2015, the urban–rural differences of the estimated RRs remained distinct.

Discussion and Conclusions

To our knowledge, this is the first study conducted in a developing country that finds an urban–rural disparity in both heat and cold mortality risks. We found that mortality risks (RRs) associated with both cold and hot temperatures were higher in rural areas than urban areas, for all types of diseases, people , and both sex groups. When we considered the number of deaths and the temperature distribution together, the urban–rural disparity of the attributable death counts to nonoptimum temperatures was much more apparent.

For heat effects, this finding challenges the general assumption in previous studies generally conducted in developed countries that urban residents are at a higher temperature due to the UHI effect and hence higher risk to extreme high temperatures (Heaviside et al. 2017). Our result is consistent with an earlier study in Jiangsu Province, China, which found an urban–rural difference in heat vulnerability and a significant modification effect due to the urbanization level (Chen et al. 2016b). Our estimates of the urban–rural differences of all-cause mortality RRs at the 99th vs. 75th percentiles of temperature [1.15 (95% CI: 1.07, 1.22) vs. 1.18 (95% CI: 1.14, 1.23); see Table S7] are lower than the estimates in the Jiangsu study [1.26 (95% posterior interval (PI): 1.23, 1.30) vs. 1.43 (95% PI: 1.36, 1.50)]. Lower heat-related mortality risks in urban areas were also found in Hubei Province, China (Zhang et al. 2017b).

Mortality risks of temperature can be described as a function of exposure and vulnerability (IPCC 2012). Even though urban populations experienced higher temperatures than rural populations due to the UHI effect, heat mortality risks were also observed to be lower in urban areas due to lower vulnerability. Our findings suggest this urban–rural vulnerability disparity might be attributable to demographic and socioeconomic factors related to the urbanization level in China, including age structure, education, GDP, health care services, type of occupations, and air conditioners (Table S2). For instance, agricultural workers, generally living in rural areas, usually work outdoors and are directly exposed to extreme temperatures. Moreover, there are less air conditioners for rural residents than urban residents (1.21 vs. 1.66 per household), which could increase the urban–rural gap in vulnerability, as reported in a previous study (Chen et al. 2016b). Recently, a heat vulnerability map aggregating factors of age, socioeconomic status, social isolation, and air conditioning for Zhejiang Province also showed that urban communities are at a markedly lower level of vulnerability than rural communities (Hu et al. 2017).

In contrast, studies in the United Kingdom (Hajat et al. 2007), Germany (Gabriel and Endlicher 2011), and Greece (Katsouyanni et al. 1993) found the mortality risks at extreme high temperatures were higher in urban areas than nonurban areas. That is possibly because in most developed countries, urban and rural populations do not have the large differences in living standards and access to health care services that are observed in China, but urban areas in developed countries might experience extremely high UHI effects (Heaviside et al. 2017). Additionally, studies in the United States reported contradictory results. For instance, Madrigano et al. (2015a) found a higher mortality increase associated with heat in urban counties than in rural counties in the northeastern United States. However, heat-related emergency room visitation rates were found higher in rural settings compared with urban areas in North Carolina (Kovach et al. 2015), as reported in a study on all-cause mortality in Ohio (Sheridan and Dolney 2003). Moreover, a nationwide study in the United States reported that heat-related mortality for the most urbanized counties were as high as those for the most rural counties and that urban–rural differences vary by regions (Berko et al. 2017). These inconsistent results indicate a considerable heterogeneity in urban–rural differences in heat health effects regionally and nationally.

Urban–rural differences in cold effects have, thus far, received very little attention. Interestingly, we found that the urban–rural disparity in RRs for cold [for all-cause mortality, 1.47 (95% CI: 1.32, 1.62) vs. 1.98 (95% CI: 1.87, 2.10)] was much higher than RRs for heat [for all-cause mortality, 1.15 (95% CI: 1.07, 1.24) vs. 1.18 (95% CI: 1.14, 1.23)], which could be ascribed to the synergy between UHI effects and vulnerability patterns. During the winter, higher temperatures due to the UHI effect in urban areas enlarge the urban–rural disparity of the risk estimates attributable to the vulnerability gap. However, during the summer, higher temperatures and lower vulnerability have the opposite effect on heat-related mortality risks in urban areas, leading to a weaker urban–rural disparity in RRs for heat.

The potential sources of the urban–rural vulnerability disparity to cold are similar to those to heat (Table S2). However, it is interesting that there is a pattern of adaptation to heat stress in hotter counties but no adaptation to cold stress in colder counties. Lower MMP in urban counties compared with rural counties also supports the fact that urban residents are more acclimated to cold weather conditions. This result is inconsistent with previous studies that reported colder regions are more adapted to cold stress than hotter regions (Guo et al. 2014), and may be due to the combined modification effects of socioeconomic and demographic factors.

In order to comprehensively estimate the mortality burden attributable to temperatures in urban and rural counties, we also estimated the attributable fractions of mortality considering the temporal relationships between human exposure and RRs (Gasparrini and Leone 2014). We found a high proportion of 16.4% (95% eCI: 14.9, 17.9) of total rural mortality attributable to cold, which is far greater than the urban estimates of 7.0% (95% eCI: 3.9, 10.1), as well as the average estimates of 15 cities in China [10.4% (95% eCI: 8.8, 11.8)] from a previous study (Gasparrini et al. 2015). We further calculated the attributable death counts, which gives useful information about the urban–rural disparity of absolute temperature-related mortality burden. All the risk measures consistently reveal the need to address the large urban–rural disparity of temperature-related mortality risks.

As one of the most developed provinces in China, the urban–rural gaps for income, education, and access to health care in Zhejiang are smaller than in other less developed regions (e.g., western and central provinces) (Pan and Shallcross 2016; Cao et al. 2010). Therefore, the urban–rural disparity in temperature–mortality relationships might be higher in less developed regions than our estimates in Zhejiang, and need further investigation. Such phenomena might also be common in other lesser developed countries, where nonurban areas lag behind in socioeconomic development and are often overlooked in health risk management.

Current and future trends suggest that the global population is expected to increasingly reside in cities, particularly in developing countries such as China (UN 2014). With vegetated areas being replaced with impervious surfaces and buildings, urbanization will continually increase urban populations’ exposure to extreme heat due to substantial UHI effects (Feng et al. 2014; Yang et al. 2017). However, according to the results of our study, urbanization-induced improvement in socioeconomic status, access to health services, etc., might enhance the capacity of growing urban populations to adapt to nonoptimum temperatures and reduce health risks. It is especially important for cold effects because both the projected warming climate, increased UHI effect, and enhanced adaptation positively contribute to the mitigation of cold-related mortality risks in urban areas. The trade-off between the health benefits and health risks brought about by urbanization should be incorporated into future risk analyses but are beyond the scope of our article.

Our study indicated that temperature effects are stronger on cardiopulmonary mortality than other causes of death. Previous studies revealed people with preexisting disease were more likely to be affected by temperatures (Basu 2009; Guo et al. 2014). It should be noted that higher cardiopulmonary prevalence in rural China compared with urban China were previously reported (Wang et al. 2015; Chen et al. 2017). Moreover, cardiopulmonary mortality was 1.5-fold higher in rural counties than urban counties in our study region, which could also partly explain stronger temperature effects on cardiopulmonary mortality in rural counties than urban counties. Therefore, improvements in cardiopulmonary health in rural areas are very important for reducing the urban–rural gaps of temperature-related mortality.

Effective multisectoral measures and policies should be implemented to reduce the urban–rural disparity of temperature health impacts. First, compared with urban residents, rural residents have relatively poorer health care services and health literacy in China (Liu et al. 2007; Zhang et al. 2017a). More efforts should be made to narrow the urban–rural gap of the access to health care, such as increasing investment in health care facilities and health care professionals in rural areas. Improving rural people’s general awareness is also needed to prevent temperature-related deaths, particularly for the elderly. Secondly, rural dwellings in China often lack thermal heating and insulation (Figure S3) (Evans et al. 2014; Yang et al. 2010). The statutory standard concerning thermal comfort is highly recommended for rural houses (Gasparrini et al. 2017); unfortunately, this has only been enforced so far for urban apartments in China. Third, due to the absence of central heating across Zhejiang Province, both urban and rural residents use air conditioners or portable heaters as heating devices in winter. Due to lower income, rural households are more likely to fall into fuel poverty compared with urban households (Bouzarovski et al. 2012). Targeted measures, such as financial assistance for paying electricity bills, will help build rural residents’ resilience (Bouzarovski et al. 2012).

Several limitations should be acknowledged in this study. First, given the county unit of analysis, many counties have both urban and rural attributes but are classified by the entire county, due to the requirements of the sample size of death counts in time-series models. This may lead to a measurement error in the exposure–response associations for both urban and rural areas, and they are likely to underestimate the existing health disparity between them. Secondly, we performed the analyses at the county level rather than at the subregional climate zone level due to the unavailability of the full addresses of deaths. This limits our research to exploring the adaptation to temperatures by the local population at a coarse resolution. Third, exposure misclassification is a well-recognized inherent limitation of environmental epidemiological studies (Zeger et al. 2000). We used population-weighted temperature rather than individual temperature to measure population exposure, which may induce individual exposure measurement errors and is likely to underestimate the temperature mortality effects (Hutcheon et al. 2010). Additionally, this study was conducted in only one developed province in China. The urban–rural disparity in other regions, particularly in other less developed regions of China and in other low- and middle-income countries, remains uncertain and needs to be further investigated.

Nonetheless, taking advantage of high-resolution temperature data, our study explored the urban–rural differences of both hot and cold mortality effects based on the two aspects of the risk determinants: exposure and vulnerability. Although urban residents are generally at a higher temperature exposure than rural residents due to the UHI effect, the mortality risks to heat are lower for urban residents. Mortality risks, attributable fractions of mortality, and attributable death counts in relation to both cold and heat were found to be higher in rural counties than in urban counties of Zhejiang Province, China. A better understanding of possible urban–rural differences is critical for global and regional mortality risk projections in the context of climate change and urbanization. Past temperature-related health risk assessments may have overlooked the important heterogeneity across subregions and underestimated the mortality risks for the population as a whole. Therefore, future projections on climate health impacts should consider the urban–rural disparity in mortality risks. The urban–rural disparity also suggests that area-specific adaptation strategies, such as narrowing the urban–rural gaps in access to health care and awareness to risk prevention, should be developed, and emergency planning should be put in place to reduce temperature-related mortality.

Click here for additional data file.

Click here for additional data file.