Optimizing future mortality rate prediction of extreme temperature-related cardiovascular disease based on skewed distribution in peninsular Malaysia

The impact of extreme heat events on human health is significant and diverse (Crimmins et al. 2016). Health is most vulnerable to extreme heat events (Patel et al. 2022). According to the Center for Disease Control and Prevention of United States, extreme heat events caused by high temperatures are the most prominent cause of heat-related human fatality, responsible for more deaths annually than hurricanes, fires, tornadoes, floods, and earthquakes (Vaidyanathan et al. 2020). Such events are expected to occur more often and are predicted to last longer because of global climate change. The World Health Organization (WHO) predicted that there will be almost 92,000 deaths annually because of heatwaves by 2030 with sub-Saharan Africa, Latin America, and southeast Asia bearing the largest burdens and extreme heat events causing more than 150,000 deaths globally from 1990 to 2018, with 2018 recorded as the warmest, particularly between June and July, through Asia, Europe, and North America, hospitalizing thousands of people with heat-related illness and recording more than 700 deaths in the immediate after-math (WHO 2019). The Intergovernmental Panel on Climate Change (IPCC) stated that southeast Asian countries and developing countries such as Malaysia will be at the greatest risk from the emergence of extreme heat events (Stocker 2014). CVDs have been recognized as the leading cause of death throughout the world. CVDs, also known as circulatory diseases, are categorized as heart and blood vessel disorders (Patel et al. 2022). According to the Department of Statistics Malaysia (DOSM), coronary heart disease (CHD), one of the CVDs, was the principal cause of death in Malaysia indicating an increasing trend between 2016 and 2019 (DOSM 2017, 2018, 2019, 2020). CVD is sensitive to climate conditions; in particular, temperature-related CVDs have become a growing public health concern. Currently, the growing CVD mortality because of global warming has become a major concern, hence, an understanding of the current pattern of seasonal variations assists in estimating the impact of future climate change. Future research in the CVD mortality rate is believed to provide valuable inputs to the stakeholders for taking precautions and preventive measures to flatten the curve trend. Health-impact projection is conducted by combining estimated health-exposure relationships with projected future exposures while considering various sources of uncertainty in the analysis. Many studies found a positive correlation between the CVD mortality rate and extreme temperature. The largest number of deaths in the United States occurred because of an increase in the summer temperature (Kouis et al. 2019). The relationship between temperature and mortality in the 2020s and 2050s was predicted by McGeehin & Mirabelli (2001). There is a positive correlation between the temperature and CVD mortality rate (Huynen et al. 2001; Hajat et al. 2002; Kynast-Wolf et al. 2010; Huang et al. 2012). Similarly, exposure to high temperatures increases the risk of CVD mortality, particularly during summers (Moghadamnia et al. 2017). An increment of 1 °C in temperature is associated with a 2–5% increase in CVD mortality rate, particularly for the elderly. Paravantis et al. (2017) found a significant correlation of daily CVD mortality rate in people over 65 years with minimum and maximum temperatures. At a remarkably high temperature, the CVD mortality rate increased by 20–35% with a 99.9% significance level. Jackson et al. (2010) estimated the future number of deaths because of heat events by considering three different future scenarios. The temperature elevated risk is apparent for CVDs and respiratory causes. For instance, the estimated number of deaths under low climate scenario yielded that the greater Seattle area may expect 68, 89, and 107 excess deaths in 2025, 2045, and 2085, respectively, because of CVDs and respiratory causes among people over 45 years. Hence, there is a strong correlation between heat events and elevated risk of mortality. An estimated rise in mortality rates with an increase of 1 °C in temperature above the heat level was observed for all causes of death, including CVDs (Hajat et al. 2014).

Baaghideh & Mayvaneh (2017) showed a strong positive association between maximum temperature and future CVD mortality rate. An increment of 1 °C is associated with a 4.27% increase in the CVD mortality rate in the following decades. Li et al. (2018) studied the projections of the burden of CVDs attributable to future temperatures. Future changes in climate are likely to cause an increment in heat-related CVDs. The monthly analysis showed that the most significant increase will occur in the summer, with percentage changes greater than 150% in the 2050s and up to 300% in the 2070s. Huang et al. (2019) further projected the relationship between future temperature and CHD mortality in the elderly. The annual temperature-related CHD mortality in the 2050s and 2070s was greater than the baseline. Liu et al. (2019) projected heat-related mortality for all causes of death including CVDs in the 2030s, 2060s, and 2090s relative to the baseline period. The projected heat-related mortality would cause a substantial increase in the majority of the GCMs with population growth, adaptation, and ageing. Gu et al. (2020) projected temperature-related mortality of CVDs under different climate scenarios with that of CHD burden showing an increase from 4.40 to 6.02% primarily because of people having heat-sensitive diseases.

