Height-diameter Models of Chinese Fir ( Cunninghamia lanceolata ) Based on Nonlinear Mixed Effects Models in Southeast China

Tree height and diameter at breast height are two important forest factors. The best model from 23 heightdiameter equations was selected as the basic model to fit the height-diameter relationships of Chinese fir with one level (sites or plots effects) and nested two levels (nested effects of sites and plots) Nonlinear Mixed Effects (NLME) models. The best model was chosen by smaller Bias, RMSE and larger . Then the best random-effects combinations for the NLME models were determined by AIC, BIC and -2LL. The results showed that the basic model with three random effects parameters φ1, φ2 and φ3 was considered the best mixed model. The nested two levels NLME model considering heteroscedasticity structure (power function) possessed with higher predictable accuracy and significantly improved model performance (LRT = 469.43, p<0.0001). The NLME model would be allowed for estimating accuracy the height-diameter relationships of Chinese fir and provided better height predictions than the models using only fixed-effects parameters.


INTRODUCTION
Chinese fir (Cunninghamia lanceolata (Lamb.)Hook) is the most important afforestation tree species, which has fast-growth and good wood qualities in southeast China.According to the National Continuous Forest Inventory, about 11.26 million ha and 734.09 million m 3 of Chinese fir were distributed more than 10 provinces in China in 2010, making it the major species in terms of volume harvested.However, as a result of great variability exists among stands in terms of silvicultural and ecological conditions, it needs proper management, especially management tools.
Knowledge of the individual tree height (h) and diameter at breast height outside bark (dbh, measured 1.30 m above ground level) is fundamental for developing growth and yield models in forest stands.However, measuring diameter is more accurate, quicker and cheaper than measuring tree height (Colbert et al., 2002).As a result, in forest inventories, diameter is measured for all the sampled trees, but height is measured only for a subsample of trees (Calama and Montero, 2004) and as h and dbh are correlated, it is common practice to fit height-diameter (h-d) models to predict h from measured dbh (Crecente-Campo et al., 2010b).
The relationship between h and dbh varies from stand to stand (Curtis, 1967), which is usually determined by measurements of both variables from trees growing in different stands or site class of the sample plot.This hierarchical structure (i.e., plots within a site class) makes the measurements dependent (Calama and Montero, 2004;West et al., 1984).As a result, there is spatial correlation among measurements in the sampling unit (Fox et al., 2001).However, the stochastic structure is often ignored and independence between measurements is assumed (Biging, 1985;Fitzmaurice et al., 2004;Gutzwiller and Riffell, 2007;Keselman et al., 1999).Therefore, the Ordinary Nonlinear Least Squares (ONLS) regression has biased estimates of the standard error of parameter estimation (Schabenberger and Gregoire, 1995).Many researchers take the Nonlinear Mixed-Effects (NLME) models into account to develop h-d models (Adame et al., 2008a;Calama and Montero, 2004;Crecente-Campo et al., 2010a;Sharma and Parton, 2007).
The fixed-effects parameters in NLME models are associated with an entire population, as in traditional regression, while the random-effects parameters are associated with individual experimental units (Pinheiro and Bates, 2000).Because of their versatility, the mixed-effect modeling approach is more useful than ONLS fitting.Use of NLME models allows h-d models to be calibrated by predicting random components from plot-or site-level covariates when a new subject is available (one not used in the fitting of the model), using the Empirical Best Linear Unbiased Predictors (EBLUPs) (Adame et al., 2008b;Calama and Montero, 2005;Nigh, 2012).
The research aimed to develop an h-d model for C. lanceolata (Lamb.)Hook in Fujian province, southeast China.A nested two levels nonlinear mixed modeling approach which included both fixed and random components was applied to the hierarchical structure of the dataset to diminish the level of variance among sampling units were included as covariates.In this study, we considered one level and nested two levels to build NLME h-d models, the first level is site and the second level is plot (nested within site).Our preliminary analysis showed that NLME models with random effects effectively removed the heteroscedasticity could provide important tools for sustainable management for this species in the study area.The predictive ability of the developed model and the applicability of NLME model technique were demonstrated on a separate validation dataset.
Data used in this study were collected from 29 sample plots of Chinese fir even trees from 4 sites which located in 4 regions, Qiantan, Shuinan, Yuhua and Yuandang work areas, which were numbered I, II, III and IV, respectively (Fig. 1).These sample plots were established from 2010 to 2013 and square in sharp but vary in size from 400 to 600 m 2 .All standing live trees (height>1.3m) on the plots were measured for breast-height diameter over bark (at 1.3 m above ground), total tree height.These data were randomly divided into two datasets with 75% used for model fitting and 25% used for model validation.Summary statistics for both fitting and validation datasets are shown in Table 1.

