Predicting crown width and length using nonlinear mixed-effects models: a test of competition measures using Chinese fir (Cunninghamia lanceolata (Lamb.) Hook.)

Abstract

Key message

Including individual-tree competition indices as predictor variables could significantly improve the performance of crown width and length models for Chinese fir ( Cunninghamia lanceolata (Lamb.) Hook.). Moreover, distance-dependent competition indices are superior to distance-independent ones when modeling crown width and length. Compared with crown width and length basic models with optimum competition indices, the performance of the two-level nonlinear mixed-effects models improved.

Context

Crown width (CW) and crown length (CL) are two important variables widely included as the predictors in growth and yield models that contribute to forest management strategies.

Aims

Individual-tree crown width and length models were developed with data from 1498 Chinese fir (Cunninghamia lanceolata (Lamb.) Hook.) trees in 16 sample plots located at Jiangle County, Fujian Province, southeastern China. Two hypotheses were proposed: (1) including individual-tree competition indices as predictor variables could significantly improve performance of both the CW—DBH and CL—DBH models; and (2) the distance-dependent competition indices would perform better than distance-independent ones.

Methods

The models were fitted using generalized linear least squares or generalized nonlinear least squares methods. In addition, to prevent correlations between observations from the same sampling unit, we introduced age classes and sample plots as random effects to develop the two-level nonlinear mixed-effects models.

Results

We found introduction of competition indices could significantly improve the performance of the CW—DBH and CL—DBH models. The distance-dependent competition index (i.e., competitor to subject tree distance) performed best in modeling the crown width and length models. Compared with crown width and length basic models with optimum competition indices, the performance of the two-level nonlinear mixed-effects models was significantly better.

Conclusion

The two hypotheses were accepted. We hope these models will contribute to scientific management of Chinese fir plantations.

1 Introduction

Forests are the largest territorial ecosystems on Earth and provide myriad ecosystem services (Buschbacher 1990; Meng et al. 2016b). Generally, ecosystem structure determines function (Warfield 2006), as is the case for forests (Meng et al. 2016a). Canopy structure plays a vital role in shaping forest ecosystems because it serves as a major layer for photosynthesis, respiration, and transpiration (Hernandezmoreno et al. 2017; Wang and Jarvis 1990). Additionally, canopy structure plays a dominant role in regulating energy flow processes, and hence has an important effect on tree growth, dynamics, and productivity (Geißler et al. 2013; Leiterer et al. 2015).

Individual-tree crowns are the basic unit of forest canopy structure (Mallinis et al. 2013). Tree crown is an explicit and direct indicator describing tree vigor and health (Zarnoch et al. 2004), long-term competitiveness (Biging and Dobbertin 1992), productivity (Stenberg et al. 1994), and wood quality (De Kort et al. 1991). It is also extensively used to estimate individual-tree growth (Leites et al. 2009), individual-tree aboveground biomass (Carvalho and Parresol 2003), crown light interception (Pukkala et al. 1991), crown fire risk (Gómez-Vázquez et al. 2014), and stand regeneration (Crookston and Stage 1999). However, labor- and time-intensive measurement of tree crowns impairs its wide application to inform forest management. As such, a number of studies have focused on developing individual crown models (Sharma et al. 2016). Crown models have been widely integrated into individual-tree-based stand simulators, such as Forest Vegetation Simulator (FVS) (Johnson 1997), SILVA (Pretzsch et al. 2006), BWINPro (Nagel and Schmidt 2006), and FORest of RUSsia − Stand (FORRUS − S) (Chumachenko et al. 2003), which have facilitated forest management decision making.

Crown width (CW) and crown length (CL) are two important features of tree crown structure (Zeng 2015). These two features are often predicted using diameter at breast height (DBH) as the independent variable with linear or nonlinear regression (Russell and Weiskittel 2011). The regression methods require the following assumptions (Ritz and Streibig 2008; Sheather 2009): (1) errors e₁, e₂, …, e_n are independent; (2) errors e₁, e₂, …, e_n have a common variance σ²; and (3) the errors are normally distributed with a mean of 0 and variance σ², that is, e|X ~ N(0, σ²). However, forestry data typically are spatially and temporally auto-correlated (Grégoire et al. 1995). The above assumptions could be violated due to the hierarchical and longitudinal structure of forestry data, resulting in biased estimation of standard error of parameters (Calama and Montero 2004; Fu et al. 2013). Because mixed-effects models provide a flexible and powerful tool for analysis of grouped data (e.g., longitudinal data, repeated measures data, blocked designs data, and multilevel data) (Zhao et al. 2013), many authors (Lindstrom and Bates 1990; Vonesh and Chinchilli 1997) have employed it to predict tree crown variables and found that this approach outperformed conventional regression methods (Fu et al. 2013; Sharma et al. 2016). For example, Sharma et al. (2016) and Sánchez-González et al. (2007) introduced sample plot as a random effect into basic CW—DBH functions. Fu et al. (2013) incorporated both site index class and sample plot as random effects into their basic CW—DBH function. Fu et al. (2017) added both block and sample plot as random effects into a basic height to crown base (HCB)—DBH function.

