Development of an Individual-Tree Basal Area Increment Model using a Linear Mixed-Effects Approach

Summary

Mixed-effects models are flexible and useful tools for analyzing data with a hierarchical stochastic structure in forestry and could also be used to significantly improve the performance of forest growth models. Here, a protocol is presented that synthesizes information relating to linear mixed-effects models.

Abstract

Here, we developed an individual-tree model of 5-year basal area increments based on a dataset including 21898 Picea asperata trees from 779 sample plots located in Xinjiang Province, northwest China. To prevent high correlations among observations from the same sampling unit, we developed the model using a linear mixed-effects approach with random plot effect to account for stochastic variability. Various tree- and stand-level variables, such as indices for tree size, competition, and site condition, were included as fixed effects to explain the residual variability. In addition, heteroscedasticity and autocorrelation were described by introducing variance functions and autocorrelation structures. The optimal linear mixed-effects model was determined by several fit statistics: Akaike’s information criterion, Bayesian information criterion, logarithm likelihood, and a likelihood ratio test. The results indicated that significant variables of individual-tree basal area increment were the inverse transformation of diameter at breast height, the basal area of trees larger than the subject tree, the number of trees per hectare, and elevation. Furthermore, errors in variance structure were most successfully modeled by the exponential function, and the autocorrelation was significantly corrected by first-order autoregressive structure (AR(1)). The performance of the linear mixed-effects model was significantly improved relative to the model using ordinary least squares regression.

Introduction

Compared with even-aged monoculture, uneven-aged mixed-species forest management with multiple objectives has received increased attention recently^1,2,3. Prediction of different management alternatives is necessary for formulating robust forest management strategies, especially for complex uneven-aged mixed-species forest⁴. Forest growth and yield models have been used extensively to forecast tree or stand development and harvest under various management schemes^5,6,7. Forest growth and yield models are classified into individual-tree models, size-class models, and whole-stand growth models^6,7,8. Unfortunately, size-class models and whole-stand models are not appropriate for uneven-aged mixed-species forests, which require a more detailed description to support the forest management decision-making process. For this reason, individual-tree growth and yield models have received increased attention throughout the last few decades because of their ability to make predictions for forest stands with a variety of species compositions, structures, and management strategies^9,10,11.

Ordinary least squares (OLS) regression is the most commonly used method for the development of individual-tree growth models^12,13,14,15. The datasets for individual-tree growth models collected repeatedly over a fixed length of time on the same sampling unit (i.e., sample plot or tree) have a hierarchical stochastic structure, with a lack of independence and high spatial and temporal correlation among observations^10,16. The hierarchical stochastic structure violates the fundamental assumptions of OLS regression: namely independent residuals and normally distributed data with equal variances. Therefore, the use of OLS regression inevitably produces biased estimates of the standard error of parameter estimates for these data^13,14.

Mixed-effects models provide a powerful tool for analyzing data with complex structures, such as repeated measures data, longitudinal data, and multi-level data. Mixed-effects models consist of both fixed components, common to the complete population, and random components, which is specific to each sampling level. In addition, mixed-effects models take into account heteroscedasticity and autocorrelation in space and time by defining non-diagonal variance-covariance structure matrices^17,18,19. For this reason, mixed-effects models have been extensively used in forestry, such as in diameter-height models^20,21, crown models^22,23, self-thinning models^24,25, and growth models^26,27.

Here, the main objective was to develop an individual-tree basal area increment model using a linear mixed-effects approach. We hope that the mixed-effects approach could be broadly applied.

Protocol

1. Data preparation

1. Prepare modeling data, which includes individual-tree information (species and diameter at breast height at 1.3 m) and plot information (slope, aspect, and elevation). In this study, the data were obtained from the 8th (2009) and 9th (2014) Chinese National Forest Inventory in Xinjiang Province, northwest China, which includes 21,898 observations of 779 sample plots. These sample plots are square-shaped with a size of 1 Mu (Chinese unit of area equivalent to 0.067 ha) and are systematically arranged over a grid of 4 km x 8 km.
   NOTE: Data for modeling (basal area) increment requires at least one growth period (i.e., two observations).
2. Randomly divide the data into two datasets, with 80% of the data from the sample plots used for model fitting (model development dataset), which consists of 17,145 observations from 623 sample plots and 20% for model validation (model validation dataset) which consists of 4,753 observations from 156 sample plots. Descriptive statistics for the key variables used are provided in Table 1.
   NOTE: This step of the modeling procedure can be omitted, and all data is used for model development.

