Non-Linear Visualization and Importance Ratio Analysis of Multivariate Polynomial Regression Ecological Models Based on River Hydromorphology and Water Quality

: Multivariate polynomial regression (MPR) models were developed for ﬁve macrophyte indices. MPR models are able to capture complex interactions in the data while being tractable and transparent for further analysis. The performance of the MPR modeling approach was compared to previous work using artiﬁcial neural networks. The data were obtained from hydromorpholog-ically modiﬁed Polish rivers with a widely varying water quality. The modeled indices were the Macrophyte Index for Rivers (MIR), the Macrophyte Biological Index for Rivers (IBMR), and the River Macrophyte Nutrient Index (RMNI). These indices measure the trophic and ecological status of the rivers. Additionally, two biological diversity indices, species richness (N) and the Simpson index (D), were modeled. The explanatory variables were physico-chemical properties depicting water quality and river hydromorphological status indices. In comparison to artiﬁcial neural networks, the MPR models performed similarly in terms of goodness of ﬁt. However, the MPR models had advantages such as model simplicity and ability to be subject to effective visualization of complex nonlinear input–output relationships, as well as facilitating sensitivity analysis using importance ratios to identify effects of individual input variables.


Introduction
Several studies acknowledge the significance of accurately accounting for multiple sources of ambiguity when modeling ecological phenomena [1]. Additionally, the assessment of large waterbodies is a challenging task [2]. Many surveys report that the reaction of aquatic communities to varying environmental conditions is often responsible for river degradation, mainly identifying eutrophication [3][4][5][6]. Several studies also focus on the quality of groundwater and of water bodies such as rivers, watersheds, and coastal environments. This plays an important role in determining its impact on public health and the environment [7,8]. Many countries in Europe have started developing river monitoring systems in their national monitoring programs based on macrophytes [9,10], According to the Water Framework Directive (WFD) [11], assessment of freshwater is based on ecological status consisting of biological indicators (fish, macroinvertebrates, phytoplankton, phytobenthos, and macrophytes), supported by water quality and physical conditions of ecosystems. The significance of macrophytes in biological river assessment is formally acknowledged under the WFD, and macrophytes are an essential element in the monitoring of ecological status and surface water quality. Macrophyte development is based on several abiotic and biotic factors consisting of nutrient concentrations, flow velocity, hydraulic using MPR to better understand and visualize the relationships between the predictors and the dependent variables, and to quantify the importance of each individual predictor on the dependent variables. This would be very cumbersome to do with ANN models. Furthermore, the MPR models themselves can be easily presented in this work for use by other investigators. The objectives of this study are to (i) develop MPR models to predict macrophyte indices in Polish rivers; (ii) perform importance ratio analysis for individual predictor variables; (iii) to visualize the non-linear relationships between the dependent and predictor variables graphically; and (iv) to compare the goodness of fit statistics and model interpretations with earlier ANN models based on the same data.

Data Management
The data were obtained from over 200 survey sites in Poland representing hydromorphologically altered rivers with a wide range of water quality conditions and comprised of 14 explanatory variables. The macrophyte indices of ecological status assessment modeled in this study were Macrophyte Index for Rivers (MIR) [30], River Macrophyte Nutrient Index (RMNI) [31], and Macrophyte Biological Index for Rivers (IBMR) [9]. Species Richness Index (N) and the Simpson Index (D) [32], which measure biological diversity were also modeled. The data was previously modeled using Artificial Neural Networks (ANNs) [18]. The three indices of the ecological status express the trophic degradation of river ecosystems. These indices are based on the coverage of macrophytes using five-grade (IBMR) or nine-grade (MIR, RMNI) scale and their indicator values. These values define the average level of water trophy of each species and the tolerance to this factor. The Simpson index (D) is one of the most often used alpha diversity indices, which includes the number of species present at the site, as well as the relative abundance of species. The species richness index (N) represents the number of species that occur at the site. The basic equations used to calculate each of these macrophyte indices are shown in the Supporting Information (Table S1). Table 1 shows the fourteen candidate predictor variables and the five indices (MIR, RMNI, IBMR, D, and N) to be predicted. The predictor variables were chosen such that they adequately represented aquatic vegetation, hydromorphological assessment, and physicochemical properties of the water body. Hydromorphological assessment was carried out following a river habitat survey (RHS) [33]. This method categorizes the physical state of a river by the calculation of certain indices. The HMS is an index extracted from the data concerning morphological modification of the river channel due to human activities (e.g., bank reinforcement, channel re-sectioning, culverting, number of weirs, etc.). The HQA index estimates the modification quality of the river through the diversity of features evaluated (e.g., number of different flow types, different substrates, and naturalness of land use). All the available fourteen independent variables were tested for predicting the five macrophyte indices.

