Developing Tree Mortality Models Using Bayesian Modeling Approach

: The forest mortality models developed so far have ignored the effects of spatial correlations and climate, which lead to the substantial bias in the mortality prediction. This study thus developed the tree mortality models for Prince Rupprecht larch ( Larix gmelinii subsp. principis-rupprechtii ), one of the most important tree species in northern China, by taking those effects into account. In addition to these factors, our models include both the tree—and stand—level variables, the information of which was collated from the temporary sample plots laid out across the larch forests. We applied the Bayesian modeling, which is the novel approach to build the multi-level tree mortality models. We compared the performance of the models constructed through the combination of selected predictor variables and explored their corresponding effects on the individual tree mortality. The models precisely predicted mortality at the three ecological scales (individual, stand, and region). The model at the levels of both the sample plot and stand with different site condition (block) outperformed the other model forms (model at block level alone and ﬁxed effects model), describing signiﬁcantly larger mortality variations, and accounted for multiple sources of the unobserved heterogeneities. Results showed that the sum of the squared diameter was larger than the estimated diameter, and the mean annual precipitation signiﬁcantly positively correlated with tree mortality, while the ratio of the diameter to the average of the squared diameter, the stand arithmetic mean diameter, and the mean of the difference of temperature was signiﬁcantly negatively correlated. Our results will have signiﬁcant implications in identifying various factors, including climate, that could have large inﬂuence on tree mortality and precisely predict tree mortality at different scales.


Introduction
Forest stands change over the years, as tree growth and mortality phenomena occur. Tree mortality, one of the main components of forest succession and dynamics, is important for the maintenance of biological and structural diversity in forest ecosystems [1]. Quantitative analysis of tree mortality is very important for understanding stand structure and dynamics. Tree mortality may be quantified through modeling that produces tree mortality models, and these are used as fundamental tools for informed decision making in forestry [2][3][4][5].
Forest stand growth is largely affected by many factors, which may be endogenous and exogenous to a stand itself. The endogenous factors are related to individual tree and stand in two main situations. Firstly, it is fully consistent with the mathematical logic, while classic statistical methods are only logical concerning making probabilistic statements about the long-term average from the hypothetical replicates of data, but not hypotheses. Secondly, relevant prior knowledge on data can be incorporated naturally into the Bayesian analysis, whereas classical methods ignore the prior knowledge other than the sample data.
Prince Rupprecht larch (Larix gmelinii subsp. principis-rupprechtii.) forests occupy approximately 65% of the forested lands in northern China, dominating the forest ecosystems in the area. Many studies have demonstrated that larch forests play critical roles in carbon storage and carbon cycling in the region [43,44]. Larch is climate sensitive, and the temperature has increased in recent years, large-scale wilting is found in larch forests in the region, although its death rate differs significantly among the regions due to differences in site conditions and environmental factors [45,46]. The climatic and ecological benefits of the intact larch forests are potentially threatened by increasing tree mortality because of changes of the endogenous and exogenous factors [47][48][49].
To address the above-mentioned issues, especially climatic effects and spatial correlations on tree mortality, this study chose the potential predictor variables at three levels: Tree-level variables, stand-level variables, and climate variables to develop the tree mortality models for larch species using the Bayesian modeling approach. The specific objectives of this study were to: (1) Select the best combination of predictor variables that could be used to develop the precise tree mortality models using the Bayesian modeling, (2) develop tree-based mortality models and evaluate their performance, and (3) evaluate the effects of the combined set of the predictor variables on tree mortality. The presented results will be useful for carbon accounting in the regions and decision making in larch forest management planning in north China.

Data
The mortality data were collected from 102 temporary sample plots (TSPs) laid out across the state-owned Guandi Mountain Larch forest (67 TSPs) and the state-owned Boqiang larch forest (35 TSPs) in northern China ( Figure 1). TSPs nested into the nine different site conditions (defined as blocks), and data were collected from July through September in 2015. Each sample plot is featured with a square shape with an area of 400 m 2 . We only considered the trees with diameter at breast height (D) > 5 cm for measurement of total height, DBH, dominant height, height to live crown base, and four crown radii. Tree height was measured with the ultrasonic altimeter; crown width was measured in four directions by the handheld laser range finder, dominant height of the stand was obtained as an average of the four tallest tree heights in each TSP. Trees were grouped into diameter classes with 4 cm interval, starting from 5 cm, i.e., [5,6), [6,10), . . . , [58,62). The median value of the intervals was assumed as a diameter class value. The distribution pattern of forest mortality is shown in Figure 2.

