Changing Susceptibility to Non-Optimum Temperatures in Japan, 1972-2012: The Role of Climate, Demographic, and Socioeconomic Factors.

BACKGROUND: Previous studies have shown that population susceptibility to non-optimum temperatures has changed over time, but little is known about the related time-varying factors that underlie the changes.
OBJECTIVE: Our objective was to investigate the changing population susceptibility to non-optimum temperatures in 47 prefectures of Japan over four decades from 1972 to 2012, addressing three aspects: minimum mortality temperature (MMT) and heat- and cold-related mortality risks. In addition, we aimed to examine how these aspects of susceptibility were associated with climate, demographic, and socioeconomic variables.
METHODS: We first used a two-stage time-series design with a time-varying distributed lag nonlinear model and multivariate meta-analysis to estimate the time-varying MMT, heat- and cold-related mortality risks. We then applied linear mixed effects models to investigate the association between each of the three time-varying aspects of susceptibility and various time-varying factors.
RESULTS: MMT increased from 23.2 [95% confidence interval (CI): 23, 23.6] to 28.7 (27.0, 29.7) °C. Heat-related mortality risk [relative risk (RR) for the 99th percentile of temperature vs. the MMT] decreased from 1.18 (1.15, 1.21) to 1.01 (0.98, 1.04). Cold-related mortality risk (RR for the first percentile vs. the MMT) generally decreased from 1.48 (1.41, 1.54) to 1.35 (1.32, 1.40), with the exception of a few eastern prefectures that showed increased risk. The changing patterns in all three aspects differed by region, sex, and causes of death. Higher mean temperature was associated (p<0.01) with lower heat risk, whereas higher humidity was associated with higher cold risk. A higher percentage of elderly people was associated with a higher cold risk, whereas higher economic strength of the prefecture was related to lower cold risk.
CONCLUSIONS: Population susceptibility to heat has decreased over the last four decades in Japan. Susceptibility to cold has decreased overall except for several eastern prefectures where it has either increased or remained unchanged. Certain climate, demographic, and socioeconomic factors explored in the current study might underlie this changing susceptibility. https://doi.org/10.1289/EHP2546.

Introduction

In response to projected climate warming, interest in the assessment of future health risks related to temperature has been growing. To evaluate future risks, knowledge about the association of temperature with mortality and how it has changed over time is critical (Arbuthnott et al. 2016). Numerous previous studies have suggested a U- or J-shaped association with increased mortality above or below a certain temperature threshold, which may take either a single value or a range of values (Basu 2009; Gasparrini et al. 2015b). This threshold is often referred to as the minimum mortality temperature (MMT) and is interpreted as the optimum temperature (Gasparrini et al. 2015b). Recently, several studies have reported changes in the temperature–mortality association over time (Arbuthnott et al. 2016). Most of these studies have indicated a reduction in heat-related risks (Bobb et al. 2014; Davis et al. 2003; Donaldson et al. 2003; Gasparrini et al. 2015a; Ha and Kim 2013; Ng et al. 2016; Petkova et al. 2014), and some have highlighted an upward shift of the MMT (Åström et al. 2016; Matzarakis et al. 2011; Todd and Valleron 2015). A few studies that have investigated the cold-related risk over time have showed conflicting results, with some reporting a reduction in risk (Carson et al. 2006; Christidis et al. 2010) and others reporting no change (Åström et al. 2013; Barnett 2007; Chung et al. 2017; Ekamper et al. 2009).

Although the evidence of temporal changes in temperature-related mortality is clear, little is known about the underlying factors that explain the changes. Identifying these factors is important to help elucidate vulnerable subgroups over time and improve the management of future risk. Previous studies have postulated many possible causes (e.g., general improvement in health care resources, increased awareness, and implementation of adaptation strategies such as air conditioning (AC), heating, and early warning systems), but quantitative evidence is very limited. Only a few studies have attempted to quantitatively attribute the changes in vulnerability to specific factors such as AC prevalence and other socioeconomic variables (Bobb et al. 2014; Ng et al. 2016); however, the results of these studies have limited interpretability because of analytical assumptions.

