Regionally Compatible Individual Tree Growth Model under the Combined Influence of Environment and Competition

To explore the effects of competition, site, and climate on the growth of Chinese fir individual tree diameter at breast height (DBH) and tree height (H), a regionally compatible individual tree growth model under the combined influence of environment and competition was constructed. Using continuous forest inventory (CFI) sample plot data from Fujian Province between 1993 and 2018, we constructed an individual tree DBH model and an H model based on re-parameterization (RP), BP neural network (BP), and random forest (RF), which compared the accuracy of the different modeling methods. The results showed that the inclusion of competition and environmental factors could improve the prediction accuracy of the model. Among the site factors, slope position (PW) had the most significant effect, followed by elevation (HB) and slope aspect (PX). Among the climate factors, the highest contribution was made by degree-days above 18 °C (DD18), followed by mean annual precipitation (MAP) and Hargreaves reference evaporation (Eref). The comparison results of the three modeling methods show that the RF model has the best fitting effect. The R2 of the individual DBH model based on RF is 0.849, RMSE is 1.691 cm, and MAE is 1.267 cm. The R2 of the individual H model based on RF is 0.845, RMSE is 1.267 m, and MAE is 1.153 m. The model constructed in this study has the advantages of environmental sensitivity, statistical reliability, and prediction efficiency. The results can provide theoretical support for management decision-making and harvest prediction of mixed uneven-aged forest.


Introduction
To obtain information on the potential future productivity of plantation forests, forest managers rely on various types of growth models [1]. The stand level growth model can be used to predict the development of even-aged forests, but for mixed forests or unevenaged forests the individual tree growth model is required [2]. The individual tree growth model is a bottom-up modeling method, starting from individual trees in the system and ending at the stand level, aiming to reveal and predict the growth mechanism of individual trees [3]. Compared with pure even-aged forests, the complexity, depth, and breadth of the forest growth simulations for mixed uneven-aged forests increased, and the scale gradually shifted from stands to individual trees [4]. And further compared with the whole stand model, the individual tree model is more applicable [5], promoting the widespread use of individual tree growth models in forest management [6].
According to the concept of mathematical statistics, the model can be divided into parametric model, non-parametric model, and semi-parametric model. A parametric of soil on plants is also well known; nitrogen, phosphorus, potassium and other nutrients, and humidity in soil vary with soil thickness and soil type [41]. Adding site factors into the growth model is of great significance for selecting reasonable afforestation tree species, increasing ecological stability, and supporting forest productivity [42,43]. Climate is an important external environment that determines the dynamic change of tree distribution and forest function at the regional scale [44,45]. Under the same site conditions, different climates will lead to different tree growth. Therefore, it is crucial to add climate factors into the growth model to improve its universality [46,47].
Close-to-nature forest management is a feasible theory and technology to improve the stability of forest ecosystems [48]. Building mixed uneven-aged forest is one of the main management measures of close-to-nature forest management. In order to build a mixed uneven-aged forest, it is necessary to conduct management operations on the forest, such as tending and thinning, adjusting the stand structure, etc. One of the important problems to be solved is to estimate the results of management operations, such as the harvest volume and biomass, etc. The solution is to build an individual tree growth model. There are many factors affecting forest growth, and it is very important to study the co-driving mechanism of competition and the environment to predict tree growth. However, there are much research on grassland and aquatic ecosystem, but there is little research on forest ecosystems [49,50] and most of them focus on the unilateral effects of stand density, climate factors, and tree species composition on tree growth. There are few individual tree growth models that include competition, site, and climate factors.
Chinese fir (Cunninghamia Lanceolata (Lamb.) Hook.) is one of the most important timber species in the south of China. The results of the ninth national continuous forest inventory (CFI) show that the area of Chinese fir plantations accounts for 1/4 of the total area of China's artificial arbor forest, and the planted area ranks the first among afforestation tree species [51]. As one of the six large forest regions in China, Fujian Province is located in the subtropical climate zone along the southeast coast of China, with complex and changeable climatic conditions and frequent extreme weather and climate events, making it a highly sensitive area to climate change [52]. Therefore, we took Chinese fir as the study tree species and Fujian Province as the study area to reveal the mechanism of competition factors (Comp), site factors (Site), and climate factors (Clim) effects on individual trees. The aims of this research study are: (1) to quantify the influence degree of factors, in order to select the relatively important factors in the individual tree growth model; (2) to explore the use of the re-parameterized (RP) method, BP, and RF algorithm to construct an individual tree growth model, which compares the accuracy and adaptability of different methods; and (3) to add the screened factors into the model and construct a regionally compatible individual tree growth model under the combined influence of environment and competition. The results of this study can provide references for the construction of individual tree growth models for different tree species in different regions, provide model support for growth harvest prediction of mixed uneven-aged forests, and also have indirect support significance for predicting forest carbon stocks and rational response to climate change.

