Modeling Height–Diameter Relationships for Mixed-Species Plantations of Fraxinus mandshurica Rupr. and Larix olgensis Henry in Northeastern China

: The mixture of tree species has gradually become the focus of forest research, especially native species mixing. Mixed-species plantations of Manchurian ash ( Fraxinus mandshurica Rupr.) and Changbai larch ( Larix olgensis Henry) have successfully been cultivated in Northeast China. Height–diameter ( H–D ) models were found to be e ﬀ ective in designing the silvicultural planning for mixed-species plantations. Thus, this study aimed to develop a new system of H–D models for juvenile ash and larch mixed-species plantations, based on competition information and tree and stand attributes. The leave-one-out cross-validation was utilized for model validation. The result showed that the H–D relationship was a ﬀ ected not only by the tree attributes (i.e., tree size and competition information) but also by stand characteristics, such as site quality and species proportion of basal area. The best model explained more than 80% and 85% variation of the tree height of ash and larch, respectively. Moreover, model validation also conﬁrmed the high accuracy of the newly developed model’s predictions. We also found that, in terms of total tree height, ash in middle rows were higher than those in side rows, while larch in the middle rows were higher in the early growth period but then became lower than those in the side rows, as the diameter increased. The newly established H–D models would be useful for forestry inventory practice and have the potential to aid decisions in mixed-species plantations of ash and larch.


Introduction
One-third of the Earth's total surface area is covered by forests [1], which act as a major sink for carbon and has an essential role in maintaining the global climate system and carbon cycle [2]. Natural forests occupy 93% of the global forest area, with a declining annual rate of 6.5 million hectares per year (2010)(2011)(2012)(2013)(2014)(2015), while the other 7% were occupied by plantations with an average annual growth of 3.1 million hectares per year   [1,3]. Plantations are gradually expanding in many countries, including China. According to the ninth forest resource survey report (2014-2018) [4], the total forest area of China has reached 220 million hectares (22.96%), where as much as 79.54 million hectares were identified as plantations. Thus, strong theoretical knowledge and comprehensive forest management planning is necessary to manage such a large plantation area.

Study Sites
This study was conducted in six forest farms (the Maoer Mountain Forest Farm, Yimianpo Forest Farm, Xiaojiu Forest Farm, Shangzhi Forest Farm, Heilonggong Forest Farm, and the Yuanbao Forest Farm) in the Shangzhi city of the Heilongjiang province, Northeast China (from 127°17′ E to 129°12′ E and from 44°29′ N to 45°34′ N, Figure 1). The Shangzhi city has a total area of 8910 km 2 with varying landscapes, from high mountains in the east to undulating low hills in the west, and the altitude ranged from 116 to 1600 m above sea level. Dominated by dark brown soil, the region had a temperate continental monsoon climate characterized by less rain and low temperature in spring, warm and rainy summer, long and cold winter, and autumns that brought frequent early frosts and a steep decline in temperature. It had a low mean annual temperature of 2.5 °C, with the highest-and lowestmonthly temperatures being 21.6 °C in July and −20.2 °C in February, respectively. The Shangzhi City had an average annual rainfall of 652.2 mm, with single rain events ranging from 5 to 577 mm. The total area of the forest land was approximately 5100 km 2 , out of which nearly 1300 km 2 were plantations. The dominant tree species for the plantations were Larix olgensis Henry, Pinus koraiensis Siebold & Zucc., Pinus sylvestris var. mongholica Litv., and Fraxinus mandshurica Rupr. In the present study, the mixed-species plantations of Manchurian ash and Changbai larch in the Shangzhi city were relatively young; their stand age ranged from 10 to 25 years. They were established at an initial density of 3300 cutting·ha −1 , with reasonably narrow spacing (1.5 m × 2.0 m) and several variations of row-wise mingling patterns ( Figure 2). The planting designs were distinguished as seven categories: (A) 1-row larch: 1-row ash (1:1) (B) 2-rows larch: 2-rows ash (2:2) (C) 3-rows larch: 2-rows ash (3:2) In the present study, the mixed-species plantations of Manchurian ash and Changbai larch in the Shangzhi city were relatively young; their stand age ranged from 10 to 25 years. They were established at an initial density of 3300 cutting·ha −1 , with reasonably narrow spacing (1.5 m × 2.0 m) and several variations of row-wise mingling patterns ( Figure 2). The planting designs were distinguished as seven categories: (A) 1-row larch: 1-row ash (1:1) (B) 2-rows larch: 2-rows ash (2:2) Forests 2020, 11, 610 4 of 22 (C) 3-rows larch: 2-rows ash (3:2) (D) 3-rows larch: 3-rows ash (3:3) (E) 5-rows larch: 3-rows ash (5:3) (F) 5-rows larch: 5-rows ash (5:5) (G) 6-rows larch: 4-rows ash (6:4) Forests 2020, 11,610 4 of 21 (D) 3-rows larch: 3-rows ash (3:3) (E) 5-rows larch: 3-rows ash (5:3) (F) 5-rows larch: 5-rows ash (5:5) (G) 6-rows larch: 4-rows ash (6:4) Based on those various planting designs, both larch and ash were divided into two categories-(1) trees in the middle of the rows (denoted as TMR); and (2) trees on the sides of the rows (denoted as TSR). TMR represents the situation in which the neighboring rows were trees of the same species, while TSR refers to the condition in which the adjacent rows were trees of different species.   Based on those various planting designs, both larch and ash were divided into two categories-(1) trees in the middle of the rows (denoted as TMR); and (2) trees on the sides of the rows (denoted as TSR). TMR represents the situation in which the neighboring rows were trees of the same species, while TSR refers to the condition in which the adjacent rows were trees of different species.

