Analysis of Crash Frequency and Crash Severity in Thailand: Hierarchical Structure Models Approach

Currently, research on the development of crash models in terms of crash frequency on road segments and crash severity applies the principles of spatial analysis and heterogeneity due to the methods’ suitability compared with traditional models. This study focuses on crash severity and frequency in Thailand. Moreover, this study aims to understand crash frequency and fatality. The result of the intra-class correlation coefficient found that the spatial approach should analyze the data. The crash frequency model’s best fit is a spatial zero-inflated negative binomial model (SZINB). The results of the random parameters of SZINB are insignificant, except for the intercept. The crash frequency model’s significant variables include the length of the segment and average annual traffic volume for the fixed parameters. Conversely, the study finds that the best fit model of crash severity is a logistic regression with spatial correlations. The variances of random effect are significant such as the intersection, sideswipe crash, and head-on crash. Meanwhile, the fixed-effect variables significant to fatality risk include motorcycles, gender, non-use of safety equipment, and nighttime collision. The paper proposes a policy applicable to agencies responsible for driver training, law enforcement, and those involved in crash-reduction campaigns.


Background
The Thailand Road Safety Center [1] stated that the reduction of road crashes is a crucial component of the country's sustainable development. Therefore, the Thai government intends to reduce annual crash fatality rates to less than 10 people for every 100,000 people by 2022. However, according to the WHO estimates, the annual death rate in Thailand remains at 32.7 people for every 100,000 people. Thailand is ranked 9th in the world in terms of crash fatality rates [2].
Given the statistical trend of road crashes under the responsibility of Thailand's Ministry of Transport, previous studies found that the highest number of crashes and deaths occurred on roads under areas supervised by the Department of Highways (DOH). Specifically, 17,045 crashes, 2651 deaths, and 15,835 injuries occurred in 2018 (2562 B.E.) [3]. Highways are main roads allocated for the transportation of people or goods. It leads to relatively high averages in daily traffic volume. Thus, a decrease in the said rates should be initiated on highways.
The significance of the frequency of road traffic crashes requires an analysis of several characteristics based on the unique features of each area. These include road designs [4], applied and proven to produce more accurate predictions than the traditional models. In another study focusing on the exposure of the parameter-influencing factors through meta regressions, Ziakopoulos and Yannis [29] conducted a spatial analysis and found that traffic volume positively contributed to the number of crashes from the first set of variables. For roadway length dimension, an increase in roadway length tended to increase the number of fatal crashes. Moreover, the variable of distance traveled would cause the crash count to increase only at the country level, as opposed to the traffic analysis zone.

Highway Organization of Thailand
The present study considers the organization of the Department of Highways (DOH) in Thailand. Here, the hierarchy of road responsibility is an interesting issue. Under the DOH, the "Office of Highways" is composed of 18 offices across the country (Figure 1), followed by "highway districts" composed of 106 districts. The different characteristics of the two sub-organizations, such as the congestion index, scope of responsibility, annual traffic volume, responsibility for road distance, and annual budget, may affect crashes and their severity levels. This notion seeks to establish a spatial model by determining the two levels of factors that influence the dependent variables: crash case level (i.e., drivers, environment, or roads) and spatial level. They pertain to the operational area (i.e., operational area under responsibility, traffic volume, and congestion index).

Methods
Studies on road traffic crash analysis at the local and global scales are divided into two groups of factors based on outcome: (1) crash occurrences and (2)   Ziakopoulos et al. [15] performed a review where they found that there was no spatial approach on a middle-income and developing country. They also suggested that the heterogeneity of crashes may be different. They are expected to produce varying results due to varying driving culture or other unobserved factors. As such, this study could fulfill the scientific knowledge of road traffic crashes in developing countries. Furthermore, if the model results are attained, they can be used to determine spatial policies. The results can likewise promote campaigns for reducing the number of crashes and death rates in Thailand, such as in other middle-income and developing countries.
The rest of the paper is structured as follows: Section 2 describes the methods used for data analysis; Section 3 presents the results of the analysis and a discussion in two parts: the crash frequency model and crash severity model; and Section 4 provides the summary and limitations of the study.