Individual Diameter Growth Model Based on RP
With DBH as the dependent variable and T as the independent variable, the Gompertz model, Logistic model, Mitscherlich model, and Richards model were selected and fitted using the least squares method, respectively. SPSS software [53] was used to fit the data, and the R 2 and RMSE were used to evaluate the fitting results. As can be seen from the Supplementary Table S1, all of the above four basic models except Richards could converge, and all model parameters could pass the significance test. The Mitscherlich model had the best fitting effect, with the largest R 2 and the smallest RMSE. Therefore, the Mitscherlich model was selected as the optimal basic model, and Comp, Site, and Clim were added to construct the individual diameter growth model based on RP (DBH-RP model). Its parameter expression is shown in Equation (1).
On the basis of the previous analysis, R software was adopted to construct a stepwise regression using the step function. The re-preference of the independent variables was completed based on the Akaike information criterion (AIC). The final model independent variables obtained include sample plot tree density (N), slope position (PW), elevation (HB), soil thickness (TRHD), landform (DM), degree-days above 18 • C (DD18), and Hargreaves reference evaporation (Eref). Therefore, parameters ω0 and ω1 could be expressed as a linear combination of the above variables. The results of parameter estimation of the DBH-RP model estimation are shown in Table 1. After processing, the DBH-RP model including T, Comp, Site, and Clim with statistically significant parameters was as shown in Equation (4).
In the formula, ω 0 and ω 1 were parameters. The formula is shown in Equations (2) and (3). The values of the coefficients in the formula are shown in Table 1.
In the test data set, R 2 of the DBH-RP model was 0.621 and RMSE was 2.813 cm; R 2 was 12.3% higher than that of the basic model, and RMSE was 9.3% lower. The predicted change in DBH with T of the DBH-RP model is shown in Figure 1. data, and the R 2 and RMSE were used to evaluate the fitting results. As can be seen from the Supplementary Table S1, all of the above four basic models except Richards could converge, and all model parameters could pass the significance test. The Mitscherlich model  had the best fitting effect, with the largest R 2 and the smallest RMSE. Therefore, the Mitscherlich model was selected as the optimal basic model, and Comp, Site, and Clim were added to construct the individual diameter growth model based on RP (DBH-RP model). Its parameter expression is shown in Equation (1).
In the formula, ω0 and ω1 were parameters. The formula is shown in Equation (2) and Equation (3). The values of the coefficients in the formula are shown in Table 1.
In the test data set, R 2 of the DBH-RP model was 0.621 and RMSE was 2.813 cm; R 2 was 12.3% higher than that of the basic model, and RMSE was 9.3% lower. The predicted change in DBH with T of the DBH-RP model is shown in Figure 1.  The residual analysis was used to evaluate the model, and statistical diagnosis was made on the model through the state of the residual graph. The residual comparison between the Mitscherlich model and the DBH-RP model is shown in Figure 2. The residual analysis was used to evaluate the model, and statistical diagnosis was made on the model through the state of the residual graph. The residual comparison between the Mitscherlich model and the DBH-RP model is shown in Figure 2. As can be seen from Figure 2, the residual was evenly distributed between −6.0 cm and 6.0 cm, but after the addition of Comp, Site, and Clim and re-parameterization, more residuals were reduced to between −4.0 cm and 4.0 cm on the 5-20 cm DBH interval, but on the DBH interval beyond greater than 20 cm, residuals did not change significantly. There was no heteroscedasticity in either the Mitscherlich model or the RP model.

