A Nonlinear Mixed-Effects Height-Diameter Model with Interaction Effects of Stand Density and Site Index for Larix olgensis in Northeast China

: Tree height is a basic input variable in various forest models, such as growth and yield models, biomass models, and carbon budget models, which serve as very important tools for the informed decision-making in forestry. The height-diameter model is the most important component of the growth and yield models and forest simulators. We developed the nonlinear mixed-effects height-diameter model with the interaction effects of stand density and site index introduced using data from 765 Larix olgensis trees in Jingouling forest farm of the Wangqing Forest Bureau in northeast China. Among the various basic versatile functions evaluated, a simple exponential growth function ﬁtted the data adequately well, and this was then expanded through the introduction of the variables describing the interaction effects of the stand density and site index on the height-diameter relationship. Sample plot-level random effects were included into this model through mixed-effects modeling. The results showed that the random effect of the stand density on the height-diameter relationship was substantially different at different classes of the site index, and the random effect of the site index was different for the different stand density classes. The nonlinear mixed-effects (NLME) height-diameter model coping with the interaction effects of the stand density and site index had a better performance than those of the NLME models with the random effect of the single variable of stand density or site index. To conclude, the inclusion of the interaction effects of stand density and site index could signiﬁcantly improve the prediction accuracy of the height-diameter model for Larix olgensis Henry. The proposed model with the interactive random effects included can be applied for the accurate prediction of Larix olgensis tree height in northeast China.


Introduction
The height-diameter model is the most important component of the growth and yield models [1] and forest simulators as tree height, which is usually estimated from the heightdiameter model, is a basic input variable for a variety of forest models, such as growth and yield models, biomass models, and carbon budget models [2]. Relative to the diameter at breast height (DBH), the measurement of the tree height is time-consuming and laborious. Therefore, in forest inventory, only a part of the sample trees have their heights measured in the sample plot, and the heights of all the other trees from the same species are predicted with the established diameter-height model [3].

1.
Quantify the interaction effects of the stand density and site index on the tree height of Larix olgensis Henry.

2.
Determine whether there would be significant interaction effects of the stand density and site index on the height-diameter relationship. 3.
Establish the nonlinear mixed-effects height-diameter model with the interaction effects of stand density and site index included for Larix olgensis in northeast China.

Study Area and Data
We used the data acquired from twenty permanent sample plots (PSPs) established in Larix olgensis plantations located in the Jingouling forest farm of Wangqing Forest Bureau of northeast China (Figure 1 and Table 1). Sample plot sizes vary from 775 m 2 to 2500 m 2 and represent a large within-sample plot variation of the height-diameter relationship. Twenty sample plots were nested within the five blocks (Block) with different site qualities, and each block contained four sample plots (Plot). PSPs were established in years between 1961 to 1964. Most of the sample plots were in the forest stands with different site conditions (Table 1) of larch forests [16]. Therefore, it is of a great significance to establish a high precision height-diameter model for Larix olgensis Henry.
With reference to the points above, this study aims to: 1. Quantify the interaction effects of the stand density and site index on the tree height of Larix olgensis Henry. 2. Determine whether there would be significant interaction effects of the stand density and site index on the height-diameter relationship. 3. Establish the nonlinear mixed-effects height-diameter model with the interaction effects of stand density and site index included for Larix olgensis in northeast China.

