Spatiotemporal Investigation of Road Traffic Accident Severity for Transportation Safety Planning in Missouri: A GIS-Based Hotspot Evolution Framework

Due to multiple factors, the surge of road traffic accidents (RTAs) has become an urgent issue, causing economic losses and property damage. The World Health Organization (WHO) declared that 1.19 million people die each year [1]. Multiple reasons, such as driver behavior, weather conditions, and characteristic design of roads, contribute to both severity and occurrence of RTAs [2], [3], [4], [5]. With the increase in both population and transportation demand, accident risk continues to grow, especially with the robust correlation between speed and RTAs [6]. RTAs comprise two major specifications, both spatial dependence and temporal variability, causing their patterns to be dynamic instead of static [7], [8], [9], [10], [11]. Thus, understanding both spatial and temporal is regarded as a crucial aspect for safety mitigation. Multiple studies have utilized both statistical and Geographic Information Systems (GIS) techniques for assessment of RTA datasets [12], [13], [14]. However, some of these studies have focused on modelling like non-spatiotemporal or traditional. These techniques tend to oversimplify conclusions and do not reflect the dynamic emergence of RTA severity. However, a clear understanding of the spatiotemporal evolution of RTA severity at a macroscopic scale of analysis remains limited.

Notably, RTA hotspot identifications across space and time are still a complex challenge due to the interaction of miscellaneous factors like traffic volume, human behavior, and environmental conditions [15], [16], [17]. Besides, variations in RTAs reporting and severity classification, as well as limitations in the quality of the dataset, complicate the spatiotemporal analysis [18]. Though multiple researchers have examined spatial, temporal, or spatiotemporal methodologies, these studies are frequently inconsistent in handling spatial heterogeneity and insufficiently represent the dynamic interactions that impact severity patterns of RTAs. Moreover, multiple studies have focused on microscopic-scale analyses, thus restricting their relevance for statewide or regional safety planning [19], [20], [21], [22]. Consequently, there is a crucial need for a robust spatiotemporal analytical framework to identify and interpret the RTA severity hotspots evolution at a macroscopic level. Therefore, the aim of this study is to investigate the spatiotemporal patterns of RTA severity in Missouri, utilizing GIS-based techniques by addressing both spatial and temporal variation of RTA severity. The significant contributions of this research are:

1. Analysing RTA severity over a large geographic area, macroscopic perspective.

2. Examining the persistent and emerging hotspot areas of RTA severity with time-varying patterns.

3. Integrating spatial statistics with traffic safety management for severity analysis.

The main relevance of this study lies in its assertion to enhance the transportation safety and planning of decision-makers by utilising a GIS-based, spatiotemporal viewpoint. Thus, in order to inspect the statistical significance of RTA severity hotspots, the G$_i^*$ technique is applied to identify the spatiotemporal evolution within a time window. This methodology enables relevant agencies and planners to move beyond conventional descriptive statistics of RTA toward a deeper knowledge of inspecting the persistence or emergence of high-risk zones. More importantly, utilising both spatial and temporal analyses for hotspot inspection enhances the objectivity and reliability of safety assessments. Hence, this complements the literature in the discipline by demonstrating how statistical methods based on GIS can assist in filling the gap between the practical transportation and spatial analysis safety applications at the statewide scale. Multi-year RTA datasets were examined by this study for not only the hot spot clusters distribution identification, but also emergence of these clusters over a selected time scale. Therefore, the main contribution of this study lies in inspecting the high-risk zones to provide a clear indication for related authorities for road safety interventions.

In order to reduce the related fatalities and injuries of road users, an active practical safety planning is crucial to identify both persistent and emerging hotspots. Therefore, the outcomes of this methodology offer decision-makers analytical tools to focus on safety measures, allocate resources properly, and design robust strategies for zones with severe hotspots. The recommended analytical methodology can be easily adapted and applied to any geographical location, time scale, or mode of transportation.