Data
In July and August of the three consecutive years (2017, 2018, and 2019), 69 sample plots from 23 plantations (3 sample plots per plantation) were established, with a permanent length of 50 m and the width of at least three repeats, which varied according to the planting designs. The so-called "repeat" was determined as a basic unit for the width of each sample plot, and one repeat consisted of some rows in accordance with the summation of each planting design's ratio. For example, for 1:1 planting design, one repeat equaled to two rows containing an equivalent proportion of both species. Meanwhile, for 3:2 planting design, one repeat was equivalent to five rows comprising three rows of larch and two rows of ash. A similar mechanism was also performed for other planting designs. Further, several basic information for all individual trees located in the sample plots was recorded, such as tree species, status (dead or alive), and type of rows (TMR or TSR). Then, D and H were measured using diameter tapes and Vertex IV Ultrasonic Hypsometer made by Haglöf Sweden, respectively.
After the basic measurement was conducted, the following stand variables were measured-(1) regardless of the tree species, such as the number of trees per hectare (N, trees·ha −1 ), quadratic mean diameter (Dq, cm), the dominant height of 100 D-based largest trees per hectare (Hd, m), basal area (BA, m 2 ·ha −1 ), and the dominant diameter (Dd, cm) of 100 D-based largest trees per hectare; (2) by taking into account the tree species, such as the number of trees per hectare (N L for larch and N F for ash), quadratic mean diameter (Dq L for larch and Dq F for ash), basal area (BA L for larch and BA F for ash), dominant height (Hd L for larch and Hd F for ash), and dominant diameter (Dd L for larch and Dd F for ash). In order to describe the specific structure of the mixed-species stand, we also calculated other simple indices at stand level [44]-(1) the proportion of each species relative to total trees (RN L = N L /N for larch and RN F = N F /N for ash); and (2) the basal area proportion of each species relative to the total basal area (RBA L = BA L /BA for larch and RBA F = BA F /BA for ash). In addition, the information on the neighboring trees represented the competitive state of the selected tree, which had an important effect on the H-D relationship. Hence, we calculated the distance-independent indices-the ratio of each individual D to Dd (RD) [45] and the basal area of the trees larger than the subject tree (BAL, m 2 ·ha −1 ). Additionally, the corresponding intra-and inter-specific indices were also calculated: RD intra , RD inter , BAL intra , and BAL inter . Summary statistics of the tree and stand variables are presented in Table 1. S.D.: Standard Deviation. N = number of trees per hectare; N L = number of trees per hectare for larch; N F =number of trees per hectare for ash; Dq = quadratic mean diameter; Dq L =quadratic mean diameter for larch; Dq F =quadratic mean diameter for ash; Hd = dominant height; Hd L = dominant height for larch; Hd F = dominant height for ash; Dd = dominant diameter; Dd L = dominant diameter for larch; Dd F = dominant diameter for ash; BA = basal area; BA L = basal area for larch; BA F = basal area for ash; RN L = relative N L for larch; RN F = relative N F for ash; RBA L = relative basal area for larch; RBA F = relative to basal area for ash.