In this research, we investigated the changing susceptibility to non-optimum temperatures among the population in the 47 prefectures of Japan, covering the entire country for the period 1972–2012 (approximately four decades). Our examination focused on three aspects of the temperature–mortality relationship: time-varying MMT and heat- and cold-related mortality risks. To study the time-varying aspects, we used a two-stage analysis with a time-varying distributed lag nonlinear model (TV-DLNM) (Gasparrini 2011; Gasparrini et al. 2015a) and multivariate meta-analysis (MV-META) (Gasparrini et al. 2012). In addition, we evaluated the extent to which these time-varying aspects could be explained by time-varying explanatory variables (e.g., climate, demographic, and socioeconomic variables, and AC prevalence). To our knowledge, this is the first study to comprehensively assess the role of various time-varying factors in changing susceptibility to temperature in these three different aspects.

Methods

Data

We collected daily time-series data on mortality and weather variables for each of the 47 prefectures in Japan for the period 1972–2012 (except for Okinawa, which was 1973–2012). Figure S1 shows the map of all prefectures, and how we divided the country into eastern and western regions. The two regions were determined such that the time-varying patterns of the three aspects of susceptibility were similar within each region but distinct between the regions. The two regions also corresponded to the two groups identified by the longitude cutoff as 136.0 E and the latitude cutoff as 35.1 N (see Table S1). Daily mortality data were extracted from a computerized death certificate database maintained by the Ministry of Health, Labor and Welfare of Japan, including daily numbers of deaths from all causes and from cardiovascular [International Classification of Diseases (ICD)-8, 390–458; ICD-9, 390–459; ICD-10, I00–I99], respiratory (ICD-8 and ICD-9, 460–519; ICD-10, J00–J99), and noncardiorespiratory conditions, respectively. For weather variables, we derived daily mean temperature and mean relative humidity based on hourly measurements provided by the Japan Meteorology Agency for a single weather station in the capital city of each prefecture. Table S1 provides the summary statistics for all-cause mortality and daily mean temperature.

Different kinds of prefecture-level meta-variables were collected over the years of the study period to assess their influences on the time-varying temperature-related mortality risk. First, we collected climate variables because many previous studies have reported that climate condition is an important factor that can explain the heterogeneity of temperature-related risks (Chung et al. 2015, 2017; Gasparrini et al. 2015b; Guo et al. 2014). We obtained annual maximum, mean, and minimum of daily mean temperature and annual mean of daily relative humidity. Second, we collected demographic information because many studies have shown that temperature-related risks vary by age (Barnett 2007; Bobb et al. 2014; Chung et al. 2015, 2017); therefore, the population’s age distribution should be a major determinant of the population’s susceptibility to temperature. We included annual data on the proportion of the population in each prefecture (), 5–9 y (), and of age () obtained from the Statistics Bureau of the Ministry of Internal Affairs and Communications of Japan. Third, we obtained socioeconomic variables because previous studies have hypothesized that the temporal change in temperature-related mortality risk may have been derived by socioeconomic improvement at either the individual or community-level (Arbuthnott et al. 2016; Åström et al. 2016; Chung et al. 2017; Gasparrini et al. 2015a; Harlan et al. 2013; Ng et al. 2016; Todd and Valleron 2015). We obtained data on annual average savings for households with two or more persons (the sum of deposits at banks and other financial/nonfinancial institutions, premiums for life and damage insurances, and various types of securities) for 1974–2009 from the National Survey of Family Income and Expenditure, which is conducted every 5 y (Statistics Bureau of the Ministry of Internal Affairs and Communications of Japan 2009). Yearly data on the annual consumer price index (CPI), which measures fluctuations in the cost of goods and services over time, were obtained for 1972–2012 (Statistics Bureau of the Ministry of Internal Affairs and Communications of Japan 2015). Yearly data on annual economic power index (EPI) were obtained from the survey of local public finance for 2003–2010 (Ministry of Internal Affairs and Communications of Japan 2009). The EPI measures the financial strength of a local government, derived as an average of the ratios of basic financial revenues to basic financial needs in the past 3 y. Finally, we collected data on the annual prevalence of AC (1972–2009) for households with two or more persons from a regional statistics database (Asahi Newspaper Publishing 2015). Data for indicators of household heating were unavailable for the current study and were therefore not considered. Table S2 shows the summary statistics and Figure S2 presents the temporal trajectories for each meta-variable.