Due to the surge in population rates, there is also an increase in the number of vehicles. With a growth rate of 2.17%, the commercial vehicle market is expected to reach around 40.29 million vehicles by 2030 [23]. As a result, it leads to an increase in car accidents. For instance, the U.S. Department of Transportation’s National Highway Traffic Safety Administration (NHTSA) revealed that around 40,901 fatalities were reported in 2023 [24]. The presence of spatial and temporal trends, correlations, and patterns, which are commonly referred to as hot spots and clusters, has been indicated by an extensive literature [12], [15], [25], [26]. Whilst the spatiotemporal analysis of accidents has been increasingly taken into consideration, the macroscopic scale of assessment has been addressed with less focus. Hence, integrating both spatial and temporal aspects for better capturing of the dynamic behavior of RTA severity patterns is crucial for safety enhancement. The RTA's severity fluctuates via spatial and temporal scales, especially when RTAs illustrate intricate spatiotemporal patterns. Therefore, it is important to understand these patterns for planning measurements and safety management. Besides, it is crucial to utilize techniques of spatiotemporal analytical models based on GIS to track the emergence of hot spots of RTA severity. For instance, Kernel Density Estimation (KDE) is often used for the spatial intensity of RTAs [27], [28], [29]. This approach is useful to produce smooth surfaces that indicate the density of incidents for high-risk areas visualization and preliminary hot spot identification. However, KDE is essentially descriptive, which lacks statistical significance assessment, hence restricting its capability for decision-making assistance. KDE is suitable when the aim is to visualise and identify the general spatial patterns without statistical inference.

In contrast, some geospatial autocorrelation methods like Moran’s I are commonly useful for the assessment of the degree of spatial dependence of the RTAs dataset [13], [26], [30], [31]. These techniques are useful for identifying the clustering patterns of RTAs by testing the spatial correlation of high and low values. Moran’s I can be useful in detecting geographical reliance, but it does not identify the accurate spatial locations of hot spots, which is regarded as unsuitable for practical applicability. Therefore, this technique is suitable when there is a need to assess the overall spatial dependence over the study area rather than identifying specific hot spot spatial locations.

Thus, to address these constraints, local statistics with spatial aspects, like G$_i^*$ method, have been widely utilized for hotspot detection [13], [28]. By utilizing neighbourhood local analysis, the identification of statistically significant clustering of both high (hot spots) and low (cold spots) RTA severity clustering can be implemented. Compared with global and density-based methods, G$_i^*$ provides more precise spatial insight and enhances safety planning [32], [33], [34]. Still, the majority of the literature is limited to a singular time period, which limits the understanding of the temporal evolution of hot spot patterns. In this case, whereas KDE offers only density visualization and Moran’s I evaluates global spatial dependency, G$_i^*$ is more preferable to detect the confined clusters, which makes it more appropriate for safety measures.

In addition, temporal dynamics play a crucial role in comprehending patterns of RTA severity [35], [36]. Multiple studies focused on integrating both spatial and temporal techniques for better capturing the evolution of the hot spot patterns over different temporal scales. In order to track whether the hot spot emerges, declines, or persists, the spatiotemporal approaches integrate both spatial and temporal dimensions. These methods introduce more comprehensive perceptions compared to spatial models. Also, these techniques are efficient for tracking hot spot patterns over various time intervals.

In general, though the performed studies have majorly contributed to comprehending the patterns of RTA, some gaps still remain. Various methodologies lack statistical reliability, fail to represent the dynamics of local hot spots, or improperly address spatial or temporal fluctuations at a macroscopic scale of analysis. Therefore, there is a crucial need to approach techniques that identify the statistically significant hot spot clusters while simultaneously inspecting their temporal emergence over broad geospatial regions.

In order to examine the RTA severity trend over the state of Missouri from 2020 to 2023, this study employs the G$_i^*$ statistic within a spatiotemporal framework. This methodology is suitable since it integrates the advantages of both local spatial and temporal assessment to address the recognized lacks in previous methodologies. By integrating both spatial statistical and temporal hot spots, this approach enables the identification of emerging, persistent, and diminishing patterns of hot spots, besides enhancing traffic safety planning.

The methodology of this study consists of three sections to examine the severity patterns of RTAs. The dataset resource and study area description are in the first part. The second part consists of the severity calculation methodology. Also, the final part of the methodology presents the spatiotemporal hot spot severity analysis to detect the statistically significant clusters of high and low severity crashes over a selected temporal scale.

Missouri is located in the Midwest of the USA, a diverse geography with an area of 68,727 m$^2$ and a population of 6,154,913 [37]. Figure 1 shows the state of Missouri with road networks. This research utilises the vehicular RTAs that have been collected by the state law enforcement and/or local agencies within the state for four successive years, from 2020 to 2023 [38]. A variable number of accidents were reported through the selected time periods, but the highest was 146,724 accidents during 2021, as illustrated in Figure 2. The dataset was first examined for accuracy and completeness. Also, reviewing whether there are any missing or unclear details, like spatial coordination or severity attributes. Removing any duplicates of entries was performed and extracted to avoid any bias during the analysis stage. Moreover, any detected inconsistencies were corrected as much as possible, and any significant errors were removed from the dataset. For the accuracy of spatial analysis, the spatial coordinates or RTA locations were checked by omitting any obvious positional errors. These preprocessing steps were crucial to ensure the reliability of the outcomes of the hot spot analysis.

