River ecosystem planktonic bacteria environment threshold identification method

Abstract

The invention relates to a river ecosystem planktonic bacteria environment threshold identification method, and belongs to the field of environmental ecology and computational biology. The method comprises the following steps: collecting a plurality of environmental factors, fusing and constructing a space-time aligned watershed ecosystem data set, and preprocessing the data set; constructing a prediction model and performing training; environmental factors are collected from a target drainage basin and input into the trained prediction model, and a nonlinear mapping relation between the environmental factors and the microbial community, namely an SHAP value, is obtained through learning; constructing a GAM response curve based on the SHAP values and the environmental factors, fitting a nonlinear response trend of the GAM response curve through a spline function, and identifying an environmental threshold corresponding to each environmental factor; through a Bootstrap resampling method, multiple times of environment threshold recognition are carried out on each environment factor, and a confidence interval is estimated based on environment thresholds obtained through multiple times of recognition, so that the stability and reliability of environment threshold recognition are improved.

Description

Technical Field

[0001] The present invention belongs to the fields of environmental ecology and computational biology, and relates to a method for identifying key environmental thresholds affecting the abundance of planktonic bacterial communities in river ecosystems by combining a neural network autoencoder and a Generalized Additive Model (GAM). Background Art

[0002] Ecological thresholds refer to the limits and critical points of the carrying capacity of an ecosystem, and are the inflection points at which an ecosystem transforms from one stable state to another. Related research on ecological thresholds has been widely carried out in ecosystems such as wetlands, forests, grasslands, and oceans. The thresholds involved include: vegetation degradation thresholds, species diversity protection thresholds, soil and water conservation thresholds, and soil safety thresholds, etc. The quantitative research on ecological thresholds can reflect the carrying limit values of an ecosystem for driving factors in terms of time or space. Therefore, based on the ecological threshold theory, managing the ecosystem and evaluating the carrying capacity of resources and the environment can scientifically and effectively manage the ecosystem and conduct safety monitoring and early warning, so as to restore or repair it in time before irreversible changes occur in the ecosystem.

[0003] With the intensification of global climate change and the increase in the intensity of human activities, the microbial community structure in freshwater ecosystems has changed significantly. Exploring its response mechanism to environmental stress and key ecological thresholds has become an important topic in environmental ecology research. Traditional ecological threshold identification methods mostly rely on macroscopic ecological indicators such as plant productivity and species diversity. These indicators have problems such as strong response lag, mismatched time scales, and low spatial resolution in responding to changes in microbial communities, and it is difficult to meet the monitoring and management requirements for the rapid dynamic changes of microorganisms. In current ecological management practices, although water quality parameters and biodiversity indices are widely used for monitoring and evaluation, they have significant deficiencies in revealing the non-linear responses and multi-factor interaction mechanisms of microbial communities. Due to their complex functional redundancy and metabolic network characteristics, microbial communities exhibit highly non-linear and high-dimensional response patterns when facing environmental changes. This characteristic poses higher requirements for threshold identification methods, that is, not only good prediction performance but also the ecological interpretable identification of key variables and their thresholds should be achieved.

[0004] Although machine learning techniques (such as random forests and deep neural networks) have been widely used in microbiome data analysis in recent years and have achieved breakthroughs in prediction accuracy, due to their "black box" nature, the application of these models in ecological mechanism analysis and policy support is limited. Especially in environmental threshold identification, traditional machine learning methods often struggle to provide ecological key explanations such as variable importance ranking and turning point location. In addition, with the rapid accumulation of high-throughput sequencing technology and multi-omics data, researchers can quantify the complex relationships between microorganisms and the environment at a more detailed level. However, there are still challenges in integrating multi-source heterogeneous data and extracting mechanistic knowledge in existing analysis frameworks. On the one hand, the non-linear coupling relationship between multidimensional environmental factors and microbial abundance is difficult to accurately capture through traditional linear models; on the other hand, there is a lack of analytical tools that combine good interpretability and non-linear modeling capabilities.

[0005] In existing research, as an unsupervised learning method, autoencoders have been used to construct a low-dimensional latent space of microbial communities, effectively improving the generalization ability of prediction. However, current research based on autoencoders mostly focuses on community structure prediction and has not effectively integrated it with interpretable models to achieve the identification of environmental thresholds and mechanism analysis.