Modeling

The modeling framework and statistical procedures are briefly described below; mathematical and algebraic details are provided in “Details for statistical methods” in the Supplemental Material. For all computations, we used R statistical software (version 3.3.3; R Development Core Team).

Two-stage modeling.

In the first stage, we modeled the time-varying association between temperature and mortality in each prefecture using a generalized linear model with quasi-Poisson distribution. We included a natural cubic B-spline of time with 8 degrees of freedom per year to control for seasonality and long-term time trend, and indicator variables to control for the day-of-week effect. The model incorporated a TV-DLNM to describe a nonlinear and delayed association that changed over time. First, we defined a cross-basis (CB) for temperature and lag as in a time-constant DLNM (Gasparrini 2011, 2014): a quadratic B-spline for temperature with three internal knots placed at the 10th, 75th, and 90th percentiles of location-specific temperature distributions, and a natural cubic B-spline for the lag with an intercept and three internal knots placed at equally spaced values in the log scale. The lag period was extended to 21 d to capture the long-delayed association. The modeling specifics were based on the justification in a previous study (Gasparrini et al. 2015b). We then added a linear interaction between the CB and time to construct the TV-DLNM (Gasparrini et al. 2015a).

From the first stage, we obtained a set of 50 coefficients (25 for the main term of the CB and 25 for the interaction between the CB and time), and reduced them over the lag space into a set of 10 coefficients (5 for the main term and 5 for the interaction) for each prefecture (Gasparrini and Armstrong 2013). The reduced number of 10 coefficients represented the time-varying temperature–mortality association cumulated over the lags. Then, in the second stage, we applied MV-META to the prefecture-specific 10 coefficients to obtain pooled estimates for the whole of Japan and the prefecture-specific best linear unbiased predictor (BLUP) with the corresponding standard error matrices. Using the pooled estimate and BLUPs, we calculated the time-varying MMT and heat- and cold-related mortality risks, as described in the following section.

We conducted the two-stage analysis using all-cause mortality in the total population and stratified by region (east/west) and sex. Additionally, we conducted the analysis using cause-specific mortality in the total population. As a sensitivity analysis, we analyzed all-cause mortality in the total population with and without adjusting for daily mean relative humidity. The relative humidity was included in the first stage model as a CB using the same specifications as the CB for temperature except that the lag period extended to 15 d. The results were nearly the same with or without the adjustment, and we reported the ones without the adjustment.

Estimating the time-varying MMT and heat- and cold-related mortality risk.

Over a grid of every 10 d for the entire study period (1972–2012), we estimated the MMT and the heat- and cold-related mortality risks as follows. The grid starts from 1 January in 1972 and includes 1,498 d in total, with 10-d intervals. At each time point of the grid, we first reduced the 10 coefficients into a set of 5 coefficients, which represented the temperature–mortality association at the specific time point. Using the 5 coefficients, we calculated a point and an interval estimate of the MMT; we simulated 1,500 values of the MMT using Monte Carlo sampling, as proposed in a previous study (Lee et al. 2017), and we then calculated the sample mean as a point estimate and the sample percentiles (i.e., the 2.5th and 97.5th percentiles) as a 95% confidence interval (CI) for the MMT. We also calculated the minimum mortality temperature percentile (MMTP). At each time point of the grid, we derived a temperature distribution using the daily mean temperatures for the 6-y moving period with the specific time point as a center, and evaluated which percentile corresponds to the estimated MMT at the specific time point. Figure S3 shows the time-varying temperature distribution obtained using the 6-y moving period. The choice of 6-y moving period seemed most appropriate to capture the smooth change in temperature distribution without showing many of the short-term fluctuations.