Differing from open-grown trees, the crown width and length for stand-grown trees is significantly influenced by competition among individual trees (Pacala et al. 1996; Sharma et al. 2016). Individual-tree competition indices (CIs), which can be classified into distance-dependent CIs and distance-independent CIs, have been being developed since the 1960s to quantify individual-tree-level competition (Schröder and Gadow 1999; Stadt et al. 2002). Unfortunately, only a few studies on crown width and length models that include distance-dependent CIs (Davies and Pommerening 2008; Purves et al. 2007; Rüdiger 2003; Rouvinen and Kuuluvainen 1997; Thorpe et al. 2010), and fewer researchers compared the performance of distance-dependent CIs and distance-independent CIs in developing crown models (Sharma et al. 2016).

In the present study, we hypothesized that including CIs as predictor variables could significantly improve performance of both CW—DBH and CL—DBH models. Moreover, since spatial arrangement of trees could have a significant influence on their growth (Pukkala and Kolstroem 1991), we proposed a second hypothesis that distance-dependent CIs perform better than distance-independent CIs. The objectives of this present study were (1) to test the hypotheses by including both distance-dependent and distance-independent CIs for development of CW and CL models, (2) identify the optimum CIs for the models, and (3) produce mixed-effects predictive models for CW and CL using optimum CIs for a Chinese fir plantation.

2 Materials and methods

2.1 Study area

Our study area was located in Jiangle County (26° 26′ N to 27° 04′ N, 117° 05′ E to 117° 40′ E) in Fujian Province, southeastern China. It has a subtropical monsoon climate with characteristics of both an oceanic and continental climate. The annual average temperature is 18.7 °C, the annual average precipitation is 1698.2 mm, the annual average sunshine time is 1730 h, and the average frost-free period is 295 days (Lin et al. 2018). The main soil type in the study area is red soil. Hills and low mountains are the main landforms, accounting for 90% of the region. The average elevation is 540 m, with the highest peak 1620 m. The main vegetation types include natural secondary forest, Chinese fir plantations, and Masson pine (Pinus massoniana Lamb.) plantations (Chen et al. 2018).

2.2 Data source

Data employed in this study were collected from 16 squared sample plots of 400 m² in even-aged, pure Chinese fir plantations. C. lanceolata accounted for 97.73% of total basal area in this plantation, with the other 2.27% represented by 15 tree species, e.g., Cinnamomum camphora (L.) Presl (0.73%), Vernicia montana Lour. (0.36%), and Magnolia officinalis Rehd. et Wils. (0.20%). These 16 sample plots were further grouped into the following categories based on age classes: 4 sample plots of young stands (6 years), 4 sample plots of middle-aged stands (16 years), 4 sample plots of near mature stands (23 years), and 4 sample plots of mature stands (30 years).

For each standing live tree, DBH and relative coordinates were measured with a precision of 0.1 cm and 0.01 m, respectively. Using a Vertex III altimeter (Vasilescu 2013), total tree heights (H) and height to crown base (HCB) were measured with a precision of 0.1 m. CW was calculated as the arithmetic mean of two crown widths, derived from the four measured crown radii at two azimuths. The first azimuth was defined as the direction from the subject tree to the center of the sample plot, and the second azimuth was perpendicular to the first (Marshall et al. 2003; Sharma et al. 2016). The number of trees measured was 1565, including 1498 Chinese fir and 67 broad-leaved trees. These 1565 trees were used to derive competition indices for each individual tree (but CW and CL models were fitted using the 1498 Chinese fir trees). The descriptive statistics for the 1498 Chinese fir trees are summarized in Table 1.

Table 1 Descriptive statistics of Chinese fir trees (n = 1498)

2.3 Methods

2.3.1 Basic model selection

The ten functions provided in Table 2 were selected as candidate basic models for the CW—DBH and CL—DBH relationships for individual trees (Sánchez-González et al. 2007; Sönmez 2009).

Table 2 Crown width or length basic models considered in the present study