Besides the spatial location (i.e., both longitude and latitude), multiple attributes are attached to each accident that can be utilized, such as the time, date, type, and reason of accident occurrence. Besides, the number of injuries and fatalities is also involved within the attributes of the dataset, as depicted in Table 1. Different temporal scales of analysis can be selected for the RTA analysis. For instance, periods of one month can be selected to capture weather patterns or seasonal impacts [39]. A smaller temporal scale of analysis could be used, such as hours, days, or weeks, depending on the dynamics of accidents. However, a one-year scale of spatial analysis is examined in this study for reducing random fluctuations of short-term scales, besides increasing the statistical stability of model estimation [40].

In order to classify the hot spot of accidents, the density of RTAs can be calculated by dividing the number of accidents of each assigned spatial area by the area itself in km$^2$. Also, the severity index is utilized in this study instead of the number of RTAs for multiple reasons, such as reflecting the consequences of accidents, supporting resource allocation priorities, and enabling advanced spatial and statistical analysis [41], [42]. This research utilized the severity index of the Belgian government with assigned factors of 1, 3, and 5 for minor injuries, severe injuries, and fatalities, respectively [43]. Eq. (1) illustrates the calculation of the accident severity index (ASI).

where, $M$ is the RTAs with minor injuries, $S$ is severe injuries, and $F$ is fatalities.

In order to determine the spatial autocorrelation of a geographical dataset, G$_i^*$ can be utilized for the identification of accident hotspots. This technique is widely used to identify the severity of high and low spatial clusters of traffic accidents, taking into consideration the local spatial autocorrelation. Rather than merely applying RTA counts, applying G$_i^*$ to severity-weighted indicators like injuries and fatalities is much better to detect spots with severe accidents and significant concentrations. Also, G$_i^*$ captures adjacent segments with similar attributes, such as traffic conditions or geometric design, then evaluates whether accidents with high severity values are surrounded by similar high severity values in adjacent areas. Hence, the results of $z$-scores and $p$-values illustrate a statistically robust recognition between the true severity hotspots and random spatial fluctuations. The spatial dependency of features as per G$_i^*$ is detailed in Eq. (2).

where, G$_i^*$ is the spatially correlated value of $i$ feature, $x_j$ is the amount of $x$ variable at spatial location $j$.

By utilizing G$_i^*$, the analytical accuracy of RTA analysis can be improved by using accident severity attributes to identify the spatial locations with excessively severe results instead of depending merely on the frequency of accidents. The determination of the spatial weight between features $i$ and $j$ can be calculated by the centred distance of the perceived spatial association. Thus, as per the distance value, the statistics of G$_i^*$ may provide different results due to the distance value. Moreover, $w_{i j}$ can be binary or non-binary values, and the sum of weights ($w_i$) is illustrated in Eq. (3).

Eq. (4) depicts that the sample mean and variance are employed for standardising G$_i^*$ values.

A fixed distance band was used for spatial relationships between accident spatial locations to ensure consistency in neighbouring identification features. The threshold distance of 0.8 km was selected based on the assessment of incremental autocorrelation. This selection identifies the optimum distance at which spatial clustering is the most apparent by identifying the optimum distance at which the $z$-score reaches its peak. This method ensures that the selected threshold distance detects the best significant spatial clustering pattern within the selected dataset and enhances the hot spot analysis. For reliability enhancement of threshold distance choosing, multiple potential lengths were assessed within a specified range before optimising the best value. The fluctuation of $z$-score values over increasing distance bands was systematically examined to identify the magnitude at which the spatial dependence is the most pronounced. Therefore, the value of 0.8 km correlated with the first stable optimum $z$-score values, representing the balance between avoiding inappropriate generalisation and capturing regional hot spot clustering. Thus, other values of distances below this threshold generated scattered clusters, while larger distances tended to produce a smooth spatial variation, reducing the sensitivity of identification of hot spot clusters. This validation stage confirms the suitability of the selected threshold for accurately depicting significant spatial correlations among RTAs spatial locations. A spatial weight was used as a binary (i.e., either 0 or 1). G$_i^*$ $w_{ij}$ is defined as 1 if the distance between features $j$ and $i$ is up to the specified threshold distance, while $w_{ij}$ = 0 otherwise. In this case, this method guarantees that merely nearby features within the designated spatial features impact the analysis.