Selection of Predictor Variables
All the predictor variables that we evaluated for their potential contributions to tree mortality were grouped into the three sets based on the information related to tree-or stand-or climate-related mortality, (1) individual predictor tree variables (I): Diameter at breast height (D), tree height, ratio of diameter to average square diameter of stand (RD) and the sum of squares of tree diameters larger than the estimated diameter (DL); (2) stand predictor variables (S): Stand density (N), stand dominant height (DH), mean quadratic diameter (SMD) and RSI; (3) Climate predictor variables (C): Sixteen bioclimatic variables computed at the spatial resolution of 1 km × 1 km grids [50]. Climate data sets  were downloaded from the Climate AP based on the longitude, latitude, and elevation of each TSP. Data acquisition interval is one year. The definitions of all the climate variables are provided in Table 1. The mean values and the amplitude according to these predictor variables are given in Table 2.

Selection of Predictor Variables
All the predictor variables that we evaluated for their potential contributions to tree mortality were grouped into the three sets based on the information related to tree-or stand-or climate-related mortality, (1) individual predictor tree variables (I): Diameter at

Individual Tree Variables
The DL was calculated as a larger value than the sum of squares of the diameter of the target tree in a stand (Equation (1)). The RD was calculated as a ratio of diameter and stand mean square diameter (Equation (2)).
where DL ijk and RD ijk are the corresponding DL and RD of the k th tree on the j th sample plot nested within the i th block, D ijk (D ijl ) is the DBH of the k th (l th ) tree on the j th sample plot nested within the i th block, and N ij is the number of observations in the j th sample plot nested within the i th block. We dealt with the multiple collinearity problems between the predictor variables using the variance inflation factor (VIF). We retained the predictor variables with VIF < 5 in our final models, for example, D 2 , 1/D and other variables with VIF > 5 were excluded. In addition, DL and RD also significantly correlated with D, which could more effectively reflect the competition intensity than D, and therefore we excluded D from our final model.

Stand Variables
We chose widely used stand-level variables, such as dominant height and relative spacing index to fit the models. The index (Equation (3)) was calculated as suggested by Wyckoff and Clark [51].
where RSI ij , N ij , and DH ij are the RSI, stand density (N, trees/ha) and dominant tree height (DH, m), respectively, of the j th sample plot nested within the i th block. The four dominant trees (four largest trees in a 400 m 2 sample plot) were identified and measured. DH is arithmetic mean of four dominant heights [52].

Climate Variables
We evaluated the potential contributions of each climate variable (Table 1) to the description of the mortality variations. Only those variables, which were uncorrelated or less correlated or had VIF < 5 and contributed significantly to the descriptions of tree mortality were retained in the final model. With this selection approach, only two predictor variables-mean of difference in temperature (DT) and mean annual precipitation (MAP)were found to significantly contribute to the mortality models, and were included as predictor variables in the extended mortality models.
As mentioned earlier, we developed the logistic mortality models through retaining only the significant predictor variables in the models. For an individual tree, it was convenient to fit the probability of survival with a binary response variable (living or dead, 1 or 0, respectively). A widely used function for precisely describing the binary data is the logistic model [13,53], which is expressed as, where P ij is the probability of survival in the j th sample plot nested in the i th block, x ij is the matrix vector of predictor variables used to fit the model, including individual, stand and climate variables, β 0 is the intercept, and β is the parameter vector. By combining predictor variables at three different scales (individual tree, stand and climate), the advantages and disadvantages of the resulting models were compared. We built the models with combinations of the three scales of variables, such as I (only individual tree variables), S (only stand variables), C (only climate variables), I + S (individual tree and stand variables), S + C (stand and climate variables), I + S + C (individual tree, stand, and climate variables), and then selected the best among these combinations using the area under receiver operating curve (AUC), which the common statistical index used for model comparison. Then, we used the best one to establish the logistic two-level models through the Bayesian modeling.