The candidate basic models were fitted using generalized linear least squares (GLS) or generalized nonlinear least squares (GNLS) method with gls or gnls function in the nlme package of R software (Pinheiro et al. 2012; R Team RDC 2013). For model evaluation, the lack-of-fit test statistics, i.e., − 2 log-likelihood (− 2LL), Akaike information criterion (AIC), Bayesian information criterion (BIC), absolute bias (Bias), root mean square error (RMSE), and coefficient of determination (R²), were calculated.

Additionally, the predictive performance of models was further validated using the ten-fold cross-validation procedure by calculating normalized mean square error for the test set (NMSE) (Levi et al. 2015) and predicted the error sum of squares (PRESS) (Quan 1988). The formulas were as follows:

$$\mathrm{NMSE}=\overset{}{\;\overline{{(y_j-\widehat{y_j})}^2}}/\overset{}{\overline{{(y_j-\overline{y_j})}^2}}=\frac{\sum{(y_j-\widehat{y_j})}^2}{\sum{(y_j-\overline{y_j})}^2}$$

(1)

$$\mathrm{PRESS}=\sum_{j=1}^{n}({{y}_{j}-\widehat{{y}_{j}})}^{2}$$

(2)

where y_j is the jth observed value in the test set, $\widehat{{y}_{j}}$ is the jth value predicted from y_j by the model fitted by the training set, and $\overline{{y }_{j}}$ is the mean of the observed values in the test set.

Based on model evaluation and validation, the final optimum basic model was determined.

2.3.2 Competition indices calculation

Ten candidate competition indices including four distance-independent indices and six distance-dependent indices (Table 3) were compared for their performance in predicting crown width and length.

Table 3 List of candidate competition indices

To calculate the competition indices, we classified trees in the sample plots into subject trees and competitor trees. Each and every tree needs be treated as both types for given data quality and amounts. Identifying competitor trees are of great importance to derive competition indices. However, competitor searching methods vary among different studies (Bella 1971; Pommerening 2008; Wang et al. 2016). Among these methods, four neighboring trees proposed by Hui and Gadow (Hui and Gadow 2003) (four neighboring trees as competitors in four different directions around the subject tree) and fixed radii (all neighboring trees as competitors around the subject trees lying around a search radius with 3.5 times the mean crown radius of canopy trees) (Mailly et al. 2003) have been extensively employed. In the data pre-analysis, we applied these two methods and found that fixed radii worked better, which is also consistent with Lorimer (1983) who recommends such a rule for most applications. In the present study, we selected the second method to search competitors.

Subject trees near sample plot edges may have competitor trees that are out of the sample plot. These are edge effects, which could bring errors into estimation of competition indices (Tang et al. 2007). We employed toroidal edge correction (also referred to as translation) to reduce such edge effects. For each sample plot, make 8 copies of the original sample plot, then one of the copies is moved from original position into up, down, left, right, upper left, lower left, lower right, and upper right direction, respectively. So a new bigger sample plot consisted of 9 original is created (Fig. 1). When calculating individual competition indices, subject trees are defined as individuals located in original sample plot (Haase 1995; Pommerening and Stoyan 2006; Diggle 2013). Scatters distribution of crown width or length and calculated CIs showed that competitions have a significant effect on the growth of crown width and length (Figs. 2, 3).

Fig. 1

Illustration of the translation edge-correction methods used in this study

Fig. 2

Scattered plots of crown width against different competition indices, respectively

Fig. 3

Scattered plots of crown length against different competition indices, respectively

2.3.3 Optimum competition index identification

To determine the optimum competition index and its optimum position, each competition index was added independently with different positions to the basic model determined in Basic model selection section. Following Schröder and Gadow (1999) and Mailly et al. (2003), we calculated the mean square error reduction (MSER) relative to a no competition index to identify the optimum combination of the competition index and its position for the crown width or length models. The formula of MSER was as follows:

$$\mathrm{MSER}=\left(1-\frac{{\mathrm{MSE}}_{2}}{{\mathrm{MSE}}_{1}}\right)\times 100$$

(3)

where ${\mathrm{MSE}}_{1}$ is the mean square error of the basic crown width or crown length model and MSE₂ is the mean square error of the models after adding the competition index.

For the CW—DBH and CL—DBH relationships, the optimum basic models determined in basic model selection section, and the optimum combination of competition indices and their positions identified in this section, were then employed as our final basic model, based on which the mixed-effects crown width or length models were further developed.

2.3.4 Two-level nonlinear mixed-effects models

General model form

Due to hierarchical structure of data (i.e., trees within a sample plot, and sample plots within an age class), we introduced two-level random effects to our basic model. The first-level random effect is age class, and sample plot (nested in age classes) serves as the second-level random effect.

The formulas of two-level linear or nonlinear mixed-effects models can be expressed as follows:

$$\begin{array}{c}{y}_{ijk}=f\left({\varphi }_{ijk},{t}_{ijk}\right)+{\varepsilon }_{ijk},i=1,\dots ,M, j=1,\dots ,{M}_{i};k=1,\dots ,{n}_{ij}\\ {\varepsilon }_{ijk}\sim N(0,{\sigma }^{2}{R}_{ij})\end{array}$$

(4)

where ${y}_{ijk}$ denotes a response value of the kth observation (tree) on the jth group (sample plot) nested within-group i (age class). M is the number of first-level groups, M_i is the number of second-level groups within-group i, n_ij is the number of observations on the jth group in group i; f(·) is a general, real-valued, differentiable function of a group-specific parameter vector φ_ijk and a covariate vector t_ijk. The within-group error ε_ijk, which contains within-group variance and correlation, is assumed to be normally distributed with zero expectation and a positive-definite variance–covariance structure R_ij, generally expressed as a function of the parameter vector λ (Fu et al. 2013).

The φ_ijk can be further written as follows:

$$\begin{array}{cc}\varphi_{ijk}=A_{ijk}\beta+B_{i,jk}\mu_i+M_{ijk}\mu_{ij},&\mu_i\sim N\left(0,\psi_1\right),\mu_{ij}\sim N(0,\psi_2)\end{array}$$

(5)

where β is a p-dimensional vector of fixed effects. ${\mu }_{i}$ and ${\mu }_{ij}$ are the first and second random effects, which are independent and normally distributed q₁- and q₂-dimensional vectors with zero mean and respective variance–covariance matrices ψ₁ and ψ₂. A_ijk, B_i,jk, and M_ijk are design matrices. ${\mu }_{i}$, ${\mu }_{ij}$, and ${\varepsilon }_{ijk}$ are mutually independent. The variance–covariance matrices of ψ₁ and ψ₂ were assumed to be unstructured.

Determining parameter effects

In this study, all possible combinations of random effects were fitted and then the best combination was selected according to lack-of-fit statistics, i.e., − 2LL, AIC, BIC, and likelihood ratio test (L.Ratio).

Determining the matrix of ${R}_{ij}$

Forestry data often exhibit autocorrelation and heteroskedasticity (Gregorie 1987), which results in incorrect standard errors and estimation intervals. Therefore, within-group variance and autocorrelation structure in ${R}_{ij}$ should be specified (Davidian 2017; Meng and Huang 2009). The following within-group variance–covariance matrix (Davidian 2017; Zhao et al. 2013) is usually used:

$${R}_{ij}={\sigma }^{2}{G}_{ij}^{0.5}{\Gamma }_{ij}{G}_{ij}^{0.5}$$

(6)

where ${\sigma }^{2}$ is a scaling factor that is equal to residual variance of the developed model, ${G}_{ij}$ is a diagonal matrix which describes heteroscedasticity, and ${\Gamma }_{ij}$ is a matrix showing the autocorrelation structure of errors.

Because our data are only from one measurement period, autocorrelation was not considered. Therefore, ${\Gamma }_{ij}$ reduced to a n_ij × n_ij identity matrix, and three variance functions, i.e., the exponential function [Eq. $($7)], the power function [Eq. (8)] and the constant plus power function [Eq. (9)], were used to reduce heteroscedasticity (Wang et al. 2019; Zhao et al. 2013).

$$\mathrm{varExp}({\varepsilon }_{ijk})={\sigma }^{2}\mathrm{exp}({2\alpha }\mathrm{u}_{ijk})$$

(7)

$$\mathrm{varPower}({\varepsilon }_{ijk})={\sigma }^{2}\mathrm{exp}({u}_{ijk}^{2\alpha })$$

(8)

$$\mathrm{varConstPower}({\varepsilon }_{ijk})={\sigma }^{2}{(\alpha +{u}_{ijk}^{\beta })}^{2}$$

(9)

where ${u}_{ijk}$ is the estimated value based on fixed parameters of the mixed-effects models and α and β are estimated parameters of variance functions.

Model prediction and evaluation

The final produced mixed-effects models for crown width and length with optimum competition indices were evaluated for their predictive performance. We performed a tenfold cross-validation procedure and NMSE and PRESS were calculated. In addition, likelihood ratio tests were conducted among the optimum basic model, optimum basic model with the best CI, and the optimum mixed-effects models for the CW and CL.

3 Results

3.1 Determination of basic models

