Incorporating competition factors in a mixed-effect model with random effects of site quality for individual tree above-ground biomass growth of Pinus kesiya var. langbianensis Mingchuan Nong 1,#, Yan Leng 1,2#, Hui Xu1,#, Chao Li1, Guanglong Ou1* 1 Key Laboratory of State Forestry Administration on Biodiversity Conservation in Southwest China, Southwest Forestry University, Kunming, China; 2 School of Biotechnology and Engineering, Dianxi Science and Technology Normal University, Lincang, China *Corresponding author: olg2007621@126.com # These authors contributed equally to this work (Received for publication 7 August 2018; accepted in revised form 21 November 2019) Abstract Background: Accurate biomass estimation has critical effects on quantifying carbon stocks and sequestration rates, and above- ground biomass (AGB) growth models are a key component of tree biomass estimation. The study objective was to develop a growth model for AGB of an individual tree by combining competition factors and site quality using a mixed-effect model. Methods: The AGB of 128 sampling trees was investigated for Simao pine (Pinus kesiya var. langbianensis) at three typical sites near Pu’er City of Yunnan Province, China. Richards’ Equation was used for the basic growth model (BM) of the AGB, and a mixed-effect model with random effect of site quality (MEM) based on BM and a mixed-effect model with fixed effect of competition factors (MEMC) based on MEM were built using S-plus. Results: Both mixed-effect models are significantly better than the basic model in fitting and predicting the individual tree AGB growth for Simao pine, but the MEM is better than the MEMC. Moreover, the mixed-effect model with competition factors and site quality is the optimal estimation model due to its highest prediction precision (P=86.08%) as well as the lowest absolute average relative error (RMA=54.34%) and average relative error (EE =6.45%). Conclusion: A model including site quality and competition factors can be used to improve the tree AGB growth estimation for the individual tree AGB growth of Simao pine. New Zealand Journal of Forestry Science Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 https://doi.org/10.33494/nzjfs492019x27x E-ISSN: 1179-5395 published on-line: 11/12/2019 © The Author(s). 2019 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Research Article Open Access as the potential climate mitigation benefits of forest management and for assessing climate change impacts on forest ecosystems (Temesgen et al. 2015; Cosmo et al. 2016). Direct measurements of biomass are time consuming and expensive, so the usual procedure is to destructively sample a subset of trees for the development of allometric models that predict biomass from commonly measured tree characteristics like diameter (DBH) and height (H) (Houghton 2005; Saint- Introduction Forests are important components of terrestrial ecosystems, and they play vital roles in the global carbon balance (Woodwell et al. 1978; Talhelm et al. 2014; Sleeter et al. 2018). Trees are an important element of forest ecosystems’ carbon storage (Goodale et al. 2002; Houghton 2005; Houghton et al. 2009; Luo et al. 2014). Accurate estimates of tree biomass are critical for quantifying carbon stocks and fluxes as well Keywords: Above-ground biomass, competition factors, growth model, mixed-effect model, Pinus kesiya var. langbianensis, Richards’ Equation, site quality Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 2 Andre et al. 2005; Somogyi et al. 2007; Sileshi 2014; Chave et al. 2014; Temesgen et al. 2015). Equations for tree biomass estimation can be categorized into regional biomass conversion factors, stand-level and tree-level biomass Equations (Di Cosmo et al. 2016; Jagodziński et al. 2018). Tree-level biomass models are used to estimate the total and the components of individual trees with easily measurable tree inventory attributes (e.g. DBH, H). The tree-level biomass model became a research priority because of its better performance and higher predicting precision compared with regional and stand-level biomass models (Temesgen et al. 2015). More than 2600 biomass models related to more than 100 tree species have been constructed all over the world (e.g. Ter-Mikaelian & Korzukhin 1997; Chojnacky 2003; Jenkins et al., 2003, 2004; Zianis 2005; Muukkonen 2007). Although many studies on tree-level biomass models have already been developed, the existing Equations are generally not good enough and need improvement (Temesgen et al. 2015). Moreover, almost all of these studies were static models and focused on the relation among the individual tree biomass and predictors, whereas few studies have considered the growth of individual tree biomass using growth Equations. Furthermore, spatial and temporal autocorrelation between individuals and groups may be neglected in modelling tree biomass growth (Li & Zhang 2010). It is necessary to develop comprehensive biomass estimation models with better prediction performance that consider differences in stand density and among sites (Temesgen et al. 2015). Forest biomass estimation variation can be observed between different ecological zones and sites (Henry et al. 2011, De-Miguel et al. 2014). Trees have different growth rates due to differences between sites, and models perform better when site quality is considered in the model (Vanninen et al. 1996; Rohner et al. 2013). Furthermore, site quality also has an important effect on tree growth, and it can be quantified as a site index (SI) which represents the potential productivity of the forestland by using the dominant height of stands at a particular standard age (Running & Mcleod 1988; Carmean et al. 1989; Sturtevant & Seagle 2004; Waring et al. 2006). Stand density can be used to reflect the effects of crowding and competition among trees on tree growth (Curtis 1985; Zeide 2005). Competition results from interactions between many biotic and abiotic factors, and this affects the forest structure (Sahney et al. 2010). With increasing stand density, individual trees can be restricted in their growth and even die (Bragg 2001). Many studies have shown that tree size is strongly affected by competition (Coomes & Allen 2007; Coates et al. 2004, 2009). By incorporating the size and distance of neighboring trees into the model, predictions of growth and mortality can be improved (Mctague & Weiskittel 2016). Biomass growth Equations are generally constructed by using forestry investigation data and contain multiple repeated observations of individual trees. The characteristics of these data led to spatial and temporal correlations among observations in the same sampling unit (Lhotka & Loewenstein 2011; Timilsina & Staudhammer 2013). A previous parameter estimation of the individual tree growth Equation was estimated by ordinary least squares (OLS). This would violate the regression assumptions of homoscedasticity of variances and independence of residuals and lead to inaccurate estimates if the growth model was developed based on non-independent data by OLS (Budhathoki et al. 2010; Njana et al. 2016). To address these problems, many researchers have tried to develop new growth models by incorporating mixed- effect modelling techniques (Budhathoki et al. 2010; Lhotka & Loewenstein 2011). A mixed-effect model, which consists of fixed effects and random effects, provides a flexible and powerful tool for the analysis of grouped data (Pinheiro & Bates 2000). Fixed effects can indicate the average relationships of a dependent variable with independent variable(s), and random effects can reflect the difference among groups (Razali et al. 2015). The advantage of mixed-effect models is that they can fit growth and yield data in forestry fairly well via multilevel random effects (Gregoire et al. 1995), and the prediction accuracy of such models can be improved through modifications of random effects (Calama & Montero 2004). In recent years, mixed-effect models have been widely applied in forestry due to their better fitting performance and prediction precision (Ramon et al. 2006), and they have been used to estimate the dominant height growth (Fang & Bailey 2001), diameter increment (Lhotka & Loewenstein 2011), tree stand basal area (Gregoire et al. 1995), basal area growth of individual trees (Budhathoki et al. 2008), wood density (Li & Jiang 2013), stand volume (Li 2012), and biomass (Zhang & Borders 2004; Fehrmann et al. 2008; Fu et al. 2012; Njana et al. 2016) by incorporating different random effects (e.g. tree-level, plot-level). Therefore, it is crucial to accurately predict biomass growth by constructing a comprehensive model with good prediction performance using a mixed-effect model that considers competition factors and site quality. Simao pine (Pinus kesiya var. langbianensis), a geographic variation of P. kesiya, is distributed in mountain areas from the northern tropical zone to the southern subtropical zone of Yunnan Province, China. Its distribution area and stocking volume account for 11% of the forestland of Yunnan Province (Compilation Committee of Yunnan Forest 1986). It has been an important tree species for afforestation in Yunnan due to its rapid growth (Southwest Forestry College & Forestry Department of Yunnan Province 1988). Moreover, Simao pine forests provide a range of goods and services and have important economic and ecological values (Wu & Dang 1992; Wen et al. 2010; Yue & Yang 2011; Li 2011). Therefore, it is important to be able to estimate biomass growth comprehensively and with high prediction accuracy by incorporating site index and competition factors to assess the potential value of these forests and to guide forest management. In this study, natural forests of Simao pine were studied by sampling the above-ground biomass (AGB) of 128 trees at three sites. Mixed-effect models incorporating site quality and a competition factor were used to construct the AGB growth model based on a transformation of Richards’ Equation. This study aimed to: (1) explore a comprehensive biomass growth model with high prediction accuracy for estimating the AGB growth of Simao pine; and (2) explain the impacts on improving AGB estimation from site quality and competition factors. Methods Study Region The study region is located in Mojiang county, Simao district and Lancang county which belong to Pu’er City, southwest of Yunnan Province, P. R. China. In this city, mountain areas comprise 98.3% of the overall region, and the study region is located between 22°02′N to 24°50′N and 99°09′E to 102°19′E. Three typical geographic areas with different climates, Tongguan Town of Mojiang County (Site I), Yunxian Town of Simao District (Site II) and Nuofu Town of Lancang County (Site III), were selected as study sites (Fig. 1), and the elevation of the sites are from 1400 m to 1600 m above sea level. The mean annual temperatures and annual precipitation both decrease from south to north (Fig. 2). Lancang county has the highest annual mean temperature (19.7 ℃) and annual precipitation (1586.5 mm), and Mojiang county has the lowest values (the annual mean temperature is 18.4 ℃ and the annual precipitation is 1306 mm) (Ou et al. 2016). Data Collection and Processing Data investigation A total sample of 128 pines in 45 plots with an area of 900 m2 were selected and investigated in the study areas (Fig. 1). Plot information including latitude, longitude, degree of slope, and aspect of slope was recorded. Tree age (A), diameter at breast height (DBH) and tree height (H) of the 128 sampled trees were recorded (Fig. 3). Tree age was determined by counting the number of annual growth rings of the stump of each felled sample trees. Then, we also recorded the basic characteristics of the surrounding trees within 5 metres of the sample trees, including tree species, DBH, H and the distance to calculate the competition index (CI). We calculated the average height of dominant trees (Ht) and the average age of stands (A) to calculate the site index (SI). The average height of the dominant trees for each plot is the mean of the three highest trees, and the average age of the stand of each plot is the mean age of the standard trees with a DBH similar to the average DBH of the plot. Biomass measurement According to the sample collection method for tree biomass modelling in China, the biomass of each tree component was measured one at a time (Zeng et al. 2011). The biomass of the stem was measured by the method of volume density. Firstly, felled trees were segmented and weighed, and the fresh weight was measured. Secondly, the segment length and diameter were measured, and the volume of the trunk was calculated. Branches and leaves were measured by the method of the graded sample branch. Dead branches and fruits were measured by the method of total weight. Thirdly, the samples from the Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 3 FIGURE 1: Location of the study sites. FIGURE 2: Monthly mean temperatures and monthly gross precipitation of three typical sites. The data are average values from 1980 to 2010 measured by the respective county weather stations. The lines are the monthly mean temperatures, and the bars are the monthly gross precipitation. different components were dried to a constant weight at 105 ℃ using an oven, and the sample density of the wood and bark was also measured using the drainage method. Finally, the biomass of wood and bark of the sampled trees was calculated using the volume and the corresponding sample density, and the branch biomass and needle biomass of sample trees were calculated Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 4 using fresh weight and the corresponding dry matter ratio. The sample trees were divided into two sets by random selection; one set with 96 sample trees was used to fit the models, and the other one was used for the test. The basic characteristics are listed in Table 1. FIGURE 3: Scatter of the diameters at breast height (DBH) and tree height (H) vs tree age. a: DBH vs age, b: H vs age. Data set SI Class (m) N Age (years) DBH (cm) Height (m) AGB (kg) Mean Standard deviation Mean Standard deviation Mean Standard deviation Mean Standard deviation Fitting 12 7 28.43 2.59 21.24 2.83 16.84 1.82 197.50 50.37 14 34 38.12 2.65 22.49 1.81 16.86 0.79 297.96 63.59 16 19 38.58 4.44 20.71 3.00 14.86 1.42 294.57 111.64 18 26 38.73 2.20 34.00 2.24 21.45 1.17 622.36 83.73 20 10 36.40 4.54 31.59 4.71 23.26 3.17 684.21 197.22 Total 96 37.49 1.50 26.11 1.30 18.37 0.67 418.06 46.68 Testing 12 2 32.00 2.00 28.75 3.55 19.30 0.10 345.37 71.94 14 5 39.00 5.67 19.66 3.20 16.78 1.52 185.89 60.91 16 13 50.00 5.59 27.68 3.55 18.71 1.57 395.02 92.28 18 10 44.70 3.68 37.88 2.57 23.80 2.26 704.72 96.04 20 2 48.00 2.00 40.50 5.60 32.25 2.75 830.54 283.31 Total 32 45.38 2.77 30.48 2.08 20.88 1.19 483.24 61.21 All 12 9 29.22 2.07 22.91 2.50 17.38 1.44 230.36 45.79 14 39 38.23 2.40 22.12 1.62 16.85 0.71 283.59 56.11 16 32 43.22 3.57 23.54 2.34 16.42 1.09 335.38 75.59 18 36 40.39 1.92 35.08 1.77 22.11 1.05 645.24 65.68 20 12 38.33 3.98 33.08 4.08 24.76 2.83 708.60 167.35 Total 128 39.46 1.35 27.20 1.11 19.00 0.59 434.35 38.18 TABLE 1. Basic characteristics of the sampled trees Competition index calculation To calculate the competition index (CI), we used the formula of Hegyi (1974), Equation 1. (1) where CIi is a competition index, Di is the diameter of sample tree i, Dj is the diameter of competition tree j around the sample tree, and DISTij is the distance from sample tree i to competition tree j. Site index calculation To derive the site index (SI), the dominant height and age of each plot were measured. SI was calculated for Simao pine using the Equation by Wang (2003), Equation 2. (2) where SI is the site index, Ht is the average height of dominant trees of each plot, A is the average age of each plot, and the base age is 20 years according to Wang (2003). The site index of the plots includes five types from 12 m to 20 m according to an interval of 2 m in this study (Table 1). Model fitting Basic model (BM) Richards’ growth Equation developed from Bertalanffy’s growth theory is used to describe biological growth changes over time. Growth Equations such as Monomolecular, Gompertz, and Logistic Equations are Richards’ growth Equations of the special form. Richards’ Equation is widespread in forestry because of its flexibility and excellent fitting performance (Richards 1959, Liu & Li 2003, Rohner et al. 2013). In the present study, Richards’ growth Equation is used as a basic biomass growth model. The general expression of Richards’ growth Equation is listed in Equation 3. (3) where parameters A, b and c are given in Equations 4–6: (4) (5) (6) where is the response variable describing the change in biomass with tree age (t), A is the asymptote of the maximum parameter for tree growth, b is the growth rate parameter that indicates the rate a tree approaches its asymptotic biomass, c is the parameter related to m, m is the power exponent of anabolism, η is the anabolism constant, and β is the catabolism constant. Parameter A in Equation 3 is the most unstable parameter; thus, a transformation of the Equation is constructed to solve this problem by using the parameter a, which is an expected-value parameter when t = t0 to replace parameter A (Fang & Bailey 2001). Therefore, we finally selected the transformation of Richards’ Equation to construct a basic biomass growth model. The Equation is shown in Equation 7. (7) where y is the organ biomass of an individual tree, a is the progressive parameter of organ biomass growth, b is the rate of growth, c is the curve shape parameter, t is the tree age, and t0 is a fixed reference age that may be fixed at any positive value according to the specific research situation (Fang & Bailey 2001). Its value was set at 20 years in this study. Mixed-effect model without fixed effects from competition factors (MEM) Forestry growth and yield data are affected by the sampling region (e.g. different regions may have different site conditions) and the correlation among the trees in the same sampling position. The data from within a sampling unit is dependent; thus, autocorrelation and heterogeneity are common in these data (Gregorie 1987). The mixed-effect modelling technique can partly remove the negative impact of heterogeneity and autocorrelation within-plots by using variance (e.g. power, exponential function) and covariance (e.g. time- autocorrelation function) structure. This technique can also explain the plot parameter variability by selecting appropriate covariates (Fang & Bailey 2001). In our study, site quality was taken as a random effect to select the mixed parameters, and power and exponential functions were used for describing the variance structures. Three time-autocorrelation functions, including the autoregressive time correlation structure with order 1 (AR(1)), continuous time AR(1) structure (CAR(1)) and autoregressive moving-average structure with both order 1 (ARMA(1,1)) function, were used to describe covariance structures. The mixed parameter selection to determinate the model forms, variance and covariance structure were all given by the research of Pinheiro & Bates (2000). Mixed-effect model with fixed effects from competition factors (MEMC) Based on the MEM, competition factors were considered fixed effects to input the Equation parameters of the MEM, and both CI and the quadratic effect (CI2) are incorporated into the models. Model evaluation In this study, Log likelihood (logLik), Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) were used to evaluate the model fitting results (Equations 8–10). Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 5 (8) (9) (10) where is the maximum likelihood estimation of θ for the Likelihood function of model , x is the random sample, q is the number of unknown parameters, and n is the sample capacity. Moreover, the sum relative error (RS), mean relative error (EE), absolute mean relative error (RMA) and prediction precision (P) were used to test the model prediction performance (Equations 11–14). (11) (12) (13) (14) where yi is the measured value, is the estimated value, is the mean value of the estimated value, ta is the distribution value of t (when confidence level is a = 0.05), N is the sample capacity, and T is the parameter number of the regression curve Equation. Results Mixed-effect model without fixed effects from competition factors (MEM) Mixed parameter selection The fitting results of different parameter combinations are listed in Table 2. Models with mixed parameters have lower AIC and BIC and higher logLik values than the basic model, referred to as the transformation of Richards’ Equation (BM). The optimal fitting result emerges when the parameter a is regarded as the mixed parameter (AIC=1305.367, BIC=1318.189, logLik= -647.684). The mixed-effect model is listed in Equation 15. (15) where y is the above-ground biomass, a is the progressive parameter of organ biomass growth, ua is the parameter for the random effect from the site index, b is the growth rate of AGB, t is the tree age, t0 is the standard age with a value of 20 years, and c is the shape parameter. Variance and covariance structure selection A comparison of the variance and covariance structure fitting results of the MEM are listed in Table 3. The optimal fitting result emerged when the power function is regarded as a variance structure but none are regarded as covariance structures (AIC=1224.4, BIC=1239.8, logLik=-606.2). Mixed-effect model with fixed effects from competition factors (MEMC) Basic model construction Based on the MEM, when parameter a is regarded as a random effect and the basic individual tree competition index is regarded as a fixed effect, the optimal model form is listed in Equation 16. Comparisons of different effects in the AGB growth mixed-effect model are listed in Table 4. The optimal fitting results emerge in the mixed-effect model with competition factors and site index (MEMC) (AIC=1298.1, logLik=-640.1). The fitting results of parameters are listed in Table 5. Parameter b is extremely significant at the p=0.01 level. (16) where y is the above-ground biomass, CI is the competition index, CI2 is the quadratic effect of CI, a is the progressive parameter of organ biomass growth, ua is the parameter for random effects from the site index, b is the growth rate of AGB, b1 and b2 are the estimated parameters of the fixed effect from and to parameter b, respectively, t is the tree age, t0 is the standard age with a value of 20 years, c is the shape parameter, and c1 and c2 are the estimated parameters of the fixed effect from CI and CI2 to parameter c, respectively. Variance and covariance structure selection The fitting results of the variance and covariance structures are listed in Table 6. The optimal fitting result emerged when the power function is regarded as a variance structure but none are regarded as covariance structures (AIC=1230.2, BIC=1250.7, logLik=-607.1). The optimal fitting results of the mixed-effect model taking the competition factor as a fixed effect are listed in Table 7, and the optimal model form is shown in Equation 17. (17) Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 6 Mixed para- meter logLik AIC BIC LRT p-value No -667.0 1342.0 1352.3 - - a -647.7 1305.4 1318.2 38.7 <0.0001 b -652.6 1315.1 1327.9 28.9 <0.0001 c No convergence a, b -647.7 1307.4 1322.8 38.7 <0.0001 a, c -647.7 1307.4 1322.8 38.7 <0.0001 b, c No convergence a, b, c -647.7 1309.4 1327.3 38.7 <0.0001 TABLE 2. Mixed parameter selection of the mixed-effects model for AGB growth where y is the above-ground biomass, a is the progressive parameter of organ biomass growth, ua is the parameter for random effects, is the growth rate of AGB, t is the tree age, t0 is the standard age with a value of 20 years, c is the shape parameter, c1 is a mixed parameter of parameter c, c2 is a mixed parameter of parameter c, CI is the basic competition index, and CI2 is the quadratic effect of CI. Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 7 Model evaluation The final fitting results of BM, MEM and MEMC are listed in Table 8. The parameter a did not differ significantly between the three models due to overlapping intervals of the estimated value, although BM has the lowest value with 61.799, and the parameters b and c have significant differences among the three models. MEM has the No. Variance structure Covariance structure logLik AIC BIC LRT p-value 1 No No -647.7 1305.4 1318.2 - - 2 Power No -606.3 1224.4 1239.8 92.839 <0.0001 3 Exponential No -617.8 1247.5 1262.9 59.889 <0.0001 4 No AR(1) -647.7 1307.3 1322.7 0.099 0.7532 5 No CAR(1) -647.7 1307.7 1322.8 0.015 0.9015 6 No ARMA(1,1) -647.7 1305.7 1322.7 0.095 0.9777 TABLE 3. Mixed-effects models considering variance and covariance structures for AGB growth TABLE 4. Comparison of different effects on the AGB growth mixed-effects model Model logLik AIC BIC LRT p-value Mixed-effects model with competition factors and regional effect -640.1 1298.1 1321.2 - - Mixed-effects model with regional effect -647.7 1305.4 1318.2 25.0 0.0001 TABLE 5. Estimated parameters of the AGB growth mixed-effects model incorporating competition factors as fixed effects Parameter Estimated value Standard deviation df t-Value p-value a 62.4101 34.2398 87 1.823 0.0718 b 0.0930 0.0225 87 4.128 0.0001 b1 -0.0015 0.0008 87 -1.812 0.0734 b2 0.00003 0.00001 87 1.956 0.0537 c 17.2559 10.4085 87 1.658 0.1009 c1 -0.5845 0.3855 87 -1.516 0.1331 c2 0.0088 0.0061 87 1.443 0.1526 TABLE 6. Mixed-effects models incorporating competition factors as fixed effects considering variance and covariance structures for AGB growth (* the model with no significant parameters b1 and b2.) No. Variance structure Covariance structure logLik AIC BIC LRT p-value 1 No No -640.1 1298.1 1321.2 2 Power No -608.1 1236.2 1261.9 63.9 <0.0001 3 Exponential No -616.2 1252.5 1278.1 47.7 <0.0001 4 No AR(1) -639.4 1298.8 1324.4 1.4 0.2394 5 No CAR(1) -639.4 1298.8 1324.4 1.4 0.2394 6 No ARMA(1,1) -639.6 1299.2 1324.8 1.0 0.3209 7* Power No -607.1 1230.2 1250.7 66.0 <0.0001 optimal fitting performance because of the lower values of AIC and BIC and the highest value of logLik (Table 8), but the difference between the MEM and MEMC is not significant because of the higher p value (p value = 0.113) according to the likelihood ratio test (LRT). While MEMC has the minimum value of EE (6.45%) and RMA (54.34%), and the highest value of prediction precision TABLE 2: Confusion matrix P (86.08%), BM has the minimum value of RS (7.71%) (Table 8). Moreover, the heteroscedasticity of the residual is not found in either mixed-effect models, but it is obvious in the BM, and the MEMC has the narrowest interval of standardised residual (Fig. 4). Therefore, the MEMC is the optimal model for the individual tree AGB growth of Simao pine. Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 8 TABLE 7. Fitting results of the three AGB growth models. Estimation parameter BM MEM MEMC Estimated value p-Value Estimated value p-Value Estimated value p-value a 61.7990±14.5400 <0.0001 76.3382±19.2527 0.0001 76.9496±11.6512 <0.0001 b 0.0516±0.0200 0.0346 0.0111±0.0115 0.3389 0.0656±0.0152 <0.0001 c 6.3479±1.0760 0.0013 2.9552±0.6778 <0.0001 9.1919±2.9250 0.0023 c1 - - - - -0.1017±0.0341 0.0037 c2 - - - - 0.0003±0.0001 0.0033 logLik -667.0 -606.2 -607.1 AIC 1342.0 1224.4 1230.2 BIC 1352.3 1239.8 1250.7 D-matrix - D=[28.1220] D=[2.9887×10-9] Heteroscedastic function value - Residual error 1.4507 1.7662 2.3904 TABLE 8. Comparison of the three AGB growth models. Model Fitting index Testing index logLik AIC BIC LRT p-value RS(%) EE (%) RMA (%) P (%) BM -667.0 1342.0 1352.3 7.71 9.82 61.93 85.81 MEM -606.2 1224.4 1239.8 121.6 <0.0001 -62.08 -37.46 55.16 79.85 MEMC -607.1 1230.2 1250.7 119.9 <0.0001 -34.26 6.45 54.34 86.08 FIGURE 4: Scatter of standardised residual vs predicted AGB for three models. a: basic model (BM); b: Mixed effect model without fixed effects from competition factors (MEM); and c: Mixed effect model with fixed effects from competition factors (MEMC). Discussion The different combinations of mixed parameters were selected, and the mixed model with mixed parameter a was considered as the optimal basic model because of the low AIC (Table 2). Parameter a represents the estimated asymptotic biomass at a standard age; however, trees have different growth potentials at the same standard age depending on the site quality. Many factors would influence tree maximum diameter at an equal basal age, such as water-holding capacity, elevation, and slope (Rohner et al. 2013). Thus, selecting parameter a as a mixed parameter in the model can improve the estimation and indicates that the site quality has a significant effect on the age-biomass relationship. Site quality had important effects on tree growth (Robichaud & Methven 1993). Past studies generally used site characteristics (a plot-level variable) as variables for biomass growth models and showed that site characteristics had a significant effect on tree growth (Lee et al. 2004). Westfall (2016) reported that mixed models, including plot random effects, can reduce prediction bias and variance for populations to a great extent compared to fixed effect models; Lhotka & Loewenstein (2011) reported analogous results. Similarly, our study found that mixed-effect models, including the random effects of site quality, are better than basic models (i.e. higher logLik and lower AIC and BIC), and incorporating random effects can improve the biomass growth model for Simao pine. Our results are consistent with previous studies for both static models (Subedi & Sharma 2011; Rohner et al. 2013; Ou et al. 2016; Chen et al. 2017; Huff et al. 2018), and growth models (Budhathoki et al. 2008; Timilsina & Staudhammer 2013; Westfall 2016). Moreover, the competition factor is an important predictor variable for the individual tree model because it intensively affects tree growth (Lhotka & Loewenstein 2011). In our study, MEMC had a better estimation due to incorporating the competition factor as a fixed effect compared with MEM, and the predictive ability was significantly improved (i.e. largest P and smallest EE and RMA). Specifically, the predictive accuracy increased to 86%. Thus, it can obviously improve the mixed model predictive accuracy if competition is regarded as an independent variable, and it also indicates that tree competition is the critical element for predicting individual tree biomass growth. Therefore, mixed-effect models have been used in forestry because of their superior fitting and prediction accuracy compared with traditional models (Huff et al. 2018). The mixed-effect model, only including the random effects of site quality, can greatly improve the fit performance of the model, and the mixed-effect model incorporating competition factors can further improve the prediction ability. In addition, autocorrelation among measurement data may result in biased estimates of the model parameter for biological data (Budhathoki et al. 2010). The mixed-effect model was then introduced to address this challenge by defining the variance and covariance structures of random effects in parameter estimations (West et al. 2007; Smith et al. 2014; De- Miguel et al. 2014; Njana et al. 2016). We also found that our prediction model had issues with heteroscedasticity which was corrected by using the power variance function (the optimal variance function) in the estimation process (Budhathoki et al. 2008). Therefore, the mixed-effect model had a good fitting performance which is consistent with previous studies (Subedi & Sharma 2011; De-Miguel et al. 2014). In contrast, all of the time-autocorrelation covariance structures in our study did not improve the model performance which indicates that the estimation models considering time autocorrelation cannot improve fitting. This might be attributable to the biomass sample data without a time series (Eisfelder et al. 2017). However, it is unrealistic and impossible to obtain the time series data of the biomass because destructive sampling was performed to collect biomass data (Temesgen et al. 2015). Therefore, we investigated sample trees with different ages (from 8 years to 80 years) in different locations to replace the time series data. The AGB change rule along with the tree ages can be explained by using the mixed-effect model considering site quality and competition factors. Thus, the investigated methods for estimating the AGB growth would be reasonably practical. Additionally, we did not consider spatial autocorrelation among trees in the same plot because Simao pine is an intolerant tree species and because sampling trees of different ages occurs at distant locations. Furthermore, the data characteristics of the independent variables are important for selecting the form of mixed effects. A random effect is applicable if the variable is a grouped one, but a fixed effect is appropriate if the data are a continuous variable in a mixed-effect model (Pinheiro & Bates 2000). In this study, the site index and competition index were incorporated into the mixed-effect models as random and fixed effects, respectively. The site index (SI) is often used as grouped data in forestry. SI tables of the dominant tree species are established to predict potential productivity of forestland; SI is based on the dominant height classes with a 2-m interval at a standard stand age (Meng 2006; Duan et al. 2009 ). Thus, it was considered as a random effect into the mixed-effect model in this study. The competition index is a continuous variable, and it was calculated using the distances between each sample tree and its neighbouring competing trees and their DBH in this study (Hegyi 1974). Thus it was incorporated into the mixed-effect model as a fixed effect. In addition, the national continuous forest inventory of China (NCFI) is being conducted using permanent plots and carried out every five years to reflect the dynamic change of forest resources at the national scale since 1987 (Kang 2011). Fang et al. (2001, 2002) calculated forest biomass carbon using the NCFI database by the variable biomass expansion factor method, with estimation error of less than 3% at the regional and national scale. At the stand level, the location of each tree in each permanent plot and the height of the stands had been recorded (Kang 2011). Therefore, the CI of each tree can be calculated and the SI of the plots can be found by looking up the SI tables of the different dominant tree species. So the comprehensive estimation method has the potential to improve biomass growth estimation and Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 9 reflect the dynamical change of stand biomass. Conclusions To improve estimation of individual tree AGB growth, a nonlinear mixed-effect model was developed by incorporating random effects of site quality and fixed effects of competition factors based on a transformation of Richards’ function. We found that the mixed-effect model was significantly better than the BM because of its better fitting and prediction performance, and the MEMC is the optimal estimation model due to its highest prediction precision. Therefore, comprehensive biomass growth estimation considering the site index and competition index can be used to predict the AGB growth of Simao pine trees, and it is a potential method for AGB growth estimation of other tree species. List of abbreviations AGB: Above-ground biomass; BM: Basic model; MEM: Mixed-effect model without fixed effect from competition factors; MEMC: Mixed-effect model with fixed effect from competition factors; DBH: Diameter at breast height; H: tree height; CI: competition index; SI: site index. Competing interests The authors declare that they have no competing interests. Authors’ contributions GO conceived the study, analyzed the data, and wrote the manuscript. MN and LY participated in field work and analyzed the data, CL participated in field work, and HX reviewed the manuscript. All authors read and approved the final manuscript. Acknowledgements We acknowledge the support from the Forestry Bureau of Mojiang County, Simao District and Lancang County for assistance with field work, and the Key Laboratory of State Forestry Administration on Biodiversity Conservation in Southwest China (Southwest Forestry University) for assistance with measurements of plant samples at the laboratory. We thank Enliang Li, Zhigang Liang and Junfeng Wang for assistance with the field work. Wenjun Wu and Mingquan Huang edited the manuscript, and Guangxing Wang (USA) reviewed and provided valuable comments. Additionally, thanks to the anonymous external reviewers for their valuable comments. Funding This study was supported by funding from the National Natural Science Foundation of China (31560209, 31760206, 31660202), and the Scientific Research Foundation of Southwest Forestry University (111416), and the Ten-Thousand Talents Program of Yunnan Province, China (YNWR-QNBJ-2018-184). Ethics approval and consent to participate Not applicable Consent for publication Not applicable Availability of data and materials Please contact author for data requests References Bragg, D.C. (2001). A local basal area adjustment for crown width prediction. Northern Journal of Applied Forestry, 18(1), 22-28. Budhathoki, C.B., Lynch, T.B., & Guldin, J.M. (2008). Nonlinear mixed modeling of basal area growth for shortleaf pine. Forest Ecology and Management, 255(8-9), 3440-3446. https://doi.org/10.1016/j. foreco.2008.02.035 Budhathoki, C.B., Lynch, T.B., & Guldin, J.M. (2010). Development of a shortleaf pine individual-tree growth Equation using non-linear mixed modeling techniques. In: Stanturf, John A. (Ed.), Proceedings of the 14th biennial southern silviculture research conference; 2007 February 26-March 1; Athens, GA. [General Technical Report SRS-121]. Asheville, NC, USA: Southern Research Station, USDA Forest Service, pp. 519-520. Calama, R., & Montero, G. (2004). Interregional nonlinear height–diameter model with random coefficients for stone pine in Spain. Canadian Journal of Forest Research, 34(1), 150-163. https://doi. org/10.1139/X03-199 Carmean, W.H., Hahn, J.T., & Jacobs, R.D. (1989). Site index curves for forest tree species in the eastern United States. [General Technical Report NC-128]. St. Paul, MI, USA: North Central Forest Experiment Station, USDA Forest Service, 142 p. Chave, J., Réjou-Méchain, M., Búrquez, A., Chidumayo, A., Colgan, M.S., Delitti, W.B.C., Duque, A., Eid, T., Fearnside, P.M., Goodman, R.C., Henry, M., Martinez-Yrizar, A., Mugasha, W.A., Muller-Landau, H.C., Mencuccini, M., Nelson, B.W., Ngomanda, A., Nogueira, E.M., Ortiz-Malavassi, E., Pelissier, R., Ploton, P., Ryan, C.M., Saldarriaga, J.G., & Vieilledent, G. (2014). Improved allometric models to estimate the aboveground biomass of tropical trees. Global Change Biology, 20(10), 3177-3190. https://doi. org/10.1111/gcb.12629 Chen, D., Huang, X., Zhang, S., & Sun, X. (2017). Biomass modeling of larch (Larix spp.) plantations in China based on the mixed model, dummy variable model, and Bayesian hierarchical model. Forests, 8(8), 268. https://doi.org/10.3390/f8080268 Chojnacky, D.C. (2003). Allometric scaling theory applied to FIA biomass estimation. In: McRoberts, Ronald E.; Reams, Gregory A.; Van Deusen, Paul C.; Moser, John W. (Eds.), Proceedings of the third annual forest inventory and analysis symposium; 2001 October 17-19; Traverse City, MI. [General Technical Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 10 Report NC-230]. St. Paul, MN: USDA Forest Service, North Central Research Station. 208 p. https://doi. org/10.2737/NC-GTR-230 Coates, K.D., Canham, C.D., & Lepage, P.T. (2004). A neighborhood analysis of canopy tree competition: effects of shading versus crowding. Canadian Journal of Forest Research, 34(4), 778-787. https:// doi.org/10.1139/x03-232 Coates, K.D., Canham, C.D., & Lepage, P.T. (2009). Above- versus below-ground competitive effects and responses of a guild of temperate tree species. Journal of Ecology, 97(1), 118-130. https://doi. org/10.1111/j.1365-2745.2008.01458.x Compilation Committee of Yunnan Forest. (1986). Yunnan Forest. Beijing: China Forestry Publishing House. Coomes, D.A., & Allen, R.B. (2007). Effects of size, competition and altitude on tree growth. Journal of Ecology, 95(5), 1084-1097. https://doi. org/10.1111/j.1365-2745.2007.01280.x Cosmo, L.D., Gasparini, P., & Tabacchi, G. (2016). A national-scale, stand-level model to predict total above-ground tree biomass from growing stock volume. Forest Ecology and Management, 361, 269-276. https://doi.org/10.1016/j. foreco.2015.11.008 Curtis, R.O. (1985). Stand density measures: an interpretation. Forest Science, 16, 403-414. De-Miguel, S., Pukkala, T., Assaf, N., & Shater, Z. (2014). Intra-specific differences in allometric Equations for aboveground biomass of eastern Mediterranean Pinus brutia. Annals of Forest Science, 71(1), 101- 112. https://doi.org/10.1007/s13595-013-0334- 4 Di Cosmo, L., Gasparini, P., & Tabacchi, G. (2016). A national-scale, stand-level model to predict total above-ground tree biomass from growing stock volume. Forest Ecology and Management, 361, 269-276. https://doi.org/10.1016/j. foreco.2015.11.008 Duan, J., Ma, L., Jia, L., Hou, Y., & Gong, N. (2009). Establishment and application of site index table for Pinus tabulaeformis plantation in the low elevation area of beijing. Scientia Silvae Sinicae, 45(3), 7-12. Eisfelder, C., Klein, I., Bekkuliyeva, A., Kuenzer, C., Buchroithner, M.F., & Dech, S. (2017). Above- ground biomass estimation based on NPP time- series? A novel approach for biomass estimation in semi-arid Kazakhstan. Ecological Indicators, 72, 13-22. Fang, J., Chen, A., Peng, C., Zhao, S., & Ci, L. (2001). Changes in forest biomass carbon storage in China between 1949 and 1998. Science, 292, 2320-2322. Fang, J., Chen, A., Zhao, S., & Ci, L. (2002). Estimating biomass carbon of China’s forests: supplementary notes on report published in science (291: 2320- 2322) by Fang et al. (2001). Acta Phytoecologica Sinica, 26(2), 243-249. Fang, Z., & Bailey, R. L. (2001). Nonlinear mixed effects modeling for slash pine dominant height growth following intensive silvicultural treatments. Forest Science, 47(3), 287-300. Fehrmann, L., Lehtonen, A., Kleinn, C., & Tomppo, E. (2008). Comparison of linear and mixed-effect regression models and a k-nearest neighbour approach for estimation of single-tree biomass. Canadian Journal of Forest Research, 38(1), 1-9. https://doi.org/10.1139/X07-119 Fu, L.Y., Zeng, W.S., Tang, S.Z., Sharma, R.P., & Li, H.K. (2012). Using linear mixed model and dummy variable model approaches to construct compatible single-tree biomass Equations at different scales - A case study for Masson pine in Southern China. Journal of Forest Science, 58(3), 101-115. Goodale, C.L., Apps, M.J., Birdsey, R.A., Field, C.B., Heath, L.S., Houghton, R.A., Jenkins, J.C., Kohlmaier, G.H., Kurz, W., Liu, S., Nabuurs, G.J., Nilsson, S., & Shvidenko, A.Z. (2002). Forest carbon sinks in the Northern Hemisphere. Ecological Applications, 12(3), 891–899. Gregorie, T.G. (1987). Generalized error structure for forestry yield models. Forest Science, 33(2), 423- 444. Gregoire, T.G., Schabenberger, O., & Barrett, J.P. (1995). Linear modelling of irregularly spaced, unbalanced, longitudinal data from permanent- plot measurements. Canadian Journal of Forest Research, 25(1), 137-156. Hegyi, F. (1974). A simulation model for managing Jack- pine stands. In: Fries, J. (Ed.), Growth models for tree and stand simulation. Stockholm, Sweden: Royal College of Forestry, pp. 74-85. Henry, M., Picard, N., Trotta, C., Raphaël J.M., Valentini, R., Bernoux, M., & Saint-André, L. (2011). Estimating tree biomass of sub-saharan African forests: a review of available allometric Equations. Silva Fennica, 45(3B), 477-569. Houghton, R.A. (2005). Aboveground forest biomass and the global carbon balance. Global Change Biology, 11(6): 945–958. https://doi.org/10.1111/j.1365- 2486.2005.00955.x Houghton, R.A., Hall, F., & Goetz, S.J. (2009). Importance of biomass in the global carbon cycle. Journal of Geophysical Research, 116(G2). https://doi. org/10.1029/2009JG000935 Huff, S., Poudel, K.P., Ritchie, M., & Temesgen, H. (2018). Quantifying aboveground biomass for common shrubs in northeastern California using nonlinear mixed effect models. Forest Ecology and Management, 424, 154-163. https://doi. org/10.1016/j.foreco.2018.04.043 Jagodziński, A.M., Dyderski, M.K., Gęsikiewicz, K., Horodecki, P., Cysewska, A., Wierczyńska, S., & Maciejzyk, K. (2018). How do tree stand parameters affect young scots pine biomass? – allometric Equations and biomass conversion and expansion factors. Forest Ecology and Management, 409, 74- 83. https://doi.org/10.1016/j.foreco.2017.11.001 Jenkins, J.C., Chojnacky, D.C., Heath, L.S., & Birdsey, R.A. (2003). National-scale biomass estimators for United States tree species. Forest Science, 49(1): 12-35. Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 11 Jenkins, J.C., Chojnacky, D.C., Heath, L.S., & Birdsey, R.A. (2004). Comprehensive database of diameter-base biomass regressions for North American tree species [General Technical Report NE-319]. Newtown Square, PA, USA: USDA Forest Service, Northeastern Research Station, 45 p. Kang, X.G. (2011). Forest Management, 4th edition. Beijing: China Forestry Publishing House. Lee, W.K., Gadow, K.V., Chung, D.J., Lee, J.L., & Shin, M.Y. (2004). DBH growth model for Pinus densiflora and Quercus variabilis mixed forests in central Korea. Ecological Modelling, 176(1-2),187-200. https:// doi.org/10.1016/j.ecolmodel.2003.11.012 Lhotka, J.M., & Loewenstein, E.F. (2011). An individual- tree diameter growth model for managed uneven-aged oak-shortleaf pine stands in the Ozark Highlands of Missouri, USA. Forest Ecology and Management, 261(3),770-778. https://doi. org/10.1016/j.foreco.2010.12.008 Li, C.M.,& Zhang, H.R. (2010). Modeling dominant height for Chinese fir plantation using a nonlinear mixed- effects modeling approach. Scientia Silvae Sinicae, 46(3), 89-95. Li, C.M. (2012). The simultaneous Equation system of total volume in fir plantation. Scientia Silvae Sinicae, 48(6), 80-88. Li, J. (2011). Dynamics of biomass and carbon stock for young and middle aged plantation of Simao pine (Pinus kesiya var. langbianensis). Beijing: Beijing Forestry University. Li, Y.X., & Jiang, L.C. (2013). Modeling wood density with two-level linear mixed effects models for Dahurian larch. Scientia Silvae Sinicae, 49(7), 91-98. https:// doi.org/10.11707/j.1001-7488.20130713 Liu, Z.G., & Li, F.R. (2003). The generalized Chapman- Richards function and applications to tree and stand growth. Journal of Forestry Research, 14(1) 19-26. https://doi.org/10.1007/BF02856757 Luo, X., He, H.S., Liang, Y., Wang, W.J., Wu, Z., & Fraser, J.S. (2014). Spatial simulation of the effect of fire and harvest on aboveground tree biomass in boreal forests of Northeast China. Landscape Ecology, 29(7),1187-1200. https://doi.org/10.1007/ s10980-014-0051-x Mctague, J.P., & Weiskittel, A.R. (2016). Individual-tree competition indices and improved compatibility with stand-level estimates of stem density and long-term production. Forests, 7(10), 238. https:// doi.org/10.3390/f 7100238 Meng, X.Y. (2006). Forest mensuration (3rd ed). Beijing: China Forestry Publishing House. Muukkonen, P. (2007). Forest inventory-based large- scale forest biomass and carbon budget assessment: new enhanced methods and use of remote sensing for verification. Helsinki, Finland: Department of Geography, University of Helsinki. Njana, M.A., Bollandsås, O.M., Eid, T., Zahabu, E., & Malimbwi, R.E. (2016). Above- and belowground tree biomass models for three mangrove species in Tanzania: a nonlinear mixed effects modelling approach. Annals of Forest Science, 73(2),353-369. https://doi.org/10.1007/s13595-015-0524-3 Ou, G.L., Wang, J.F., Xu, H., Chen, K.Y., Zheng, H.M., Zhang, B., Sun, X.L., Xu, T.T., & Xiao, Y.F. (2016). Incorporating topographic factors in nonlinear mixed-effects models for aboveground biomass of natural Simao pine in Yunnan, China. Journal of Forestry Research, 27(1), 119-131. https://doi.org/10.1007/s11676- 015-0143-8 Pinheiro, J.C., & Bates, D.M. (2000). Mixed effects models in S and S-plus. New York: Springer Verlag. Ramon, C.L., George, A.M., Walter, W.S., Russell, D.W., & Oliver, S. (2006). SAS for mixed models, 2nd ed. Cary, NC, USA: SAS Institute Inc. Razali, W. W., Razak, T. A., Azani, A. M., & Kamziah, A.K. (2015). Mixed-effects models for predicting early height growth of forest trees planted in Sarawak, Malaysia. Journal of Tropical Forest Science, 27(2), 267-276. https://doi.org/10.2307/43582392 Richards, F.J. (1959). A flexible growth function for empirical use. Journal of Experimental Botany, 10(29), 290-300. https://doi.org/10.1093/ jxb/10.2.290 Robichaud, E., & Methven, I.R. (1993). The effect of site quality on the timing of stand breakup, tree longevity, and the maximum attainable height of black spruce. Canadian Journal of Forest Research, 23(8), 1514-1519. https://doi.org/10.1139/x93- 191 Rohner, B., Bugmann, H., & Bigler, C. (2013). Estimating the age-diameter relationship of oak species in Switzerland using nonlinear mixed-effects models. European Journal of Forest Research, 132(5-6),751- 764. https://doi.org/10.1007/s10342-013-0710- 5 Running, S.W.,& Mcleod, S.D. (1988). Comparing site quality indices and productivity in ponderosa pine stands of western Montana. Canadian Journal of Forest Research, 18(3), 346-352. https://doi. org/10.1139/x88-052 Sahney, S., Benton, M.J., & Ferry, P.A. (2010). Links between global taxonomic diversity, ecological diversity and the expansion of vertebrates on land. Biology Letters, 6(4), 544-547. https://doi. org/10.1098/rsbl.2009.1024 Saint-Andre, L., M’Bou, A.T., Mabiala, A., Mouvondy, W., Jourdan, C., Roupsard, O., Deleporte, P., Hamel, O., & Nouvellon, Y. (2005). Age-related Equations for above- and below-ground biomass of a Eucalyptus hybrid in Congo. Forest Ecology and Management, 205(1-3), 199-214. https://doi.org/10.1016/j. foreco.2004.10.006 Sleeter, B.M., Liu, J., Daniel, C., Rayfield, B., & Loveland, T.R. (2018). Effects of contemporary land-use and land-cover change on the carbon balance of terrestrial ecosystems in the united states. Environmental Research Letters, 13(4). https://doi. org/10.1088/1748-9326/aab540 Sileshi, G.W. (2014). A critical review of forest biomass estimation models, common mistakes and corrective measures. Forest Ecology and Management, 329, 237-254. https://doi. Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 12 org/10.1016/j.foreco.2014.06.026 Smith, A., Granhus, A., Astrup, R., Bollandsås, O.M., & Petersson, H. (2014). Functions for estimating aboveground biomass of birch in Norway. Scandinavian Journal of Forest Research, 29(6), 565-578. https://doi.org/10.1080/02827581.201 4.951389 Somogyi, Z., Cienciala, E., Mäkipää, R., Muukkonen, P., Lehtonen, A., & Weiss, P. (2007). Indirect methods of large-scale forest biomass estimation. European Journal of Forest Research, 126(2), 197-207. https://doi.org/10.1007/s10342-006-0125-7 Southwest Forestry College, & Forestry Department of Yunnan Province. (1988). Iconographia Arbororum Yunnanicorum. Kunming, China: Yunnan Science and Technology Press. Sturtevant, B.R., & Seagle, S.W. (2004). Comparing estimates of forest site quality in old second-growth oak forests. Forest Ecology and Management, 191(1-3), 0-328. Subedi, N., & Sharma, M. (2011). Individual-tree diameter growth models for black spruce and Jack pine plantations in northern Ontario. Forest Ecology and Management, 261(11), 2140-2148. https:// doi.org/10.1016/j.foreco.2011.03.010 Talhelm, A.F., Pregitzer, K.S., Kubiske, M.E., Zak, D.R., Campany, C.E., Burton, A.J., Dickson, R.E., Hendrey, G.R., Isebrands, J.G., Lewin, K.F., Nagy, J, & Karnosky, D.F. (2014). Elevated carbon dioxide and ozone alter productivity and ecosystem carbon content in northern temperate forests. Global Change Biology, 20(8), 2492-2504. Temesgen, H., Affleck, D., Poudel, K., Gray, A., & Sessions, J. (2015). A review of the challenges and opportunities in estimating above ground forest biomass using tree-level models. Scandinavian Journal of Forest Research, 30(4), 326-335. https:// doi.org/10.1080/02827581.2015.1012114 Ter-Mikaelian, M.T., & Korzukhin, M.D. (1997). Biomass Equations for sixty-five North American tree species. Forest Ecology and Management, 97(1), 1–24. https://doi.org/10.1016/S0378- 1127(97)00019-4 Timilsina, N., & Staudhammer, C.L. (2013). Individual tree-based diameter growth model of slash pine in Florida using nonlinear mixed modeling. Forest Science, 59(1), 27-37. https://doi.org/10.5849/ forsci.10-028 Vanninen, P., Ylitalo, H., Sievänen, R., & Mäkelä, A. (1996). Effects of age and site quality on the distribution of biomass in Scots pine (Pinus sylvestris L.). Trees, 10(4): 231-238. https://doi.org/10.1007/ BF02185674 Wang, H.L. (2003). Study on stand growth models for Pinus kesiya var. langbianensis natural secondary woodland. Kunming, China:Southwest Forestry College. Waring, R.H., Milner, K.S., Jolly, W.M., Phillips, L., & Mcwethy, D. (2006). Assessment of site index and forest growth capacity across the pacific and inland Northwest USA with a MODIS satellite-derived vegetation index. Forest Ecology and Management, 228(1-3), 1-291. Wen, Q.Z., Zhao, Y.F., Chen, X.M., Yang, Z.X., Ai, J.L., & Yang, X.S. (2010). Dynamic study on the values for ecological service function of Pinus kesiya forest in China. Forest Research, 23(5), 671-677. West, B., Welch, K.B., & Galecki, A.T. (2007). Linear mixed models: a practical guide using statistical software. Boca Raton FL, USA: Chapman and Hall/CRC. Westfall, J.A. (2016). Strategies for the use of mixed- effects models in continuous forest inventories. Environmental Monitoring and Assessment, 188(4), 1-11. https://doi.org/10.1007/s10661-016-5252- 0 Woodwell, G.M., Whittaker, R.H., Reiners, W.A., Likens, G.E., Delwiche, C.C., & Botkin, D.B. (1978). The biota and the world carbon budget. Science, 199(4325), 141-146. https://doi.org/10.1126/ science.199.4325.141 Wu, Z.L.,& Dang, C.L. (1992). The biomass of Pinus kesiya var. langbianensis stands in Pu-er district, Yunnan. Journal of Yunnan University: Natural Science Edition, 14(2), 161-167. Yue, F.,& Yang, B. (2011). Study on carbon sink of Pinus kesiya forests. Jiangsu Agricultural Sciences, 39(5), 467-469. Zeide, B. (2005). How to measure stand density. Trees,19(1), 1-14. Zeng, W.S., Zhang, H.R., & Tang, S.Z. (2011). Modeling method on the standing trees. Beijing: China Forestry Publishing House. Zhang, Y., & Borders, B.E. (2004). Using a system mixed-effects modeling method to estimate tree compartment biomass for intensively managed loblolly pines - an allometric approach. Forest Ecology and Management, 194, 145-157. https:// doi.org/10.1016/j.foreco.2004.02.012 Zianis, D., Muukkonen, P., Mäkipää, R., & Mencuccini, M. (2005). Biomass and stem volume Equations for tree species in Europe. Silva Fennica Monographs, 4(4). Nong et al. New Zealand Journal of Forestry Science (2019) 49:11 Page 13