Methods
Studies on road traffic crash analysis at the local and global scales are divided into two groups of factors based on outcome: (1) crash occurrences and (2) crash fatalities. Figure 2 presents the procedure followed in this study.

Crash Frequency Models
Count data consist of non-negative integer values. They are normally used for transportation models, such as number of routes, drivers' change per week, and number of crashes on roads per year. For the probability equation of the Poisson model, α takes a value of zero. This implies that the choice between the two models (i.e., Poisson and NB) is dependent on the value of α, which is relatively frequently over-dispersed model. Washington, et al. [30] provided a test for overdispersion. If the overdispersion parameter is significant, the NB model would be more appropriate than the Poisson model.
The probability of the number of collisions occurring on a road segment is considered an NB distribution. Predicting the annual number of crashes per year should consider that no crashes may take place on certain roads. Thus, this variable can be categorized as normal-count and zero-count (Figure 3).

Crash Frequency Models
Count data consist of non-negative integer values. They are normally used for transportation models, such as number of routes, drivers' change per week, and number of crashes on roads per year. For the probability equation of the Poisson model, α takes a value of zero. This implies that the choice between the two models (i.e., Poisson and NB) is dependent on the value of α, which is relatively frequently over-dispersed model. Washington, et al. [30] provided a test for overdispersion. If the overdispersion parameter is significant, the NB model would be more appropriate than the Poisson model.
The probability of the number of collisions occurring on a road segment is considered an NB distribution. Predicting the annual number of crashes per year should consider Sustainability 2021, 13, 10086 5 of 19 that no crashes may take place on certain roads. Thus, this variable can be categorized as normal-count and zero-count ( Figure 3).
Washington, et al. [30] provided a test for overdispersion. If the ove is significant, the NB model would be more appropriate than the P The probability of the number of collisions occurring on a road an NB distribution. Predicting the annual number of crashes per ye no crashes may take place on certain roads. Thus, this variable can mal-count and zero-count (Figure 3).  The general model may not be comprehensive for separating parts. Therefore, the most suitable model for the dual state is the zer called the zero-inflated negative binomial (ZINB) when established Vuong statistics were used to compare between the ZINB and the N the density function of the ZINB and f 1(.) as the density of function decision guidance suggested that if the Vuong statistic is 1.96, and NB overdispersion parameter is α > |1.96|, the zero-inflated models because it is evident that some simultaneous zero-altered process III :Count and discrete dependent variable models "of Washington, cusses this analysis. The models have crash frequency as the respo highway characteristics or traffic volume as explanatory variables.
In certain cases, the reasons for the prediction may include The general model may not be comprehensive for separating the analysis into two parts. Therefore, the most suitable model for the dual state is the zero-inflated model. It is called the zero-inflated negative binomial (ZINB) when established on the NB model [31]. Vuong statistics were used to compare between the ZINB and the NB model with f 1 (.) as the density function of the ZINB and f 1 (.) as the density of function of the NB model. The decision guidance suggested that if the Vuong statistic is 1.96, and the t-statistics of the NB overdispersion parameter is α > |1.96|, the zero-inflated models are plausible. This is because it is evident that some simultaneous zero-altered process is at work [32]. "Part III: Count and discrete dependent variable models" of Washington, et al. [30] further discusses this analysis. The models have crash frequency as the response variable and the highway characteristics or traffic volume as explanatory variables.
In certain cases, the reasons for the prediction may include a mutual correlation among observable data. This relationship is due to spatial data, such as data on crashes within the same area. Thus, the effect should be determined and not just observed. As such, the present study classifies roads into areas according to the designation of the subdepartments of the Department of Highway (DOH). In light of this relationship, model application should generate random and fixed effects, especially when an unpredictable factor is taken as an indicator variable.
The equation for the random effects is modified using the equation for log-linear function: In this equation, λ ij refers to the expected number of crashes for road segment i in group j (i.e., the regulatory space of the sub-departments of the DOH, where unobserved heterogeneity is expected). On the other hand, X ijm pertains to the vector of the explanatory variable mth, while β mj stands for the parameter prediction vector. ε i is the random effect representing within-area variance. µ 0j and µ mj denote the unique effect associated with sub-areas of the DOH. ε i , µ 0j , and µ mj are usually assumed at normal distribution with the mean zero. Spatial correlation is calculated based on the proportion of spatial variation from the total variation [33][34][35][36]. The equation is as follows: where σ 2 η represents the variance obtained from the prediction of i or variance within the same sub-DOH area and σ 2 e refers to the variance obtained from the estimation of the fixed effects or the variance between the sub-areas of the DOH.
Program R was used for analysis. The dplyr package (a grammar of data manipulation) was used for data input [37]. The ggplot2 package was used to create the graphic of the data distribution. The traditional crash frequency models include Poisson that was analyzed by glm (stats) (fitting generalized linear models) and NB that was analyzed by glm.nb (MASS) (fit a NB generalized linear model). On the other hand, the ZINB was analyzed by zeroinfl (pscl) (fit zero-inflated regression models for count data via maximum likelihood) [38]. The Vuong (pscl) was used to check the plausibility of the zero state. With regard to the final model, the spatial zero-inflated negative binomial models (SZINB) were analyzed by glmmTMB (fit a generalized linear mixed model using Template Model Builder) [39].