For a better capturing of localized clustering patterns of accident severity, the utilized fixed distance band technique was selected instead of K-nearest neighbours for a better maintenance of a consistent spatial scale across the selected study area, besides a better capturing of localized clustering patterns of accident severity.

The implementation of G$_i^*$ within the GIS framework, it enables a clear vision of the spatial pattern severity of RTAs. It is worth mentioning that a local spatial autocorrelation analysis can reveal a set of distributed target attributes in the case of using the correct spatial autocorrelation between these specifications. Besides, it provides a robust foundation for prioritising high-risk locations and optimising planning for road safety enhancement to reduce fatalities and injuries.

During the planar spatial analysis, Euclidean or Manhattan distance can be employed. In this study, Euclidean distance was utilized for spatial proximity measurement among accident spatial locations. This method is useful for spatial analysis to provide a clear representation of spatial distance among features. The relationship between features can be examined by the investigation of spatial autocorrelation, which might provide information on how to comprehend the correlation among RTAs severity occurrences. The connection between events can be spatially examined in several perspectives, such as fixed distance, K-nearest neighbours, or the space-time window technique. The selected distance threshold in this research was established taking into consideration the spatial autocorrelation variability of inherent [25]. The segmentation of road networks has been utilized at this study as a spatial analytical unit instead of individual RTA locations. Each road segment is treated as a geospatial unit where the RTA data are allocated. One segment of road would comprehend multiple incidents that were consolidated by merging their severity attributes into a single metric utilizing ASI. This aggregate allows for representing the total impact of RTAs on each segment of the road network. The analysis of hot spots was conducted over these aggregated units of link-based data, enabling the detection of high and low severity levels of road segments. This approach illustrates a network-oriented perspective for clustering of hot spots, fitting closely with the linear characteristics of transportation infrastructure rather than point-based modeling.

It is crucial to know that extremely high G$_i^*$ $z$-score values can occur, especially in network-based RTA analyses, whenever using thorough segmentation of roads and vast datasets or accidents. The statistic of G$_i^*$ can be influenced by geospatial intensity besides the total number of investigated spatial units. Hence, multiple adjacent road segments with high aggregated values of severity may generate high positive values of $z$-score. Therefore, in this research, the determined $z$-scores illustrate a robust, statistically significant intensity of RTA severity clusters versus any methodological or computational errors.

By utilising a GIS-based framework, ArcGIS Pro (version 3.2), the G$_i^*$ statistics were applied for $z$-scores, $p$-values, and confidence level bins (G$_{i,\mathrm{Bin}}$) calculations for identifying the statistically significant clusters of both high and low RTA severity values. Moreover, a one-year interval can be selected as a temporal scale for this study. It is worth mentioning that the yearly temporal scale was selected to ensure both reliability and consistency of results. Hence, the annual aggregation of the RTA dataset was utilized to mitigate the impact of short-term fluctuation, inaccurate reporting, and irregular variation. These issues can be pronounced over selecting smaller temporal scales, like seasonal or monthly. Furthermore, the aim of this study is to identify long-term patterns in RTA severity rather than temporary fluctuations. Thus, the annual temporal scale is more suitable for persistent patterns detection over a selected time scale. In addition, the preliminary dataset inspection demonstrated the existence of probable sparsity and unreliability over certain time periods. This case can adversely impact the statistical efficacy of the analysis when smaller temporal scales are used.

It is worth mentioning that the selected methodology of this study focuses on Missouri as a case study, but it is applicable to any other geographical region. The application of both ASI and G$_i^*$ statistics can be used over various traffic and geographical conditions. However, the effectiveness of this technique may differ according to some patterns’ characteristics, like traffic volume, road network, and urban development. Hence, the application of this approach over other regions would require model parameter alteration and include some related factors to those areas. For repeatability guarantee, the analysis procedure was executed as a series of clearly specified logical steps. The dataset was processed by omitting repeated records, correcting errors, and validating geospatial coordinates. Later, the occurrence of RTAs was spatially assigned to their related road segments, and severity measurements were aggregated for each segment using the ASI. A fixed distance threshold technique was then utilized for assuring uniformity in neighbourhood definitions across all datasets. The G$_i^*$ statistic was then applied using uniform parameter values for all temporal scales, the annual scale of analysis. All outcomes, like $p$-values and $z$-scores, were generated using the same computational settings in ArcGIS Pro (version 3.2). Finally, the standardised methodology implies that each phase of the methodology might be duplicated using similar preprocessing rules, spatial linking techniques, statistical parameters, and aggregated logic to any similar set of data over any region.