Individual Diameter Growth Model Based on BP
The individual diameter growth model based on BP (DBH-BP) was trained. Table 2 shows the fitting results of the sub-model and the whole model based on BP.
When the independent variable only includes T, the model can explain 63.8% of the DBH variable, and it is concluded that T was the most important factor affecting the growth of DBH. However, after Comp, Site, and Clim were introduced into the model, the predictive ability of the model improved. R 2 increased by 7.9%, 12.6%, and 4.1%, respectively, while MAE and RMSE decreased. The predictive accuracy of the whole model composed of all factors reached the maximum, which could explain more than 80.8% of the DBH variable. These results indicate that the introduction of Comp, Site, and Clim could improve the prediction accuracy of the individual tree growth model. The predicted change in DBH with T of the DBH-BP model is shown in Supplementary Figure S1. Meanwhile, it can be seen from the distribution diagram of the residuals of the four models (Figures 3 and S2) that they were evenly distributed without obvious heterogeneity. However, with the addition of Comp, Site, and Clim to the model, the residual error of the model was gradually reduced from −6.0 cm and 6.0 cm to −3.0 cm and 3.0 cm, and the model accuracy was improved. As can be seen from Figure 2, the residual was evenly distributed between −6.0 cm and 6.0 cm, but after the addition of Comp, Site, and Clim and re-parameterization, more residuals were reduced to between −4.0 cm and 4.0 cm on the 5-20 cm DBH interval, but on the DBH interval beyond greater than 20 cm, residuals did not change significantly. There was no heteroscedasticity in either the Mitscherlich model or the RP model.

Individual Diameter Growth Model Based on BP
The individual diameter growth model based on BP (DBH-BP) was trained. Table 2 shows the fitting results of the sub-model and the whole model based on BP. When the independent variable only includes T, the model can explain 63.8% of the DBH variable, and it is concluded that T was the most important factor affecting the growth of DBH. However, after Comp, Site, and Clim were introduced into the model, the predictive ability of the model improved. R 2 increased by 7.9%, 12.6%, and 4.1%, respectively, while MAE and RMSE decreased. The predictive accuracy of the whole model composed of all factors reached the maximum, which could explain more than 80.8% of the DBH variable. These results indicate that the introduction of Comp, Site, and Clim could improve the prediction accuracy of the individual tree growth model. The predicted change in DBH with T of the DBH-BP model is shown in Supplementary Figure S1.
Meanwhile, it can be seen from the distribution diagram of the residuals of the four models (Figures 3 and S2) that they were evenly distributed without obvious heterogeneity. However, with the addition of Comp, Site, and Clim to the model, the residual error of the model was gradually reduced from −6.0 cm and 6.0 cm to −3.0 cm and 3.0 cm, and the model accuracy was improved.

Individual Diameter Growth Model Based on RF
RF model was used to build an individual diameter growth model (DBH-RF). Table 3 shows the sub-model and whole model fitting results using RF model.

Individual Diameter Growth Model Based on RF
RF model was used to build an individual diameter growth model (DBH-RF). Table 3 shows the sub-model and whole model fitting results using RF model.  Table 3 shows that when the independent variable contained only T, the model could explain 67.9% of the DBH. However, after Comp, Site, and Clim were introduced into the model, the predictive ability was improved, R 2 was increased by 8.0%, 12.7%, and 2.8%, respectively, and MAE and RMSE were decreased. The prediction accuracy of the whole model composed of all factors reached the maximum, which could explain more than 84.9% of the DBH variable, and MAE and RMSE were the minimum. The predicted change in DBH with T of the DBH-RF model is shown in Supplementary Figure S3.
At the same time, it can be seen from the distribution diagram of the prediction and residuals of the four models (Figures 4 and S4) that the residuals of the four models were evenly distributed without obvious heterogeneity. However, with the addition of Comp, Site, and Clim, it can be seen that the model residuals gradually shrunk from −4.0 cm and 6.0 cm to −3.0 cm and 2.0 cm, and the model residuals decreased significantly.   Table 3 shows that when the independent variable contained only T, the model could explain 67.9% of the DBH. However, after Comp, Site, and Clim were introduced into the model, the predictive ability was improved, R 2 was increased by 8.0%, 12.7%, and 2.8%, respectively, and MAE and RMSE were decreased. The prediction accuracy of the whole model composed of all factors reached the maximum, which could explain more than 84.9% of the DBH variable, and MAE and RMSE were the minimum. The predicted change in DBH with T of the DBH-RF model is shown in Supplementary Figure S3.
At the same time, it can be seen from the distribution diagram of the prediction and residuals of the four models (Figures 4 and S4) that the residuals of the four models were evenly distributed without obvious heterogeneity. However, with the addition of Comp, Site, and Clim, it can be seen that the model residuals gradually shrunk from −4.0 cm and 6.0 cm to −3.0 cm and 2.0 cm, and the model residuals decreased significantly.

Individual Height Growth Model
The conditions and tree species of the study area remained unchanged, and it was convenient and practical to build a machine learning model for a specific research object in the nntool toolbox of MATLAB [54]. Referring to the method of individual diameter

