Identifying Key Factors of Hazardous Materials Transportation Accidents Based on Higher-Order and Multilayer Networks

This paper focuses on the application of higher-order and multilayer networks in identifying critical causes and relationships contributing to hazardous materials transportation accidents. There were 792 accidents of hazardous materials transportation that occurred on the road from 2017 to 2021 which have been investigated. By considering time sequence and dependency of causes, the hazardous materials transportation accidents causation network (HMTACN) was described using the higher-order model. To investigate the structure of HMTACN such as the importance of causes and links, HMTACN was divided into three layers using the weighted k-core decomposition: the core layer, the bridge layer and the peripheral layer. Then causes and links were analyzed in detail. It was found that the core layer was tightly connected and supported most of the causal flows of HMTACN. The results showed that causes should be given hierarchical attention. This study provides an innovative method to analyze complicated accidents, which can be used in identifying major causes and links. And this paper brings new ideas about safety network study and extends the applications of complex network theory.


Introduction
Hazardous materials transportation safety, especially road transportation, is closely tied to people's lives and property and getting a growing concern. For example, the National Roadway Safety Strategy, published by the U.S. Department of Transportation, proposed zero roadway fatalities and serious injuries [1]. The State Council of China printed the National Road Traffic Safety Plan for the "Fourteenth Five Year Plan" and put it forward to decrease fatalities and frequency of traffic accidents [2]. Road transportation of hazardous materials is considered to be the most dangerous to public health [3]. When it happened, nearby residents usually suffer the most. And its accident causes are complex by the coupling effect of many unsafe factors including humans, vehicles, the environment, etc., and the inherent danger of hazardous materials. At present, hazardous materials transportation accidents occurring on the road are maintaining a high frequency of growth [4,5]. For example, during the period 2013-2019, 2777 accidents involving the road transportation of dangerous goods occurred in China, except occurring in loading, unloading, storage and maintenance [5]. Therefore, it is urgent to curb these accidents.
Identifying key causes and relationships is quite important to prevent accidents. Accident causes often involve nonlinear relationships and complex dynamic processes. Various factors, such as human, equipment, vehicle, environment and management, interact with each other and increase the complexity of transportation accidents [6]. For example, the combination of heavy fog weather and fatigue driving can lead to serious accidents. However, existing studies on accident causation mostly analyzed the influence of singlefactors without considering the interactions of multiple factors [7][8][9]. The popular methods considering influences of accident forms, 9 accident types are classified into single factors, which contain collision, scrape, roll over, fall over, fire, leakage, explosion, poisoning and others. Based on our previous research on hazardous materials transportation accidents [37], Figure 1 shows the single factors of HMTAs. Multi-factors emphasize the combined actions of two or more single-factors, for example, a multi-factor E2H6 means combined actions of E02 (slippery road) and H06 (improper operation). Ultimately, 129 accident causal factors, containing 67 single-factors and 62 multi-factors, are found. Then, the 792 accidents are transformed into accident chains, which corresponded to the codes. The specific process is illustrated by two examples, as shown in Table 1.  Unsafe distance → multicar collision → leakage H09 → E07 → A06

Modeling HMTACN
Since time sequence and dependence are key characteristics of HMTAs, this paper builds HMTACN with the idea of non-Markov chains in higher-order network theory. With HMTAs' causal factors and the chains aforementioned, path dependency should be extracted, then a network can be constructed. Based on the BuildHON+ algorithm [27],

Modeling HMTACN
Since time sequence and dependence are key characteristics of HMTAs, this paper builds HMTACN with the idea of non-Markov chains in higher-order network theory. With HMTAs' causal factors and the chains aforementioned, path dependency should be extracted, then a network can be constructed. Based on the BuildHON+ algorithm [27], the process of this network contains two steps.
Step 1 is used to illustrate the path dependency extraction process and step 2 is used to illustrate the network construction process.
where, W(i → j) is the frequency from node i to j, and h denotes the set of nodes that node i may go at time t + 1.
the process of this network contains two steps.
Step 1 is used to illustrate the path dependency extraction process and step 2 is used to illustrate the network construction process.
where, ( → ) is the frequency from node to , and ℎ denotes the set of nodes that node may go at time t + 1. Step 2: The network construction process is illustrated by an example of five path dependencies, shown in Figure 3. First, the first-order nodes and their corresponding edges are formed by converting the first-order path dependencies. As shown in Figure 3b, H12 → H02 is attached to a network and converted into H12 and H02. Second, the secondorder node H02|H12 and two outgoing edges are formed by converting H02|H12 → E07 and H02|H12 → A01 in Figure 3c. Similarly, the second-order node E07|H02 and two outgoing edges are formed by E07|H02 → A05 and E07|H02 → A03. Finally, the edges of all higher-order nodes are reconnected. Node H02 is derived from H12 to form the secondorder node H02|H12, and reconnects H12 → H02 to "H02|H12" in Figure 3d. Similarly, the edge H02|H12 → E07 is reconnected to the second-order node "E07|H02" in Figure 3e.

