Development of a Mixed-E ﬀ ects Individual-Tree Basal Area Increment Model for Oaks ( Quercus spp.) Considering Forest Structural Diversity

: In the context of uneven-aged mixed-species forest management, an individual-tree basal area increment model considering forest structural diversity was developed for oaks ( Quercus spp.) using data collected from 11,860 observations in 845 sample plots from the 7th (2004), 8th (2009), and 9th (2014) Chinese National Forest Inventory in Hunan Province, south-central China. Since the data was longitudinal and had a nested structure, we used a linear mixed-e ﬀ ects approach to construct the model. We also used the variance function and an autocorrelation structure to describe within-plot heteroscedasticity and autocorrelation. Finally, the optimal mixed-e ﬀ ects model was determined based on the Akaike information criterion (AIC), Bayesian information criterion (BIC), log-likelihood (Loglik) and the likelihood ratio test (LRT). The results indicate that the reciprocal transformation of initial diameter at breast height (1 / DBH), relative density index (RD), number of trees per hectare (NT), elevation (EL) and Gini coe ﬃ cient (GC) had a signiﬁcant impact on the individual-tree basal area increment. In comparison to the basic model developed using least absolute shrinkage and selection operator (LASSO) regression, the mixed-e ﬀ ects model performance was greatly improved. In addition, we observed that the heteroscedasticity was successfully removed by the exponent function and autocorrelation was signiﬁcantly corrected by AR(1). Our ﬁnal model also indicated that forest structural diversity signiﬁcantly a ﬀ ected tree growth and hence should not be neglected. We hope that our ﬁnal model will contribute to the scientiﬁc management of oak-dominated forests.


Introduction
The more than 500 extant species of oak (Quercus spp.) are widely distributed across the Northern Hemisphere, including Mesoamerica [1,2]. Oak-dominated forests are of great ecological and economical importance [3,4]. For instance, oaks are extensively used in soil and water conservation and restoration efforts since they have strong, adventitious root systems and exhibit a high tolerance to waterlogging [5][6][7]. Additionally, oaks are well adapted to fire and hence are frequently used for constructing forest fire belts [8]. Oaks also play an important role in maintaining biodiversity. For instance, Caprio and Ellena [9] reported that the retention of native oaks is the key factor for the conservation of winter bird diversity in local deciduous woods. In addition to the ecological functions, oak timber is distinguished for its great strength, durability and beauty [10]. As a high-grade We randomly divided the data into two datasets, with 80% of the data used for model fitting (9739 observations from 676 sample plots) and the other 20% for model validation (2121 observations from 169 sample plots).

Development of Basic Model
The individual-tree basal area increment (BAI) can be described as a function of tree size (SIZE), competition pressures (COMP) and site conditions (SITE) in a stand [26,56,58,59]. We also considered the effects of forest structural diversity (DIVE) on BAI. Strictly speaking, structural diversity includes three categories-tree species diversity, tree size diversity, and tree position diversity. Unfortunately, We randomly divided the data into two datasets, with 80% of the data used for model fitting (9739 observations from 676 sample plots) and the other 20% for model validation (2121 observations from 169 sample plots).