Two-Level Models
Due to the hierarchical structure of data (trees in sample plots that are nested in a block). We established both the one-and two-level models: Block (level 1) and sample plot (level 2) level models. Let i(i = 1, 2, . . . , M) represent the level 1 unit (block) and j(j = 1, 2, . . . , n i ) represent the level 2 unit (plot); thus, a two-level random intercept logistic model is expressed as: where µ i is a two-level random effect assumed to have a normal distribution with mean 0 and variance σ 2 i . v j is a two-level random effect assumed to have a normal distribution with mean 0 and variance σ 2 ij . u i and v ij represent the random effects of the i th block and j th sample plot, respectively. β 1 -β q are the fixed-effect parameters and x 1ij -x qij are q selected predictors.
To develop the models, fixed effects only, block random effects, and two-level random effects were included into the intercept, in addition to the three-level variables (individual, stand, and climate variables). Both the one-and two-level models were fitted to account for the autocorrelations of the observations in the block (level 1) with sample plot (level 2).
All the parameters of the Bayesian logistic models are estimated through the Markov chain Monte Carlo (MCMC) simulation using the R package MCMCglmm, which uses the combination of Gibbs sampling, slice sampling, and Metropolis-Hastings sampling [54]. We used "non-informative" priors for all the regression coefficients, i.e., a normal distribution with zero mean and a large variance (104). It is important to note that the posterior mean depends largely on the choice of the non-informative priors for the variance component, i.e., uniform (0, 1) and Inverse Gamma (0.001, 0.001) [53]). We set some parameters, such as the total number of iterations: 60,000 with a burn-in of 10,000 for tree mortality models. Additionally, thinning parameters were all set to 10 to reduce the autocorrelations and we chose the Inverse Gamma (0.001, 0.001) as the variance component.
The Akaike Information Criteria (AIC) was used to evaluate the fitting performance of models [55]. The smaller the AIC, the better the model fitting effect would be.

Model Evaluation
The Bayesian two-level models were compared with the Bayesian fixed-effect (nonhierarchical) models using the deviance information criterion (DIC) [55], which is calculated as: where Q Ω is some set of the model parameters and it is the effective number of parameters in the model, mean deviance (Q) is calculated over all the iterations, and it is the posterior mean of the deviance, Q(Ω) = 2Q − Qhat. Qhat is a point estimate of deviance given by Qhat = −2 log(Prob(y Ω)) . The advantage of DIC over other criteria in Bayesian model selection is that DIC can be easily calculated from samples generated by an MCMC simulation. Models with lower DIC values indicate a better fit to the data in which differences ≥5 are regarded as substantial evidence and differences ≥10 are regarded as very strong evidence in favor of the model with the lowest DIC [54]. The value of the area under the receiver operating curve (AUC) is a thresholdindependent indicator to differentiate live and dead trees. The models with larger AUC values indicate a better fit to the data [56]. In general, an AUC value of 0.5 suggests no discrimination, 0.6-0.7 indicates poor discrimination, 0.7-0.8 suggests acceptable discrimination, and 0.8-0.9 stands for an excellent discrimination between live and dead trees [8]. The threshold is calculated by the ROCR package in R [57], which is used to determine that the predicted values of the developed models were 1 or 0. If the predicted values of tree mortality were greater than or equal to the corresponding threshold, the tree mortality would be equal to 1; otherwise, it would be equal to 0.

Model Evaluation
The results showed that I, S, C, I + S, S + C, and I + S + C were the best model forms for simulating tree mortality through the Bayesian modeling with a combination of three levels of the predictor variables ( Table 3). The models including individual variables (e.g., I, I + S, I + C, and I + S + C) showed the superior fit indices (higher AUC values) compared to other models (S, C, S + C), which had smaller AUC values. Thus, we considered only four models (I, I + S, I + C, I + S + C) with the individual predictor variables to develop both the one-and two-level models.

Two-Level Mortality Model
As mentioned above, we only selected four alternatives (I, I + S, I + C, I + S + C) to include the random effects at one level and two levels into the tree mortality models. We chose the optimal threshold from the base models, the one-level model, and the two-level model. Then, we used the optimal threshold to calculate the AUC value for each model. For example, we chose a better threshold from the basic, one-level, and two-level models using individual level variables, and then we applied this threshold to the three models to calculate the AUC value. The parameter estimates, associated standard deviations, DIC, and AUC for Bayesian fixed-effect (non-hierarchical), one-level, and two-level models are presented in Table 4. The fixed-effect parameter estimates in the multilevel models were larger than those in the fixed-effect models; therefore, ignoring the random effects underestimated most of the fixed-effect parameters. Table 4. Parameter estimates of tree survival models from fixed-effect (non-hierarchical: Base, one-level, and two-level random effects) models using the Bayesian method. * Indicates the difference between parameter estimates and 0 using t-test. '***': p < 0.0001, '**': p < 0.001, and '*': p < 0.05, '.': p < 0.10. The definitions of these parameters are the same as those in Table 3.
The DIC values of the Bayesian one-level and two-level models were much smaller than the base model, which indicates that both the one-level and two-level models are superior to the base model. Compared to the base model, the two-level model showed better-fit indexes (bigger AUC and smaller DIC). Compared to the different model by using different level predictor variables, the Bayesian two-level model showed superior fit indexes (smaller DIC and bigger AUC).
For the model constructed by the combination of predictor variables, for all forms of the components, AUC values of all models were greater than 0.8; thus, all the models were considered adequate in terms of their precisions.
In accordance with the base, one-level and two-level models of I (I + S, I + C, I + S + C), tree mortality significantly negatively correlated with RD and SMD, but positively correlated with DL. Two-level models showed the largest AUC and the smallest DIC.