Individual Height Growth Model
The conditions and tree species of the study area remained unchanged, and it was convenient and practical to build a machine learning model for a specific research object in the nntool toolbox of MATLAB [54]. Referring to the method of individual diameter growth model construction, the H data set was imported, and the RP, BP, and RF models were directly used to train the H data set to verify the applicability of the model.

Individual Height Growth Model Based on RP
The Gompertz model, Logistic model, Mitscherlich model, and Richards model were selected to fit H models using the least squares method. The R 2 and RMSE were used to evaluate the fitting results. As shown by the fitting results in Supplementary Table S2, the Mitscherlich model had the best fitting effect. The parameter estimations of ω0 and ω1 are shown in Table 4. The re-parameterized model including Comp, Site, and Clim is as follows.
In the formula, ω 0 and ω 1 are parameters. The formulae are shown in Equations (2) and (3). The values of the coefficients in the formula are shown in Table 4.
The evaluation index of the H-RP model was calculated, and the R 2 and RMSE values were 0.683 and 2.102 m. As can be seen from the results, compared to the Mitscherlich model, the R 2 of the H-RP model improved by 4.1%, and the RMSE decreased by 5.0%. The predicted H changes with T and the residual diagram of the H-RP model are shown in Figure 5, indicating that the H-RP model could reflect the growth process of tree height. From the residual of the model, it can be seen that the residuals of the RP model were uniformly distributed between −3 m and 3 m, and the residual diagram did not have heteroscedasticity.  Figure S5. When calculating the evaluation index, the R 2 was 0.731, RMSE was 1.857 m, and MAE was 1.686 m, which met the accuracy requirements in practice. On the H test set, compared to the RP model, the R 2 of the BP model improved by 7.0% and the RMSE decreased by 11.7%. From the residual of the model, it can be seen that the residuals of the BP model were uniformly distributed evaluation index, the R 2 was 0.731, RMSE was 1.857 m, and MAE was 1.686 m, which met the accuracy requirements in practice. On the H test set, compared to the RP model, the R 2 of the BP model improved by 7.0% and the RMSE decreased by 11.7%. From the residual of the model, it can be seen that the residuals of the BP model were uniformly distributed between −4 m and 3 m, and the residual was evenly distributed between −3 m and 3 m; there was no heteroscedasticity.

Individual Height Growth Model Based on RF
The RF model was used to construct an individual tree height model (H-RF) containing Comp, Site, and Clim, and the predicted H changes with T and model residuals were obtained, as shown in Supplementary Figure S6. When calculating the evaluation index, the R 2 of the H-RF model for tree height was 0.845, RMSE was 1.267 m, and MAE was 1.153 m. On the H test set, compared to the RP model, the R 2 of the RF model improved by 23.7%, and the RMSE decreased by 39.7%. These results clearly demonstrate the advantages and potential of machine learning algorithms over traditional models. It can be seen that the residuals of the RF model were uniformly distributed between −2 m and 2 m, model accuracy was significantly improved, and there was no heteroscedasticity in the model.

Performance Evaluations of Three Optimal Models
In order to compare the prediction accuracy of different whole models in individual tree growth modeling, the prediction accuracy of the three models was compared and analyzed (Supplementary Table S3). Figure 6 shows a comparison between the predicted values and the actual values of the three models for DBH and H. It can be seen that the predicted results of the three models were close to the actual values. The BP and RF models were effective for the use of data sets containing a large number of environmental predictors with complex interactions and nonlinear relationships. Among them, the RF model had the best prediction effect with good generalization ability and statistical reliability.

Discussion
Existing studies on individual tree growth models are constructed using different methods for Chinese fir. Using multiple stepwise regression methods to build the linear mixed-effects model with sample plots as random effects, the R 2 of the model is 0.676 [55]. Considering climate and competition, an individual tree growth model at DBH was constructed using the Bayesian model with an R 2 of 0.6121 [56]. The results of our study