Crash Severity Modeling
The study determined the factors influencing crash severity at the second level in terms of spatial characteristics. Such was carried out under the supervision of each district divided by the boundary of the province or district for large provinces. Logistic models applying spatial concepts are first used to determine driver factors (driver i) in the district group j. Each crash consisted of the driver (ith driver) in the jth road segment. The fatality of common drivers was calculated by using the log odds in general areas. In other words, various characteristics of the areas have indifferently affected the likelihood of fatality, which determines the two severity levels including the fatal and non−fatal crashes. Here, Y ij p ij ∼ Bernouilli p ij , p ij = Pr Y ij = 1 represents the probability of driver i from road j that has fatal crashes, coming from the relationship of estimated parameters of explanatory variables such as driver factors (e.g., gender, drunkenness, vehicle type), road factors (e.g., median type, slope), and environmental factors (e.g., visibility). Such a relationship can be derived as follows: Here, η ij represents the log odds of driver i on road j that has fatal crashes. β 0j is the constant or the average value of log odds 1 that has fatal crashes on road j only. X kij is the kth explanatory variable at the individual level for predicting the odds of fatality j. On the other hand, β kj are parameter values that indicate the relationship (slope) between the driver level variables and the log odds that possibly causes the fatalities from the collisions. Equation (6) assumes that crashes on each road result in different degrees of severity. Thus, the equations are different from those of common logistic models. Each road has a constant value β 0j . However, the slopes β kj are different. Considering a spatial model as an explanation for the variation in the regression coefficient and variables at the upper level based on road characteristics, these qualities are assumed to be constant. Slope values, such as the median types and intersections, exert their effects at the driver level (Equation (6)). These can vary according to the explanatory variables of each road segment Z mj . The equations are formulated as follows: γ 00 indicates the log odds of fatal crashes on the road. On the other hand, Z mj represents the mth explanatory variable for road j, and M denotes the number of road level variables. γ 0m indicate the slope or relationship with Z mj corresponding to β 0j and β kj , respectively. γ 00 is the intercept for or β 0j and γ k0 is the intercept for β kj . µ 0j and µ kj stand for the prediction errors in information at the road level and demonstrate the unique effect of road j. Kim et al. [19], Huang et al. [20], Xiao et al. [40], Se et al. [41], Champahom et al. [42] provide for extended details on the calculation presented above.
Before the spatial analysis, the test results should be obtained from the models without parameter estimation (unconditional model). This may be done by considering the proportion of variance of dependent variables (outcome) within (crash on the same road) and between (each road) groups. The latter are called intra-class correlation coefficients (ICC). They were calculated by assuming a logistic distribution for the dependent variables. For errors at the driver level with a variance value of π 2 /3 [43], the ICC values can be calculated as follows: where σ 2 µ 0 is the variance of the dependent variables between roads (level 2). The ICC values calculated for multilevel modeling should not be zero. In the case where the ICC value is zero, then no variation exists between the data of crashes on each road. According to previous research, the ICC values for second-level variables, such as intersection, should be higher than 0.09 or 9% [19].
The following step assesses the model suitability by comparing the intercept-only model with all parameters set to zero. Moreover, the convergence model is compared with parameter vector β. The values used for comparison were the log-likelihood of both models. This is also called pseudo-R 2 or McFadden's ρ 2 : where LL(β) is the log-likelihood value of the model with parameter estimation. On the other hand, LL(0) represents that without or with parameter estimation but only with a constant value. If ρ 2 has a value close to 1, then the model can predict values close to the actual data [30]. R program was used for analysis. Moreover, glm (stats) with the family is binomial, and it was used to analyze the traditional model. The glmer (lme4) (fitting generalized linear mixed-effect models) was used for coefficients using the spatial analysis technique [44].