Climate Effects on Tree Mortality
A summary of the parameter estimates obtained in the climate sensitive mortality models is provided in Table 4. Two climatic variables DT and MAP were significantly correlated with tree mortality. In the two-level models, DT had a negative effect on the mortality of Chinese fir. While the climate effects, including DT and MAP on the tree mortality were significant (Table 4), the effects were relatively small (Figure 3).
A summary of the parameter estimates obtained in the climate sensitive mortality models is provided in Table 4. Two climatic variables DT and MAP were significantly correlated with tree mortality. In the two-level models, DT had a negative effect on the mortality of Chinese fir. While the climate effects, including DT and MAP on the tree mortality were significant (Table 4), the effects were relatively small (Figure 3).

Discussion
The Bayesian modeling method could adequately describe the phenomenon of variations at multiple clustering levels of data, such as the block or the sample plot level in the tree mortality analysis. Compared to the traditional modeling method, such as the ordinary least square regression and the maximum likelihood method, the Bayesian method has independent prior distributions to build the models [58]. Both the traditional and Bayesian methods have their own features in the modeling tree and stand characteristics. We did not intend to compare the Bayesian models against the traditional methods in details, as the latter modeling approach is not suitable for the hierarchical data structure. In particular, when the Bayesian prior is uninformative, the results of these two methods would be almost similar [59][60][61]. Thus, in this study, we chose the Bayesian modeling method, which is considered the best.
Several researchers have also reported that tree-level, stand-level, and climatic variables such as DBH, stand basal area, and mean annual temperature significantly contribute to tree mortality [10,23,62]. We also evaluated these variables; however, prediction precision of the resulting model did not significantly improve. The insignificant contributions of such variables to the mortality models may be due to the inherent collinearity among other variables that are included into the final models. Due to the un-