Next, we calculated the heat-related mortality risk as the relative risk (RR) of mortality for the 99th percentile of the time-varying temperature distribution defined by the 6-y moving period versus the time-varying MMT. The cold-related mortality risk was calculated as the RR for the first percentile of the time-varying temperature distribution versus the time-varying MMT. In many prefectures, the MMT was estimated as exceeding the 99th percentile of the temperature distribution at later time periods. In such cases, we let the heat-related risk be null.

Investigating the relationship with prefecture-specific meta-variables.

We applied linear mixed effects models (LMEMs) to examine the association between the three aspects (MMT and heat risks and cold risks) and prefecture-specific meta-variables. For each aspect, we obtained the prefecture-specific estimates at the centering point (1 July) for the years when each meta-variable was available, and the corresponding standard errors. Depending on the meta-variable, the number of yearly measurements ranged from 8 to 41 (see Table S2 and Figure S2). Then, for each aspect and each meta-variable, we fitted the LMEM with the prefecture-specific random intercepts and with the prefecture-specific time trend adjustment, using a natural cubic spline with two internal knots placed at 1985 and 2000 to avoid temporal confounding (for EPI, we adjusted for linear time trends because the yearly measurements for EPI were available only after 2003). The residual serial and spatial correlations were accounted for using the first-order serial autoregressive structure and including the longitude and latitude of each prefecture in the model, respectively. The uncertainty for each aspect estimate was incorporated using the inverse of the squared standard errors as a weight. We estimated the association between each aspect and each meta-variable as the change in MMT and heat and cold risks (as log of RR) per unit change in each meta-variable. We conducted the LMEM analysis for the total population and by sex. In all models, we excluded Okinawa because of its outlying pattern. Statistically significant associations were determined with , to account for multiple testing.

Results

Table S1 shows that the temperature distribution is dispersed with lower values of the minimum and up to 75th percentiles in eastern prefectures than in western prefectures, although the 95th percentile and maximum were similar among all prefectures. Figure S2 shows that the annual mean and minimum temperatures were generally lower in eastern prefectures and gradually increased over the years, whereas the annual mean relative humidity decreased over time. Figure S2 also indicates that and consistently decreased, whereas increased over time, and annual savings and CPI increased at faster speeds in earlier period. The increasing pattern in AC prevalence was clearly different between the east and the west.

Figure 1A shows the time-varying MMT (black) and MMTP (blue) in Japan. Overall, the MMT and MMTP increased over the study period, at a higher increasing rate especially after 1990. The pattern of increase was, however, distinct between the eastern and western region; it was roughly linear in the eastern prefectures but nonlinear with an abrupt increase in the mid-2000s in western prefectures (see Figure S4 for region-specific pooled results). The prefecture-specific results presented in Figure S5 show that the MMT and MMTP increased in almost all prefectures. Figure 1B,C shows the prefecture-specific MMTs in the first and last year (1972 and 2012) of the study period. The MMT increased throughout the country between the 2 y, and it was generally higher in the western prefectures in both years.

Figure 1.

Time-varying minimum mortality temperature (MMT), minimum mortality temperature percentile (MMTP), heat- and cold-related mortality risks. (A) Time-varying MMT and MMTP in Japan; (B and C) prefecture-specific MMT in 1972 and 2012; (D) time-varying heat-related mortality risk in Japan; (E and F) prefecture-specific heat-related mortality risks in 1972 and 2012; (G) time-varying cold-related mortality risk in Japan; (H and I) prefecture-specific cold-related mortality risks in 1972 and 2012. The MMT is the daily temperature associated with the lowest daily mortality; heat-related mortality risk is the relative risk (RR) comparing mortality at the 99th percentile of temperature vs. the MMT; cold-related mortality risk is the RR comparing mortality at the first percentile of temperature vs. the MMT. Black solid and dashed lines in (A), (D), and (G) indicate the central estimate and 95% pointwise confidence intervals.

Figure 1D shows that overall heat risk in Japan decreased over the study period and was close to null at the end of the period. Figure S6 shows the generally decreasing risk in most prefectures (except for two western prefectures, Miyazaki and Okinawa). Despite the general decrease, a “hump” in the trajectory was visible in many prefectures and overall in Japan during the 1990s. In many prefectures, the risk nearly approached null at the end of the period, except for several eastern prefectures where it remained significant. Region-specific pooled results show that heat risk remained in the east but not in the west (see Figure S4). Figure 1E,F shows that heat risk decreased between 1972 and 2012, and it was generally lower in the western prefectures in both years.