Projection of climate change impact including health impact is often studied using the GCM (Noor et al. 2018). However, the GCM outputs cannot provide reliable information on spatial scales below 200 km, and therefore, cannot be used for local impact studies (Amirabadizadeh et al. 2016; Noor et al. 2018; Adeyeri et al. 2020). Additionally, there are several errors and uncertainties in GCM outputs (Amirabadizadeh et al. 2016; Noor et al. 2018; Adeyeri et al. 2020). To assess climate change on local scales, an approach known as downscaling is used to bridge the gap between the resolution of GCM outputs and the local climatic process (Noor et al. 2018). Murphy (1999) concluded that both dynamical and statistical downscaling exhibit similar skill levels and produce hydrologically plausible results with the former being more sensitive to climate change with respect to future climate scenarios. However, these two approaches can nonetheless differ significantly in their projections of future climate conditions. Murphy (1999) and Wilby et al. (2000) also highlighted the need to bias correct GCM outputs to assure meaningful results in meteorological applications. This has been supported by Shrestha et al. (2014) who stated that there is good agreement between statistically downscaled GCM simulated temperature with the observations, whereas the dynamically downscaled simulated results differ substantially from observations. However, the differences between the 2050s and 1970s obtained from the two approaches are qualitatively similar. Shaaban et al. (2011) concluded that dynamical downscaling can produce more reliable results than statistical downscaling as it does not incorporate the time evolution of dynamical features on the climate at continental and global scales. However, Wilby et al. (2000) argued that statistical downscaling has shown good performance to reduce biases for temperature and it has the advantage of efficiently producing ensembles of surface climate variables given by a limited set of gridded predictor variables. In contrast, producing dynamically downscaled outputs was the most time consuming and computationally demanding. Despite the higher level of sophistication and physical realism of the model, dynamical downscaling was not generally as efficient as statistical downscaling and does not depend on static regression relationships among physical processes. Overall, statistically downscaled GCM outputs remain the preferred approach provided that biases are explicitly accounted for. To date, there are no truly objective techniques for identifying which particular downscaling approach is the best or most competent in generating future climate projections for a given region and climate variables (Conlon et al. 2016; Smid & Costa 2018; Bettolli et al. 2021; Wang et al. 2022). Various aspects of the atmospheric system, which each model handles differently, may be altered as the climate changes. Conlon et al. (2016) and Jose et al. (2022) suggested that to assess the best climate projection, an ensemble mean value of GCMs projection needs to be produced. In addition to considering the mean change across a set of models, it is useful to identify the most and least extreme projections within a suite of climate-model outputs. Smid & Costa (2018) suggested that the most cutting-edge approach to provide future localized climate information is to combine dynamical downscaling with further statistical downscaling when assessing the climate evolution. However, the skills, expertise, access to funding and time computational power resources were arguably outstanding.

