Xiao Zhou · Liyong Fu · Ram P. Sharma ·Peng He · Yuancai Lei · Jinping Guo

Abstract Tree mortality models play an important role in predicting tree growth and yield, but existing mortality models for Larix gmelinii subsp. principis-rupprechtii,an important species used for regeneration and af forestation in northern China, have overlooked potential regional inf luences on tree mortality. This study used data acquired from 102 temporary sample plots (TSPs) in natural stands of Prince Rupprecht larch in the state-owned Guandi Mountain Forest ( n = 67) and state-owned Boqiang Forest ( n = 35)in northern China. To model stand-level tree mortality, we compared seven model forms of county data. Three continuous (dominant height, plot mean diameter, and basal area per hectare) and one dummy variable with two levels(region) were used as f ixed ef fects variables. Tree morality variations caused by forest blocks were accounted for using forest blocks as a random ef fect in selected models. Results showed that tree mortality signif icantly positively correlated with stand basal area and dominant height, but negatively correlated with stand mean diameter. Incorporating both the dummy variables and random ef fects into the tree mortality models signif icantly increased the f itting improvements,and Hurdle Poisson mixed-ef fects model showed the most attractive f it statistics (largest R 2 and smallest RMSE) when employing leave-one-out cross-validation. These mixedef fects dummy variable models will be useful for accurately predicting Larix tree mortality in dif ferent regions.

Keywords Base models · Regional mortality models ·Mixed-ef fects modeling · Model validation · Forest management

Introduction

Larix gmeliniisubsp.principis-rupprechtii(Mayr) A. E.Murray (Pinaceae) is the main af forestation tree species in the mountains of northern China (Fu 2017) because of its faster growth, excellent wood materials, stronger resistance to bad weather and wind, and contributions to soil conservation. It is thus important in mountainous regions. In recent years, however, large-scale wilting has been found inLarixforests in these regions, although the death rate ofLarixspecies dif fers signif icantly among regions due to dif ferences in site conditions and environmental factors (Chen and Hua 1991; Ban et al. 1997).

Predicting tree mortality is one of the important parts of forest growth and yield models (Clutter and Jones 1980;Knoebel and Burkhart 1986). Information on tree mortality and potential causes are very important for understanding forest dynamics (Das and Nathan 2015) because mortality may strongly inf luence future stand status (Bircher et al.2015). Mortality has a signif icant inf luence on the prediction accuracy of changes of global stand structures (Dietze and Jaclyn 2014). Because we lack a comprehensive understanding of tree mortality is incomprehensive, as tree mortality is rarely observed and its cause is unclear (Das and Nathan 2015; Vanoni et al. 2016). As a result, tree mortality is still the most dif ficult part to incorporate in models to predict forest growth and harvest (Hamilton and Edwards 1976).

Tree mortality results from the combined ef fect of environmental factors, stand factors and genetic characteristics of tree and thus dif fers among stands; trees most often gradually decline in vitality until they die. Tree mortality is highly complex, with a multifactor synergism and has considerable degree of the randomness; therefore, the underlying mechanisms are dif ficult to elucidate (Sala et al. 2010),limiting modeling capability (Galbraith et al. 2010; Adams et al. 2013).

The many factors and their interactions of the factors af fecting mortality at the same time in the same place make it dif ficult to describe variability of the observed mortality data using traditional modelling approaches such as ordinary least square regression. Alternatively, mixed-ef fects modelling makes it possible to describe mortality more ef fectively than the traditional modeling approach (Zhang et al. 2014,2017b). Thus, mixed-ef fects tree mortality models should be developed to improve predictions of forest damage at stand levels.

Stand level mortality is a count variable because there are often no dead trees in a stand. The least squares method implicitly presumes that the data are Gaussian distributed with constant variances or at least satisfy Gauss–Markov assumptions. If the least squares method is applied to data with a large proportion of zero counts, the estimated results would be biased. Thus, linear models are not appropriate to describe mortality. A generalized linear model is often used to describe the variability of tree mortality, in which a dependent variable follows an exponential distribution,which may not be appropriate for potential mortality patterns. Other probability distributions such as negative binomial distribution, binomial distribution and Poisson distribution are commonly used for mortality modeling with better results (Zhang et al. 2014).