Figure 1G shows that overall cold risk also decreased in Japan. The risk was higher in the 1970s and 1980s, but began to decline in the late 1980s before leveling off in the 1990s and 2000s. The temporal patterns differed between eastern and western prefectures (see Figure S7 for prefecture-specific results and Figure S4 for region-specific results); cold risk decreased in the western region, whereas it increased or roughly leveled-off in the eastern region. Nevertheless, in both regions, cold risk was still significant at the end of the period in all prefectures. Figure 1H,I shows that the higher cold risk in western prefectures in 1972 reduced to levels comparable with those of the eastern prefectures in 2012.

Figure 2 shows the lag-cumulative RR curves in four different years (approximately the midpoint in each decade) for Japan, Tokyo (an eastern prefecture with the largest population in Japan), and Fukuoka (a western prefecture with the third largest population in the west). These two prefectures were chosen because the region-specific overall pattern was well-reflected in the two prefectures. In Japan, both the right and left sides of the U-shape curve lowered over the decades and the MMT shifted slowly to the right. In Tokyo, the right side of the RR curve lowered but the left side was elevated over time. In Fukuoka, both sides of the U-shape curve lowered and the right side nearly flattened over time. Figure S8 shows the RR curves in the first and last years of the study period for each of the 47 prefectures. The right side of the curve lowered in most of prefectures, whereas the change of the curve on the left side was more variable.

Figure 2.

Lag-cumulative relative risk (RR) curves at four time points (1976, 1985, 1996, 2006) (i.e., approximate midpoints in every decade of the study period) in Japan (the entire country), Tokyo (an eastern prefecture), and Fukuoka (a western prefecture). The RRs represent the relative risk of mortality for daily mean temperature relative to the MMT for each location and time point. Solid lines with shaded areas indicate the estimated RR curves with 95% confidence regions.

Figure S9 presents the time-varying MMT and the heat- and cold-related mortality risks by sex. The figure shows that the MMT and MMTP increased and the heat risk decreased over time by larger amount for females than for males. The cold risk-change pattern was similar between males and females. Figure S10 displays the time-varying trends by different causes of death. The increase in MMT and the decrease in heat and cold risks were much larger for respiratory mortality than for cardiovascular mortality. For noncardiorespiratory mortality, the heat and cold risks were initially low and either decreased by only a small amount or remained at similar levels.

Table 1 shows that some climate variables were significantly associated () with heat or cold risk; higher mean temperature was associated with lower heat risk, whereas higher mean relative humidity was associated with higher cold risk. Demographic variables were strongly associated with MMT and both risks; higher levels of were associated with lower MMT and lower heat risk but with increased cold risk, and higher was associated with lower MMT and lower cold risk. Also, higher value of was associated with higher MMT and higher cold risk. Among the socioeconomic variables, higher EPI was associated with lower cold risk. AC was not strongly associated with any of the three aspects. Tables S3 and S4 show that some of the significant meta-variables differed between males and females. For example, a positive association between mean humidity and cold risk was shown only in females, whereas higher minimum temperature was associated with higher cold risk only in males. In addition, higher values of annual savings and EPI were associated with increased MMT and decreased cold risk only in females, whereas higher CPI was associated with increased cold risk only in males.

Table 1

Associations between meta-variables and each outcome (minimum mortality temperature, heat-related mortality, and cold-related mortality).

Meta-variable	MMT		Heat-related RR		Cold-related RR
	Slope estimatea (95% CI)	p-Value	Slope estimatea (95% CI)	p-Value	Slope estimatea (95% CI)	p-Value