The Comparison of Total Tree Height (H) Between TMR and TSR
Previous studies revealed that in ash and larch mixed-species plantations, H of TMR and TSR were significantly different under the same specific stand age [27,28]. Prior to the development of H-D models, we utilized the probability density curves of H in different rows, to compare the differences of H between TMR and TSR. Furthermore, the differences were examined statistically by using analysis of variance (ANOVA) or nonparametric tests, at a significance level of 0.05. The choice of using either ANOVA or nonparametric tests, depended on the normality test of each sample. If the normality test was successfully passed, ANOVA would be used; otherwise, nonparametric tests would be more appropriate.

Generalized H-D Model
The H-D relationship was strongly influenced by the stand characteristics and neighboring trees' information [32]. Thus, a generalized H-D (GHD) model was necessary. Once the best H-D model was selected, we evaluated the potential of several stand variables and competition indices contributions to the H-D relationship. First, the selected basic model was fitted to the data for each plot to obtain the corresponding parameter estimation values. Furthermore, the relationships between model coefficients and stand variables or competition indices were scrutinized by scatter plot and correlation analysis [58]. The stand and competition variables needed to be significant, and the fitting statistics of the GHD model needed to be improved.
Dummy variable models are frequently used in forestry modeling, to describe the effect of species-, region-, and even origin-specific differences [40,59,60]. In this study, a dummy variable ( ) was included in the GHD model to describe the row-specific effect ( = 0 for TMR and = 1 for TSR), which was called a GHDR model. An F-test was then used to test whether separate models were required for the different type of rows, from a statistical point of view [59]. Parameters of the GHD model were expressed by a linear function , which was defined as follows: where is the parameter in the GHD models; βj and βk are the estimated parameters in the GHDR model.

Functions Number and Forms
Reference Richards, 1959 Curtis, 1967 Bates, 1980 Wykoff, 1982 Schumacher, 1939 Ratkowsky, 1990 Huang, 1992 The H-D relationship was strongly influenced by the stand characteristics and neighboring trees' information [32]. Thus, a generalized H-D (GHD) model was necessary. Once the best H-D model was selected, we evaluated the potential of several stand variables and competition indices contributions to the H-D relationship. First, the selected basic model was fitted to the data for each plot to obtain the corresponding parameter estimation values. Furthermore, the relationships between model coefficients and stand variables or competition indices were scrutinized by scatter plot and correlation analysis [58]. The stand and competition variables needed to be significant, and the fitting statistics of the GHD model needed to be improved.
Dummy variable models are frequently used in forestry modeling, to describe the effect of species-, region-, and even origin-specific differences [40,59,60]. In this study, a dummy variable (S 1 ) was included in the GHD model to describe the row-specific effect (S 1 = 0 for TMR and S 1 = 1 for TSR), which was called a GHD R model. An F-test was then used to test whether separate models were required for the different type of rows, from a statistical point of view [59]. Parameters of the GHD model were expressed by a linear function S 1 , which was defined as follows: where β i is the parameter in the GHD models; β j and β k are the estimated parameters in the GHD R model.

Mixed-Effects Model
The mixed-effects model is a model that includes both fixed-and random-effect variables. The application of a mixed-effect model has increased in the field of forest modeling [40,41]. In this study, a nonlinear mixed-effects (NLME) modeling framework was used for the hierarchical structure of H-D data, by setting the sample plots as random effects in HD B , GHD, and GHD R . Numerous studies have applied the mixed-effects models to describe H-D relationships and have improved the model fitting and prediction accuracy [32,43]. At a plot-level, the NLME model of the jth tree height in the ith sample plot was modeled as [61]: where H ij represents the height of the jth tree in the ith plot; m and n i are the number of plots and observations in the ith plot, respectively; f is a nonlinear function of a plot-specific parameter vector φ ij and a covariate vector ν ij ; and ε ij is a within-plot error term, which subjects to a multivariate normal distribution with a mean value vector of 0 and variance-covariance matrix of R. φ ij is given as: where β is the fixed-effect parameter vector; b i is the random-effect parameter vector of the ith plot, which was assumed to have a multivariate normal distribution with a mean value vector of 0 and variance-covariance matrix of G, A ij , and B ij are the incidence matrices of the appropriate dimensions, consisting of 0 or 1. In the weighted regression analysis, a power function was used to correct the heteroscedasticity of variance exhibited in the H-D models, and to find the most effective variance functions [58]. All regressions were calculated in the R software 3.5.3., by using both base and nlme packages [57,62].