Development of Basic Model
The individual-tree basal area increment (BAI) can be described as a function of tree size (SIZE), competition pressures (COMP) and site conditions (SITE) in a stand [26,56,58,59]. We also considered the effects of forest structural diversity (DIVE) on BAI. Strictly speaking, structural diversity includes three categories-tree species diversity, tree size diversity, and tree position diversity. Unfortunately, due to the limited coordinate information available in the data, only tree species diversity and tree size diversity were considered. The relationship is expressed as follows: The following independent variables were evaluated in our study: (1) Tree size effects: initial DBH, square of initial DBH (DBH 2 ), reciprocal transformation of initial DBH (1/DBH) and natural logarithm of initial DBH (logDBH). (2) Competition pressures effects: number of trees per hectare (NT), stand basal area per hectare (BA), the sum of the basal area (m 2 /ha) in trees with DBHs larger than the subject tree's DBH (BAL), the ratio of BAL and DBH (BAL/DBH), and the relative density index (RD = DBH/QMD). (3) Site conditions effects: elevation (EL), slope (SL) and aspect (ASP). Additionally, slope and aspect were combined using Stage's transformation, i.e., SLCos = SL × cos(ASP), SLSin = SL × sin(ASP) [60]. (4) Structural diversity effects: 1) tree species diversity: the Shannon-Wiener index (SHI), Pielou index (PI) and Simpson's index (SII); 2) tree size diversity: the Gini coefficient (GC) and the standard deviation of the DBHs (SDDBH) [61].
where p i is the proportion of basal areas in the ith species.
where SHI is the Shannon-Wiener index and S is the total number of species in a sample, across all samples in a dataset.
where p i is the proportion of basal areas in the ith species and n is the number of species observed.
where ba t is basal area for the tree in rank t (m 2 /ha) and t is the rank of a tree in order from 1, . . . , n.
Each of these variables has been documented to have a significant effect on individual-tree basal area increment. All variables used in the model development and their detailed descriptions are provided in Table 1.
Different forms of dependent variables can result in different model performances [26,62,63]. The most frequently used dependent variables are the diameter increment (DBH 2 -DBH, in which DBH 2 denotes the DBH measured after a 5-year growth), squared diameter increment (DDS = DBH 2 2 − DBH 2 ), squared diameter increment plus a constant value of one (DDS + 1), and the natural logarithmic transformation of each [25,47,62]. These potential dependent variables were evaluated using the following lack-of-fit statistics, i.e., absolute bias (Bias), root mean square error (RMSE), and the coefficient of determination (R 2 ). Additionally, residual plots were also produced to inspect the homogeneity of the variance and normality of the residuals. After comparing and testing these dependent variables, log(DDS + 1) was selected as the dependent variable for our basal area increment model since it performed best in terms of the normality and homogeneity of the residuals, and lack-of-fit statistics. Similar results were also reported by Calama and Montero [47], Adame et al. [63] and Lhotka and Loewenstein [25]. In this study, the basic individual-tree basal area increment model was developed by least absolute shrinkage and selection operator (LASSO) regression. The R package used for LASSO regression was glmnet [64]. For this method, the λ tuning parameter was determined by cross-validation using the fitting data with the mean square error as the criterion to obtain the candidate model. In order to check for multicollinearity amongst the variables, the variance inflation factor (VIF) was calculated. All variables with a VIF larger than 5 were removed to minimize overfitting [47,52,65]. Additionally, the relative importance refers to the quantification of an individual regressor's contribution to a multiple regression model. In order to examine the relative effects of independent variables entering the final basic model, we calculated the relative importance values using the R Package "relaimpo".

Development of Linear Mixed-Effects Model
Linear Mixed-Effects Model Since our data was longitudinal (repeated measurement for each tree) and had a nested structure (tree nested within sample plots), we used a linear mixed-effects approach using sample plot as a random effect to produce our individual-tree basal area increment model. The general form of a linear mixed-effects model is as follows: where y i represents the vector of the observations for the dependent variable, x i is a design matrix including covariates and derivative terms associated with fixed effects parameters, β is the vector of fixed parameters, z i is a design matrix for the random effects, and b i is the vector of random parameters. e i is error term, and R i is a positive definite variance-covariance matrix for the error term which includes the variance and correlation.

Construction of Parameter Effects
The parameters are divided into mixed-effects and fixed effects parameters. This is important when developing the mixed-effects model. In this study, all possible combinations were fitted and then the best combination with convergence was selected according to the Akaike information criterion (AIC), Bayesian information criterion (BIC), log-likelihood (Loglik) and the likelihood ratio test (LRT) [66,67]. The AIC, BIC and LRT are defined as follows: where d is the number of estimated parameters of the model (including the intercept), n is the number of observations, and L 1 and L 2 are the likelihood value of different models.

Determination of Random Effects Variance-Covariance Structure
The variance-covariance matrices for the random effects are common to all sample plots. These matrices account for the variability between sample plots. In this study, we selected three structures to describe the random effects variance-covariance matrix. They were a diagonal matrix, a compound symmetry and general positive-definite matrix. The hypothetical 3 × 3 variance-covariance matrices are shown in Equations (10)- (12).
Compound symmetry = is the covariance between the b u th and the b v th random effects, intercepting at σ 2