Study Area and Data
We used the data acquired from twenty permanent sample plots (PSPs) established in Larix olgensis plantations located in the Jingouling forest farm of Wangqing Forest Bureau of northeast China ( Figure 1 and Table 1). Sample plot sizes vary from 775 m 2 to 2500 m 2 and represent a large within-sample plot variation of the height-diameter relationship. Twenty sample plots were nested within the five blocks (Block) with different site qualities, and each block contained four sample plots (Plot). PSPs were established in years between 1961 to 1964. Most of the sample plots were in the forest stands with different site conditions (   Table 1. A total of 765 Larix olgensis trees in 20 PSPs with 14 types of M * S and 2 altitude levels. Number of Larix olgensis trees (Numbers), stand density class (M), and site index class (S) of each sample plot in model fitting data and model testing data can be seen in this table. (Note: A = altitude of the sample plot; A = 1 when the altitude of the sample plot is higher than 700 m; A = 2 when the altitude of the sample plot is less than 700 m. M= stand density class of the sample plot; S= stand index class of the sample plot. The ranges of M and S are shown in Table 2 1  301  3  2  30  301  3  2  15  1  302  6  1  30  302  6  1  15  1  303  4  2  30  303  4 2 15 A total of 765 observations of height-diameter pairs of Larix olgensis were used in this study. Each observation point contained five variables, such as D, H, altitude (A), stand density class (M), and site index class (S). The environmental factors, such as A, M, and S are categorical variables. As altitude ranges from 600 m to 780 m, we divided this into two levels: >700 m and <700 m, denoted as level 1 and level 2, respectively. Table 1. A total of 765 Larix olgensis trees in 20 PSPs with 14 types of M * S and 2 altitude levels. Number of Larix olgensis trees (Numbers), stand density class (M), and site index class (S) of each sample plot in model fitting data and model testing data can be seen in this table. (Note: A = altitude of the sample plot; A = 1 when the altitude of the sample plot is higher than 700 m; A = 2 when the altitude of the sample plot is less than 700 m. M = stand density class of the sample plot; S = stand index class of the sample plot. The ranges of M and S are shown in Table 3 6  3  15  1  308  1  5  30  308  1  5  15  2  309  5  5  30  309  5  5  15  2  310,311,318  2  3  30  310,311,318  2  3  15  2  312  3  4  30  312  3  4  15  2  313  1  4  30  313  1  4  15  2  314  2  1  30  314  2  1  15  2  315  6  1  30  315  6  1  15  2  316  4  1  30  316  4  1  15  2  317  5  2  30  317  5  2  15  2  319  2  4  30  319  2  4  15  2  320  6  3  30  320  6  3  15 Different stand density levels significantly differently affect the tree growth [25][26][27]. Thus, it was necessary to divide stand density into 6 classes (denoted as 1, 2, . . . , 6) called stand density classes (M) according to the interval of 150 trees ha −1 . The effect of Meyer's site index (Meyer's SI) on tree height was analyzed by Du [28]. The Meyer's SI was estimated by equation supplied by Du [28]. The site index with 15.5 m as the starting point and 0.5 m as the interval was divided into 5 classes (denoted as 1, 2, . . . , 5), termed site index class (S). Out of 765 Larix olgensis trees, data from 510 trees was used as model fitting data and data from 255 trees used as model testing data. Data partition was done randomly; we randomly split the model fitting data and model testing data from the sample plots with same M * S according to the proportion 2:1. Summary of statistics is presented in Table 2, and height-diameter relationship plotted against D is shown in Figure 2.  [28]. The site index with 15.5 m as the starting point and 0.5 m as the interval was divided into 5 classes (denoted as 1, 2, ..., 5), termed site index class (S). Out of 765 Larix olgensis trees, data from 510 trees was used as model fitting data and data from 255 trees used as model testing data. Data partition was done randomly; we randomly split the model fitting data and model testing data from the sample plots with same M * S according to the proportion 2:1. Summary of statistics is presented in Table 3, and height-diameter relationship plotted against D is shown in Figure 2.

Base Models
Since Figure 2 shows a strong nonlinear relationship of height with diameter, we considered nine basic versatile nonlinear functions (Table 4) to describe such a relationship. We fitted these functions to the entire data set and selected the best performing one according to the AIC (Akaike information criterion, Equation (1) [29] and BIC (Bayesian information criteria, Equation (2) [30]. As the BIC is slightly more stringent for the number of parameters, we prefer choosing the model with the lowest BIC as the best model. We used the R-nls function to estimate the basic height-diameter models. The best fitting model was then expanded with introduction of the interactive effects of stand density and site index, and the sample plot-level random effects.
where k is the number of model parameters, n is the number of samples, L is the likelihood function value. Table 4. Candidate base models we considered.

The NLME Models
For the parameters with fixed effects in the nonlinear mixed-effects model, the most important thing is to determine what random effects each parameter needs to include. There are two ways to achieve this [40]. One method is to add all random effects for each parameter with AIC and BIC as main criteria to evaluate the fitting performance. Another method is to judge whether the mixed-effects model is properly parameterized based on the correlation between the estimated random effects. In this paper, we used the former method to choose the random effects for each parameter. There were six combinations of the random factors M, S, and M * S for each parameter. However, we excluded the random factor M + S because model did not converge when we added this to the model.

Parameter Estimation
The parameters of the NLME models were estimated by "nonlinear mixed-effects" module in Forstat2.2 [23]. A general NLME model was defined as: where φ i is formal parameter vector and includes the fixed effect parameter vector β and random effect parameter vector u i of the ith sample plot; symbols A i and Z i are the design matrices for β and u i , respectively. H ij and x ij are total height and the predictor vector of the jth tree on the ith sample plot, respectively. The estimated random effect parameter vector u i would be:û whereψ is the estimated variance-covariance for the random effects,R i is the estimated variance-covariance for the error term within the sample plot i. In this study, no structure covariance type BD (b) [41] was selected as the covariance type of ψ, and R(ψ = L T L, L is an upper triangular matrix). We assumed that the variances of random effects produced by structural variables were independent equal variances and there was no heteroscedasticity in our model; therefore, variance-covariance of sample plot i isR i = σ 2 I(σ 2 is the variance of the residual; I is the identity matrix.). The value of variance matrixψ or covariance matrixR i was calculated by restricted maximum likelihood with the sequential quadratic algorithm [21]. The f (·) is an interactive NLME model, andẐ i is an estimated design matrix:Ẑ where x i is a vector of the predictor on the sample plot i.

Model Evaluation
We used five statistical indicators to evaluate the performance of the interactive NLME height-diameter models such as MPSE,RMSE, and R 2 calculated by Equations (6)-(9) using both the model fitting and model testing data sets.
where H k andĤ k are the observed and estimated values of total height of the kth tree, respectively; H is the average value of the observed tree height; n is the total number of trees; MPSE is the mean percent standard error and it is a precision index reflecting the estimated value of individual tree height; RMSE is root mean squared error and it is used to measure the deviation between the estimated value and the observed value; and R 2 was used as a main criterion for model evaluation.
All the random effects of site index and stand density of the interactive NLME heightdiameter model were estimated using the Forstat software of the 2.2 version. All the interactive random effects of site index and stand density for each sample plot were used to evaluate the performance of the interactive NLME model using both the model fitting and model testing data sets.

Base Models
Nine basic diameter-height models were used as candidate base models ( Table 4). The entire data of 765 larch plants were used to fit the base models. According to the fit statistics produced by fitting nine candidate models (Table 5), the BIC of M2 was the lowest and the AIC of M2 was in the middle; therefore, M2 in Table 4 (Equation (10) was selected as the best base model in this study.
where H is the total height of Larix olgensis, D is the diameter at breast height of Larix olgensis, φ 1 and φ 2 are the formal parameters of this model. ε (ij)k is the error term of kth tree on the sample plot (ij).

The NLME Models
A total of 36 combinations of the random effects (Table 6) were derived from the formal parameters (φ 1 and φ 2 of the base model Equation (10)   were applied using model fitting data. The best performing combination was then selected based on the AIC and BIC scores.  Table 6. Thirty-six alternatives of the random effect constructions in the nonlinear mixed-effects models. Both the AIC and BIC of the interactive NLME model 15 are the smallest (Table 7). Therefore, the random effects constructions M * S and M * S were selected as the random variables on φ 1 and φ 2 , respectively. The AIC and BIC of the single effect model (model 1,2,3, and 4) were lower than those of the interactive effects model (model 15), which indicated that a substantial proportion of the tree height variations was better explained by the crossed random effects of the stand density and site index than the random effects of the stand density alone or a single site index alone.

The Interactive NLME Height-Diameter Model
Model 15 (Table 7) with the random construction variables [M * S, M * S] that had the best performance (AIC = 1771.64, BIC = 1797.02) was chosen as a final model. Model 15 with the random effects were the interaction effects of M and S and was defined as the interactive NLME height-diameter model in this article. The expression of the best interactive NLEM height-diameter model for further analysis is shown in Equation (11).
where M 1 (in this study, M 1 = 6) is the total classes of stand density, M 2 (in this study, M 2 = 5) is the total classes of site index, n ij is the number of observation points contained in the ith stand density class and the jth site index class. (ij) is the sample plot with the ith stand density class and the jth site index class. The ε (ij)k is the error term of the kth tree in sample plot (ij), which we assumed as R (ij) = σ 2 I (σ 2 (σ 2 > 0) is the variance of the residual). H is the tree height measurement value.

Parameter Estimates
All the parameter estimates of the basic model and interactive NLME model were significant (p < 0.05). The parameter estimates of the basic model obtained by an ordinary nonlinear least square (ONLS) function are as Equation (12), which we termed as the NLS model.Ĥ where D is the measured value of DBH,Ĥ is the estimated total tree height with the NLS model, ε ∼ N(0, 1.952I). The interactive NLEM height-diameter model with the estimated parameters obtained with the linearization approximation-sequential quadratic algorithm implemented in the "nonlinear nixed-effects" module of the Forstat software of the 2.2 version is given by: with where (ij) is the sample plot with the ith stand density class and the jth site index class, and k is the kth observation on the (ij) sample plot. D (ij)k is the measured value of the DBH of tree k on (ij) sample plot, whileĤ (ij)k is the estimated H of tree k on (ij) sample plot with the interactive NLME model. The fixed effect parameter estimates for the interactive NLME model are: β 1 = 3.1138, β 2 =−6.8180. The fixed effect parameter estimates of the NLME model were different from the parameter estimates of the NLS model (Equation (12), which also showed that different stand density classes and site index classes resulted in the different parameter estimates of the height-diameter model.
In addition, in order to compare the performance of the single-level mixed-effects model and the interactive mixed-effects model, we established a single-level mixed-effects model based on Model 2 with the best performance (with the smallest AIC and BIC) among the four single-level NLME models ( Table 7). The single-level NLME height-diameter model with the estimated parameters obtained with the linearization approximationsequential quadratic algorithm implemented in the "nonlinear mixed-effects" module of the Forstat software of the 2.2 version is given by: with where i is the sample plot with the ith stand density class and the j is the sample plot with the jth site index class, and k is the kth observation on the sample plot with the ith stand density class and the jth site index class. D (ij)k is the measured DBH of tree k on the (ij) sample plot, whileĤ (ij)k is the estimated H of tree k on the (ij) sample plot with a single-level NLME model.

Random Effects of the Interactive NLME Height-Diameter Model
The interactive random effects estimates obtained by fitting the NLME height-diameter model using the model fitting data are presented in Table 8. For two sample plots with the same M (stand density class) and altitude but different S (site index class), the estimated value of u1 and u2 of the interactive NLME height-diameter model Equation (13) are different. In addition, as for two sample plots with the same S and A but different M, the value of u1 and u2 of the interactive NLME height-diameter model Equation (13) are also different. This indicated that the random effect of the stand density on the Larix olgensis tree height was different at the different classes of site index, and the random effect of the site index on the Larix olgensis tree height was different at the different classes of stand density. Hence, there were strong interaction effects of the stand density and site index on the height-diameter relationship. Note: M * S means the type of stand density and site index for a sample plot (e.g., M * S = 3 * 2 means a sample plot with stand density class = 3 and site index class = 2). The u1 and u2 are the random effects on the parameters φ 1 and φ 2 of the interactive NLME height-diameter model (Equation (13)), respectively.

Model Evaluation
We used the likelihood ratio test to test the significant difference between the different models. Specifically, we compared a complex model with a simple model to test whether the complex model can significantly fit our data. If the difference between them is significant, the complex model could be used in future data analysis. The formula for the likelihood ratio test is: where LR is the likelihood ratio, L1 is the maximum likelihood value of the complex model, L2 is the maximum likelihood value of the simple model, d f 1 is the degree of freedom of the complex model, and d f 2 is the freedom of the simple model. The likelihood ratio test results are shown in Table 9 The difference between the interactive NLME model and the NLS model was significant. The difference between the interactive NLME model and the single-level NLME model was also significant. Therefore, the interactive NLME model we developed could be used for further data analysis. Table 9. Results of likelihood ratio test. LR is the likelihood ratio.

Complex Model vs. Simple Model LR p-Value
Interactive NLMEM (Equation (13)) vs. NLS (Equation (12) The estimated random effects of the interactive NLME height-diameter model (Equation (13)) in Table 8 were used for further evaluation. The height-diameter curves produced corresponding to different M * S are shown in Figure 3. We compared the statistical indicators, such as MPSE, RMSE and R 2 (Equations (6)-(9)), of three height-diameter models (Equations (12)-(14)) ( Table 10). Regardless of whether the model was used to predict the model testing data or model fitting data, the indicator values of the interactive NLME model were lower than those obtained for the NLS model and also lower than those of the single-level NLME model, which indicated that the crossed random effect of the stand density and site index could significantly improve the prediction accuracy of the model. Figure 4 shows the random distributions of the residuals produced with the three models.    (12), the single-level NLME model (Equation (14)), and the inter-   (12)), which is a non-random effect or fixed effect model, the single-level NLME model (Equation (14)), and the interactive NLME model (Equation (13)).  Figure 3. Height-diameter curves produced with Equation (13) corresponding to 14 combinations of M * S (stand density * site index; 3 * 2= stand density is class 3 and site index is class 2). The left is the height-diameter curves superimposed on the model fitting data, and the right is the height-diameter curves superimposed on the model testing data.  (12), the single-level NLME model (Equation (14)), and the interactive NLME model (Equation (13) for both the model fitting data and model testing data, respectively.  (12)), the single-level NLME model (Equation (14)), and the interactive NLME model (Equation (13)) for both the model fitting data and model testing data, respectively.

Discussion
We evaluated nine basic nonlinear functions of different forms proposed in the previous height-diameter modeling studies to identify the most suitable for our data. Based on the best-fitted base model (M2 -Table 4), the interactive NLME height-diameter model for Larix olgensis Henry was established through the integration of the variables describing the effects of the stand density and site quality on the height-diameter relationship. This expanded model described significantly larger variations of the height-diameter relationship than its base model counterpart did. When we tested the effect of the model with the model testing data, the R 2 of the interactive NLME height-diameter model increased by 5.2% compared with that of the base model. The interaction random effects of the stand density and site index were thus able to improve the prediction accuracy of the height-diameter model for Larix olgensis.
The height-diameter relationship can be influenced by the stand density [5][6][7]42], site index [8,13], and the interaction effect of the stand density and site index [14]. In our study, we assumed that the substantial proportion of the tree height variations was better explained by the interaction effect of the stand density and site index than the random effects of single stand density or single site index. This may be because there are interaction effects between the stand density and site index. In the previous model establishment [8,9], it is usually assumed that the stand density and site index are independent of each other; however, Ritchie [43] documented that the stand density affected the dominant height (the site index is estimated using dominant height and stand age). Under the assumption that the estimated value of the site index is indeed affected by density management, Ritchie [43] proposed a conjecture that, if the site index is used to drive the growth model, it may lead to the underestimation of the height growth. This conjecture was proved to be correct by our study as the AIC and BIC scores of the interactive NLME model (model 15 in Table 7) were higher than the single variable-based random effect models (model 1, 2, 3, and 4 in Table 7).
In addition, the random effect of the stand density on the tree height is different at different classes of the site index, and the random effect of the site index on the tree height is different at different classes of the stand density (Table 8). Hence, the data we analyzed further confirmed that there are significant interaction effects of stand density and site index on the height-diameter relationship. Moreover, in the low-density forest stands (in the forest stand with M = 1 and M = 2, respectively (Figure 5a), trees mainly grew radially and the DBH of the tree was larger with the increased site index of the stand. In the early development of the medium-density stands (in the forest stand with M = 3 and M = 4, respectively (Figure 5b), the higher the site index of the stand was, the larger the DBH of the tree in the stand was. After the canopy of the forest was closed, trees began to grow vertically due to resources constraint. It can be seen from the second half of Figure 5b that the height of the tree was taller with the increased site index of the stand. In high-density forest stands (in the forest stand with M = 5 and M = 6, respectively (Figure 5c), tree growth occurred vertically due to the fierce competition among the trees and resource constraints, so the higher the site index of the stand forest was, the taller the tree in the stand was. Our opinions were consistent with Xiang's conclusions [14] that they thought that the larger the initial planting density of the forest, the slower the growth of the DBH, and the better the site quality is, the faster the trees grow.
The interaction effects of the stand density and site index can significantly improve the prediction accuracy of the height-diameter model as the predicted performance of interactive NLEM height-diameter model (Equation (13) with R 2 = 0.6234) was better than the NLS model (Equation (12) with R 2 = 0.5717). In practical application, it is usually necessary to consider the impact of the interaction effects between the predictors on the estimated variable. Therefore, researchers in other disciplines can refer to the interactive NLME height-diameter model presented in this article to build the interactive NLME model. However, it can be seen from the residual plots (Figure 4) of the interactive NLME height-diameter model that, compared with large-diameter trees, the height of smalldiameter trees varies largely, which may be because the trees first focus on the height growth to win glory and then use their resources to increase the diameter [7]. To a certain extent, the altitude may affect the impact of the interaction effects between the stand density and site index on tree height, which requires further experiments to prove. In our data, the fitted height curves of Larix olgensis were obviously different at different altitudes ( Figure 6), and this indicates that the altitude could affect the parameter estimates of the height-diameter model to a certain extent. In order to further improve the accuracy of the height-diameter model of Larix olgensis, based on the interactive NLME height-diameter model presented in this article, an interactive NLME height-diameter model considering the group variable altitude could be further developed. The interaction effects of the stand density and site index can significantly improve the prediction accuracy of the height-diameter model as the predicted performance of interactive NLEM height-diameter model (Equation (13) with 2 R = 0.6234) was better than the NLS model (Equation (12) with 2 R = 0.5717). In practical application, it is usually necessary to consider the impact of the interaction effects between the predictors on the estimated variable. Therefore, researchers in other disciplines can refer to the interactive NLME height-diameter model presented in this article to build the interactive NLME model.
However, it can be seen from the residual plots (Figure 4) of the interactive NLME height-diameter model that, compared with large-diameter trees, the height of small-diameter trees varies largely, which may be because the trees first focus on the height growth to win glory and then use their resources to increase the diameter [7]. To a certain extent, the altitude may affect the impact of the interaction effects between the stand density and site index on tree height, which requires further experiments to prove. In our data, the fitted height curves of Larix olgensis were obviously different at different altitudes ( Figure 6), and this indicates that the altitude could affect the parameter estimates of the height-diameter model to a certain extent. In order to further improve the accuracy of the height-diameter model of Larix olgensis, based on the interactive NLME height-diameter model presented in this article, an interactive NLME height-diameter model considering the group variable altitude could be further developed.
In addition, there should be 30 types of M * S in the cross combination of six stand density levels and five site index levels, but our data only contained 14 types of M * S. Therefore, we planned the study to validate and verify the current results with more sample plots with extra M * S types and an increased amount of data in the future.  (13) corresponding to the different M * S superimposed on full data set. Plot (a) described the height-diameter relationship in low-density stands, plot (b) described the height-diameter relationship in medium-density stands, and plot (c) described the height-diameter relationship in high-density forest stands. (Symbol M * S means stand density * site index; 1 * 4= stand density of the stand is class 1 and site index of the stand is class 4). H= tree height (m); D = diameter at breast height (cm).

Conclusions
Through the height-diameter modeling with the two random components introduced in the model, it is concluded that there were strong interaction effects of the stand density and site index on the height-diameter relationship for Larix olgensis Henry. Furthermore, the random effects of the stand density on the Larix olgensis tree height were In addition, there should be 30 types of M * S in the cross combination of six stand density levels and five site index levels, but our data only contained 14 types of M * S. Therefore, we planned the study to validate and verify the current results with more sample plots with extra M * S types and an increased amount of data in the future.

Conclusions
Through the height-diameter modeling with the two random components introduced in the model, it is concluded that there were strong interaction effects of the stand density and site index on the height-diameter relationship for Larix olgensis Henry. Furthermore, the random effects of the stand density on the Larix olgensis tree height were substantively significant at the different classes of site index. In addition, the interaction between the stand density and site index significantly improved the prediction accuracy of the heightdiameter model when such an effect was included in the model. In practical forestry application, it is usually necessary to consider the impact of the interaction effects between the predictors on the estimated response variable. Researchers in other disciplines can refer to the interactive nonlinear mixed-effects height-diameter model presented in this article to build the interactive nonlinear mixed-effects models for other tree species.