H12→H02→E07→A05 H09→E07→A06
H12 Step 2: The network construction process is illustrated by an example of five path dependencies, shown in Figure 3. First, the first-order nodes and their corresponding edges are formed by converting the first-order path dependencies. As shown in Figure 3b, H12 → H02 is attached to a network and converted into H12 and H02. Second, the secondorder node H02|H12 and two outgoing edges are formed by converting H02|H12 → E07 and H02|H12 → A01 in Figure 3c. Similarly, the second-order node E07|H02 and two outgoing edges are formed by E07|H02 → A05 and E07|H02 → A03. Finally, the edges of all higher-order nodes are reconnected. Node H02 is derived from H12 to form the second-order node H02|H12, and reconnects H12 → H02 to "H02|H12" in Figure 3d. Similarly, the edge H02|H12 → E07 is reconnected to the second-order node "E07|H02" in Figure 3e.
Finally, HMTACN contains 243 nodes and 545 edges. Nodes represent causal factors with 124 first-order nodes and 119 higher-order nodes; the maximum order is four. Edges represent interrelationships of nodes. The HMTACN is visualized by Gephi, shown in Figure 4a. Gephi is an open-source network analysis and visualization software that provides a visual way to analyze complex networks. Figure 4b shows part of HMTACN, which clearly displays detailed interactions of factors. There are five first-order nodes (A02|, E07|, A6HM3|, H06| and V04|) representing causal factors and twenty higher-order nodes representing path dependencies. For example, the second-order node A01|H06 represents a second-order dependency with H06 → A01, and the third-order node V04|E07.A03 represents a path dependency with A03 → E07 → V04.  Finally, HMTACN contains 243 nodes and 545 edges. Nodes represent causal factors with 124 first-order nodes and 119 higher-order nodes; the maximum order is four. Edges represent interrelationships of nodes. The HMTACN is visualized by Gephi, shown in Figure 4a. Gephi is an open-source network analysis and visualization software that provides a visual way to analyze complex networks. Figure 4b shows part of HMTACN, which clearly displays detailed interactions of factors. There are five first-order nodes (A02|, E07|, A6HM3|, H06| and V04|) representing causal factors and twenty higherorder nodes representing path dependencies. For example, the second-order node A01|H06 represents a second-order dependency with H06 → A01, and the third-order node V04|E07.A03 represents a path dependency with A03 → E07 → V04.  Finally, HMTACN contains 243 nodes and 545 edges. Nodes represent causal factors with 124 first-order nodes and 119 higher-order nodes; the maximum order is four. Edges represent interrelationships of nodes. The HMTACN is visualized by Gephi, shown in Figure 4a. Gephi is an open-source network analysis and visualization software that provides a visual way to analyze complex networks. Figure 4b shows part of HMTACN, which clearly displays detailed interactions of factors. There are five first-order nodes (A02|, E07|, A6HM3|, H06| and V04|) representing causal factors and twenty higherorder nodes representing path dependencies. For example, the second-order node A01|H06 represents a second-order dependency with H06 → A01, and the third-order node V04|E07.A03 represents a path dependency with A03 → E07 → V04.

Degree and Strength
Degree: k i The degree of node i is the number of edges connected to it and is calculated as shown in Formula (2): where: N-the total number of nodes in the network; a i,j -a i,j = 1, there is an edge between node i and j, otherwise, a i,j = 0.

