Increase in traffic injury crashes following the 2016 Kumamoto earthquake in Japan: A model comparison

Abstract Objective Following natural disasters, the number of motor vehicle crashes may increase as drivers are often forced to drive under stressful conditions. This study aims to analyze the changes in motor vehicle crashes that resulted in injury or death (injury crash) following the 2016 Kumamoto earthquake in Japan. An existing study reported that the increased crashes resulted in property damage following the earthquake; however, the effects on injury crashes remain unreported. Methods Interrupted time series analysis is employed to investigate the changes in injury crashes following the earthquake. The results are compared based on several time series models, including negative binomial and autoregressive integrated moving average models. Monthly injury-crash data from 2011 to 2020 in Kumamoto and Fukuoka city is used. Results The results reveal a 1,642-count or 20% increase (1.20-times increase, 95% confidence interval: 1.12, 1.27) in injury crashes due to the earthquake in Kumamoto city, where the earthquake damage was heavy. In contrast, statistically significant change is not detected in Fukuoka city, where the earthquake damage is negligible. Conclusion The results indicate that the earthquake has increased the motor-vehicle-crash risk and that traffic crash alerts are important following disasters.


Introduction
Driving under stressful conditions following the occurrence of disasters may cause motor vehicle crashes. The disaster victims may not sleep well because they may have to live in shelters or temporary housings. The victims who lost their houses will be worrying about their life recovery for the future. The workers in the affected areas require to work harder for the recovery and may be more tired than in a usual situation. The anxiety, stress, distraction, and sleep deprivation in these conditions may induce careless driving, resulting in motor vehicle crashes (Casey et al. 2019).
The literature reports have indicated increased motor vehicle crashes following the disasters; for example, increased crashes following natural earthquakes in Japan (Hino et al. 1997;Hagita and Yokozeki 2019) and the human-induced earthquakes in Oklahoma, USA (Casey et al. 2019). Possible driving performance reductions following the 2010 Christchurch earthquake in New Zealand were also reported (Helton and Head 2012). Moreover, terrorist attacks have also changed the drivers' mental condition, resulting in increased crashes (Gigerenzer 2004;Stecklov and Goldstein 2004;Su et al. 2009). Although these studies have consistently reported the possible increase in motor vehicle crashes owing to natural or human-made disasters, further studies are required to generalize the findings. Maruyama and Taguchi (2021) have analyzed the motor vehicle crashes following the 2016 Kumamoto earthquake in Japan, and the present study attempts to extend their study. Specifically, the authors have attempted to overcome the following three limitations of the existing research. First, in the study conducted by Maruyama and Taguchi (2021), only the motor vehicle crashes with property damage (property damage crashes) were analyzed; the changes in motor vehicle crashes resulting in injury or death (injury crashes) have remained unreported. In the Kumamoto prefectural level data, an evident increase in the property damage crashes, rather than in injury crashes, following the earthquake was found; therefore, they mainly focused on the property-damage-crash analysis. As shown in this study, a detailed analysis can reveal the increase in injury crashes as well as property damage crashes, following the earthquake.
Second, in the study conducted by Maruyama and Taguchi (2021), an interrupted time series (ITS) analysis based on Poisson and negative binomial regression models considering seasonality using Fourier terms, was reported; however, the results obtained based on the other time-series models were unreported. Comparing the results estimated by several models would be useful in validating the robustness of the results. Therefore, this study employs several time-series models and compares the results.
Third, in the study conducted by Maruyama and Taguchi (2021), the time periods considered for the data were limited. Only the property-damage-crash data from 2015 to 2018 were used because the data before 2015 were unavailable. This limitation hinders the smooth analysis of crash trends and the development of models prior to the earthquake. This study used the data from 2012 to 2020 and analyzed the trend changes in injury crashes following the earthquake.
To summarize, this study aims to analyze the injury crash changes following the 2016 Kumamoto earthquake using data from 2012 to 2020 and compare the results across several ITS models. Particularly, the authors have compared the results employing a Poisson and a negative binomial regression model, considering seasonality, using Fourier terms and monthly dummy variables, as well as autoregressive integrated moving average models.
ITS analysis is a valuable tool for evaluating the longitudinal effect of unexpected events and has been used in various research fields (Lopez Bernal et al. 2017;Calvert and Erickson 2020;Khan et al. 2022) (additional references can be found in the bibliography in Online Appendix A1). However, their application in the examination of the earthquake effects on traffic crashes is limited (Casey et al. 2019;Maruyama and Taguchi 2021). The aim of this study is not to judge the merits and demerits of ITS, but to apply and compare several ITS models.
The results indicate that the injury crashes in Kumamoto city have increased significantly during the 46 months following the earthquakes. These findings are useful for establishing countermeasures against traffic crashes, following the earthquake disasters. This paper is organized as follows. The next section describes the data used, and the methods section summarizes the several models used. The subsequent section shows the results following the discussion in the final section.