Methods:
Nested two levels NLME model: Available data were based on a sample of multiple measurements (diameter at breast height-height) taken from different plots located in different regions.As a consequence of this nested structure, we detected high correlation among measurements taken from the same plot.To alleviate this, a mixed-model approach has been proposed by other authors (Crecente-Campo et al., 2010b;Adame et al., 2008a, b;Sharma and Parton, 2007).A general expression for a NLME model can be defined as (Lindstrom and Bates, 1990;Vonesh and Chinchilli, 1997): ( )   ( ) ϕ1, ϕ2 and ϕ3: The parameters of models where, M = The number of sites M i = The number of plots within the i th site n ij = The number of measurements h ijk = The height of the k th tree in the j th plot within the i th site d ijk = The diameter at breast height ϕ ijk = A parameter vector r×1 = (where r is the number of parameters in the model) f = A nonlinear function of the predictor variables and the parameter vector ε ijk = The within-group error including the withingroup variance (Davidian and Giltinan, 1995) which is assumed normally distributed with zero expectation and a positive-definite variancecovariance structure R ijk , generally expressed as a function of the parameter vector δ (Meng and Huang, 2009): Moreover, ϕ ijk can be expressed as: ( , where, λ is the p×1 vector of fixed population parameters (where p is the number of fixed parameters in the model), µ i and µ ij are the q 1 ×1 and q 2 ×1 vectors of random effects associated with the first level and the second level, respectively (where q 1 and q 2 are the numbers of random parameters of two levels in the model), which are assumed to be normal (or Gaussian) with a mean of 0 and have the variance-covariance matrixes ψ i and ψ ij , which are the q 1 ×q 1 and q 2 ×q 2 variance-covariance matrixes of random effects in two levels, respectively.A ijk , B i,jk and B ijk are design matrices of size r×p, r×q 1 and r×q 2 for the fixed and random effects of two levels specific to each sampling unit, respectively.
Height-diameter equations: Several functions to model the relationship between dbh and h of the trees in a stand are available (Curtis, 1967;Fang and Bailey, 1998;Huang et al., 2000Huang et al., , 1992)).In this study, 23 functions (Table 2) were compared to determine which showed the best predictive ability and fit to the heightdiameter data (Huang et al., 2000).All functions tested have two or three parameters and are nonlinear.
Using the R nls function to fit 23 above-mentioned equations by ONLS regression.Different initial values for the parameters were used to ensure that a global minimum was achieved.The best function was selected by applying three statistical criteria, which are Bias, Root Mean Square Error (RMSE) and the adjusted coefficient of determination (ܴ ௗ ଶ ) (Zhang et al., 2011).
The calculation formulas of these statistics were listed as follow: ( ) where, y ijk and ‫ݕ‬ ො = The height measurements and predictions, respectively ‫ݕ‬ ത = The average of measurements Mixed parameter ascertainment: Which parameters in the model should be considered as random effects and which can be treated as purely fixed effects is the key question to fit the NLME models.Pinheiro and Bates (2000) suggested that it could start with a model with random effects for all parameters and then examine the fitted object to decide which, if any, of the random effects can be eliminated from the model by Akaike's information criterion (Akaike, 1974), Bayesian information criterion (Weiss, 2005) and -2 logarithm likelihood (-2LL).An appropriate variance function structure for NLME models were determined by Likelihood Ratio Test (LRT) (Fang and Bailey, 2001;Pinheiro and Bates, 2000).All NLME models presented in this research were calibrated using the NLME function in the R statistical software (Ihaka and Gentleman, 2004;Pinheiro et al., 2005).

