Multilevel Nonlinear Mixed-Effect Crown Ratio Models for Individual Trees of Mongolian Oak (Quercus mongolica) in Northeast China.

In this study, an individual tree crown ratio (CR) model was developed with a data set from a total of 3134 Mongolian oak (Quercus mongolica) trees within 112 sample plots allocated in Wangqing Forest Bureau of northeast China. Because of high correlation among the observations taken from the same sampling plots, the random effects at levels of both blocks defined as stands that have different site conditions and plots were taken into account to develop a nested two-level nonlinear mixed-effect model. Various stand and tree characteristics were assessed to explore their contributions to improvement of model prediction. Diameter at breast height, plot dominant tree height and plot dominant tree diameter were found to be significant predictors. Exponential model with plot dominant tree height as a predictor had a stronger ability to account for the heteroskedasticity. When random effects were modeled at block level alone, the correlations among the residuals remained significant. These correlations were successfully reduced when random effects were modeled at both block and plot levels. The random effects from the interaction of blocks and sample plots on tree CR were substantially large. The model that took into account both the block effect and the interaction of blocks and sample plots had higher prediction accuracy than the one with the block effect and population average considered alone. Introducing stand density into the model through dummy variables could further improve its prediction. This implied that the developed method for developing tree CR models of Mongolian oak is promising and can be applied to similar studies for other tree species.

Introduction

Tree crown size is an important variable that is commonly involved in growth and yield models used as decision-support tools in forest management [1-3]. Tree crown size is usually characterized using crown length (CL), crown width (CW), and crown ratio (CR). Tree CR is defined as the percentage of crown length from the base of live crown to the tree top to its total height. The value ranges from 0 (for the trees without crown, such as dead or defoliated trees) to 1 (for the trees with crowns extending the entire tree bole). Tree CR is considered to be an indirect measure of photosynthetic capacity of trees [2]. Moreover, tree CR is also a useful indicator of vigor [4-7] and stand density [8]. Additionally, it is a variable of interest in the management of many non-timber resources including recreation forests and wildlife habitats [9]. In forest inventory, however, measuring tree CR for all sampled trees is time- and money-consuming (e.g., [5, 10]). Therefore, developing accurate tree CR models is necessary, which allows forest managers to accurately predict tree CR.

One common approach used to obtain values of tree CR is to develop either deterministic or stochastic models from tree characteristics [1–3, 5, 10, 11]. So far, most of the obtained models (see Table 1) are simply linear or nonlinear (e.g., Exponential, logistic and Weibull models) and often developed as a function of diameter at breast height (1.3 m above ground) (D), total tree height (H), age (A) and so on, using ordinary least squares techniques [12, 13]. When the relationships of tree CR with other tree variables are modeled, measurements of the variables are often collected within sample plots that are allocated in different stands that represent different site conditions, called blocks in this study). This hierarchical structure leads to the fact that the selected trees are located within the same plots and the plots selected within the same blocks and that measurements are likely correlated with each other significantly [14-16]. In this case, using ordinary least square regression technique to fit the data often leads to bias as the standard errors of parameter estimates would get inflated [15]. To deal with this problem, several approaches have been proposed [14, 17, 18].

Table 1

Existing tree crown ratio (CR) models.

Model	Predictors	Model form	Source	Model no
CR = [1 + exp(βX)]-1	HT, D, HT/D,BAL, CCF, ELEV,SL,AZ	Logistic	Hasenauer and Monserud, 1996	II.1
CR = [1 + exp(βX)]-1	D,HT,TSC,CCF,BAL, ELEV, SL, ASPECT	Logistic	Temesgen et al., 2005	II.1
CR = [1 + exp(βX)]-1	CCR, D	Logistic	Toney and Reeves, 2008	II.1
CR = [1 + exp(βX)]-1	D, HT, BA	Logistic	Leites et al., 2009	II.1
CR = [1 − exp(βX)]-1	MHT, TSC	Logistic	Popoola and Adesoye, 2012	II.2
CR = a/[1 + b exp(−βX)]	Age, SD, DH, D	Logistic	Soares and Tomé, 2001	II.3
CR = a + exp(−βX)	MHT, TSC	Exponential	Popoola and Adesoye, 2012	II.4
CR = exp(βX)	D, HT, BA, CCF,PCT	Exponential	Leites et al., 2009	II.5
CR = a[1 – b exp(−βX)]	Age, SD, DH, D	Exponential	Soares and Tomé, 2001	II.6
CR = [1 − exp(−ϕ βX)]	Age, D, HT,	Exponential	Dyer and Burkhart, 1987	II.7
CR = [1 − exp(−ϕ βX)]	BA, DH, D, HT	Exponential	Hynynen, 1995	II.7
CR = [1 + exp(−βX)]−1/2	MHT, TSC	Richards	Popoola and Adesoye, 2012	II.8
CR = a/[1 + b exp(−βX)]1/c	Age, SD, DH, D	Richards	Soares and Tomé, 2001	II.9
CR = a/[1 − b exp(−βX c)]	Age, SD, DH, D	Weibull	Soares and Tomé, 2001	II.10