Data
Two significant earthquakes with magnitudes of 6.5 and 7.3 hit Kumamoto, Japan (JMA 2018) on April 14 and 16, 2016, respectively. The earthquakes led to more than 270 direct and indirect fatalities, 2,800 injuries, and 206,000 instances of residential damage, including more than 8,600 total collapses (FDMA 2019). The maximum number of evacuees in the Kumamoto prefecture was greater than 180,000 (CO 2019). Figure 1 shows the location of the study area. Kumamoto is located in the center of Kyusyu, the southern island of Japan.
This study analyzed and compared the injury crash data from Kumamoto and Fukuoka city. The epicenter of the 2016 Kumamoto earthquake is located in the town next to the Kumamoto city, and the town and city were damaged heavily due to the earthquake. Fukuoka city is located approximately 100 km away from the epicenter, where almost no damage was reported. Therefore, the data from Fukuoka can function as a control group in this study. These cities are the prefectural capitals of Kumamoto and Fukuoka prefectures. According to the 2015 census, Kumamoto and Fukuoka cities have a population of 741 thousand and 1.539 million, respectively. These two cities have been chosen for the present study because these are major cities in Kyusyu, and the injury-crash statistics are publicly available. No cities in Kyusyu were found to provide the injurycrash data for free for the period of investigation in this study.
The injury-crash data from Kumamoto and Fukuoka cities were retrieved from a statistical report (FC 2022;KC 2022). This study used the monthly injury-crash data from January 2012 to December 2020, and the comparisons of cities with/without earthquake damage were used to reveal the effects clearly.
Online Appendix A2 demonstrates the basic demographic characteristics of the two cities during 2012 and 2020. In 2016 (year of the earthquake), in the Kumamoto and Fukuoka cities, the population densities were 1894.87 and 4530.51 person/km 2 , average household sizes were 2.33 and 2.00, percentage of seniors were 24.6% and 21.0%, average incomes were 2,753-and 3,329-thousand-yen, number of passenger cars per population were 0.323 and 0.283, respectively. In these cities, a significant change was not detected in the population, number of vehicles, or other demographic data after the earthquake. Therefore, although some demographic characteristics are different, the two cities appear to be suitable for the ITS study.

