Study framework
In this study, MCDA was applied to assess the risk areas for FMD in China based on FMD outbreak data between 2010 and 2020 (Figure 1). We provided detailed data on outbreaks during this period from the Ministry of Agriculture of the People's Republic of China (MAPRC) and OIE.
1. We used Web of Science, PubMed and Google Scholar to select and evaluate the risk factors associated with the occurrence and spread of FMD. The five-Likert scale is used to select the ten most important spatial risk factors.
2. Analytic hierarchical process (AHP) was applied to assess the weights occupied by ten risk factors in FMD outbreaks.
3. Weighted linear combination (WLC) was used to transform the empirical data into risk map after standardization of all risk factors.
4. Sensitivity and uncertainty analysis was performed on the risk map.
Data transformation and standardization
We selected fifty risk factors from the literature and self-assessment (Supplement 1), these factors were obtained from actual surveys and analyses. The selected risk factors were classified and scored using the five-Likert scale [23]. Ten professional veterinarians participated in this work, including three Ph.D. and seven M.S. Ten most important risk factors for FMD outbreaks were selected, including livestock distribution, environment, transportation and location conditions.
Density data for buffalo, cattle, pigs, goats and sheep were collected to the extent of mainland China, and resampling was applied to transform them into map layers with a spatial resolution of 2.5 min. Proximity analysis was used to calculate the distance of each region in mainland China from outbreak points, national boundaries, livestock markets and slaughterhouses. The data were transformed into map layers with distances in km to each location. The kernel density estimation (KDE) was used to calculate densities of major railroads and roads in mainland China. Complex distance was used to measure local density changes in this method.
In this study, five types of relationships were proposed (linear, increasing sigmoidal, decreasing sigmoidal, symmetrical, and bi-directional). Based on previous literature for the relationship between these risk factors and FMD, sigmoidal was finally determined [15]. Standardized spatial layer maps of all risk factors were calculated by fuzzy membership functions at a scale of 0–1 (unsuitable to perfectly suitable, Supplement 2).
All geographic data were calculated and transformed using ArcGIS 10.2 (ESRI, Redlands, CA, USA). The standardization of the risk maps was performed using IDRISI (Clark Labs, Worcester, USA).
Risk factors weight
AHP was used to determine the weights of risk factors. The factors were grouped according to their characteristics. The elements of the same level were dominated by the elements of the previous level. It dominated certain elements of the next level. The judgment matrix was constructed by comparing the importance of all factors at each level using a 9-level scale. "1" means that the two elements are of equal importance; "3" means that the former is moderately important than the latter; "5" means that the former is strongly important than the latter; "7" means that the former is very strongly important than the latter; "9" indicates that the former is extremely important than the latter; 2, 4, 6 and 8 are the intermediate values of the above judgments. The importance of all factors was determined according to the scores in 2.2. For weight calculation and hierarchical single ranking, we applied the square root to calculate the eigenvectors and eigenvalues. The resulting eigenvectors are the weight ranking of the factors. Equation 1 is given as follows:
Where Mi is the product of the elements in the i-th row in the judgment matrix P; is the geometric mean of Mi; Wi is the weight of the i-th factor; W is the eigenvector; λmax is the maximum eigenvalue. The judgment matrix was constructed by comparing the scores of two elements. Equation 2 is given as follows:
Where CI is the consistency of the indicator; CR is the relative consistency of the indicator; λmax is the maximum eigenvalue; RI is the average random consistency indicator.
CR < 0.10 indicates a high level of consistency in the pairwise comparisons. However, if CR ≥ 0.10, the original values in the pairwise comparison matrix should be reconsidered and revised. R software version 4.0.3 (R Foundation for Statistical Computing, Vienna, Austria) was used for the application of AHP, and we provided the R code in this study (Supplement 3).
Risk map
A risk map was produced from the selected spatial risk factor layers and weighted using WLC. Equation 3 is given as follows:
Where n is the number of risk factors; w is the weight of factor I; v is the value function of risk factors at level i(ai). R is the total value of risk factors at each cell level in the end [24]. ArcGIS 10.2 was used to map and classify the risk map.
Sensitivity analysis and uncertainty
The one-at-a-time (OAT) was used to conduct sensitivity analysis. The change in the weight value of only one factor at a time (other factors remain unchanged as much as possible) reflects the degree of influence and regularity of the single factor weight change for results. For the adjustment weight (wa) of each main factor, the rate of change was between -20% and 20%, and the step size was 1%. The weights of other factors (wi) were calculated as follows (Equation 4):
Where Wi0 is the initial weight of each risk factor; Wa0 is the initial weight of the main change risk factor. The sum of all weights is equal to 1 [25].
The mean of absolute change rates (MACRs) of the adjusted-weight risk maps were calculated as follows (Equation 5):
Where Rk is the adjusted-weight risk map; R0 is the initial risk map; N is the number of pixels [15].
The uncertainty map was the standard deviation of the risk maps generated after all factors changed the weights [26]. The weights were randomly selected for each iteration through Monte Carlo sampling (between ±20%), and the process was repeated 400 times. The total number of output risk maps was calculated by multiplying the number of iterations with the number of criteria. Finally, 400 weight-adjusted values per criteria throughout this range were combined to compute 4000 risk maps [27].
Risk map validation
The FMD outbreak data in mainland China from 2018 to 2020 were used to test the accuracy of the map. Areas with coefficients greater than 0.7 in the risk map were considered to have an outbreak of FMD. Combining predicted and actual data, the receiver operating characteristic (ROC) analysis was used to assess the predictive ability of the risk map. The R package-InformationValue was used to calculate the ROC.