Model Assessment and Validation
In the present study, six types of models were evaluated-(1) HD B model; (2) GHD model; (3) GHD R model; (4) mixed-effect HD B model (HD BM ); (5) mixed-effect GHD model (GHD M ); and (6) mixed-effect GHD R model (GHD RM ). These models were validated using the leave-one-out cross-validation [63], in which the H-D models were constructed with all-but-one sampled plot data (fitting data included m-1 sample plots) and the fitted model was used to predict the excluded sampled plot. Statistical indices were used to assess and compare model performance, and the best model was selected, based on the following criteria: Statistical significance of the estimated parameters; Statistics of model assessment, such as the adjusted coefficient of determination (R a 2 ), root mean square error (RMSE), and Akaike information criterion (AIC) [64]. Performance of model validation, such as mean absolute error (MAE) and mean absolute error percent (MAE%) [65].
The marginal adjusted coefficient of determination (R am 2 ) only involved fixed effect, while the conditional adjusted coefficient of determination (R ac 2 ) concerned both fixed-and random-effects, which was further used to evaluate the mixed-effects models [40]. However, the mixed-effects model validation needed prior information to calculate the random-effects parameters. The parameters of random effects were calculated using the best linear unbiased predictions (BLUPs) method [61]. A vector of random effects parameters of sampled plot k was calculated with Equation (4).
Forests 2020, 11, 610 where b k is a vector of random effects parameters of sampled plot k calculated by all sample trees; G is the variance-covariance matrix estimated in the modeling process; R k is the corresponding variance-covariance matrix of within-group errors; Z k is the matrix of the partial derivatives of the nonlinear function, with respect to its random parameters; and ε k is error terms of the sampled plot k using the fixed effects parameters of the mixed-effects models.

Comparison of Sample Designs
Although validation of the mixed-effects models used all sample tree height observations to calculate the random effect parameters, different numbers of available measured trees were then used to compare the prediction performance. In calculating the random effect parameters, researches have pointed out that the increasing number of sampled trees per plot would decrease the bias between the observations and the predicted values [64]. However, a large number of sample trees per plot would lead to a higher inventory cost, which only gives a slight improvement in the model's accuracy. Thus, four trees per sample plot seemed to be adequate for the calculation of random effects [41]. Simulations were performed using subsamples of tree height measurements, where predefined sample sizes of 2, 3, 4, 5, 6,7,8,10,12,14,16,18,20,24, and 26 trees were subsampled without replacement and the calculations were repeated continuously for 100 times. Within each subsample, we estimated the random effect parameter values and obtained the final prediction bias. These simulations were also performed in R software 3.5.3. using the nlme package [62].

The Variation of H in Different Rows
The probability density plots of the two species were used to analyze the distribution of H across different statuses (TMR and TSR, Figure 4). Preliminary visual inspection showed that the two species' heights were obviously different across both TMR and TSR, in which, the data centered on different values. As the non-normality (p < 0.01) of H in different rows of two species showed in Figure 4, a nonparametric test (Mann-Whitney test [66]) was applied to compare whether the value of the average height of TMR and TSR were significantly different between larch and ash. The result of the two-way Mann-Whitney test supported that the average height of TMR and TSR were significantly different (p < 0.01). However, limited evidence confirmed the different height of TMR and TSR by the data collected so far; hence, it was also necessary to use a modeling method for further analysis and verification. Although validation of the mixed-effects models used all sample tree height observations to calculate the random effect parameters, different numbers of available measured trees were then used to compare the prediction performance. In calculating the random effect parameters, researches have pointed out that the increasing number of sampled trees per plot would decrease the bias between the observations and the predicted values [64]. However, a large number of sample trees per plot would lead to a higher inventory cost, which only gives a slight improvement in the model's accuracy. Thus, four trees per sample plot seemed to be adequate for the calculation of random effects [41]. Simulations were performed using subsamples of tree height measurements, where predefined sample sizes of 2, 3, 4, 5, 6,7,8,10,12,14,16,18,20,24, and 26 trees were subsampled without replacement and the calculations were repeated continuously for 100 times. Within each subsample, we estimated the random effect parameter values and obtained the final prediction bias. These simulations were also performed in R software 3.5.3. using the nlme package [62].

The Variation of H in Different Rows
The probability density plots of the two species were used to analyze the distribution of H across different statuses (TMR and TSR, Figure 4). Preliminary visual inspection showed that the two species' heights were obviously different across both TMR and TSR, in which, the data centered on different values. As the non-normality (p < 0.01) of H in different rows of two species showed in Figure 4, a nonparametric test (Mann-Whitney test [66]) was applied to compare whether the value of the average height of TMR and TSR were significantly different between larch and ash. The result of the two-way Mann-Whitney test supported that the average height of TMR and TSR were significantly different (p < 0.01). However, limited evidence confirmed the different height of TMR and TSR by the data collected so far; hence, it was also necessary to use a modeling method for further analysis and verification.

