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

Abstract

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 fitted 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 significantly 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.

1. 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 height-diameter 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].

Stand density is an important factor affecting tree growth [4]. Low stand density may reduce the competition impact, increase soil moisture availability, and produce favorable physiological responses, thereby significantly affecting tree growth [5]. Various studies have shown a significant influence of the stand density on the height-diameter relationship. For example, Zeide and Vanderschaaf [6] documented that Pinus taeda L. in higher density stands were taller than those in lower density stands. Thanh and Dinh [7] documented that the height-diameter relationship of Pinus koraiensis appears to differ between stand density classes. In addition, the authors of these studies [8,9,10,11] used stand density as an independent variable to construct the height-diameter models and were able to significantly improve the prediction accuracy.

The height-diameter relationship varies from one site type to another, and, even within the same site type, the relationship would not remain the same over time [12]. Feldpausch and Banin [13] documented that these differences in site equality in different seasons would effectively amplify the dry season effect, and may cause the forests in the Amazon to have larger tree height and diameter. According to Krisnawati [8], the individual tree growth in a stand may take longer on a poorer site than in a site with a higher index class. In addition, the site index can improve the prediction accuracy of the height-diameter model when it is included as an independent variable; for example, Krisnawati [8] and Missanjo [9] included the site index as an independent variable to establish the height-diameter model.

Summing up, the height-diameter relationship is significantly influenced by the stand density and site index. However, there may be interaction effects of the stand density and site index, which has been disregarded in the previous height-diameter modeling studies. When the density remains the same, the better the site quality, the better the tree growth [14]. When the site index remains the same, the kurtosis and peak value of the small-diameter timber yield could decrease with an increase of the stand density [14]. Therefore, exploring the methods of quantifying the impact of the interaction effects of stand density and site index on the tree height-diameter relationship is very important.

Stand density and site index are random factors that produce random effects on tree height. The nonlinear mixed-effects model (NLMEM) can adequately describe the influence of the random effects on the height-diameter relationship. In addition, forestry data are generally characterized by hierarchical structures (multiple measurements within the same subject, such as sample plots within the same site); the NLMEM can effectively deal with the observation dependence and correlations generated by hierarchical data structure [11,15,16,17]. Therefore, the NLME modeling is very popular in forestry; for example, many modelers have applied NLME modeling to develop height-diameter models [2,15,18,19], height to crown base models [10], and other applications [20]. All these studies show that the accuracy of the NLMEM is remarkably higher than any other modeling approach, such as ordinary least squares regression. However, as mentioned earlier, the existing NLME height-diameter models do not cope with the interaction effects of stand density and site index. In this study, the NLME height-diameter model was developed through the inclusion of the interaction effects of stand density and site index. The parameter values of this model were obtained by using the linearization approximation-sequential quadratic algorithm proposed by Fu [21,22] and implemented on Forstat2.2 [23].

Larix olgensis Henry, also known as Changbai larch in China, has economic and ecological values, such as fast growth, most desirable wood properties, and strong resistance to diseases and insect pests. Larix olgensis is widely distributed in Northern China, North Korea, Russia, and other places [16]. Larix sp. is of a great interest to foresters because they vary greatly, grow rapidly, produce high-quality wood for many purposes, and provide exciting opportunities for short rotation forestation [24]. Obtaining a series of tree or stand attribute data in a timely and effective manner is essential for the reasonable management 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:

• Quantify the interaction effects of the stand density and site index on the tree height of Larix olgensis Henry.
• Determine whether there would be significant interaction effects of the stand density and site index on the height-diameter relationship.
• 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.

2. Materials and Methods

2.1. 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² to 2500 m² 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). Within each of the sample plots, all standing and living trees (total height, H > 1.3 m and diameter at breast height, D > 5 cm) were measured. The Research Institute of Forest Resources Information Techniques, Chinese Academy of Forestry established the PSPs and maintained the database.

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. Range of stand density class and site index class for Larix olgensis. SD is stand density; SI is site index.

Variable	Class	Class Range
M (Stand density class)	1	SD < 700
trees ha−1	2	700 ≤ SD < 850
	3	850 ≤ SD < 1000
	4	1000 ≤ SD < 1150
	5	1150 ≤ SD < 1300
	6	≥1300
S (Site index class)	1	SI < 15.5
m	2	15.5 ≤ SI < 16.0
	3	16.0 ≤ SI < 16.5
	4	16.5 ≤ SI < 17.0
	5	SI ≥ 17.0

).

A	Model Fitting Data				Model Testing Data
	Sample Plot	M	S	Number	Sample Plot	M	S	Number
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
1	304	3	3	30	304	3	3	15
1	305,307	4	1	30	305,307	4	1	15
1	306	6	3	30	306	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

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).