When sample plots have little or no dead trees, a large amount of zero data may be possible, and the structure of this data is discrete (Eid and Tuhus 2001). The Poisson model is often used for counting; however, the Poisson regression must have equality of the mean and the variance, so the negative binomial model is sometimes used for counting (Rashid 2016; Zhang et al. 2017a). However,in practical situations, some data are too discrete, and the negative binomial model is not suitable. Sometimes, if the model is implemented, interpretation of the results can distort (Ping et al. 2008). In this situation, researchers have applied the zero-inf lated model and Hurdle model to f it the mortality data because these methods can ef fectively solve the heterogeneity problems of the data (Hu et al. 2011; Yang 2014; Fang et al. 2016). Bayesian estimation methods have also been used (Alspach and Sorenson 1972) for mortality modeling. In forestry, counting models are mainly used for counting forest f ire incidents (Kwak et al. 2012; Xiao et al.2015; Susaeta et al. 2016) and rarely applied to tree mortality modeling (Af fleck 2006; Li et al. 2019; Zhang et al.2014). A two-step method has also been used to build standlevel mortality models (Woollons 1998; Eid and Oyen 2003)and involves f itting a logistic function to tree mortality data,and stand-level mortality is obtained by summing the number of dead trees in the stand. Because detailed information on individual trees is needed, with the chance of errors accumulating and thus reducing the accuracy of the stand-level mortality information.

Considering all these issues, here we used seven commonly used counting model (Poisson model and negative binomial model, zero-inf lated Poisson model, zero-inf lated negative binomial model, Hurdle Poisson model, Hurdle negative binomial model, and logistic regression model) to f it the tree mortality data. Considering dif ferent candidate models to f it data provides a good opportunity to select the most suitable model according to data patterns. The bestf itted model was then selected to describe the phenomenon of the regional random death ofLarix. The presented mortality models will be useful for estimating comprehensive growth processes forLarixforests in northern China for developing more ef fective silvicutural strategies and forest management plans.

Materials and methods

1Data collection

We established 102 temporary sample plots (TSPs) in stateownedLarixforests in Shanxi Province, China to collect mortality data (Fig. 1): 67 in the Guandi Mountain Forest and 35 in the Boqiang Forest. These 102 TSPs did not have any have obvious damage due to disease and pest. Each TSP was square-shaped and 0.04 ha. The TSPs were selected to provide representative information for a variety of stand structures and densities, tree heights and ages, and site productivity. Data was collected from July through September in 2015. For each stand structure, stand origin was recorded,stand age was determined, and canopy density, and height ofLarixwith DBH larger than 5 cm were measured. All 102 TSPs originated from natural forests. Tree height was measured with an ultrasonic altimeter; the crown was measured in four directions using a hand-held laser range f inder;the age of each dominant tree was determined by counting rings obtained from cores drilled at breast height the dominant height of the stand was obtained as an average of the f ive tallest trees in each quadrat within the sample plot. Summary statistics are presented in Table 1. The climate in the studied area is temperate continental. In Guandi Mountain Forest, mean annual temperate ranges from 3 °C to 7 °C,and mean annual rainfall is 822.6 mm. In Boqiang Forest,mean annual temperate ranges from ‒1 °C to 8 °C, and mean annual rainfall is about 400 mm.

Fig. 1 Study area showing the sample plot locations

Methods

We mainly choose stand factors to assess their af fected on tree mortality. Stand factors include stand density, competition index, stand productivity, stand structure, etc. In most cases, these factors may be considered simultaneously or several of them are considered (Af fleck 2006; Zhang et al.2014; Das and Nathan 2015). Based on the research data, the main factors af fecting stand level mortality were calculated,including stand density, stand mean diameter, stand dominant height, basal area per hectare, relative spacing index.Tree mortality patterns are shown in Fig. 2.