The lack-of-fit statistics of the candidate CW and CL models are summarized in Tables 4 and 5, respectively. Although the lack-of-fit statistics were almost identical for most of the functions, the function [CR4] (cubic form) showed a slightly better performance for both the crown width and length. NMSE and PRESS produced by the tenfold cross-validation also indicated slightly better performance for the function [CR4]. Moreover, the scattered plots of the CW or CL against DBH (Fig. 4) were suitfully fitted by the function [CR4]. Therefore, this function was selected as the basic model for modeling the CW—DBH and CL—DBH relationships. The optimum basic models are as follows:

$${\mathrm{CW}}_{ijk}=0.5403{{\mathrm{DBH}}_{ijk}}^{0.5794}+{\varepsilon }_{ijk}$$

(10)

$${\mathrm{CL}}_{ijk}=0.4356{{\mathrm{DBH}}_{ijk}}^{0.8829}+{\varepsilon }_{ijk}$$

(11)

where CW_ijk, CL_ijk, and DBH_ijk are crown width, crown length, and diameter at breast height, respectively, of the kth tree in the jth plot in the ith age class; ε_ijk is a model error term.

Table 4 Lack-of-fit statistics for CW—DBH basic models

Table 5 Lack-of-fit statistics for CL—DBH basic models

Fig. 4

Scattered plots of crown width or length against DBH for [CR4], respectively

3.2 Inclusion of competition index

Based on Eq. (10) and Eq. (11), effects of adding different competition indices at varying positions on the improvement of CW—DBH and CL—DBH models were compared (Tables 6 and 7).

Table 6 Comparison of different competition indices added to the CW—DBH basic model [CW4]

Table 7 Comparison of different competition indices added to CL—DBH basic model [CL4]

In terms of the CW—DBH model, it can be seen from Table 6 that JiQi has the largest MSER of 32.38, but it was not selected since its insignificant coefficient. Therefore, the Alem was identified as the optimum competition index since it yielded the second largest MSER of 20.44. The corresponding model form was written as:

$${\mathrm{CW}}_{ijk}=(0.9533{+0.0516\mathrm{Alem})DBH}_{ijk}^{0.2474}+{\varepsilon }_{ijk}$$

(12)

Regarding crown length, the OP exhibited a better performance than the other competition indices with the largest MSER of 14.31 (Table 7). The corresponding model form was as follows:

$${\mathrm{CL}}_{ijk}=0.5991{{\mathrm{DBH}}_{ijk}}^{\left(0.5122+0.2897{\mathrm{OP}}_{ijk}\right)}+{\varepsilon }_{ijk}$$

(13)

where JiQi_ijk and OP_ijk are JiQi and OP, respectively, of the kth tree in the jth sample plot in the ith age class.

In addition, the CW and CL curves overlaid on the measured data (DBH and the optimum competition index, i.e., Alem and OP) were presented in Fig. 5, which showed that the scattered plots were suitfully fitted by the Eqs. (12–13).

Fig. 5

Scattered plots of crown width or length against DBH and CI for Eqs. (12–13)

3.3 Development of the two-level mixed-effects model

There are 49 (${C}_{3}^{1}{C}_{3}^{1}$ + ${C}_{2}^{1}{C}_{3}^{2}{C}_{3}^{1}$ +${C}_{2}^{1}{C}_{3}^{3}{C}_{3}^{1}$+ ${C}_{3}^{2}{C}_{3}^{2}$ + ${C}_{2}^{1}{C}_{3}^{2}{C}_{3}^{3}$ + ${C}_{3}^{3}{C}_{3}^{3}$  = 49) potential combinations of random effects for Eq. (12) when considering all independent variables and the intercept in the basic model. When fitted to the data, all combinations reached convergence. Among these 49 mixed-effects models, Eq. (14) yielded the smallest AIC (1571), BIC (1613) and − 2LL [1554.91 (8)]. In addition, the L.Ratio (L.Ratio = 105.72, p < 0.0001) also indicated the best performance of Eq. (14). Therefore, Eq. (14) was the resulting mixed-effects model.

$${\mathrm{CW}}_{ijk}=[({\phi }_{1}+{u}_{1ij})+{\phi }_{3}{\mathrm{Alem}}_{ijk}]{DBH}_{ijk}^{{\phi }_{2}+{u}_{1i}+{u}_{2ij}}+{\varepsilon }_{ijk}$$

(14)

where ${\phi }_{1}$–${\phi }_{3}$ are the fixed-effects parameters, u_1i is a random-effects parameter generated by age class for the CW model, and u_1ij and u_2ij are random-effects parameters generated by interaction of age class and sample plot for the CW model.

