Simultaneous Compatible System of Models of Height, Crown Length, and Height to Crown Base for Natural Secondary Forests of Northeast China

Abstract

Individual trees are characterized by various sizes and forms, such as diameter at breast height, total height (H), height to crown base (HCB), crown length (CL), crown width, and crown and stem forms. Tree characteristics are strongly related to each other, and studying their relationships is very important. The knowledge of the compatibility and additivity properties of the major tree characteristics, such as H, CL, and HCB, is essential for informed decision-making in forestry. H can be used to represent site quality and CL represents biomass and photosynthesis of crown, which is the performance of individual tree vigor and light interception, and the longer the crown length (or shorter HCB) is, the more vigorous the tree would be. However, none of the studies have uncovered their inherent relationships quantitatively. This study attempts to explore such relationships through the application of appropriate modeling approaches. We applied seemingly unrelated regression, such as nonlinear seemingly unrelated regression (NSUR), which is commonly used for exploring the compatibility and additivity properties of the variables, for the proposes. The NSUR involves the variance and covariance matrices of the sub-models that are used for the interpretation of the correlations among the variables of interest. The data set acquired from Mongolian oak forest and spruce-fir forest in the Jingouling forest farm of the Wangqing Forest Bureau in the Northeast of China were used to construct two types of model systems: a compatible model system (the model system of H, CL, and HCB can be estimated simultaneously) and an additive model system (the sum of HCB and CL is H, the form of the H sub-model equals the sum of the HCB and CL sub-models) from the individual models of H, CL, and HCB. Among the various tree-level and stand-level variables evaluated, D (diameter at breast), Dg (quadratic mean diameter), DT (dominant diameter), CW (crown width), SDI (stand density index), and BAS (basal area of stand) contributed significantly highly to the variations of the response of the variables of interest in the model systems. Modeling results showed the existence of the compatibility and additivity of H, CL, and HCB simultaneously. The additive model system exhibited better fitting performance on H and HCB but poorer fitting on CL compared with the simultaneous model system, indicating that the performance of the additive model system could be higher than that of the simultaneous model system. Model tests against the validation data set also confirmed such results. This study contributes a novel approach to solving the compatibility and additivity of the problems of H, CL, and HCB models through the application of the robust estimating method, NSUR. The results and algorithm presented will be useful for constructing similar compatible and additive model systems of multiple tree-level models for other tree species.

1. Introduction

Informed decision-making is of great importance in forestry [1,2]. Forest managers usually make effective plans using quantitative information of past and present stand conditions, such as inventory data, growth and yield models, for prediction of the future stand conditions [3]. Individual tree size, such as diameter at breast height (D), total tree height (H), height to crown base (HCB), crown length (CL), ratio of D to quadratic mean D (Dq), basal area (BA), as well as crown width (CW), play an indispensable role for constructing tree growth and yield models. D is usually used as a predictor of the D-H model because of its measurement convenience and comparable lower measurement error (degree of accuracy degree 0.1 cm) [4,5,6,7]. D is also used as one of the main predictors in many other tree-based models, including the HCB model [8], the CW model [9], the CL model, the crown ratio model, the mortality model [10,11], the recruitment model, the site index and the dominant height model [5,12], the biomass model [13,14], and so on.

As mentioned above, D is one of the crucial tree variables for growth and yield models. Besides this, HCB, H, and CL are equally important for these models. Tree height (H), especially dominant height, is the reflection of site quality [5,12]; the position of HCB is the indicator of the crown ratio, while the crown ratio is the indicator of tree vigor and photosynthesis [15], competition, individual tree foliage and branch biomass [16,17,18], wood quality [19,20], and wind firmness [21]. H and HCB are inherently correlated in a tree, and thus height growth would always be a factor to influence HCB. HCB can be estimated from H, and H can also be estimated from HCB [22]; the latter may be more acceptable due to less cost and time requirement. Moreover, HCB is a member of the crown, attributed with CW, crown profile, crown biomass, crown volume, crown surface area, CL, light intercepting surface, site of photosynthesis, and so on [2,9,23,24,25,26,27,28,29,30]. Furthermore, H is simply a sum of CL and HCB, and CL is as essential as HCB, even more essential in estimating crown measures and tree vigor [31,32], crown volume, biomass, and taper [33,34]. The correlation between H, CL, and HCB is very high and inherent. All in all, H, HCB, and CL are key measures to describe growth, yield, mortality, and recruitment dynamics of forest stands. However, so far, there is no research on the compatibility and additivity of H, CL, and HCB.

Most of the existed studies with regard to compatibility are concentrated on the compatibility and additivity of components and the total parts. For example, the compatibility and additivity of biomass for slash pine (Pinus elliottii Engelm. var. elliottii) trees were carried out [35] in the southeastern United States with three modeling approaches, and they found no single system to predict biomass that is the most accurate for all the components and total tree biomass. These problems in terms of compatibility and additivity almost always concentrate on regression models rather than classification models. When it comes to classification problems, the machine learning methods would be the best choice rather than other regression methods [36,37,38,39,40,41] in forest fire supervision, forest biomass mapping, parameters extraction of remote sensing figures, and other aspects. However, when it comes to regression problems, traditional methods would usually be chosen first, such as NSUR (nonlinear seemingly unrelated regression method), OLS (ordinary least square method), and QR (quantile regression method), et al.