Discussion
The Bayesian modeling method could adequately describe the phenomenon of variations at multiple clustering levels of data, such as the block or the sample plot level in the tree mortality analysis. Compared to the traditional modeling method, such as the ordinary least square regression and the maximum likelihood method, the Bayesian method has independent prior distributions to build the models [58]. Both the traditional and Bayesian methods have their own features in the modeling tree and stand characteristics. We did not intend to compare the Bayesian models against the traditional methods in details, as the latter modeling approach is not suitable for the hierarchical data structure. In particular, when the Bayesian prior is uninformative, the results of these two methods would be almost similar [59][60][61]. Thus, in this study, we chose the Bayesian modeling method, which is considered the best.
Several researchers have also reported that tree-level, stand-level, and climatic variables such as DBH, stand basal area, and mean annual temperature significantly contribute to tree mortality [10,23,62]. We also evaluated these variables; however, prediction precision of the resulting model did not significantly improve. The insignificant contributions of such variables to the mortality models may be due to the inherent collinearity among other variables that are included into the final models. Due to the uncertainties and complexity in the tree mortality process, we lacked an effective tool or method to determine a reasonable combination of tree-level, stand-level, and climatic variables and their interactions with the stand-specific conditions for predicting tree mortality.
A ratio of the diameter to the average square diameter of the stand (RD) had a negative correlation with tree mortality and the sum of squares of tree diameters greater than the estimated diameter (DL) had a positive correlation. With smaller RD and bigger DL, the probability of tree mortality would be smaller. This may be due to smaller RD and bigger DL in the same stand, indicating that the individual tree DBH is small. As the tree mortality rate generally decreases with increasing tree size, mortality rates can be higher for the juvenile trees, decrease with increasing tree size, and start to increase again with further increase of tree size [63]. Similarly, the effect of SMD (stand quadratic mean diameter) (size index) on tree mortality was negative (i.e., greater SMD, smaller stand mortality), indicating that tree mortality is more likely in the forests with many small trees compared to the forests with larger trees [64]. This finding is also consistent with the study [65], which shows that the mortality rate of young trees decreases with the increase of tree size and begins to increase when the mature age is reached. According to the local data survey [66], the mature age of Prince Rupprecht larch is 85 years. In this study, the data we collected were from natural secondary forests intervened in by human activities in the 1970s and all the trees almost did not reach mature age. In addition, Prince Rupprecht larch is a light-loving tree species. If young trees grown under the main forest layer obtain less light resource due tree crowding and intense competition, they would die easily.
Several existing studies show the climate warming-induced drought that might be the main driver of a widespread increase of tree mortality [13,[67][68][69][70][71]. Our study confirms that DT and MAP, which were selected among the 16 bioclimatic predictor variables, have significant influences on tree mortality. A significant positive correlation with DT suggests that low mean temperature in the warmest month and high mean temperature in the coldest month increases tree mortality. Lower mortality is related to higher DT, and higher DT indicates a higher probability of the high temperature of the study area within the growing season, and hence the lower mortality rate [51,72]. This is also consistent with the findings from the Table 4 and Figure 3 that an increase in temperature increases the probability of survival (i.e., positive trend in Figure 3, and positive coefficients in Table 4).
Temperature affects tree growth and mortality by changing photosynthesis, respiration, cell division and elongation, chlorophyll synthesis, enzyme activity, water absorption, and transpiration [73]. Generally speaking, when the available water is not limited during the growing season, the increase of temperature difference physiologically reduces the probability of tree mortality. In this study, the temperature difference between day and night in the state-owned Guandi Mountain Larch forest and the state-owned Boqiang larch forest within the growing season of the Prince Rupprecht larch is small, and this means that the Prince Rupprecht larch grows in the environment of extreme high temperature for a long time, and the long-term extreme high temperature may increase the forest mortality. The main reasons for this may be that: (1) Increasing water shortage of trees and increasing drought pressure may directly or indirectly lead to tree mortality; (2) promotion of growth and reproduction of insects and pathogens attacking trees; (3) Prince Rupprecht larch trees situated close pores to prevent hydraulic failure (i.e., water column cavitation), which may lead to carbon starvation because high respiratory costs lead to depletion of carbon reserves [72]. The effect of drought caused by climate change on tree mortality is consistent with the recent drought in the withering event of the subtropical monsoon evergreen broad-leaved forest in southern China. In temperate forests in the western United States, a positive correlation was observed between the short-term change of forest loss and drought caused by climate change [19].
The significant positive correlations between the MAP and tree mortality (Table 4, Figure 3) may be due to the co-effects of temperature and precipitation on tree mortality [72]. In northern China, overall, the temperature of the wettest quarter is quite high, and there is the greatest amount of precipitation during the year [74]. However, Prince Rupprecht larch is a typically drought-resistant species in northern China [66]. Therefore, too much precipitation and too high temperature could inhibit physiological processes of the Prince Rupprecht larch and growth [75]. In addition, a large amount of precipitation also causes excessive nutrient loss by flooding [76]. These are the reasons that higher mortality is related to a greater MAP.
Compared to the mortality models built with the tree-level predictor variables, the improved method with the model built with the tree-level predictor variables could show better performance. With the smallest DIC and the biggest AUC, the Bayesian two-level model (I + S + C) was the best for predicting tree mortality of larch. This is similar to the results with the random intercept term included in the multilevel tree mortality models, which largely improved the fit index compared to the fixed-effect model [77]. Moreover, to effectively account for spatial correlations of our data, we preferred using the Bayesian modeling method. All three forms of the model (base, one-level, and two-level models), which were constructed using the Bayesian modeling approach, showed the excellent statistical indicators (AUC > 0.8). The influence of climate factors on the tree mortality are shown to be significant through modeling, and resulting mortality models can effectively assess the forest mortality under climate change.

Conclusions
We built the tree mortality models with the inclusion of tree-level, stand-level, and climatic predictor variables using the Bayesian modeling approach for Prince Rupprecht larch in northern China. The main conclusions are: (1) The best model included the predictor variables at three levels: Individual tree-and stand-, and environmental (climate)-levels in the Bayesian logistic models. (2) The Bayesian two-level model, which includes tree-level, stand-level, and climatic predictor variables, outperformed all the other forms of the models, describing larger variations of tree mortality and accounting for multiple sources of the unobserved heterogeneities. (3) Tree mortality significantly positively correlated with the sum of squares of tree diameters larger than the estimated diameter, and mean annual precipitation, but negatively correlated to the ratio of the diameter to the average square diameter of stand, the stand arithmetic mean diameter, and the mean of difference in temperature. (4) Presented mortality models will have significant implications for identifying different factors affecting tree mortality and precise prediction of the mortality. (5) With the mortality data collected from a wider distribution of the tree species of interest and advanced modeling techniques, the prediction performance of the tree mortality models may be improved, which we aim for in the future.
Author Contributions: L.X., X.C. and X.Z. collected data; L.X. and X.Z. analyzed data; L.X., X.C., R.P.S., X.Z. and J.L. wrote manuscript and contributed critically to improve the manuscript, and gave a final approval for publication. All authors have read and agreed to the published version of the manuscript.