Strength: S i
The strength of node i is the sum of the weights with the edges directly connected to it, which is called "accident causal flows" in this paper. It is calculated as shown in Formula (3): where: w i,j -the weights of the edges between node i and j.

Weighted k-Core Decomposition
Since the weighted k-core decomposition was proposed by Antonios Garas et al. [38], it has often been applied in multilayer to analyze the multilayer structure of networks [33,34]. And the specific stratification depends on the study subjects and data characteristics. To further analyze the causal factors and reveal the multilayer characteristics, the weighted k-core decomposition is used to analyze HMTACN in this paper. The calculation is shown in Formula (4): where: k i -the degree of node i; w ij -the weight between node i and its neighboring j; When the parameters α = β = 1, the degree and weight have a similar influence. The steps of the weighted k-core decomposition are as follows.
Step 1: Calculate the total degree k i of all nodes.
Step 2: Gradually remove nodes with k i = 1, which belongs to the same layer of the network.
Step 3: Gradually remove nodes with k i = 1, which may be completely disconnected from the main network.
Step 4: Remove the new nodes with k i < k i and keep repeating this process until k i ≥ k s + 1 for all nodes in the network. Then mark the removed nodes as k s and k s = k s + 1.
Step 5: The algorithm stops when all the nodes are marked with k s .

Measurements in HMTACN
In HMTACN, several nodes of different orders may correspond to one causal factor in accident chains. When the HMTACN is layered using the weighted k-core decomposition, nodes of different orders may be assigned to different layers. The following are the corresponding measurements.
(1) Layer-degree: k il The layer-degree of the node i is the sum of the degrees of all nodes of different orders corresponding to causal factor i in layer l. It is calculated as Formula (5).
It should be emphasized that HMTACN is a directed network, the layer-in degree and layer-out degree are proposed. The layer-in degree of node i is the sum of the entry degrees of all nodes of different orders corresponding to causal factor i in layer l. It is calculated as Formula (6).
The layer-out degree of node i is the sum of the exit degrees of all nodes of different orders corresponding to causal factor i in layer l. It is calculated as Formula (7): where: M-the total number of all nodes of different orders corresponding to node i in layer l; k m il -the degree of the mth node corresponding to node i in layer l.
The layer-strength of the node i is the sum of the strengths of all nodes of different orders corresponding to causal factor i in layer l. It is calculated as Formula (8): where: S m il -the strength of the mth node corresponding to node i in layer l.

(3) Layer-ks: k sil
The layer-k s of the node i is the sum of the k s of all nodes of different orders corresponding to causal factor i in layer l. It is calculated as Formula (9): (9) where: k m sil -the k s of the mth node corresponding to node i in layer l.

(4) The ratio of causal node connection and causal flow within and between the layers
In order to deeply analyze network topology and node characteristics for different layers, the ratio of causal node connection and the ratio of causal flow are defined from intra-layer and inter-layer aspects, as shown in Table 2. Table 2. The ratio of causal node connection and causal flow within and between the layers.

The Ratio of Causal Node Connection The Ratio of Causal Flow
Inter-layer Where: N layer -the number of neighbors of node i within the layer, for example, N core means the number of neighbors in the core layer, N bridge means the number of neighbors in the bridge layer, and N total means the total number of neighbors of node i within and between the layers. F layer -the number of causal flows of node i with neighbors in the layer, for example, F core means the number of causal flows with the neighbors in the core layer of this causal factor, F bridge means the number of causal flows with the neighbors in the bridge layer of this causal factor, F total means the total number of causal flows that start or end within this causal factor. R a in -the ratio of causal nodes within the layers; R a out -the ratio of causal nodes between layers; R f in -the ratio of causal flow within the layers; R f out -the ratio of causal flow between the layers; R a cb -the ratio of causal node connection between the core layer and bridge layer, R a cp , R a bc , R a bp are similar; R f cb -the ratio of causal flow between the core layer and bridge layer,