Basic H-D Model Results
By fitting the 12 candidate functions listed in Table 2, Ra 2 ranged from 0.5780 to 0.6011 for ash and from 0.6295 to 0.6513 for larch (Table 3), which showed a relatively slight fluctuation in a small range. The other two fitting statistics (RMSE and AIC) also displayed a similar pattern. Table 3 indicated that the number 10th function from Ratkowsky was found to be better than others for both

Basic H-D Model Results
By fitting the 12 candidate functions listed in Table 2, R a 2 ranged from 0.5780 to 0.6011 for ash and from 0.6295 to 0.6513 for larch (Table 3), which showed a relatively slight fluctuation in a small range. The other two fitting statistics (RMSE and AIC) also displayed a similar pattern. Table 3 indicated that the number 10th function from Ratkowsky was found to be better than others for both larch and ash, with a higher value of R a 2 and a lower value of RMSE and AIC. The function in bold is the best basic model. R a 2 = adjusted coefficient of determination; RMSE = root mean square error; AIC = Akaike information criterion.

Generalized H-D Model
Correlation analysis of each parameter in the HD B model indicated that the dominant height considering species and species proportion by basal area had a relatively greater impact on the parameters. Different competition indices had different effects on species; therefore, the final forms of the GHD model were: Ash: Larch: There was a substantial improvement in model fits when the covariates were added, (Tables 4 and 5). The R a 2 increased from 0.6011 to 0.7800 for ash and from 0.6513 to 0.8418 for larch. The corresponding RMSE and AIC also showed a significant reduction. Stand variables and competition indices explained more variations on the basis of HD B . Then, a dummy variable S 1 was added into the different parameters, and the fitting performance was evaluated by a smaller AIC value and the statistical significance of parameters. The results of the F-test between the GHD and GHD R models turned out to be statistically significant for both ash and larch (p < 0.01). The final forms of the GHD R models for both species were as follows:    Ash: Larch: Based on the GHD R models' formula, the effects of S 1 (TMR and TSR) on the H-D relationship for larch and ash were simulated at the average values of other predictors with three random sampled plots ( Figure 5). The GHD R model results showed that there was a significant difference in H between TMR and TSR for both species. The H of TMR was found to be higher than those of TSR for both ash and larch. However, it is worth noting that for larch, the H of TMR was higher in the early growth period, but then the difference decreased with time, before finally being lower than the H of TSR.

Mixed-Effects H-D Model
The plot-level random effect was introduced in the HD B , GHD, and GHD R models; and was found to be significantly different among the three models, in which a variation of the H-D models was observed among the different plots. However, employing a plot-level random effect on both the GHD and GHD R models for the two species showed a slight improvement. The HD BM model described nearly 20% more variations than the HD B model in the H-D relationship, and there was an apparent difference between the fit performance evaluated by only considering the fixed-effect parameters and the mixed-effect parameters (R am 2 and R ac 2 were 0.5367 and 0.8073, respectively, for ash, and 0.5849 and 0.8590 for larch). The R am 2 values of GHD M and GHD RM models were slightly smaller compared to the corresponding R ac 2 values. The final form of the HD BM , GHD M , and GHD RM models were: HD BM models: Ash: Larch: GHD M models: Ash: Larch: GHD RM models: Ash: Larch: where u i1 and u i2 are the random parameters, while the other symbols are defined above. The parameter estimates, variances, and fit statistics for all H-D models are listed in Tables 4 and 5. The scatterplots between the standardized residuals and the estimated height by utilizing weighted regression are shown in Figure 6. Within the sampled plot, heteroscedasticity in the residuals was reduced by applying a weighted regression in the form of a power function, at the base of the predicted height. The values of the power variance functions are also listed in Tables 4 and 5.
The effects of Hd, RBA, and BAL considering tree species along with Rd on the H-D relationship were simulated with H-D curves by species-specific GHD RM model of ash (Figure 7) and larch (Figure 8). The variables of interest were roughly divided into four equal intervals, where other predictors were at the mean value of our data. The simulation showed that all variables, except Hd F , had a negative effect on H for both species, which was similar to others [40]. As shown in Figures 7 and 8, Hd gave a more considerable contribution to the H-D model, among all variables. The differences of H between RBA F and RBA L were larger as the D increased, and the species competition significantly affected the initial H-D relationship for ash and increased the influence of the competition to larch.
where and are the random parameters, while the other symbols are defined above. The parameter estimates, variances, and fit statistics for all H-D models are listed in Tables 4 and 5. The scatterplots between the standardized residuals and the estimated height by utilizing weighted regression are shown in Figure 6. Within the sampled plot, heteroscedasticity in the residuals was reduced by applying a weighted regression in the form of a power function, at the base of the predicted height. The values of the power variance functions are also listed in Tables 4 and 5.  Figure  8). The variables of interest were roughly divided into four equal intervals, where other predictors were at the mean value of our data. The simulation showed that all variables, except HdF, had a negative effect on H for both species, which was similar to others [40]. As shown in Figures 7 and 8, Hd gave a more considerable contribution to the H-D model, among all variables. The differences of H between RBAF and RBAL were larger as the D increased, and the species competition significantly affected the initial H-D relationship for ash and increased the influence of the competition to larch.   Table 4. Figure 7. The effects of the corresponding variables were simulated for ash. The H-D curves were generated using the fixed-effect parameter of the GHD RM model in Table 4. Figure 7. The effects of the corresponding variables were simulated for ash. The H-D curves were generated using the fixed-effect parameter of the GHDRM model in Table 4.  Table 5.