[0006] In summary, the current technology mainly has the following deficiencies: 1. Lack of an environmental threshold identification method that can simultaneously consider non-linear modeling capabilities and ecological interpretability; 2. Although neural networks can model complex relationships, they lack a quantitative evaluation mechanism for feature contribution degrees and are difficult to reveal key regulatory factors; 3. The potential threshold response relationship between environmental factors and microbial communities has not been systematically identified and verified. Summary of the Invention

[0007] In view of this, the purpose of the present invention is to provide a method for identifying environmental thresholds of planktonic bacteria in river ecosystems, which integrates the feature extraction advantages of neural networks and the interpretability of generalized additive models, and interprets the model results through the SHAP (Shapley Additive Explanations) algorithm, so as to achieve the accurate identification of environmental response thresholds of microbial communities and ecological mechanism analysis, and provide scientific support for ecological early warning and precise management of river ecosystems.

[0008] To achieve the above object, the present invention provides the following technical solutions:

[0009] A method for identifying environmental thresholds of planktonic bacteria in a river ecosystem, the method comprising:

[0010] S1. Collect environmental factors including satellite remote sensing inversion data, physical and chemical parameters of hydrological monitoring stations, and OTU abundances of planktonic bacteria, and fuse them to construct a spatio-temporal aligned dataset of the watershed ecosystem. After normalizing and applying centered log-ratio transformation to the dataset, divide it into a training set and a test set;

[0011] S2. Construct a prediction model, train the prediction model using the test set, and verify the prediction accuracy of the prediction model using the test set;

[0012] S3. Collect environmental factors from the target watershed and input them into the trained prediction model to learn the non-linear mapping relationship between environmental factors and microbial communities, i.e., SHAP values;

[0013] S4. Construct a GAM response curve based on SHAP values and environmental factors, and fit the non-linear response trend of the GAM response curve through a spline function to identify the environmental thresholds corresponding to each environmental factor;

[0014] S5. Use the Bootstrap resampling method to identify the environmental thresholds of each environmental factor multiple times, and estimate the confidence interval based on the environmental thresholds obtained from multiple identifications to improve the stability and reliability of environmental threshold identification.

[0015] Further, in step S1, the centered log-ratio transformation of the OTU abundance data of planktonic bacteria includes: for the data collected from the target watershed, replace all OTU abundances that are zero with a constant that is much smaller than the minimum non-zero OTU abundance; after replacement, for each OTU abundance, first take its natural logarithm, and then subtract the natural logarithm of the geometric mean of all OTU abundances, that is, the centered log-ratio transformation result of each OTU abundance is obtained.

[0016] Further, in step S2, the constructed prediction model includes an encoder, a decoder, and an interpretability analysis module. Among them, the structures of the encoder and the decoder are symmetric; the encoder is used to extract the potential ecological gradient features of the input data, and the decoder is used to reconstruct the predicted abundance distribution of the microbial community; the interpretability analysis module is used to calculate the SHAP values corresponding to each input environmental factor, quantify the contribution of each environmental factor to the prediction result, so as to enhance the interpretability of the prediction model.

[0017] Among them, the encoder includes three fully connected hidden layers, and a tanh activation function is connected after each hidden layer; among them, along the direction of data transmission, the number of neurons in each hidden layer decreases in turn. The decoder is symmetric with the encoder. In the decoder, along the direction of data transmission, the number of neurons in each hidden layer increases in turn.

[0018] The interpretability analysis module calculates the SHAP value through the following formula:

[0019]

[0020] In the formula, F represents the complete set of environmental factors in the input prediction model environment, S represents the subset excluding environmental factor i, f(S) represents the prediction output of the prediction model under subset S, and φ i represents the SHAP value of environmental factor i, that is, the average marginal contribution of environmental factor i to the model output.

[0021] Furthermore, in step S4, constructing the GAM response curve based on the SHAP value and environmental factors includes, for each environmental factor, taking the value of the environmental factor as the horizontal axis and the SHAP values corresponding to different values of the environmental factor as the vertical axis to construct the GAM response curve;

[0022] Identifying the environmental threshold by fitting the non-linear response trend of the GAM response curve with a spline function includes using the spline function to fit the non-linear trend of the GAM response curve, and then performing numerical derivative calculation on the fitted curve to identify the turning point where the SHAP value changes from negative to positive or from positive to negative on the fitted curve, that is, the point where the SHAP value is 0, and this point is the environmental threshold corresponding to the environmental factor.