There are 49 (${C}_{3}^{1}{C}_{3}^{1}$ + ${C}_{2}^{1}{C}_{3}^{2}{C}_{3}^{1}$ +${C}_{2}^{1}{C}_{3}^{3}{C}_{3}^{1}$+ ${C}_{3}^{2}{C}_{3}^{2}$ + ${C}_{2}^{1}{C}_{3}^{2}{C}_{3}^{3}$ + ${C}_{3}^{3}{C}_{3}^{3}$  = 49) potential combinations of random effects for Eq. (13) when considering all independent variables and the intercept in the basic model. When fitted to the data, a total of 38 combinations reached convergence. Among these 38 mixed-effects models, Eq. (15) yielded the smallest AIC (4869), BIC (4912), and − 2LL [4853.39 (8)]. In addition, the L.Ratio (L.Ratio = 9.05, p = 0.0108 < 0.05) also indicated the best performance of Eq. (15). Therefore, Eq. (15) was the resulting mixed-effects model.

$${\mathrm{CL}}_{ijk}=({\phi }_{1}+{v}_{1ij}){\mathrm{DBH}}^{[{\phi }_{2}+{v}_{1i}+\left({\phi }_{3}+{v}_{2ij}\right)\mathrm{OP}]}+{\varepsilon }_{ijk}$$

(15)

where ${\phi }_{1}$–${\phi }_{3}$ are fixed-effects parameters, ${v}_{1i}$ is a random-effects parameter generated by age class for the CL model, and ${v}_{1ij}$ and ${v}_{2ij}$ are random-effects parameters generated by interaction of age class and sample plot for the CL model.

We used three variance functions to reduce the heteroscedasticity of the residuals and their performances were compared in terms of AIC, BIC, − 2LL, and L.Ratio. The results showed that the mixed-effects models with variance function exhibited better performance than the one without considering the variance function. Performance also differed among the mixed-effects models with different variance functions (Tables 8 and 9). According to AIC, BIC, − 2LL, and L.Ratio, the exponent function was determined as the best variance function for CW and CL models.

Table 8 The lack-of-fit statistics of the crown width mixed-effects model using different variance functions

Table 9 The lack-of-fit statistics of the crown length mixed-effects model using different variance functions

After the determination of parameter effects and error variance–covariance structure, the final forms of the CW and CL nonlinear mixed-effects models for Chinese fir in Fujian Province, southeastern China, were proposed in Eq. (16) and Eq. (17). The residual plots and QQ plots are shown in Figs. 6 and 7, respectively.

$${\mathrm{CW}}_{ijk}=[(0.4229+{u}_{1ij})-(8.730e-4{)\mathrm{Alem}}_{ijk}]{DBH}_{ijk}^{0.7212+{u}_{1i}+{u}_{2ij}}+{\varepsilon }_{ijk}$$

(16)

where

$$\begin{array}{l}{u}_{i}=\left[{u}_{1i}\right]\sim \mathrm{N}\left\{\left[0\right],{\psi }_{i}=\left(0.0014\right)\right\}\\ {u}_{ij}=\left[\begin{array}{c}{u}_{1ij}\\ {u}_{2ij}\end{array}\right]\sim \mathrm{N}\left\{\left[\begin{array}{c}0\\ 0\end{array}\right],{\psi }_{i}=\left(\begin{array}{cc}0.0704& -0.0507\\ -0.0507& 0.0415\end{array}\right)\right\}\\ {\varepsilon }_{ij}\sim \mathrm{N}(0, {R}_{ij}=0.1052{G}_{ij}^{0.5}{\Gamma }_{ij}{G}_{ij}^{0.5}), {G}_{ij} =\mathrm{ exp }(0.0732{y}_{i}), {\Gamma }_{ij} ={I}_{ni}\end{array}$$

$${\mathrm{CL}}_{ijk}=(1.3129+{v}_{1ij}){\mathrm{DBH}}^{[0.1669+{v}_{1i}+\left(0.3562+{v}_{2ij}\right)\mathrm{OP}]}+{\varepsilon }_{ijk}$$

(17)

where

Fig. 6

Residual plot and QQ plot of the selected crown width mixed-effects model [Eq. (16)]

Fig. 7

Residual plot and QQ plot of the selected crown length mixed-effects model [Eq. (17)]

$$\begin{array}{l}{v}_{i}=\left[{v}_{1i}\right]\sim \mathrm{N}\left\{\left[0\right],{\psi }_{i}=\left(0.0029\right)\right\}\\ {v}_{ij}=\left[\begin{array}{c}{v}_{1ij}\\ {v}_{2ij}\end{array}\right]\sim \mathrm{N}\left\{\left[\begin{array}{c}0\\ 0\end{array}\right],{\psi }_{i}=\left(\begin{array}{cc}0.0555& -0.0164\\ -0.0620& 0.0057\end{array}\right)\right\}\\ {\varepsilon }_{ijk}\sim \mathrm{N}(0, {R}_{ij}=0.6124{G}_{ij}^{0.5}{\Gamma }_{ij}{G}_{ij}^{0.5}), {G}_{ij} =\mathrm{ exp }(0.0958{y}_{i}), {\Gamma }_{ij} ={I}_{ni}\end{array}$$

