Variable-Exponent Taper Equation Based on Multilevel Nonlinear Mixed Effect for Chinese Fir in China

: A variable-exponent taper equation was developed for Chinese ﬁr ( Cunninghamia lanceolate (Lamb.) Hook.) trees grown in southern China. Thirty taper equations from different groups of models (single, segmented, or variable-exponent taper equation) were compared to ﬁnd the excellent basic model with S-plus software. The lowest Akaike information criteria (AIC), Bayesian information criteria (BIC), and -2loglikelihood (-2LL) was chosen to determine the best combination of random parameters. Single taper models were found having the lowest precision, and the variable-exponent taper equations had higher precision than the segmented taper equations. Four variable-exponent taper models that developed by Zeng and Liao, Bi, Kozak, Sharma, and Zhang respectively, were selected as basic model and had no difference in ﬁt statistics between them. Compared with the model without seldom parameter, the nonlinear mixed-effects (NLME) model improves the ﬁtting performance. The plot-level NLME model was found not to remove the residual autocorrelation. The tree-level and two-level NLME model had better simulation accuracy than the plot-level NLME model, and there were no signiﬁcant differences between the tree-level and two-level NLME model. Variable-exponent taper model developed by Kozak showed the best performance while considering two-level or tree-level NLME model, and produced better predictions for medium stems compared to lower and upper stems.