[0023] Furthermore, step S5 is specifically as follows: for each environmental factor i, perform n times of resampling with replacement from its values and the corresponding SHAP values, construct the GAM response curve based on the samples obtained each time, fit the non-linear response trend of the GAM response curve with a spline function, then perform numerical derivative calculation on the fitted curve to identify the point where the SHAP value is 0 on the fitted curve, that is, the environmental threshold; obtain n independent environmental thresholds of environmental factor i through n times of resampling, and calculate the confidence interval with a 95% confidence level for these n environmental thresholds.

[0024] Among them, if there are multiple points where the SHAP value is 0 on the fitted curve of a certain environmental factor, select the turning point with the largest absolute value of the derivative in the fitted curve as the main environmental threshold.

[0025] The beneficial effects of the present invention are as follows: By constructing a prediction model integrating an encoder-decoder and an interpretability analysis module, the present invention breaks through the bottleneck of insufficient interpretability of traditional machine learning models. The encoder-decoder architecture effectively extracts the potential ecological gradient features of microbial communities and environmental factors. Combining with the quantitative evaluation mechanism of SHAP values, it not only realizes the non-linear mapping modeling of high-dimensional data, but also quantifies the contributions of various environmental factors to the prediction results, enhancing the interpretability of the constructed model. In addition, the present invention establishes a rapid response mechanism between the dynamics of microbial communities and environmental factors, solving the problem of strong lag of traditional macroscopic indicators. Through the contribution analysis of SHAP values, the key environmental regulation factors and their interaction networks can be accurately located, providing a scientific basis for the targeted governance of watershed ecosystems.

[0026] In summary, the environmental threshold identification method provided by the present invention breaks through the dependence of traditional statistical methods on linear relationships and single-factor responses, and has the ability to handle complex ecological relationships such as high-dimensional, non-linear, and multi-factor interactions. It can automatically identify key environmental factors and their response thresholds without relying on expert experience rules, providing a theoretical basis and technical support for the dynamic monitoring and precise management of microbial communities in river ecosystems, and having broad application prospects.

[0027] Other advantages, objectives, and features of the present invention will be described to some extent in the subsequent specification, and to some extent, will be obvious to those skilled in the art based on the study of the following text, or can be taught from the practice of the present invention. The objectives and other advantages of the present invention can be realized and obtained through the following specification. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] In order to make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be described in detail with reference to the accompanying drawings, where:

[0029] Figure 1 is a schematic flow chart of a method for identifying environmental thresholds of planktonic bacteria in a river ecosystem provided by an embodiment of the present invention;

[0030] Figure 2 is a schematic diagram of the training results of a machine learning model in an embodiment of the present invention;

[0031] Figure 3 is a scatter plot of temperature-SHAP values of Actinobacteria;

[0032] Figure 4 is a schematic diagram of the GAM fitting results of Actinobacteria. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0033] The following describes the embodiments of the present invention through specific examples. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments. Various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the drawings provided in the following embodiments only illustrate the basic concept of the present invention in a schematic manner. Without conflict, the following embodiments and the features in the embodiments can be combined with each other.

[0034] Among them, the drawings are only for illustrative purposes, showing only schematic diagrams rather than actual diagrams, and should not be construed as a limitation on the present invention; in order to better illustrate the embodiments of the present invention, some components in the drawings will be omitted, enlarged or reduced, which does not represent the size of the actual product; for those skilled in the art, it is understandable that some well-known structures and their descriptions in the drawings may be omitted.

[0035] In the drawings of the embodiments of the present invention, the same or similar reference numerals correspond to the same or similar components; in the description of the present invention, it should be understood that if there are terms such as "upper", "lower", "left", "right", "front", "rear", etc. indicating the orientation or positional relationship, they are based on the orientation or positional relationship shown in the drawings. This is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation. Therefore, the terms describing the positional relationship in the drawings are only for illustrative purposes and should not be construed as a limitation on the present invention. For those of ordinary skill in the art, the specific meanings of the above terms can be understood according to specific circumstances.