Model Development
The dataset was divided into three parts randomly when the data were modeled using ANN [18]. To make comparisons meaningful, the same fit, cross, and validation datasets were retained for MPR modeling. The fit dataset (n = 140) was used for determining the regression coefficients. The cross dataset (n = 30) was used to determine global goodness-offit for cross-validation. This means the coefficients were computed using the fit dataset, but only those terms were added to the model that improved the prediction of the cross dataset and having a p-value < 5%. The validation dataset (n = 30) was used to generate the same global goodness-of-fit statistics for final validation after the models were generated i.e., after the best model is finalized based on fit and cross datasets. This approach minimizes the possibility of overfitting and helps to ensure the resulting model is generalizable, i.e., applicable to new data.
We produced 5 regression models for each index depending on 14 predictors, such as The exponents could be either integers or fractions. In this study, only integer exponents between −3 and +3 were considered. Preliminary examination of fractions and integers outside this range did not result in significantly better models. Given this and the Occam's Razor principle, higher absolute exponents, multipliers, and fractional exponents were excluded from detailed consideration. TaylorFit software gives a conventional mathematical description of MPR in the form of a general equation (Equation (1)). Basically, it is an extension of multilinear regression (MLR) that includes interactions and other nonlinearities. Based on Taylor Theorem, it describes a polynomial series equation with a finite number of terms for each of the index to yield a definite degree of accuracy.
Further, the number of terms increases with the number of exponents and the number of multiplicands. The maximum number of multiplicands for each candidate term used in this study was 4. A stepwise algorithm, as follows, was used to build the model: The model always started with an intercept (average of the dependent variable values). The software generated terms that best combined with existing terms of a model, depending on the allowable exponents and multiplicands set by the user. The terms were sorted in terms of best t-statistics of the fit data.

2.
A candidate term was chosen from among the statistically significant variables of fit data, which also improved the R 2 of a separate cross-correlation dataset. Both these criteria had to be met for a term to be added into the model. This procedure reduced the possibility of overfitting and improved the generalizability of the model.

3.
After any term was added to the model, the other terms previously added were tested for statistical significance and removed if not. 4.
The above process was repeated iteratively for additional candidate terms. 5.
The model was thus built by an iterative process by adding and removing candidate terms from among statistically significant terms based on fit dataset that also improved the R 2 of the test dataset, until the model could not be improved by addition or removal of any single term. 6.
After the model was complete, it was tested against a 3rd independent validation dataset. The performance with the validation dataset was used for comparison with the ANN models.

Sensitivity Analysis
Sensitivity analysis was conducted based on select variables that existed in each of the model indices using the MPR. The local sensitivity index, δ, is given by [24,34] where f is the model prediction of the dependent variable (MIR, RMNI, IBMR, D, and N, respectively, in our case), and X is the independent variable whose effect on the dependent variable is being examined (e.g., HMS, HQA, P tot , etc.). However, local sensitivity indices of different variables cannot be compared with each other because of the varying units and ranges of different predictor variables, and therefore, they cannot be used to obtain their relative importance on the output. Consequently, a more useful indicator of sensitivity considers the spread in the variable involved and normalizes the units. The importance ratio [34] as defined in Equation (3), accomplishes both goals.
where s X and s f are the standard deviations of the independent and dependent variables in the fit dataset, respectively. The interpretation of the importance ratio is that it indicates the sensitivity relative to the spread of the variables. For example, if the importance ratio for y versus x is 0.30, this means that if we increased x by, say, 1% of its standard deviation, then y would increase by about 0.3% of its standard deviation. The importance ratio (and the sensitivity) for nonlinear models is not a single value, as for linear models. It will be different for each combination of independent variables. Thus, it is better to examine it as a distribution. Thus, we reported tenth, fiftieth, and ninetieth percentiles as well as a graphical display of the cumulative distribution. For a single central tendency, an appropriate measure would be the root-mean-square (RMS) value. We used the RMS importance ratio as a singular way to compare the impact of different independent variables on the prediction. These are examples of ways that representative models such as MPR can provide information that is difficult to obtain with ANN models.
MPR modeling using TaylorFit can produce models in the same form as the evolutionary polynomial regression (EPR) approach [35]. The difference is in the type of optimization approach used to select model terms, resulting in different constraints on the model. Plus, TaylorFit has built-in tools for sensitivity and importance ratio analysis as described above.