Variables	Fitting data					Validation data
	Min	Max	Mean	S.D.		Min	Max	Mean	S.D.
DBH1 (cm)	5	124.8	19.9	13.2		5	101.5	19.5	13.4
QMD (cm)	6.7	82.3	22.5	8.5		9.2	73.3	21.8	9.2
ID (cm)	0.1	14.4	1.1	1		0.1	16.9	1	1.1
BAL (m3)	0	5.2	1.7	0.9		0	5.4	1.7	1
NT (trees/ha)	14.9	3642	1072	673.7		14.9	3418	1205	829.3
BA (m2/ha)	0.1	77.5	34.2	13.9		0.1	80.6	34.5	15.3
EL (m)	2	3302	2189	340.3		1441	3380	2256	308.3

Table 1. Summary statistics for fitting and validation data. DBH₁: initial diameter at breast height at 1.3 m (DBH), DBH₂: DBH measured after 5 years of growth, QMD: quadratic mean diameter, ID: diameter increment for 5 years (DBH₂ – DBH₁), BAL: the basal area of trees larger than the subject tree (the subject tree: the tree which was calculated the competition indices), NT: the number of trees per hectare, BA: basal area per hectare, EL: elevation, S.D.: standard deviation.

2. Basic model development

1. Consult references to identify variables that affect individual-tree basal area increments.
2. Select and compute variables based on the data. Generally, the individual-tree basal area increment is affected by three groups of variables: tree size, competition, and site condition^27,28,29,30.
   1. Consider tree-size effects such as DBH₁, square of DBH₁ (), the inverse transformation of DBH₁ (1/DBH₁), and the common logarithm of DBH₁ (logDBH₁) or combinations of them.
   2. Consider competitive effects such as both one- and two-sided indices of competition to more comprehensively quantify the level of competition experienced by a tree, as well as its social position within the stand. One-sided competition include BAL and the relative density index (RD=DBH₁/QMD); two-sided competition include NT, and BA.
      NOTE: The distance dependent competition indices should be considered if data is available.
   3. Consider site effects such as aspect (ASP), slope (SL), and EL. SL and ASP should be included using Stage’s transformation³¹.
3. Select log( - +1) ( denotes square of DBH₂) as the dependent variable.
4. Develop the basic model using the stepwise regression method. Ensure that the model is biologically reasonable and exhibits significant differences between independent variables. Utilize the variance inflation factor (VIF) to check for multicollinearity.
5. Leave the independent variables with p < 0.05 and VIF < 5 in the basic model.
6. Output the basic model results and the residual plot. The basic model produced here serves as a foundation for the further development of a mixed-effects model.

3. Linear mixed-effects model development with the package “nlme” in R software

1. Read the model development dataset and load the Package “nlme”.
   >model.development.dataset=read.csv("E:/DATA/JoVE/modelingdata.csv", 
   header=TRUE)
   >library(nlme)
2. Select sample plots as random effects to develop the mixed-effects model.
3. Fit all possible combinations of random effects with the maximum likelihood (ML) method and output the results.
   >Model<-lme(Y~1/DBH₁+BAL+NT+EL,data=model.development.dataset, 
   method="ML", random =~1|PLOT)
   >summary(Model)
   1. Set random =~1 is the intercept to random parameters. Change the random statements until all combinations are fitted. For example, to set 1/DBH₁ and BAL as random parameters, the code is as follows: random =~1/DBH₁+BAL-1. In addition, in the process of fitting, the codes may report errors due to the nonconvergence of the fitted model.
4. Select the best model by Akaike’s information criterion (AIC), the Bayesian information criterion (BIC), the logarithm likelihood (Loglik), and likelihood ratio test (LRT).
   >anova(Model.1, Model.6)
   >anova(Model.6, Model.23)
   >anova(Model.23, Model.30)