Determining the structure of variance-covariance:
The eigenvalues of variance-covariance matrixes ψ i and ψ ij are strictly positive (Pinheiro and Bates, 2000).
A hypothetical 3×3 variance-covariance matrix is shown as follow (Calama and Montero, 2004;Fu et al., 2013): where, ߪ ௨ ଶ , ߪ ௩ ଶ and ߪ ௪ ଶ are the variance for the random effects u, v and w, respectively, σ uv = σ vu is the covariance among random effects u and v, σ uw = σ wu is the covariance among random effects u and w, σ vw = σ wv is the covariance among random effects v and w.
Determining the structure of R ijk : To account for within-plot heteroskedasticity in R ijk , which includes weighting factors (Davidian and Giltinan, 1995;Meng and Huang, 2009), the approach is given by Calama and Montero (2004)  where, σ 2 is a scaling factor for the error dispersion (Grégoire et al., 1995), G ij is a n ij ×n ij diagonal matrix within-plot error heteroscedasticity variances and I ij is an n ij ×n ij matrix showing the within-plot autocorrelation structure of error.In our case, because no correlation patterns emerged between measurements from the same plot, I ij reduced to a n ij ×n ij identity matrix.
In this study, the variance heterogeneity was removed by two variance functions: the power function and the exponential function (Pinheiro and Bates, 2000): where, δ is estimated parameter.

Parameter estimation:
The parameters in equation were estimated using the R NLME function (Lindstrom and Bates, 1990;Pinheiro and Bates, 2000).
In this study, the random effects parameters can be calculated with the validation data of height and diameter, by the Empirical Best Linear Unbiased Predictors (EBLUPs) (Vonesh and Chinchilli, 1997): where, D = The estimated variance-covariance matrix for the random-effects ܾ R ijk = The estimated variance-covariance matrix for the error term Z ijk = The estimated partial derivatives matrix with respect to random effects parameters

CONCLUSION
Twenty three height-diameter equations were considered in this study for describing the heightdiameter relationship for Chinese fir in pure plantation stands in Fujian province, southeast China.Then the best function Eq. ( 10) as the basic model was determined by Bias, RMSE and ܴ ௗ ଶ .The ONLS regression approach is commonly used to build the h-d models, but to a certain extent, there is a restriction of this method because of its assumption of independence and ignorance of heteroscedastiity.In this study, one level (site or plot) and nested two levels (plot nested within site) NLME models with variance function approach were used to estimate the height-diameter relationship based on the basic model because of the hierarchical structure of measurements data.We used the power function and the exponential function as the variance functions to eliminate the heteroscedasticity.The random effects parameters can be calculated with the measurements of the height and diameter (data not used in the fitting of the model), by EBLUPs Eq. ( 9).The results showed that the nested two levels NLME models Eq. ( 19) provided better model fitting and more precise estimations than the others Eq. ( 10), ( 17) and (18) (Table 5, Fig. 2 and 3).Therefore, NLME models could improve the fitting performance in modeling height-diameter relationship for Chinese fir, which is important for owners and managers of plantations to make appropriate sustainable management plans for different sites, even plots.

Fig. 1 :
Fig. 1: Twenty nine sample plots located in four regions in Fujian province, southeast China and Crecente-Campo et al. (2010a): (6) model for predicting the heightdiameter relationship of Chinese fir trees in a pure plantation in the study area.

Table 2 :
Mathematical expressions of the height-diameter equations

Table 3 :
Performance criteria of NLME models for combinations of random-effects

Table 5 :
Evaluation indices of each model