Multilayered Structure of HMTACN
According to the weighted k-core decomposition mentioned in Section 3.3.2, the multilayer structure of HMTACN is proposed. Figure 5 shows the result of the layers. Where: -the number of neighbors of node within the layer, for example, means the number of neighbors in the core layer, means the number of neighbors in the bridge layer, and means the total number of neighbors of node within and between the layers.
-the number of causal flows of node with neighbors in the layer, for example, means the number of causal flows with the neighbors in the core layer of this causal factor, means the number of causal flows with the neighbors in the bridge layer of this causal factor, means the total number of causal flows that start or end within this causal factor.
-the ratio of causal nodes within the layers; -the ratio of causal nodes between layers; -the ratio of causal flow within the layers; -the ratio of causal flow between the layers; -the ratio of causal node connection between the core layer and bridge layer, , , are similar; -the ratio of causal flow between the core layer and bridge layer, , , are similar.

Multilayered Structure of HMTACN
According to the weighted k-core decomposition mentioned in Section 3.3.2, the multilayer structure of HMTACN is proposed. Figure 5 shows the result of the layers. Theoretically, the nodes corresponding to each k s can be used as a separate layer, but too many nodes are not conducive to the study. Inspired by the division of multilayer networks in the traffic domain [29,33,34], HMTACN can be divided into three layers: the core layer, the bridge layer and the peripheral layer by the distribution characteristics of k s . As Figure 5 shows, the first gap appears and k − shell shows a continuous decreasing trend when k s ≥ 15. And the corresponding nodes within this range are assigned to the core layer, which indicates 19 nodes are the core nodes in HMTACN; when k s ≤ 2, the k − shell reach the maximum value, and the corresponding nodes are assigned to the peripheral layer; and when 3 ≤ k s < 15, these corresponding nodes are assigned to the bridge layer.
The multilayer structure of HMTACN is shown in Figure 6 by Gephi. The node size represents the number of factors connected to it. Here green represents the core layer, orange represents the bridge layer, and purple represents the peripheral layer. The color of the edge is the color of the target node. Since the existence of higher-order nodes in the multilayer structure of HMTACN, one cause may correspond to several different nodes and be assigned to different layers. Also, multiple nodes in a layer may correspond to a cause.
For example, H02 corresponds to the H02|H12, H02|E07 and H02|V10 in HMTACN, where H02 is part of the core layer, H02|H12 and H02|E07 are part of the bridge layer and H02|V10 is part of the peripheral layer. Therefore, node H02 belongs to three layers. It is found that 19 causal factors belong to two layers accounting for 15.32%, and 6 causal factors belong to three layers, accounting for 4.84%.

layer.
The multilayer structure of HMTACN is shown in Figure 6 by Gephi. The node size represents the number of factors connected to it. Here green represents the core layer, orange represents the bridge layer, and purple represents the peripheral layer. The color of the edge is the color of the target node. Since the existence of higher-order nodes in the multilayer structure of HMTACN, one cause may correspond to several different nodes and be assigned to different layers. Also, multiple nodes in a layer may correspond to a cause. For example, H02 corresponds to the H02|H12, H02|E07 and H02|V10 in HMTACN, where H02 is part of the core layer, H02|H12 and H02|E07 are part of the bridge layer and H02|V10 is part of the peripheral layer. Therefore, node H02 belongs to three layers. It is found that 19 causal factors belong to two layers accounting for 15.32%, and 6 causal factors belong to three layers, accounting for 4.84%.  In different layers, the nodes and edges present different characteristics, as shown in Table 3. In the core layer, there are 16 causal factors. The value of S il has a wide range from 34 to 846. Actually, 14 causal factors have an average layer strength between 34 and 247, and 2 causal factors have much higher S il with 704 and 846, respectively. Obviously, the core layer has the largest k sl value. In addition, 70 edges in the core layer bear 1511 causal associations, with an average of 21.59 per edge, accounting for 63% of the network. In the bridge layer, there are 40 causal factors and the value of S il have a small variation, from 3 to 37. The 67 edges in the bridge layer bear 143 causal associations, with an average of 2.13 per edge. There are 98 causal factors in the peripheral layer and the value of S il is only between 1 and 5, indicating that each causal factor in the peripheral layer appears less frequently. The 57 edges in the peripheral layer bear 61 causal associations, with an average of 1.07 per edge. The edges between the three layers can reflect the inter-layer relationship of causal factors. Clearly, the relationship between the core and bridge layers is prominent in Table 3.