Figure 3 depicts the monthly variation of RTAs in Missouri from 2020 to 2023. In general, the data reflects a seasonal trend, when the number of RTAs is recorded highest during October of each year. However, at the beginning of 2020, the number of RTAs was relatively high, but during April 2020, there was a sharp recorded drop. However, there was a gradual recovery over the successive months. A more stable pattern was recorded over both 2021 and 2022, when the RTA numbers fluctuated between 11,000 and 13,000. Also, the year 2021 records the highest number of RTAs during October, with around 14,200, due to potential impacts of traffic-related variations. The year 2023 follows the same pattern as the previous years, but with a lower level than 2021. The variation among different years reflects the impact of some external factors such as traffic volume, weather conditions, and mobility. It is worth mentioning that these seasonal fluctuations have been addressed by multiple researchers, where the frequency of RTAs increases during fall months due to the weather change conditions, reduction of visibility, and increase in traffic demand due to the schools and work activities [2], [17]. Nevertheless, the recorded drop during April 2020 is not predictable as normal seasonal patterns, which can be defined by the mobility reduction during the restrictions of COVID-19.

Figure 4 illustrates the age impact on the fatality of road users in 2022 in Missouri. The vast majority, the adults, contributed with 59% of fatalities. This percentage indicates that this age group bears the highest weight. However, seniors contributed for 21% of fatalities, which is also regarded as a significant contribution. Both youth and teen groups reflect the lowest recorded percentage of fatality, with both recorded 16% and 4%, respectively.

These indicated patterns are consistent with the previous literature, when middle-aged adults represent the maximum percentage of fatalities due to the highest exposure to traffic conditions and commuting activities [2], [5]. Moreover, due to the reduced reaction time, besides the increase in physical fragility, these high percentages of senior fatalities have also been reported by multiple studies. On the contrary, the lower recorded percentages among younger road users might reflect less driving exposure and obeying driving regulations.

The contributing circumstances that would cause the fatality of RTAs can be due to multiple factors, with different percentages. For instance, the speed contributes to the highest share of fatality with 29%. This highlights that excessive or unsafe driving speed can be regarded as the main cause. At 26%, the aggressive behavior of drivers reflects a substantial impact on the fatality of RTA occurrences. Both alcohol and drugs rank third in causing fatalities of RTAs, with 20%, which reflects the issue of driving while intoxicated. Lack of compliance with traffic regulations resulted in around 16% of severity during 2022. Though they reflect the lowest proportion at 9%, both distraction and fatigue are still contributing factors. Generally, the percentages illustrate that factors related to the behavior of drivers are significantly impacted by RTA fatalities, as in Figure 5.

The outcomes of this study are in line with previous studies, which referred to the main role of speeding and aggressiveness of the drivers in causing fatalities [6], [41]. The findings across different studies emphasise the major impact of driver behavior on traffic safety. Furthermore, the accidents related to drugs and alcohol might indicate variances due to social behavior and enforcement intensity compared to other areas.

For the spatial distribution of RTA density over the counties of Missouri, Figure 6 depicts multiple ranks of density. For instance, areas encompassing main metropolitan areas like St. Louis City and St. Louis County, followed by Jackson County (Kansas City area), endure the highest rates of accidents, particularly. This is due to the fact that these areas are recognised with high population, extensive road networks, and high traffic volumes. The second rank, such as St. Charles, Greene, Jasper, and Jefferson, illustrated the impacts of major cities besides regional transit routes. Also, most rural counties have the lowest accident rates. In general, the outcomes reflect a spatial clustering of RTAs rather than being uniformly distributed.

Due to some major factors like the increase in traffic congestion, interaction among road users, and the presence of traditional centres, metropolitan areas regularly have the highest concentration of hot spot clusters of RTAs. This pattern of spatial clustering is widely demonstrated in the literature of transportation safety [9], [13], [30]. In contrast, some literature illustrates that rural regions have a noticeable frequency of RTAs, when the severity may be elevated over these areas due to factors like the high driving speed of some road users, inadequate lighting, and lack of proper signage and road markings [18].

The geospatial RTA distribution is majorly impacted by the functional and hierarchy classification of road network. For instance, counties with substantial routes of interstate, freight corridors, economic activities, and metropolitan arterials illustrated a much higher density of RTAs than other rural areas. Any sites operating as administrative or economic centers typically generate higher trip generations, hence raising the occurrence of traffic incidents. Due to the presence of significant network connections, certain counties may indicate a high density of RTAs, where traffic flows interact. These outcomes significantly indicate that the distribution of RTAs is strongly correlated to the functional classification of roads in addition to the intensity of land-use, rather than merely spatial location.