Methods
In the present study, the ITS models present the changes in the level and slope of traffic-crash data in the periods impacted by the earthquake (that is, impact period). The authors have used the following five ITS models to analyze the monthly injury-crash data: four time-series count models and an autoregressive model. A Poisson regression model is given as follows.
where y t is the number of monthly crashes at time t, and Poisson k t ð Þ is the discrete random variable following Poisson distribution with mean k t : crisisðtÞ is 1 if time t falls within the impact period, and 0, otherwise. b 0 , b 1 , b 2 , and b 3 are the parameters to be estimated. b 1 is the slope in the unimpacted period, b 2 represents the level change at the beginning of the impact period, and b 3 represents the slope change during the impact period. Figure 2 illustrates the changes in the level and slope of the crash data during the impact periods. The seasonality(t) represents the seasonality in time t.
The negative binomial regression model is given as follows: where the distribution is parameterized with its mean, k t , using Eq.
(2) and dispersion parameter, /: A limiting case of the negative binomial distribution with / ! 1 corresponds to the Poisson distribution. The mathematical expressions of the Poisson and negative binomial regression are described in detail in Online Appendix A3.
In the model with Fourier terms, the seasonality(t) is set as follows: where a j and b j are the parameters to be estimated. To prevent overfitting, the maximum J ¼ 2 is set in this study.
In the model with monthly dummy variables, seasonality(t) is set as follows: where month ðt, kÞ is the monthly dummy and is 1 if time (t) is month k, and 0, otherwise. m k is the parameter to be estimated. The data in January is set as the base by setting m 1 ¼ 0: Hereafter, the Poisson regression and negative binomial regression model with Fourier terms are referred to as Model A and B, respectively, and those with the monthly dummy as Model C and D, respectively.
Subsequently, the seasonal autoregressive integrated moving average model with exogenous factors (seasonal ARIMAX or SARIMAX) is considered (Hyndman and Athanasopoulos 2018). The model is an extension of the ARIMAX (p,d,q) model, which can be represented as follows: where e t is the white noise that follows a normal distribution with mean zero and variance r 2 (that is, e t $ N 0, r 2 ð Þ). ; i , h i , b 0 , b 2 , and b 3 are the parameters to be estimated. SARIMAX(p,d,q)(P,D,Q)[s] is formulated similarly; however, its formulation is not detailed here (Hyndman and Athanasopoulos 2018). The following setting is used in this study: s ¼ 12.
Hereafter, the SARIMAX model is referred to as Model E. The ITS analysis using SARIMAX follows the method given by Schaffer et al. (2021), which is described in detail in Online Appendix A4. In summary, the possible models are summarized in Table 1.
The selection of impact model is important factor in obtaining unbiased results (Lopez Bernal et al. 2018), and the following procedure is used to determine the impact periods and possible level and slope change, as established by Maruyama and Taguchi (2021). The impact period is determined such that it minimizes the Akaike information criteria (AIC). The procedure for handling the possible structural changes in the model should be focused upon. The level and slope change of the model during the impact period (Figure 2) implies the structural change of the model. The standard AIC for the model with structural should not be directly compared with the AIC without structural change. Kurozumi and Tuvaandorj (2011) proposed a modified AIC that can address this issue. They demonstrated that a value of six times the number of structural changes should be added to the modified criteria. This modified AIC is used to determine the impact period.
Notably, the (modified) AIC should be used for selecting the model with the same sample size. Thus, it cannot be used to compare the SARIMAX model, with differencing, with the models without differencing. For such comparison, root mean squared error (RMSE) is used (Hyndman and Athanasopoulos 2018).
In addition, the statistically insignificant estimate of b 3 indicates the no-slope change for the impact period. Furthermore, the statistically insignificant estimates of b 2 and b 3 indicate the no-level and no-slope changes, implying no earthquake impact on the crashes. The impact period and impact model are set through the systematic procedure.
The report on the autocorrelation of error terms in ITS is also important. The Ljung-Box test is used in this study to investigate the autocorrelation. If the null hypothesis of no autocorrelation of error terms is rejected for a model with an impact period, it is not considered to be the final model. The text statistic Q is as follows.
where r k is the autocorrelation for lag k, h is the maximum lag being considered, and T is the sample size. The statistic Q follows an v 2 distribution with ðh À pÞ degrees of freedom under the null hypothesis of no autocorrelation, where p is the number of model parameters. Hyndman and Athanasopoulos  The estimation procedure can be summarized as follows. First, all the five models (Model A-E) were estimated. The impact periods for the models A-D were determined by the modified AIC, and that for Model E was determined by the method demonstrated in Online Appendix A4 (Schaffer et al. 2021). The model with the lowest modified AIC while satisfying the residual error test (Ljung-Box test) was selected as the best model among the models A-D. The best model was compared with the Model E by RMSE and the final model was selected.

Results
The annual injury-crash changes in both Kumamoto ( Figure 3) and Fukuoka (Figure 4) cities reveal a gradual decreasing trend. Following the Kumamoto earthquake in April 2016, a small increase in crashes can be observed for Kumamoto city. The ITS method has been employed to examine whether the increase is statistically significant.