Determining the Matrix of R i
Heteroscedasticity and autocorrelation should be corrected when performing a regression analysis; if they are not, the computed standard errors and estimation interval will be incorrect, which will result in spurious regression problems [56,65]. Because the modeling data in this study was obtained from repeated measurements, the heteroscedasticity and autocorrelation of R i should be addressed [68]. The heteroscedasticity could be observed through the residual plot and quantile-quantile (QQ) plot. Additionally, we also conducted the Breusch-Pagan test to detect heteroscedasticity using the R Package "lmtest." Because our data was of a longitudinal structure, i.e., two measurement periods, autocorrelation might have occurred. We therefore performed the Durbin-Watson test to investigate the autocorrelation using the R Package "car." If there were heteroscedasticity and autocorrelation, we could introduce variance functions and autocorrelation structures to remove the heteroscedasticity of the residuals and correct autocorrelation among the observed values for the same tree obtained in different growth periods. The R i was written as follows [69]: where σ 2 is a scaling factor which is equal to the residual variance of the developed model, G i is a diagonal matrix which describes the heteroscedasticity and Γ i is a matrix showing the autocorrelation structure of errors. Three different variance functions, i.e., the exponential function (Equation (14)), the power function (Equation (15)) and the constant plus power function (Equation (16)), were tested to remove heteroscedasticity [70]. varExp varPower where u i is the estimated value based on fixed parameters of the linear mixed-effects model and α, β are estimated parameters of variance functions. A first-order autoregressive structure (AR(1)), the compound symmetry structure (CS), and a combination of first-order autoregressive and moving average structures (ARMA(1,1)) were selected to correct autocorrelation [70]. The Equations (17)- (19) are the detailed matrix forms.
where ρ is the autoregressive parameter, υ is a moving average component and σ 2 is the residual variance.

Random Parameters Prediction
The random parameters in Equation (20) were predicted using the best linear unbiased predictor (EBLUP) [46,71] whereb i is a vector of random effects parameters of sample plots;D is the m × m variance-covariance matrix for random effects (m is the number of random effects parameters);Ẑ i is the design matrix for the random effects parameters for the observations;ê i is the residual vector whose components are given by the difference between the observations of individual-tree basal increment and predicted values from the mixed-effects model with only the fixed effects.

Model Prediction and Evaluation
As stated in our data and pre-analysis section, 20% of the plot data was used for model prediction. Using a mixed-effects model for prediction, two different predictive situations can be considered, i.e., the population average (PA) response, and the subject-specific (SS) response [47,67,72]. (a) Population average responses (also termed marginal mean responses) are predicted using only the fixed effects part of a mixed-effects model if no additional information on random effects is available [47,72]. In this case, the expected value of the EBLUP's for the random effects are zeroes, and the population average responses could hence be expressed as E(y) = Xβ. (b) Subject-specific responses (also called conditional mean responses) can be obtained using both fixed effects part and random effects part of a mixed-effects model, when the random parameters vector b can be predicted [47,65]. The subject-specific responses can be written as E(y b) = Xβ + Zb. In this study, the random parameters for the evaluation data (the 20% of plot data) was predicted using the same method in Section 2.2.3.
where ob j i is the observed value, est i is the estimated value, ob j i is mean observed value, and N is the number of observations.

Basic Model
The basic model using LASSO regression originally had nine independent variables. The tuning parameter λ, obtained by cross-validation, was 3.208 × 10 −3 . To avoid over overfitting and multicollinearity, we eliminated the variables of the candidate model with VIF > 5. This left 1/DBH, RD, NT, EL and GC to be entered into our basic model. Equation (24) was determined as our final basic model, to which random effects would be added to produce our linear mixed-effects individual-tree basal area increment model. The parameter estimates and their statistical description are shown in Table 2.

Random Parameter Effects
There are 63 (C 1 6 + C 2 6 + C 3 6 + C 4 6 + C 5 6 + C 6 6 = 63) potential different combinations of random effects for Equation (24) while considering all independent variables and the intercept in the basic model. A total of 49 combinations reached convergence when fitted to the data. Amongst these 49 mixed-effects models, Equation (25) yielded the smallest AIC (4836.349) and BIC (4994.392) and largest LogLik (-2396.174). In addition, the LRT also indicated the best performance of Equation (25). Equation (25) was the resulting mixed-effects model.
where β 1 -β 6 are the fixed effects parameters, and b 1 -b 5 are the random effects parameters.