Climateb
Tmax	−0.0004 (−0.0058, 0.0049)	0.8727	−0.0002 (−0.0005, 0)	0.0645	−0.0002 (−0.0004, 0.0001)	0.1812
Tmean	0.0052 (−0.0057, 0.0161)	0.3478	−0.0009 (−0.0014, −0.0004)	0.0006*	−0.0004 (−0.001, 0.0001)	0.1210
Tmin	0.005 (0.0009, 0.0092)	0.0173	0 (−0.0002, 0.0002)	0.6377	0.0001 (0, 0.0003)	0.1492
RHmean	0.0004 (−0.003, 0.0038)	0.8010	0 (−0.0001, 0.0002)	0.6240	0.0003 (0.0001, 0.0005)	0.0022*
Demographic, age (%)c
≤4y	−11.2218 (−17.5004, −4.9432)	0.0005*	−0.5718 (−0.8039,−0.3397)	0.0000*	0.3213 (0.1239, 0.5187)	0.0014*
5–9 y	−10.2994 (−16.2165, −4.3822)	0.0007*	−0.0964 (−0.3204, 0.1277)	0.3991	−0.3363 (−0.5171, −0.1556)	0.0003*
≥65y	8.7305 (4.2986, 13.1624)	0.0001*	−0.222 (−0.4897, 0.0458)	0.1041	0.4578 (0.2148, 0.7009)	0.0002*
Socioeconomicd
Savings	−0.0274 (−0.0689, 0.0142)	0.1965	0.0011 (−0.0013, 0.0034)	0.3860	0.0004 (−0.0015, 0.0023)	0.6797
EPI	0.7935 (−0.3674, 1.9543)	0.1826	−0.0212 (−0.042, −0.0003)	0.0480	−0.0807 (−0.1152, −0.0462)	0.0000*
CPI	0.0057 (0, 0.0114)	0.0499	−0.0001 (−0.0004, 0.0001)	0.2138	0.0002 (0, 0.0004)	0.0190
Air conditioninge
AC	0.0005 (−0.0009, 0.0019)	0.4836	0.0001 (0, 0.0001)	0.1152	−0.0001 (−0.0002, 0)	0.0727

Note: AC, air conditioning; CI, confidence interval; CPI, consumer price index (a measure of the cost of goods and services by household); EPI, economic power index (a measure of financial strength of a local government, higher value represents more strength); MMT, minimum mortality temperature; RHmean, average daily mean relative humidity; RR, relative risk; Savings, annual average savings per household with two or more persons; Tmax, maximum daily mean temperature); Tmean, average daily mean temperature; Tmin, minimum daily mean temperature.

.

Change in each outcome [MMT in degrees Celsius (); heat- and cold-related mortality in ln(RR)] per 1-unit increase in each meta-variable. Each meta-variable is scaled such that the slope estimate is not a value very close to zero.

Climate variables for 1972–2012 (annual values, ): Tmax, Tmean, Tmin, and RHmean.

Demographic variables for 1972–2012 (annual values for percent of population, ).

Socioeconomic variables: Savings (annual average, available every 5 y from 1974 to 2009, ); EPI (2003–2010, ); CPI (1972–2012, ).

Air conditioning prevalence for households with two or more persons, 1972–2009 ().

.

Discussion

In this study, we investigated temporal changes in three different aspects (MMT, heat and cold risks) of the temperature–mortality association in Japan over the four decadal period (1972–2012), and we examined the role of climate, demographic, and socioeconomic factors in the changing association. We found that MMT increased and heat risk decreased overall in the country. Cold risk decreased in general, except for in several eastern prefectures where the risk either increased or remained flat throughout the study period. Some of the studied climate, demographic, and socioeconomic factors were strongly associated with MMT as well as with heat and cold risks.

The MMT and MMTP increased over time in almost all prefectures of Japan. These results are consistent with previous studies that have reported an increasing trend in MMT and MMTP in European countries during the last decades (Åström et al. 2016; Ekamper et al. 2009; Todd and Valleron 2015). In our study, the time-varying MMTP was calculated based on the time-varying temperature distribution; therefore, the increased MMTP implies that MMT increased faster than the progressive shift in temperature distribution. The observed increasing patterns in the MMT and MMTP differed between the eastern and western prefectures; the MMT took a longer time to shift upward in the western region, possibly because of the already higher MMTs during the 1970s and 1980s compared with the eastern region (see Figure S4).