A	Model Fitting Data				Model Testing Data
	Sample Plot	M	S	Number	Sample Plot	M	S	Number
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
1	304	3	3	30	304	3	3	15
1	305,307	4	1	30	305,307	4	1	15
1	306	6	3	30	306	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

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.

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⁻¹. 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 3. Summary statistics of stand variables from sample plots used for model fitting and model validation, respectively. D= diameter at breast height (cm); H = tree height (m); SD = stand density (tree ha⁻¹); SI = site index (m); STD = Standard deviation of the variable.

Data	Variable	Min	Max	Mean	STD
Model fitting data
	Area (m2)	775	2500	1500	600
	D (cm)	5.60	35.50	17.69	5.30
	H (m)	5.50	22.00	16.13	2.52
	SD (trees ha−1)	640	1819	1048	266
	SI (m)	15.03	17.30	16.04	0.67
Model testing data
	D (cm)	7.10	33.30	17.61	5.18
	H (m)	7.00	22.00	16.25	2.79
	SD (trees ha−1)	640	1490	1031	277
	SI (m)	15.02	17.30	16.09	0.69

, and height-diameter relationship plotted against D is shown in Figure 2.

2.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.

$$AIC=2k-\ln (L)$$

(1)

$$BIC=k\ln (N)-2\ln (L)$$

(2)

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.

Basic Model	Function Expression	Function Form	Source
M1	$H=1.3+{\varphi }_1{D}^{{\varphi }_2}$	Power function	[31]
M2	$H=1.3+\exp ({\varphi }_1+\frac{{\varphi }_2}{D})$	Growth	[32]
M3	$H=1.3+{\varphi }_1(1-\exp (-{\varphi }_2{D}^{{\varphi }_3}))$	Weibull	[33]
M4	$H=1.3+{\varphi }_1{(1-\exp (-{\varphi }_{2}D))}^{{\varphi }_3}$	Chapman-Richards	[34]
M5	$H=1.3+{{\displaystyle \varphi }}_1+\frac{{{\displaystyle \varphi }}_2}{D+{{\displaystyle \varphi }}_3}$	Richards	[35]
M6	$H=1.3+{{\displaystyle \varphi }}_1\exp (-{{\displaystyle \varphi }}_2\exp (-{{\displaystyle \varphi }}_{3}D))$	Gompertz	[36]
M7	$H=1.3+\frac{{{\displaystyle \varphi }}_1}{1+\frac{1}{{{\displaystyle \varphi }}_2{{\displaystyle D}}^{{{\displaystyle \varphi }}_3}}}$	Hossfeld IV	[37]
M8	$H=1.3+{{\displaystyle \varphi }}_1\exp (-{{\displaystyle \varphi }}_2{{\displaystyle D}}^{-{{\displaystyle \varphi }}_3})$	Korf	[38]
M9	$H=1.3+\frac{{{\displaystyle \varphi }}_1}{1+{{\displaystyle \varphi }}_2\exp (-{{\displaystyle \varphi }}_{3}D)}$	Logistic	[39]

Note: H = tree height (m); D = diameter at breast height (cm). ${\varphi }_1$, ${\varphi }_2$, and ${\varphi }_3$ are the formal parameters to be estimated.

Table 4. Candidate base models we considered.

Basic Model	Function Expression	Function Form	Source
M1	$H=1.3+{\varphi }_1{D}^{{\varphi }_2}$	Power function	[31]
M2	$H=1.3+\exp ({\varphi }_1+\frac{{\varphi }_2}{D})$	Growth	[32]
M3	$H=1.3+{\varphi }_1(1-\exp (-{\varphi }_2{D}^{{\varphi }_3}))$	Weibull	[33]
M4	$H=1.3+{\varphi }_1{(1-\exp (-{\varphi }_{2}D))}^{{\varphi }_3}$	Chapman-Richards	[34]
M5	$H=1.3+{{\displaystyle \varphi }}_1+\frac{{{\displaystyle \varphi }}_2}{D+{{\displaystyle \varphi }}_3}$	Richards	[35]
M6	$H=1.3+{{\displaystyle \varphi }}_1\exp (-{{\displaystyle \varphi }}_2\exp (-{{\displaystyle \varphi }}_{3}D))$	Gompertz	[36]
M7	$H=1.3+\frac{{{\displaystyle \varphi }}_1}{1+\frac{1}{{{\displaystyle \varphi }}_2{{\displaystyle D}}^{{{\displaystyle \varphi }}_3}}}$	Hossfeld IV	[37]
M8	$H=1.3+{{\displaystyle \varphi }}_1\exp (-{{\displaystyle \varphi }}_2{{\displaystyle D}}^{-{{\displaystyle \varphi }}_3})$	Korf	[38]
M9	$H=1.3+\frac{{{\displaystyle \varphi }}_1}{1+{{\displaystyle \varphi }}_2\exp (-{{\displaystyle \varphi }}_{3}D)}$	Logistic	[39]