Table 1 Summary of stand-level variables; D: stand mean diameter;N: stand density (number of stems per hectare), DH: stand dominant height; S: basal area per hectare; and RSI: relative spacing index;MC: stand-level mortality in a sample plot; SD: standard deviation

Fig. 2 Tree mortality distribution patterns by number of dead trees per sample plot

Variable selection

Stand variables characterized by a meaningful biological explanation were selected as predictor variables in the tree mortality models. The dominant height (DH), which describes the combined ef fects of stand development and site productivity calculated and evaluated its potential contribution to the tree mortality models. Similarly, sample plot mean diameter (D), number of trees (N) and basal area per hectare (S), and relative spacing index (RSI), which were assumed to describe stand density and competition, were also evaluated for their potential contributions to the tree mortality variations. Multicollinearity among the independent variables was verif ied with the variance inf lation factor(VIF). According to a common rule-of-thumb, multicollinearity among variables was considered to occur when VIF > 5(Akinwande et al. 2015). Thus, the variance inf lation factor (VIF) was used to examine whether variables would be signif icantly correlated with each other, and variables with VIF < 5 were retained in our f inal models. We retained only three stand-level predictor variables in our tree mortality models, and they are DH,D, andS.

Model development

We considered seven commonly used versatile functions to develop the tree mortality models, such as Poisson model and negative binomial model (NB), which refer to as the standard function, zero-inf lated Poisson model (ZIP), zeroinf lated negative binomial model (ZINB), Hurdle Poisson model (HP), Hurdle negative binomial model (HNB), and logistic regression model. We expanded each of these functions through the inclusion of important stand-level variables (D, DH,S), random component and dummy variable.More details of the expanded models were given in Table 2.

When a dummy variable describing regional variations in tree mortality was added to parameterβ1in all the seven models, dummy variable tree mortality models were formed(Table 3).

We formulated the tree mortality models using each of the seven base models by incorporating dummy variable describing regional mortality variation and random ef fects accounting for forest block ef fects in the tree mortality models. The mixed-ef fects tree mortality models with dummy variables we formulated are given in Table 4 .

In all the models in Table 4, the vectors of errors and block-level random ef fects (u i1,u i2) are def ined byζ i～N(0,R) andμ i～N(0,D), respectively, meaning that error vector is assumed to have a normal distribution with zero mean and within-block variance–covariance matrixR i, def ined by Eq. 22.

In this study, all parameter vectors can be estimated through the maximum likelihood method. Parameter estimation was implemented using the glmmTMB package (Brooks et al. 2017) in R 3.6.3 (R Core Team 2020).

Model selection and goodness of f it

Various statistical indicators were used to compare the f itting performance of the candidate models presented above. Tocontrast the goodness of f its between these models, coef ficient of determination (R2 ), mean residual error (MD), total relative error (TRE), and root mean square error (RMSE)were used. The expressions of these indicators are given below:

Table 2 Summarized forms of mortality functions, which were expanded through inclusion of three predictor variables

Table 3 Tree mortality models with dummy variable

Using MD, TRE, RMSE andR2 alone does not ensure whether the models f itted data optimally. The validity of the tree mortality models developed from the seven dif ferent base models can be evaluated using an independent data set.However, such a validation procedure was not feasible in this study because of the limited availability of data. Instead,the predictive performance of the tree mortality models was evaluated using the leave-one-out cross-validation (LOOCV)approach (Nord-Larsen et al. 2009; Timilsina and Staudhammer 2013).

Results

Basic tree mortality models

The parameter estimates and f it statistics of all the seven candidate models using three stand-level variables (DH,D,S) as predictors, which we have def ined as basic tree mortality models, are presented in Table 5, and their models in Table 2.

All the parameter estimates for each base model were signif icant at the 0.05 level. Model 3 and Model 5 showed better-f it statistics compared to the other models. Model 1 and Model 2 were inferior compared to the other models.The complex models provided better f its than the simpler ones did. These results also conf irmed the superiority of the complex models or discrete data. Model 3 and Model 5 provided better f its than all other models did, suggesting that they were the best suited to the data structure and the sample plots with no mortality data (zero data).

