Data source and study population
Our study drew upon a tailored database within the Health Care Bigdata provided by the National Health Insurance Service in South Korea. This data encompasses both outpatients and inpatient admissions diagnosed with heat-related illnesses, as per the International Classification of Diseases, 10th Edition (ICD-10: T67), spanning from 2011 to 2020. The cities covered in our study are Seoul, Busan, Daegu, Incheon, Gwangju, Daejeon, and Ulsan.
Table 1 presents the number of patients affected by heat-related illnesses (HRI). Across the decade, 296,829 individuals were diagnosed with HRI. Among them, 46,393 were aged 65 or older, while 250,436 were under 65. The dataset also included 13,735 disabled patients, segmented into 4,106 with severe disabilities and 9,629 with mild disabilities. The low-income and outdoor worker comprised 51,653 and 33,159 individuals, respectively.
Ethical approval
Our study protocol underwent thorough review and received approval from the Institutional Review Board (IRB) of Kongju National University (Confirmation No. 2022-46). Given the observational nature of this study and the utilization of anonymized statistical data, the IRB committee dispensed with the need for informed consent. All methodologies conformed to the guidelines stipulated by the Korean government concerning health and medical data usage.
Socio-economic classification
Participants were categorized based on socio-economic determinants into elderly, youth, outdoor workers, low-income individuals, and disabled groups. Those aged 65 and above were labeled as elderly, while those below this threshold were termed youth. The outdoor worker category encompasses professionals across sectors, including agriculture, hunting, forestry, fishing, mining, electricity, construction, and manufacturing. The low-income bracket was constituted of individuals in the lowest quintile for insurance premiums and medical aid beneficiaries. The type of disabilities defined officially by National Health Insurance Service includes people with physical disabilities, brain lesions, blind people, other disabilities (language, intellectual, autistic, mental, kidney, heart, respiratory tract, liver, face, ostomy, urinary tract, epilepsy) and national merit with related disabilities. Disability ratings follow the disability rating system and the criteria for determining the degree of disability in Korea. Korea is classified into severe (grades 1–3) and mild (grades 4–6) from 2019.
Weather data
The weather data were sourced from ASOS (The Automated System Observing System) for each city provided by the Korea Meteorological Administration. Our analyses incorporated daily maximum temperature (Tmax) and daily average humidity.
Empirical methodology
The statistical analysis was divided into two-stages. In the first stage, time-series regression was applied to each city in order to derive estimates of location-specific associations between heat-related illnesses and ambient temperature, reported as relative risk (RR). DLNM (Distributed Lag Nonlinear Model) was employed to analyze the relationship between heat-related illnesses and ambient temperature. DLNM has been widely used in epidemiological26,27,28,29. The key features of the DLNM is nonlinearity, time lagged effects, and control for confounders. The model allows capturing nonlinear relationships between the exposure and outcome variables. The model account for delayed effects of exposures on outcomes. The model can incorporate adjustments for potential confounding factors or time-varying covariates. The variable composition of the model for this study is as follows:
$$\textlny_t=\alpha +\beta \times Tmax_t-3+NS_1\left(sn_t\right)+NS_2\left(rhavg_t\right)+NS_3\left(doy_t\right)+\gamma \times weekday_t+\varepsilon _t,$$
where \(y_t\) denotes the number of patients with heat-related illnesses, \(rhavg_t\) is the daily average humidity, \(doy_t\) is the day in a year, and \(sn_t\) is the serial number. The degree of freedom of \(sn_t\) is set to 6 times the spanning time of the data, 10-year. \(weekday_t\) is the day of the week, Natural cubic spline (NS) shows a nonlinear relationship between dependent and independent variables. Other than the variables at time \(t\), a daily maximum temperature, \(Tmax_t-3\), is the 3-day lagged times series. The use of the time-lagged series is in line with Hess et al., Heo et al., and Royé et al., which provide evidence that the impact of heat waves strongly remains on the period26,27,28.
Once the DLNM is fitted, the estimated coefficients obtained from the model can be used to calculate the relative risk. The DLNM estimates a log-linear relationship between exposure and outcome, the relative risk can be derived by exponentiation the coefficient associated with the exposure variable, as follow.
β is the estimated coefficient from the DLNM representing the log-relative risk. The relative risk measures the ratio of patient occurrence over certain temperature while holding all other confounding factors.
To compare relative risks between groups or cities, we converted nonlinear RRs to the average RR above the 90th percentile daily maximum temperature. The daily maximum temperature at which the number of patients begins to increase, in other word, the threshold temperature, varies for each city depending on location-specific temperature range. We found that the threshold temperature was similar to approximately the 90th percentile daily maximum temperature for each city (You can see this in the results). Therefore, to express the nonlinear RR as a quantitative effect size to facilitate comparison between cities or groups, and to use it as an effect size for heat waves in the next stage, we averaged the RR above the 90th percentile daily maximum temperature.
In the second stage, meta-analysis was applied to synthesize and analyze information across multiple cities. In meta-analysis, we pooled the RRs from the first stage to calculate a more accurate or comprehensive relative risk. The weight assigned to each effect size (means the RR above the 90th percentile Tmax in this study) from location-specific time-series regression analysis reflects the contribution to the overall effect estimate. The weight for an effect size was determined by the inverse variance method in this study. In the meta-analysis, a fixed-effect model and a random-effect model are employed to estimate the average effects. The fixed-effect model is for data with characteristics of homogeneity on population while the random-effect model is for heterogeneity. Q statistics and Higgins’ I2 are applied for the homogeneity test for the population. The homogeneity tests were conducted using Q statistics and Higgins’s I2. Due to the limited number of cases in the meta-analysis, homogeneity was considered significant if the p-value in the Q-test was less than 0.1 and if I2 exceeded 50% (Table 2).
The first-stage regression was performed with the R software (version 4.1.2), using functions in the package dlnm (version 4.3.0). The package meta (version 6.5.0) of the R was used for the second-stage meta-analysis.
link