Importance of Causes and Links in HMTACN
4.2.1. Importance of Causes in HMTACN As aforementioned, HMTACN is divided into the core layer, bridge layer and peripheral layer based on k s . The k s value of a node can reflect the number of connected nodes, the frequency of related edges and the ability to connect other layers. Usually, the larger the k s is, the more important the node is. Therefore, causes are classified into seven levels (I-VII) by dividing the nodes within different layers, shown in Table 4. In Table 4, there are 5 causes in level I with nodes belonging to core, bridge and periphery layers simultaneously. And the k s value in level I is the largest. This indicates the causes in level I are the most important causal factors in HMTACN. Among them, the three higher-order nodes that represent dependencies A01|H09, A01|H06 and A03|H06 belong to the core layer, which cannot be captured by traditional complex networks. In level II, only A05 belongs to core and bridge layers and its k s value is 42. However, in level III, the k s value is 68.83 which is relatively larger than level II. This is because the k s value of V04 is 157 and V06 is 145 which increases the average k s value in level III. In levels I, II, III and IV, the k s value is large and node number is small relatively. Therefore, the 16 causes within levels I to IV are core causal factors in the HMTAC system.
From level V to level VII, k s values decrease and node numbers increase successively. Further, the causes in level V are mainly vehicle factors, and nodes in level VI cover almost all causal factor types. In level VII, multi-factors account for 71.62%, and 47.17% are related to environmental factors.

Importance of Links in HMTACN
The links importance of HMTACN is measured by w l (the average weight of the edges in layer l) shown in Table 5. The greater w l , the greater accident causal flow, the more important the link is. As shown in Table 5, at level I, the values of w l is the most with 70 links accounting for 12.84%. Thus, it is important to focus on these links. From level II, the values of w l decrease sharply. This means these links undertake fewer cause connections. At level II, the number of links is the largest, accounting for 32.29%. Level III has fewer links and lower weights. The values of w l from levels IV to VI are similar, while the number at level IV is much larger than at levels V and VI. This indicates most of the bridge and peripheral layer causes tend to connect with the core layer, while the links between the bridge and peripheral layer are sparse.
In order to reveal the significant links of HMTACN, the study explores the links within level I, shown in Table 6. In the type of vehicle causes, V04, V03 and V02 are easily leading to A06 accidents. Especially, the link frequency of V04 → A06 is the highest, accounting for 24.42%. In the type of accident causes, A01 → V04 is the most important link, the value of weight is 216. In addition, A03 is easily leading to V04 and A06. In the type of human cause, H12 → A01 is clearly the most important link, the value of weight is 210. In addition, H02 and H06 are easily leading to A03. In the type of environmental causes, the weights of E09 → E02, E07 → V04 and E02 → E07 are similar and smaller, indicating that the environmental cause alone is not the most important causal factor.