3.4 Model prediction and evaluation

The results of L.Ratio (Table 10) indicated that in comparison to Eq. (10) and Eq. (12), the addition of random parameters in Eq. (16) could significantly improve the predictive ability for crown width. A similar improvement for crown length can also be found in Eq. (17) relative to Eq. (11) and Eq. (13). In addition, these improvements were reflected by NMSE and PRESS of the ten-fold cross-validation. Figure 8 also showed that the final models were unbiased, and model fits were good.

Table 10 Comparisons among optimum basic model, optimum basic model with optimum competition index and optimum mixed-effects model for crown width and length, respectively

Fig. 8

Plots of predicted values against observed values of the final models

4 Discussion

We found that including CIs as predictor variables could significantly improve the performance of CW—DBH and CL—DBH basic models (Tables 6 and 7). The first hypothesis proposed in the Introduction was therefore accepted. Similar results were reported by Sharma et al. (2016), Thorpe et al. (2010), and Davies and Pommerening (2008), all who suggested that competition could significantly affect crown dimension and should be included into crown models. The CIs are comprehensive indicators describing the competition status of individual trees and hence can well represent an individual tree’s space occupation and allocation of resources (light, water, or nutrients), which in turn determines growth (Kahriman et al. 2018). In addition to CW and CL models, CIs have been extensively used to develop individual-tree growth models. Many authors (Brand 2011; Collet and Chenost 2006) have found individual-tree growth models including CIs could provide for more accurate estimation.

The distance-dependent competition index, Alem, performed best in modeling the CW—DBH relationship (Table 6). Similarly, the distance-dependent competition index, OP, showed best performance in the CL—DBH model (Table 7). Our second hypothesis, i.e., the distance-dependent CIs perform better than distance-independent CIs for CW—DBH and CL—DBH relationship, was also accepted, which suggested that the distance-dependent CIs more appropriately and adequately describe competitive situations among individual trees in a stand than the distance-independent ones (Contreras et al. 2011; Sharma et al. 2016). Therefore, the interactions of the trees in a spatial manner over the restricted distances is essential for growth dynamics of the tree crowns. Many authors have also documented that distance-dependent CIs are superior to distance-independent ones, especially in mixed-species stands (Contreras et al. 2011; Quinonez-Barraza et al. 2018; Sharma et al. 2016). Additionally, the model showed that CW significantly increases with increasing Alem (Fig. 9). Alem, proposed by Alemdag (1978), is an index defined on the basis of the growing space for the subject tree (Burkhart and Tomé 2012). Therefore, trees with higher Alem or trees with larger growing space and area potentially available is expected to have larger crown sizes. OP was significantly positively related to CL (Fig. 9). Since OP reflects the degree to which the subject tree in the spatial structural unit is not shaded by competitor trees (Hu 2020). OP is an indicator of individual-tree light environment status in a forest stand and its relationship with CL suggested that trees with larger OP have stronger competitive ability for light which is considered as the major limiting resource for individual-tree growth. Therefore, a high OP will increase light supply, light capture, and light use efficiency, resulting in a larger crown.

Fig. 9

Effect of competition indices on the crown width or length respectively for Chinese fir. The curves were produced using parameter estimates of optimum basic models with optimum competition index [Eqs. (12–13)]

A disadvantage of distance-dependent competition index is that it requires tree attributes and tree locations, which are expensive and labor intensive to acquire. However, with the development of remote sensing and geographic information systems (such as Light Detection and LiDAR) have facilitated the acquisition of tree-level spatial and dimensional data for the entire stands (Packalen and Maltamo 2006; Rowell 2009; Suratno et al. 2009), and are increasingly being used for forest and natural resource applications (Contreras et al. 2011; Suratno et al. 2009). The crown width and length models that depend on DBH and distance-dependent CIs in this study can be easily integrated with LiDAR inventory data and will be a useful tool to evaluate the influences of alternative management measures over time.

