*2.2. Methods*

In the "Spatial Autocorrelation" tool of ArcMap, the Global Moran's *I* statistic, based on feature locations and attribute values [36], was adopted to measure the degree of global spatial autocorrelation (or, roughly speaking, to decide whether there was a spatial cluster effect) for different kinds of war in historical China. Before the tool can be applied, a grid containing the point data of battlefields was generated to facilitate the statistical analysis. In this study, 100 km was selected as the length of each square cell (Figure 2, and the reason for this is detailed in Supplementary Materials). Then, the number of battles in each cell of the grid was counted. Next, the Global Moran's *I* statistic was calculated according to 2020

$$I = \frac{n}{s\_0} \frac{\sum\_{i=1}^{n} \sum\_{j=1}^{n} w\_{i \cdot j} (\mathbf{x}\_i - \overline{\mathbf{x}}) \left(\mathbf{x}\_j - \overline{\mathbf{x}}\right)}{\sum\_{i=1}^{n} \left(\mathbf{x}\_i - \overline{\mathbf{x}}\right)^2} \tag{1}$$

where *n* is the total number of cells, *xi* and *xj* are the counted battle number in cells *i* and *j* respectively, *x* is the mean value of battle number in all cells, and *wi,j* denotes the proximity between *i* and *j*. When *i* and *j* are adjacent, *wi,j* = 1; otherwise, *wi,j* = 0. *S0* denotes the aggregate of all spatial proximity as

$$S\_0 = \sum\_{i=1}^{n} \sum\_{j=1}^{n} w\_{i,j} \tag{2}$$

The *z*-score for the Global Moran's *I* statistic is computed as

$$z = \frac{I - \mathbb{E}[I]}{\sqrt{\mathbb{V}[I]}} \tag{3}$$

where E[*I*] and V[*I*] are the expectation and variance, respectively, with the formulas

$$\mathbb{E}[I] = -1/(n-1) \tag{4}$$

$$\text{V}[I] = \text{E}[I^2] - \text{E}[I]^2 \tag{5}$$

**Figure 2.** Figure 2. Battle location and the grid (100 × 100 km per cell) adopted in this study.

The index ranges between −1 and 1. The spatial autocorrelation is positive (negative) if the value is larger (smaller) than 0. The higher (lower) the value, the more clustered (dispersed) the war. There is no spatial autocorrelation (i.e., random distribution) when the value is equal to 0.

As the Global Moran's *I* becomes positive, it is possible that there are some hot spots in which battles are clustered in space. Thus, in this research, EHSA was used to visualize the foci of war. The operational procedure of EHSA is realized in two major steps. First, a space–time cube with network common data form (.nc format) was created by running the tool "Create Space Time Cube" in the package "Space Time Pattern Mining Tools" of ArcMap. Instead of real 'points' (e.g., battle coordinates), the hot spots are an array of grids in which point data are aggregated and counted. Each cell (or bin) must have the same size. In this study, 100 km was chosen as the length of the side for each bin, because a bin less than 100 km (e.g., 50 km) would lead to innumerable cells, each of which only covers a small area. In this case, the bin number that contains battlefields may be small, whereas others may be largely blank. By contrast, if the length of the side is more than 100 km (e.g., 200 km), the resultant patterns may be too coarse to decipher, since only limited amounts of bins reside in the study area (Figure S1). Therefore, the medium size of 100 × 100 km was suitable for this task.

During the first step, the Mann–Kendall trend test (or M–K test) was automatically conducted to evaluate the overall trend for all bins in the cube. It is a non-parametric test used to analyze whether data are consistently increasing or decreasing over time [37,38]. More details about the M–K test are provided in Supplementary Materials. Second, the cube was input in the Emerging Hot Spot Analysis tool, and the Getis-Ord *Gi\** statistic (i.e., the traditional hot spot analysis) was calculated for each bin. The formula for obtaining the *Gi\** is given as [39]

$$\mathbf{G}\_{i}^{\*} = \frac{\sum\_{j=1}^{n} w\_{i,j} \mathbf{x}\_{j} - \overline{\mathbf{x}} \sum\_{j=1}^{n} w\_{i,j}}{\mathbf{S} \sqrt{\frac{n \sum\_{j=1}^{n} w\_{i,j}^{2} - \left(\sum\_{j=1}^{n} w\_{i,j}\right)^{2}}{n-1}}} \tag{6}$$

where

$$\overline{\mathbf{x}} = \frac{\sum\_{j=1}^{n} \mathbf{x}\_{j}}{n} \tag{7}$$

$$S = \sqrt{\frac{\sum\_{j=1}^{n} x\_j^2}{n} - (\overline{x})^2} \tag{8}$$

Once the analysis was completed, each bin had an associated *z*-score, *p*-value, and hot spot classification. With the resultant trend, *z*-score, as well as *p*-value for each location, the EHSA tool categorized eight kinds of hot spots, which are listed in Table 2.

**Table 2.** Definitions of various hot spots in EHSA.


(http://desktop.arcgis.com/en/arcmap/latest/tools/space-time-pattern-mining-toolbox/learnmoreemerging.htm).

In terms of the spatiality of war, EHSA surpasses the traditional hot spot analysis because bins with high frequencies of battles surrounded by high values can be determined and various hot spots are categorized based on the variation trends over time. Accompanying hot spots, however, the spatiotemporal patterns of cold spots (i.e., clusters of low values) are generated, which may be

meaningless since war hot spots were prioritized in this study. Only significant hot spots are presented (*p* < 0.05).