Discussion
Existing studies on individual tree growth models are constructed using different methods for Chinese fir. Using multiple stepwise regression methods to build the linear mixed-effects model with sample plots as random effects, the R 2 of the model is 0.676 [55]. Considering climate and competition, an individual tree growth model at DBH was constructed using the Bayesian model with an R 2 of 0.6121 [56]. The results of our study clearly demonstrated the advantages and potential of machine learning algorithms in the individual growth modeling of Chinese fir. Compared with the formula model, the machine learning model can model the complex nonlinear relationship without being restricted by statistical assumptions. Especially when the model contains independent variables and a large amount of data, it is more convenient to use machine learning in the modeling process and has obvious advantages.
The influence of environmental factors on the growth of an individual tree was analyzed through the results of 10-fold cross-validation. Supplementary Figure S7 shows that, at the level of individual trees, the age of individual trees is the main factor affecting the growth of trees. The relative importance of individual trees was extremely high, exceeding 60%, and the relative importance of competition factors to the model was about 8%.
In Site, PW had the most significant effect. In the BP and RF models, the model accuracy was increased by about 4%. Followed by HB and TRHD, the model accuracy increased by about 3%, slope aspect (PX) and DM had the least effect, and the contribution rate was about 1%. This result was not entirely consistent with those of previous studies [57]. The main reason for this is the difference in the study area. Fujian Province has more hills and mountains, and the vertical variation in heat and precipitation is more pronounced, which determines the growth difference of the study subjects; therefore, PW and HB become important factors.
PW represents the soil erosion and accumulation capacity. Generally speaking, soil temperature gradually decreases from top to bottom, while moisture gradually increases. Chinese fir is an acidic positive tree species that likes deep, fertile, and moist soil and has good drainage conditions, and a down slope position is more conducive. HB represents temperature, elevation increase, and temperature decrease, which is not conducive to the growth of Chinese fir at DBH [58]. TRHD represents soil fertility and water holding capacity. The thicker the soil, the higher the comprehensive fertility and water holding capacity, promoting the growth of DBH. This conclusion is consistent with the findings of Monserud et al. [59]. PX represents light, and light is more sufficient on sunny slope than on shady slopes, and the photosynthesis of Chinese fir is stronger, promoting DBH growth; this conclusion is consistent with the findings of LU [60]. DM represents the soil nutrient distribution and temperature and humidity conditions in the woodland. Compared with hills, Chinese fir has a higher material accumulation capacity in mountains with abundant rainfall and high air humidity [61].
We also found that the partial dependence of a variable on DBH growth is highly correlated with the relative importance of the variable (by comparing the ranking results of characteristic importance). Specifically, when the relative importance of a variable is greater, the average annual diameter growth of individual trees changes more sharply with the change in the variable. When the relative importance of a variable is less, the average annual diameter growth of individual trees changes more smoothly with the change in the variable. At the same time, competition factors, site factors, and climate factors have interactive effects on the growth of individual trees.

Study Area
The study area is Fujian Province, China, and the geographical location map is shown in Figure 7. Fujian Province is located on the southeast coast of China, 115 • 50 ~120 • 40 E, 23 • 33 ~28 • 20 N. The average annual temperature of the province is 17-21 • C, and it has an annual rainfall of 1400-2000 mm. The altitude of the province is high in the northwest and low in the southeast. Most of the province is hilly, with few plains, and hills and valleys account for more than 80% of the total area. The southeastern coastal region has a south subtropical climate, while the northeast, north, and west parts of the province have a central subtropical climate. The climate characteristics vary greatly between regions in the province. With a forest area of 81,158 km 2

Study Area
The study area is Fujian Province, China, and the geographical location map is shown in Figure 7. Fujian Province is located on the southeast coast of China, 115°50′~120°40′ E, 23°33′~28°20′ N. The average annual temperature of the province is 17-21 °C, and it has an annual rainfall of 1400-2000 mm. The altitude of the province is high in the northwest and low in the southeast. Most of the province is hilly, with few plains, and hills and valleys account for more than 80% of the total area. The southeastern coastal region has a south subtropical climate, while the northeast, north, and west parts of the province have a central subtropical climate. The climate characteristics vary greatly between regions in the province. With a forest area of 81,158 km 2

Climate Data
The climate factors' data for the study were generated using the ClimateAP (http: //ClimateAP.net, accessed on 22 March 2023) [63] software of Microsoft, providing climate data covering the entire Asia-Pacific region. By entering the longitude, latitude, and elevation annual-scale, climate data for the sample plots of the survey year can be obtained. A total of 17 climate factors were used, and the basic overview is shown in Supplementary Table S4.