Data Collection
The highway crash reports on cases of road crashes from 2011 to 2017 that took place on Thai highways. During collisions, a highway officer investigates and records each case. The record contains information about the scene of the crash (i.e., specific name of the road and location in kilometers), date, physical characteristics of the road (i.e., median type and intersection), environmental factors (i.e., weather, time, and lighting), vehicle type, and crash type. The data collection was supported by the highway Bureau of Highway Safety, Department. Before developing the model, the correlations between explanatory variables were checked. It was found that there is no pair with correlation greater than 0.7, which was an appropriate indicator when collinearity begins to severely distort model estimation and subsequent prediction [45]. Table 1 provides a summary of the crash frequency models concerning road segment data. These were divided according to the physical characteristics of roads. The study found that the mean value of crash occurrence is 6.39 times, the segment length is 3.08 km, the lane width is 3.47 m, and the road with median is 33%. In terms of traffic volume, findings show that an average of 14,465 vehicles per day are on the streets, of which 16.33% are trucks. The maximum annual average daily traffic (AADT) is 339,248 vehicles per day. This value refers to the district of Thanya Buri and Klongluang of Prathum Tani province. This route is road no. 1 (Phahonyothin Road), which has 10-14 lanes with a barrier in the middle. This area is an economic and complex area. This means that there are many trip production and trip attractions including industry zones, complex residences, universities, and shopping malls. In addition, the road connects the biggest regions of Thailand including the "North" and the "Northeast." The official report of the AADT on this road could be assessed by the Thailand Department of Highways [46]. In Table 2, the overdispersion parameter is significant (0.143, SD = 0.002, t-statistic > 1.96). Therefore, it can be concluded that the NB fits that Poisson model. The result of Vuong statistic is 18.916 (t-statistic > 1.96). It could be stated that the ZINB is proved to be more appropriate. In summary, the study found that the Akaike information criterion (AIC) values continuously decrease until the SZINB reaches the smallest value. Here, a small AIC value indicated a model with better fit [12,47]. However, this study was analyzed based on the spatial analysis model, which used correlation spatial data. Vaida and Blanchard [48] suggested that the conditional Akaike information criterion (cAIC) was a more appropriate model than the ordinary AIC. The penalty term in cAIC is related to the effective degrees of freedomρ [49]. The cAIC allows for the comparison of models with different random effect structures. It likewise allows for the comparison of mixed-effect models with clusterspecific models where the parameters are fixed. The result of the cAIC shows that the SZINB is the smallest. This finding can be interpreted in such a way that the SZINB is the most-fitted model [50].