Note: H = tree height (m); D = diameter at breast height (cm). ${\varphi }_1$, ${\varphi }_2$, and ${\varphi }_3$ are the formal parameters to be estimated.

2.3. 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.

2.4. 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:

$$H_{ij}=f({\varphi }_i,x_{ij})$$

(3)

with

$${\varphi }_i=A_i\beta +Z_i{u}_i,$$

where ${\varphi }_i$ is formal parameter vector and includes the fixed effect parameter vector $\beta$ and random effect parameter vector $u_i$ of the ith sample plot; symbols $A_i$ and $Z_i$ are the design matrices for $\beta$ 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:

$${\widehat{u}}_i=\widehat{\psi }{{\widehat{Z}}_i}^T{({\widehat{Z}}_i\psi {{\widehat{Z}}_i}^T+{\widehat{R}}_i)}^{-1}(y_i-f(\widehat{\beta },u_i*,x_i)+{\widehat{Z}}_i{\widehat{u}}_i)$$

(4)

where $\widehat{\psi }$ is the estimated variance–covariance for the random effects, ${\widehat{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 $\psi$, and $R$($\psi ={{\displaystyle L}}^{\mathrm{{\rm 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 is ${\widehat{R}}_i={\sigma }^{2}I$(${\sigma }^2$ is the variance of the residual; I is the identity matrix.). The value of variance matrix $\widehat{\psi }$ or covariance matrix ${\widehat{R}}_i$ was calculated by restricted maximum likelihood with the sequential quadratic algorithm [21]. The $f(·)$ is an interactive NLME model, and ${\widehat{Z}}_i$ is an estimated design matrix:

$${\widehat{Z}}_i=\frac{\partial f(\beta ,u_i,x_i)}{\partial u_i}$$

(5)

where $x_i$ is a vector of the predictor on the sample plot i.

2.5. Model Evaluation

We used five statistical indicators to evaluate the performance of the interactive NLME height-diameter models such as $MPSE$,$RMSE$, and ${{\displaystyle R}}^2$ calculated by Equations (6)–(9) using both the model fitting and model testing data sets.

$$MPSE={\displaystyle \sum _{k=1}^N\frac{\left|\left(H_k-{\widehat{H}}_k\right)/{\widehat{H}}_k\right|}{N}}\times 100$$

(6)

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{k=1}^N{{\displaystyle (H_k-{\widehat{H}}_k)}}^2}}{N}}$$

(7)

$$\overline{H}=\frac{1}{N}{\displaystyle \sum _{k=1}^N{{\displaystyle H}}_k}$$

(8)

$${{\displaystyle R}}^2=\frac{{\displaystyle \sum _{k=1}^N({\widehat{H}}_k-\overline{H})}}{{\displaystyle \sum _{k=1}^N({{\displaystyle H}}_k-\overline{H})}}$$

(9)

where $H_k$ and ${\widehat{H}}_k$ are the observed and estimated values of total height of the $k\mathrm{th}$ tree, respectively; $\overline{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 ${{\displaystyle 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 height-diameter 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.

3. Results

3.1. 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. AIC (Akaike information criterion) score and BIC (Bayesian Information Criterion) score of the base models (Table 4).

Model	AIC	BIC
M1	2818.65	2832.57
M2	2721.78	2735.70
M3	2721.78	2737.93
M4	2718.63	2737.19
M5	2718.85	2737.41
M6	2720.33	2738.89
M7	2718.05	2736.61
M8	2718.29	2736.85
M9	2723.98	2742.54

), 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.

$$H_{(ij)k}=1.3+\exp \left[{\varphi }_1+\frac{{\varphi }_2}{D_{(ij)k}}\right]+{\epsilon }_{(ij)k}$$

(10)

where H is the total height of Larix olgensis, D is the diameter at breast height of Larix olgensis, ${\varphi }_1$ and ${\varphi }_2$ are the formal parameters of this model. ${{\displaystyle \epsilon }}_{(ij)k}$ is the error term of kth tree on the sample plot (ij).

3.2. The NLME Models

A total of 36 combinations of the random effects (

Table 6. Thirty-six alternatives of the random effect constructions in the nonlinear mixed-effects models.

Model	${\mathit{\varphi }}_1$	${\mathit{\varphi }}_2$	Model	${\mathit{\varphi }}_1$	${\mathit{\varphi }}_2$
1	M	M	19	M + M ∗ S	M
2	M	S	20	M + M ∗ S	S
3	S	M	21	M + M ∗ S	M ∗ S
4	S	S	22	M + M ∗ S	M + M ∗ S
5	M	M ∗ S	23	M + M ∗ S	S + M ∗ S
6	M	M + M ∗ S	24	M + M ∗ S	M + S + M ∗ S
7	M	S + M ∗ S	25	S + M ∗ S	M
8	M	M + S + M ∗ S	26	S + M ∗ S	S
9	S	M ∗ S	27	S + M ∗ S	M ∗ S
10	S	M + M ∗ S	28	S + M ∗ S	M + M ∗ S
11	S	S + M ∗ S	29	S + M ∗ S	S + M ∗ S
12	S	M + S + M ∗ S	30	S + M ∗ S	M + S + M ∗ S
13	M ∗ S	M	31	M + S + M ∗ S	M
14	M ∗ S	S	32	M + S + M ∗ S	S
15	M ∗ S	M ∗ S	33	M + S + M ∗ S	M ∗ S
16	M ∗ S	M + M ∗ S	34	M + S + M ∗ S	M + M ∗ S
17	M ∗ S	S + M ∗ S	35	M + S + M ∗ S	S + M ∗ S
18	M ∗ S	M + S + M ∗ S	36	M + S + M ∗ S	M + S + M ∗ S

Note: columns ${\varphi }_1$ and ${\varphi }_2$ show the random effects formulation variables that are acting on parameter ${\varphi }_1$ and parameter ${\varphi }_2$, respectively. M ∗ S denotes the crossed random effects of M (stand density) and S (site index). Symbol ∗ is not a simple multiplication; it means the crossed effects of variables.

) were derived from the formal parameters (${\varphi }_1$ and ${\varphi }_2$ of the base model Equation (10) affected by the random variables M (stand density class), S (site index class), and the crossed random effect of the stand density and site index (M ∗ S). All the possible combinations of the random effects were applied using model fitting data. The best performing combination was then selected based on the AIC and BIC scores.

Both the AIC and BIC of the interactive NLME model 15 are the smallest (

Table 7. AIC (Akaike information criterion) score and BIC (Bayesian Information Criterion) score of 36 models derived from the base model.

Model	AIC	BIC	Model	AIC	BIC
1	1798.34	1823.72	19	1787.53	1817.14
2	1796.58	1817.73	20	1792.20	1817.58
3	1793.00	1814.15	21	1773.64	1803.25
4	1804.01	1829.34	22	1777.64	1815.71
5	1783.41	1804.56	23	1789.41	1823.26
6	1785.13	1814.74	24	1779.64	1821.94
7	1785.41	1810.80	25	1789.22	1814.61
8	1787.12	1820.97	26	1794.21	1823.82
9	1783.41	1804.57	27	1773.64	1803.25
10	1785.37	1810.76	28	1775.64	1809.48
11	1787.41	1817.03	29	1777.64	1815.71
12	1789.37	1823.22	30	1779.64	1821.94
13	1787.23	1808.39	31	1789.53	1823.37
14	1790.22	1811.37	32	1796.19	1830.03
15	1771.64	1797.02	33	1775.64	1809.48
16	1773.64	1803.25	34	1779.64	1821.94
17	1773.64	1803.25	35	1779.64	1821.94
18	1775.64	1809.48	36	1783.64	1834.40

). Therefore, the random effects constructions M ∗ S and M ∗ S were selected as the random variables on ${\varphi }_1$ and ${\varphi }_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.

3.3. 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).

$$\{\begin{array}{l}{H}_{(ij)k}=1.3+\exp \left[{\beta }_1+u_{1(ij)}^{(M*S)}+\frac{{\beta }_2+u_{2(ij)}^{(M*S)}}{D_{(ij)k}}\right]+{\epsilon }_{(ij)k}\\ u_{(ij)}~N(0,\psi ),{\epsilon }_{(ij)}~N(0,R_{(ij)})\\ \\ i=1,\dots ,M_1, j=1,\dots ,M_2, k=1,\dots ,n_{(ij)}\end{array}$$

(11)

where M₁ (in this study, M₁ = 6) is the total classes of stand density, M₂ (in this study, M₂ = 5) is the total classes of site index, ${{\displaystyle 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 ${{\displaystyle \epsilon }}_{(ij)k}$ is the error term of the kth tree in sample plot (ij), which we assumed as $R_{(ij)}={\sigma }^{2}I$ (${{\displaystyle \sigma }}^2$ (${{\displaystyle \sigma }}^2>0$) is the variance of the residual). H is the tree height measurement value.

3.4. 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.

$$\widehat{H}=1.3+\exp \left[3.1194+\frac{-6.9169}{D}\right]+\epsilon$$

(12)

where $D$ is the measured value of DBH, $\widehat{H}$ is the estimated total tree height with the NLS model, $\epsilon \sim N\left(0,1.952I\right)$.

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:

$${\widehat{H}}_{(ij)k}=1.3+\exp \left[3.1138+u_{1(ij)}^{(M*S)}+\frac{-6.8180+u_{2(ij)}^{(M*S)}}{D_{(ij)k}}\right]+{\epsilon }_{(ij)k}$$

(13)

with

$${{\displaystyle u}}_{ij}=\left[\begin{array}{l}{{\displaystyle u}}_{1ij}\\ {{\displaystyle u}}_{2ij}\end{array}\right]\sim \left\{\left[\begin{array}{l}0\\ 0\end{array}\right],\psi =\left(\begin{array}{ll}0.0053& -0.0981\\ -0.0981& 1.9956\end{array}\right)\right\}, {\epsilon }_{ij}\sim N\left(0,1.6996I\right),$$

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 ${\widehat{H}}_{(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: ${\beta }_1$ = $3.1138$, ${\beta }_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 approximation-sequential quadratic algorithm implemented in the “nonlinear mixed-effects” module of the Forstat software of the 2.2 version is given by:

$${\widehat{H}}_{(ij)k}=1.3+\exp \left[ 3.1165+u_{(i)}^{(S)}+\frac{-6.8545+u_{(j)}^{(M)}}{D_{(ij)k}}\right]+{\epsilon }_{(ij)k}$$

(14)

with

$$u_{(i)}~N(0,0.0003),u_{(j)}~N(0,0.1216),{\epsilon }_{ij}~N(0,1.8505I),$$

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 ${\widehat{H}}_{(ij)k}$ is the estimated H of tree k on the (ij) sample plot with a single-level NLME model.

3.5. 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. Estimated values of the random effects.

A	M ∗ S	u1	u2	A	M ∗ S	u1	u2
1	3 ∗ 2	−0.0592	1.4189	2	3∗ 4	0.0544	− 0.9706
1	6 ∗ 1	−0.0666	0.6652	2	1∗ 4	− 0.0261	0.9751
1	4 ∗ 2	0.0679	−0.8770	2	2 ∗ 1	0.1206	−2.3081
1	3∗ 3	0.0161	−0.2747	2	6 ∗ 1	−0.0666	0.6652
1	4 ∗ 1	−0.0274	0.3191	2	4 ∗ 1	−0.0274	0.3191
1	6 ∗ 3	−0.0349	1.1108	2	5 ∗ 2	−0.0111	−0.2363
1	1 ∗ 5	0.0118	−0.0328	2	2 ∗ 4	−0.1043	2.1220
2	5 ∗ 5	0.0044	−0.3663	2	6 ∗ 3	−0.0349	1.1108
2	2 ∗ 3	0.0543	−1.5452

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 ${\varphi }_1$ and ${\varphi }_2$ of the interactive NLME height-diameter model (Equation (13)), respectively.

. 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.

3.6. 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:

$$LR=2(\ln L1-\ln L2)~{\chi }^2(df1-df2)$$

(15)

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, $df1$ is the degree of freedom of the complex model, and $df2$ is the freedom of the simple model.

The likelihood ratio test results are shown in

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))	27.81	p < 0.001
Interactive NLMEM (Equation (13)) vs. Single-level NLMEM (Equation (14))	23.36	p < 0.001

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.

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 ${{\displaystyle R}}^2$ (Equations (6)–(9)), of three height-diameter models (Equations (12)–(14)) (

Table 10. Evaluation indicators for the NLS model (Equation (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)).

Data Set	Model	MPSE	RMSE	${{\displaystyle \mathit{R}}}^2$
Model fitting data	NLS model	6.8618	1.1820	0.6880
	Single-level NLMEM	6.6037	1.3450	0.6991
	Interactive NLMEM	6.3076	1.1306	0.7189
Model testing data	NLS model	7.4793	1.2306	0.5717
	Single-level NLMEM	7.2896	1.4591	0.5799
	Interactive NLMEM	6.8234	1.1803	0.6234

). 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.

4. 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² 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 ${{\displaystyle R}}^2$ = 0.6234) was better than the NLS model (Equation (12) with ${{\displaystyle 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 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.

5. 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 height-diameter 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.