Model Validation
The GHDRM model was found to have the best performance among the six models, followed by GHDM, HDBM, GHDR, GHD, and HDB. The GHDRM model gained the lowest value of MAE (0.8095 m for ash and 0.8111 m for larch) and MAE% (6.4927% for ash and 6.7883% for larch) ( Table 6). The results showed that the fit statistics and the performance of model validation of all models did not have a significant increase (Tables 4 and 5). For larch, apart from the HDB model, the value of MAE ranged from 0.8095 to 0.8299, and MAE% decreased from 6.6421% to 6.4927%, and a similar regular pattern also appeared for ash. Overall, the prediction performances were within a reasonable range and improved as the model's complexity increased.  Table 5.

Model Validation
The GHD RM model was found to have the best performance among the six models, followed by GHD M , HD BM , GHD R , GHD, and HD B . The GHD RM model gained the lowest value of MAE (0.8095 m for ash and 0.8111 m for larch) and MAE% (6.4927% for ash and 6.7883% for larch) ( Table 6). The results showed that the fit statistics and the performance of model validation of all models did not have a significant increase (Tables 4 and 5). For larch, apart from the HD B model, the value of MAE ranged from 0.8095 to 0.8299, and MAE% decreased from 6.6421% to 6.4927%, and a similar regular pattern also appeared for ash. Overall, the prediction performances were within a reasonable range and improved as the model's complexity increased. The residuals of the six models were plotted across four D classes (Figure 9). The average values of MAE and MAE% were relatively stable in all D classes and six models, except the HD B models for both small size (5 < D < 10 cm) ash and larch, which indicated a stable prediction accuracy within the range of D. This could improve the prediction accuracy by applying both the join of covariates and NLME modeling approach for larch and ash. Compared to the ordinary least square, the NLME method showed a significant improvement from the perspective of MAE and MAE%.
The residuals of the six models were plotted across four D classes (Figure 9). The average values of MAE and MAE% were relatively stable in all D classes and six models, except the HDB models for both small size (5 < D < 10 cm) ash and larch, which indicated a stable prediction accuracy within the range of D. This could improve the prediction accuracy by applying both the join of covariates and NLME modeling approach for larch and ash. Compared to the ordinary least square, the NLME method showed a significant improvement from the perspective of MAE and MAE%.

Result of the Sampled Designs Comparison
The prediction performance of the GHDRM model for each species with different numbers of prior measured trees are shown in Figure 10. The value of MAE and MAE% were decreased with increasing numbers of sampled trees, which implicated that the prediction accuracy increased when more prior height information was used to calculate the random effect parameters. There was no significant difference between six, seven, or even more sampled trees, based on the significant test of MAE and MAE% in the different numbers of previously measured trees.

Result of the Sampled Designs Comparison
The prediction performance of the GHD RM model for each species with different numbers of prior measured trees are shown in Figure 10. The value of MAE and MAE% were decreased with increasing numbers of sampled trees, which implicated that the prediction accuracy increased when more prior height information was used to calculate the random effect parameters. There was no significant difference between six, seven, or even more sampled trees, based on the significant test of MAE and MAE% in the different numbers of previously measured trees.