Core Layer
The core layer is the crucial layer in HMTACN. To discover the characteristics of the core layer, 16 causal factors are further analyzed by calculating the ratio of causal node connection and causal flow within the core layer and between the core layer and other layers. Table 7 shows the detailed calculations with R a in , R a out , R a cb , R a cp , R f in , R f out , R f cb , R f cp , which have been defined in Section 3.3.3.
As shown in Table 7, accident types acting as intermediaries remain larger K il in core layer, including A03 (roll over), A01 (collision), A06 (leakage) and A05 (fire). This indicates that these accidents occur most frequently and connect more factors, especially roll over, collision and leakage of hazardous materials. The unsafe status of the vehicle is relatively more in the core layer, including V04 (tank damaged), V03 (valve damaged), V02 (valve loosed), V06 (pipe rupture) and V14 (oil tank damaged), which mainly focus on equipment loading hazardous materials. It indicates that they connect more factors and the damage to vehicle facilities are key reasons in the core layer, especially V04. Unsafe behaviors of humans, like H12 (non-hazmat personal reasons), H02 (improper avoidance) and H06 (improper operation) are the main causes, especially H12. Also, environmental causes like E07 (multi-car collision), E02 (slippery road) and E09 (rain and snow weather)  To visualize the characteristics of the core layer, calculations of Table 6 are drawn in Figures 7 and 8. Figure 7 shows that R a and R f are not direct relation to the node higher-order degree (k il ) in the core layer. Except for E09 (rain and snow weather), A05 (fire) and V03 (valve damaged), the causal factors show that R a in is much smaller than R a out , but R f in is much larger than R f out . This suggests that interactions between causal factors within the core layer are more likely to lead to accidents than with the bridge and periphery layers. In addition, E09 and A05 show that R a in > R a out and R f in > R f out , indicating that they tend to interact with causal factors in the core layer. It is worth noting that only the V03 R a in < R a out but R f in < R f out indicates that it is more closely connected to other layers. This can be attributed to both causal factors of E11 (high temperature) and E13 (other environmental reasons) in the bridge layer that are more likely to cause valve damage. It is found that the core layer has the strongest effect in HMTACN.  Figure 7 shows that and are not direct relation to the node higher-order degree ( ) in the core layer. Except for E09 (rain and snow weather), A05 (fire) and V03 (valve damaged), the causal factors show that is much smaller than , but is much larger than . This suggests that interactions between causal factors within the core layer are more likely to lead to accidents than with the bridge and periphery lay-  Figure 8 shows that the causal factors have the characteristics of R a cp < R a cb and R f cp < R f cb on the whole, indicating that the interaction between the core layer and the bridge layer is stronger than that between the periphery layer. Further, there are 4 edges with strength greater than 10 between the core layer and bridge layer, and 4 higher-order nodes that represent dependencies are generated: A01|H02, A03|E02, A01|A03, A01|E02. This illustrates two findings. Firstly, H02 (improper avoidance) → A01 (collision), E02 (slippery road) → A03 (roll over), A03 (roll over) → A01 (collision), E02 (slippery road) → A01 (collision) are the crucial causal links between the core layer and bridge layer. Similar findings have been found in previous studies [4]. Secondly, the dependency may have been missed by previous studies using complex networks. Such as A01|H02, collision usually depends on the generation of improper avoidance leading to chain accidents. Figure 7. (a) and of casual factors in core layer vs. layer-degree. (b) and of casual flows in core layer vs. layer-degree within the layer. Figure 7 shows that and are not direct relation to the node higher-order degree ( ) in the core layer. Except for E09 (rain and snow weather), A05 (fire) and V03 (valve damaged), the causal factors show that is much smaller than , but is much larger than . This suggests that interactions between causal factors within the core layer are more likely to lead to accidents than with the bridge and periphery layers. In addition, E09 and A05 show that > and > , indicating that they tend to interact with causal factors in the core layer. It is worth noting that only the V03 < but < indicates that it is more closely connected to other layers. This can be attributed to both causal factors of E11 (high temperature) and E13 (other environmental reasons) in the bridge layer that are more likely to cause valve damage. It is found that the core layer has the strongest effect in HMTACN.  on the whole, indicating that the interaction between the core layer and the bridge layer is stronger than that between the periphery layer. Further, there are 4 edges with strength greater than 10 between the core layer and bridge layer, and 4 higher-order nodes that represent dependencies are generated: A01|H02, A03|E02, A01|A03, A01|E02 This illustrates two findings. Firstly, H02 (improper avoidance) → A01 (collision), E02 (slippery road) → A03 (roll over), A03 (roll over) → A01 (collision), E02 (slippery road) → A01 (collision) are the crucial causal links between the core layer and bridge layer. Similar findings have been found in previous studies [4]. Secondly, the dependency may have been missed by previous studies using complex networks. Such as A01|H02, collision usually depends on the generation of improper avoidance leading to chain accidents.

Bridge Layer
Similar to the core layer analysis, Figure 9 shows the ratios within the bridge layer and between the bridge layer and the other two layers. In the bridge layer, most of the causal factors show R a in < R a out and R f in < R f out , indicating that the bridge layer is more inclined to interact with outside the layer. Also, most causal factors in the bridge layer have the characteristics of R a bc > R a bp and R f bc > R f bp on the whole (in Figure 10), indicating that they tend to interact with the core layer.