Tree mortality models with dummy variable

When a dummy variable describing regional variations in tree mortality was added to parameterβ1in the seven models, f it statistics obtained were substantially better than those of their base model counterparts (Table 2) (see model forms in Table 3).

The parameter estimates of the dummy variable and all other parameters of each tree mortality model were signif icant (p< 0.05), except forβ4 andβ5 in the NB, ZINB and HNB models (Table 6). Model 10 and Model 12 provided a better f it than all the other models, with the greatestR2 and smallest RMSE and TRE. Model 9 was inferior to the other models, with the smallestR2 and greatest RMSE and TRE.

Table 4 Mixed-ef fects models with dummy variable and random ef fects

Table 5 Parameter estimates and f it statistics of all the base models. R 2 : coef ficient of determination; RMSE: root mean square error; TRE: total relative error

Table 6 Parameter estimates and f it statistics of tree mortality models with dummy variable included

Mixed-ef fects models with dummy variable and random ef fect

We formulated the tree mortality models using each of the seven base models by incorporating a dummy variable describing regional mortality variation and random ef fects accounting for forest block ef fects into the tree mortality models (see model forms in Table 4).

Except for Model 18, which did not converge with global minimum, the f it statistics of all the other mixedeffects dummy variable models were significantly

improved, and all the parameter estimates of each model were signif icant (Table 7). Model 15 f itted better than all the other models, with the greatestR2 , and the smallest RMSE and TRE. Model 16 showed an inferior f itting to other models.

Model evaluation with LOOCV

We used only those predictor variables in the tree mortality models, which had VIF < 5, to insure no collinearity occurred among them. We used the selected variables inall model types: basic models, dummy variable models,and mixed-ef fects dummy variable models. We evaluated all these model types using the LOOCV approach. The prediction improvement was substantial through adding the dummy variable and random ef fects to the basic models (Table 8). For Model 19,R2 is the largest and RMSE is 7.1% lower than that of Model 15. Among the basic models(Eqs. 1–7), Model 3 and Model 5 had the most attractive prediction statistics. Model 10 and Model 12, which are the dummy variable models and the mixed-ef fects dummy variable model, Model 19, appeared to be the best in their prediction performance.

Table 7 Parameter estimates and f it statistics of the mixed-ef fects models with dummy variable and random ef fects included

Table 8 Prediction statistics obtained with the leave-one-out crossvalidation (LOOCV) approach

When we compared the observed and predicted tree mortality distribution patterns (Fig. 3), base Model 5, dummy variable Model 12, and mixed-ef fects Model 19 showed better f itting ef fects.

Discussion

In this study, seven different mortality functions were considered for fitting the mortality data collected from two dif ferent regions of northern China, and their f itting performances were evaluated using common statistical measures.

Sample plot mean diameter (D) and stand basal area (S),which ref lect the stand diameter growth, may describe the morality caused by stand density and competition. Dominant height (DH) may ref lect the combined ef fects of site quality and stand development on the tree mortality.

Stand variablesSandDcan be calculated simply and accurately using diameter at breast height, which was the most reliably measurable variable in f ield survey data. In our models, tree mortality is signif icantly related toS, DH andDin the sample plot. VariableShad a positive correlation with tree mortality, indicating thatSraised the tree mortality rate, which may be due to resource limitations in the stand. An increase in basal area per hectare may cause crowding, and the intense competition may increase the mortality rate in the forest (Dieguezaranda et al. 2005; Wiegand et al. 2006; Moustakas et al. 2008; Zhang et al. 2015). DH was also positively correlated with tree mortality, and with larger DH, tree mortality increased. Site conditions dif fer among regions and may be responsible for dif ferences in tree mortality between the two regions with trees of the same dominant height.

Conversely, the ef fect ofDon tree mortality was negative;that is, as the stand mean diameter became smaller, the tree mortality increased. This result indicates that tree mortality was more likely in forests with many small trees compared to forests with larger trees (Juknys et al. 2006; Larson and Franklin 2010).