Future climate change exposures are typically derived using simulations from computationally expensive climate models. Even though these models reflect state-of-the-art knowledge on the climate system, their outputs are known to exhibit complex spatial–temporal biases when compared to observations. Factors contributing to this bias include errors in parameters describing physical and chemical processes, incorrect representation of the underlying processes with mathematical equations, and discretization of meteorological fields in space and in time. Hence, climate-model simulations for the projection period need to be bias-corrected prior to estimating future health impacts (Holthuijzen et al. 2021). BCM is renowned in statistical downscaling because it can reduce biases between GCM outputs and observations, and it skillfully reduces the coarse resolution of GCMs outputs into a finer resolution of observation scale (Teutschbein & Seibert 2012; Cannon et al. 2015; Fang et al. 2015; Ringard et al. 2017; Hong et al. 2022; Ngai et al. 2022). Wood et al. (2004) employed BCM to reduce the statistically downscaled biases of GCM outputs. BCM was successful in reproducing the main features of the observed hydrometeorology from the retrospective climate simulation, when applied to statistically downscaled GCM outputs, and competent to preserve the statistical moments of the hydrological series (Sennikovs & Bethers 2009). Generally, BCM for future simulation is accomplished in two steps. First, the bias between observations and simulations during the historical/estimated period is assessed. Then a correction algorithm is applied to future simulations by assuming the bias can be extrapolated to future periods. Olsson et al. (2015) stated that quantile mapping was able to provide temperature data that were sufficiently close to observed discharges in the control period and produced more realistic projections for mean annual and seasonal changes compared to the uncorrected GCM data. Putra et al. (2020) have found that BCM could reduce biases in the GCM outputs relative to observations with a coefficient of determination of 0.81 for spatial distribution. Tadese et al. (2020) also showed that BCM and variance scaling performed well in correcting the biases in GCM outputs despite the corrected maximum and minimum temperature being slightly overestimated for the mean and standard deviation. Despite the fact that BCM has shown considerable skill in reducing the biases, there are a few limitations that need to be addressed. Firstly, the remaining uncertainty regarding how well a calibrated BCM performs for conditions different from those used for calibration. Although a good performance of BCM during the calibration period has been shown, that does not guarantee a good performance under changed future conditions due to the stationary assumption of BCM. Secondly, BCM is less efficient in downscaling both high and low temperature extremes (Maurer & Hidalgo 2008; Abatzoglou & Brown 2012). The most commonly used BCMs were based on shifting or scaling climate-model simulations that have been shown to perform poorly for removing bias at extremes (Räisänen & Räty 2013). Thus, BCM has a tendency to overfit on calibration data, especially at extremes where data are scarce and highly variable (Piani et al. 2010; Grillakis et al. 2013; Lafon et al. 2013; Mamalakis et al. 2017; Holthuijzen et al. 2021). The limitation of BCM is mainly because it assumes that the temperature follows the normal distribution (Hempel et al. 2013; Gaitán et al. 2019; Lange 2019). Some studies recommended the normal to be the most appropriate distribution for describing the statistical properties of extreme temperatures (Zhan et al. 2020). However, Semenov (2008) concluded that normal is not a good fit. The means of yearly maxima for daily maximum temperature were reproduced less accurately compared to those outside the 95% confidence interval (CI) of observed means. The study suggested replacing the normal distribution with a more flexible distribution. Similarly, Barrow & Hulme (1996) proposed that normal distribution could yield accurate results when variations in the temperature were less than 0.1 °C, but produced errors for variations greater than 10 °C. BCM with gamma (GAM) distribution further decreased the simulation biases. Similarly, Pastén-Zapata et al. (2020) tested BCM in simulating the temperature and observed that the existing BCM with normal distribution did not improve the representation of daily temporal variability, especially the extremes. These findings are significant for tackling the main part of the statistical modeling process which is the selection of appropriate distribution for the BCM that could be overlooked in the past research. Various studies were conducted to fit a suitable statistical distribution on extreme temperatures and different conclusions were drawn depending on factors such as the area and duration of study. Particularly, the right distribution may enhance the BCM's performance in simulating the extremes. Several studies suggested GAM and GUM owing to their capability to capture the mean and variance of daily extreme temperatures (Folland et al. 1999; Jones et al. 1999; Horton et al. 2001). Other studies concluded that Weibull (WEI) is a good fit for both maximum and minimum temperature data (Seimela 2021), whereas several studies fitted daily extreme temperatures with generalized extreme value (GEV) distribution (Kharin & Zwiers 2005; Underwood 2013; An et al. 2020). Wilson (2016) revealed that the GEV distribution was better in modeling the extreme temperature behavior based on the lowest values of the Akaike information criterion (AIC) and Bayesian information criterion (BIC). Two comparative tests, namely, Anderson–Darling and Kolmogorov–Smirnov, confirmed that the GEV distribution is adequate for the data. Similarly, in Malaysia, Amin et al. (2018) concluded that Frèchet (FRE) and GEV are a better fit for monthly maximum temperature. GEV was observed to be suitable for a daily and annual maximum of daily maximum temperature data (Salleh & Hasan 2018). Results indicated that the nonstationary model is preferable for Kuala Terengganu, Muadzam Shah, and Senai stations. Hasan et al. (2013) and Hasan et al. (2014) fitted the annual extreme temperatures in Malaysia with the GEV distribution. Additionally, nonstationary GEV has been observed to perform better than the stationary GEV for several stations such as Subang, Bayan Lepas, Kota Kinabalu, Sibu, and Labuan. Hasan et al. (2014) applied cluster analysis on the studied stations, thereby revealing similarities that indicate the best performance of the nonstationary GEV distribution. Contrarily, Supian & Hasan (2021) compared several distributions to find the best distribution for fitting the annual extreme temperatures. The results indicated that a majority of stations preferred the generalized skew logistic distribution as the best fitted probability distribution. Some studies indicated excellent performances of the normal, LGNORM, and WEI distributions to describe the annual extreme temperature. In summary, no single distribution can be concluded to perform the best. Therefore, it is crucial to thoroughly analyze various aspects while selecting the most appropriate distribution for the extreme temperatures to achieve an accurate projection of the CVD mortality rate. Thus, in this study, the two issues in BCM will be solved by incorporating the skewed distribution to fit the extreme values of annual temperatures and to include the linear covariate to take into account nonstationarity (i.e., linear trend) in the data series. Hence, the main objectives are to identify a better statistical distribution among the skewed distributions (i.e., LGNORM, GAM, Pearson type 3 (P3), GEV, GUM, WEI, and FRE distributions), formulate the BCM-QM_skewed with better statistical distribution and linear covariate in reducing the biases between the historical/estimated and observed annual extreme temperatures in peninsular Malaysia and project the mortality rate of CVD based on annual extreme temperatures projection in peninsular Malaysia using the annual attributable death (AAD) equation.