Bridge Layer
Similar to the core layer analysis, Figure 9 shows the ratios within the bridge layer and between the bridge layer and the other two layers. In the bridge layer, most of the causal factors show < and < , indicating that the bridge layer is more inclined to interact with outside the layer. Also, most causal factors in the bridge layer have the characteristics of > and > on the whole (in Figure  10), indicating that they tend to interact with the core layer.    To further explore important characteristics of causal factors in the bridge layer, these top 14 causal factors with larger have been analyzed, shown in Table 8. These 14 causal factors have 76 edges with the core layer, accounting for 59.38% of all edges between the bridge layer and core layer. In addition, the value of , , and are much larger than 0.5. Except for E13 (other environmental reasons) and E2H1 (slippery road and fast driving), M07 (inadequate safety check and maintenance) has the largest , followed by M08 (other management reasons). Causal factors connecting to M07 and M08 are V02 (valve loosed), V03 (valve damaged), V04 (tank damaged), V06 (pipe rupture), A05 (fire) and A06 (leakage), indicating that M07 and M08 need only one step to interact with these causal factors. But other causes require two or more steps. This indicates that the influence of the management factor is greater and the problem is more serious. In short, the interaction between the bridge layer and core layer is far closer than the periphery layer. To further explore important characteristics of causal factors in the bridge layer, these top 14 causal factors with larger R a bc have been analyzed, shown in Table 8. These 14 causal factors have 76 edges with the core layer, accounting for 59.38% of all edges between the bridge layer and core layer. In addition, the value of R a out , R f out , R a bc and R f bc are much larger than 0.5. Except for E13 (other environmental reasons) and E2H1 (slippery road and fast driving), M07 (inadequate safety check and maintenance) has the largest R a bc , followed by M08 (other management reasons). Causal factors connecting to M07 and M08 are V02 (valve loosed), V03 (valve damaged), V04 (tank damaged), V06 (pipe rupture), A05 (fire) and A06 (leakage), indicating that M07 and M08 need only one step to interact with these causal factors. But other causes require two or more steps. This indicates that the influence of the management factor is greater and the problem is more serious. In short, the interaction between the bridge layer and core layer is far closer than the periphery layer. In the bridge layer, the value of S il is less than 28 for most causal factors, except V16 (packaging issues), V12 (tire overheating), and H09 (unsafe distance). Therefore, these causal factors have been analyzed furthermore in Table 9. It is found that V16 only belongs to the bridge layer and has the largest value of K il and S il in the bridge layer, so the packaging issue is the most significant causal factor in the bridge layer. Moreover, only node V16 has a higher K il in than K il out , which indicates that many causal factors can lead to the packaging issues, such as H04 (improper braking), A03 (roll over), A01 (collision), V15 (other vehicle reasons) and E04 (poor road). In addition, the value of R f out and R f bc of V16, V12 and H09 are relatively larger, which illustrates that they are more closely connected with the core layer.

Periphery Layer
Similar to the core layer analysis, Figure 11 shows the ratio within the peripheral layer and between the periphery layer and the other two layers.
causal factors have been analyzed furthermore in Table 9. It is found that V16 only belongs to the bridge layer and has the largest value of and in the bridge layer, so the packaging issue is the most significant causal factor in the bridge layer. Moreover, only node V16 has a higher than , which indicates that many causal factors can lead to the packaging issues, such as H04 (improper braking), A03 (roll over), A01 (collision), V15 (other vehicle reasons) and E04 (poor road). In addition, the value of and of V16, V12 and H09 are relatively larger, which illustrates that they are more closely connected with the core layer. Similar to the core layer analysis, Figure 11 shows the ratio within the peripheral layer and between the periphery layer and the other two layers. As shown in Figure 11, the trends of R a and R f in the periphery layer are basically the same without obvious patterns. This is because these causal factors within the periphery layer appear less frequently in accident chains and most of them are multi-factors. Moreover, 68% of multi-factors are related to environmental factors. This explains that environmental factors often work in conjunction with other causal factors to cause accidents, so the abnormal environment should be alerted in time.
Furthermore, there are 44 causal factors with R a out = R f out = 1 in the periphery layer, which illustrates that these causal factors are only correlated with other layers. And among them, 31 causal factors only belong to the periphery layer, 19 of them are multi-factors, and most of them interact with E01 (downhill road), E02 (slippery road) and E03 (turning road).