5. Determine the structure of R_i. Address the heteroscedasticity and autocorrelation of R_i³². The R_i is written as follows:
      (1)
   where σ² is an unknown scaling factor that is equal to the model residual variance, G_i is a diagonal matrix describing heteroscedasticity, and Γ_i is a matrix describing autocorrelation.
   1. Observe whether the residuals have heteroscedasticity from the residual plot. If there is heteroscedasticity (the residuals have a clear pattern or trend), introduce three frequently used variance functions—the constant plus power function, the power function, and the exponential function—to model the errors variance structure.
      >Model.30.1<-lme(Y~1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",random=~1/DBH₁+BAL+NT|PLOT,
      weights=varConstPower(form=~ fitted(.)))
      >summary(Model.30.1)
      >Model.30.2<-lme(Y~1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",random=~1/DBH₁+BAL+NT|PLOT,
      weights=varPower(form=~ fitted(.)))
      >summary(Model.30.2)
      >Model.30.3<-lme(Y~1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",random=~1/DBH₁+BAL+NT|PLOT,
      weights=varExp(form=~ fitted(.)))
      >summary(Model.30.3)
   2. Determine the best variance function for the model according to the AIC, BIC, Loglik, and LRT.
      >anova(Model.30, Model.30.1)
      >anova(Model.30, Model.30.2)
      >anova(Model.30, Model.30.3)
   3. Introduce three commonly used autocorrelation structures—the compound symmetry structure (CS), first-order autoregressive structure [AR(1)], and a combination of first-order autoregressive and moving average structures [ARMA(1,1)]—to account for autocorrelation.
      >Model.30.3.1<-lme(Y~1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",
      random=~1/DBH₁+BAL+NT|PLOT, weights=varExp(form=~fitted(.)), corr= corCompSymm())
      >summary(Model.30.3.1)
      >Model.30.3.2<-lme(Y~1/DBH₁+BAL+NT+EL,data=model.development.dataset,  method="ML",
      random=~1/DBH₁+BAL+NT|PLOT,weights=varExp(form=~ fitted(.)), corr=corAR1())
      >summary(Model.30.3.2)
      >Model.30.3.3<-lme(Y~1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",
      random=~1/DBH₁+BAL+NT|PLOT,weights=varExp(form=~ fitted(.)), corr=corARMA(q=1,p=1))
      >summary(Model.30.3.3)
   4. Determine the best autocorrelation structure according to the AIC, BIC, Loglik, and LRT.
      >anova(Model.30.3, Model.30.3.2)
      NOTE: The G_i and Γi cannot be defined if there is no heteroscedasticity and autocorrelation.
   5. Output the final results of the mixed-effects model using the restricted maximum likelihood (REML) method.
      >Mixed.model<-lme(Y~1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="REML",random=~1/DBH₁+BAL+NT|PLOT,
      weights=varExp(form=~ fitted(.)), corr=corAR1())
      >summary(Mixed.model)

4. Bias correction

1. Transform the predicted values of basal area increment using the final model on a logarithmic scale to the original scale. However, such a linear back transformation of predicted value from a log-transformed model produces an associated log-transformation bias. To deal with the log-bias, a correction factor was derived and integrated into the prediction equation, which estimates actual predicted basal area increment for a given tree [Equation (2)]:
      (2)
   where is predicted logarithmic value of basal area increment from the model, while is the predicted back transformed value of basal area increment for 5 years after correcting for log-transformation bias. is variance from random effects at plot and σ² is residual variance.
2. Convert basal area increment ( ) to the diameter increment.

5. Model prediction and evaluation

1. Prepare the model validation dataset produced in section 1.2 for prediction.
2. Use the linear mixed-effects model to predict individual-tree basal area increment. The random components were calculated using the following best linear unbiased predictor:
      (3)
   where is a vector for the random components; is the variance-covariance matrix for between-plots variability; is the design matrix for the random components acting at the complementary observations; is the residual vector whose components are given by the difference between the basal area increments and the predicted increments using the fixed-effects model.
3. Evaluate and compare the predictive ability of the basic model and the linear mixed-effects model using the following three statistical indicators^23,33.
      (4)
      (5)
      (6)
   where obj_i is the basal area increments, est_i is the predicted basal area increments, is the mean of observations, and N is the number of observations.

Results

The basic basal area increment model for P. asperata was expressed as Equation (7). The parameter estimates, their corresponding standard errors, and the lack-of-fit statistics are shown in Table 2. The residual plot is shown in Figure 1. Pronounced heteroscedasticity of the residuals was observed.
   ...

Discussion

A crucial issue for the development of mixed-effects models is to determine which parameters can be treated as random effects and which should be considered fixed effects^34,35. Two methods have been proposed. The most common approach is to treat all parameters as random effects and then have the best model selected by AIC, BIC, Loglik, and LRT. This was the method employed by our study³⁵. An alternative is to fit basal area increment model...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Computer	acer
Microsoft Office 2013
R x64 3.5.1