The G$_i^*$ Statistics were utilized for spatial clustering evaluation by comparing the local sum of ASI within the neighbourhood road segments to the expected sum under spatial randomness. These levels of confidence (i.e., 90%, 95%, and 99%) illustrate the thresholds of statistical significance derived from $z$-scores and $p$-values. Higher confidence levels (i.e., 99%) depict statistically more robust clustering of hot spots, which resulted in larger $z$-score values and smaller $p$-values, reflecting a lower probability when the observed trend of spatial pattern is due to random chance. However, 90% reflects the lowest levels of confidence, which illustrates the weakest conditions but is still statistically significant. This classification enables recognition between the statistically significant clusters of severity with high or low values, in addition to the areas with identified patterns that do not significantly diverge from randomness. The results of the hot spot maps show a clear statistically significant spatial clustering of both high and low RTA severity across the study area based on the analysis of G$_i^*$. The red clusters depict high-severity hot spots, while the blue clusters show cold spots with low-severity levels. The intensity color variation indicates the magnitude of $z$-scores calculated from the statistics of G$_i^*$, where the higher values of $z$-score illustrate stronger spatial correlation of clustering for high or low values of severity. The segments of the road network were utilized as spatial units of analysis instead of relying on individual RTA points. Thus, each road segment was treated as a spatial unit for RTA points to be assigned. By using a spatial association technique, severity values of RTAs were spatially assigned to the nearest road segment (i.e., buffering within 50 m). The utilized aggregation illustrates the entire impact of RTAs on each related road segment. Afterwards, G$_i^*$ Statistics were performed for statistically significant identification of hot spot clusters of high and low severity road segments.

Instead of depending on distinct accident locations, road segments of road networks were utilized as geospatial units for analysis. As a result, each segment of road was identified as a spatial unit for the assignment of RTA incidents. Within a 50-meter buffering technique, the calculated severity level for each incident was geospatially assigned to the adjacent road segment utilizing a spatial approach of association. The utilized segmentation technique illustrates the thorough impact of accidents on each relevant road segment. Then, G$_i^*$ statistical methodology was performed to identify significant clustering of high and low severity levels of road sections. It is worth mentioning that the reported hot spot clusters reflect the statistically significant road segments instead of geospatial hot spot locations. The linked and linear characteristics of road networks enable various continuous segments of roads located on the same route that simultaneously demonstrate clusters with statistically significant differences. Thus, a single geospatial hot spot cluster may consist of multiple statistically significant portions, and this leads to highly reported hot spot counts. The segmentation resolution that has been selected determines the number of recognized hotspot cluster segments. Therefore, more precise segmentation produces a larger number of road units for analysis, which enhances the probability of additional statistically significant segment detections within the road network.

For the selected period of analysis (i.e., from 2020 to 2023), the outcomes reveal statistically significant evidence of spatial clustering of RTA severity across the selected temporal scale of analysis. For instance, the main metropolitan areas like St. Louis and Jackson counties depicted a highly persistent confidence hot spot (from 95 to 99%). This pattern of elevating RTA severity illustrates a fundamental high risk due to such factors as high traffic volumes, dense population, and complicated roadway designs. It is worth noting that the frequency of persisting hot spots in main metropolitan areas like St. Louis and Jackson counties is strongly correlated to the presence of main interstate corridors, complex interchanges, high-density population, and increased traffic demand. These areas have dense freeway systems with constant traffic interactions, lane changes, and congestion. Consequently, this increases the probability of spatially severe RTA clustering along surrounding road segments.

Besides, some minor metropolitan areas and some significant major routes also had persistent hot spots. On the other hand, due to factors such as low volumes of traffic, the results reflected some cold spots in some rural areas. The robust spatiotemporal analysis shows that the RTA severity clustering in overall Missouri is spatially consistent and statistically reliable, as depicted in Figure 7, Figure 8, Figure 9, and Figure 10. The persistence of hot spot clusters over multiple years reveals that the detected patterns are not random but essentially correlated within the system of transportation. This persistence is typically related to persistent functional limitation, like high traffic demand, road capacity, and conflict scenarios in some locations. This implies that such geospatial locations are frequently areas with high risk and have persisted in enduring RTA severity unless certain measures of safety are taken.