[0036] The present invention proposes a method for identifying environmental thresholds of planktonic bacteria in river ecosystems based on the combination of neural network autoencoders and generalized additive models, which closely combines SHAP values with GAM. By integrating multi-source heterogeneous data and constructing an interpretable machine learning model, it dynamically identifies the key environmental thresholds that affect the changes in the planktonic bacterial community structure in river ecosystems and reveals its potential ecological regulation mechanism. GAM is not only used to visualize the response trend between environmental factors and SHAP values, but also used to model and explain their non-linear relationship at the statistical modeling level. This method takes the value of each environmental factor as the independent variable and the corresponding SHAP value as the dependent variable, and constructs a GAM model based on spline functions to achieve smooth fitting of the local change trend. Through the GAM model, the "non-linear mutation points" in the response curves of planktonic bacteria to each environmental factor can be identified, and these mutation points reflect the sensitive interval of the ecosystem to external disturbances, that is, the key threshold area.

[0037] Based on this, an embodiment of the present invention provides a method for identifying the environmental thresholds of planktonic bacteria in river ecosystems, which mainly includes three steps: ecological big data processing, interpretable machine learning modeling, and dynamic identification of environmental thresholds.

[0038] Step 1: Data collection and processing

[0039] Integrate multi-source heterogeneous data, including satellite remote sensing inversion data (such as chlorophyll a concentration distribution), physical and chemical parameters (such as water temperature WT, dissolved oxygen DO, pH value) collected in real time by hydrological monitoring stations, and metagenomics data obtained by high-throughput sequencing (such as 16S rRNA amplicon sequences), to construct a spatio-temporally aligned watershed ecosystem database. Specifically, in this embodiment, data from 417 sampling points in the upper reaches of the Yangtze River (from Batang, Sichuan to Yichang, Hubei) between 2015 and 2023 are used, covering the OTU abundance information of planktonic bacteria and related hydrological and water environment parameters. Environmental factors cover climate conditions (air temperature AT, precipitation Pre) and water body physical and chemical characteristics (total nitrogen TN, total phosphorus TP, dissolved organic carbon DOC, chlorophyll Chla, conductivity Cond, pH, dissolved oxygen DO, water temperature WT).

[0040] In terms of data preprocessing, first remove missing values and outliers, and then perform Min-Max normalization on environmental variables to unify the dimension range. The specific formula is:

[0041]

[0042] For planktonic bacteria abundance data, use centered log-Ratio (CLR) transformation to eliminate the closure effect problem in microbial community composition data, and correct the uncertainty brought by zero values by adding a small constant. The specific formula is:

[0043]

[0044]

[0045] where x i is the relative abundance of the i-th OTU in the sample, and g(x) is the geometric mean of the abundances of all OTUs.

[0046] In this process, first use the geometric mean method to standardize zero values, that is, in each sample, uniformly replace all zero OTU abundances with a very small constant (such as 10 -6) Thus, mathematical errors during logarithmic transformation are avoided. Subsequently, the natural logarithm of all OTU abundances in each sample is taken, and the natural logarithm of the geometric mean of all OTU abundances in that sample is subtracted, thereby obtaining the CLR transformation result, as shown in Equation (2). Through centered log-ratio transformation, the error caused by component relativity is effectively solved, improving the accuracy and stability of subsequent modeling.

[0047] Finally, the preprocessed data is divided into a training set and a test set at a ratio of 7:3 for subsequent model training. The original environmental dataset is shown in Table 1, and the relative abundances of bacteria at the phylum level are shown in Table 2.

[0048] Table 1

[0049]

[0050] Table 2

[0051]

[0052] Step 2: Interpretable machine learning modeling and SHAP value calculation

[0053] An autoencoder structure is adopted to construct a prediction model. The encoder embeds environmental factors into the latent space to extract ecological gradients, and the decoder realizes the non-linear prediction of the composition of the planktonic bacterial community. To achieve model interpretability, the KernelExplainer method in the SHAP library is further used to explain the entire process of the input environmental factors. KernelExplainer quantifies the contribution of each environmental factor to the prediction result in a specific sample by simulating the influence of different feature combinations on the model output, and outputs the corresponding SHAP values, laying a foundation for threshold identification.

