Biomass Modeling of Larch ( Larix spp . ) Plantations in China Based on the Mixed Model , Dummy Variable Model , and Bayesian Hierarchical Model

With the development of national-scale forest biomass monitoring work, accurate estimation of forest biomass on a large scale is becoming an important research topic in forestry. In this study, the stem wood, branches, stem bark, needles, roots and total biomass models for larch were developed at the regional level, using a general allometric equation, a dummy variable model, a mixed effects model, and a Bayesian hierarchical model, to select the most effective method for predicting large-scale forest biomass. Results showed total biomass of trees with the same diameter gradually decreased from southern to northern regions in China, except in the Hebei province. We found that the stem wood, branch, stem bark, needle, root, and total biomass model relationships were statistically significant (p-values < 0.01) for the general allometric equation, linear mixed model, dummy variable model, and Bayesian hierarchical model, but the linear mixed, dummy variable, and Bayesian hierarchical models showed better performance than the general allometric equation. An F-test also showed significant differences between the models. The R2 average values of the linear mixed model, dummy variable model, and Bayesian hierarchical model were higher than those of the general allometric equation by 0.007, 0.018, 0.015, 0.004, 0.09, and 0.117 for the total tree, root, stem wood, stem bark, branch, and needle models respectively. However, there were no significant differences between the linear mixed model, dummy variable model, and Bayesian hierarchical model. When the number of categories was increased, the linear mixed model and Bayesian hierarchical model were more flexible and applicable than the dummy variable model for the construction of regional biomass models.