Discussion
In this study, the method combining higher-order and multilayer networks is applied to the field of hazardous materials transportation accidents for the first time. A total of 792 accidents occurring on roads from 2017 to 2021 are analyzed. The hazardous materials transportation accident causation network (HMTACN) is constructed using a higher-order network and divided into three layers using weighted k-core decomposition: the core layer, bridge layer, and peripheral layer.
As a result, 16 key causes covering five types of causes are identified through the analysis of the core layer. In the bridge layer, the management factors including M07 (inadequate safety check and maintenance) and M08 (other management reasons), the vehicle factors including V16 (packaging issues) and V12 (tire overheating), the human factor including H09 (unsafe distance) are critical causes of accidents. The analysis of the peripheral layer indicated that environmental factors often contributed to accidents in conjunction with other factors. The important results of each layer will be discussed from the perspective of several cause types, including vehicle, human, environment, accident types and management factors.
(i) Vehicle factors such as V04 (tank damaged), V03 (valve damaged), V02 (valve loosed), V06 (pipe rupture) and V14 (oil tank damaged) contribute to hazardous material accidents, as the damage to vehicle facilities can lead to leaks and increase the probability of accidents. These factors are associated with the maintenance condition, aging, and manufacturing quality of the vehicles. Therefore, it is necessary to implement regular maintenance programs to ensure the integrity of vehicle tanks, valves, and pipelines. Also, promptly repair or replace any damaged components to prevent leaks, and strengthen quality control measures during the manufacturing process to ensure vehicles are built to high standards. In addition, V16 (packaging issues) and V12 (tire overheating) are also important causes of accidents. This is consistent with the findings of Ma et al. [39] that packaging problems are an important cause of hazardous materials transportation accidents. Therefore, the packaging of hazardous materials should be improved.
(ii) Human factors, including H12 (non-hazmat personnel reasons) by non-hazmat personnel, H02 (improper avoidance), H06 (improper operation) and H09 (unsafe distance) by hazardous material transport personnel, have a high likelihood of causing accidents. These are consistent with the findings of Oggero et al. [40] that valve failure in machinery and improper human operation are the major causes of accidents. Human factors are influenced by driver experience, training, awareness, and attitude. Therefore, developing comprehensive training programs for both non-hazardous personnel and hazardous material transport personnel is important. Also, it is necessary to emphasize emergency response protocols and defensive driving techniques and implement stringent screening and qualification processes to ensure that personnel have the appropriate competencies.
(iii) Environmental factors such as E07 (multi-car collision), E02 (slippery road) and E09 (rain and snow weather) often work in conjunction with other factors to cause hazardous material accidents. Poor road conditions increase the probability of human error. Therefore, it is necessary to establish effective communication channels to provide real-time weather updates and road condition alerts to drivers.
(iv) Different accident types, such as A03 (roll over), A01 (collision), A06 (leakage) and A05 (fire), have different causes and potential impacts. For example, roll over can cause traffic congestion and road blockages, collisions can result in vehicle damage and personal injuries, leakage can lead to environmental pollution and health risks, and fire can trigger explosions. Therefore, developing and enforcing strict safety standards covering the handling, storage and transportation of hazardous materials is crucial.
(v) Management factors play a significant role in HMTACN, and the lack of management leads to an increased probability of all factors and accidents [20]. Therefore, it is necessary to implement comprehensive safety checks, maintenance and quality control measures to ensure the integrity of vehicles and their components.

Conclusions
In conclusion, this study proposes a new accident causation method that combines higher-order and multilayer networks. Unlike previous methods, this approach considers the interactions and dependencies of factors, allowing for more accurate identification of key causes by considering the interactions and dependencies of factors. The method also conducts a quantitative analysis of complex accident causation systems, particularly when dealing with a large amount of accident data. Based on the findings, recommendations are discussed to mitigate the risks of hazardous material road transportation accidents. Implementing these measures can help reduce accidents and enhance the safety of hazardous material transportation.
According to this study, the significant causes and links should be given hierarchical attention due to their different roles in the hazardous materials transportation accident causation network. Overall, this study contributes to the understanding of hazardous material transportation accidents and provides valuable insights for developing effective preventive measures.
This paper considered the static indicators of hazardous materials transportation accidents causation network. Further research should focus on analyzing the complex interactions and propagations of causes and links. Additionally, the application of surrogate safety measures, such as traffic conflict analysis, for dynamic safety evaluation of hazardous materials transportation deserves further study. In the future, the major causes and relationships may be changed, so more detailed and updated dynamic data on hazardous materials transportation accidents should be studied continuously.