Results
The MPR models obtained for five indices are shown in Table 2. The t-statistics and associated probability that the coefficient is different from zero (p(t)) are shown in Table 3. Macrophyte communities from different river systems make comparisons between different types of rivers difficult [36]. Each of the fitted models cannot be applied as-is to any other river because these locations will have different physicochemical properties and different hydromorphological status. However, a similar approach can work in other geomorphological regions and other river types by generating new models following the procedure presented in this paper.  Table 1).

MIR
3.34780 The MPR model for MIR had validation R 2 and MSE values of 0.580 and 10.86, respectively ( Figure S1, Table 4). In comparison, the corresponding values of the ANN model were 0.702 and 9.790, respectively. Figure 1a shows the nonlinear response of MIR with respect to HQA and conductivity. For the nonlinear plots, all other variables in the model, except those represented on the plot, were held at their average values.  Compared to ANN, the use of MPR enabled a much better understanding of the predictors effect on different macrophyte indices. Sensitivity analysis, or other explanatory analysis that are utilized in ANN, showed relations between variables [37], but it did not give precise information about existing relationships. Previous studies [18] showed the general importance of parameters on MIR and other macrophyte indices, but MPR revealed a precise interaction between the modeled indices and individual predictor variables. For the first index, as conductivity increased, MIR decreased at a decreasing rate. As HQA increased, MIR increased as well. These results are like other studies showing the importance of conductivity and hydromorphological conditions to macrophytes [38]. The cumulative frequency distributions of the importance ratios for the dependent variables for MIR (Figure 1b) showed that conductivity had a negative importance ratio for all observed values, indicating its negative correlation with MIR. The importance ratios were calculated using Equation (3) for all combinations of observed values, and a cumulative distribution was built by rank ordering all calculated importance ratios. For MIR, the variable with the most impact on its variability was conductivity, as evident by its relatively high root mean square importance ratio (1.358 from Table 5). HQA had mostly positive importance ratios, indicating its positive correlation with MIR for all observed conditions. pH, HMS, and BOD had relatively lower importance ratios than HQA and conductivity. It is evident that the macrophyte increased with the HQA index depicting a broad measure of diversity and "naturalness" of the river. The HQA score is determined by the presence and extent of habitat features of the known wildlife examined in the survey. On the other hand, HMS index, which assesses the human modification activities around the river increased while the MIR index decreased.  Compared to ANN, the use of MPR enabled a much better understanding of the predictors effect on different macrophyte indices. Sensitivity analysis, or other explanatory analysis that are utilized in ANN, showed relations between variables [37], but it did not give precise information about existing relationships. Previous studies [18] showed the general importance of parameters on MIR and other macrophyte indices, but MPR revealed a precise interaction between the modeled indices and individual predictor variables. For the first index, as conductivity increased, MIR decreased at a decreasing rate. As HQA increased, MIR increased as well. These results are like other studies showing the importance of conductivity and hydromorphological conditions to macrophytes [38]. The cumulative frequency distributions of the importance ratios for the dependent variables for MIR (Figure 1b) showed that conductivity had a negative importance ratio for all observed values, indicating its negative correlation with MIR. The importance ratios were calculated using Equation (3) for all combinations of observed values, and a cumulative distribution was built by rank ordering all calculated importance ratios. For MIR, the variable with the most impact on its variability was conductivity, as evident by its relatively high root mean square importance ratio (1.358 from Table 5). HQA had mostly positive importance ratios, indicating its positive correlation with MIR for all observed conditions. pH, HMS, and BOD had relatively lower importance ratios than HQA and conductivity. It is evident that the macrophyte increased with the HQA index depicting a broad measure of diversity and "naturalness" of the river. The HQA score is determined by the presence and extent of habitat features of the known wildlife examined in the survey. On the other hand, HMS index, which assesses the human modification activities around the river increased while the MIR index decreased. The RMNI model showed validation R 2 and MSE values of 0.650 and 0.050, respectively ( Figure S2, Table 4). Consequently, the ANN model had the values of 0.715 and 0.050, respectively. RMNI is extremely sensitive to nutrient composition, as is evident in the nonlinear plot ( Figure 2a). As the total phosphorus and nitrogen increase, the RMNI index increase gradually [39]. As seen, the increase in the index was gradual until it attained a constant trend. Figure 2b depicts the importance ratio plots consisting of positive implications on RMNI, except nitrate nitrogen, which had a negative impact. In contrast to MIR, RMNI did not show any relationship to HMS and did not seem to rely significantly on HQA either. These findings contrast with the results of modeling using neural networks but were consistent with the assumptions of the index [31]. For RMNI, the variable with the most influence on its variance was ammonia nitrogen, as perceptible by its relatively high root mean square importance ratio (14.640 from Table 5).  Table 4), compared to the ANN model, which had an R 2 value of 0.532 and an MSE of 0.410. Figure 3a shows the nonlinear response of IBMR with respect to HMS and organic nitrogen. An increase in HMS and organic nitrogen had a negative effect on IBMR because of implications of human activities attributed to the HMS index. The clear relation with these two environmental variables is reliable with the IBMR assumption [9], but so far rarely indicated in other studies. The importance ratio of independent variables (Figure 3b) showed negative impacts in terms of importance ratio for IBMR, as we observed in MIR. Nitrite nitrogen HMS seemed to have a positive correlation with IBMR. Organic nitrogen seemed to have negative as well as positive impacts with the index. However, biochemical oxygen demand did not seem to have any significant impact of importance ratio on IBMR. The highest influence of a variable on IBMR was shown by nitrite nitrogen, which had an RMS of 14.607 (Table 5). This index has weight given to trophic level and heavy organic pollutants [40], the effect of which was effectively captured in the MPR model through organic nitrogen and nitrite nitrogen. Unlike other indices, IBMR has not been proven to have shown effects of HMS and HQA indices.    Table 4), compared to the ANN model, which had an R 2 value of 0.532 and an MSE of 0.410. Figure 3a shows the nonlinear response of IBMR with respect to HMS and organic nitrogen. An increase in HMS and organic nitrogen had a negative effect on IBMR because of implications of human activities attributed to the HMS index. The clear relation with these two environmental variables is reliable with the IBMR assumption [9], but so far rarely indicated in other studies. The importance ratio of independent variables (Figure 3b) showed negative impacts in terms of importance ratio for IBMR, as we observed in MIR. Nitrite nitrogen HMS seemed to have a positive correlation with IBMR. Organic nitrogen seemed to have negative as well as positive impacts with the index. However, biochemical oxygen demand did not seem to have any significant impact of importance ratio on IBMR. The highest influence of a variable on IBMR was shown by nitrite nitrogen, which had an RMS of 14.607 (Table 5). This index has weight given to trophic level and heavy organic pollutants [40], the effect of which was effectively captured in the MPR model through organic nitrogen and nitrite nitrogen. Unlike other indices, IBMR has not been proven to have shown effects of HMS and HQA indices.  Table 4), compared to the ANN model, which had an R 2 value of 0.532 and an MSE of 0.410. Figure 3a shows the nonlinear response of IBMR with respect to HMS and organic nitrogen. An increase in HMS and organic nitrogen had a negative effect on IBMR because of implications of human activities attributed to the HMS index. The clear relation with these two environmental variables is reliable with the IBMR assumption [9], but so far rarely indicated in other studies. The importance ratio of independent variables (Figure 3b) showed negative impacts in terms of importance ratio for IBMR, as we observed in MIR. Nitrite nitrogen HMS seemed to have a positive correlation with IBMR. Organic nitrogen seemed to have negative as well as positive impacts with the index. However, biochemical oxygen demand did not seem to have any significant impact of importance ratio on IBMR. The highest influence of a variable on IBMR was shown by nitrite nitrogen, which had an RMS of 14.607 (Table 5). This index has weight given to trophic level and heavy organic pollutants [40], the effect of which was effectively captured in the MPR model through organic nitrogen and nitrite nitrogen. Unlike other indices, IBMR has not been proven to have shown effects of HMS and HQA indices.   The importance ratio plot (Figure 4b) for D showed nitrate nitrogen having a negative effect on the index, whereas HQA had a significant positive correlation with D. Biological oxygen demand (BOD 5 ) showed the highest RMS importance ratio value of 7.324 (Table 5), depicting the dominance of the variable on D. The species richness index (N) had an R 2 value of 0.330 and an MSE value of 8.39 ( Figure S5). The ANN derived model had 0.415 for validation and MSE of 9.29. HMS, HQA, and conductivity importance ratio plots (Figure 5b) showed a slight positive correlation with N index for all the observed values. Though nitrite nitrogen seemed to have a significant positive effect, but it also seemed to have a negative correlation with N. The nonlinear plot (Figures 4a and 5a) for D and N were both dependent on HMS and HQA. These plots captured the effect of human activities on the area effectively, which was clearly shown by the increasing rate of HMS index that decreased both D as well as N. Simultaneously, with increasing HQA, which demonstrated natural habitat around the region, appeared to result in surges of values for both the indices. Modeling quality for biodiversity indices, both D and N, supported the previous attempts on the relationship between biological diversity patterns and environmental factors in various types of ecosystems [41].
Water 2021, 13, x FOR PEER REVIEW 10 of 13 The importance ratio plot (Figure 4b) for D showed nitrate nitrogen having a negative effect on the index, whereas HQA had a significant positive correlation with D. Biological oxygen demand (BOD5) showed the highest RMS importance ratio value of 7.324 (Table  5), depicting the dominance of the variable on D. The species richness index (N) had an R 2 value of 0.330 and an MSE value of 8.39 ( Figure S5). The ANN derived model had 0.415 for validation and MSE of 9.29. HMS, HQA, and conductivity importance ratio plots (Figure 5b) showed a slight positive correlation with N index for all the observed values. Though nitrite nitrogen seemed to have a significant positive effect, but it also seemed to have a negative correlation with N. The nonlinear plot (Figures 4a and 5a) for D and N were both dependent on HMS and HQA. These plots captured the effect of human activities on the area effectively, which was clearly shown by the increasing rate of HMS index that decreased both D as well as N. Simultaneously, with increasing HQA, which demonstrated natural habitat around the region, appeared to result in surges of values for both the indices. Modeling quality for biodiversity indices, both D and N, supported the previous attempts on the relationship between biological diversity patterns and environmental factors in various types of ecosystems [41].  As can be seen, data can be modeled with different approaches, each producing its own performance results; but how can one determine if the differences in performance are significant? The approach described here will work not only for comparing two MPR models, but to compare MPR with MLR, or MPR with ANN, etc. This can be carried out using the data from TaylorFit, but the calculation itself will have to be done in external software such as Excel. The importance ratio plot (Figure 4b) for D showed nitrate nitrogen having a negative effect on the index, whereas HQA had a significant positive correlation with D. Biological oxygen demand (BOD5) showed the highest RMS importance ratio value of 7.324 (Table  5), depicting the dominance of the variable on D. The species richness index (N) had an R 2 value of 0.330 and an MSE value of 8.39 ( Figure S5). The ANN derived model had 0.415 for validation and MSE of 9.29. HMS, HQA, and conductivity importance ratio plots (Figure 5b) showed a slight positive correlation with N index for all the observed values. Though nitrite nitrogen seemed to have a significant positive effect, but it also seemed to have a negative correlation with N. The nonlinear plot (Figures 4a and 5a) for D and N were both dependent on HMS and HQA. These plots captured the effect of human activities on the area effectively, which was clearly shown by the increasing rate of HMS index that decreased both D as well as N. Simultaneously, with increasing HQA, which demonstrated natural habitat around the region, appeared to result in surges of values for both the indices. Modeling quality for biodiversity indices, both D and N, supported the previous attempts on the relationship between biological diversity patterns and environmental factors in various types of ecosystems [41].  As can be seen, data can be modeled with different approaches, each producing its own performance results; but how can one determine if the differences in performance are significant? The approach described here will work not only for comparing two MPR models, but to compare MPR with MLR, or MPR with ANN, etc. This can be carried out using the data from TaylorFit, but the calculation itself will have to be done in external software such as Excel. As can be seen, data can be modeled with different approaches, each producing its own performance results; but how can one determine if the differences in performance are significant? The approach described here will work not only for comparing two MPR models, but to compare MPR with MLR, or MPR with ANN, etc. This can be carried out using the data from TaylorFit, but the calculation itself will have to be done in external software such as Excel.
Disparate models can be compared by using the ratio of their MSE values. Note that MSE = SSE/df. The df is unambiguous with MLR or MPR models. In the case of ANN models, there is disagreement in the literature about what the df should be. Some references indicate that df is much less than the number of weights in the ANN (which can be quite high); others indicate it should be higher. A conservative approach would be to just assume that df is equal to the number of weights in the ANN.
In any case, once the MSE is computed for each of the two models, the ratio of the model with the larger MSE (call it "model 1") and the model with the smaller MSE (model 2) will be the F-statistic for the comparison: Then, the probability that that value of F or larger can occur by chance can be computed from the F-statistic and the degrees of freedom for the two models: In Excel, for example, the p(F) can be computed using the built-in function "=FDIST(F,df1, df2)". If p(F) ≤ α, one may accept that model 2 has a significantly smaller (better) MSE. Table 4 shows the results of comparing each of the MPR models with the corresponding ANN models. The F-statistic shown is the ratio of the MSE for the MPR model to that of the ANN model. The lowest value of p(F) is 0.377 (for MIR), indicating that there is a 37.7% probability that the difference in MSE could occur by random error. Since this is much larger than the typical threshold α = 0.05, we must accept the null hypothesis that MSE for the MPR models may be equal to the MSE of the ANN models.
In other words, although the ANN models produce higher R 2 , they are not significantly higher, and we may conclude that the MPR models are as good as the ANN models in fitting these data.
The Pearson correlation [42] between independent variables is shown in Figure S5 (Supporting Information). It can be seen from the figure that most independent variables are not correlated, though this does not seem to be of concern, as MPR eliminates terms with high collinearity from the final model.