Except DBH and the competition index, which affect the crown width and length significantly, the crown structure is also largely influenced by site index, which is the most widely used measures of potential forest productivity and has strongly piratical utilized as a key input to both forest processes-based models and empirical growth and yield models (Amponsah et al. 2004; Fish et al. 2011; Yang 1998). Therefore, in this study, the SI was calculated respectively for the 16 sample plots by SI model (base age 20) developed by Duan and Zhang (2004) for the Chinese fir plantation to estimate the site quality. Unfortunately, because the results were not significant, and the site index was not left in the model. It might be attributed to the small study area, resulting in little difference in site quality. Additionally, many authors included other stand- or individual tree-level variables to model crown structure. For instance, Fu et al. (2017) evaluated 18 variables, including stand age [A (years)], height [H (m)], height to crown base [HCB (m)], site index [SI (m)], plot dominant height [DH (m)], and total diameter [LDTD (cm)]; four variables, i.e., DBH, DH, HCB, and H, were left the final model. For our model, we did not include other variables because more independent variables would reduce model generality (Hasenauer 2006; Vanclay 1994). Furthermore, additional independent variables also entail more labor and time, making data collection more expensive. Additionally, more independent variables impact model convergence and computational speed of parameter estimation (Montgomery et al. 1982).

Generality, precision, and reality are three important properties of a model that researchers strive to maximize (Levins 1966). However, Levins (1966) argued one of the three properties has to be sacrificed in order to achieve a higher level of the other two. Normally, many models have sacrificed generality to enhance reality and precision (Burkhart and Tomé 2012). Burkhart and Tomé (2012) documented that there is an increased interest in improving generality of models in the context of rapidly changing management and environmental conditions. In our study, we only included two variables to model CW and CL to enhance generality.

We produced two-level nonlinear mixed-effects models for predicting CW and CL. The performance of models significantly improved after introducing random effects. For instance, in comparison to the optimum fixed-effects CW model [Eq. (12)] and CL model [Eq. (13)], the AIC of the optimum mixed-effects models for CW and CL dropped by 653 and 330 respectively, and the adjusted R² increased by 0.1952 and 0.0577. Moreover, the model prediction and evaluation using the tenfold cross-validation further supported the conclusion that the mixed-effects approach had significantly improved the models’ predictive performance. NMSE and PRESS were significantly reduced. The advantage of using a mixed-effects modeling approach to model crown structure has been also documented by other authors (Fu et al. 2013, 2017; Sharma et al. 2016). For instance, using a nested, two-level nonlinear mixed-effects approach, Fu et al. (2013) improved the prediction accuracy of crown width by including a site index and sample plot as random effects.

Mixed-effects models are tools used to deal with repeated measures and spatially-correlated data (Pinheiro and Bates 2006). Many authors documented that introducing random effects could correct or reduce autocorrelation and heteroscedasticity (Hasenauer 2006; Montgomery et al. 1982; Pommerening and Stoyan 2006; Wang et al. 2016). A similar result was found in our study. In addition, we further introduced three variance functions to refine our model. Finally, based on the lack-of-fit statistics, the exponent function was determined as the optimum option for reducing the heteroscedasticity of the residuals. Similar results were also reported by Zhao et al. (2013), Calama and Montero (2005) and Wang et al. (2019).

Determining parameters is of great importance when developing a mixed-effects model (Calama and Montero 2005). In previous studies, a common approach has been to first determine parameter effects and then incorporate additional predictor variables as fixed effects into a model (Budhathoki et al. 2008; Calama and Montero 2005). The additional predictor variables could introduce random effects (Fu et al. 2017), which often was not accounted for in previous studies. In this study, we determined parameter effects after including CIs in the models, and the random effects that the CIs introduced could be therefore be well represented.

The traditional method of computing CIs, we applied in this study, migth result in significant bias on the predicted variable of interest. Compared to the traditional method, more modern and more appropriate methods of computing CIs simultaneously with the model parameters have been used in a few studies, which can be able to describe competition impact on growth and development of individual-tree characteristics and produce more accurate output (Pommerening et al. 2011; Sharma and Brunner 2017). Therefore, the improving methods of computing CIs are strongly recommended in the future. Additionally, climate change has been reported to affect crown structure (Carnicer et al. 2011; Solberg 2004). For example, Carnicer et al. (2011) documented that drought resulting from climate change has a long-lasting chronic effect on crown structure, including increased crown defoliation rates for all tree species examined. It is imperative to integrate climatic variables into crown width and length models when data is available.

5 Conclusion

We produced CW and CL models using a two-level nonlinear mixed-effects approach. Our models showed that including CIs as predictor variables could significantly improve the performance of the models. Similarly, the second hypothesis, that distance-dependent CIs are superior to distance-independent ones when modeling CW and CL, was also supported. Compared with the CW and CL models with no random effects, the performance of the two-level nonlinear mixed-effects models improved. It is hoped that the final CW and CL models will contribute to the scientific management of the Chinese fir in other sites and under different environmental conditions.