Basic H-D Model and Generalized Model
In order to develop the H-D model for mixed-species plantations of ash and larch, 12 alternative models based on D as the only predictor were evaluated in this study. The function number 10 in Tables 3 and 4, which were previously used to predict tree height in British Columbia [67], gave the best performance to be fitted in our data. In this function, was considered to be the asymptotic height, and represented steepness, known by its derivative function that would be extended to the formula of stand variables and competition information, to develop a generalized model.
Hd mostly represented site quality and was frequently used in the H-D models for both monocultures and mixed-species plantations [32,40]. Species proportion by basal area is an important attribute in mixed-species plantations [44], which had an effect on stand structure and tree growth.

Basic H-D Model and Generalized Model
In order to develop the H-D model for mixed-species plantations of ash and larch, 12 alternative models based on D as the only predictor were evaluated in this study. The function number 10 in Tables 3 and 4, which were previously used to predict tree height in British Columbia [67], gave the best performance to be fitted in our data. In this function, β 0 was considered to be the asymptotic height, and β 1 represented steepness, known by its derivative function that would be extended to the formula of stand variables and competition information, to develop a generalized model.
H d mostly represented site quality and was frequently used in the H-D models for both monocultures and mixed-species plantations [32,40]. Species proportion by basal area is an important attribute in mixed-species plantations [44], which had an effect on stand structure and tree growth. Beside stand variables, competition information expressed by some indices, including distance-independent (i.e., BAL) [32,68], is necessary for the H-D model. The distance-dependent competition indices might show a better performance in describing the competition status. However, it is time-consuming and labor-intensive, since each neighboring tree's coordinates need to be measured precisely. Due to the limitation of both fieldwork duration and economic resources, we chose the indices of BAL and RD to explain the competition variations in the H-D models. When these covariates were added into the basic H-D models of ash and larch, R ac 2 increased by 34.39% and 31.89%, while RMSE reduced by 30.94% and 36.41%, respectively. Hence, the extended basic H-D model was preferable, since it gave a relatively high prediction accuracy. The results of the model validation by utilizing the leave-one-out cross-validation also showed that the generalized H-D models were better than the basic models, in terms of prediction. Thereby, H of the mixed-species plantations of ash and larch would be accurately estimated with the newly developed models, which were found to be more flexible and effective. GHD R used a dummy variable to analyze the differences between TMR and TSR for both species, which proved to be significantly different by the F-test result between GHD R and GHD model (p < 0.01). Wang et al. [26] found that the mixture of ash and larch improved intraspecific competition of ash, which might be the reason why H of ash in middle rows was higher than those in the side rows, in three different random sample plots, for the same value of D ( Figure 5). In contrast to ash, H of larch was found to be higher in TMR than TSR in small D; then the H of larch was found to be higher in TSR than TMR as the D increased, which was similar to others [27,28]. Larch is classified as shade intolerant-species. In mixed-plantations of larch and ash, larch in TSR was often shaded by ash, primarily the middle and lower part of their crown. As a result, larch in TMR got more lights and grew faster, mainly in the early phase after planting [69]. In the present study, we found that there was a large number of dead larches in mixed-species plantations, especially for those that grew in the side rows. One possible reason behind this phenomenon was the effect of interspecific relationship, since the remaining larches in TSR were found to grow healthier with broader space.

H-D Model with NLME
The NLME model considers the potential differences between plot-levels and is widely applied in H-D modeling [40,41,64]. The performances of the GHD, GHD R , and HD BM models were quite similar, which indicated that predicting a random effect can be thought of as an empirical surrogate for stand-level variables that affect the H-D relationship [68]. Generally, the convergence of mixed-effects models become harder to achieve as the number of random parameters increase. It was almost impossible to achieve the convergence when the models contained more than two random parameters, in which further improvement on the algorithm was required [70]. Actually, the mixed-effects model by only considering the fixed parameters could be used to predict height [40], but the prediction accuracy would be significantly lower and was not recommended to be applied. In our study, the R am 2 of HD BM , GHD M , and GHD RM were lower than the same models estimated by nonlinear least squares. Therefore, the predictions of the NLME models containing the random effects were better and had relatively high precision.
As shown in Table 6, the prediction ability of the mixed-effects model was better by using the random effect parameters calculated by all measurements of H. However, it was very difficult or even unrealistic to measure the H of all trees in forest resource inventory, which implied that choosing several trees for H measurement would be more favorable. Many studies have analyzed the impact of different numbers of available measured trees on the prediction of mixed-effects model [40,41,43].
The prediction performance of mixed-effects model validation with a priori samples of different sizes was compared, and the result showed that the model's performance increased as the prior information increased. Based on the significant test of MAE and MAE% in Figure 10, six trees would be a reasonable choice and close to the number of dominant trees per sample plot, as slightly different from other studies [41,43].