Introduction
Forest biomass plays an important role in regulating global carbon balance to help mitigate the effects of climate change.With the development of large regional-scale biomass monitoring work, developing a biomass model for large-scale application has become an important job.In recent years, many researchers have explored biomass models at national, regional, and global levels [1][2][3][4].Allometric equations were often used to estimate forest biomass in these studies [1,[5][6][7][8].The allometric equations usually had good fit performance and high values of R 2 , and the variables for prediction of aboveground biomass, such as tree diameter and height, are easy to measure in the field.However, a main drawback of these equations is that they produce different results when applied to sites outside the geographic allocation where the equations were originally developed [6,9,10].Several authors have noted that the problem with the scaling coefficient in allometric equations is variation, which is dependent on species, tree stage, and site [7,9,[11][12][13].Case and Hall [3] also found that there was an increase in prediction error when using allometric equations to estimate tree biomass at different geographical scales.The crucial problem is how to make biomass models have a wider range of applicability, and not only meet the accuracy requirements for large-scale biomass estimates, but also improve local biomass prediction accuracy.The three models for estimating forest biomass, including the dummy variable model, the mixed effects model, and the Bayesian hierarchical model, provide possible solutions to this problem.
The mixed effects model is an improvement on the practical statistical techniques when estimating fixed effect and random effect parameters to reflect the difference between tree and plot, and to improve the precision of the prediction model [14].In recent years, the mixed effects model has been widely used in the field of forestry.Fang and Bailey [15] applied a modified Richards growth function with nonlinear mixed effects for modeling slash pine dominant height growth with different silvicultural treatments.Daniel and Michael [16] developed a multilevel nonlinear mixed effects model system for describing dominant height, basal area, stand density, and volume.Zhang and Borders [17] used the mixed effects model to improve the precision of tree biomass estimation in the intensively managed loblolly pine plantations.Trincado and Burkhart (2006) [18] explored a nonlinear mixed effects model to describe stem taper, which showed that the method increased their prediction flexibility and efficiency.Fehrmann et al. [19] compared the mixed model with k-nearest neighbor methods in the establishment of Norway spruce and Scots pine tree biomass models in Finland.The relative root mean square errors of linear mixed models and k-nearest neighbor methods estimates are slightly lower than those of an ordinary least squares regression model.The results show that nonparametric methods are suitable in the context of single-tree biomass estimation.Derek et al. [20] applied the mixed effects model to develop branch models for white spruce with hierarchical structure data.In short, mixed effects models have become important tools for modeling tree height, stem taper, branch diameter, crown, and tree biomass.The Bayesian method is usually used to evaluate ecological models [4,21].The results of Bayesian analyses with non-informative priors were similar to using parameters and statistics to fit allometric biomass equations with a classic statistical approach, however, the Bayesian method with informative priors performed better than the non-informative priors and the classic statistical approach [22][23][24][25].In recent years, the Bayesian method has been widely used in forestry, for applications such as diameter distribution [26,27], tree growth [28], individual tree mortality [29,30], and stand-level height and volume growth models [11].Zapata-Cuartas et al. [31] used Bayesian methods to estimate aboveground tree biomass, and obtained a reliable result.Zhang et al. [25] applied Bayesian methods to estimate Chinese fir tree biomass.When Bayesian methods are applied to quantify the forest biomass of a geographical location, the problems that occur often involve variation between different locations, which instead quantifies the hierarchical characteristics of the source data.The Bayesian hierarchical approach can resolve the problems in fitting the data.Zapata-Cuartas et al. [31] also concluded that the Bayesian hierarchical model may be an effective approach to improve biological estimation at a large regional scale.
In regression analysis, the dummy variable was considered as a virtual variable, which usually takes the values 0 or one [32].The dummy variable method is commonly used to deal with indicator or categorical variables, which are involved in all quantitative methods (Tang et al., 2008;Fu et al., 2012) [33,34].When the generalized model for a large-scale range is established, the dummy variable model is also a reliable model.Wang et al. [35] applied the mixed effects model and the dummy variable method to establish the dominant height model, and concluded that the two methods could be used to establish a specific parameter and obtained similar results.Fu et al. [34] used the linear mixed effects model and the dummy variable model to construct compatible single-tree biomass equations at different scales for Masson pine (Pinus massoniana) in southern China, and found that there were no significant differences between the results obtained from the two models.
Larch (Larix spp.), a fast-growing deciduous coniferous tree, is one of the main tree species used for timber production in China.Because it has a straight stem, better quality timber, and high commercial value, the planting area and volume for larch plantations in China are approximately 3.14 million ha, 18.4 million m 3 , respectively [36].Currently, the larch plantations, not only provide a large volume of commercial timber, but also play an important role in forest carbon fixation.Therefore, the goal of this research was to build larch stem wood, branches, stem bark, needle, root, and total biomass models at the regional level, using a general allometric equation, the dummy variable model, the mixed effects model, and the Bayesian hierarchical model, to select the most effective method for predicting large-scale tree biomass.

Site Description
The data were collected from six main growth regions of larch plantations in China (Table 1, Figure 1).The experiment sites in this study covered the main timber production regions of larch plantations in China.
commercial value, the planting area and volume for larch plantations in China are approximately 3.14 million ha, 18.4 million m 3 , respectively [36].Currently, the larch plantations, not only provide a large volume of commercial timber, but also play an important role in forest carbon fixation.Therefore, the goal of this research was to build larch stem wood, branches, stem bark, needle, root, and total biomass models at the regional level, using a general allometric equation, the dummy variable model, the mixed effects model, and the Bayesian hierarchical model, to select the most effective method for predicting large-scale tree biomass.

Site Description
The data were collected from six main growth regions of larch plantations in China (Table 1, Figure 1).The experiment sites in this study covered the main timber production regions of larch plantations in China.

Biomass Data
The biomass data of 360 sample trees from 209 plots were used in this study, which were collected via destructive sampling from Larch plantations in China between 2009 and 2013.The sample plots were located in Hubei, Gansu, Hebei, Liaoning, and Heilongjiang provinces and the Inner Mongolia autonomous region (Figure 1).The data covered the complete period of stand development for timber production, with stand age ranging from 3 to 45 years.The stem diameter at breast height (DBH), total height (HT), height to crown base (HCB), and crown width (CW) were recorded for each tree within the plots.The sample trees were chosen according to the dominant, average, and inferior tree outside the plot in the stand.Then the sample trees were felled as carefully as possible to minimize damage to their crowns.The tree height, live crown length, and height of lowest living branch of the trees were recorded, and subsequently the trees were stripped of their branches.The live crown was divided into three parts: top, middle, and bottom.From each part, two branches were selected for further sub-sampling in order to determine fresh and oven-dried biomass ratio.For about half of the sub-samples, shoots were further sub-sampled by removing all needles from the branches in order to determine needle to branch biomass ratio.All of the branches were removed and the tree length was recorded.Each stem was then divided into 1m sections.The fresh weight of each stem section was measured.From each section, a stem disc (2 cm wide) was taken as a sub-sample in order to determine the fresh to oven-dried biomass ratio for each stem section.The bark of each stem disc was removed and the bark to wood mass ratio was determined in the lab.The whole roots were manually excavated, and were then sorted into root stump, large roots (≥5 cm), medium roots (2-5 cm) and small roots (≤2 cm).The aboveground stumps that remained after tree felling were included in the root stump biomass.The fresh biomass of each root diameter class was determined in the field.Sub-samples of each class were taken and their fresh biomass was determined with an electronic balance.All sub-samples of branches, needles, stem wood, stem bark, and roots were oven-dried at 80 • C until a constant weight was reached (usually within 24-48 h).According to the ratios of fresh to oven-dried biomass, the total branch, needle, stem wood, stem bark, and root biomass were calculated.The biomass data of each component in different regions was summarized in Table 2.

General Model
Allometric equations were widely used for estimating biomass [1,2,5,7].The form of the equation provides a good balance between accurate predictions and low data requirements by using the most commonly and easily measured variable, which is DBH.We also used the following allometric equation as a base model to construct different biomass equations in this study: where Y is the oven-dry weight of the biomass component of a tree, DBH is the diameter at breast height in cm, δ represents the multiplicative error term, and a and b are parameters.Equation ( 1) can be changed into the following linear form by logarithmic transformation: where: y = lnY, x = lnDBH, a 0 = lna, ε = ln δ.

Dummy Variable Model
The general form of the dummy variable model based on Equation ( 2) is expressed as follows: where z i is the dummy variable, a i is the corresponding specific parameter of region i, other parameters are the same as in Equation (2).

Linear Mixed Effects Model
The general form of the linear mixed effects model is as follows [14,32,33]: where y i is the n i -dimensional response vector for the ith region; β is the p-dimensional parameter vector for fixed-effects; X i is the n i × p design matrix for fixed-effects variables; u i is the q-dimensional region-specific vector of parameters for random effects (only one random effect); Z i is the n i × q design matrix for random-effect variables (random effect of the ith region); ε i is the n i -dimensional error vector.It is assumed that: Further, it is common to assume that ε i and u i are normally distributed and thus:

Bayesian Hierarchical Model
Equation ( 2) could be changed into the following form when the Bayesian hierarchical approach is used: where u 0i , u 1i are cluster-specific random variables to be predicted and assumed to be N 0, σ 2 0 , and N 0, σ 2 1 , respectively.ε i is assumed to be N 0, σ 2 .Variance, σ 2 0 or σ 2 1 , measures the between-subject variability, while σ 2 accounts for the within-subject variability for all the regions.
The relevant prior knowledge about the data can be incorporated into Bayesian analyses, which contrast the classical statistical approach [37].In the study, we need to choose appropriate prior distributions for all parameters, including a and b.The prior distributions of a and b (total tree, root, stem wood, stem bark, branch, and needle) were obtained for 36 biomass equations from six Chinese larch publications (Appendix A Table A1) [38][39][40][41][42][43].We assume the parameters a and b are distributed as a bivariate normal distribution.The Bayesian hierarchical model was fitted using the MCMC method in the R package R2WinBUGS [44].

Model Fitting and Evaluating
The general model, dummy variable model, linear mixed model, and Bayesian hierarchical approach were used to establish the stem wood, stem bark, branch, needle, root, and total biomass models.All the models were fitted using the "R3.2.1" software.The following statistics were calculated to compare the effect of the region-level biomass for larch plantations: prediction determination coefficient (R 2 ), root mean square error (RMSE), and absolute bias (MAB).
where y i is the observed biomass values, is the arithmetic mean of all observed biomass values, ŷi is the estimated biomass values based on models, and n is the sample number.

Results
In this study, the general biomass model, Equation (2), the dummy variable model, Equation (3), linear mixed effects model, Equation (4), and the Bayesian hierarchical model, Equation (5), were used to fit the biomass data on 360 sample larch trees from six regions (S1, S2, S3, S4, S5, S6) in China, developing the model of total tree, root, stem wood, stem bark, branch, and needle biomass respectively.Parameter estimates of the models, Equations (2)-( 5) are presented in Tables 3-6.The parameters of each model were statistically significant (p < 0.01), and the parameters of the dummy variable model, linear mixed effects model, and Bayesian hierarchical model may predict the biomass of different regions.The residual distributions of the dummy variable model, linear mixed effects model, and Bayesian hierarchical model were better than the general biomass model, and they were more uniform and had no obvious residual trends (Figure 2).The special and random parameters of region type indicate that the biomass component of a tree in six regions are obviously different.In addition, the biomass of a tree with the same diameter gradually decreased from southern to northern regions in China, except in region S3, Weichang County in Hebei Province, in a temperate region.However, the differences were not statistically significant.For example, when the tree diameter at breast height was 25 cm, the total biomass of single tree was 356 kg on the Changlinggang farm in Hubei (S1) in the northern subtropical region; 333 kg in Tianshui City in Gansu (S2) in the warm-temperate region; 293 kg on the Dagujia farm in Liaoning (S4) in a temperate region; 286 kg on the Mengjiagang farm in Heilongjiang (S5) and 256 kg in the Wuerqihan Forestry Bureau in the Inner Mongolia autonomous region (S6), which are both located in a cold-temperate region; and 239 kg in Hebei (S3) (Figure 3).In terms of the three fit statistics (R 2 , MAB and RMSE) for the general biomass model (Equation ( 2)), the dummy variable model (Equation (3)), the linear mixed effects model (Equation ( 4)) and the Bayesian hierarchical model (Equation ( 5)) are presented in Table 7.It was found that the R 2 values of Equations ( 3)- (5) were larger than Equation (2), and MAB and RMSE were lower than Equation (2), when they were used to estimate stem wood, stem bark, branch, needle, root, and total biomass.This indicated that the dummy variable model, mixed effect model, and Bayesian hierarchical model were better in the prediction of biomass than the general biomass model.Additionally, the fitting optimization index of the F-test also verified the results (Table 7).However, there were no significant differences among the results of fitting for the mixed effect model, dummy variable model, and Bayesian hierarchical model.In other words, the accuracy of prediction of biomass using the mixed effects model, dummy variable model, and Bayesian hierarchical model were similar.For the three models, the Bayesian hierarchical model was slightly better for predicting stem wood biomass.The predicted results of the biomass component indicated that the generalized R 2 is highest for stem wood, bark, root, total biomass, and lowest for needle and branch biomass.This may be due to the difficulties with available sample crown.In short, in this study, each component biomass model and total biomass model were well developed in different regions.Based on RMSE, MAB, and R 2 values, it could be concluded that the mixed model, dummy variable model, and Bayesian hierarchical model performed better than the general model, and the results of the three methods were approximate in their ability to describing the fit of the data.

Discussion
Larch is one of the most important afforestation tree species in China.The accurate estimation of larch biomass and carbon stock is critical.Allometric equations have been widely used to estimate forest biomass in many studies, but the use of allometric equations in regions outside the area in which they were developed or for different species, has been strongly debated in forestry literature.Some authors have even recommended using general allometric equations [12,13].In this study, the scaling coefficient in allometric equations was varied in different biomass components (stem wood, branches, stem bark, needle, root and total biomass), and across different larch species and regions.Although it is generally known that the scaling exponent b is between two and three [12], it is also usually realized as a species-specific equation, with different coefficients a and b, because trees may differ in architecture as well as wood density [6,45].West et al. [46] applied a process-based model which is called the WBE model, to estimate values of scaling exponents using a functional relationship; it indicated that the aboveground biomass of a tree species should scale against stem diameter with b = 8/3 (2.67), regardless of species, site, and age.However, Zianis and Mencuccini [12] estimated another empirical exponent (b = 2.36) based on a worldwide list of 279 biomass allometric equations, which is different from the scaling exponent (b ≈ 2.67).They indicated that the ratio of biomass and DBH for trees growing in different environmental conditions cannot be constant.In the study, the average value of parameter b was 2.42 for the total biomass model, and it did not differ from the theoretical value (2.67) and empirical value (2.36), because the scaling exponent (b = 2.67 or 2.36) is the average of all species.For larch plantation, the b value (2.42) may be more appropriate.
The larch data was collected from the six main ecological regions in China for larch plantation distribution and growth.The biomass equations suitable for application to regional forest biomass prediction are the linear mixed model, dummy variable model, and Bayesian hierarchical model.Based on the three evaluation statistics (R 2 , RMSE, MAB), the four models were compared: the general biomass model (Equation ( 2)), the dummy variable model (Equation (3)), the linear mixed model (Equation ( 4)), and the Bayesian hierarchical model (Equation ( 5)).There were no significant differences among the three statistics of fit in the dummy variable model, linear mixed model, and Bayesian hierarchical model.However, when compared with the general biomass model, dummy variable model, linear mixed model, and Bayesian hierarchical were significantly better; F-test and residual plots also verified the significant differences.
Some studies only compared the dummy variable model with the mixed effect model in solving large-scale forestry biomass and growth estimation accuracy problems [34,35].In this study, an alternative method, the Bayesian hierarchical model, was used to model the large-scale tree biomass.Bayesian inference has been developed in forest biomass estimates [25,31].The main difference between Bayesian and classical approaches is how to define their prior knowledge of the sample data.The classical approach assumes the parameters are a fixed and unknown constant, whereas the Bayesian approach assumes that the parameters follow some statistical distribution, especially for a bivariate normal distribution [25,31].Due to known prior data, the Bayesian method has advantages in the prediction of small sample data.Zapata-Cuartas et al. [31] built the biomass model using the Bayesian and the classical approach to compare model precision of different sample sizes of 6, 10, 20, 30, 40, 50 and 60, and found that model efficiency (RMSE) of two approaches was almost similar when the sample size was larger than 60, otherwise the Bayesian approach was superior.Huang et al. [47] also indicated that the fitted precision for tree biomass using the Bayesian approach was better than the classical approach when the sample size was smaller than 50.Thus, using the Bayesian method, with small sample data to estimate forest biomass will significantly reduce the time and costs of sampling.The emphasis of this study is mainly on methodology.Choosing between the dummy model and mixed model for prediction of forest biomass has been a hot research topic with the aim of improving the accuracy of models for application in large regional scales [34,35].For practical applications, the choice of model depends on a number of categories (i.e., regions) and the sample number in each category.If the number of categories is small (less than 10), the dummy model is preferred; if the number of categories is large, the mixed model is more appropriate [35].However, the Bayesian hierarchical method has advantages in the prediction of small sample data (numbers < 60) [31,47].In this paper, we defined the categories by the six geographical regions, and the number of samples per region was equal to 60.Furthermore, the sample trees from each region came from various larch plantations which covered different site conditions, stand ages, and stand densities.Thus, it is necessary to analyze the random effects, taking the specific region's variables as random parameters in the mixed model.Because the sample number is equal to 60, the effect of the Bayesian hierarchical model is the same as the linear mixed model.Based on this knowledge, we tend to recommend the mixed model and the Bayesian hierarchical model for developing larch stem wood, branch, stem bark, needle, root, and total biomass models.
The three methods (dummy variable model, linear mixed model, and Bayesian hierarchical model) have no differences in regional biomass estimation ability.The modeling results showed that the total biomass of a tree with the same diameter gradually decreased from the southern to northern regions in China, except in Hebei Province (S3) (Figure 3).This may be because the water and heat conditions in the southeastern region are better and the trees have enough growing space, but from south to north, the water and heat conditions worsen, which impacts the growth and development of the trees.In addition, the sample plots in Hebei Province (S3) were located in a plateau region with high winds, sandy soil, and limited recipitation.The climate conditions may have resulted in larch having the lowest biomass in this region (S3).

Conclusions
Based on the biomass data for larch plantations in China, the generalized single tree biomass model suitable for regional scale forest biomass estimation was developed using the dummy variable model, linear mixed model, and Bayesian hierarchical model, which solved the problem of forest biomass estimate accuracy amongst different geographical scales.Our results indicate that using the mixed model, dummy variable model, and Bayesian hierarchical model improved the model-fitting results, providing an effective method for estimating larch biomass at the regional scale.Note that if the model fitting process accounts for differences in species and origins as a nested factor based on regional differences, the model fitting results should improve, and the mixed model, dummy variable model, and Bayesian hierarchical model would be more effective than the general allometric equations.The choice of the dummy variable model, mixed model, or Bayesian hierarchical model mainly depends on number of categories and samples.If the number of categories is larger, the mixed model and Bayesian hierarchical model should be the better choices.When the sample numbers are less than 60, the Bayesian hierarchical model may be more appropriate.In general, the mixed model and Bayesian hierarchical model are more flexible and applicable.The modeling results showed that the total biomass of a tree with the same diameter gradually decreased from southern to northern regions in China, except in Hebei Province (S3).

Figure 2 .
Figure 2. Residual boxplots of the biomass component of all models (A-D) represent the general biomass model, linear mixed effects model, the dummy variable model, and the Bayesian hierarchical model respectively).

Figure 2 .
Figure 2. Residual boxplots of the biomass component of all models (A-D) represent the general biomass model, linear mixed effects model, the dummy variable model, and the Bayesian hierarchical model respectively).

Figure 3 .
Figure 3.The regression curves of the linear mixed model for different regions.

Figure 3 .
Figure 3.The regression curves of the linear mixed model for different regions.

Table 1 .
The code and basic featuresof study regions.

Table 2 .
Biomass of each larch component (kg tree −1 ) in six regions in China.

Table 7 .
Evaluation statistics of the general model, dummy variable model, mixed effects model, and Bayesian hierarchical model for biomass modeling.

Table 7 .
Evaluation statistics of the general model, dummy variable model, mixed effects model, and Bayesian hierarchical model for biomass modeling.