Sample Plot Data
The sample plot data were collected from the CFI of Fujian Province from 1993 to 2018. According to the "Technical regulations for continuous forest inventory" [64], the sample plots were generally square, with a side length of 25.82 m and an area of 667 m 2 . The measured attributes included individual DBH, H, age, sample plot tree density, and site factors. Site factors include landform, elevation, slope aspect, slope position, slope gradient, soil name, and soil thickness. The summary of sample plot data and individual tree data obtained after data collation is shown in Table 5.
The DBH data set covered 24 plots in different counties from 1993 to 2008. The plots were repeatedly measured four times with an interval of 5 years, and a total of 1432 individual trees were collected. According to the rules of the CFI, H is generally determined by selecting 3-5 trees nearest to the center of the sample plot as a sample for height determination. The study selected 36 plots in different counties; the plots were repeated 3 times with an interval of 5 years, and a total of 357 individual trees were collected. The H data set covered the years from 2008 to 2018. We obtained the original data set. Due to the old age of some data in the original data set and errors in tree status records, this study followed the Pauta Criterion [65], taking the sample plot as the unit and eliminating non-relevant data according to the triple standard deviation method. Finally, we obtained 1189 DBH data and 307 H data of individual trees. In order to facilitate model construction, DBH data set and H data set were stored as Excel tables in the same format. The DBH data set was taken as an example (Supplementary Table S5).
The train_test_split function was used in the scikit-learn machine learning library in Python language to randomly partition the data set [66], with 80% as the training set for model training and 20% as the test set for validating and evaluating the reliability and generalization ability of the model.

