Nonlinear Simultaneous Equations for Individual-Tree Diameter Growth and Mortality Model of Natural Mongolian Oak Forests in Northeast China

A nonlinear equation system for individual tree diameter growth and mortality of natural Mongolian oak forests was developed based on 13,360 observations from 195 permanent sample plots in Northeast China. Weighted regression was used in a distance-independent diameter growth equation for dealing with heterocedasticity. Since diameter growth and mortality models have common predictors including the diameter at breast height (DBH), stand basal area (BA), basal-area-in-larger trees (BAL), and site index (SI), parameters were estimated using nonlinear three-stage least squares (N3SLS) and seemingly unrelated regression (SUR) which accounts for correlations of errors across models. The system equation provided better projection than individual fitting of the equation based on maximum likelihood estimation. Compared with the separate tree growth model, the simultaneous equations using N3SLS and SUR produced more efficient parameter estimation and smaller bias. Furthermore, N3SLS had more accurate projection. Overall, the simultaneous model will facilitate the growth and yield projection for better management of Mongolian oak forests in the region. OPEN ACCESS Forests 2015, 6 2262


Introduction
Oak is the largest tree species in China [1].Mongolian oak (Quercus mongolica Fisch.) is a main tree species of the secondary broad-leaved forest that grows mainly throughout Northeast China.It not only has economically important value for timber production, but also has the capacity to improve environment, soil and water conservation.Unfortunately, the wood quality of most Mongolian oak in Northeast China is poor because of short clear bole lengths, frequent occurrence of water sprouts and crooked stems.However, the wood quality of Mongolian oak growing in other regions, such as Japan, is excellent, and its wood has been widely used and also exported internationally [2].These facts indicate that the Mongolian oak growing in Northeast China has a potential to grow well like those species in Japan if adequate silvicultural treatments are applied.Forest management requires growth and yield models that enable prediction of the development of forest stands under different silvicultural systems.To overcome the negative effects of clearcutting, uneven-aged selection management has been practiced in Northeast China [3].The selection management needs the individual tree growth and yield models because they provide good simulation of growth and mortality for short-term projections, produce detailed information about stand structure and allow consideration of a wide variety of silvicultural treatments.However, there are no individual tree growth models of Mongolian oak forests for good forest management with the exception of the integrated stand growth model developed by using system equations [4].
Diameter growth and mortality models are basic components of individual tree growth modeling.Diameter at breast height (DBH) is the most important and the easiest tree factor to measure.Based on DBH, we can directly calculate stand basal area.In addition, because there is a strong relationship between DBH and tree height and volume, tree height and volume can be calculated through tree height-DBH curve and volume equation.Therefore, the establishment of a diameter growth model for forest growth forecast is of great significance.On the other hand, mortality is an important growth process in maintaining the population dynamic changes, species diversity, and structure diversity.A mortality model is important to understand forest evolution, and thus an essential component of the growth and yield models.Both individual tree diameter growth equations and mortality models have been largely used to predict annual or periodic growth and mortality of trees [5][6][7][8][9][10][11][12][13][14][15][16][17][18].Furthermore, from a biological perspective, the diameter growth and mortality equations of individual-tree growth model are related.Due to the highly internal correlation of these two equations, we can obtain more accurate and compatible parameter estimation through the establishment of joint estimation of simultaneous equations using nonlinear three stage least squares (N3SLS) and seemingly unrelated regression (SUR), which has been widely applied in forest growth and yield modeling [11,[19][20][21][22][23][24].
Therefore, the objective of this study was to develop nonlinear equations for individual tree diameter growth and mortality through simultaneous equations, in order to obtain more accurate and compatible parameter estimation for reasonable management of natural Mongolian oak forest in Northeast China.N3SLS and SUR were examined and compared with the traditional individual tree diameter growth and mortality equations based on maximum likelihood estimation (MLE).

Data
The data used to develop individual-tree diameter growth and mortality models for predicting natural Mongolian oak were collected from 195 re-measured permanent sample plots (PSP) in Wangqing Forestry Bureau, Jilin Province, Northeast China (Figure 1).The area of each plot was 0.06 hm 2 and all of the plots were natural pure forests.These PSPs were selected using the following criteria: tree-and stand-level characteristics of these sample plots have been remeasured periodically with an interval of ten years from 1997 to 2007, and plots are located in forests without any evidence of silvicultural treatments or any other forms of human interference, such as cutting and artificial regeneration.One hundred and sixty plots were selected for model calibration, and the rest 35 plots were for model validation.On each sample plot, physical site attributes were collected including geographic coordinates, slope, aspect, and elevation.Trees with DBH above 5 cm were measured with species, diameter, and status (recruitment, live, or dead).In total, we received 13,360 observations.For each plot, we calculated the following stand or tree factors: stand basal area (BA), number of trees per hectare (N), stand density index (SDI), stand quadratic diameter and its increment (ΔDq), and competition index.
Summary statistics of tree and plot level variables for model calibration are shown in Table 1.These plots are good representative for the research site.Mortality occurred in 152 plots out of 160 plots.The distribution of total trees and dead trees by diameter classes is shown in Figure 2, which indicates that mortality mainly occurred in the small diameter classes.

Model Development
The usually identified components of individual-tree growth models are diameter (or basal area) increment, mortality, and recruitment [25].Recruitment was not considered in the present study.At first, diameter growth and mortality models are established separately by maximum likelihood estimation, then due to the internal correlation of these two equations, simultaneous equations are built to obtain more accurate and compatible parameter estimation by three-stage least squares (N3SLS) [22] and seemingly unrelated regression (SUR) [23].

Diameter Growth Model
The diameter growth model predicts increment as a function of different variables describing tree size, competition and site conditions [13].After the test of a variety of forms of diameter growth equations, the individual tree diameter growth adding one in ten years (Δd+1) was selected as the predicted variable.The potential explanatory variables and their transformation were tested: Tree size: tree diameter at breast height (DBH).Distance-independent competition indices (COMP): basal area larger than subject trees (BAL), the ratio of objective tree's DBH and average DBH from stand basal area (RD), stand basal area per hectare (BA), number of trees per hectare (N), stand density index (SDI), the square sum of tree's diameter bigger than objective tree (DL), the ratio of objective tree's DBH and the biggest tree's DBH (DDM).
Site: site condition which is a composited function of slope, aspect, altitude [26] and site index which is the stand mean height at reference age and computed as Equation ( 1) [27] where H is stand mean height, T is stand mean age.The model of periodical diameter growth within ten years (Δd+1) was expressed as Equation ( 2): Non-significant variables were removed in the modeling process.In order to avoid multicollinearity, the value of the variance inflation factor (VIF) was computed for the groups of variables that best predicted.Those alternatives with a VIF larger than 10 were rejected, as suggested by Myers and Kutner et al. [28,29].

Mortality Model
The possibility of a tree dying in next 10 years was modeled by using logistic function Equation (3), which is widely applied for tree mortality [14][15][16][17].The candidate variables for the mortality model were numerous and diverse.Hamilton classified such variables in four groups: measures of individual tree size, tree competition, tree vitality and measures of stand density [18].Several possibilities for transformation of this variable exist and have been applied [15,18,[30][31][32].Only significant variables (p < 0.05) with VIF less than 10 were selected.

(
) where Pi is the probability of tree mortality, b0-bn are parameters to be estimated and x1-xn are explanatory variables.

Simultaneous Systems of Individual Diameter Growth and Mortality Models
Because a strong concomitant relationship exists between (Δd+1) and Pi, they used some same tree-and stand-level variables as predictors in both Equation (2) and Equation (3).It is likely to improve projection accuracy and ensure logical consistency.Treated in this manner, Equation (2) and Equation (3) were fitted as a nonlinear simultaneous system.Fitting this simultaneous system using maximum likelihood estimation would result in biased estimates of the coefficients, termed simultaneity bias [33][34][35].To remove simultaneity bias, the simultaneous system was firstly fitted using nonlinear three-stage least squares (N3SLS) regression, with the following three stages [24]: (1) Tree-and stand-level variables were used to estimate a first-stage estimate of (Δd+1) ((Δd+ 1)first) with a nonlinear model.Then a first-stage estimate of Pi (Pi first) was developed by the same basic procedure with a binary logistic model.(2) Some same or different tree-and stand-level variables were taken as predictor variables in both Equation (2) and Equation (3) in the second stage.Simultaneity bias across simultaneous systems of equations was removed by this.Simultaneous correlation and the correlation of errors, which arise in the two equations of the system, were then estimated by the error terms of these second-stage equations.(3) In the third stage, the estimated simultaneous correlation was used in which a seemingly unrelated regressions approach was taken to fit the two equations as a system.
For comparison, a consistent and asymptotically efficient estimator also was obtained if the system of nonlinear equations was jointly estimated under seemingly unrelated regression SUR directly [23].
Initial model residual diagnosis showed that there was heterocedasticity for diameter growth equation.To deal with the problem, an additional weighting factor was used.Heterocedasticity residual variance was modeled as a power function of an independent variable [36,37]: The q value was subsequently included in the weighting factor (w) as: 1/ q w DBH = .

Model Evaluation
For individual tree diameter growth model, the following statistics were calculated: mean bias (MB) and mean relative bias (MB%), which showed the bias of model prediction, root mean square error (RMSE), and relative root mean square error (RMSE%) which is expressed in the same units as the dependent variable, giving an idea of the mean error when using the model.
Mean bias: Mean relative bias: Root mean square error: Relative root mean square error: where i y is observed value, i y ˆ is predicted value, y is average value, n is the number of trees, and p is the number of parameters.
For individual tree mortality model, a probability threshold has to be determined to assign mortality.If the estimated probability of mortality exceeds the threshold then the tree is considered dead.The principle of using the maximum summary threshold (MST) was used in the study [38], namely taking the probability threshold when the sum of sensitivity and specificity is maximum (Table 2).Specificity is the ratio of the number of correctly predicted survived trees and observed survived trees.Sensitivity is the ratio of the number of correctly predicted dead trees and observed dead trees.Where a 11 is the number of correctly predicted dead trees, a 12 is the number of incorrectly predicted dead trees, a 21 is the number of incorrectly predicted survived trees, a 22 is the number of correctly predicted survived trees, A 1m is observed dead trees, A 2m is observed survived trees, A n1 is predicted dead trees, A n2 is predicted survived trees, A nm is total trees.
The area under the receiver operating characteristic (ROC) curve was calculated for the mortality model.According to Hosmer and Lemeshow [39], the area under the ROC curve is a threshold independent measure of model discrimination, where a value of 0.5 suggests no discrimination, 0.7-0.8suggests acceptable discrimination, and 0.8-0.9suggests excellent discrimination.To carefully check the bias among diameter classes, the deviations between predicted and observed values were tested by means of Pearson Chi-Squared Statistics [40,41].Chi-Squared values for classes within an attribute were calculated as: ( ) where NDobs, NDpred are the number of observed and predicted dead trees in a diameter class, respectively.Large chi-squared values provided evidence of lack of fit.When the chi-squared values we calculated are less than the critical chi-squared values with significant level p = 0.05, it means there is no significant difference between the probability of observed and predicted dead trees.
Model calibrations were carried out in R with the lme4 and systemfit package [42].

Individual-Tree Diameter Growth Model
Equation ( 9) was the final model and its parameter estimates were listed in Table 3.The most important explanatory variables were the log of tree DBH (lnDBH), stand basal area (BA), basal area larger than subject trees (BAL) and site index (SI).The standard errors of all explanatory variables are close to zero and the values of the variance inflation factor (VIF) are lower than 5.

Mortality Model
Estimated parameters for mortality model Equation (10) are shown in Table 4.The most important explanatory variables were tree size: inverse of diameter in breast height (1/DBH); competition factors: basal area larger than subject trees (BAL), stand basal area per hectare (BA); Site factor: site index (SI).All the parameter estimates were highly significant and all VIFs were smaller than 5.The threshold of mortality probability is 0.3 according to MST principle.

Simultaneous Systems of Equations Using N3SLS and SUR
Individual-tree diameter growth and mortality model are interrelated.For these above two models, the final instrumental variables used to obtain first stage estimates were DBH, BA, BAL, and SI.Therefore, the simultaneous equations are developed using seemingly unrelated regression N3SLS and SUR.
( ) where ( 1) d ∆ + is the first-stage estimate of (Δd+1) , i P ˆ is the first-stage estimate of Pi, and a0-a4 and b0-b4 are sets of species-specific parameters for the first-stage equations of (Δd+1) and Pi, respectively.Parameter estimates for the system of equations using N3SLS and SUR are presented in Tables 5  and 6, respectively.Both methods produce similar parameter estimates, but SUR produced smaller standard error than N3SLS.For the individual-tree diameter growth model, higher values of lnDBH and SI led to greater increments, while higher values of any of the other variables (BA, BAL) resulted in lower increments.It is consistent with the separate individual-tree diameter growth model.In the case of the mortality model, higher values of SI resulted in lower mortality, while higher values of any of the other variables (1/DBH, BAL, BA) resulted in greater mortality.It is not consistent with the separate mortality model but biologically reasonable since small trees and dense stands with high competition normally resulted in high mortality.Thus, the simultaneous systems of equations removed simultaneity bias and had more accurate estimated parameters.
The plots of weighted residuals against predicted values for the final individual-tree diameter growth equations (Figure 3) showed no evidence of heterogeneous variance over the full range of the predicted values and no systematic pattern in the variation of the residuals after the correction for heterocedasticity and non-normality.Moreover, the plots of weighted residuals against predicted values for the system of simultaneous equations using N3SLS (Figure 4) and SUR (Figure 5) showed, more obviously, no evidence of heterogeneous variance over the full range of the predicted values.It indicated that simultaneous equations using N3SLS and SUR were more precise than estimates of the separate individual-tree diameter growth model, N3SLS had the most accurate projection.The ROC curve for both fittings were very similar, and the area under the curve indicated that the model discriminated well [39] between live and dead trees (Figures 6-8).The AUC values of simultaneous equations and mortality model were 0.776, 0.761, 0.740, respectively.
The chi-squared statistics in Table 8 gave no evidence of lack of fit between predicted and observed values for both models.Examining the predicted and observed values by diameter class showed that simultaneous equations using N3SLS and SUR produced more accurate projection than separate mortality model (Figure 9), but there are larger errors in big diameter classes.It indicated that estimates of simultaneous equations were more precise than estimates of the separate mortality model for the range of data tested, and N3SLS had the most precise projection.

Discussion
Diameter growth and mortality models have long been recognized as a critical component in individual tree growth modeling.This study firstly presented a separate individual-tree diameter growth model and mortality model for Mongolia oak in Northeast China, based on permanent sample plots.Potential predictor variables were selected based on their biological importance for tree growth and mortality rather than simply fitting statistics because if a model does not make biological sense, it will not perform well for any data set other than that used for model development [18].The predictor variables for both models are related to tree size, tree competition and site conditions.After testing different independent variables, the following were included in the two models: the initial diameter at breast height (DBH), basal area of larger trees (BAL), stand basal area (BA), and site index (SI).Owing to same predictors and error correlation of both models, the simultaneous fitting was done through N3SLS and SUR.Possible correlated error components in the equations could be explained by the simultaneous fitting, because some tree-and stand-level variables depended on parameter estimates and were used in both equations.In contrast to individual fitting of the equations, this approach results in a decrease in the approximate standard errors of the coefficients and more accurate projections as shown in this study.Both stand-and tree-level system equations were reported to be superior to separate growth models [22,23,[43][44][45].However, the simultaneous fitting using N3SLS provided the more precise projection and similar parameter estimates than those obtained using SUR in the region.Because the bias from nonlinear N3SLS were smaller than those from SUR, which means N3SLS was more efficient than SUR in this region.Therefore, N3SLS is the most appropriate method for system equation estimate in the study.Recently, a system of equations was developed by Sattler and LeMay to predict crown length and crown radius for trees in structurally complex stands [24].Compatible branch and foliage biomass simultaneous equations were constructed based on the crown biomass mixed models, in which the crown biomass was divided into the two components of branch and foliage biomass [46].
In the simultaneous fitting using N3SLS, the individual tree diameter growth model showed that higher values of lnDBH and SI led to greater increments, while higher values of any of the other variables (BA and BAL) resulted in lower increments.The results confirmed that bigger trees and trees in better sites representing more vigor had higher diameter growth.The basal area of larger trees (BAL) confirmed that an increase in competition leads to a reduction in the diameter growth.BAL has been considered as a variable to capture competition for light [47].Stand basal area (BA) also has negative effects, which showed a decrease in the growth rate as the stand becomes more crowded.It is consistent with the separate individual tree diameter growth model [13,[48][49][50].In addition, the crown characteristics were also included in some individual tree diameter growth models [9,51].However, crown measurements were not used in this study, and no crown models were available for the species in the study area.Therefore, it will be interesting to test whether crown attributes would be beneficial to assess Mongolia oak stand dynamics in future research.
In the case of the mortality model from simultaneous fitting using N3SLS, higher values of SI representing good site condition resulted in lower mortality, while higher values of any of the other variables (1/DBH, BAL, BA) resulted in greater mortality.The results implied that the mortality probability of big trees is lower than that of small trees.Similar results have been reported for this species from Austria [15], Norway [41], and Finland [52].However, Monserud and Sterba [15] hypothesized an increased mortality rate for very large diameters due to senescence of over-mature trees, included DBH 2 in the model, and found it significant in predicting mortality for Norway spruce.The basal area of larger trees (BAL), stand basal area (BA) indicated that an increase in competition of stand leads to a higher possibility of mortality.This means that the denser the stand and the more suppressed the tree, the greater is the mortality rate [41,53,54].The size of the tree affects one-sided or asymmetric competition and occurs when larger trees have a disproportionate share of resources [55], so smaller trees are more likely to die due to stagnation.Light is the most common example of an asymmetric resource because of its directionality, but it is perhaps not the only one, since the degree of symmetry or asymmetry of mineral nutrients and water are not enough studied [47].The basal area of larger trees and the ratio of basal area of individual tree and stand have been considered as a variable to capture competition for light, and, hence, it is seen as a surrogate of one-sided competition [56], while stand basal area and number of trees is used to express two-side competition [53].This makes sense as these trees grow in relatively sparse stands with high availability of light, where competition seems to be for belowground resources [57].
Owing to lack of stand dominant height, site productivity in the study was realized by the stand mean height at a given age, which is a practical measure of site productivity in the past [58].It had significant effects on diameter and mortality in the both models as expected.There are still unexplained variances in both individual diameter growth and mortality models, and further model improvement will be considered in the future with the application of a mixed-effect model with plot as a random effect, for example.Height growth and recruitment models will also be needed.

Conclusions
A system of nonlinear simultaneous equations for individual tree growth and mortality in 10 years for natural Mongolian oak forests in Northeast China was developed.Parameters were estimated using nonlinear three-stage least squares (N3SLS) and seemingly unrelated regression (SUR), which provided better predictions than individual fitting of the equation, and N3SLS had the best prediction results.The model will facilitate the growth and yield prediction for better management of Mongolian oak forests in the region.

Figure 1 .
Figure 1.Map of Mongolian oak plots distribution in Wangqing Forestry Bureau, Northeast China.

Figure 2 .
Figure 2. Distribution of total trees and dead trees by diameter classes of Mongolian oak.

Figure 3 .
Figure 3. Weighted residuals and predicted values of diameter growth model.

Figure 4 .
Figure 4. Weighted residuals and predicted values of simultaneous equations for N3SLS.

Figure 5 .
Figure 5. Weighted residuals and predicted values of simultaneous equations for SUR.

Figure 6 .
Figure 6.ROC curve of separate mortality model.

Figure 8 .
Figure 8. ROC curve of simultaneous equations for SUR.

Figure 9 .
Figure 9. Observed and predicted mortality trees for simultaneous equations using N3SLS and SUR and separate mortality model.

Table 1 .
Summary statistics of tree-and plot-level variables for model calibration of Mongolian oak in Wangqing Forestry Bureau, Northeast China (n = 10,698).
since no dominant height was available in both class A and B inventory in China.

Table 2 .
Contingency table of the observed values and predicted values.

Table 3 .
The statistical analysis results for first-stage equations of individual-tree diameter growth model for Mongolian oak (n = 10,698).

Table 4 .
The statistical analysis results for first-stage equations of individual tree mortality model (n = 10,698).

Table 5 .
The parameter estimates from N3SLS regression of (Δd+1) and Pi in the system of equations for Mongolia oak (n = 10,698).

Table 6 .
The parameter estimates from SUR of (Δd+1) and Pi in the system of equations for Mongolia oak (n = 10,698).

Table 7 .
Estimated errors of separate diameter growth model and simultaneous equation using N3SLS and SUR.

Table 8 .
Test in predicting and observing the mortality probability for calibration and validation data sets.