Modelling Individual Tree Diameter Growth of Quercus mongolica Secondary Forest in the Northeast of China

Abstract

Quercus mongolica secondary forest is widely distributed in the northeast of China, but it usually has low productivity, unstable structure, poor health, and low biodiversity. Diameter is a tree variable that is commonly used for forest growth measurement, to provide the basis for forest management decision. Two level generalized linear mixed effects individual tree diameter growth model were developed using data from two times surveys of 12 Q. mongolica secondary forest permanent plots that were distributed among Wangqing forest farms. Random effects of 14 tree species and 12 plots were introduced into the basic model consisting of three factors: tree size, competition of surrounding trees, and site quality. The results showed that initial diameter at breast height(DBH) was the most important variable affecting diameter growth, followed by competition, while the effect of site quality on diameter growth was not significant. Compared with the basic model, the prediction accuracy of the mixed effect model was improved by 17.69 %, where R² reached to 0.6805, indicating that it is suitable for the individual-tree diameter growth prediction of the secondary forest of Q. mongolica.

1. Introduction

As forest management is a long-term activity, the response of growth indicators to management measures may not be reflected in the short term. Forest growth models can estimate the long-term production by using the growth law of trees or stands, and then evaluate the impact of management measures to provide a theoretical basis for the formulation of forest management decision, so it is highly valued. The research of the forest growth model mainly includes model construction, parameter determination, high-precision estimation, and the prediction of individual tree level and stand level [1,2]. A large number of growth and harvest models have been developed and applied to forestry production practice [3]. The stand growth model may not able to represent the growth dynamics of individual trees and is only applicable to regular and even-aged stands, as tree sizes and competition largely varies in structurally complex stands [4,5,6]. The individual tree growth model is more flexible, detailed and accurate for the prediction of individual tree growth dynamics, so compared with the stand level model, it can more accurately evaluate the impact of different management measures on growth [7]. In addition, it is more suitable for uneven-aged mixed forest, because it can more accurately show the competition relationship between different tree species or different sizes of trees [2,8]. The individual tree growth models usually consist of diameter increment, height increment, crown ratio, and mortality models as core models [9] which can be used as input in forest simulators for decision making in sustainable forest management [10,11].

In the past few hundred years, the focus of forest management was mostly around the cultivation of pure plantation, and the series of models for pure forest were relatively perfect [12]. However, compared with pure forest, mixed forest has higher ecological and social values, so it has attracted more and more attention, and models used for the growth prediction of mixed forest are also increasingly needed [13,14,15,16]. The least squares method is the most commonly used method for estimating model parameters. Due to the lack of independence of the observed data [17] and the existence of time correlation and spatial heterogeneity among the data [18], the least square method leads to a large error in the prediction results of the model in the mixed forest with multiple tree species. With the development of computer technology, model estimation methods are also developing. For example, the introduction of measurement error model and mixed effect improves the application range of the model and the accuracy of simulation and prediction [19,20,21]. Mixed effect model has been widely used in the forestry field [22] as it can solve the above problems effectively by setting fixed effect parameters and random effect parameters [23].

According to the data of China’s ninth forest inventory, Quercus stands are ranked first in terms of volume and area [24], moreover, oak management has also changed from mainly producing fuelwood to producing high-value wood. Q. mongolica is one of the main tree species of coniferous and broad-leaved mixed forest in Northeast and North China. It has well-developed root system, strong adaptability, and stronger resistance to natural disasters, therefore, it plays an important role in maintaining regional ecological balance and ecosystem restoration [25]. Because it is usually distributed in poor site conditions, it has not been paid attention in management. However, with the increasing environmental threats caused by climate change, oak species have been paid more and more attention because of their good environmental adaptability and high-quality wood production capacity [26,27,28,29].

Despite its higher economic and ecological benefits compared to other tree species, the individual tree-based growth model by introducing random effects for Q. mongolica in China is still lacking. This study aims to develop a tree diameter growth model using variables describing the tree size, stand competition, and site quality. The linear mixed-effect modelling method was applied to reflect the influence of tree species and the difference between sample plots. Moreover, compared with introducing variables as fixed effects, it is more concise to take them as random effects. The model will be used to predict the growth dynamics of Q. mongolica secondary forest, while the factors that influence the diameter growth are analyzed to provide the basis for guiding the management of this kind of forest.

2. Materials and Methods

2.1. Data Materials