Model Independent Variable Selection
The independent variables of the individual tree growth model included four groups of factors: age (T), Comp, Site, and Clim.
According to existing studies [67,68], N is the most important factor affecting the level of competition in a sample plot. N affects the diameter growth of trees, as shown by the decrease in diameter with increasing N. At the same time, N changes light conditions and the soil environment, which also affects the content of soil organic matter. Therefore, N was chosen to quantify competition in the sample plot.
We used RF feature importance evaluation to sort Site and Clim. The basic idea of RF feature importance evaluation is [69]: calculate the reduction value of the Gini coefficient (DGi) for Site and Clim at node splitting, sum the DGi of all nodes in the random forest, then after averaging over all trees to calculate the importance of Site and Clim, and finally normalize and sort according to the importance size. The results as shown in Figure 8.
According to Figure 8, the top 5 factors with the highest relative importance of features in Site were PW, HB, TRHD, PX, and DM. The top 5 factors with the highest relative importance of features in Clim were DD18, MAP, Eref, MCMT, and MAT.
Secondly, we calculated the Pearson correlation coefficient (Pearson's r) between each of the Site and the Clim. If r > 0.7, the two factors are strongly correlated; to avoid multicollinearity violating statistical assumptions, only one of the two factors with a higher contribution rate was selected. A thermal map of the correlation coefficient was obtained, as shown in Figure 9.
According to Figure 9, there was no strong correlation between Site variables. In Clim variables, MAT was strongly correlated with Eref, MCMT, and MAP; meanwhile, MCMT was strongly correlated with DD18 and Eref. Therefore, MAT and MCMT were excluded. feature importance evaluation is [69]: calculate the reduction value of the Gini coefficient (DGi) for Site and Clim at node splitting, sum the DGi of all nodes in the random forest, then after averaging over all trees to calculate the importance of Site and Clim, and finally normalize and sort according to the importance size. The results as shown in Figure 8.
According to Figure 8, the top 5 factors with the highest relative importance of features in Site were PW, HB, TRHD, PX, and DM. The top 5 factors with the highest relative importance of features in Clim were DD18, MAP, Eref, MCMT, and MAT. Secondly, we calculated the Pearson correlation coefficient (Pearson's r) between each of the Site and the Clim. If r > 0.7, the two factors are strongly correlated; to avoid multicollinearity violating statistical assumptions, only one of the two factors with a higher contribution rate was selected. A thermal map of the correlation coefficient was obtained, as shown in Figure 9.
According to Figure 9, there was no strong correlation between Site variables. In Clim variables, MAT was strongly correlated with Eref, MCMT, and MAP; meanwhile, MCMT was strongly correlated with DD18 and Eref. Therefore, MAT and MCMT were excluded. After feature importance evaluation ranking and Pearson's r screening, we discarded unimportant factors and variables with strong autocorrelation. The final Site included PW, HB, TRHD, PX, and DM. The final Clim included DD18, MAP, and Eref. Therefore, the individual tree growth model driven by competition and environment can be expressed as Equation (6).
where y represents the DBH model or H model, f ( ) represents modeling method, T represents age, Comp represents competition factors, Site represents site factors, and Clim represents climate factors.
To analyze the contribution of competition and environment to the dependent variables of DBH and H step-by-step, the following four types of models were, respectively, compared and analyzed: (1) Sub-model 1, that is, the independent variable only contains T; (2) Sub-model 2, that is, the independent variable contains T and Comp; (3) Sub-model 3, that is, the independent variable contains T, Comp and Site; (4) Whole model, that is, the independent variable contains T, Comp, Site, and Clim; Therefore, the individual tree growth model driven by competition and environment can be expressed as Equation (6).
where y represents the DBH model or H model, f( ) represents modeling method, T represents age, Comp represents competition factors, Site represents site factors, and Clim represents climate factors.
To analyze the contribution of competition and environment to the dependent variables of DBH and H step-by-step, the following four types of models were, respectively, compared and analyzed: (1) Sub-model 1, that is, the independent variable only contains T; (2) Sub-model 2, that is, the independent variable contains T and Comp; (3) Sub-model 3, that is, the independent variable contains T, Comp and Site; (4) Whole model, that is, the independent variable contains T, Comp, Site, and Clim;

Re-Parameterized Model
In the base model, stand growth is completely age dependent. To quantify the effects of competition, stand, and climate on individual tree growth, an RP method was adopted, in which the parameters (ω0 and ω1) in the base model were expressed as linear combinations of Comp, Site, and Clim. The RP model construction consisted of the following main steps: Step 1: Base model selection. The Gompertz, Logistic, Mitscherlich, and Richards models were selected as alternative base models. The four base models were fitted using the study data set, and the optimal base model was selected by solving the model parameters using the least squares method. The alternative base models are shown in Table 6. Table 6. Alternative basic models.

ID
Model Equation 1 1 Gompertz y = a * e −b * e −c * t 2 Logistic Step 2. Stepwise regression screening variables. After the basic model was selected, to avoid model multicollinearity, the stepwise regression method was used to remove the largest variable of VIF (Variance Inflation Factor) from the model and then re-fit it. The process was repeated until the VIF < 5 (in general, no collinearity between independent variables can be considered if VIF < 5).
Step 3. Construct the parameterized model. The selected variables were added into the basic model as covariates and fitted through the NLS (Nonlinear Least Squares) function of R software [70] to obtain the form with the best fitting effect.

Back Propagation Neural Network
The essence of a neural network is the approximation of an algorithm or function. It is comprises multiple layers of neurons. Each layer of neurons receives the original input or input from the neuron of the previous layer and then performs mathematical operations and outputs the operation results to the next layer of neurons. The calculation formula of the neuron output value n is shown in Equation (11).
where t is the input number of neurons in this layer, ω i is the weight of neuron input in layer i, m i is the input value of neuron in layer i, and A is a predefined nonlinear function. The BP model was selected in this study, the model was composed of an input layer, output layer, and several hidden layers; each layer contained several neurons, and the neurons between layers were connected by a weight or threshold value. Each layer of neurons affected only the next layer of neurons, and the neurons in the same layer were not connected with each other. The model structure is shown in Figure 10.
layer , is the input value of neuron in layer , and A is a predefined n tion.
The BP model was selected in this study, the model was composed of output layer, and several hidden layers; each layer contained several neu neurons between layers were connected by a weight or threshold value. Eac rons affected only the next layer of neurons, and the neurons in the same connected with each other. The model structure is shown in Figure 10. The BP model construction process was as follows: Step 1. Determine the training sample. The data set was divided, with 80% as the training set and 20% as the test set. T, Comp, Site, and Clim were selected as the input variables of the model, and DBH and H were selected as the output variables, respectively.
Step 2. Determine the network structure and parameters. We chose three layers of neural networks, including one input layer, one hidden layer, and one output layer. Equation (12) was used to calculate the number of hidden layer nodes, the three-layer network structure is 10:6:1. Set the transfer function between the input layer and the hidden layer as Log-sigmoid, and set the transfer function between the hidden layer and the output layer as Purelin.
where n is the number of neurons in the input layer, o is the number of neurons in the output layer, and m is any integer between 1 and 10.
Step 3. Network training. The BP model was trained using the Python language [71]. Select the training data set, call the BP model for training, and compare the influence of different inputs on the fitting accuracy of DBH and H.
In this study, the model results were optimal when the learning rate was set to 0.01, the target accuracy was 0.001, and the maximum number of iterations was 1000.

Random Forest
The RF model is composed of multiple decision trees, and each node in the decision tree is a condition about a feature. The Bootstrap resampling method was used to extract multiple samples from the original data, conduct decision tree modeling for each sample, and combine the prediction of multiple decision trees to jointly predict the results. To compare the difference, the RF model was set up to use the same data set partitioning ratio as the BP model, with 80% as the training set and 20% as the test set, and the mean using the same input variables as the BP model, including T, Comp, Site, and Clim. The RF model was trained using the Python language. The training process of the RF model is as follows.
Step 1. The original data sample content is N. Randomly generate M variables for a binary tree on N nodes. The choice of binary tree variables satisfies the principle of minimum Gini impurity.
Step 2. Use the bootstrap combination method to sample with replacement n sample sets in M to form n decision trees; then, the unsampled samples were used for the prediction of a single decision tree.
Step 3. According to the random forest composed of n decision trees, the final result is the average output of each decision tree. The RF model has two important parameters to set: the number of decision trees (ntree) and the number of variables randomly selected by tree nodes (mtry). Generally speaking, the overall error rate tended to be stable after ntree reached 500. Meanwhile, Breiman [72] suggested that for regression problems, the default value of mtry should be set to 1/3 of the number of all arguments. Taking this as a reference, we take ntree = 1, 20 . . . 500 to tune ntree. This study used at most 10 independent variables; therefore, 1 ≤ mtry ≤ 10, mtry = 1, 2 . . . 10. The results show that the model works best when ntree = 160 and mtry = 4.

Model Evaluation
The following fitting statistical indicators were selected as the model selection criteria. The calculation formulas of each test indicator are shown in Equations (13)-(15): (1) Coefficient of Correlation (R 2 ); (2) Root Mean Square Error (RMSE); (3) Mean Absolute Error (MAE) where the y i is the actual value of i,ŷ i is the predicted value of i, y is the dependent variable of the actual value of the mean, and n is the number of samples. In addition to the above statistical indicators, we also chose residual plots for the model evaluation. As an important regression diagnostic quantity, the residuals imply important information about the model assumptions. The analysis of residuals can examine the following issues [73]: the feasibility of the linearity assumption of the regression function; the reasonableness of the assumption of equal variance of random errors; the reasonableness of the assumption of independence of errors; the feasibility of the assumption of normal distribution of errors; the presence of outliers in the observations; and the omission of certain important independent variables in data collection or model fitting.

Conclusions
The re-parameterized model and machine learning model were used to construct the DBH model and H model of Chinese fir. The input variables included T, Comp, Site, and Clim. The results showed that: (1) The addition of competition and environmental factors could improve the prediction accuracy of the individual tree growth model. In terms of site factors, PW had the most significant effect, followed by HB and PX. In terms of climate factors, the DD18 contribution rate was the highest, followed by MAP, and Eref. (2) Among the three models, the RF model fits best. The model constructed in this study could well reflect the growth process of individual DBH and H of Chinese fir in Fujian Province and has reference significance for the growth models of other provinces and other tree species.
Climate factors were added into the individual tree growth model, which fully reflected the spatial differences of climate factors and solved the key problem of applicability of the individual tree growth model in different regions. The model was applied to the forest management decision support system. The results showed that the individual tree growth model, including T, Comp, Site, and Clim, could be used to predict the stand harvest in different sites and climates and provide decision support for harvest prediction and management method selection of mixed uneven-aged forests.
Since the area of Chinese fir plantation forest is the largest in terms of plantation forests in China, this paper only took Chinese fir as an example to construct the individual tree growth model. In future research, the modeling approach in this paper could be applied to the construction of individual tree growth models for other tree species, such as Masson pine (Pinus massoniana Lamb.) and broadleaf species. We also recommend using longer study years to investigate the variables to provide more reliable results.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants12142697/s1, Figure S1: DBH-BP model predicts the change in DBH with T. Figure S2: Residual diagram of sub-model 2 and sub-model 3 based on the BP model. Figure S3: DBH-RF model predicts the change of DBH with T. Figure S4: Residual diagrams of sub-model 2 and sub-model 3 based on the RF model. Figure S5: H-BP model predicted H changes with T and residual distribution diagrams. Figure S6: H-RF model predicted H changes with age and residual distribution diagram. Figure S7: Accuracy changes of BP and RF models with the number of variables. Table S1: DBH basic model parameter fitting value results and fitting index results.

Data Availability Statement:
The data is not publicly available because it relates to the copyright of the data provider.

Conflicts of Interest:
The authors declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.