The heat-related mortality risk in Japan decreased over the whole study period. The present results were consistent with those of many previous studies conducted in other locations (Bobb et al. 2014; Chung et al. 2017; Davis et al. 2003; Donaldson et al. 2003; Gasparrini et al. 2015a; Ha and Kim 2013; Ng et al. 2016; Petkova et al. 2014). A major driver for the reduction was shown to be an MMT increase that was faster than the increase of the 99th temperature percentile, leading to a narrowing gap between the two temperatures over time. More importantly, another driver was the decreased slope above MMT over the past decades. Figure 2 displays the RR curves estimated in four different years of the study period, where the slope above MMT decreased over time in Japan, Tokyo, and Fukuoka. Figure S8 shows similar patterns in almost all prefectures. Our results imply that population susceptibility to heat decreased in both contexts—with the threshold and the slope above the threshold.

Despite the generally decreasing heat risk, a temporary increase was observed during the 1990s in a number of prefectures and overall in Japan. These humps in the trend were induced by jumps in the 99th temperature percentile during the 1990s (see Figure S3) in the corresponding prefectures, and heat risk was reduced to a decreasing trajectory in tandem with the stabilized 99th percentile after the 1990s. At the end of study period, the heat risk reached approximately the null value in many prefectures, whereas it remained significant in several eastern prefectures. This suggests that the mechanism of adaptation to heat possibly differs between the two regions. For example, in the eastern prefectures with a cooler climate, populations are less exposed to hot weather and may be less well adapted physiologically to heat. Moreover, the generally lower AC prevalence in the east (see Figure S2) might reflect a lower perceived need to cope with heat, and different strategies to manage heat exposure might have contributed to the observed risk difference between the regions.

Cold risk generally decreased in Japan, but in several eastern prefectures the risk either increased or remained about the same, reiterating the results of previous studies that also varied depending on location (Åström et al. 2013; Barnett 2007; Carson et al. 2006; Christidis et al. 2010; Chung et al. 2017; Ekamper et al. 2009). A major driver for the cold risk change was shown to be the changing slope below the MMT. Figure 2 shows the slope below MMT increased in Tokyo, whereas it decreased in Fukuoka over time. Figure S8 also shows that in most southern prefectures, the slope decreased between the first and last years (1972 and 2012), whereas it either increased or remained the same between these two time points in many eastern prefectures. In addition to the slope, other factors like the MMT movement and first temperature percentile change following the warming climate might underlie the decrease or increase in cold risk. Regardless of the direction of change, the cold risk was still significant even in very recent periods.

Our results showed that higher mean temperature was associated with lower heat risk and higher minimum temperature was associated with higher cold risk in males, which is consistent with the results of previous studies (Chung et al. 2015, 2017; Gasparrini et al. 2015b; Guo et al. 2014), supporting the notion that populations in warmer climate are less vulnerable to heat and more vulnerable to cold. Interestingly, higher relative humidity was associated with higher cold risk. One recent study reported a similar finding that the cold-related risk of cardiovascular mortality was increased with higher humidity in China (Zeng et al. 2017). Many studies have investigated humidity as an independent exposure or a confounder for temperature but little research has focused on the joint effect of temperature and humidity. Our results suggest that humidity may be a potential modifier of temperature-related risks; further research is merited to confirm this hypothesis.

Demographic structure were shown to be important determinants of the different aspects of susceptibility. Lower heat and cold risks were associated with a higher proportion of children, which may imply that if a population is relatively young, it is at lower risk of mortality associated with heat and cold. However, higher proportion of infants/preschoolers () was associated with increased cold risk, which may be explained by the strong seasonality of infant mortality with peak in winter because of infections such as influenza and pneumonia. Greater proportions of elderly adults showed links with higher MMT and higher cold risk. All these findings suggest that future change in a population’s age distribution must be carefully considered in health risk assessment of climate change. In societies like Japan, aging has been progressing, which suggests increased susceptibility of the population. This social phenomenon may have dampened the reduction of cold risk attributable to factors other than age over time. Socioeconomic factors were also found to be related to MMT and cold risk. Interestingly, a sex difference was observed in significant socioeconomic variables. Higher annual savings and EPI were linked with increased MMT and decreased cold risk only in females, whereas higher CPI was associated with increased cold risk only in males. Based on these findings, we hypothesize that females are more likely than males to utilize their socioeconomic resources to cope with health risks; this hypothesis requires further investigation in future research.