Conclusions
Based on our research, the macrophyte indices were modeled using the simple technique of multiple polynomial regression, and then sensitivity analysis was performed to demonstrate the influence of independent variables on the indices. The ecological indices (MIR, RMNI, and IBMR) performed better than the diversity indices (N and D). The MPR modeling can facilitate interpretation of ecological indices (MIR, RMNI, and IBMR) as well as the biological diversity indices (D and N). The same observations when modeled using artificial neural networks (ANNs) showed higher R 2 but had the complexity and difficulty in understanding the patterns of explanatory variables with the indices. Moreover, none of the R 2 values for the ANNs were statistically significantly higher than those for the MPR models, based on use of the F-statistic for comparing MSE values for the two modeling approaches. Consequently, inclusion of importance ratio analysis improved the comprehensive performance of the status of river water quality. Further, our MSE (mean square error) as well as df (degrees of freedom) were much lower. This demonstrates that MPR has added advantages in visualizations and sensitivity analysis.
Similar models can be developed to predict the eutrophication of other rivers where macrophytes play a major role. The simplicity of MPR models could be easily adopted in the analysis of "big" ecological data, and nonlinear sensitivity analysis could be used to develop relationships between the variables that are difficult to discover using black box methods.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/w13192708/s1, Figure S1: Comparison of predicted and observed MIR, Figure S2: Comparison of predicted and observed RMNI, Figure S3: Comparison of predicted and observed IBMR, Figure S4: Comparison of predicted and observed D, Figure S5: Comparison of predicted and observed N, Figure S6: Correlation between Independent Variables, Table S1: Equations for Macrophyte indices calculations.

Conflicts of Interest:
The authors declare no conflict of interest.