Random Effects Variance-Covariance Structure
The three different random effects variance-covariance structures, i.e., a diagonal matrix, a compound symmetry and general positive-definite matrix, were introduced to Equation (25), respectively, and their performances were compared. The general positive-definite matrix showed the best performance in terms of AIC, BIC, Loglik and LRT (Table 3). Therefore, the model with general positive-definite matrix was selected for our linear mixed-effects individual-tree basal area increment model.

Error Variance-Covariance Structure
We used three variance functions to remove the heteroscedasticity of the residuals and their performances were compared in terms of AIC, BIC, Loglik and LRT. The results showed that the mixed-effects models with the variance function exhibited better performance than the one without considering the variance function (Table 4). Performance also differed amongst the mixed-effects models with different variance functions. According to AIC, BIC, Loglik and LRT, the exponent function was determined as the best variance function. The three correlation structures, i.e., AR(1), ARMA(1,1) and CS, were then used to correct the autocorrelation. According to the AIC, BIC, Loglik and LRT, AR(1) performed the best and hence was selected as the best autocorrelation structure for our mixed-effects model.

Final Form of the Linear Mixed-Effects Individual-Tree Basal Area Increment Model
After the determination of parameter effects, random effects variance-covariance structure and error variance-covariance structure, the final form of the linear mixed-effects individual-tree basal area increment model for oaks in Hunan Province, south-central China, was proposed in Equation (26) and Equation (27). The statistical analysis results of variables of the final mixed-effects model for oaks are shown in Table 5.  The residual and QQ plots of our final mixed-effects model (Equation (26)) were shown in Figure 3. Compared with the basic model (Equation (24)), we observed a significant improvement in the residuals in terms of homogeneity and normality.
Forests 2019, 10 FOR PEER REVIEW 11 The residual and QQ plots of our final mixed-effects model (Equation (26)) were shown in Figure  3. Compared with the basic model (Equation (24)), we observed a significant improvement in the residuals in terms of homogeneity and normality.

Model Prediction and Evaluation
The basic model (Equation (24)) and mixed-effects model (Equation (26)) were used for prediction using the validation data. Table 6 lists the three lack-of-fit statistics of the basic model and the mixed-effects model. We found that our mixed-effects model outperformed the basic model regardless of PA and SS. Additionally, Figure 4 shows predicted DBH increments over DBH with all

Model Prediction and Evaluation
The basic model (Equation (24)) and mixed-effects model (Equation (26)) were used for prediction using the validation data. Table 6 lists the three lack-of-fit statistics of the basic model and the mixed-effects model. We found that our mixed-effects model outperformed the basic model regardless of PA and SS. Additionally, Figure 4 shows predicted DBH increments over DBH with all other predictors set to their mean. We could observe that the model shows the sigmoidal curve, i.e., decreasing DBH increment for trees with very large DBH (e.g., 59 and 61 cm in this study), though the DBH 2 did not enter into our model. Basic model (Equation (24)) Mixed-effects model (Equation (26)) (24)) and mixed-effects model (Equation (26)) over DBH with all other predictors set to their mean for oaks in Hunan Province, south-central China.

Discussion
In this study, we developed a linear mixed-effects individual-tree basal area increment model for oaks considering forest structural diversity in Hunan Province, south-central China, using data from 845 sample plots of CNFI measured three times each. The model was described as a stochastic process, where the fixed component explained the mean value for the basal area increment and the random part incorporated unexplained residual variability acting at the level of sample plot (e.g., moisture, soil parameters, nutrient content, etc.) [47,63]. Our final mixed-effects model included the variables 1/DBH, RD, NT, EL and GC in the fixed component. Additionally, the relative importance values for 1/DBH, RD, NT, EL and GC were 33.15%, 38.78%, 13.94%, 12.71% and 1.42%, respectively. From the perspective of forest management, the independent variables that have higher relative importance value should be given higher priority in forest management planning.
Consistent with results reported by Cao et al. [73], Uzoh and Oliver [56], Lei et al. [74], Pokharel and Dech [59], and Timilsina and Staudhammer [75], 1/DBH was significantly negatively related to basal area increment. It had the second largest effect on basal area increment. RD, a distanceindependent individual-tree-level competition index, was found to be positively related to basal area increment, which had more effect on basal area increment than any other variable. Similar findings were documented by Lei et al. [74] and Yan [32]. Both 1/DBH and RD were indicators of individual  (24)) and mixed-effects model (Equation (26)) over DBH with all other predictors set to their mean for oaks in Hunan Province, south-central China.