We did not find any association between the heat risk and AC prevalence, which was consistent with previous studies (Bobb et al. 2014; Ng et al. 2016). A potential reason has been provided that the prevalence may not reflect the actual usage of AC owing to high electricity costs. We also did not find a strong association between heat risk and socioeconomic indicators, which is counterintuitive because the heat risk and the socioeconomic indicators consistently decreased and increased over time, respectively. Although the two time-series indicate a correlation under a long-term framework, our analysis focused on a short-term association (i.e., the current year), adjusting for the long-term trend. More refined analysis considering any delayed or temporally changing association may reveal significant associations between improvement in socioeconomic indicators and a reduction in heat risk.

Our study has several methodological strengths. First, previous studies that relied on a subperiod approach (e.g., annual or decadal analysis or comparison between two periods) to study temporal changes in the temperature–mortality association were unable to capture continuous change and had reduced statistical power as a result of data stratification. We used time as a modifying factor covering the entire period of data in a unified analysis, which is useful to estimate a smooth change as well as to increase power. We adopted the TV-DLNM proposed in a previous study (Gasparrini et al. 2015a) and further enhanced the second stage meta-analysis by pooling the coefficients for both the main and interaction terms, to avoid conducting meta-analyses at multiple time points. In addition, previous studies used fixed thresholds or temperature values (e.g., the average MMT or temperature percentiles for the whole study period) (Chung et al. 2017; Gasparrini et al. 2015a), thereby disregarding the possibility of changing MMT and the reference temperature percentiles over time. We applied a time-varying definition for heat- and cold-related risks, and calculated the risks based on time-varying MMT and temperature percentiles. Finally, we used an LMEM to investigate the association with meta-variables, which allowed for the use of the information of spatio-temporal variations in each aspect of susceptibility and each meta-variable.

We acknowledge several limitations in our study. First, the spatial unit was prefecture, which is an administrative unit with relatively large resolution. A more refined resolution would allow for finer investigation of spatial variation. Second, each meta-variable was investigated separately in the LMEM. Multiple variables can be included simultaneously in the model, but in such multi-variable modeling, there is a potential multi-collinearity and interpretation may become complex. Future research should consider using a factor-analysis approach to identify a few meaningful factors from the highly correlated meta-variables (e.g., one factor for each category of climate condition, demographic structure, and socioeconomic status) and include those factors simultaneously in the model. In addition, we assumed the association between each aspect and each meta-variable does not change over time and space for simplicity, and we did not consider any delayed association. In future research, the role of meta-variables can be more carefully examined for delayed and spatio-temporally varying influence. Finally, cause-specific mortality risks can be investigated further, linked with different factors.

Our results have several important implications for public health policy and practice. First, whereas the population showed it has adapted to heat and heat risk has approached the null toward the end of the study period, cold risk remains significant; therefore, more public health resources should be allocated to manage the cold risk. Second, the observed sex differences imply that health risk management should be strategized differently by sex. In particular, the negative association between socioeconomic variables and cold risk observed only among females suggest that a linking strategy may be needed between socioeconomic resources and the health risks of weather stressors for males. Also, management of risk that recognizes region-specific factors related to climate and demography might be more effective. Finally, the assessment of future human health risks under a changing climate should consider the changes in the temperature–mortality association, demographic structure, and socioeconomic situation given that the current projection is based on limited assumptions of no change in these components, which may lead to underestimation.

Conclusions

Our findings suggest that the population of Japan has adapted to heat over the last four decades, with some differences by regions, sex, and causes of death. Our results also suggest that the cold risk has decreased overall in Japan, although it has increased or remained flat in several eastern prefectures. Climate, demographic, and socioeconomic factors seem to have different roles that affect population susceptibility to non-optimum temperatures. The assessment of future health risks under projected climate change should take into consideration the potential impact of the different factors.

Click here for additional data file.