[0054] Among them, the constructed prediction model is a hybrid architecture model based on an artificial neural network autoencoder, which is used to learn the non-linear mapping relationship between environmental factors and microbial communities. The structure of this model includes an encoder, a decoder, and an interpretability analysis module. The encoder part consists of three fully connected hidden layers with the number of neurons decreasing successively (such as 128 - 64 - 16), and a tanh activation function is connected after each layer to compress the environmental factor information layer by layer and extract potential ecological gradient features; the decoder part is a symmetric structure of the encoder (16 - 64 - 128) for reconstructing the predicted abundance distribution of the microbial community. This model takes various environmental factors (such as air temperature AT, precipitation Pre, dissolved organic carbon DOC, total phosphorus TP, etc.) as inputs, extracts potential ecological gradient features through the encoder module composed of a multi-layer perceptron, and reconstructs the abundance distribution of the microbial community through the decoder to effectively predict the microbial structure. The prediction results are as Figure 2 shown.

[0055] A SHAP value close to zero generally indicates that the factor makes a relatively small contribution to the current prediction. It may be in the critical zone of ecosystem regulation and has the ability to identify potential ecological thresholds. The calculation formula of the SHAP value is as follows:

[0056]

[0057] where F is the complete set of input environmental factors, S is the subset without environmental factor i, f(S) is the predicted output of the model under subset S, and φ i is the SHAP value of environmental factor i, that is, the average marginal contribution of environmental factor i to the model output.

[0058] Step 3: GAM fitting and environmental threshold identification

[0059] To identify the key regulation intervals of environmental factors on microbial abundance, the method described in this embodiment constructs a generalized additive model (GAM) based on the relationship curve between the SHAP value and environmental factors, captures the non-linear response trend through spline functions, and systematically analyzes the non-linear relationship between environmental factors and microbial community responses.

[0060] Specifically, for each environmental factor, with the value of the environmental factor as the horizontal axis and the SHAP value corresponding to different values of the environmental factor as the vertical axis, a GAM response curve is constructed, and its non-linear trend is fitted using spline functions. Subsequently, by calculating the numerical derivative of the fitted curve, the turning points where the SHAP value changes from negative to positive or from positive to negative in the fitted curve are identified, thereby defining the key inflection points where the microbial community responds to specific environmental factors, that is, the critical threshold points of ecological balance (where the SHAP value is 0).

[0061] To improve the robustness of threshold identification, the Bootstrap resampling method is further introduced: based on the response data of each environmental factor, n = 1000 resamplings with replacement are performed. In each sampled sample, GAM fitting and derivative analysis are carried out to obtain the corresponding critical threshold points, and finally 1000 independent threshold estimates are obtained. The distribution of the main thresholds in all resampling results is statistically analyzed, and the confidence interval (Confidence Interval, CI) with a confidence level of 95% is calculated. In the case where there are multiple SHAP zero-crossing points for an environmental factor, the turning point with the largest absolute value of the GAM curve derivative is preferentially selected as the main threshold to ensure that the identified critical points have significant ecological response significance and regulatory value.

[0062] Taking Actinobacteria as an example in this embodiment, its environmental threshold is identified, and the results are shown in Table 3:

[0063] Table 3

[0064]

[0065] As can be seen from Table 3, there are clear "ecological response thresholds" in the regulation of the planktonic bacterial community structure by different environmental factors. For example, when the water temperature (WT) exceeds 15.808 °C, the abundance of Actinobacteria changes significantly, indicating that this temperature value may be a key node for the ecosystem to shift from one stable state to another. Similarly, total nitrogen (TN) and total phosphorus (TP) show threshold effects near 0.738 mg / L and 0.058 mg / L respectively, suggesting that small fluctuations in the concentrations of nitrogen and phosphorus nutrients may also trigger changes in the community structure.

[0066] The scatter plot of the SHAP values of the temperature of Actinobacteria is as Figure 3 shown. The GAM fitting curve of Actinobacteria is as Figure 4 shown. On the fitting curve, the red dashed line indicates the position where the SHAP value is close to 0, the blue area indicates that the environmental factor has a negative impact on the abundance of planktonic bacteria, and the yellow area indicates that the environmental factor has a positive impact on the abundance of planktonic bacteria. It should be noted that for some environmental factors (such as TN, DOC, WT, pH), the GAM fitting results show multiple intersections with a SHAP value of zero, indicating that their ecological responses may have multiple sensitive intervals.

[0067] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that the technical solutions of the present invention can be modified or equivalently replaced without departing from the spirit and scope of the present technical solution, and they should all be covered by the scope of the claims of the present invention.