The efficiency is ordered by SUR2 (jointly fitting component biomass equations using weighted NSUR) > SUR1 (weighted nonlinear seemingly unrelated regression (NSUR)) > DRM (an alternative disaggregative approach). In addition, all these three model systems provided more accurate biomass predictions than the previously published biomass equations. Similarly, Dong et al. [42] compared the nonadditive against the additive biomass models for larch (Larix olgensis Henry) trees in Northeast China. Jiang et al. [43] proposed a simultaneous compatible system of taper and volume equations. Jose et al. [44] constructed compatible height and site index models for five pine species to ensure the simultaneous estimation of height and site index. The nonlinear seemingly unrelated regression (NSUR) is a commonly used method to estimate parameters of these compatible-simultaneous systems of the models of biomass, taper, volume, and site index. With the model system estimated, the inherent correlations of the components are effectively accounted for and guaranteed as compatible, which is not possible with the ordinary least square regression. The accuracy produced using the NSUR is relatively higher, as errors of each individual model involved in a system are simultaneously reduced. There exist other methods of estimating simultaneous systems of equations besides NSUR; for example, Yang et al. [45] utilized full information maximum likelihood (FIML), two-stage least square (2SLS), and three-stage least squares (3SLS) methods to construct simultaneous individual tree diameter and height-to-crown base, which sound great for data coming from remote sensors of LiDAR; these methods provide a convenient way to estimate these indexes, which can be applied to remote sensing and can shrink the cost to acquire them. However, the application conditions are constrained by some prerequisites. Nevertheless, no studies have been done on the compatibility and additivity of H, CL, and HCB with these methods.

Thus, based on the knowledge of the previous modeling studies, this study attempts to develop a compatible system and an additive system of H, CL, and HCB models to simultaneously estimate H, CL, and HCB with the data acquired from permanent sample plots from the northeast Jilin province of China. Then, in addition to NSUR, FIML [46], two-stage least squares (2SLS), and three-stage least squares (3SLS) will also be used to try to estimate the parameters of models with different forms of models. The study develops the compatible and additive simultaneous systems of H, CL, and HCB models and compares the discrepancy of NSUR and other methods using the statistical indices. The presented system of the models will be useful for the effective management of natural secondary forests.

2. Materials and Methods

2.1. Study Area

The study was carried out in the Jingouling experimental forest farm (129°97′–130°22′ E, 43°32′–43°49′ N) and the Tazigou experimental forest farm (130°5′–130°20′ E, 43°17′–43°25′ N) of the Wang Qing forest bureau in Jilin province, China. The area of the Wang Qing forest bureau belongs to the Xueling branch of the Laoyeling mountain range of the Changbai Mountain system, which has a temperate continental monsoon climate, with four distinct seasons, long winter and short summer, and a cold climate. The first frost period is generally in the middle of September, the last frost period is generally in the middle of May, the frost-free period is 110–130 days, and the accumulated temperature is 2110 °C. The altitude is about 360–1477 m. The soil is brown forest soil, with neutral pH and 45 cm thickness.

The stand type is spruce-fir coniferous broad-leaved mixed natural secondary forest and Mongolian oak broad-leaved mixed natural secondary forest. Spruce-fir coniferous broad-leaved mixed natural secondary forest is a forest vegetation type left over from the primitive zonal forest vegetation (spruce-fir coniferous and broad-leaved mixed forest and broad-leaved Korean pine forest) of Chang bai Mountain after many “plucking” cuttings. Mongolian oak broad-leaved mixed natural secondary forest is widely distributed in northeast China, mostly in the secondary succession stage of forest communities formed after the over-cutting of broad-leaved Korean pine forest, with low forest productivity and degraded ecological function.

2.2. Sampling and Measurements

The data used in this study were collected from 384 permanent sample plots (PSPs) with the same area of 625 m² established across the Tazigou forest farm and the Jingouling forest farm out of 28 forest classes between July and August in 2013. Two forest types were included, where 192 PSPs were located in Mongolian oak forest (Mongolian oak is the main tree species) and another 192 PSPs were located in spruce-fir forest (spruce-fir is the main tree species). Sample plots are representative of the various conditions, stand density, and stand stocking. All standing living trees with diameter at breast height (D) ≥ 1 cm and total tree height ≥ 1.3 m were measured for D, total height (H), height-to-crown base (HCB), and four crown radii. Two azimuths were used to determine the accurate four-crown radii positions; one of the azimuths is perpendicular to the other one, and their direction was from south to north and west to east. In each quadrant, the crown radii were measured from the center of the target tree bole to the great extent of the crown from the target tree bole. Crown width was computed as the arithmetic mean of two-crown widths, and stand density was calculated by the sum of trees in every plot; then, the stand density index would be calculated: SDI = N(D_g/D₀)^1.605, where D_g and D₀ are the stand quadratic diameter and standard average diameter; in China, D₀ is 20 cm, N is the number of trees in the forest stand, and 1.605 represents the self-thinning slope [47]. Dominant diameter at breast height (DT) was calculated according to the average value of the top five sample trees in terms of D in every plot. Other stand variable calculation methods were used commonly. Out of the total measured trees (27,574), trees with DBH < 5 cm and height < 1.3 m were excluded, which retained 15,301 trees to be used in the modeling. This dataset was randomly partitioned into model fitting and model validation datasets, consisting of 70% and 30% of the total dataset, respectively. A data summary is presented in

Table 1. Model-fitting and model-validation data (Std: standard deviation; D: diameter at breast height (cm); H: total tree height (m); HCB: height to live crown base (m); CW: crown width (m); CL: crown length (m); BAS: stand basal area (m² ha⁻¹); DT: stand dominant D (cm); SDI: stand density index).

Stand Type	Variable	Model-Fitting Data				Model-Validation Data
		Max	Min	Mean	Std	Max	Min	Mean	Std
Quercus mongolica Fisch. ex Ledeb.	D (cm)	59.40	5.20	14.68	8.75	59.10	5.00	14.58	8.47
	Dg (cm)	22.40	13.58	18.34	1.54	26.67	12.82	16.71	2.25
	H (m)	36.60	1.60	11.19	4.73	29.40	1.50	11.28	4.66
	HCB (m)	19.80	0.30	4.62	2.86	17.40	0.20	4.76	2.92
	CW (m)	15.55	0.45	3.86	1.72	14.40	0.40	3.87	1.72
	CL (m)	31.40	0.20	6.58	3.69	22.40	0.20	6.52	3.59
	BAS (m2 ha−1)	24.10	5.75	13.24	3.54	30.74	8.04	20.24	4.74
	DT (cm)	45.04	20.96	30.50	4.71	44.14	25.22	32.64	4.75
	SDI	172.86	13.56	72.15	30.77	319.13	11.24	160.22	64.61