The detected spatial hot spot clustering of high severity of RTAs parallels previous studies that demonstrated hot spot evolution over metropolitan areas. The observed clustering of RTA severity over these areas reflects the high consistency with previous findings of the literature, which reported that hot spot clusters of RTAs are mainly concentrated along main urban roads and highways due to factors like dense population and high traffic volumes [26], [31]. The great similarity proposes that the patterns of spatial RTA are mainly impacted by urban structure and characteristics of the road network. However, due to the local conditions of traffic management and differences in urban intensity development, this study revealed higher hot spot cluster concentrations compared to some other literature.

The number of statistically significant hot spot clusters (confidence levels at 90%, 95%, and 99%) increased from 29,889 in 2020 to around 107,601 in 2021. The resulting values should not be regarded as independent or discrete hot spots across the selected study area. During the road segmentation process, they signify the combined number of statistically related segments of the highway that have been produced. It is worth noting that parallel segments along the same route direction might simultaneously achieve the criteria of statistical significance, resulting in much higher counts of hot spots compared with conventional area-based techniques of clustering. As a result, the resulting numbers reveal the geospatial resolution of network-based analysis rather than multiple discrete clusters of hot spots. A superior spatial resolution of a network-based segmentation methodology is demonstrated by high results, especially when each segment of road is evaluated separately as a unique analytical unit. This increase is attributed to the usage of the road segmentation technique as a spatial unit of analysis. The selected scale of segmentation enabled the reliable detection of regional versus geospatial changes of severity across the highway networks. However, a larger resolution of segmentation tends to enhance the total count of investigated spatial units, which leads to an increase in the frequency of statistically significant segments of hot spots compared with finer methods of aggregation.

Each segment is evaluated individually after the aggregation process of RTA severity, allowing the identification of differences over the road networks. As a result, more occurrences of localized events can be revealed rather than generalised patterns of the region. At the same time, the cold spots of 2021 decreased by around 50% compared to the previous year of 2020. The identified hot spots were correlated with high positive $z$ scores (i.e., $z >$ 2.58) with regard to the $p$-values less than or equal to 0.01($p \leq$ 0.01), signifying statistically significant clustering at the confidence level of 99%. For 2020, the high value of $z$-score ($z$ = 136) illustrates a robust severity clustering of RTAs, with $p <$ 0.01 confirming that the reflected pattern is unlikely to be due to random variation, which is statistically significant. The considerably high $z$-score values recognized in this methodology are due to the utilization of high-resolution road segments, in addition to the vast amount of aggregated road networks during G$_i^*$ analysis. The associated $p$-values were less than or equal to 0.01($p \leq$ 0.01), confirming the statistical significance at 99% confidence level. The situation of 2021 is similar to the previous year, when the highest z -score values were up to around 202, and that reflects an intense spatial intensification of high severity of RTAs, depicting a highly non-random clustering pattern. Also, $p$-values are less than or equal to 0.01, identifying robust evidence against the null hypothesis of spatial arbitrariness. The average ASI decreased from 3 in 2020 to 2.5 during 2021. The maximum ASI value in 2020 was 38, while the maximum ASI value in 2021 was 52.

The situation during 2022 revealed a noticeable increase in hot spot clusters (confidence levels of 90%, 95%, and 99%), with 11.1% with 119,556 clusters compared with the previous year. These results reflected further spatial clustering of high-severity clustering. At the same time, the cold spot clusters decreased by around 25% compared to 2021. Both results of $z$-scores and $p$-values for 2022 revealed around 279 combined $p <$ 0.01, which indicated a robust deviation from spatial randomness and confirmed the high severity of RTA hot spot clusters. This also indicates a robust and statistically significant clustering of high-severity RTAs. The average ASI for 2022 was 2.47, with a slight decrease compared to the previous year, 2021. However, the maximum ASI value during 2022 was the highest among all four years, with 111, reflecting the presence of the highest critical accident clusters.

It is worth noting that the observed increase in RTA severity over the study period is strongly consistent with previous findings of the literature, which link with the increasing levels of severity due to the high traffic volumes, in addition to the complexity of roadways [35], [41]. However, some other literature reported contrasting results, when the severity can be reduced by improving safety measures. These conclusions may be related to the road infrastructure variations. Besides the overall correlation with the demand of traffic, previous literature reveals that the severity of RTA is highly impacted by spatiotemporal factors such as composition of traffic volumes, environment of roads, and fluctuation of driving conditions. These factors would intensify RTA impact in complicated metropolitan areas [16], [30]. In addition, literature investigating spatiotemporal hot spot patterns of RTA revealed that the correlation between the dynamics of traffic flow and attributes of infrastructure majorly influences the severity distribution among various locations [33], [36]. Meanwhile, some studies reported a noticeable reduction in severity due to the improvements in traffic control measurements, enforcement strategies, and infrastructure enhancements [6], [16], [25]. This trend reveals the major impact of local safety policies and implementation levels for severity mitigation over time [12], [18].