Influence on Species Mixing in H-D Relationship
One of the advantages of mixed-species over monoculture plantings is that it yields a higher productivity [7], and the mixed-species plantations of ash and larch in Northeast China is thought to be a more productive plantation than their corresponding monocultures [28,30]. Zhang et al. [28] concluded that H of both species in mixed-species plantation were higher than in monoculture plantations. The final H-D model for larch in this study differed from the model formula developed by Zang et al. [32] using data from monoculture plantations of larch across a broader geographic area; both included stand variables and competition indices that improved fitting performance and increased the precision of predictions. The predictions of H utilizing H-D model developed by Zang et al. [32] were significantly lower compared to both the observed and predicted H by using the model developed in this study ( Figure 11). Therefore, preliminary analysis showed that H of larch in mixed-species plantation of ash and larch were higher than that in monoculture plantations of larch. However, a further comparison of the H-D model using more detailed data between ash in mixed-species and monoculture plantations is essential, since an H-D model in monoculture plantations is still lacking.
Generally, the H-D models developed in this study show strong practicability and operability in predicting height by relying on diameter and dominant tree height data, which could be used to calculate other covariates. Therefore, the developed models in this research were useful for predicting tree height and in future work, it will be an essential tool for analyzing the differences in total biomass, stand volume, wood quality, etc. of mixed-species plantations and monoculture plantations.
parameters, in which further improvement on the algorithm was required [70]. Actually, the mixedeffects model by only considering the fixed parameters could be used to predict height [40], but the prediction accuracy would be significantly lower and was not recommended to be applied. In our study, the Ram 2 of HDBM, GHDM, and GHDRM were lower than the same models estimated by nonlinear least squares. Therefore, the predictions of the NLME models containing the random effects were better and had relatively high precision.
As shown in Table 6, the prediction ability of the mixed-effects model was better by using the random effect parameters calculated by all measurements of H. However, it was very difficult or even unrealistic to measure the H of all trees in forest resource inventory, which implied that choosing several trees for H measurement would be more favorable. Many studies have analyzed the impact of different numbers of available measured trees on the prediction of mixed-effects model [40,41,43]. The prediction performance of mixed-effects model validation with a priori samples of different sizes was compared, and the result showed that the model's performance increased as the prior information increased. Based on the significant test of MAE and MAE% in Figure 10, six trees would be a reasonable choice and close to the number of dominant trees per sample plot, as slightly different from other studies [41,43].

Influence on Species Mixing in H-D Relationship
One of the advantages of mixed-species over monoculture plantings is that it yields a higher productivity [7], and the mixed-species plantations of ash and larch in Northeast China is thought to be a more productive plantation than their corresponding monocultures [28,30]. Zhang et al. [28] concluded that H of both species in mixed-species plantation were higher than in monoculture plantations. The final H-D model for larch in this study differed from the model formula developed by Zang et al. [32] using data from monoculture plantations of larch across a broader geographic area; both included stand variables and competition indices that improved fitting performance and increased the precision of predictions. The predictions of H utilizing H-D model developed by Zang et al. [32] were significantly lower compared to both the observed and predicted H by using the model developed in this study ( Figure 11). Therefore, preliminary analysis showed that H of larch in mixedspecies plantation of ash and larch were higher than that in monoculture plantations of larch. However, a further comparison of the H-D model using more detailed data between ash in mixedspecies and monoculture plantations is essential, since an H-D model in monoculture plantations is still lacking.

Conclusions
Tree height has important implications on stem volume, biomass, carbon storage, and others, with a growing interest in mixed-species plantations. The main objective of the developed H-D models in this study was to predict the missing height data in mixed-species plantations, which might reduce field workload and improve the prediction accuracy. Compared to previous research work on H-D models, we considered the row-specific effect along with intra-and inter-specific competition effects to describe the H-D relationship in mixed-species plantation of ash and larch. The GHD model included D as a predictor, and two covariate predictors (stand variables and competition information), which showed a significant improvement as compared to the basic models. According to the actual condition, we could use the different H-D models to predict tree height. The three mixed-effects models HD BM , GHD M , and GHD RM performed well, but the use of GHD RM was preferable. However, the newly developed H-D models was not only important to predict the tree height but was also essential to comprehend the vertical structure of mixed-species plantation of ash and larch in Northeast China.
Author Contributions: L.X. collected the field data, performed data analysis, and wrote the paper; F.R.A.W. collected the field data and wrote the paper; L.D. and F.L. supervised and coordinated the research project, designed and installed the experiment, took some measurements, and contributed to writing the paper. All authors have read and agreed to the published version of the manuscript.