Spruce fir	D (cm)	55.50	5.00	16.37	8.41	63.87	14.28	43.31	9.72
	Dg (cm)	23.48	14.22	18.38	1.89	26.18	12.22	16.92	2.35
	H (m)	34.20	2.30	15.63	6.07	55.50	5.00	16.37	8.41
	HCB (m)	21.50	0.40	6.72	3.85	34.20	2.30	15.63	6.07
	CW (m)	12.85	0.77	4.27	1.67	21.50	0.40	6.72	3.85
	CL (m)	26.20	0.40	8.91	4.44	12.85	0.77	4.27	1.67
	BAS (m2 ha−1)	32.16	6.48	23.18	5.25	26.20	0.40	8.91	4.44
	DT (cm)	40.98	21.96	33.56	3.69	8.96	0.61	2.41	1.27
	SDI	67.39	15.04	47.72	10.58	26.18	12.22	16.92	2.35

. Relationships of H, CL, and HCB with D, CW, BAS, Dg, DT, and SDI are shown in Figure 1 and Figure 2.

2.3. Selecting the Base Model

In this study, the monomolecular model for H and HCB and the logistics model for CL were chosen as base models to develop the systems of H, CL, and HCB models. We also evaluated other model forms in our preliminary analysis; however, the monomolecular model for H and HCB and the logistics model for CL performed better in terms of the convergence and reduction of the sum of squared errors of both the model system and individual models of H, CL, and HCB. Our other evaluated models are the Richards model, the Korf model, and the Weibull model. Bronisz and Mehtätalo [48] chose a monomolecular model and a logistics model as the best models to fit H and CW for young silver birch stands on post-agricultural lands in central Poland and individual trees of Chinese fir (Cunninghamia lanceolata) in south-central China, as well as a height-to-crown base model (HCB) for Mongolian oak in northeast China, based on covariates including D, H, HCB. The following are the forms of base models that showed better fitting performance than many of the other models evaluated.

Base models:

$$H=a_0(1-\exp (-{\mathrm{a}}_{1}D))+{\epsilon }_H$$

(1)

$$HCB=c_0(1-\exp (-c_{1}D))+{\epsilon }_{HCB}$$

(2)

$$CL=b_0/(1+b_1\exp (-b_{2}D))+{\epsilon }_{CL}$$

(3)

where a₀, a₁, b₀, b₁, b₂, c₀, and c₁ are parameters to be estimated, and ${\epsilon }_H$, ${\epsilon }_{CL}$, and ${\epsilon }_{HCB}$ are residual errors of the models H, CL and HCB, respectively.

2.4. Selecting Additional Independent Variables

D was chosen as the main independent variable for all three individual models because its correlations with dependent variables (H, CL, HCB) were highly significant. As mentioned earlier, H, CL, and HCB are inherently correlated to each other, and H is the sum of the other two. We consider the tree variables and stand variables that could improve the fitting accuracy of the base model (Equations (1)–(3)). They were included in the base models as covariate predictors; the final model forms with the added covariate predictors were determined through some statistical measures to be defined later, and in order to avoid multicollinearity among predictor variables, only those predictors that have insignificant correlations or correlation less than 0.3 were chosen as the predictors of models. Meanwhile, the constraint of the simultaneous system was defined by the matrix of variance-covariance of H, CL, and HCB. The individual tree-level variables and stand-level variables were finally chosen after many experiments. The variables that negatively affect the fitting accuracy and convergence of the models were eliminated.

2.5. Introducing a Dummy Variable

Since our data represent two forest types (spruce-fir forest, Mongolian oak forest), the effect of species-specific difference was described by one dummy variable that was included in the simultaneous model system (i.e., if species is spruce-fir, then sp = 1, otherwise 0). The error and trail approach was applied to locate the dummy variable to the right positions of the parameters in each base model.

2.6. Structure of Simultaneous Compatible and Additive Model Systems

A model system (MS) is composed of H, CL, and HCB base models (Equation (4)), which is also known as simultaneous MS. The method had been identified by several studies.