It is worth mentioning that by 2021, hot spot areas expanded beyond main metropolitan spatial areas and covered some adjacent suburban counties, besides some secondary urban areas. For instance, in 2022, spatial hotspot areas emerged along major interstate highways, arterial roads, and commercial roads, especially in the central and southern areas of the state. These outcomes reflect the statistical significance of RTA severity patterns. Moreover, the outcomes of 2023 reflected intense and persistent hot spots, which expanded across the entire state of Missouri, with emerging new significant clusters, signifying augmented severity of RTA compared to previous years’ conditions, as in Figure 10.

This emergence of hot spot clusters signifies a shift in the spatial distribution of traffic risk over some spots that were considered less critical previously. This trend might be accompanied by the dispersion of urban activities, where population growth introduces new travel patterns with new trip generations. Moreover, the increased dependence on regional connectivity routes may intensify traffic flow out of conventional metropolitan areas, which contributed to the emergence of new high-risk areas. The rise of hot spot clusters may also indicate conflicts between capacity and the increasing demand for secondary road networks. The detected geospatial hot spot expansion onto suburban and secondary urban areas might be associated with the increased demand for regional transportation modes, besides the overall expansion of urban development exceeding metropolitan regions. Multiple recently identified segments of hot spots were inspected along interstate highways in addition to some major arterials. This indicates the spatial expansion of traffic activity, which leads to a high generation of fatal RTAs.

In 2023, the hot spot clusters (confidence levels of 90%, 95%, and 99%) increased to around 62.4% compared with the previous year, 2022, indicating a substantial evolution of high-severity clustering over the analysis area. This increase is attributed to the presence of around 194,286 hot spot clusters. Still, the $z$-scores recorded values more than 2.58 and $p$-values less than 0.01 ($p <$ 0.01). This confirms strong clustering of high-severity RTAs. The average ASI value recorded was the lowest at 2 during 2023. Also, the maximum ASI value was 39.

It is worth mentioning that the hot spot spatial expansion of RTA severity in this study aligns with previous literature, indicating that hot spot clusters extend along major corridors due to increasing mobility demand on RTA severity distribution [25]. Also, this study illustrates an ongoing adjustment in the hot spot level of intensity, attributed to the continuous population growth, besides traffic pattern evolution. The outcomes of the study reveal not only the spatial clustering of RTA severity, but also its temporal variation over an annual temporal scale. This technique enables the emergence of hot spots, which add an insight different from the literature that considers only spatial assessment during a single time period. Generally, the outcomes of this study demonstrate robust agreement with previous studies related to spatiotemporal patterns assessment of RTAs, besides extending the literature by emphasising the dynamic development of hot spots over time. In contrast to other literature, this study offers a more thorough comprehension of RTA severity distribution and its emergence over an annual scale of analysis [44], [45], [46].

The main objective of this study was to capture the temporal variations of RTA hotspot severity and determine statistically spatiotemporal clusters in Missouri. By using GIS-based geospatial analysis, this study illustrates the limitations of traditional accident research, emphasising the risk of injuries and fatalities. Driver behavior, specifically speeding, aggressively driving, and substance influence, all of which, was identified as a major contributor to the emergence of accident severity. The spatiotemporal analyses revealed that the RTAs in Missouri are significantly clustered, with persistent high-severity hotspots over major metropolitan areas like St. Louis and Jackson counties, alongside major interstates. The hotspot expansions into suburban areas beyond 2021 indicated a geographical increase in risk of severity. These trends of persistent patterns reveal critical insights for long-term safety management and transportation planning. Besides the statistically significant and persistent patterns of hotspot RTA severity, the statistics of G$_i^*$ can be beneficial for long-term hotspot analysis for traffic safety enhancement. The findings of this study can be useful for engineers, planners, and decision-makers to prioritise safety measures. All metropolitan areas that have persistent hotspots can benefit from targeted enforcement, such as increasing police patrols and enforcing speed regulations, or even structural enhancements such as redesign of junctions and/or signal optimization. Furthermore, the variation of temporal scales allows for better seasonal planning. However, some limitations can be inferred from this study. For instance, due to dataset constraints, the correlation between roadway geometry and RTA spatial locations was not inspected. Moreover, the uncertainty of the geospatial location of incidents may affect cluster identification. Therefore, further research should consider diverse smaller temporal scales, severity-weighted analysis, and machine learning frameworks to predict RTA severity dynamics better and investigate the impact of spatial parameter variance.