When the dummy variable accounting for mortality variations due to regional dif ferences was added, the model f it statistics slightly improved. However, basal area per hectare(S) appeared insignif icant in the dummy variable models(negative binomial model (NB), zero-inf lated negative binomial model (ZINB), and Hurdle negative binomial model(HNB)). The reason may be due dif ferences in regions, such as Guandi Mountain and Wutai Mountain. In the dif ferent regions,Dwould be signif icantly dif ferent due to dif ferent sites in the two regions. However, the models had higher f itting accuracy when tree mortality models including the dummy variables and random ef fects were considered based on the basic models. Because the random ef fects were added asDand the intercept, the dif ferences wre explained as the ef fects of stand mean diameter in dif ferent blocks.

For discrete data, the Poisson model and the negative binomial model have poor prediction accuracy in the basic model, while the zero-inf lated Poisson model (ZIP) and Hurdle Poisson model (HP) had unique advantages for prediction. The prediction accuracy of the zero-inf lated negative binomial model is lower than that of the zero-inf lated Poisson model, which may be due to the numerous zero data that would not be in the expansion state (Long and Freese 2006).

Fig. 3 Observed and predicted distributions of tree mortality for Larix gmelinii subsp. principis- rupprechtii

Compared to the basic models, parameters of dummy variable models increased to a certain extent, and models index have improved. The dummy variable models can well integrate dif ferent areas ofLarixstands, improve the accuracy of the tree mortality model, and expand the compatibility of the model.

If the logistic model is based on the maximum likelihood estimation, it can only estimate the probability of the dependent variable. Thus, it is not appropriate to use maximum likelihood estimation in our study. Any of the maximum likelihood-based criteria for model selection, such as the Akaike information criterion (Strawderman et al. 2000)cannot be applied to the logistic model. Alternatively, the leave-one-out cross-validation (LOOCV) was applied to evaluate prediction performance of the models.

The Poisson model performed well in principle with tree mortality, but was unable to account for the large zero fraction; the zero-inf lated negative binomial, zero-inf lated Poisson, and Hurdle negative binomial models overestimated the count part, but underestimated the zero part, resulting in a low prediction accuracy. Because of the complexity of tree mortality, it is dif ficult to interpret the f itted mortality function when the zero part was added. However, when the dummy variable (region) and random ef fects (block) were included into each of the seven base models, the prediction accuracy shown by LOOCV signif icantly improved, suggesting that there were signif icant variations in tree mortality caused by regional conditions and subject (forest block).This result justif ies applying the mixed-ef fects dummy variable modeling approach in our study.

According to the f ield survey, altitude dif ferenced inthe block is large, which is expected to greatly inf luence mortality ofLarixspecies. However, we did not include any site variables such as altitude, aspect and slope into our models.In the future, these variables need to be considered in mortality modeling. Similarly, climate change also contributes to tree mortality on a large scale (Mantgem and Stephenson 2007; Kurz et al. 2008; Allen et al. 2010; Yang 2014; Hartmann et al. 2018), so climatic factors also need to incorporated into tree mortality models.

Conclusions

Fitting and comparison of the seven basic models through incorporation of the dummy variable describing regional effects and random components describing the forest block ef fects on the tree mortality led us to the following conclusions:

(1) The models f itted with the dummy variable and random ef fects signif icantly improved f it statistics and prediction statistics compared with the basic models.

(2) Among the various model formulations (basic, dummy,mixed models), the random ef fects for the Hurdle Poisson model described the largest variation in the tree mortality.

(3) Tree mortality was signif icantly positively correlated with stand basal area and stand dominant height, but negatively correlated with sample plot mean diameter.

(4) The models only considered mortality that was caused by competition; however, the impact of different regional climate scenarios on tree death is not yet clear and needs to be studied in the future.

AcknowledgementsWe thank Dr. Guangshuang Duan for his kind help on the seven methods with R.

Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source,provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.