Note: MHT: merchantable height, TSC: tree slenderness coefficient, D: diameter at breast height, H: total tree height, BA: stand basal area, CCF: stand crown competition factor, PCT: The percentile in the stand basal area distribution, CCR: compacted crown ratio, BAL: basal area per ha for trees larger than the subject tree, ELEV: elevation, SL: slope, ASPECT: aspect, SD: stand density, DH: dominant tree height, AZ: azimuth of aspect in radians, Age: tree age, a, b, c: model parameters, β: parameter vector, x: vector of stand or tree variables.

One of the approaches is the use of nonlinear mixed-effect (NLME) models [19]. This approach provides an efficient way to analyze correlated hierarchical structure data and make accurate local predictions ([19, 20]. NLME models contain both fixed- and random-effect parameters. The fixed parameters provide the potential to account for covariate or treatment effects as in traditional regression, while the random parameters offer the capacity to explain various sources of heterogeneity and randomness in the data caused by known and unknown factors [20-25]. Therefore, the NLME models have been widely applied in forest growth and yield modelling over the past decades (e.g., [21-25]). They can be used not only to account for stochastic variability in tree CR models, but also to provide the potential of increasing accuracy of tree CR prediction. To the authors’ knowledge, however, only a few studies have been conducted to model tree CR using the NLME models, especially multilevel NLME models [26, 27].

The objective of this study was to develop an individual tree CR model for Mongolian oak (Quercus mongolica) natural forests where sample plots were selected and the trees within the plots were measured. A two-level NLME CR model with random effects at both block and plot levels was first built to account for autocorrelation of hierarchical structure data. Various tree and stand characteristics were then tested for their contributions to the improvement of the CR model fits and predictions. The predictive ability and applicability of the obtained NLME CR model were validated and demonstrated using an independent data set.

Materials and Methods

Study area and data

The used data were obtained from a total of 118 permanent sample plots (PSPs) established in Mongolian oak natural forests allocated in Wangqing Forest Bureau of northeast China (123°56′–131°04′ E, 43°05′–43°40′ N) (Fig 1, Table 2). These PSPs were nested within 15 blocks, out of which 12 blocks with 100 PSPs are located in Tazigou forest farm and the other 3 blocks with 18 PSPs in Jincang forest farm, that were randomly located in the study area to represent the stands with different site conditions. When the data were collected, there were no regulations setup to limit the scientific research in this forestry bureau and thus no specific permission was required for the field work. In addition, there were no any endangered or protected species involved in the field study.

Fig 1

Location (upper left) of study area: Wangqing Forest Bureau (upper right) in northeast China and spatial distribution of 15 blocks and 118 sample plots (bottom).

Table 2

Numbers of sample plots and trees grouped into classes of stand density index (SDI) for Mongolian oak.

Variable	Class	Class midpoint	Class range	Number of plots	Number of trees
Stand density index	1	150	0 < SDI ≤ 300	17	253
(trees ha-1)	2	350	300 < SDI ≤ 400	40	705
	3	450	400 < SDI ≤ 500	29	554
	4	550	500 < SDI ≤ 600	18	916
	5	600	SDI > 600	8	706

All the PSPs established in 2010 (100 plots) and 2013 (18 plots) by the Research Institute of Forest Resources Information Techniques, Chinese Academy of Forestry, were square in shape with an average size of 719 m2 and varying from 400 to 2500 m2. All standing live trees (H > 1.3 m and D > 5cm) within each of the plots were measured for D, H and trunk height to the crown base (HCB). Tree HCB was defined as the height from the ground to the base of the first normal green branch as a part of the crown; this excluded the secondary branches (epicormic and adventitious) [5]. Furthermore, a single green branch was not the base of tree crown if there were at least three dead whorls above it [5]. Three to five dominant or sub-dominant trees within each plot were chosen to measure plot dominant tree height (DH) and dominant tree diameter at breast height (DD) [28]. Stand age (A) was determined by the mean values of three average sample trees in each plot [29]. Canopy density (CD) for each plot was obtained using a moose-horn method [30]. The measurements of the tree variables were collected from a total of 3685 trees (Table 2).

Stand density (SD) affects tree growth significantly [3, 31, 32]. Our preliminary analyses revealed that the differences of tree CR among different SD values of the PSPs were significant at a risk level of 0.05. All the PSPs were thus stratified into five stand density classes by stand density index (SDI) [33] (Table 2) and the effect of SDI on tree CR was then analyzed based on the inventory data. The expression of SDI is given as [33]: where SD is stand density (trees ha-1), MD is stand arithmetic mean diameter (cm), D 0 is standard mean diameter (cm), β is an estimated parameter, and ε is an error term. The values of D 0(D 0 = 20) and β (β = 1.3798) in this study were obtained based on the study by Du et al. (2000) in which stand density index models for 11 tree species groups for Wangqing Forest Bureau, including Mongolian oak, were developed. The effect of Meyer’ site index (Meyer’ SI) on tree CR was also analyzed in this study. The Meyer’ SI was given by Du et al. [29]: where α is a parameter, A is stand age (years) at time when DH measurements were collected, and A 0 is the standard age (base age) (years) of SI, equal to 40 for Mongolian oak [29]. A multivariate analysis was carried out to detect outlier data based on the distribution of Mahalanobis distance between the observations and their expectations [15]. Only 6 PSPs were rejected from this technique as they represented special cases, namely very low dense stands or over mature stands. The remaining 112 PSPs were randomly divided into two groups: one for model fitting and the other for model validation. The data used for model fitting consisted of 2166 trees from 74 PSPs and the data for model validation was composed of 968 trees from 38 PSPs (Table A in S1 File). Summary statistics of the data and relevant stand characteristics are listed in Table 3.

Table 3

Summary statistics of stand variables from the sample plots used for model fitting and validation, respectively.

Data	Variable	Min	Max	Mean	SD	CV%
Plots for model fitting
	Area (m2)	400	2500	719	581	80.83
	CR	0.06	1.00	0.61	0.17	28.49
	HCB (m)	0.40	13.30	4.60	2.39	51.86
	D (cm)	1.40	70.10	15.42	8.40	54.49
	H (m)	1.80	25.60	12.18	4.18	34.29
	SD (trees ha-1)	275	1863	877	515	59
	CD	0.46	0.90	0.79	0.09	10.97
	DH (m)	12.54	23.78	17.04	2.57	15.06
	DD (cm)	16.75	38.90	24.53	5.08	20.69
Plots for model validation
	Area (m2)	0.04	0.25	0.08	0.07	92.32
	CR	0.01	0.93	0.60	0.18	29.15
	HCB (m)	0.50	14.00	4.76	2.43	51.04
	D (cm)	1.50	48.20	16.48	8.65	52.50
	H (m)	2.20	25.60	12.49	4.14	33.16
	SD (trees ha-1)	300	1575	662	381	58
	CD	0.46	0.90	0.79	0.10	12.17
	DH (m)	12.43	22.98	17.51	1.90	10.88
	DD (cm)	16.88	34.55	25.48	5.32	20.89

Note: Area: sample plot area, CR: crown ratio, HCB: trunk height to crown base, D: diameter at breast height, H: tree height, SD: stand density, CD: canopy density, DH: dominant tree height, DD: dominant tree diameter.

Modeling approach

Two-level NLME models

Tree CR mixed-effect models with two-level random effects were formulated according to multilevel nonlinear mixed-effect model techniques in this study [34]. That is, blocks were used as the first level random effect and sample plots were nested in each of block, showing interaction, as the second level random effect: where CR is the crown ratio of the k tree on the j plot nested within the i block, M is the number of blocks, M is the number of plots within the i block, n is the number of observations in the j plot of the i block, and f(·) is a real-valued and differentiable function of a group-specific parameter vector and a covariate vector t . The within-group error ε that accounts for within-group variance and correlation was assumed to follow a normal distribution that has an expectation of zero and a positive-definite variance-covariance structure R [35]. R is generally expressed as a function of the parameter vector λ [24]:

Moreover, can be expressed as [16]: where β is a p-dimensional vector of fixed effects, meaning the first-level random effects, u is independently normal distributed q 1-dimensional vector with zero means and variance- covariance matrix Ψ 1, indicating the second-level random effects, u is independently normal distributed q 2-dimensional vector with zero means and variance-covariance matrix Ψ 2, and A , B and M are design matrices, and u , u and ε are independent of each other.

Predictor variables

Growth of individual trees is potentially affected by three groups of variables: tree size and vigor effects, site condition effects and competition [36-39]. In this study, a total of 14 variables, including 2 tree size and vigor effects related variables, 1 site condition effects related variable, and 11 competition effects related variables (Table 4), were used to evaluate their effects on tree CR.

Table 4

A total of 14 candidate variables used for developing crown ratio model.

Groups of variables	Variables
Tree size and vigor effects	diameter at breast height (D), total tree height(H)
Site condition effects	Meyer’ site index
Competition effects	stand density (SD), canopy density (CD), dominant tree height (DH), plot arithmetic mean diameter (AMD), plot dominant tree diameter (DD), plot quadratic mean diameter (QMD), number of trees with diameter larger than the target tree (LDN), total diameter of all trees with diameter larger than the target tree (LDTD), mean diameter of all trees with diameter larger than the target tree (LDMD), total basal area of all trees with diameter larger than the target tree (LDTBA), and mean basal area of all trees with diameter larger than the target tree (LDMBA)

The tree and stand variables were selected using a stepwise regression procedure. The measurements of the variables were first analyzed by graphical examination and correlation statistics [36]. Different combinations of the variables and their logarithmic transformations in linear models were then tested based on the coefficient of determination R . All the calculations were performed using R/S-Plus nls function [40].

Model selection

The base model used to develop the NLME model of tree CR was determined from a total of 10 candidate models (Table 1) with selected predictors based on the performance of model fitting and prediction. The models were first fit to the data from 2166 trees of 74 PSPs and their predictions were then compared with the observations from 968 trees of 38 PSPs for the validation dataset. Nonlinear regressions were carried out using ordinary nonlinear least square (ONLS) technique with the R/S-Plus nls function [40]. The following four statistical criteria were used to select the model that had the highest accuracy for both fitting and prediction [25, 41]: where CR and are the observed and predicted values of tree CR for the i observation (i = 1,…,N), N is the total number of observations, is the mean of the observed tree CR values, is the mean prediction error, ξ is the variance of prediction errors, δ is the root mean square error, R 2 is the coefficient of determination. δ combines the mean prediction error () and the variance of prediction errors (ξ), giving a robust measure of the overall model accuracy, and therefore was selected as a primary criterion for model evaluation [25]. The selected model with predictor variables was then used as a base model to construct NLME model of tree CR.

Construction of parameter effects

Several approaches were proposed to determine which parameters should be modeled as mixed effects in the base model [21, 42, 43]. In this study, we fitted all possible combinations of random effects on parameters including intercept and slope terms in the base model and selected the one with the smallest Akaike’s information criterion (AIC) and the largest log-likelihood (Loglik) [34]. To avoid over parameterization, the likelihood ratio test (LRT) was also used [16, 43].

Determining Ψ1 and Ψ2 structure

The variance-covariance matrices for the random effects, Ψ and Ψ , define the variability that exists among the sample plots and blocks, respectively. As in the study of Calama and Montero [44], both Ψ and Ψ are assumed unstructured in this study.

Determining R structure

To account for the within-plot heteroscedasticity and autocorrelation in R [24, 35], we applied the approach suggested by Davidian and Giltinan [35]: where σ 2 is a scaling factor for error dispersion [21], given by the value of residual variance of the estimated model, G is a n × n diagonal matrix explaining the variance of within-plot heteroscedasticity and Γ is a n × n matrix accounting for the within-plot autocorrelation structure of errors.

Three widely used variance models, including exponential model Eq (8), power model Eq (9) and constant plus power model Eq (10) [16] each with different predictors, were used to account for heterogeneity in variance. The LRT and AIC were employed to determine an appropriate variance model. where x is one of the selected predictors; γ, γ 1 and γ 2 are the estimated parameters; is an error term. In addition, two autocorrelation structures, AR (1) (autoregressive process of order one) and ARMA (autoregressive moving average process) for matrix Γ , were evaluated, and the one that had a higher fitting accuracy (i.e., smaller AIC) and provided the expected residual pattern was selected.

Parameter estimation

The parameters in Eq (1) were estimated by maximum likelihood (ML) estimation method using the Lindstrom and Bates (LB) [19] algorithm which alternates between two steps: a penalized nonlinear least square (PNLS) step and a linear mixed-effect (LME) step [16, 19]. This algorithm was implemented using the R/S-Plus nlme function [16] and for its details, readers can refer to Lindstrom and Bates [19] and Pinheiro and Bates [16].

Model prediction and evaluation

When two-level NLME models are used to predict tree CR, both population average (PA) and subject-specific (SS) responses [15, 34, 45] are often considered. The former is related to the fixed effect response and prediction of tree CR in the stands where independent stand or tree variables required by the models are measured, but tree CR measurements are not collected. The latter is related to the predictions of tree CR in the stands where in addition to the independent variables, the dependent variable tree CR is also measured in a sub-sample of trees. In order to reduce both measurement cost and potential errors, as Calama et al [15] and Temesgen et al [46] suggested, four randomly selected trees within each plot were measured for estimation of the random effect parameters in this study [15, 46].

Population average model

A population average (PA) model means that it contains the fixed effects as global parameters and the random effects as zero and has the following general form [34]: Where are the predicted values of tree CR for the k th tree in the j th plot in the i th block, t are other stand and tree variables for the corresponding tree, and is the estimates of PA model parameters and contains .

Subject-specific response at block level

When a model accounts for the effects of block on tree CR and ignores other effects such as the interaction of blocks and sample plots, its parameter estimates in Eq (11) contain both and . Usually, the random effect parameter estimates are obtained by an empirical best linear unbiased prediction (EBLUP) approach [20, 25, 34]. where the superscript T denotes the matrix transpose operation, is the q 1 × 1 dimensional random effects from blocks, is the q 1 × q 1 dimensional variance-covariance matrix, is the n i × q 1 dimensional design matrix of the partial derivatives of the nonlinear function f(·) with respect to random effect parameters u , . The values of blocks in the data used for both model fitting and validation were identical. Thus, the values of do not need to be recalculated in prediction. That is, they would be equal to the estimates obtained when the model is fitted.

Subject-specific response at both plot and block levels

In addition to , in Eq (11) also contains both and . That is, it simultaneously incorporates the block effects and the interaction of blocks and sample plots. The random effect parameters were calculated by Eq (12). The specific structure of , and (i = 1,…,M) in the Eq (12) in this case are presented in Appendix 1.

Model assessment

Statistics , ξ, δ and R 2 calculated by Eqs (3)–(6), and LRT were applied to assess the predictive ability of the developed tree CR models using both fitting and validation datasets. The obtained mixed-effect models were compared based on Loglik, AIC and LRT.

Results

Selection of base model

To avoid over-parameterization and collinearity in the models, only those variables that had statistically significant contributions to improving the quality of the models fit to the data were selected. The obtained variables were: D, DH and plot dominant tree diameter (DD). Especially, the variables D and DH were strongly correlated with tree CR and included in most of the candidate models in Table 1. All the candidate models with selected predictors converged except II.9 and II.10 in Table 1. The statistics for accuracy measures Eqs (3)–(6) to assess the performance of converged candidate models for both fitting and validation are calculated (Table B in S1 File). It is found that the values of the same statistics were very similar to each other among the models. For the data used to fit the models, the values of for all the candidate model were very close to 0, and the values of both ξ and δ had the same range of 0.1681 to 0.1687. The coefficients of determination varying from 0.2469 to 0.2554 were very small. For the data used to validate the models, the values of , ξ and δ for all the models had small ranges of -0.0174 to -0.0147, 0.1712 to 0.1717, and 0.1718 to 0.1725, respectively. But, the model validation results showed that model II.1 had the smallest values of mean prediction error, , variance of prediction error, ξ = 0.1712, and root mean square error, δ = 0.1718, suggesting a slightly higher prediction accuracy than the others. Additionally, another advantage for model II.1 is that its predictions are constrained to be between 0 and 1, i.e., 0% to 100% CR. Thus, the base model selected to construct the CR mixed-effect model is: where D are the diameter at breast height (cm) of the k th tree in the j th plot in the i th block stand density index class, DH and DD are the plot dominant tree height (m) and plot dominant tree diameter at breast height (cm), respectively, of the j th plot of the i th block stand density index class, and ϕ 0, ϕ 1, ϕ 2 and ϕ 3 are the model parameters and other variables as defined previously.

Mixed-effect model

Considering four model parameters (ϕ 0-ϕ 3) involved and both block effects and interaction of blocks and sample plots, a total of 15 different combinations of random effects were obtained for Eq (13). The mixed-effect model with each of the combinations was fitted to the data. It was found that only 10 mixed-effect model alternatives converged (Table C in S1 File). The following model Eq (14) with block effects and the interaction showed the smallest AIC, -1919.58, and the largest Loglik, 970.79 (Table C in S1 File): where β 0, β 1, β 2 and β 3 are fixed-effect parameters, u 1i and u 3i are random-effect parameters caused by blocks on ϕ 1 and ϕ 3, respectively, u 1ij and u 3ij are random-effect parameters caused by the interaction of blocks and sample plots on ϕ 1 and ϕ 3, respectively.

The residuals from Eq (14) were graphed against the predicted values for the data used for model fitting in Fig 2. The absolute values of the residuals tended to decrease as the observed tree CR values increased. It indicated that some heteroscedasticity remained even after incorporating the effects of blocks and the interaction effects of blocks and sample plots in the nonlinear mixed-effects CR model (14). Furthermore, the empirical autocorrelation function (ACF) for Eq (14) indicated that the autocorrelation was significant at a risk level of 0.05 among residuals within plots.

Fig 2

The residuals of predicted crown ratio values from Eq (14) graphed against the predicted values for Mongolian oak in Wangqing Forest Bureau of northeast China.

Within-plot variance-covariance (R) structure

The assessment statistics based on three variance models with each selected predictor (D, DH, and DD) for Eq (14) were shown in Table 5. For comparison purpose, the assessment results of this model were also derived and given when the variances of the error term ε was assumed homogeneous. Power model Eq (8) with DH as a predictor failed to converge. In the case of homogeneous variance, whether D or DD or DH was used as a predictor, the values of AIC and Loglik were significantly different from those using power model Eq (8), exponential model Eq (9) and constant plus power model Eq (10) at a risk level of 0.05. Even with random effects in the parameters, heteroscedasticity still existed in the mixed-effect model Eq (14) for tree CR. Among these variance models, exponential model Eq (9) with DH as a predictor showed the highest accuracy for model fitting (Table 5) (AIC = -1987, Loglik = 1006). Eq (14) with the Autocorrelation structure AR (1) produced a smaller AIC (AIC = -1950) than that with the ARMA (1, 1) (AIC = -1931). Therefore, Eq (14) was fitted with an AR (1) autocorrelation structure and an exponential variance model. The final nonlinear mixed-effects CR model was: where ρ is an estimated parameter in the matrix Γ ij with the autocorrelation structure AR (1), and other variables and coefficients in this model defined as above.

Table 5

Performance assessment of mixed-effect model Eq (14) using measurements of crown ratio with different variance models.

Variance	D				DD				DH
model	AIC	Loglik	LR	p value	AIC	Loglik	LR	p value	AIC	Loglik	LR	p value
1	-1920	971			-1920	971			-1920	971
PF	-1965	994	47.26	<0.0001	-1965	994	46.95	<0.0001	None		F	F
EF	-1970	997	52.02	<0.0001	-1957	991	39.72	<0.0001	-1987	1006	70.03	< 0.0001
CPF	-1963	994	47.25	<0.0001	-1963	994	46.97	<0.0001	-1962	994	46.97	< 0.0001

Note: D: diameter at breast height, DD: dominant tree diameter, DH: dominant tree height, AIC: Akaike information criterion, Loglik: log-likelihood, LR: likelihood ratio, Variance model 1 means that the variances are homogeneous, PF: power model—Eq (8), EF: exponential model—Eq (9), CPF: constant plus power model—Eq (10), F denotes a model failing to converge.

Parameter estimation

Table 6 shows the parameter estimates and performance results for fixed effects and variance-covariance for Eqs (13), (14) and (15). Eq (13) had only one variance parameter (σ 2) because its error variance was assumed to be homogeneous. Eqs (14) and (15) had much smaller AIC values and larger Loglik values than Eq (13), which implied that the effects of blocks and the interactions of blocks and sample plots on tree CR were significant. Among these models, Eq (15) had the smallest AIC, -1993, and the largest log-likelihood, 1009. After the parameter estimates were put into Eq (15), the tree CR model for Mongolian oak in northeast China is: Where

Table 6

Parameter estimates and performance assessment results of models.

	Parameter estimates	Eq (13)	Eq (14)	Eq (15)	Eq (17)
Fixed-effect	β^_0	0.6827	0.3650	0.3990	0.6026
parameters	β^_1	-0.0126	-0.0166	-0.0168	-0.0169
	β^_2	-0.0097	-0.0006	0.0007	-0.0046
	β^_3	-0.0413	-0.0325	-0.0361	-0.0386
Canopy density	β^0(1)	—	—	—	-0.1054
effects	β^0(2)	—	—	—	0.0199
	β^0(3)	—	—	—	-0.0662
	β^0(4)	—	—	—	0.0344
Variance	σ^Block02	—	0.1627	0.1643	0.1733
parameters	σ^block12	—	0.0071	0.0070	0.0067
	ρ^_Block01	—	0.324	0.173	0.185
	σ^Block*Plot02	—	0.2251	0.2058	0.2086
	σ^block*Plot12	—	0.0073	0.0042	0.0049
	ρ^_Block*Plot01	—	-0.94	-1	-1
	σ^^2	0.1683	0.1514	0.3622	0.3620
	γ^	—	—	-0.0521	-0.0521
	ρ^	—	—	0.0618	0.0613
Model	AIC	-1568	-1920	-1993	-2063
assessment	Loglik	789	971	1009	1045

Note: AIC, Akaike information criterion; Loglik, log-likelihood; , parameter estimate for autoregressive process of order one.

Model prediction

Table 7 presents the statistics of model performance based on the data used for both model fitting and validation for the base model Eq (13), and two mixed-effect models Eqs (14) and (15) in three cases: PA, block, and block plus interaction effects of blocks and sample plots. In Table 7, Eq (13) itself was regarded as a population average (PA) model and its parameters were directly estimated by fitting to the data without random effects, while Eq (14) PA model only consisting of the components that account for the fixed effects was indirectly obtained by first fitting the model that contained both fixed effects and random effects to the data and then removing the components related to the random effects. This was also applied to Eq (15) PA model. The results showed that although all the models produced the mean prediction errors that were not significantly different from zero at a significance level of 0.05, the prediction accuracy of the models with both the block effects and interaction was the highest, followed by the models with block effects alone, and by the PA models. The major difference between Eqs (14) and (15) was seen in the case of considering both the block effects and interaction. Eq (15) had much smaller statistics of , ξ and δ than Eq (14) although both models slightly under-predicted tree CR and resulted in the mean prediction errors that were not significantly different from zero at a risk level of 0.05.

Table 7

Performance assessment results of the models.

Model No.	Model fitting				Model validation
	e¯	ξ	δ	R2	e¯	ξ	δ
Eq (13)	0	0.1681	0.1681	0.2554	-0.0147	0.1712	0.1718
Eq (14)
PA	-0.0016	0.1686	0.1686	0.2444	-0.0162	0.1702	0.1709
B	0.0044	0.1547	0.1547	0.4579	-0.0095	0.1647	0.1650
B+B*Plot	0.0017	0.0861	0.0861	0.6359	0.0034	0.1143	0.1144
Eq (15)
PA	-0.0012	0.1687	0.1687	0.2439	-0.0159	0.1704	0.1712
B	0.0038	0.1549	0.1549	0.4568	-0.0094	0.1648	0.1650
B+B*Plot	0.0012	0.0814	0.0814	0.6435	0.0025	0.1054	0.1054
Eq (17)	0.0005	0.0412	0.0412	0.6827	0.0012	0.0653	0.0653

Note: : mean prediction error, ξ: variance of prediction error, δ: root mean square error, R2: coefficient of determination between the observed and predicted values, B: block, and B*Plot: interaction of block with plot.

Compared to that from Eq (13), the prediction accuracies of Eqs (14) and (15) with both the block effects and the interaction of blocks and sample plots were much higher. Without the interaction effects, moreover, Eqs (14) and (15) led to very small coefficients of determination (R2 varying from 0.2439 to 0.4579). Adding the interaction effects into both models increased the coefficients to 0.6359 and 0.6435 for Eqs (14) and (15), respectively. Based on the results of model validation, Eq (15) with the interaction effects decreased the root mean square error by 38.65% and 38.33% compared to Eqs (13) and (14) (PA model). These implied that the random effects from the interaction of blocks and sample plots on tree CR were substantially large.

The residuals from Eq (15) with the interaction effects of blocks and sample plots were graphed against the predicted values for the data used for model fitting in Fig 3. Compared to those in Fig 2, the residuals from Eq (15) with the interaction effects, especially the under-predictions, were greatly decreased. This further indicated that the exponential variance model with HD as a predictor accounted for heteroscedasticity effectively. Similar results of the residuals from the Eq (15) were also found for the data used for model validation. Additionally, the error structures of Eq (15) did not show any signs of autocorrelation. Therefore, Eq (15) was promising for predicting tree CR of Mongolian oak.

Fig 3

The residuals of predicted crown ratio values from Eq (15) graphed against the predicted values for Mongolian oak in Wangqing Forest Bureau of northeast China.

Model extension

Our preliminary analyses revealed that the correlations between tree CR and the predictors varied with stand density class. This implied that the accuracy of predicting tree CR could be potentially increased by taking into account the stand density classes. To account for the variation among stand density classes, dummy variables (K 1, K 2, K 3 and K 4) with values of 1 and 0 were introduced into the model. That is, K 1 = 1, K 2 = 1, K 3 = 1, and K 4 = 1 meant the first, second, third and fourth stand density class, respectively, and K 1 = K 2 = K 3 = K 4 = 0 implied the fifth stand density class. Adding the dummy variables into Eq (15) led to following model:

The estimates of model parameters for Eq (17) were listed in Table 6 and its statistics of performance to predict tree CR were shown in Table 7. In both tables, Eq (17) was compared with Eqs (13)–(15) based on the values of the performance statistics. The results showed that compared to Eqs (15) and (17) decreased the value of AIC by 3.51% and increased the value of Loglik by 3.57% (Table 6). Eq (17) also increased the coefficient of determination by 167.31%, 7.36%, and 6.09% compared to those from Eqs (13), (14) and (15), respectively. Based on the results from the validation data, Eq (17) led to a root mean square error of 0.0653, decreasing by 61.99%, 42.92%, and 38.05% compared to those from Eqs (13), (14) and (15). The residuals from Eq (17) were graphed against the predicted values of tree CR for five stand density classes (Fig 4). The results implied that if the stand density of each sample plot was known, Eq (17) in which stand density class was introduced through dummy variables can result in further improvement of predictions.

Fig 4

Residuals of predicted crown ratio values from Eq (17) for each of five stand density classes graphed against the predicted values for Mongolian oak in Wangqing Forest Bureau of northeast China.

Discussion and Conclusions

A nonlinear mixed-effects model was appealing for the analysis of correlated hierarchical structure data because of its flexibility to account for the covariance structures that are not taken into account in traditional regression approaches. A modified logistic model (Eq (15)) with nonlinear mixed-effects model approach at both block and plot levels was recommended for modelling CR of Mongolian oak trees. The results (Tables 6 and 7) showed that a random effects model provided higher accuracy of prediction for the datasets used for both model fitting and validation compared to a similar model without random effects. It is also worth pointing out that the variance-covariance matrices for the random effects and within-plot errors play a more important role in prediction. If we ignore the random effects, within-plot heterogenity and correlation, and apply the ordinary least square approach to the final chosen model (Eq (15)), the statistics (, ξ and δ) for the prediction CR are obviously much larger than those obtained by appropriate variance-covariance structure (see Table 7).

In this study, three tree and stand variables, including D, DD and DH, were selected and involved in the mixed-effect model of tree CR. This finding was supported by other tree CR model studies in which D was often used an independent variable [1, 2, 10, 11]. As Hasenauer and Monserud [5] and Soares and Tomé [1] concluded, the parameter β 1 in our Eq (15) was negative, meaning that increasing D would result in a larger value of tree CR (Table 6). This implied that the effect of D on tree CR was significant. In fact, D is one of the most important tree variables and usually applied to account for stand structure, tree vigor and competition capacity. Soares and Tomé [1] also included DH in their tree CR model for Eucalyptus globules Labill plantations. In most forest growth studies, DH often implies development of trees and stands [47]. On the other hand, DH is a measureable stand characteristic and indicates the site quality in terms of stand growth and yield [48]. In the study of Soares and Tomé [1], greater values for DH resulted in smaller values of CR for Eucalyptus globulus Labill in the north and central coastal regions of Portugal. The finding in this study was consistent with this conclusion, that is, tree CR decreased as DH increased because the parameter β 2 in the model (15) was positive. In addition, the effect of dominant tree diameter DD on CR was also found to be significant in this study and the negative β 3 implied that tree CR increased as DD increased.

Monserud and Sterba [49], Temesgen et al. [10], and Leites et al. [2] found that many other stand or tree variables had significant effects on tree CR, such as tree height (H), and H/D, named as tree slenderness coefficient [3]. In this study, these variables were also examined and it was found that their contributions to the prediction accuracy of the tree CR model were significant. However, measuring H of each tree is prohibitively high cost and time-consuming [15, 47]. Therefore, H and H/D were not selected as predictors in the model (15) to enhance the feasibility and practicality of the model, Site index (SI) affects tree growth significantly ([34, 36, 45, 50], but the effects of Meyer’ SI on CR is insignificant in this study. This is mainly because (i) DH contained in the proposed model could reflect the site quality effectively; (ii) in Wangqing Forest Bureau of northeast China, the environmental and climatic factors greatly vary from sample plot to sample plot due to the complication of environmental and climatic conditions so that the relationship between tree CR and SI greatly differs from place to place; and (iii) introducing SI would complicate the model for a very small portion of gain in predictability. Thus, the stand and tree variables such as SI were not selected.

The tree CR is widely used as a predictor variable to predict growth and yield of trees and forests and also used as a key factor to determine target trees in the management of close-to-nature forests [51-53]. In practice, forest inventory data do not necessarily include the measurements of CR. The mixed-effect CR model (such as Eq (15)) proposed in this study can be used to predict tree CR in those cases. The method that is used to predict CR based on the mixed-effect CR model in practical application consists of two steps: parameter estimation and prediction, which is similar as the approach that proposed in this study. For different tree species or the same tree species but it comes from different area, the values of parameters in the mixed-effect CR model may be different. Therefore, the parameters in the model should be estimated before prediction. In addition, it is also noted that when predicting tree CR using the mixed-effect models such as Eq (15), the random-effect parameters should be estimated from the prior information obtained from the observations of the dependent variable [15, 25, 43, 44, 54]. In addition to the predictors D, DD and DH, tree CR has to be measured from a small sub-sample of trees so that the random-effect parameters can be estimated using EBLUP theory [20, 25, 34]. The appropriate number of trees to be measured for estimating the random-effect parameters has been discussed in the literature [34, 44, 46, 47]. Temesgen et al. [46] and Yang et al. [34] suggested that the more trees used for estimation of random-effect parameters, the higher the prediction accuracy. However, as the number of trees increases, the decrease of prediction error will become insignificant, while the cost of inventory will greatly increase. In fact, the analysis of data using 1–6 trees revealed the accuracy of prediction with four trees was similar to that with five or six trees. As Calama and Montero [44] suggested, therefore, in this study four randomly selected trees were measured for estimation of random-effect parameters.

The prediction of tree CR would have some potential random error. The value of the error varies mainly depending on the prediction accuracy of the mixed-effect CR model (such as (15)). Therefore, the predictions from the forest growth and yield models that use predicted tree CR as a predictor are associated with uncertainty and subject to the random error in the predicted tree CR. One direct method to reduce the uncertainty is using an accurate CR model to predict CR. This is why we recommended the use of the nonlinear mixed-effect modelling approach to develop a CR model for improving the prediction accuracy in this study. In addition, several other approaches, such as simultaneous equations [55], errors-in-variables (EIV) models [56] and Bayesian estimation [57], could also be used to reduce the uncertainty in the case of using the predicted CR as a predictor variable in the forest growth and yield modelling.

Measuring tree CR is subject to error, even though stand and tree attributes are commonly assumed to be measured without error [58, 59]. The measurement errors due to mistakes of field crew and faulty instrument can be substantial [58, 60]. For example, tree CR is generally measured with standard height measurement instrument, but, when the crown is uneven, one often visually rearranges the crown branches to obtain the values of tree CR. There is ample evidence in the literature about the ambiguity of visual estimation of tree characteristics. The studies of Nicholas et al. [61] and Ghosh et al. [62] highlighted the degree of variation that can arise from subjective measurements of tree and stand variables. All the existing tree CR models (Table 1) including the models developed in this study are assumed that (i) tree CR is a random variable; and (ii) the independent variables are fixed and observed without error. It is well known that the violation of the second assumption may lead to biased estimation of regression coefficients and the standard errors of the coefficients and consequently, to misleading test of hypothesis [63]. If the predictors in Eq (15) are considered to have measurement errors, a new approach for the nonlinear mixed-effect models should be developed. However, so far, no algorithms and corresponding computational programs to implement this approach have been available. We are in the process of developing such an appropriate algorithm to deal with this problem, which will certainly be very challenging.

Appendix 1 Parameters and Matrix Structures for Calculation of Random Effects at Both Plot and Block Levels

At both plot and block levels, Is (q 1 + M q 2) × 1 dimensional augmented random effect vector for the i block. is (q 1 + M q 2) × (q 1 + M q 2) dimensional block diagonal positive definite matrix, is design matrix. Z and E (j = 1,…,M ) are n × q 1 and n × q 2 dimensional matrices of the partial derivatives of the nonlinear function f(·) with respect to random effect parameters u and u , respectively. M is the number of plots within the i block. n is the number of observations in the j plot of the i block. q 1 and q 2 are the number of dimensions of u and u , respectively. Ψ and Ψ are the variance-covariance matrixes of u and u , respectively. As mentioned in modeling approach section, the values of in prediction are the same as the estimates obtained in modeling.

Contains the following files: Table A, Observed tree attributes of the fitting and validation data sets. Table B, Results of performance for the used crown ratio models. Table C, Detailed results for the 10 converged mixed-effect models using R/S-Plus nlme function.

(RAR)

Click here for additional data file.