Discussion
In this study, we developed a linear mixed-effects individual-tree basal area increment model for oaks considering forest structural diversity in Hunan Province, south-central China, using data from 845 sample plots of CNFI measured three times each. The model was described as a stochastic process, where the fixed component explained the mean value for the basal area increment and the random part incorporated unexplained residual variability acting at the level of sample plot (e.g., moisture, soil parameters, nutrient content, etc.) [47,63]. Our final mixed-effects model included the variables 1/DBH, RD, NT, EL and GC in the fixed component. Additionally, the relative importance values for 1/DBH, RD, NT, EL and GC were 33.15%, 38.78%, 13.94%, 12.71% and 1.42%, respectively. From the perspective of forest management, the independent variables that have higher relative importance value should be given higher priority in forest management planning.
Consistent with results reported by Cao et al. [73], Uzoh and Oliver [56], Lei et al. [74], Pokharel and Dech [59], and Timilsina and Staudhammer [75], 1/DBH was significantly negatively related to basal area increment. It had the second largest effect on basal area increment. RD, a distance-independent individual-tree-level competition index, was found to be positively related to basal area increment, which had more effect on basal area increment than any other variable. Similar findings were documented by Lei et al. [74] and Yan [32]. Both 1/DBH and RD were indicators of individual tree competition status in a forest stand and their relationship with basal area increment suggested that larger trees have stronger competitive ability for resources, especially for light which is considered as the major limiting resource for individual-tree diameter growth [65,76]. For instance, Ma et al. [77] reported the nutrient use efficiency of the China fir in mature stands was almost twice that in young stands. RD has frequently been selected to represent individual tree competition effects due to its performance in models [74,78]. In comparison to the individual-tree-level competition indices, i.e., 1/DBH and RD, NT can be thought to represent stand-level competition, which had the third largest effect on basal area increment. In our study, a negative relationship was observed between NT and basal area increment, indicating that stand-level competition reduces individual tree growth.
There are two different types of competition, i.e., one-and two-sided competition [18,79]. In one-sided competition, larger trees are not affected by their smaller neighbors. In two-sided competition, by contrast, resources are shared by all trees either equally or proportionally to size [80]. It is commonly assumed that one-sided competition is driven by the availability of aboveground resources, whereas two-sided competition is more reflective of belowground competition [81]. Therefore, many growth and yield models consider 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 [18,79]. Thus, we also considered one-sided (BAL, BAL/DBH, RD) and two-sided (NT, BA) competition when developing the basic model in the present study. RD and NT were included in the model to represent the comprehensive effects of aboveground and belowground competition.
Although we included three site-level independent variables, i.e., EL, SLSin, and SLCos, only EL was left in our final model. EL had the fourth largest effect on basal area growth. We observed a negative correlation between EL and basal area increment, which suggested that individual trees exhibited faster growth in the lower elevation. This negative correlation could be attributed to the shorter growing season at higher elevations [60,65,82]. Additionally, many authors argued that EL could indirectly affect tree growth by altering temperature, moisture, light, soil nutrient availability and other chemical and physical agents in a forest stand [83,84].
GC had the smallest influence on basal area increment. The negative correlation between GC and basal area increment suggested that variation in large tree size reduces basal area increment. Similar results were also reported by Liang et al. [34], Cordonnier and Kunstler [85] and Bourdier et al. [86]. Bourdier et al. [86] attributed the negative tree size inequality effect to the reduced total light interception of the stand, which reduces light use efficiency. Their detailed explanations are as follows: for stands with a comparable basal area and mean diameter, the increase in GC indicates greater difference in the DBHs between the biggest and small trees. Because tree canopy width and depth increase asymptotically with tree size, the light interception efficiency, which is determined by the canopy characteristics, also increases asymptotically. Therefore, a decrease in the amount of light intercepted might occur when the gain in light intercepted by the bigger trees is unlikely to compensate for the loss of light intercepted by the smaller trees. In comparison, the reduced light use efficiency with increases in the GC may be attributed to the fact that large individuals already intercept more light and an increase in light only slightly improves their growth, whereas small individuals live in low light conditions where supplementary light has a stronger effect on growth. Therefore, from a forest management perspective, extremely large trees, which have significant relative dominance and absolute advantages in light competition, should be cut to reduce tree-size inequality in a stand. Consequently, the total light interception and light use efficiency of a stand might increase.
Due to the hierarchical structural of our data, we used the sample plot as a random effect to produce the individual-tree basal area increment model. Our results indicated a significant improvement in model performance after introducing the random effects. For example, the AIC dropped from 6317.446 to 4705.766, and BIC from 6367.733 to 4878.177. Similar results have been extensively documented by other authors [65,67,75].
Although introducing random effects could correct or reduce autocorrelation and heteroscedasticity [65,68,87,88], our mixed-effects model still exhibited autocorrelation and heteroscedasticity. We therefore further introduced three variance functions and three correlation structures to refine our model. Finally, based on the lack-of-fit statistics, the exponent function and AR(1) were determined as the optimum option for correcting the heteroscedasticity and autocorrelation of the residuals. Similar results were also reported by Calama and Montero [47] and Fabian et al. [56]. Additionally, the model prediction and evaluation using validation data further supported the conclusion that the mixed-effects approach had significantly improved the model's predictive performance, either in PA or SS in comparison to the basic model without random effects. For instance, the Bias decreased from 0.2239 to 0.1716, the RMSE decreased from 0.2776 to 0.2167 and the R 2 increased from 0.3663 to 0.6140.
Although the model performance was improved after we used the mixed-effects method to develop the individual-tree basal area increment model for oaks, there are some limitations of our final mixed-effects model. Firstly, we only considered the autocorrelation structure to correct autocorrelation. However, Zhao et al [65] reported that modelling the autocorrelation structure alone could not completely match the possible common effects and suggested introducing random period effect to describe autocorrelation, if data is sufficient. Pokharel and Dech [59] directly introduced random period effects to develop their diameter growth model. In the present study, unfortunately only three period observations were available for the same trees and hence it makes no sense to employ period as a random effect.
Secondly, climatic change has been widely reported to influence forest growth, productivity, tree species composition and distribution [89][90][91]. For instance, Battles et al. [92] assessed the impact of climate change on the productivity and health of the mixed-conifer forest in California and found that the productivity of stem volume increment was reduced by 19% in the most extreme climate change. Unfortunately, we did not include climatic variables due to lack of data. It is imperative to integrate climatic variables into forest growth models when data is available.
There are many individual-tree forest simulators throughout the world, for instance, the FVS, TASS, PROGNAUS in the US, the SILVA, and BWINPro in Germany and CAPSIS in France [93][94][95]. For a comprehensive prediction of forest dynamics, this individual-tree simulation system not only included individual-tree basal area or diameter increment models, but also integrated an individual-tree mortality and recruitment model. Therefore, the development of an individual-tree mortality and recruitment model of oak species is strongly recommended in the future.

Conclusion
In this study, we produced a linear mixed-effects individual-tree basal area increment model for oaks (Quercus spp.) that considers forest structural diversity. Our final model showed that it is necessary to consider forest structural diversity in model development. Additionally, compared with the basic model built by LASSO regression, the performance of the mixed-effects model was greatly improved. The model is a useful tool to predict the basal area growth of individual trees and propose relevant forest management strategies in uneven-aged mixed-species forests.
Author Contributions: All authors made significant contributions to the manuscript: J.M. and W.Z. conceived, designed and performed the experiments; W.W. analyzed the data and results; X.C., W.Z. and J.W. contributed reagents/materials/analysis tools; J.M. and W.W. are the main authors who developed and revised the manuscript.