$$\{\begin{array}{l}H=a_0(1-\exp (-a_{1}D))+{\epsilon }_H\\ CL=b_0/(1+b_1\exp (-b_{2}D))+{\epsilon }_{CL}\\ HCB=C_0(1-\exp (-c_{1}D))+{\epsilon }_{HCB}\\ \epsilon ={\epsilon }_H+{\epsilon }_{CL}+{\epsilon }_{HCB}\end{array}$$

(4)

where the H, CL, and HCB denote the observed H, CL, and HCB values of measured trees, respectively; ${\epsilon }_H$, ${\epsilon }_{CL}$, and ${\epsilon }_{HCB}$ are error terms of each model in the simultaneous equation system, D is the diameter at breast height. a₀, a₁, b₀, b₁, b₂, c₀, and c₁ are the parameters to be estimated. $\epsilon$ is the three-dimensional error vector that is assumed to be normally distributed, with zero means and the variance–covariance matrix Σ. The matrix Σ contains variances and covariances to interpret the inherent correlations of H, CL, and HCB.

We assumed that the error terms (${\epsilon }_H$, ${\epsilon }_{CL}$, and ${\epsilon }_{HCB}$ commonly compose the error vector $\epsilon$) would be unrelated among the observations; however, they would be contemporaneously correlated between the sub-models. Observations of H, CL, and HCB would have a variance and covariance matrix (Equation (5)):

$$\Sigma =\left(\begin{array}{lll}{\sigma }_{H\times H}& {\sigma }_{H\times CL}& {\sigma }_{H\times HCB}\\ {\sigma }_{CL\times H}& {\sigma }_{CL\times CL}& {\sigma }_{CL\times HCB}\\ {\sigma }_{HCB\times H}& {\sigma }_{HCB\times CL}& {\sigma }_{HCB\times HCB}\end{array}\right)$$

(5)

As is known, H, CL, and HCB exist with the inherent correlations, and the sum of CL and HCB results in H. The form of an additive model system is shown in Equation (6):

$$\{\begin{array}{l}CL=b_0/(1+b_1\exp (-b_{2}D))+{\epsilon }_{CL}\\ \\ HCB=c_0(1-\exp (-c_{1}D))+{\epsilon }_{HCB}\\ \\ H=CL+HCB+{\epsilon }_H\\ \\ \epsilon ={\epsilon }_H+{\epsilon }_{CL}+{\epsilon }_{HCB}\end{array}$$

(6)

Our assumed error terms and structure of the variance-covariance matrix would be the same as defined in Equations (4) and (5).

2.7. Methods of Estimating Models

2.7.1. Ordinary Least Squares with Separating Regression

Base models (Equations (1)–(3)) and their modified variants were separately fitted by nonlinear ordinary least squares with separating regression (OLSSR).

2.7.2. Nonlinear Seemingly Unrelated Regression (NSUR)

Three base models (Equations (1)–(3)) were linked to form a single equation system through the inclusion of the noise terms. The hypothesis was that the noises or disturbances between the three models would be uncorrelated among all the observations. The covariance matrix of the disturbance terms, however, is generally unknown. A feasible generalized least squares (FGLS) algorithm is applied to estimate the matrix, and it is asymptotically efficient. The algorithm used for this propose is generally known as a “seemingly unrelated regression” (SUR) estimator. Parameters of the simultaneous equation system could be estimated by NSUR, with FGLS associated. The calculation steps we employed are below:

Step 1: The OLS method was used to fit the three base models to obtain residuals $\epsilon ={({\epsilon }_H,{\epsilon }_{CL},{\epsilon }_{HCB})}^T$. Then, model residuals were used to estimate the variance and covariance of the matrix Σ with a 3 × 3 dimension (Equation (5)). This matrix was used to interpret the inherent correlations of H, CL, and HCB.

Step 2: The estimate of covariance matrix $\widehat{\Sigma }$ was used to estimate the covariance matrix R in the form of covariance matrix $\widehat{R}=\widehat{\Sigma }\otimes I_N$. The FGLS is used to estimate the parameter vector (Equation (7) [49].

$$B={({\underset{K\times 3N}{X}}^T(\underset{3\times 3}{{\widehat{\Sigma }}^{-1}}\otimes I_N)\underset{3N\times K}{X})}^{-1}{\underset{K\times 3N}{X}}^T(\underset{3\times 3}{{\widehat{\Sigma }}^{-1}}\otimes I_N)\underset{3N\otimes 1}{y}$$

(7)

2.7.3. Attempting Other Estimation Methods

We also tried fitting simultaneous and additive MS using FIML (full information maximum likelihood), 2SLS (two-stage least square), and 3SLS (three-stage least square) with various start values from OLSSR but failed as none of the MS converged.

In the FIML, 2SLS, and 3SLS model systems, the endogenous variables (error-in variable), exogenous variable (error-out variable), and instrument variables (usually exogenous variables) are required to be set in advance. The instrument variables, endogenous variables, and exogenous variables were all needed in the 2SLS and 3SLS methods; the endogenous variables and exogenous variables were needed in the FIML method, but instrument variables were not required. Out of the prerequisites of these methods, the difficulty of model system convergence would be enlarged. When these methods were utilized to estimate the model system, their complexity will hinder the iteration and convergence.

2.8. Model Evaluation

Many statistical indices can be used to evaluate individual models and model systems constructed from them [6]. We applied the most common statistical measures for this purpose, such as coefficient of determination (R²), bias ($\overline{e}$), bias variance ($\delta$), root mean square error (RMSE), and total relative error (TRE) [50].

$$R^2=1-{\displaystyle \sum _{i=1}^N{(y_i-{\widehat{y}}_i)}^2}/{\displaystyle \sum _{i=1}^N{(y_i-{\overline{y}}_i)}^2}$$

(8)

$$\overline{e}={\displaystyle \sum _{i=1}^N{e}_i}/N={\displaystyle \sum _{i=1}^N(y_i-{\widehat{y}}_i)/N}$$

(9)

$$\delta ={\displaystyle \sum _{i=1}^N{(e_i-\overline{e})}^2/N-1}$$

(10)

$$RMSE=\sqrt{{\overline{e}}^2+\delta }$$

(11)

$$TRE={\displaystyle \sum (y_i-{\widehat{y}}_i)/{\displaystyle \sum {\widehat{y}}_i}}\times 100$$

(12)

where $y_i$ and ${\widehat{y}}_i$ are the observed and predicted values of H, CL, or HCB for the ith observation value, respectively, for the ith observation; ${\overline{y}}_i$ is mean value of samples; N is the number of samples.

Both the model systems (compatible and additive model systems) were tested against the 30% data set allocated for validation. The same statistical indices (Equations (8)–(12)) that were used to evaluate the fitting performance of the individual models and model systems were also used to evaluate the prediction performance of the simultaneous and additive model systems. Unless otherwise stated, we used a 5% confidence level to evaluate the individual models and model systems.

The basic data processing and calculation of individual tree variables and stand variables were carried on with R-3.6.1 software. SAS 9.4 software was used for estimating base models separately and two model systems with the NLIN and MODEL procedure in the SAS/ETS modular system [51].

The flowchart of the paper (Figure 3) is shown below:

3. Results

3.1. Base Models

OLS was applied to estimate the parameters of all the base models (1), (2), and (3) separately for H, CL, and HCB. The parameters estimates of the H, CL, and HCB models were all significant (p < 0.05). The model for H appeared to be better than those for CL and HCB (

Table 2. Evaluation indices of the base models (Equations (1)–(3)), as well as the expanded models (Equations (13)–(18)) (R²: coefficient of determination; $\overline{e}$: bias; $\delta$: bias variance; RMSE: root mean square error; TRE: total relative error).

Model	R2	$$\overline{\mathit{e}}$$	$$\mathit{\delta }$$	RMSE	TRE
Equation (1)	0.6123	−0.0943	14.3249	3.8119	6.9947
Equation (2)	0.5131	−0.0027	8.7434	2.9569	12.4000
Equation (3)	0.2599	−0.0572	9.4194	3.0696	26.4562
Equation (13)	0.6196	0.0747	13.4882	3.6733	6.5960
Equation (14)	0.5645	−0.0007	7.8187	2.7962	10.9460
Equation (15)	0.3062	−0.0439	8.8323	2.9722	24.4007
Equation (16)	0.6641	0.0999	11.9039	3.4516	5.7788
Equation (17)	0.6243	−0.0068	6.7464	2.5973	9.3043
Equation (18)	0.3443	−0.0406	8.3472	2.8894	22.7553

).

3.2. Inclusion of Covariate Predictor Variables

All the chosen predictor variables after the multicollinearity control were included in the base models, which resulted in more complex models (Equations (13)–(15)). The individual tree-level variables and stand-level variables were finally chosen via many experiments; the variables that negatively affect the fitting accuracy and convergence of the models were eliminated. The forms of H, CL, and HCB models are shown below:

$$H=\left(a_0{D}_g+a_{1}DT+a_{2}BAS+a_{3}CW\right)\left(1-\exp (-a_{4}D)\right)+{\epsilon }_H$$

(13)

$$CL=\frac{b_0}{(1+b_1\exp (-b_{2}D-b_{3}CW+b_4{D}_g+b_{5}SDI+b_{6}HCB))}+{\epsilon }_{CL}$$

(14)

$$HCB=\left(-c_{0}CW+c_1{D}_g+c_{2}BAS+c_{3}SDI+c_{4}DT\right)\left(1-\exp (-c_{5}D)\right)+{\epsilon }_{HCB}$$

(15)

We estimated Equations (13)–(15) with OLSSR, and the estimated parameters were all significant (p < 0.05). After the inclusion of additional covariate predictors to the base models, the expanded models were complex, but their parameter estimates (Equations (13)–(15)) were all significant (p < 0.05) (Table 2.) The result indicated that the inclusion of individual tree-level and stand-level variables to the CL and HCB models were a better way to enhance the accuracy and reduce the errors of these models compared with the original models, while in the H model, it was indicated more variables maybe limit the fitting ability of the model.

Then, one dummy variable, SP, was introduced to the expanded models (Equations (13)–(15)) to differentiate the two forest types: Mongolian oak and spruce-fir forest; spruce-fir sp = 1, otherwise 0. Equations (16)–(18) are the expanded forms with dummy variables included, and the parameter estimates and fit statistics are presented in Table 2:

Dummy H, CL, and HCB model:

$$H=\left((a_0+a_{01}SP)D_g+a_{1}DT+a_{2}BAS+a_{3}CW\right)\left(1-\exp \left(a_{4}D\right)\right)+{\epsilon }_H$$

(16)

$$CL=\frac{b_0+b_{01}SP}{(1+b_1\exp (-b_{2}D-b_{3}CW+b_4{D}_g+b_{5}SDI+b_{6}HCB))}+{\epsilon }_{CL}$$

(17)

$$HCB=\left(-c_{0}CW+\left(c_1+c_{11}SP\right)D_g+c_{2}BAS+c_{3}SDI+c_{4}DT\right)\left(1-\exp (-c_{5}D)\right)+{\epsilon }_{HCB}$$

(18)

The fitting improvement, obtained with the dummy variable included, was compared to the models without a dummy variable (Equations (13)–(15)). The simultaneously compatible model system and the additive model system were developed from Equations (16)–(18). The results of these are given in the following sub-sections.

3.3. Simultaneous H, CL, and HCB Model System Estimated with NSUR

When the equations of H (16), CL (17), and HCB (18) were formed as a system (19) to be estimated, HCB was used as a covariate in the CL model. The form of simultaneous MS of the models is shown in Equation (19). This MS was estimated using NSUR.

$$\{\begin{array}{l}H=\left((a_0+a_{01}SP)D_g+a_{1}DT+a_{2}BAS+a_{3}CW\right)\left(1-\exp \left(a_{4}D\right)\right)+{\epsilon }_H\\ \\ CL=\frac{b_0+b_{01}SP}{(1+b_1\exp (b_{2}D+b_{3}CW+b_4{D}_g-b_{5}SDI+b_{6}HCB))}+{\epsilon }_{CL}\\ \\ HCB=\left(-c_{0}CW+\left(c_1+c_{11}SP\right)D_g+c_{2}BAS+c_{3}SDI+c_{4}DT\right)\left(1-\exp (-c_{5}D)\right)+{\epsilon }_{HCB}\\ \\ \epsilon ={\epsilon }_H+{\epsilon }_{CL}+{\epsilon }_{HCB}\end{array}$$

(19)

All the parameters, except for a₂ and b₅, were significant (p < 0.05) (

Table 3. Parameter estimates of a system of the models of H, CL, HCB (Equation (19)) (H: total tree height; CL: crown length; HCB: height to live crown base; values in the parenthesis are standard errors).

Simultaneously Model System (Equation (19))
H Model		CL Model		HCB Model
Parameter	Estimates	Parameter	Estimates	Parameter	Estimates
a0	0.7292 (0.0258)	b0	13.2480 (0.1139)	c0	0.1437 (0.0293)
a01	0.2566 (0.0071)	b01	3.3560 (0.0981)	c1	0.2399 (0.0263)
a1	0.1705 (0.0179)	b1	22.8787 (1.2407)	c11	0.1282 (0.0062)
a2	−0.0306 (0.0184)	b2	−0.1222 (0.0017)	c2	−0.1965 (0.1579)
a3	0.0304 (0.0288)	b3	−0.0444 (0.0052)	c3	0.1123 (0.0274)
a4	0.0706 (0.0012)	b4	−0.0982 (0.0031)	c4	0.1032 (0.0190)
		b5	−0.0013 (0.0010)	c5	0.0780 (0.0028)
		b6	0.0917 (0.0019)

). The estimated parameters of dummy variables a₀₁, b₀₁, c₁₁ for the H model, the CL model, as well as the HCB model were significant (p < 0.0001), indicating the differences of H, CL, and HCB between Mongolian oak and spruce-fir. The magnitudes and signs of these parameters are biologically plausible. For the H model, the effect of stocking (BAS) negatively affects the asymptote, while D_g, DT, and CW have a positive effect. For the CL model, the effect of D, D_g, and nCW positively affect the asymptote, while SDI and HCB have a negative effect. For the HCB model, the effect of CW and stocking (BAS) negatively affects the asymptote, while D, D_g, DT, and SDI have a positive effect.

The fitting statistics of the simultaneous MS (Equation (19)) of H, CL, and HCB models are listed in

Table 4. Evaluation indices of (Equation (19)) (H: total tree height; CL: crown length; HCB: height to live crown base; R²: coefficient of determination; $\overline{e}$: bias; $\delta$: bias variance; RMSE: root mean square error; TRE: total relative error).

System	Sub-Model	R2	$$\overline{\mathit{e}}$$	$$\mathit{\delta }$$	RMSE	TRE
Compatible model (Equation (19))	H	0.6636	0.9268	11.9689	3.4554	7.0813
	CL	0.5923	0.1026	7.3113	2.7059	10.4018
	HCB	0.2860	0.0269	9.0904	3.0152	25.3271

. All statistics for the MS of the H, CL, and HCB models apparently had lower performance than those for the H, CL, and HCB models fitted separately, i.e., the values of R² and RMSE for the H model in the system of the model with NSUR were a little lower than those estimated with OLS. The OLS estimation was based on only the H model and never took correlations of H with the other two dependent variables, CL and HCB. There was no existing constraint between the three dependent variables. It indicated that when simultaneous MS did not exist, OLS was not necessarily worse than NSUR.

The estimates of the variance-covariance matrix R (Equation (4)) in the estimated system of the H, CL, HCB models by NSUR are shown below:

$$\sum =\left(\begin{array}{lll}11.92& 5.94& 4.71\\ 5.94& 6.75& -0.73\\ 4.71& -0.73& 8.35\end{array}\right)$$

for the system of the H, CW, and HCB models (19).

All the elements in the covariance matrix for (Equation (19)) were significant (p < 0.05), which suggested that the correlations among the H, CL, and HCB were definitely large and could not be ignored. The statistical indices of the simultaneously MS system fitted by NSUR were lower than those of separately fitted models by OLSSR, but the H, CL, and HCB were simultaneously estimated, and the correlation between them was interpreted clearly by the variance-covariance matrix.

The simultaneous MS (Equation (19)) residual plots (Figure 4) of fitting data against predictors are shown above; no obvious heteroscedasticity occurred in Figure 3.

3.4. Additive H, CL, and HCB Model System Estimated with NSUR

The sum of CL and HCB should be H in order to ensure the additivity of H MS (Equation (20)). We adopted the additive MS containing three equations (Equation (20)), with constraints on H-CL-HCB structural parameters and correlations of the error terms among the equations.

The parameter estimates by NSUR for the additive Equation (20) are listed in

Table 5. Parameter estimates of the systems of H, CL, and HCB (Equation (20)) (H: total tree height; CL: crown length; HCB: height to live crown base; values in the parenthesis are standard errors).

Compatible Additive Model System (Equation (20))
H Model
CL Model		HCB Model
Parameter	Estimates	Parameter	Estimates
b0	11.6479 (0.1146)	c0	0.3462 (0.0342)
b01	2.6157 (0.0983)	c1	0.5571 (0.0302)
b1	7.4269 (0.7177)	c11	0.1648 (0.0073)
b2	−0.1503 (0.0035)	c2	−2.3110 (0.0998)
b3	0.0085 (0.0091)	c3	1.1529 (0.0484)
b4	−0.0177 (0.0054)	c4	−0.0358 (0.0186)
b5	−0.0075 (0.0016)	c5	0.0597 (0.0023)
b6	−0.0389 (0.0035)

. The parameters were highly significant (p < 0.05), except for b₃ (p = 0.2289) for the CL sub-model of Equation (20). The estimated parameters of dummy variables b₀₁ and c₁₁ were also significant (p < 0.001). For the CL model, the effect of D, Dg, and HCB positively affect the asymptote, while SDI and CW have a negative effect. For the HCB model, the effect of CW, stocking (BAS), and DT negatively affect the asymptote, while D, Dg, and SDI have a positive effect. In terms of H, it is hard to interpret how these variables affect H. As H is the sum of CL and HCB, the comprehensive effect of CL and HCB will finally influence the H. All of these parameters are biologically plausible.

All statistics for the additive MS of H, CL, and HCB (Equation (20)) apparently possessed much higher accuracy and much lower errors than those for the sub-model of H, and HCB for the model of Equation (19), except for CL, which occurred the reverse condition (

Table 6. Evaluation indices of Equation (20) (H: total tree height; CL: crown length; HCB: height to live crown base; R²: coefficient of determination; $\overline{e}$: bias; $\delta$: bias variance; RMSE: root mean square error; TRE: total relative error).

System	Sub-Model	R2	$$\overline{\mathit{e}}$$	$$\mathit{\delta }$$	RMSE	TRE
Additive model (Equation (20))	H	0.7749	−0.0194	8.1403	2.8267	3.8886
	CL	0.5250	0. 0214	8.5295	2.9206	11.9447
	HCB	0.3430	−0.0407	8.3631	2.8922	22.8642

), i.e., the values of R², and δ, RMSE, and TRE for the CL model in the MS of the model (20) were 11.31% smaller and 16.66%, 7.93%, and 14.83% larger than those from the CL model in simultaneous MS (19), respectively. The $\overline{e}$ occurred in contrast to other error indices. Overall, the improvement of MS (20) was greater than in the model system (19).

$$\{\begin{array}{l}H=\frac{b_0+b_{01}SP}{(1+b_1\exp (b_{2}D+b_{3}CW+b_4{D}_g-b_{5}SDI+b_{6}HCB))}\\ +\left(-c_{0}CW+\left(c_1+c_{11}SP\right)D_g+c_{2}BAS+c_{3}SDI+c_{4}DT\right)\left(1-\exp (-c_{5}D)\right)+{\epsilon }_H\\ \\ CL=\frac{b_0+b_{01}SP}{(1+b_1\exp (b_{2}D+b_{3}CW+b_4{D}_g-b_{5}SDI+b_{6}HCB))}+{\epsilon }_{CL}\\ \\ HCB=\left(-c_{0}CW+\left(c_1+c_{11}SP\right)D_g+c_{2}BAS+c_{3}SDI+c_{4}DT\right)\left(1-\exp (-c_{5}D)\right)+{\epsilon }_{HCB}\\ \\ \epsilon ={\epsilon }_H+{\epsilon }_{CL}+{\epsilon }_{HCB}\end{array}$$

(20)

The estimates of the variance-covariance matrix R (Equation (4)) in the estimated MS of the additive H, CL, and HCB models by NSUR are shown below:

$$\sum =\left(\begin{array}{lll}8.15& 4.15& 3.99\\ 4.15& 8.53& -4.38\\ 3.99& -4.38& 8.37\end{array}\right)$$

for the system of the H, CL, and HCB models (20).

All the elements for Equation (20) were also significant (p < 0.05); the difference between Equations (19) and (20) lies in the values of these elements of the variance-covariance matrix. The values in Equation (20) were smaller than those of Equation (19), except that the covariance of CL and HCB was increased, which may be caused by the correlation between model CL and HCB being enhanced when HCB was taken as a predictor in the CL sub-model. The overall variance and covariance were weakened, which may be caused by the shared parameters of CL and HCB with H.

The additive MS (Equation (20)) residuals plots (Figure 5) of fitting data against predictors were shown above; no obvious heteroscedasticity occurred in Figure 4.

3.5. Comparison of Compatible and Additive Model Systems

Both the MS (compatible and additive model systems) were tested against the 30% dataset allocated for model validation. Validation statistics R² and $\overline{e}$, $\delta$, RMSE, and TRE of the additive MS estimated with NSUR were superior to the simultaneous MS estimated with NSUR (

Table 7. Evaluation indices of Equations (19) and (20) based on the validation data (H: total tree height; CL: crown length; HCB: height to live crown base; R²: coefficient of determination; $\overline{e}$: bias; $\delta$ : bias variance; RMSE: root mean square error; TRE: total relative error).

System	Sub-Model	R2	$$\overline{\mathit{e}}$$	$$\mathit{\delta }$$	RMSE	TRE
Simultaneous model system (Equation (19))	H	0.6382	0.9486	11.3382	3.4983	6.9084
	CL	0.5716	0.1987	7.5201	2.7495	11.3537
	HCB	0.2906	−0.1787	8.8704	2.9837	22.9350
Compatible additive model system (Equation (20))	H	0.7531	−0.0877	8.3422	2.8896	4.0012
	CL	0.4847	0.2431	9.0334	3.0153	13.7654
	HCB	0.3172	−0.3309	8.4588	2.9271	21.0731

).

The boxplot of residuals of Equations (19) and (20) (Figure 5) showed slight heteroscedasticity in the HCB sub-models, no matter in which model systems, compared to the residuals of the H and CL models without heteroscedasticity.

The range of prediction residuals (Figure 6) produced with the simultaneous MS was wider than those of the additive MS in terms of H and HCB and the ability to predict the CL of the simultaneous MS prior to the compatible additive MS.

The prediction of Equation (19) (Figure 7) and Equation (20) (Figure 8) are shown above; the scatterplots of prediction with validation data are differentiated between two types of forest types.

The residual VS predictors of Equation (19) (Figure 9) and Equation (20) (Figure 10) are shown above; no obvious heteroscedasticity occurred in Figure 9 and Figure 10; the fitting efficiency is great.

4. Discussion

In the previous studies [14,24,42,50], the simultaneous MS and additive MS almost always revolved around individual tree and stand biomass (stem, branch, foliage, bark, root, and total individual tree biomass), where the sum of the components’ biomass equals to the total biomass or the sum of four direction radii equals to twice the crown width. Biomass additivity and crown additivity are all about the additivity; their compatibility has been seldom mentioned. The NSUR is commonly used to solve the compatibility and additivity of these problems. Moreover, AP (adjustment in proportion) is also a common method to estimate biomass systems and crown radii systems. These methods [31,52,53,54,55,56] all ensure the additivity and compatibility of models. i.e., Fu et al. [8] developed a system of nonlinear additive crown width models utilizing seemingly unrelated regression for Prince Rupprecht larch in northern China. The compatibility and additivity of H, CL, and HCB is a desirable aspect for a system of models used for predicting H, CL, and HCB with these predictors, including individual tree-level and stand-level variables. So far, however, no study has analyzed the compatibility and additivity properties of H, CL, and HCB. In this study, the simultaneous MS and additive MS of H, CL, and HCB were developed with the most widely applied estimation method, NSUR, due to its superiority to other methods (FIML, 2SLS, 3SLS, and OLS) as NSUR accounts for the inherent correlations between H, CL, and HCB. The additive MS was superior to the simultaneous MS. It may be because the parameters in additive MS are shared by H with CL, and HCB that the parameter estimation by the variance-covariance matrix R (Equation (5)) is simpler than that of the simultaneous MS.

The base models of H, CL, and HCB are of great importance for the simultaneous MS of H, CL, and HCB model development. Albert [57] developed a generalized nonlinear mixed-effects height–diameter model for Norway spruce in mixed-uneven aged stands, and 21 (H-D) functions were evaluated for their fit performance; (Fu et al.) [25,32] developed a crown additive model fitted by NSUR and a generalized nonlinear mixed-effects height-to-crown base model with a logistics model; Yang [58] developed a generalized nonlinear mixed-effects height-to-crown base model for Picea crassifolia trees in northwest China. In this study, due to the data attribute and the model attribute, when several forms of models would be simultaneously estimated, the convergence would be difficult to meet. The monomolecular model has two parameters; in addition, the fitting accuracy is relatively high and easily estimated when estimated together with CL in the logistics model with many experiments; finally, the monomolecular model was the best base model for H and HCB, and the logistic model for CL was selected, finally, as the base model. The accuracy of these models was definitely enhanced, and the errors were reduced to a great extent when the tree-level and stand-level variables were introduced as such variables significantly affect the growth of H, CL, and HCB. The innovation of this study differs from previous studies: first, the developed simultaneous MS and developed additive MS of H, CL, and HCB; second, the model forms change—the allometric growth biomass models are the best choice for biomass, but the base models of this study need to be tested by many model forms because the variation discrepancy between H, CL, HCB and biomass is enormous.

Our study includes two types of forest stands; the growth differences between the two types of forest stands are native, so the introduction of dummy variables will be reasonable. The accuracy was improved for every sub-model when a dummy variable was introduced. The dummy variable model parameters were all significant (p < 0.05) except for b₃ (p = 0.2289) for CL sub-model of Equation (20), indicating the discrepancy of forest two-kinds stand types is significant. The reason why only several individual tree- and stand-level variables were introduced into our base model was that introducing too many predictors into the systems could impede the computation convergence and result in over-parameterization that causes biased parameter estimates, as reported by Tang [59]. Furthermore, including too many predictor variables will increase forest inventory costs. Therefore, D, CW, HCB, Dg, and SDI are appropriate as covariate predictors in the system of the H, CL, and HCB models.

The simultaneous MS of H, CL, and HCB and additive MS could be used for estimating part of the values of H, CL, and HCB of individual trees for Mongolian oak and spruce-fir in northeast China, where H, CL, and HCB are not simultaneously known because of some inconvenient reasons of measuring these predictors in regular forest inventories. D and CW can be measured to predict H by calculating the DT, Dg, and BAS of the stand. Then, D, DT, Dg, BAS, and SDI can be used to predict HCB, and, finally, CL will be calculated. The simultaneous MS of H, CL, and HCB models will greatly reduce the cost and time for the collection of required input variables in the field. The shortcoming of our model lies in the dearth of H as a predictor of HCB, according to many studies, because H is admitted as a robust variable to predict HCB. Therefore, the fitting accuracy of our HCB sub-model was comparatively lower than the other two models of H and CL. When it comes to the additive MS, H cannot be a covariate in the model. Hence, H would be excluded.

NSUR and OLSSR were utilized to estimate the parameters of the simultaneous MS of H, CL, and HCB, respectively, for simultaneously compatible model MS and additive MS. When NSUR was used, the correlation between H, CL, and HCB was taken into consideration by the variance-covariance matrix; the variance-covariance matrix is an excellent interpretation of their inner relationships. When OLSSR was used, the fitting accuracy was higher than NSUR, but it could not estimate H, CL, and HCB simultaneously. Compared with NSUR, the OLSSR never considers the correlation of sub-models, so the results produced from OLSSR would usually be biased, while NSUR considers the deeper relations of H, CL, and HCB and interprets the relations by R matrix, so the results from NSUR are unbiased. There will be more studies on applying the TSEM (two-stage error-in-variable method) and other methods to the simultaneous system of the H, CL, and HCB models in order to offer more efficient methods to solve these problems. We also tried fitting simultaneous and additive model systems using FIML, 2SLS, and 3SLS; however, none of the model systems converged. The methods of using FIML, 2SLS, and 3SLS to make the simultaneous and additive model systems converge will be our next study point.

5. Conclusions

In this study, NSUR and OLSSR approaches were used to develop the simultaneous MS of H, CL, and HCB models. The data was collected from two types of forest stand; the dummy variables were introduced as indicator variables to explain the discrepancy of the two forest types. The variance-covariance matrix of NSUR was used to interpret the correlations among H, CL, and HCB, which ensured the compatibility and additivity of these sub-models, while OLSSR failed to do it. The logistic model could be used as the base model for estimating CL, and the monomolecular model was the best base model for H and HCB by multiple tests. The logistics model may have been the better choice when CL was estimated in this study. Compared with OLSSR and other methods, NSUR was used widely and is convenient for coming through to the additivity and compatibility of stand and individual biomass and crown components. The compatibility and additivity of H, CL, and HCB can be solved by NSUR; the additivity MS had better performance than simultaneous MS. NSUR would be the best choice when the compatible MS and additive MS of several models need developing. The priority of FIML, 2SLS, and 3SLS had not been found compared to NSUR in this study, but whether these methods have their own advantages over other problems needs to be studied deeply.