The study area is located in Tazi Gou Forest farm (129°56′–131°04′ E, 43°05′–43°40′ N), Wangqing Forestry Bureau, Jilin Province, China (Figure 1). This region is dominated by low mountains and hills, with an altitude between 300–1200 m and a slope between 5° and 25°. It belongs to a temperate continental monsoon climate with an annual average temperature of about 3.9 °C and annual precipitation 600–700 mm, mostly in summer. The lowest temperature is about −32 °C in January and the highest temperature is about 22 °C in July. The area is mainly dark brown soil with fertile soil texture, acidic pH, and high humus content. The main tree species are Q. mongolica, Betula platyphylla, Pinus koraiensis, Acer pictum, Tilia mandshurica, Populus ussuriensis, Larix olgensis, Betula dahurica, Fraxinus mandshurica, Abies nephrolepis, Ulmus davidiana var. japonica, Picea jezoensis, and Phellodendron amurense, among others.

Twelve one-hectare permanent sample plots (100 m × 100 m each) were established in Q. mongolica secondary forest in August 2013. After the establishment of the sample plot, measurements were made and then the thinning operation was taken. The measurements of plots include height (H), diameter at breast height (DBH), crown width and coordinates of all living trees (D ≥ 5.0 cm) as well as slope, elevation, and other variables. The thinning operation was designed for comparing the influence of traditional thinning and crop tree releasement. The felled trees of traditional thinning are of poor quality, low vitality, and high density, while the felled trees of crop tree releasement are the competitors of crop trees.

In August 2018, 12 sample plots were remeasured. The data for establishing the diameter growth model are from the two measurements of investigation data of 12 Q. mongolica sample plots in 2013 and 2018. The calculation of independent variables in 2013 is based on the stand status after cutting. A total of 9365 trees with measurements taken twice were used for the model establishment. The basic information at the beginning of the establishment of the sample plots is shown in

Table 1. General information of oak secondary forest.

Plot	Thinning Intensity of Volume/%	Aspect	Slope Position	Gradient/°	Altitude /m	Basal Area /(m2/ha)	Volume /(m3/ha)	DBH/cm			Height/m
								Min.	Max.	Mean	Min.	Max.	Mean
1	6.37	West	Middle	8	705	19.29	136.10	5.0	66.0	14.6	2.0	29.6	11.0
2	6.67	West	Middle	8	738	19.21	134.57	5.0	59.1	13.7	0.5	24.1	11.1
3	4.04	West	Middle	7	721	20.05	140.85	5.0	55.7	12.9	1.3	28.5	10.2
4	2.88	West	Middle	8	741	21.34	153.49	5.0	50.2	13.6	1.5	25.8	10.9
5	0	West	Up	10	794	27.01	192.62	5.0	53.8	15.1	1.4	21.8	10.0
6	0	West	Middle	8	639	13.65	107.11	5.0	50.7	15.0	1.8	27.1	13.3
7	3.46	West	Middle	8	637	19.18	140.37	5.0	50.9	13.1	1.5	21.8	8.9
8	9.76	West	Middle	8	635	21.19	159.18	5.0	58.5	12.8	3.2	31.1	10.9
9	1.74	West	Middle	7	677	18.93	137.41	5.0	52.2	14.3	1.4	18.9	9.8
10	4.22	West	Middle	7	685	22.57	175.08	5.0	70.1	14.3	3.1	29.4	11.7
11	4.16	West	Middle	6	703	23.38	175.86	5.0	58.2	14.1	2.2	26.5	10.9
12	0	West	Middle	10	755	20.49	143.85	5.0	65.3	15.9	1.4	23.4	11.0

.

2.2. Basic Model

Naturally, the main factors affecting the diameter growth of individual trees are tree size, competition of surrounding trees, and site quality [30,31,32,33,34]. Combined with the biological importance and the correlation with the diameter growth of the dependent variable, a total of 14 variables were selected, among which the diameter at breast height (D), basal area (BA), tree height (H), average crown width (CW), the coefficient of slenderness (CS), and crown ratio (CR) were related to the tree size. The distance independent individual tree competition index include basal area of the trees larger in diameters than a subject tree (BAL, m² ha⁻¹), total basal area per ha (m²/ha) for all species in the stand (BASUM), the stand density index (SDI = N(D_g/20)^1.605) [35], the ratio of basal area of target tree species to stand basal area(SC) and the ratio of DBH of the target tree to the average DBH of the stand(CI). The distance-dependent individual tree competition indexes include Hegyi [36] and aCI based on the intersection angle [37], which can be depicted as:

$$Hegy{i}_i={\displaystyle \sum _{j=1}^4\left(\frac{D_j}{D_i}\cdot \frac{1}{d_{ij}}\right)}$$

(1)

where Hegyi_i is the competition index of the object tree i (ith), d_j is the DBH of the competitive tree j (jth), D_i is the DBH of the ith object tree, D_ij is the horizontal distance between the object tree and the competitive tree.

$${\mathrm{aCI}}_{\mathrm{i}}=\frac{1}{4}{\displaystyle {\displaystyle \sum }_{\mathrm{j}=1}^4}\frac{({\mathsf{\alpha }}_1+{\mathsf{\alpha }}_2\cdot {\mathrm{c}}_{\mathrm{ij}})}{180°}$$

(2)

where ${\alpha }_1=\left\{\begin{array}{l}\mathrm{arctg}(H_i/L_{ij}),H_j{>H}_i\\ \mathrm{arctg}(H_j/L_{ij}),\mathrm{otherwise}\end{array}\right.$, ${\alpha }_2=\mathrm{arctg}\left(\frac{H_j-H_i}{L_{ij}}\right)$, $c_{ij}=\left\{\begin{array}{l}1,H_j>H_i\\ 0, \mathrm{otherwise}\end{array}\right.$, where aCI_i is the competition index of the ith tree based on the intersection angle, H_i is the initial height of of the ith object tree, H_j the height of the jth neighboring tree. L_ij is the horizontal distance between the object tree and the neighboring tree. The spatial structure of forest can be quantified by describing the spatial relationships between the reference tree and its four immediate neighboring trees [38]. Therefore, the competition trees selected in this study are the 4 trees directly adjacent to the target trees.

Site index is most commonly used for site productivity, however, it is suitable for regular stands [39]. Therefore, in this study, the relationship between tree height and diameter was used to represent site productivity (SPI) [40]. It is common to use the relationship between tree height and diameter to evaluate the site productivity of mixed forest [4], as it is difficult to precisely know the site index and age.

$$SPI=1.3+\left(H_i-1.3\right)\cdot {\left(\frac{1-\exp \left(-b{D}_0\right)}{1-\exp \left(-b{D}_i\right)}\right)}^{kb}$$

(3)

where k = 14.4204, b = 0.1657, D₀ = 20 [41], H_i is the dominant height of the ith plot, D_i is the average diameter of dominant trees in the ith plot.

Therefore, the individual tree diameter growth model of Q. mongolica natural secondary forest was constructed by regression analysis, as shown below (Equation (4)).

$$\ln \left(D_2^2-D_1^2\right)=a+Size+Comp+Site$$

(4)

$$Size=b_1{D}_1+b_{2}BA+b_{3}H+b_{4}CW+b_{5}CS+b_{6}CR$$

(5)

$$Comp=c_{1}Hegyi+c_{2}aCI+c_{3}BAL+c_{4}BASUM+c_{5}SDI+c_{6}SC+c_{7}CI$$

(6)

$$Site=d_{1}SPI$$

(7)

where D₂ is the final DBH, D₁ is the initial DBH.

After determining the form of the basic model, these candidate factors are gradually introduced to measure the accuracy improvement of the basic model by different combinations. The variables with strong collinearity were excluded by combining the Variance inflation factor [42,43]. We start with graphs and relevant statistical analysis to see whether the variables are suitable for model fitting [44,45]. The variables with significant coefficients and higher fitting accuracy are selected primarily. Moreover, optimal basic model is determined by the stepwise regression method of linear model. Finally, four variables (D, CW, CS and CR) reflecting the size of the forest and four variables (SDI, aCI, SC and Hegyi) reflecting the competition were used as the independent variables to establish the individual tree diameter growth model of Q. mongolica natural secondary forest.

2.3. Two Level Nested Mixed Effect Model of Individual Tree Diameter Growth

Based on linear basic model of individual tree diameter increment, estimation models with sample plot nested in tree species two level random effects were formulated according to the linear mixed estimation modeling technique (Equation (8)):

$$\left\{\begin{array}{c}\ln {(D_2^2-D_1^2)}_{ijk}=f({\phi }_{ijk},x_{ijk})+{\epsilon }_{ijk},\\ {\phi }_{ijk}=A_{ijk}\beta +B_{ijk}{u}_{ij},\\ u_{ij}~N(0,\Psi ),\\ i=1,\dots ,M,j=1,\dots ,n_i,k=1,\dots ,n_{ij}.\end{array}\right.$$

(8)

where $\ln {(D_2^2-D_1^2)}_{ijk}$ was the diameter increment of the kth tree on the jth sample plot of the ith tree species, $M$ was the number of tree species, $n_i$ was the number of sample plots in the ith tree species, $n_{ij}$ was the number of individual tree in the jth sample plot of the ith tree species, and $f(\cdot )$ was a linear function of a parameter vector ${\phi }_{ijk}$ and a covariate vector $x_{ijk}$. The within-group error vector ${\epsilon }_{ij}={({\epsilon }_{ij1},\dots ,{\epsilon }_{ij{n}_{ij}})}^T$ that accounted for the within-group variance and correlation was assumed to follow a normal distribution that had an expectation of zero and an identity matrix variance-covariance structure. where $\beta$ was a parameter vector of fixed effects, $u_{ij}$ was an independently and normally distributed vector of random effects with zero means and variance-covariance matrix $\Psi$, $A_{ijk}$ and $B_{ijk}$ were design matrices, and $u_{ij}$ and ${\epsilon }_{ijk}$ were independent of each other.

2.4. Model Estimation and Evaluation

Generally, data sets are divided into a modeling set and a verification set. The modeling set is used to build models and estimate parameters, while the verification set is used to evaluate the extrapolation ability of models. The cross validation method was used to evaluate the prediction accuracy of the two-level mixed effect model to make better use of the existing data. The data structure of this study is based on sample plots, including 12 sample plots. In the cross-validation method, 1 sample plot was taken as the validation set each time, and the remaining 11 sample plots were taken as the modeling set for model construction. Then, the constructed model was used to predict the 1 sample plot. This cycle was repeated 12 times, and the predicted values of the 12 sample plots were obtained. The Akaike information criterion (AIC, Equation (9)), root mean square error (RMSE Equation (10)), total relative error (TRE, Equation (11)) and maximum coefficient of determination (R², Equation (12)) were calculated based on the observed and predicted values of these dependent variables. The larger R² is, the smaller AIC, RMSE and TRE are, the stronger prediction ability of the model is.

$$AIC=-2\ln \left(1\right)+2p$$

(9)

$$RMSE=\sqrt{\frac{{\displaystyle \sum {(y_i-{\widehat{y}}_i)}^2}}{n}}$$

(10)

$$TRE=\frac{{\displaystyle \sum (y_i-{\widehat{y}}_i)}}{{\displaystyle \sum {\widehat{y}}_i}}\times 100$$

(11)

$$R^2=1-\frac{{\displaystyle \sum {(y_i-{\widehat{y}}_i)}^2}}{{\displaystyle \sum {(y_i-\overline{y})}^2}}$$

(12)

where n is the total number of trees, $y_i$ is the diameter growth of the ith tree, ${\widehat{y}}_i$ is the estimated value of the ith tree diameter growth, $\overline{y}$ is the average value of all tree diameter growth, p is the number of parameters in the model, and l is the maximum likelihood function value of the model.

The parameters of the basic model and mixed effect model were estimated in R 3.6.2, and the mixed effect model was calculated by Nonlinear Mixed Effect (NLME) package [46].

3. Results

3.1. Determination of the Optimal Basic Model

The fitting evaluation indexes of each basic model with different variables are shown in

Table 2. Evaluation index for basic model of individual-tree diameter increment.

Number of Variables	Selected Variables by Stepwise Regression	RMSE	R2	TRE	AIC
1	lnD	0.7867	0.5081	4.909	21618
2	lnD, lnCW	0.7697	0.5291	4.689	21220
3	lnD, lnCW, SDI	0.7584	0.5429	4.546	20950
4	lnD, lnCW, SDI, aCI	0.7460	0.5577	4.392	20649
5	lnD, lnCW, SDI, aCI, CS	0.7398	0.5651	4.316	20498
6	lnD, lnCW, SDI, aCI, CS, CR	0.7364	0.5691	4.275	20416
7	lnD, lnCW, SDI, aCI, CS, CR, SC	0.7322	0.5739	4.225	20314
8	lnD, lnCW, SDI, aCI, CS, CR, SC, Hegyi	0.7297	0.5768	4.194	20253
9	lnD, lnCW, SDI, aCI, CS, CR, SC, Hegyi, D2	0.7285	0.5782	4.180	20226

. With the gradual introduction of variables, the evaluation indexes RMSE, TRE, and AIC all decreased, and the decreasing rate slowed down when the number of independent variables exceeded 5. R² increased with the increase of independent variables, and the increase rate also slowed down when the number of independent variables exceeded 5. When all the 9 independent variables were taken into account, the fitting effect of the basic model was the best, with a maximum R² of 0.5782 and a minimum RMSE of 0.7285. The optimal basic model is as follows (Equation (13)):

$$\ln (D_2^2-D_{}^2)=a_1+a_2\ln D+a_3{D}^2+a_4\ln CW+a_{5}SDI+a_{6}aCI+a_{7}CS+a_{8}CR+a_{9}SC+a_{10}Hegyi$$

(13)

where D and D₂ are the DBH values from the first and second measurements, respectively, a₁, a₂, a₃.... a₁₀ are parameters.

Table 3. Fitting results of the optimal basic model for individual-tree diameter growth.

Variables	Parameters	t Value	p Value	Variance Inflation Factor
Intercept	2.2350	16.356	<0.01	-
lnD	0.9899	29.439	<0.01	5.914
lnCW	0.4898	17.196	<0.01	2.404
SDI	−0.0014	−16.731	<0.01	1.031
aCI	−1.2310	−15.521	<0.01	2.092
CS	−0.3823	−10.878	<0.01	1.609
CR	0.4880	10.858	<0.01	1.083
SC	−0.4778	−10.559	<0.01	1.288
Hegyi	−0.0170	−7.553	<0.01	1.491
D2	−0.0002	−5.395	<0.01	3.384

shows that all the selected variables have a significant impact on the basic model, and their variance inflation factors are less than 6, so there is no multicollinearity problem [43]. In the aspect of tree size factor, the coefficient of the logarithm of D (lnD) is positive, while the coefficient of D² is negative, which indicates that the diameter growth of tree increases with the increase of diameter, and to a certain extent, its speed will decrease, showing a parabolic form, which is in line with the law of tree growth. The coefficients of CW and CR are both positive, indicating that the larger the crown width is, the higher the ratio of crown length is, the faster the growth of tree diameter is. While the coefficient of CS is negative, reflecting that the larger the ratio of height to DBH is, the smaller the growth of tree diameter is, and trees allocate more resources to the growth of tree height. The coefficients of SDI, aCI, SC and Hegyi are all negative in reflecting the competition factors around trees, which indicates that the competition among natural secondary forest trees of Q. mongolica already exists and inhibited the diameter growth of trees.

3.2. Mixed Effect Model of Diameter Growth

The model (13) has 10 parameters (a₁–a₁₀), therefore, there will be 58,236 different combinations by introducing random effects of tree species and sample plots into the parameters, and each combination forms a mixed effect model of diameter growth. In order to find the best combination of random effects of tree species and sample plots, the covariance matrix [47] of random effects of tree species and sample plots was set as diagonal matrix due to the convergence of mixed effects model.

According to the evaluation indexes of mixed effect model (see

Table 4. Mixed effect model evaluation indexes of different combinations of random effects of tree species and sample plots.

Species-Level Random Effects	Plot-Level Random Effects	Structure of Ψ	RMSE	R2	TRE	AIC
a1	a1	DM	0.6600	0.6539	3.407	18,829
a1, a2	a1, a2	DM	0.6449	0.6695	3.249	18,530
a1, a2, a3	a1, a2, a3	DM	0.6389	0.6756	3.187	18,437
a1, a4	a1, a4	DM	0.6512	0.6630	3.315	18,731
a1, a5	a1, a5	DM	0.6600	0.6539	3.407	18,833
a1, a6	a1, a6	DM	0.6527	0.6614	3.330	18,686
a1, a7	a1, a7	DM	0.6519	0.6622	3.323	18,717
a1, a8	a1, a8	DM	0.6567	0.6573	3.373	18,788
a1, a9	a1, a9	DM	0.6600	0.6539	3.407	18,832
a1, a10	a1, a10	DM	0.6576	0.6563	3.382	18,792
a1, a2, a4	a1, a2, a4	DM	0.6400	0.6745	3.198	18,458
a1, a2, a6	a1, a2, a6	DM	0.6430	0.6714	3.229	18,491
a1, a2, a7	a1, a2, a7	DM	0.6390	0.6755	3.188	18,455
a1, a2, a3, a4	a1, a2, a3, a4	DM	0.6349	0.6796	3.147	18,381
a1, a2, a3, a7	a1, a2, a3, a7	DM	0.6341	0.6805	3.138	18,377
a1, a2, a3, a7	a1, a2, a3	DM	0.6368	0.6778	3.165	18,401
a1, a2, a3, a7	a1, a2, a3, a7	UN	---	---	---	---
a1, a2, a3, a7	a1, a2, a3, a7	CS	0.6568	0.6571	3.374	18,952

Note: DM means diagonal matrix, UN means unstructured matrix, CS means compound symmetric matrix, same below.

), combined with the variance of random effect of tree species and sample plot on different parameters, the greater the variance, the greater the impact of random effect on this parameter. It is determined that the random effects of tree species and sample plot are added to parameters a₁, a₂, a₃ and a₇, and the evaluation index of diameter growth mixed effect model is the best. Then, considering the structure of variance covariance matrix of random effects of tree species and sample plot, the likelihood ratio test [48,49] is used to avoid over fitting. Finally, it is found that when the variance and covariance matrices of random effects of selected tree species and sample plots are diagonal matrices, each evaluation index of the model is the best.

Therefore, the two-level mixed effect model of individual tree diameter growth of Q. mongolica secondary forest was finally determined as (Equation (14)):

$$\begin{array}{l}\ln {(D_2^2-D_{}^2)}_{ijk}=2.098+u_{1i}+v_{1ij}+(1.081+u_{2i}+v_{2ij})\ln D_{ijk}+(-0.00021+u_{3i}+v_{3ij})D_{ijk}^2+0.5263\ln C{W}_{ijk}\\ -0.00144SD{I}_{ijk}-1.169aC{I}_{ijk}+(-0.2170+u_{4i}+v_{4ij})C{S}_{ijk}+0.2354C{R}_{ijk}-1.292S{C}_{ijk}-0.00981Hegy{i}_{ijk}+{\epsilon }_{ijk}\end{array}$$

(14)

$$\begin{array}{c}{u}_i=(u_{1i},u_{2i},u_{3i},u_{4i}{)}^{\prime }~N\left\{{(0,0,0,0)}^{\prime },{\Psi }_1=diag(1.057,0.4386,0.000519,0.2608)\right\}\\ v_{ij}=(v_{1ij},v_{2ij},v_{3ij},v_{4ij}{)}^{\prime }~N\left\{{(0,0,0,0)}^{\prime },{\Psi }_2=diag(0.1675,0.0518,0.000194,0.1988)\right\}\\ {\epsilon }_{ijk}\sim N(0,0.4096)\end{array}$$

(15)

See model (13) for the meaning of variables in the formula.

The parameters of model (14) were significant (p < 0.001). Compared with the model (13), each evaluation index was significantly improved, in which RMSE decreased by 14.89%, R² increased by 17.69%, TRE and AIC decreased by 24.93% and 9.14%, respectively. It is shown that the introduction of random effects of tree species and sample plots greatly improves the fitting accuracy of the model, and the model (14) could be used to describe the diameter growth of individual tree in Q. mongolica secondary forest.

3.3. Model Evaluation

Based on the cross-validation method, the prediction evaluation of the three prediction methods (average response, random effect prediction at tree species level, and sample plot level) of model (13) and model (14) are shown in

Table 5. The results of evaluation metrics for model (13) and mixed random effects model (14) using cross-validation.

Model		e	RMSE	R2	TRE
Model (13)		0.0073	0.7375	0.5678	4.299
Model (14)	Mean response	0.0671	0.7676	0.5317	4.827
	Species level	0.0338	0.6948	0.6165	3.852
	Plot level	0.0454	0.6907	0.6210	3.831

. The prediction results of the model (14) based on different levels of random effects are better than the prediction of its mean response and the prediction of the model (13), with smaller RMSE and TRE and larger R². In particular, the prediction accuracy of the model (14) based on the sample plot level nested in the tree species level was the highest. Compared with the average response and tree species level prediction of the model (13) and (14), its RMSE decreased by 6.35%, 10.02% and 0.59%, and its R² increased by 9.37%, 16.80%, and 0.73%, respectively. The results showed that the random effect of tree species level was significant, but the random effect of sample plot level was not as obvious as tree species. Diameter growth curve of Q. mongolica natural secondary forest (Figure 2) was produced by model (14). It can be seen that the diameter growth curve well passes through the observed value, which verifies the validity of the mixed effect model (14).

Therefore, the mixed effect model (14) is recommended to predict the individual tree diameter growth of Q. mongolica secondary forest.

The results showed that the variable D was the most important factor affecting the diameter growth, and the influence of the other 7 variables on the diameter growth from large to small order is as follows: CW, aCI, CS, CR, SDI, SC, and Hegyi. CW, CS, and CR were positively correlated with the diameter growth, while the other 5 variables were negatively correlated. The regression coefficient of SPI is not significant, so this index was not considered in the model construction. The results show that the mixed effect model has better prediction results, and the prediction accuracy based on the sample plot level nested in the tree species level is the highest, with RMSE of 0.6907 and R² of 0.6210. It can be seen from the distribution map of the predicted values and residuals (Figure 3) that there is no obvious heterogeneity in the distribution of residuals, which indicates that the two models have good prediction results, and model (14) is obviously better than model (13).

4. Discussion

Tree growth is affected by many factors, including physical, physiological, and ecological aspects [50]. Both biological importance and statistical fitting were considered when variables were selected in this study. If the model has no biological significance, it is just for model development and has no guiding significance for production and research [51]. This study showed that the variables significantly related to the diameter growth of Q. mongolica natural secondary forest were size and competition. The effect of tree size factor on diameter growth is higher than that of competition factor. R² of the model reached 0.5081 after introducing the first variable lnD. The introduction of the remaining eight variables only increased R² by 13.80%, indicating that initial DBH had the greatest impact on individual tree diameter growth. Many related studies also show this, but most of these related studies are based on traditional regression models (least square method), which is not suitable for data with complex structure [34,47,52,53]. Meanwhile, it is found that the basal area of the stand seems to be more important than the relative size of the tree in other studies [50,54,55]. CW (Figure 4a) has the greatest effect on diameter growth in addition to diameter, which may be because there is a close positive correlation between crown and diameter growth. The crown characteristics were also included in some individual tree diameter growth models [56,57].

In this study, the effect of site condition on diameter growth is not significant, it may because the 12 sample plots are in the similar site conditions, where the microenvironment does not change much. The 12 plots were all in the west direction, with a gradient between 6 and 10 degrees, and were all in the middle slope (except for the No.5 plot, which was uphill) (see Table 1). Another reason why site conditions have little effect on diameter growth may be that Q. mongolicus, which is mainly distributed in arid and poor soil in China, does not have a high demand for site conditions. The above results are consistent with the conclusions of relevant scholars [58,59]. Related studies also show that site quality has a greater impact on tree height growth, while competition factor has a greater impact on diameter, which is also verified by this study [58]. With the increase of independent variables and the introduction of random effects, the coefficient of determination of the model has been improved, likely due to the variety of species in the stands. Many related studies have also improved the model performance by introducing multiple predictors and random effects [60,61].

The independent variables in the model have different influences on diameter growth, among which the coefficients of lnD, CW and CR are positive. It can be seen that the higher the values of these three items are, the more vigorous the tree is and the faster the growth rate is. However, the coefficient of D² is negative, which indicates that the growth rate of trees will decrease when the diameter grows to a certain extent. In this study, CS is regarded as a variable representing the tree size, which can also be treated as a variable representing competition. The coefficients of CS and four competition variables SDI, aCI, SC and Hegyi are all negative. Many related studies also show that diameter growth decreases with the increase of competition [47,53,61,62]. Both stand density and competition have significant effects on diameter growth, as trees compete for limited resources such as space, nutrients, water, and light requirements during growth. In case of excessive competition, trees may grow taller, but with smaller crown and smaller diameter. Especially for the trees in the sub-forest layer, more energy will be allocated for tree height growth in order to obtain more lights [63,64]. As for the dominant trees, because the canopy can receive enough sunlight, more energy will be distributed to the diameter growth, and the stability is better than that of the compressed trees, which also explains the positive coefficient of lnD [65]. This study is based on the observation data from two times surveys of 12 sample plots in a forest farm. Long term and large-scale observation data are needed in order to predict the growth dynamics of Quercus mongolica secondary forest more accurately.