Kumamoto city
In Kumamoto city, the negative binomial model with monthly dummy (Model D) was determined to be the best model (Table 2). It produced the smallest modified AIC among the general linear models (Models A-D) and smaller RMSE than that of the SARIMAX model (Model E). The SARIMAX model rejected the null hypothesis of no autocorrelation of error terms by the Ljung-Box test and cannot be the final model. Model D reveals no statistical significance in the Ljung-Box test, and the (partial) autocorrelation function of the residuals reveals no spikes (See Figure A2 in Online Appendix A5). This supports the no autocorrelation among the residuals and modeling assumption. Model D estimated the impacted period as from April 2016 to February 2020. All the general linear models produced similar impacted periods. See Online Appendix A6 for the details of parameter estimates of the models. Figure 5 demonstrates the actual, modeled, and counterfactual crash trends for Kumamoto city based on Model D. It clearly demonstrates the level changes in crashes during the impact period. The results reveal the 1,642-count or 20% increase (1.20-times increase, 95% confidence interval (CI): 1.12, 1.27) in injury crashes. The increase was derived from the change in the level and slope.

Fukuoka city
The procedure, same as that for Kumamoto city, revealed that SARIMAX (Model E) is the best model in Fukuoka city. However, no impacted periods were detected in Model E. The general linear model (Models A-E) revealed crash decrease after the earthquake; however, these phenomena are unlikely to be caused by the earthquake. Therefore, the results indicate that the earthquake did not affect the injury crash in Fukuoka city (See Online Appendix A7 for the details).

Discussion
The results of Kumamoto city revealed a significant increase in injury crash; whereas, Fukuoka city, as a control group, revealed no such increase; this clearly indicates the increase in injury crash in the earthquake-affected area.
The result is consistent with those reported in the previous studies. Hino et al. (1997) reported a 10.4-17.0% increase in injury crashes in the affected area, compared with the previous year in the 1995 Great Hanshin-Awaji earthquake in Japan. Hagita and Yokozeki (2019) reported a 11.9-15.3%   increase in injury crashes compared with those in the previous year in 1995 in Kobe city following the earthquake. They also examined the changes in injury crashes following seven natural earthquakes during the period of 2000-2008 in Japan. They reported a 7.3-16.3% increase in the injury crashes, compared with that in the previous year in urban areas.
Major parts of Kumamoto city are urban areas, and the present result of a 20% increase in Kumamoto city is consistent with their results. Note that their method is based on descriptive statistics; however, the results in the present study statistically support the increase in the injury crashes based on the ITS method and successfully estimate the impact period. Casey et al. (2019) reported a 4.6% increase in motor vehicle crashes in the highly exposed counties in Oklahoma, USA, in the month following the earthquakes, induced by wastewater injection (that is, human-induced earthquakes), based on the SARIMAX model. The extent of damage caused by the human-induced earthquakes is evidently less than that caused by the natural earthquakes in Japan analyzed in this study and the literature. Therefore, the small impact of the earthquake might explain the small increase and short impact period (that is, one month) reported by Casey et al. (2019). Maruyama and Taguchi (2021) analyzed the property damage crashes following the Kumamoto earthquake based on the ITS method. They classified Kumamoto city as the affected area, and the relative risk (RR) was 1.25 (95% CI: 1.14, 1.35). Kumamoto city consists of five wards and the RR range is 1.18-1.32. April 2016 to December 2018 was considered as the impact period. Therefore, the earthquake appears to impact property damage crashes more significantly than injury crashes.
The reason for the larger increase in property damage crashes compared with injury crashes following the earthquake remains unclear. Possible reasons may be related to the reason behind the increased crashes: anxiety, stress, distraction, and sleep deprivation due to earthquakes (Casey et al. 2019), which may induce careless driving, thereby resulting in small crashes. Drivers may be under a high-alert condition following the earthquake, which might have prevented heavy crashes. The difference between the changes in property damage and injury crashes following the earthquake is not extensively examined in the literature. Therefore, further studies are required, including the data that encompasses the cause of the crash, which is unavailable to the present study.
Two limitations of the present study are noted here. One limitation being that the monthly traffic flow data for the study periods were unavailable. The changes in traffic crashes may be related to the change in the traffic flow. The traffic flow data were monitored and made public only recently. Future studies may use these data to obtain more valid results. The second limitation is that monthly injury-crash data for the study periods for cities other than Kumamoto or Fukuoka were unavailable. Comparison with several earthquake damages will strengthen the results of this study in the future work.
Even with these limitations, the present study contributes significantly to the literature. This study is the first to demonstrate the statistically significant increase in the injury crashes following natural earthquakes based on the ITS method. The evidence of the increased traffic crashes suggests that crash alerts are important after earthquakes.