Crash Frequency Model
, p is the number of column or the number of fixed effects and N = number of road segments,ρ is the effective degree of freedom. Pseudo-R 2 is calculated as Pseudo-R 2 = 1 − ((LL(m0)/LL(m1)).
The pseudo-R 2 of SZINB is relatively low because there are many variables (e.g., human factors, weather, lighting, traffic signs, and land use factors) that are not measured [51]. However, the pseudo-R 2 of SZINB was still within the acceptable criteria of Mohammadi et al. [52], Abdul Manan et al. [53]. These authors all accepted the goodness of fit at 0.06.
The model results of SZINB (Table 3) were presented in three parts: 1.
The random effect conditional model refers to the prediction of the relationship between the independent and dependent variables. These variables were grouped according to the sub-areas of the Department of Highway (DOH).

2.
The normal state denotes the factors influencing collision number as non-zero state.

3.
Lastly, zero state pertains to the factors that do not lead to crashes. The zero sate is modeled to find the significant variables that could be specifically led to be the zero-crash segment. These significant variables can be developed to be effective road safety policy.
The inverted signs between states mean that the factors positively or negatively affected the crash frequency. These factors include road length and LN_AADT. It is positive in a normal state and negative in a zero state. Thus, these variables were clearly interpreted to have increased crash frequency. This was expounded in the study of Dong et al. [54].
According to the intercept of random effect, it is statistically significant and different from zero [55]. This could mean that the crash frequencies vary across office highway districts. The variances of random effects for variables of road segment levels were found to be insignificant. In other words, the parameter estimation values were lower than those of the standard errors. Thus, the results indicate that such variables do not differently affect crash frequency in the sub-areas of the DOH. This finding is consistent with the study done by Aguero-Valverde [28]. He found that the random effects of spatial error were not significant. The fact that the variance of random effect was not significant could be explained by the highway district office having no different control power on distance. As such, the traffic volume and the proportion of trucks were not significantly different in each area under each highway district office's responsibility. Therefore, a management policy that aims to reduce the number of crashes can be uniformly applied throughout the country. Traffic volume and the presence of road medians were the top factors that lead to crash occurrences. The reason behind this notion is that these two variables are clearly related. Roads with high volumes of traffic and that are built with medians subsequently lead to increased chances of crashes. Aguero-Valverde [56], Liu et al. [57] provided relevant study results. The reason why the existence of traffic islands caused more crashes could be explained by the fact that the Department of Thai Highways designed roads according to their traffic volume [58]. Roads with traffic islands are usually roads with high traffic volume [57]. Road shoulder and the traffic island width are other factors that contribute to car crashes. These results were inconsistent with those of Joon-Ki et al. [8], which suggested that roads with narrow shoulders led to more frequency of trailing. However, for developing countries, Mahmud et al. [31] found that shoulders did not contribute to more frequent crashes. As for the relationship of these two factors with traffic volume, it was found that the two factors positively correlated with traffic volume (Figure 4). Nevertheless, these two factors should be further investigated to discover why increasing road width of traffic islands and shoulders increased the crash frequency in Thailand.
For zero states, the majority of negative values indicated that significant factors typically lead to an increase in crash frequency. For example, the increasing distance between road segments can lead to the non-zero state, followed by the existence of medians and increase in traffic volume [56]. In view of the intercept of zero states, it was clearly found that it had a significant positive effect on a non-crashes. This implied that the mean frequency of crashes on each road segment was zero. Meanwhile, most of these factors led to a negative effect in zero states, which means an increase in crash frequency.  For zero states, the majority of negative values indicated that significant factors typically lead to an increase in crash frequency. For example, the increasing distance between road segments can lead to the non-zero state, followed by the existence of medians and increase in traffic volume [56]. In view of the intercept of zero states, it was clearly found that it had a significant positive effect on a non-crashes. This implied that the mean frequency of crashes on each road segment was zero. Meanwhile, most of these factors led to a negative effect in zero states, which means an increase in crash frequency. Table 4 displays the data on crash severity, which indicates the number of road crash victims. After screening, data on 37,685 victims (2014-2018) were collected. Results showed that in terms of driver characteristics, age was mostly a factor for middle-aged people. Moreover, it was shown that most drivers do not use seat belts and helmets. Drunk driving is found to be very rare, which accounts for only 2% of the total number of crashes. For road factors, the study found that the overall picture points to roads with medians of 33%. Median types consist of raised median (26%), depressed median (23%), and colored median (5%). Environmental factors imply that most crashes occur during daytime (60%). On the other hand, nighttime crashes and crashes in poorly illuminated areas reached 10%. In terms of crash-type factor, the study further found that single crashes obtained the highest proportion at 49%. This was followed by rear-end crashes at 25%, side rear-end crashes at 13%, crashes against people traveling, and head-on collisions.    Table 4 displays the data on crash severity, which indicates the number of road crash victims. After screening, data on 37,685 victims (2014-2018) were collected. Results showed that in terms of driver characteristics, age was mostly a factor for middle-aged people. Moreover, it was shown that most drivers do not use seat belts and helmets. Drunk driving is found to be very rare, which accounts for only 2% of the total number of crashes. For road factors, the study found that the overall picture points to roads with medians of 33%. Median types consist of raised median (26%), depressed median (23%), and colored median (5%). Environmental factors imply that most crashes occur during daytime (60%). On the other hand, nighttime crashes and crashes in poorly illuminated areas reached 10%. In terms of crash-type factor, the study further found that single crashes obtained the highest proportion at 49%. This was followed by rear-end crashes at 25%, side rear-end crashes at 13%, crashes against people traveling, and head-on collisions. Table 5 provides the parameter estimation results given that frequency data are considered with empirical data. The study found that the mixed-effect model produced more accurate prediction results than the fixed-effect model. This is evidenced by the AIC and pseudo-R 2 values. This result is in line with many past studies [19,24,59]. An ICC value of 0.24 can be interpreted as the variability of severity in highway districts at approximately 24%. Thus, the dataset of the current study is suitable for spatial model analysis [19,20]. Table 6 illustrates the comparison of the parameter estimation between traditional logistic regression and the mixed-effect model. The overall numbers of the factors influencing fatal crashes are found along the same direction. In terms of fixed parameters, the mixed-effect model has a smaller number of significant variables. This finding is relevant to that of Yu et al. [47], Wali et al. [60]. Few studies found that the mixed effect model has significant variables greater than or equal to a fixed-effect model [19]. However, there is no study necessarily comparing the number of significant variables between mixed and fixed-effect models. Contradictorily, it was found that the estimation is insignificant. This result is relevant to the study by Ahmed et al. [26], which found that the parameter estimation of the road grade between the Poisson and spatial model (−0.530 and 0.273, respectively) was insignificant.

Crash Severity Model
In terms of vehicle size, the findings show that small cars, such as motor vehicles, incurred the largest fatality risk. This point is understandable due to the physical characteristics of small cars (i.e., mostly motorbikes), which are relatively less secure than mediumand large-sized cars. Thus, high fatality risks are the result [21,22,61]. Similar to many results, females were more likely to die than males. This notion is possibly caused by the females' decision-making phase, which takes more time than that of males [5]. The outstanding result, which indicates that seat belt use can lead to decreased instances of death, is similar to that of Delen et al. [62]. In their study, they found that using seat belts was the first factor in assessing crash severity. Moreover, drunk drivers were more likely to die than those who did not drink [7]. With regard to the driver's age, the present study found that the elderly were more likely to die. This tendency is consistent with the work of Abdel-Aty [5], who posited that death may be common among the elderly due to their physical condition. They are more prone to injuries than younger drivers [63]. In terms of environmental factors, nighttime collisions were most likely to cause fatality, whereas collisions at night and in poorly illuminated areas increased the likelihood of death [18,64,65].
Road factors pointed to the main traffic lane as a source of low fatality risk [66]. A subsequent variable is the areas where vehicles decelerated before entering. The model results revealed that crashes at median openings and intersections also led to less cases of death. However, rear-end crashes frequently occurred in such areas. An explanation for this finding is that auxiliary lanes are located in areas with median openings. In addition, clear warning signs are likely posted before junctions. Therefore, every car may possibly recognize these areas and slows down. This result is consistent with that of Chen et al. [21].   Lastly, in terms of crash types, the parameter estimation still goes on the same side in both the traditional and random effect model. According to the random effect model, there were two types of collisions in general. These are the pedestrian and head-on crashes. These two types of crashes had non-significant effects on the fixed coefficient. In the random effect, the pedestrian crash had not been found to be significant, while a head-on crash was significant, which implied that when a crash pattern was determined to vary along the road segment, there would be a significant effect on crash fatality. This is relevant to the study by Zeng et al. [63] who found that head-on crashes are significant on the variance of a random term. In Thailand, head-on crashes are slightly more dangerous. This is because in some districts (e.g., industrial zones, university zones), small vehicles are often reversing. For sideswipe crashes, it was found that they tended to have a significant effect on non-fatal collisions. Such finding was inconsistent with that of Chen et al. [21]. Nevertheless, when considering the constant values from the model of Xiao et al. [40], it was found that sideswipe collisions (same direction) negatively affected fatality crashes. Meanwhile, rear-end crashes were often the collisions that resulted in no fatalities. These results were consistent with the findings of Kim et al. [19]. They found that when the intercept was compared with crash fatality among crash types, it yielded the most negative values of rear-end crashes. In other words, the average rear-end crashes would cause the least fatality rates. The reason for this is that rear-end crashes often occurred frequently, especially in an urban area where vehicle speed was not different. This resulted in only few instances of death. The studies concluded that rear-end collisions in urban areas were less violent than those in rural areas [42]. That is, a single crash led to fewer instances of death. Moreover, if there were no fixed objects, the likelihood of death was also lower [65].
The variance of random effect was found to be insignificant at the 0.05 level (i.e., intercept, meaning that the mean chance of fatal crashes was not different between road segments under the control of the highway district office and the other variables such as the divided median, flush median, raised median, depressed median, and median opening). This implied that the relationship between these factors and crash fatality did not vary with road segments [12]. However, many variables did (i.e., intersection and sideswipe crashes). This meant that each crash at an intersection in each road segment was different [55]. The causes of the difference may come from the control power of the highway district office over most of the urban roads. This directly affected the speed of vehicles and may lead to different levels of severity of intersection crashes. The variance of head-on collisions was also significant, indicating that the roads under the control of the highway district office caused different severity levels of head-on crashes. The results depended on different types of traffic islands. These results suggested that barriers should be built by highway district offices who have the control over several rural roadways.

Conclusions and Implementation
The study aimed to focus on the factors influencing the frequency and severity of road traffic crashes. The abovementioned factors were divided into two levels: collision case and highways. These factors were analyzed using the spatial model. The model results predicted crash risks for road segments under the responsibility of the highway districts.
For the crash frequency model, a Vuong z-statistic test found that a zero-inflated model is more appropriate than an NB model. According to the ICC, it was found that 9.7% of the crash frequency varies correspondingly with the responding area of the Department of Highway district (DOH). Moreover, the SZINB is the most suitable based on the conditional cAIC. According to the results of the SZINB, in terms of a fixed parameter that consists of a normal state and a zero state, these states are consistent. Meanwhile, if a factor is positive in the normal state, that factor would be negative in the zero state. For example, road length and traffic volume (i.e., for the random effect model, it was found that only the intercept is significant) could be interpreted as the crash frequency had unobserved heterogeneities. Thus, spatial analysis should be applied to analyze the crash frequency data. The crash severity model made an analysis based on the logistic regression models for understanding the relationship between the contributing factors and crash fatality (1 = fatal crash, 0 = non-fatal crash). A comparison between the fixed-effect and mixedeffect models was made. This was indicated by the AIC and pseudo-R 2 . It was found that the mixed-effect model is fit. Besides, based on the ICC of crash fatality, there was a 24% variance among road segments. Thus, the crash severity model should be accounted for by the unobserved heterogeneity model or a mixed effect based on the spatial model. The results of the fixed coefficient estimation between the two models are a consistent sign. The number of significant parameters in the mixed-effect model is fewer. The parameter estimation of random coefficients indicated many significant factors. These factors affected crash fatality varying among the road segments (e.g., intersection, sideswipe collision, and head-on collision). The overall fixed coefficient estimation makes sense. For example, a small vehicle crashing may more likely lead to a fatal crash than a large vehicle. Moreover, a crash at night, without lights, increased chances of a fatal crash. The study subsequently proposed policies that target the decrease in the crash frequency and fatal crashes.
The results confirmed that applying spatial correlation to the sub-areas of the DOH is more suitable than the traditional model in terms of crash frequency and severity in Thailand. The proposed guidelines are summarized as follows: The DOH can implement uniform policies throughout the country with special consideration for roads with median and high traffic volumes. In the case of such roads, a road safety audit may be implemented for a safety assessment. Authorities responsible for road supervision should regularly inspect the safety of such roads as an extra safety measure. For example, road safety audit (RSA) and identification of hazardous locations (IHL) should be thoroughly implemented to improve road safety in the case of limited funding. Road safety policy in Thailand with regards to RSA and IHL is not routine performed. When the DOH wants to perform RSA or IHL, they will hire a consultancy company. Thus, the RSA or IHL will proceed a few years at a time. Previous research about road safety in Thailand used to suggest a policy about RSA and IHL such as Jomnonkwao et al. [61], Champahom et al. [67], Se et al. [68]. However, they suggest only generally but did not specify. A new policy from the result of this study suggests the DOH should begin to do the RSA and IHL specifically on the road segment with the median with high AADT and more than four lanes road. Those road segments usually connect between large areas such as provinces or districts, since the fixed effect of AADT, median, and number of lanes are significant. Moreover, the random effect confirms significant variance among the department of highways offices.
Guidance policy reduces crash fatalities, which is selected based on significant parameters from the model. The guidances include four parts, as follows:

•
A survey for poorly lighted crash-prone areas should be first conducted. There is previous research that suggested this point such as Champahom et al. [42], Se et al. [68], Se et al. [69]. These works mentioned that the DOH should survey the dark road segments and consider installing lighting poles. According to the road safety policy in Thailand, there are regularly surveys of lighting conditions on the road segments. However, some processes are slow and limited in terms of budget. To conserve the budget for installing lighting-poles, this study recommends surveying a light condition of more than 3% at the slope segment first, since the factors significantly positively affect likelihood of fatality. • Followed by the posting of additional signs for anti-drunk driving campaigns, to reduce fatal crashes. Law enforcement and a campaign against drunk drivers is a regular proceeding, mentioned in previous works, namely, Jomnonkwao et al. [70], Se et al. [71]. To establish the guidelines for determining campaigns developed by Phillips et al. [72], healthy organizations such as the Thai Health Promotion Foundation should use personal communication or roadside media as part of their campaign delivery strategy. Lastly, checkpoints inspecting for alcohol use and seat belt use should likewise be reinforced. Recommendation of using safety equipment has been mentioned in terms of drivers' attitudes [73][74][75] and in terms of law enforcement by using CCTV for detection [76,77].

•
The Department of Land Transport can help promote such campaigns and topic training because head-to-head collisions and rear-end crashes are the most frequently occurring crash types that lead to high rates of fatality. Otherwise, short-term campaigns can be used as guidelines [72]. Therefore, it is necessary to consider areas prone to crashes, especially where driving on the wrong side of the road is a likely tendency. In Thailand, the campaign or the training topic of awareness of collision types are a low priority. The previous research recommended that the DOH should highlight some crash types such as rear-end crashes, single crashes and head-on crashes [42,67,68]. These studies suggest some policies, for example, the driver should be aware of the appropriate headway with the front vehicle. Regarding the results of this study, rear-end crashes and head-on collisions were found to have significant variance among the road segments. Thus, a new policy from this study recommended that the road safety government does not necessarily attempt implementation of safety guidelines throughout the country. However, they should focus on some areas, for example, roads with mixed vehicle types, such as industrial areas and communities without limited accessibility of trucks.

•
In addition, officials of the DOH should investigate the blind spots and risk spots, especially undivided roads and roads divided by flush medians. Currently, the safety policy of Thailand is trying to reduce the length of flush medians and to increase the length of barriers. It is relevant to the road safety studies in Thailand, which suggest reducing flush medians because it is more likely to sustain fatal injuries in a crash [42,69]. This study confirms significant spatial variance of median-openings and intersections. Therefore, the recommendation of this study is to begin considering installation of barriers to divide intersections and median opening points.

Limitations and Further Research
This study has its limitations. First, the independent variables have rarely been introduced into the crash frequency model due to the limited data in Thailand. Thus, the policies for reducing crash rates should be interpreted with caution. This is important, especially as there are other developing countries that will apply this idea in the future. Moreover, the pseudo-R 2 seems low because the relative factors were not included in the model. This study could not be fully accounted for due to the limitation in the dataset used. Other variables, such as the ratio of curves and the number of intersections, should also be considered [78].