Introduction
Stem profile equations, commonly known as taper equations, which refer to the general decrease in the regular outline of a solid body from its base to its tip. The taper equation is important mainly because it can predict the follows: (1) diameter at any height of the tree stem or height at any diameter; (2) volume at any height of the tree stem or total volume of the tree stem and total peeling volume; and (3) merchantable timber volume and dry timber height with any diameter restriction or relative height restriction [1]. The taper equation can be generally divided into three categories in the history of development. A single taper equation [2][3][4][5] is composed of simple mathematical equations, and it has disadvantage that the significantly swollen part at the bottom of the tree stem is inconsistent with the model fitting. The segmented taper equation [6][7][8][9] refers that several polynomials representing different parts of the tree stem are connected through the inflection point, and the tree stem is assumed as concave, paraboloid and conical respectively from bottom to top. However, more complex iterative operation as well as more fitting data and experience are required for the estimation of segmented taper equation inflection point. The variableexponent taper equation [10][11][12][13] refers that the tree stem shape can be better estimated through the change of the independent variable exponent in the continuous function, which can be applied for theoretically describing the tree stem of any shape. Variable-exponent taper equation shows better imitative effect and application prospect [10,[13][14][15] than other two taper equations.
Taper equation fitting data are often based on repeated measurements of diameters in the same tree at different heights. A hierarchical structure is available. There is a correlation between measured data for many times. The following assumption of error term independent identical distribution of least square method is violated, thereby possibly leading to deviation in variable parameter estimation [16]. Mixed effect model has been applied in the research of the taper equation in recent years [17,18]. Nonlinear mixed effect model is an effective way of processing stratified data and repeated measurement data, including fixed parameters and random parameters, which can meet the assumption of error term independent identical distribution. Randomness and elasticity due to known or unknown factors can be explained [19]. A nonlinear mixed model method has been utilized for fitting Max and Burkhart segmented taper equation, and the results shows that mixed effect model can improve the fitting precision of the model [20].Chinese fir is a native coniferous timber tree species in China with a long cultivation history, which is widely distributed in 18 provinces and regions in southern China. Chinese fir is favored by people because of its rapid growth, excellent material, and straight and full stem forms [21]. The result of the national continuous forest inventory shows that [22] the area of Chinese fir man-made forest has reached 9.90 million hectares, the standing volume has reached 755 million cubic meters, and they respectively account for 27.23% and 32.57% of the main dominant tree species of artificial forest in China. Chinese fir has become the most important timber plantation tree species in China. It is lack of systematic comparative analysis among many basic models while building Chinese fir taper equation. Mixed effect model is rarely used for constructing variable-exponent taper equation particularly, while the autocorrelation of the fitting data is not deeply considered.
The objective of the present study was to evaluate the performance of some wellknown taper functions and develop a mixed-effects variable-exponent taper equation for Chinese fir trees in southern China.

Data
A total of 183 trees sampled from even-aged Chinese fir stands were used in the present study. The trees were taken from stands located in Jiangxi province of southern China, and included the typical range of diameters and heights of Chinese fir stands. The total tree height (H: m) and diameter at breast height (DBH: 1.3 m above ground level) were measured. Diameter outside bark (cm) was also measured at heights of 0.2, 1, 1.3, and 2 m and then at intervals of 1 m along the remainder of the stem. A total of 123 trees selected randomly were used as fit data, and the rest of the 60 trees were used as validation data. The tree height and DBH data statistics were shown in Table 1. Figure 1 shows the trend of the relative diameter along every tree stem with the relative height. The relative height is the ratio of the height of the off-ground tree stem and total tree height. The relative diameter is the ratio of the diameter at the height of the off-ground tree stem to DBH.

Establishment of Mixed Effect Model
Mixed effect model contained both fixed parameters and random parameters. Random parameters were fixed parameters added to nonlinear model, which varied with changes of different blocks. A mixed effect model was established according to plot-level effect, tree-level effect and two-level effect respectively in the study. The model expression is as follows: where Y ijk refers to the k th observed value of diameter d at the height h of the j th tree in the i th plot. f (·) refers to differentiable function containing parameter vector Φ ijk and concomitant variable vector t ijk , ε jik refers to a error term obeying normal distribution, it is assumed that the expected value is zero with positive variance covariance structure R ij , wherein Φ ijk can be expressed as follows: where β refers to n × 1-dimension fixed parameters vector, u i and u ij , respectively, refer to random parameter vectors of plot level and tree level, which obey independent normal distribution, and expected value is zero with variance covariance structures ψ 1 and ψ 2 . A ijk , B i,jk and C ijk refer to design matrix.

Determination of Random Parameters
One key for establishing a mixed effect model lies in determination of fixed parameters and random parameters in the model. Too many random parameters may lead to excessive parameterization or convergence [42], The possibility of 1-2 random parameter combinations in the selected basic model was fitted. The optimal combination of fixed parameters and random parameters was determined by comparing the minimum values of statistics AIC (Akaike information criteria), BIC (Bayesian information criteria), and -2LL(-2loglikelihood).

Variance-Covariance Structure of Stochastic Effect
Variance covariance structure can reflect the changes among plots and among trees. Two random parameters are adopted as examples according to the study of Calama [43], variance covariance structure can be set as a positive definite structure matrix, and the structure is shown as follows: where σ 2 i (i = 1, 2) refers to the variance of random parameters, σ ij (i, j = 1, 2, i = j) refers to the covariance of random parameters 1, 2.

Intra-Group Variance-Covariance Structure
The intra-group variance-covariance structure should be determined to eliminate the correlation of intra-group error. The previous studies show that the mixed effect model with an unstructured covariance-variance matrix can be applied for explaining the autocorrelation among most observed values [16,44]. The unstructured variancecovariance matrix of random effects was used in this study. For the random error term, constant variance and independence of observations was assumed because the emphasis is on prediction rather than hypothesis testing.

Model Evaluation
The prediction of random-effects parameters was accomplished by an approximate Bayes estimator [41]. Statistical indicators for evaluating the model include adjusted determination coefficient (R 2 adj ), mean variation (Bias) and root-mean-square error (RMSE). These expressions are shown as follows: where y i refers to the observed value,ŷ i refers to the predicted value, y i refers to the average value of the observed value, n refers to the number of observed values, and p refers to the model parameter number.

Selection of Basic Model
After estimating the fit statistics of each one of 30 models fitted, six models were selected for further analysis. In general, the variable-exponent taper models showed the best results; therefore, the four best variable-exponent taper models (Zeng and Liao [37]; Bi [12]; Kozak [1] and Sharma and Zhang [13]), the best single taper function (Cervera [26]) and the best segmented model (Max and Burkhart [6]) were selected. The parameter estimates and fit statistics of each model are shown in Table 3.

Construction of Mixed Effect Model
A mixed effect model was built on the basis of four selected variable-exponent taper models in the study. The possibility of 1-2 random parameter combinations in the selected basic model was fitted. Wherein, Zeng and Liao [37] had 10 pairs of combinations. Bi [12] had 28 pairs of combinations. Sharma and Zhang [13] had 10 pairs of combinations. Kozak [1] had 45 pairs of combinations. The parameter combination of AIC, BIC, and -2LL with the minimum statistical indicators was selected as the optimal mixed effect model. The fitting results are shown in Table 4. It can be seen from Table 4 that the AIC value based on the plot-level effect mixed model is decreased by 3.3-5.8% compared with the basic model aiming at four variable-exponent taper equations. The AIC value based on tree-level effect mixed model is reduced by 40.7-46.6% compared with the basic model, and the AIC value nested with two level effect mixed models is decreased by 40.6-46.5% compared with the basic model. BIC and -2LL values had similar decrease rate. Wherein, AIC, BIC and -2LL values of the model Bi [12] in the basic model and the mixed model with plot-level effect were the minimum. The AIC, BIC, and 2LL values of model Kozak [1] in the mixed model with tree-level effect and two-level effects were the lowest. Figure 2 is the lag residual diagram of the mixed effect model added with intra-group variance-covariance structure of Kozak [1] model after determination of mixed parameters. A, B, C, and D are model lagged residual diagrams based on common nonlinear least square method, plot-level effect, tree-level effect, and nested two-level effect mixed effect respectively in the figure. It can be seen that the common least square method and mixed models based on plot-level effect had obvious positive correlation. The residual diagrams of mixed model based on tree-level effect and nested two-level effect had no prominent trend. It is obvious that the structure-free intra-group variance-covariance matrix mixed effect model can effectively eliminate autocorrelation.  Figure 2 is the lag residual diagram of the mixed effect model added wi tra-group variance-covariance structure of Kozak [1] model after determination of parameters. A, B, C, and D are model lagged residual diagrams based on common linear least square method, plot-level effect, tree-level effect, and nested two-level mixed effect respectively in the figure. It can be seen that the common least s method and mixed models based on plot-level effect had obvious positive corre The residual diagrams of mixed model based on tree-level effect and nested two effect had no prominent trend. It is obvious that the structure-free intra-group var covariance matrix mixed effect model can effectively eliminate autocorrelation.  Table 5 respectively shows fixed parameter and random parameter estimated v of four variable-exponent taper equations based on tree-level effect mixed model. be found that all parameters are statistically significant. Table 6 shows the fit sta comparison of four variable-exponent taper equations in the traditional least s method and those based on plot-level effect, tree-level effect, and nested two-level The simulation precision in mixed mode based on plot, tree-level effect, or n two-level effect is higher than that of the traditional least square method. The simu precision based on tree-level effect and nested two-level effect is higher than that on plot-level effect mixed model. The determination coefficient of mixed effect based on plot-level effect is higher than that of the basic model by 0.0016-0.002 determination coefficient of mixed model based on tree-level effect and nested two    Table 5 respectively shows fixed parameter and random parameter estimated values of four variable-exponent taper equations based on tree-level effect mixed model. It can be found that all parameters are statistically significant. Table 6 shows the fit statistics comparison of four variable-exponent taper equations in the traditional least square method and those based on plot-level effect, tree-level effect, and nested two-level effect. The simulation precision in mixed mode based on plot, tree-level effect, or nested two-level effect is higher than that of the traditional least square method. The simulation precision based on tree-level effect and nested two-level effect is higher than that based on plot-level effect mixed model. The determination coefficient of mixed effect model based on plot-level effect is higher than that of the basic model by 0.0016-0.0020. The determination coefficient of mixed model based on tree-level effect and nested two-level is higher than that of the mixed model based on plot-level effect by 0.0104-0.0117. The simulation precision of the mixed model based on the tree-level effect and nested two-level had no prominent difference. Bi [12] model has the highest precision in the least square method basic model and the mixed model of nested two-level. Kozak [1] model has the highest determination coefficient and the lowest root-mean-square error in the mixed model based on tree-level effect and nested two-level effect. It is consistent with the goodness of fit indicators AIC, BIC, and -2LL. Table 5. Parameter estimates and variance components of the tree levels nonlinear mixed-effects (NLME) (standard errors in parentheses).

Predicted Diameters and Model Calibration
The established Kozak [1] mixed effect model based on tree-level effect was used for fitting the validation data to test the model's prediction ability. The relative height of each tree in validation data was divided into 10 parts, and the prediction mean variation and SD of the diameter in each part was segmented for statistics. The result was shown in Table 7. It can be seen from Table 7 that the mean variation of the Kozak [1] mixed effect model based on the tree-level effect was smaller in the 10 segmented parts. The mean variation was slightly larger in the stem base (0.0 < h/H ≤ 0.1) and at the top of the stem (0.7 < h/H ≤ 1.0), both variations were larger than 0.8cm. However, the mean variation at the relative height of 0.3 < h/H ≤ 0.4 was the smallest, namely, reaching 0.0448 cm. The prediction accuracy of the Kozak [1] mixed effect model based on the tree-level effect is higher in the middle part of the trunk than that in the base and at the top of the stem, and SD also has the same trend. Figure 3

Discussion
Single taper equations have been found to describe stem shape quite accurat [2,45], although more flexible equations were developed to provide a better descript of the stem profile [6,10]. It was found that the model of Cervera [26] was optimal in gle taper equations and only slightly lower than the other two types of stem form eq tions. However, the single taper equations often fail to characterize the entire stem pro and produce bias, especially for the area near the butt and at the very top sections of t Taper profiles for three different diameter at breast height (DBH) and the similar height tree, three different height and the similar DBH tree in Kozak [1] with tree-level effects NLME model.

Discussion
Single taper equations have been found to describe stem shape quite accurately [2,45], although more flexible equations were developed to provide a better description of the stem profile [6,10]. It was found that the model of Cervera [26] was optimal in single taper equations and only slightly lower than the other two types of stem form equations. However, the single taper equations often fail to characterize the entire stem profile and produce bias, especially for the area near the butt and at the very top sections of tree [6,8]. The model of Max and Burkhart showed the same excellent accuracy as previous studies [6]. A variable-exponent taper equation describes the bole shape with a changing exponent or variable from ground to top to represent the neiloid, paraboloid, conic, and several intermediate forms [10]. This approach is based on the assumption that the stem form varies continuously along the length of a tree [39]. In comparison with single and segmented taper models, this approach provides the lowest bias and highest precision in taper prediction [10,11,14]. In this study, the fitting precision of the variable-exponent taper equation is higher than other two groups of taper equations, which is consistent with the conclusions of previous studies [10,13]. The four optimal variable-exponent taper equations including Zeng and Liao [37], Bi [12], Kozak [1], Sharma and Zhang [13] almost provide similar results in terms of the goodness of fit statistics.
Mixed models have also been used successfully in the development of taper equations for several species throughout the world [46][47][48][49][50], showing that estimates of the mean response (fixed-effects model) can be significantly improved by including random effects parameters. This was confirmed in the present study. Further, the fitting precision of the mixed model based on tree-level effect and nested two-level effect were proved to be higher than the mixed model based on plot-level effect. The mixed models based on nested two-level effect and plot-level effect had no significant difference. Kozak [1] model had the highest adjusted determination coefficient among the four optimal variable-exponent taper equations while considering tree-level effect and nested two-level effect. The inclusion of an intra-group variance-covariance matrix in the mixed effect model can eliminate the autocorrelation, which is consistent with the results of previous studies [16,44,51,52]. It can be found that the fitting residual scatterplot of taper model based on common least square method and plot-level effect shows significant positive correlation, while the residual diagram of mixed model based on tree-level effect and nested two-level effect both have removed this trend in a large part. Although the mixed model accounting for autocorrelation does not improve the prediction performance, it improves interpretation of the statistical properties [53].
In addition to this model and its simulation methods selection, some tree or stand factors (e.g., crown height, ratio, and planting density) that have deep relationship with tree form and quality [13,54] have been introduced into taper equations to improve modeling performance [8,13,55]. Additionally, for a taper equation, its aim is to develop a volume equation and calculate stand merchantable timber volume. The stem taper model should be integrable. Unfortunately, the model of Kozak [1] is not analytically integrable. Therefore, numerical integration methods and iterative procedures to estimate height at aspecific end diameter should be used.

Conclusions
The variable-exponent taper equations were developed for Chinese fir, the most important commercial tree species in southern China. Through the comparison of thirty taper equations, single taper models were found having the lowest precision, and the variable-exponent taper equations had higher precision than segmented taper equations. Variable-exponent taper models developed by Zeng and Liao [37], Bi [12], Kozak [1], Sharma and Zhang [13] almost provide similar results in terms of the goodness of fit statistics, and can be selected as taper models of Chinese fir trees. The fitting precision of the mixed model based on tree-level effect and nested two-level effect proved to be higher than the mixed model based on plot-level effect. Variable-exponent taper model developed by Kozak [1] produced better predictions for medium stems compared to lower and upper stems.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.