In this study, observed extreme temperatures (1976–2005), historical/estimated extreme temperatures (1976–2005), future extreme temperatures (2006–2100), and observed mortality (1990–2017) data recorded annually were used. The secondary data of observed daily temperature data series between 1976 and 2005 in peninsular Malaysia were obtained from the Malaysian Department of Irrigation and Drainage (Wong et al. 2011). Then, the annual extreme temperatures were extracted using the block maxima approach. Similarly, the historical/estimated (1976–2005) and future (2006–2100) GCMs of annual extreme temperatures were obtained under moderate and worst-case future scenarios of RCPs known as RCP4.5 (moderate case) and RCP8.5 (worst case). The model for interdisciplinary research on climate (MIROC5) under the Coupled Model Intercomparison Project 5 (CMIP5) with resolution 1.406° × 1.4° (Williams et al. 2012) was used. The definitions of both RCPs are presented in Appendix 1. The secondary data of annually observed mortality rate of CVD between 1990 and 2017 in peninsular Malaysia were obtained from the Global Burden of Disease Collaborative Network (http://ghdx.healthdata.org/gbd-results-tool) (Roth 2018). This database is openly available and can be downloaded from the global health data exchange web portal.

Model 1 (M1): and are constants (i.e., stationary model)

Model 2 (M2): and are constants (i.e., nonstationary model)

Model 3 (M3): and are constants (i.e., nonstationary model)

Model 4 (M4): is constant (i.e., nonstationary model)

This study focuses on clustering the stations to study the impact of mortality rate because of CVD in different regions (Jolliffe & Philipp 2010). The identification of clusters is based on the similarity of distributions and states (see Appendix 6). This study defines region 1 as Perlis, region 2 as Kedah, region 3 as Pulau Pinang, region 4 as Kelantan and northeast of Perak, region 5 as Perak, region 6 as Pahang, region 7 as Terengganu and northeast of Pahang, region 8 as southwest of Perak, region 9 as Kuala Lumpur, Malacca, Negeri Sembilan and Selangor, and region 10 as Johor. Cluster analysis can potentially contribute to the health sector by identifying key targets for these applied initiatives, thereby enabling the allocation of resources in an optimally effective manner and designing efficient adaptation and mitigation strategies in future (Clatworthy et al. 2005). The minimum mortality temperature is 28 °C. RR values under RCP8.5 are higher than those under RCP4.5 for all regions. CI describes the inherent uncertainty in this estimate and the range of values within which the true effect lies with reasonable certainty. The effect size is precisely known for small intervals, whereas the uncertainty is greater for larger intervals, i.e., large sample size tends to give more precise estimates of effects, and a small CI corresponds to a small sample size. The projected estimated extreme temperatures are 30–32 °C and 33–42 °C under RCP4.5 and RCP8.5, respectively. Arias et al. (2021) stated that there is a high confidence that the temperature change might increase to 1.5 °C and 4 °C toward the late century under RCP4.5 and RCP8.5, respectively. The RR values produced lower uncertainties (high confidence) for extreme temperatures between 30 and 32 °C for regions 1–4, 7–8 under both RCPs, whereas RR values produce lower uncertainties (high confidence) for regions 5–6, 9–10 with extreme temperatures between 33–36 °C and 33–42 °C under RCP4.5 and RCP8.5, respectively. Notice that the RR value for the point estimate of 61 °C produced high uncertainties with very wide interval estimation (extremely low confidence). Thus, the point estimate can be assumed to be an imprecise estimate (Appendix 7). Notably, regions 5–6 and 9–10 project higher extreme temperatures compared to other regions because of the anthropogenic climate change in urban areas caused by human activities, such as burning fossil fuels including coal, oil, and natural gas which leads to severe environmental damage linked to climate change. The high warming rate, particularly in large areas such as Kuala Lumpur, Selangor, and Johor could be enhanced by the urban heat island effect. The short period of data coverage and higher frequency of the El-Niño Southern/Oscillation (ENSO) occurrences in the last few decades may contribute to higher warming rates in some Malaysian stations. Some studies based on shorter data records of the last four or five decades reported higher warming rates (Tangang et al. 2007). Furthermore, higher projected extreme temperatures in regions 5 and 6 (⁠36 °C) could also be influenced by geographical factors because both regions are situated near the Indian and Pacific oceans in which the inter-annual variability variation in sea surface temperature could influence the large-scale climate variability, such as ENSO and Indian Ocean Dipole (Isa et al. 2020).

This study developed BCM-QM_skewed by incorporating the skewed distribution and a linear covariate into BCM to better fit the extreme temperatures in peninsular Malaysia. LGNORM and GUM are the most appropriate distributions to model the annual extreme temperatures at stations throughout peninsular Malaysia. LGNORM and GUM are the most suitable models to fit the annual extreme temperatures in the western and eastern parts of peninsular Malaysia, respectively, based on the AIC and BIC diagnostic tests. It can also be inferred that a BCM-QM_skewed with time-varying covariate is a good fitted model for annual extreme temperatures at most stations across peninsular Malaysia which suggests that the trend of extreme temperatures has been linearly increasing during 1976–2005. BCM-QM_skewed exhibits considerable skills in reducing the biases from the GCM model for all stations. The corrected annual extreme temperatures are comparable to the observed values based on the RMSE, MAE, PBIAS, and PP values. The projection trends for extreme temperatures increase for both RCPs. Higher projection is more pronounced under RCP8.5 with precise estimates ranging between 33 and 42 °C compared to that under RCP4.5 with precise estimates ranging from 30 to 32 °C. The CVD mortality rate falls in the range of 0.002–0.014. Regions 5, 8–9 should be given extra attention because the trend of extreme temperatures has been shown to increase yearly, which may cause higher CVD mortality rates in future, with a median rate that can attain a value of 0.014 (i.e., 14 deaths per year) . Overall, the CVD mortality rate across all regions in peninsular Malaysia is significant, specifically under RCP8.5. However, high CVD mortality rates coincide with high extreme temperatures. The significance of this study rests on the fact that the health sector in Malaysia heavily depends on climatic conditions. Thus, studying future scenarios may provide important information to researchers and local governments, to propose appropriate adaptation measures and increase resilience. Furthermore, appropriate policies and plans can be sanctioned to prepare the public for changes because of extreme temperatures, particularly, on CVD mortality. Using the accurate predictions of future mortality rates of CVD, the government can play a crucial role by enacting policies that encourage the provision of preventive healthcare services and promotion of a healthy lifestyle. However, the following improvements may be considered in the future. Multi-model ensembles can be used to improve the performance of GCM simulations instead of a single GCM to reduce the uncertainties. The selection of suitable distribution continues to improve with time and hybrid distribution (i.e., a combination of two or more distributions) may be associated with BCM-QM_skewed. Furthermore, other causes/factors of CVD and various categories of patients, such as gender, age, and race, could be potentially examined to achieve more meaningful outcomes to assist the respective local authorities, leaders, and policy makers in making action plans and preparing preventive measures. Integrating various climate covariates indifferent to spatio-temporal scale in the BCM-QM_skewed model may also be considered for future studies.

The daily mean temperatures of 58 stations in peninsular Malaysia were obtained from Ir. Dr Wong Chee Loong. We would also like to thank Editage (www.editage.com) for English language editing. This research was funded by the Malaysian Ministry of Higher Education Fundamental Research Grant Scheme (FRGS/1/2018/STG06/UPM/02/5) grant number 5540130.

All relevant data are available from http://ghdx.healthdata.org/gbd-results-tool.

The authors declare there is no conflict.