Claims

1. A method for identifying the environmental threshold of planktonic bacteria in a river ecosystem, characterized in that The method includes: S1. Collect environmental factors including satellite remote sensing inversion data, physical and chemical parameters of hydrological monitoring stations, and OTU abundances of planktonic bacteria, and fuse them to construct a spatio-temporally aligned dataset of the watershed ecosystem. After normalizing and performing centered log-ratio transformation on the dataset, it is divided into a training set and a test set; S2. Construct a prediction model, train the prediction model using the test set, and verify the prediction accuracy of the prediction model using the test set; S3. Collect environmental factors from the target watershed and input them into the trained prediction model to learn the non-linear mapping relationship between environmental factors and microbial communities, i.e., SHAP values; S4. Construct a GAM response curve based on SHAP values and environmental factors, and fit the non-linear response trend of the GAM response curve through a spline function to identify the environmental thresholds corresponding to each environmental factor; S5. Use the Bootstrap resampling method to identify the environmental thresholds of each environmental factor multiple times, and estimate the confidence interval based on the environmental thresholds obtained from multiple identifications to improve the stability and reliability of environmental threshold identification.

2. The method according to claim 1, wherein In step S1, the centered log-ratio transformation of the OTU abundance data of planktonic bacteria includes: for the data collected from the target watershed, replace all OTU abundances that are zero with a constant much smaller than the minimum non-zero OTU abundance; after replacement, for each OTU abundance, first take its natural logarithm, and then subtract the natural logarithm of the geometric mean of all OTU abundances, that is, the centered log-ratio transformation result of each OTU abundance is obtained.

3. The method according to claim 1, characterized in that, In step S2, the constructed prediction model includes an encoder, a decoder, and an interpretability analysis module. Among them, the structures of the encoder and the decoder are symmetric; the encoder is used to extract the latent ecological gradient features of the input data, and the decoder is used to reconstruct the predicted abundance distribution of the microbial community; the interpretability analysis module is used to calculate the SHAP values corresponding to each input environmental factor, quantify the contribution of each environmental factor to the prediction result, so as to enhance the interpretability of the prediction model.

4. The method according to claim 3, wherein The encoder includes three fully connected hidden layers, and a tanh activation function is connected after each hidden layer; among them, along the direction of data transmission, the number of neurons in each hidden layer decreases sequentially; The decoder is symmetric with the encoder structure. In the decoder, along the direction of data transmission, the number of neurons in each hidden layer increases sequentially.

5. The method according to claim 3, wherein The SHAP value is calculated in the interpretability analysis module by the following formula: Wherein, F represents the complete set of input prediction model environmental factors, S represents the subset excluding environmental factor i, f(S) represents the prediction output of the prediction model under subset S, and φ i represents the SHAP value of environmental factor i, that is, the average marginal contribution of environmental factor i to the model output.

6. The method according to claim 1, wherein In step S4, the construction of the GAM response curve based on SHAP values and environmental factors includes, for each environmental factor, taking the value of the environmental factor as the horizontal axis and the SHAP values corresponding to different values of the environmental factor as the vertical axis to construct a GAM response curve; The non-linear response trend of the GAM response curve is fitted by a spline function, and the identification of the environmental threshold includes using the spline function to fit the non-linear trend of the GAM response curve, and then calculating the numerical derivative of the fitted curve to identify the turning point where the SHAP value changes from negative to positive or from positive to negative on the fitted curve, that is, the point where the SHAP value is 0, and this point is the environmental threshold corresponding to the environmental factor.

7. The method according to claim 1, characterized in that, Step S5 includes, for each environmental factor i, performing n times of resampling with replacement from its values and the corresponding SHAP values, constructing a GAM response curve based on the samples obtained from each sampling, fitting the non-linear response trend of the GAM response curve by a spline function, then calculating the numerical derivative of the fitted curve, and identifying the point where the SHAP value is 0 on the fitted curve, that is, the environmental threshold; obtaining n independent environmental thresholds of the environmental factor i through n times of resampling, and calculating the confidence interval with 95% confidence level of these n environmental thresholds.

8. The method according to claim 6 or 7, characterized in that If there are multiple points where the SHAP value is 0 on the fitted curve of a certain environmental factor, then select the turning point with the largest absolute value of the derivative in the fitted curve as the main environmental threshold.

Images