- Email: [email protected]
Predicting tree recruitment with negative binomial mixture models
Forest Ecology and Management 270 (2012) 209–215
Contents lists available at SciVerse ScienceDirect
Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco
Predicting tree recruitment with negative binomial mixture models Xiongqing Zhang a, Yuancai Lei a,⇑, Daoxiong Cai b, Fengqiang Liu c a
Research Institute of Forest Resource Information Techniques, Chinese Academy of Forestry, Beijing 100091, PR China Experimental Center of Tropical Forestry, Chinese Academy of Forestry, Pingxiang 532600, Guangxi, PR China c Southern Swedish Forest Research Centre, Swedish University of Agricultral Sciences, 232 38 Alnarp, Sweden b
a r t i c l e
i n f o
Article history: Received 4 October 2011 Received in revised form 31 December 2011 Accepted 19 January 2012 Available online 18 February 2012 Keywords: Tree recruitment Negative binomial model Zero-inflated negative binomial model Hurdle negative binomial model Chinese pine
a b s t r a c t Tree recruitment models play an important role in simulating stand dynamic processes. Periodic tree recruitment data from permanent plots tend to be overdispersed, and frequently contain an excess of zero counts. Such data have commonly been analyzed using count data models, such as Poisson model, negative binomial models (NB), zero-inflated models, and Hurdle models. Negative binomial mixture models (zero-inflated negative binomial model, ZINB; Hurdle negative binomial model, HNB) including NB model were used in this study to predict tree recruitments of Chinese pine (Pinus tabulaeformis) in Beijing. ZINB model and HNB model were suitable for dealing with excess zero counts, for which two equations are created: one predicting whether the count occurs (logistic function) and the other predicting differences on the occurrence of the count (NB model). Based on the model comparisons, the results showed that negative binomial mixture models performed well in modeling tree recruitment, and ZINB model was the best model of negative binomial mixture models. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction Forest dynamic process includes three components: tree mortality, tree growth, and recruitment of new trees. For these decades, foresters have made many researches on the first two components (e.g., Monserud and Sterba, 1999; Zhao et al., 2004; Cao, 2006; Zhang and Lei, 2009; Subedi and Sharma, 2011, etc.). Dynamic models may neglect the simulation of the third component, tree recruitment, because of the nature of data or difficulty in modeling. Recruitment, sometimes referred to as ingrowth (Shifley et al., 1993) is usually quantified by means of number of trees or saplings that reach or exceed a specific threshold size (e.g., 1.3 m height, 5 cm DBH over bark, etc.) over a certain period as a result of different regeneration processes like establishment, growth, and mortality of saplings (Lexerød, 2005; Lexerød and Eid, 2005). Recruitment contributes significantly to volume production in multi-layered forests (Andreassen, 1994). Recruitment model is part of the forest management tools that enable prediction of the development of forest stands, and disregarding recruitment will give biased prediction of future forest growth and yield. Because of the importance of tree recruitment, many researchers have developed models to predict recruitment over a certain period, including linear, nonlinear, and Tobit models (e.g., Adams and Ek, 1974; Shifley et al., 1993; Liang et al., 2005).
⇑ Corresponding author. Tel.: +86 10 62889199; fax: +86 10 62888315. E-mail addresses: [email protected] (X. Zhang), [email protected] (Y. Lei). 0378-1127/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.foreco.2012.01.028
Considering the fact that in permanent sample plots a relatively high number of the plots have no occurrences of recruitment even over periods of several years, it means that recruitment data are bounded and characteristically exhibit varying degrees of dispersion and skewness in relation to the mean. Additionally, the data often contain an excess number of zero counts. Yet least squares method implicitly presumes that the data are Gaussian distributed with constant variance, or at least satisfy the Gauss–Markov conditions. If the method is still used to deal with the data with large proportion of zero counts, the estimated results will be biased. Alternatively, if only recruitment observations are used for model development, recruitment will be overestimated. To overcome this problem, Ferguson et al., (1986) used a twostep approach to predict recruitment in the PROGNOSIS model (Stage, 1973; Wykoff, 1986). With a two-step approach, the probability of zero recruitment in the subsequent measurement interval, given current stand conditions, is obtained from a logistic regression in the first step. Then the amount of recruitment, given that some is known to occur, is described by a recruitment function fit by least squares in the second step. The method is notable in its recognition of excess number of zero fraction as an aspect of the response distribution not well solved by standard nonlinear regression models and has been utilized with some modifications in other researches (e.g., Woollons, 1998; Fridman and Ståhl, 2001; Álvarez-González et al., 2004; Lexerød, 2005; Fortin and DeBlois, 2007). But the method implies that a large zero fraction reflects additional structure in the mean function, structure that should carry some ecological relevance if it is to be recognized.
210
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
Additionally, this method does not admit a well-defined stochastic structure that can serve as a basis for inference (Affleck, 2006). In recent years, there has been considerable interest in models for count data, such as manufacturing defects (Lambert, 1992), patent applications (Crepon and Duguet, 1997), species abundance (Welsh et al., 1996; Hall, 2000; Cunniningham and Lindenmayer, 2005), medical consultations (Gurmu, 1997), use of recreational facilities (Gurmu and Trivedi, 1996; Shonkwiler and Shaw, 1996). There are also a few studies so far have addressed this issue in forestry (e.g., Keefe, 2004; Rathbun and Fei, 2006; Eskelson et al., 2009; Flores et al., 2009; Li et al., 2011). Fortin and DeBlois (2007) applied zero-inflated Poisson model (ZIP) and zero-inflated discrete Weibull (ZIdiW) to model tree recruitment of hardwood stands. But the base condition for the positive counts in the ZIP model is that the variance equals the expectation, and to ZidiW model, it is difficult for estimating an additional vector of parameters. Negative binomial model (NB) is another basic count-data model. But, it is not suitable to deal with the data with excess of zero counts (Zeng et al., 2009). Zero-inflated negative binomial model (ZINB) and Hurdle-negative binomial model (Hurdle-NB) were common methods to analyze this type of data (Affleck, 2006; Zeileis et al., 2008). The objective of this study was to develop and compare the ZIP model, Hurdle-Poisson model, NB model, ZINB model, and HurdleNB model for predicting tree recruitment over a period of time. The modelings were restricted to include available independent-variables directly from forest inventories for forest management planning.
count data, and the probability mass function (PMF) for Poisson is characterized by as
PðY ¼ yÞ ¼
ek ky y!
ð1Þ
where y refers to random variable of recruitment count, y = 0, 1, 2. . ., and k > 0. A Poisson regression model is obtained by relating the mean k to a vector of independent variables X, by k ¼ ExpðXbÞ; where b is a vector of regression coefficients to be estimated. A character of the Poisson probability function (Eq. 1) is that the mean and the variance are equal, that is Var½YjX ¼ E½YjX ¼ k. When data are not suitable for the Poisson distribution, it is typically because of overdispersion, meaning the model’s variance exceeds the mean value. The Negative binomial regression (NB) can be considered a generalization of Poisson regression and assumes that the conditional mean k of Y is not only determined by X but also a heterogeneity component unrelated to X, which addresses the issue of overdispersion by including a dispersion parameter to accommodate the unobserved heterogeneity in the count data (MacNeil et al., 2009). The PMF for the negative binomial is given by
!h1
Cðy þ h1 Þ h1 PðY ¼ yÞ ¼ 1 1 Cðy þ 1ÞCðh Þ h þ k
k
h1 þ k
y ð2Þ
where h represents the dispersion parameter. The mean is k, as the Poisson model, but the variance is k þ hk2 , thus allowing the variance to exceed k. NB regression model is obtained by relating the mean k to a vector of independent variables, k ¼ ExpðXbÞ. The pre^ corresponds to the expectation k of the distribution. diction y
2. Data 3.2. Zero-inflated models The data, collected by the Inventory Institute of Beijing Forestry, consisted of a systematic sample of permanent plots with a 5-year re-measurement interval. The plots (squared plots), 0.067 ha each, were in Chinese pine (Pinus tabulaeformis Carr.) plantations situated on upland sites of Beijing (Fig. 1). The data consisted of measurements from 132 plots obtained between 1986 and 2001. Summary statistics of stand variables are presented in Table 1.
3. Method Tree recruitment (a threshold size of diameter at breast height for tree recruitment in China is DBH = 5 cm) was assumed to relate with location variables, site conditions and stand characteristics. The models should be biologically realistic to ensure logical behavior when applied outside the range of the modeling data. Thus, an appropriate biological interpretation was the main criterion with respect to variable selection (Lexerød, 2005). Independent variables, with an appropriate biological interpretation, were selected among those describing location (altitude Alt, slope Slo), site conditions (dominant height H), and stand characteristics (stand density N, stand arithmetic mean diameter Dm, stand basal area B, relative spacing, etc.). In this study, these independent variables above would be used to predict tree recruitment. Moreover, we analyzed the multicollinearity between independent variables through variance inflation factor (VIF) test. A common rule of thumb is that if VIF > 5 then multicollinearity exists. The VIF test was performed through DAAG package in R software.
Although NB model is better than Poisson model for dealing with dispersed data (Liu and Cela, 2008), there are some situations where a major source of overdispersion is a preponderance of zero counts. This suggests consideration of more flexible count data models such as zero-inflated models. They are mixture models combining a count distribution with a point mass at zero. In zero-inflated models there are two sources of zeros: zeros may come from both the point mass and from the count component (Lambert, 1992; Eskelson et al., 2009). Actually, two regression equations are created in zero-inflated models: one predicting whether the count occurs (logistic function) and the other predicting differences on the occurrence of the count (Poisson or NB) (Karazsia and van Dulmen, 2008). Supposing a discrete variable Y (recruitment counts) obeys to zero-inflated distribution, the PMF for zero-inflated distribution is given by:
PðY ¼ yÞ ¼
p þ ð1 pÞf ð0Þ
y¼0
ð1 pÞf ðyÞ
y>0
ð3Þ
Where p is the probability of belonging to the point mass at zero and (1 p) is the probability of belonging to the count component. Logistic regression was commonly used to fit point mass compop nent: logitðpÞ ¼ Logð1p Þ ¼ Xd, where d is a vector of parameters to be estimated. In the count component, it was the same to Poisson or NB model k ¼ ExpðXbÞ.
3.1. Negative binomial model
3.2.1. Zero-inflated Poisson model (ZIP) In count component, if Y obeys to the Poisson distribution with parameter k in Eq. (3), then we get the ZIP distribution. The PMF for ZIP is given by:
Count data models are a subset of discrete-response regression models and aim to explain the number of occurrences or counts of an event. Poisson regression is the simplest regression model for
8 < p þ ð1 pÞExpðkÞ PðY ¼ yÞ ¼ ð1 pÞExpðkÞky : y!
y¼0 y>0
ð4Þ
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
211
Fig. 1. The spatial location of Chinese pine in Beijing. Lines in the figure are the boundaries of administrative region in Beijing.
Table 1 Summary statistics of stand-level variables, including tree recruitment counts in a stand Tr, stand age A, stand dominant height H, stand density N, stand basal area B, stand arithmetic mean diameter Dm, altitude Al, slope Slo and the relative spacing Rs. The relative spacing is a function of dominant height and stand density: Rs = (10,000/ N)0.5/H. Variables
Min
Max
Mean
SD
Tr (trees/plot) A (years) H (m) N (trees/ha) B (m2/ha) Dm (cm) Rs (trees/ham) Alt (m) Slo (°)
0 12 2.5 132 0.3580 5.8 0.1712 140 5
59 55 17.4 2164.1790 28.3105 16.1260 3.4512 1400 40
7.9242 28.3561 6.6402 899.2537 8.0792 9.9097 0.6823 443.2500 25.4318
12.0061 8.5132 2.9006 453.5848 5.6711 2.3951 0.4675 266.9939 8.1082
SD, standard deviation.
3.2.2. Zero-inflated negative binomial model (ZINB) If the NB distribution is used as count component, the mixture distribution is zero-inflated negative binomial, and the PMF for ZINB is given by:
PðY ¼ yÞ ¼
8 > < p þ ð1 pÞ > : ð1 pÞ
h1 h1 þk 1
1=h
Cðyþh Þ Cðyþ1ÞCðh1 Þ
h
y¼0 1
h
1
1=h
þk
y
k h þk 1
ð5Þ y>0
When h ! 0, the ZINB is changed into ZIP model. 3.3. Hurdle models Hurdle models, originally proposed by Mullahy (1986) in the econometrics literature (e.g., Cameron and Trivedi, 1998, 2005) are another model class capable of dealing with excess zero counts. They are also called two-part models (Heilbron, 1994), thus similar to zero-inflated models. Hurdle models consider count outcome generated by two systematically different statistical processes, a binomial distribution determining if a count outcome is zero or nonzero and a truncated-at-zero distribution for count data governing all positive counts conditional on count component. However, they are slightly different from zero-inflated model with all zero counts from two different sources and assume that zero counts might come from a single statistical process (Liu and Cela, 2008). Similar to the zero-inflated models, logistic function often was used to model point mass at zero. About count component, Poisson model or negative models were used in this case. 3.3.1. Hurdle-Poisson model (HP) One of the most popular Hurdle model is called HP model, which is the combination of a logit regression modeling point mass at zero and a truncated Poisson regression modeling positive counts conditional on count component. Its PMF is given as:
212
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
8
y¼0 ð6Þ
y>0
3.3.2. Hurdle-negative binomial model (HNB) Similar to ZINB model, HNB model is another important Hurdle model, which is the combination of a logit regression modeling point mass at zero and a truncated negative binomial regression modeling count component. Its PMF is given as:
( PðY ¼ yÞ ¼
p
h i1 y 1 ðyþh1 Þ k ð1 pÞ CðhC1 ð1 þ hkÞh 1 kþh1 ÞCðyþ1Þ
y¼0 y>0 ð7Þ
^ yi y ffiffiffiffiffi i pr ¼ p v^ i
ð11Þ
^i is the estimated mean and v^ i is the variance of y ^i of the ith where y observation. To evaluate the adequacy of the models for predicting the overall counts of tree recruitments, an approximate v2 goodness-of-fit was calculated using observed and predicted counts (Eskelson et al., 2009):
v2 ¼
P 2 m X ½#ðyi ¼ jÞ Pðyi ¼ jÞ P Pðyi ¼ jÞ j¼1
ð12Þ
where # represents the frequency of observations yi in count class j, m represents the number of count classes, and P(y = j) is the predicted probability that an observation belongs to count class j.
3.4. Parameter estimation
5. Results and discussion
All vectors of parameters b, d, and h can be estimated through the optimization of the likelihood function, i.e., the maximum likelihood method (ML). Standard errors of model parameters were estimated with two methods. The first one, based on standard ML theory from the inverse of the k k (k is the number of parameters estimated in the model) observed Fisher information matrix (Affleck, 2006):
According to the VIF test, we found that there was no multicollinearity among stand age, stand dominated height, relative space, altitude, slope, and any two of the combinations of stand density, stand arithmetic mean diameter, and stand basal area (VIF < 5). Table 2 lists the parameter estimates and their standard deviation errors for NB model, ZIP model, ZINB model, HP model, and HNB model. In the positive count component of the models, the majority of variables were significant at the 0.05 level; and for ZINB model, parameters of dominant height and slope were significant at the 0.1 level. Dominant height reflects the site situation. Vanclay (1989) found that the number of recruits became larger on good sites. In this study, the relation between dominant height and tree recruitment was positive (i.e. the greater H, the greater the number of tree recruitment). Stand basal area is a measure of stand competition, reflecting both size, and density of trees in a stand, and as it increased, tree recruitments decreased. This was consistent with the results of previous studies (Vanclay, 1989, 1992; Lexerød, 2005). Similarly, the effect of stand arithmetic mean diameter (size index) was negative (i.e. the greater Dm, the smaller number of tree recruitment). This meant that tree recruitment was more likely in forests with many small trees compared to forests with larger trees (Lexerød, 2005). In contrast, as the altitude became larger, tree recruitments increased. And the parameter estimates of altitude for these five models were nearly same. On the other hand, as the slope increased, tree recruitments decreased. In the zero component of the models, the parameters of ZINB model were all significant at the 0.1 level (Table 2).Two model fit statistics are gi^ ; yÞ for ZINB model was 737.5109, as ven in Table 2. The 2Lðu compared to 1241.5210, 1240.2000, 744.4256, 748.6767 for ZIP model, HP model, NB model, and HNB model, respectively; AIC was 761.5109, versus 1257.5210, 1264.1680, 762.4734, 768.6767. ^ ; yÞ and AIC of the ZINB These numbers showed that the 2Lðu model were much smaller than ZIP model and HP model, and slightly smaller than NB model and HNB model. It meant that the ZINB regression provided the ‘‘best’’ model. Pearson residuals are shown in Fig. 2. All five regression models overestimated for small recruitment counts and underestimated for large counts. The residuals ranged from 2 to 6 for ZIP model and HP model, while ranged from 1 to 3 for NB model, ZINB model, and HNB model. The v2 value of ZINB model was 2.9944, as compared to 3.8499, 3.9821, 44.2936, 46.1286 for NB model, HNB model, ZIP model, and HP model, respectively, and p-value >0.05 for negative binomial mixture models, p-value <0.0001 for Poisson mixture models (For the v2 -test, classes were grouped together to obtain expected frequencies of 5 or greater). The lower the v2 statistic, the better the model fit. ZIP model and HP model accommodated the large zero fraction, but their Poisson component fit the positive count poorly. Fig. 3 showed the empirical
" ^ m ¼ va^r ½u
@ 2 Lð/; yÞ ju^ @/@/0
#1
"
n X @ 2 Lð/; yi Þ ¼ ju^ @/@/0 i¼1
#1 ð8Þ
^ yÞ is the log-likelihood of the developed model, u ^ is the where Lð/; ^ includes the vector of estik 1 vector of estimated parameters (u mated parameter b, d, and h depending on the model fitted), n is number of observations, and y is a column vector whose elements are the responses yi. The variance estimator is consistent if the model is correctly specified. But it is biased if the data are overdispersed. Thus, standard errors were also calculated from the matrix product.
^ r ¼ va^r½u ^ m va^r ½u
" # n X @Lð/; yi Þ @Lð/; yi Þ ^ m r ½u j ^ va^ u @/ @/0 i¼1
ð9Þ
The standard errors obtained with Eq. (9) are consistent in model specification and require only the unbiased ML estimating equations (Cameron and Trivedi, 1998). Estimation of parameters was implemented with pscl package in R software (Jackman, 2011). 4. Model selection and goodness-of-fit Negative binomial model, zero-inflated models and Hurdle models calibrated with the same data set can be compared through ^ ; yÞ) as well as Akaike information crilog-likelihood values (2Lðu ^ ; yÞ indicate that modterion (AIC). Smaller values of AIC and 2Lðu el is better.
^ ; yÞ þ 2k AIC ¼ 2Lðu
ð10Þ
Vuong test is considered as a better method to compare two non-nested models for count data (Vuong, 1989). The Vuong tests were performed in R software. Goodness-of-fit was also evaluated from predictive and residual diagnostics. Residuals measure the departure of predicted values from observed values of the dependent variable. The most commonly used residual is the raw resid^. For the classical linear model with normally ual: r ¼ y y ^Þ N½0; r2 , so that in large distributed homoscedastic error ðy y samples the raw residual has the desirable properties of being symmetrically distributed around zero with constant variance. However, for count data, the raw residual is heteroscedastic and asymmetric. The obvious correction for heteroscedasticity is the Pearson residual (Cameron and Trivedi, 1998):
213
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215 Table 2 Parameter estimations and fit statistics of ZIP model, HP model, NB model, ZINB model, and HNB model. Parameter
ZIP model
HP model
Estimation
Sig.
Positive count component of the model **** Intercept 5.3011 (0.1861) **** H 0.0721 (0.0184) **** B 0.0384 (0.0102) **** Dm 0.3592 (0.0282) **** Alt 0.0014 (0.0001) **** Slo 0.0233 (0.0038) Zero component of the model Intercept – – **
Estimation
Sig.
Estimation
Sig.
Estimation
Sig.
5.2921 (0.1881) 0.0681 (0.0186) 0.0372 (0.0103) 0.3562 (0.0286) 0.0014 (0.0001) 0.0231 (0.0038)
****
4.8353 (0.6504) 0.1220 (0.0545) 0.0559 (0.0210) 0.3515 (0.0808) 0.0014 (0.0004) 0.0235 (0.0137)
****
4.8954 (0.6197) 0.1321 (0.0507) 0.0567 (0.0317) 0.3631 (0.0792) 0.0014 (0.0004) 0.0228 (0.0136)
****
5.3638 (0.6455) 0.0875 (0.0525) –
****
–
–
–
–
*
–
–
**
–
–
*
–
–
*
–
–
0.7300 (0.1090) 744.4256 762.4734
****
247.7297 (140.0640) 8.6282 (5.1010) 3.9387 (2.0332) 284.6551 (164.2425) 0.0350 (0.0205) 0.7944 (0.1515) 737.5109 761.5109
N h
–
^ ; yÞ 2Lðu AIC
1241.5210 1257.5210
1240.2000 1264.1680
Rs
– *
HNB model
Sig.
1.1600 (0.4901) –
B
0.1042 (0.0502) –
ZINB model
Estimation
0.1487 (0.0761) 0.1327 (0.0598) 0.7956 (0.4708) 0.0009 (0.0005) –
H
NB model
****
****
****
****
****
**
**
****
***
**
***
*
****
****
*
–
–
*
0.1487 (0.0761) 0.1327 (0.0598) 0.7956 (0.4708) 0.0009 (0.0005) 0.7641 (0.2877) 748.6767 768.6767
*
***
*
*
*
*
****
20
30
40
50
60
0
30
40
HP
20
30
40
50
60
Observed recruitments
50
60
50
60
-2 0 2
4 6
ZINB 1 2 3 4
10
20
Observed recruitments
-1 0
10
Observed recruitments
Pearson residual
10
-2 0 2 4 6
Pearson residual
1 2 3 4
ZIP
-1
Pearson residual Pearson residual
0
0
10
20
30
40
Observed recruitments
1 2 3 -1
Pearson residual
HNB
0
10
20
30
40
50
–
0.4210 (0.0745) 0.0013 (0.0004) 0.0224 (0.0134)
Values in parentheses are standard deviation errors, and Sig. is significance. * P < 0.1. ** P < 0.05. *** P < 0.01. **** P < 0.001.
NB
*
60
Observed recruitments Fig. 2. Pearson residual plots of tree recruitments for ZIP model, HP model, NB model, ZINB model and HNB model.
***
*
**
*
*
***
214
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
Fig. 3. Comparison between the empirical distribution of observed counts (black solid line) and the empirical distribution of predicted counts (according to the ZIP model, HP model, NB model, ZINB model, and HNB model).
distribution of the observations and the empirical distribution of the predictions of the five models. It was clear that NB model fit the observed count outcomes as well as ZINB model and HNB model, which was much better than ZIP model and HP model. A ZIP or HP will reflect data accurately when overdispersion is caused by a preponderance of zeros. If overdispersion is attributed to factors beyond the inflation of zeros, a ZINB or HNB model is more appropriate (Long and Freese, 2006). And the ZIP model or HP model accommodated the excess of zero component, but its Poisson component fit the positive component poorly (Fig. 3). Poisson is highly restrictive as the variance of the outcome is assumed to equal its mean. Count data often exhibits overdispersion. NB offers a dispersion parameter which explains the overdispersion of positive component (MacNeil et al., 2009). Therefore, NB mixture models (NB, ZINB, and HNB) performed better than Poisson mixture models (ZIP, and HP). However, we should know that the dispersion parameter h of the NB, ZINB, or HNB model is scale dependent when dealing with real data (Pielou, 1957; Eberhardt, 1967). For the Hurdle models, the zero hurdle component describes the probability of a zero count while for the zero-inflated models, the zero-inflation component predicts the probability of observing a zero count from the point mass component and count component (Zeileis et al., 2008). Overall, both models led to the same qualitative results and very similar model fits (Table 2; Figs. 2 and 3). NB model provides flexible structures for describing empirical count data. The model showed closer agreement with the tree recruitment data than the ZIP model and HP model on the basis of AIC and in terms of their ability to reproduce the observed count distribution (Table 2; Fig. 3). The improvement fit of the NB model could be ascribed to unobserved heterogeneity in tree recruitment across and with plots over time. The NB model could also be attributed to contagion, wherein the occurrence of one event augments the probability of successive events (Boswell and Patil, 1970). However it has only one mode, like the Poisson, does not deal with the extra-zero problem even though it allows for extra dispersion (Cunniningham and Lindenmayer, 2005). And ZINB model or HNB has three advantages. Firstly, by allowing separate models for the two components, tree recruitment can be modeled com-
pared with simply assuming the data is approximately Poisson or NB. Secondly, the zero component of the model and positive count component of the model of the two models are directly interpretable in terms of the empirical features of the data. Thirdly, the two models will better model the variability seen in this data. This will ensure that the information in the data is used appropriately (Barry and Welsh, 2002). Additionally, in this study, Vuong test-statistic V was 4.9617, pvalue <0.0001 for comparing HNB model with HP model, and V was 5.0648, p-value <0.0001 for comparing ZINB model with ZIP model. This showed that ZINB model or HNB model was much better than ZIP model and HP model. In the same way, for comparing the ZINB model with NB model, V was 1.6963, p-value was 0.0449, which indicated that ZINB model performed better than NB model at the significant level 0.05. We also compared ZINB model with HNB model, V was 1.4883, p-value was 0.0683. It showed that ZINB model was better than HNB model at significant level 0.1. Overall, ZINB model provided more opportunities for exploring when the probability of tree recruitment and the number of tree recruitment may be driven by different factors. Although the ZINB model is proved to have a good fit for recruitment data with an excess of zero counts, the extra zeros are not the only feature of the recruitment patterns. This typically occurs when considering a rare species, a mast species, counting recruitment of each different tree species in mixed forests, recruitment on small plots, or recruitment in short census interval for slow-growing tree species. With increasing census interval, the probability of zero occurrence gradually decreases (Yao et al., 2001).
6. Conclusion In this paper we have reviewed several count data models for analyzing tree recruitment. When the data is overdispersed, the negative binomial model is a good choice. The NB model is more flexible as it allows for the variance to be greater than mean and considers the observed heterogeneity in the model (Yaacob et al., 2010). However, the phenomenon of tree recruitment was always characterized by large number of zero counts. In this case, zero
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
negative binomial model and Hurdle negative binomial model might give a more satisfactory fit to the data. Based on the model comparisons, ZINB model performed better than other four models for modeling tree recruitments. These techniques can also be extended to other applications, such as stand mortality, snag abundance, species abundance, tree generation, etc. (e.g., Barry and Welsh, 2002; Affleck, 2006; Eerujäinen et al., 2007; Eskelson et al., 2009). Acknowledgments The authors express their appreciation to State Forestry Administration (SFA), the Ministry of Science and Technology (MOST) and National Natural Science Foundation of China (NSFC) for fiscal support (Project Research Grants No. 201204510, No. 2005DIB5JI42 and No. 31170588) and the Inventory Institute of Beijing Forestry for its data. The authors also thank the two anonymous reviewers for improving the scientific quality of this manuscript and the editors for their careful work. References Adams, D.M., Ek, A.R., 1974. Optimizing the management of uneven-aged forest stands. Can. J. For. Res. 4, 274–287. Affleck, D.L.R., 2006. Poisson mixture models for regression analysis of stand-level mortality. Can. J. For. Res. 36, 2994–3006. Álvarez-González, J.G., Dorado, F.G., Gonzalez, A.D.R., López-Sánchez, C.A., Gadow, K.V., 2004. A two-step mortality model for even-aged stands of Pinus radiate D.Don in Galicia (Northwestern Spain). Ann. For. Sci. 61, 439–448. Andreassen, K., 1994. Development and yield in selection forest. Medd. Skogforsk 47, 1–37. Barry, S.C., Welsh, A.H., 2002. Generalized additive modeling and zero inflated count data. Ecol. Mod. 157, 179–188. Boswell, M.T., Patil, G.P., 1970. Chance mechanisms generating the negative binomial distributions. In: Patil, G.P. (Ed.), Random Counts in Scientific Work. The Pennsylvania State University Press, University Park, Penn, pp. 3–22. Cameron, A.C., Trivedi, P.K., 1998. Regression Analysis of Count Data. Cambridge University Press, Cambridge. Cameron, A.C., Trivedi, P.K., 2005. Microeconometrics: Methods and Applications. Cambridge University Press, Cambridge. Cao, Q.V., 2006. Predictions of individual-tree and whole-stand attributes for loblolly pine plantations. For. Ecol. Manage. 236, 342–347. Crepon, B., Duguet, E., 1997. Research and development, competition and innovation-pseudo-maximum likelihood and simulated maximum likelihood methods applied to count data models with heterogeneity. J. Econometrics 79, 355–378. Cunniningham, R.B., Lindenmayer, D.B., 2005. Modeling count data of rare species: some statistical issues. Ecology 86, 1135–1142. Eberhardt, L.L., 1967. Some developments in distance sampling. Biometrics 23, 207– 216. Eerujäinen, K., Miina, J., Valkonen, S., 2007. Models for the regeneration establishment and the development of established seedlings in uneven-aged, Norway spruce dominated forest stands of southern Finland. For. Ecol. Manage. 242, 444–461. Eskelson, B.N.I., Temesgen, H., Barrett, T.M., 2009. Estimating cavity tree and snag abundance using negative binomial regression models and nearest neighbor imputation methods. Can. J. For. Res. 39, 1749–1765. Ferguson, D.E., Stage, A.R. and Boyd, R.J., 1986. Predicting regeneration in the grand fir-cedar-hemlock ecosystem of the northern rocky mountains. For. Sci. Monog. 26, 41 p. Flores, O., Rossi, V., Mortier, F., 2009. Autocorrelation offsets zero-inflation in models of tropical saplings density. Ecol. Model. 220, 1797–1809. Fridman, J., Ståhl, G., 2001. Models for prediction of basal area mean diameter and number of trees forest stands in southeastern Norway. Scand. J. For. Res. 16, 455–466. Fortin, M., DeBlois, J., 2007. Modeling tree recruitment with zero-inflated models: the example of hardwood stands in Southern Quebec. Canada. For. Sci. 53, 529– 539. Gurmu, S., 1997. Semi-parametric estimation of hurdle regression models with an application to Medicaid utilization. J. Appl. Econom. 12, 225–242. Gurmu, S., Trivedi, P., 1996. Excess zeros in count models for recreationaltrips. J. Bus. Econ. Stat. 14, 469–477.
215
Hall, D.B., 2000. Zero-inflated Poisson and binomial regression with random effects: a case study. Biometrics 56, 1030–1039. Heilbron, D., 1994. Zero-altered and other regression models for count data with added zeros. Biometrical J. 36, 531–547. Jackman S. 2011. pscl: Classes and Methods for R Developed in the Political Science Computational Laboratory, Stanford University. Department of Political Science, Stanford University, Stanford, California. R package version 1.04.1. Available from:
Contents lists available at SciVerse ScienceDirect
Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco
Predicting tree recruitment with negative binomial mixture models Xiongqing Zhang a, Yuancai Lei a,⇑, Daoxiong Cai b, Fengqiang Liu c a
Research Institute of Forest Resource Information Techniques, Chinese Academy of Forestry, Beijing 100091, PR China Experimental Center of Tropical Forestry, Chinese Academy of Forestry, Pingxiang 532600, Guangxi, PR China c Southern Swedish Forest Research Centre, Swedish University of Agricultral Sciences, 232 38 Alnarp, Sweden b
a r t i c l e
i n f o
Article history: Received 4 October 2011 Received in revised form 31 December 2011 Accepted 19 January 2012 Available online 18 February 2012 Keywords: Tree recruitment Negative binomial model Zero-inflated negative binomial model Hurdle negative binomial model Chinese pine
a b s t r a c t Tree recruitment models play an important role in simulating stand dynamic processes. Periodic tree recruitment data from permanent plots tend to be overdispersed, and frequently contain an excess of zero counts. Such data have commonly been analyzed using count data models, such as Poisson model, negative binomial models (NB), zero-inflated models, and Hurdle models. Negative binomial mixture models (zero-inflated negative binomial model, ZINB; Hurdle negative binomial model, HNB) including NB model were used in this study to predict tree recruitments of Chinese pine (Pinus tabulaeformis) in Beijing. ZINB model and HNB model were suitable for dealing with excess zero counts, for which two equations are created: one predicting whether the count occurs (logistic function) and the other predicting differences on the occurrence of the count (NB model). Based on the model comparisons, the results showed that negative binomial mixture models performed well in modeling tree recruitment, and ZINB model was the best model of negative binomial mixture models. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction Forest dynamic process includes three components: tree mortality, tree growth, and recruitment of new trees. For these decades, foresters have made many researches on the first two components (e.g., Monserud and Sterba, 1999; Zhao et al., 2004; Cao, 2006; Zhang and Lei, 2009; Subedi and Sharma, 2011, etc.). Dynamic models may neglect the simulation of the third component, tree recruitment, because of the nature of data or difficulty in modeling. Recruitment, sometimes referred to as ingrowth (Shifley et al., 1993) is usually quantified by means of number of trees or saplings that reach or exceed a specific threshold size (e.g., 1.3 m height, 5 cm DBH over bark, etc.) over a certain period as a result of different regeneration processes like establishment, growth, and mortality of saplings (Lexerød, 2005; Lexerød and Eid, 2005). Recruitment contributes significantly to volume production in multi-layered forests (Andreassen, 1994). Recruitment model is part of the forest management tools that enable prediction of the development of forest stands, and disregarding recruitment will give biased prediction of future forest growth and yield. Because of the importance of tree recruitment, many researchers have developed models to predict recruitment over a certain period, including linear, nonlinear, and Tobit models (e.g., Adams and Ek, 1974; Shifley et al., 1993; Liang et al., 2005).
⇑ Corresponding author. Tel.: +86 10 62889199; fax: +86 10 62888315. E-mail addresses: [email protected] (X. Zhang), [email protected] (Y. Lei). 0378-1127/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.foreco.2012.01.028
Considering the fact that in permanent sample plots a relatively high number of the plots have no occurrences of recruitment even over periods of several years, it means that recruitment data are bounded and characteristically exhibit varying degrees of dispersion and skewness in relation to the mean. Additionally, the data often contain an excess number of zero counts. Yet least squares method implicitly presumes that the data are Gaussian distributed with constant variance, or at least satisfy the Gauss–Markov conditions. If the method is still used to deal with the data with large proportion of zero counts, the estimated results will be biased. Alternatively, if only recruitment observations are used for model development, recruitment will be overestimated. To overcome this problem, Ferguson et al., (1986) used a twostep approach to predict recruitment in the PROGNOSIS model (Stage, 1973; Wykoff, 1986). With a two-step approach, the probability of zero recruitment in the subsequent measurement interval, given current stand conditions, is obtained from a logistic regression in the first step. Then the amount of recruitment, given that some is known to occur, is described by a recruitment function fit by least squares in the second step. The method is notable in its recognition of excess number of zero fraction as an aspect of the response distribution not well solved by standard nonlinear regression models and has been utilized with some modifications in other researches (e.g., Woollons, 1998; Fridman and Ståhl, 2001; Álvarez-González et al., 2004; Lexerød, 2005; Fortin and DeBlois, 2007). But the method implies that a large zero fraction reflects additional structure in the mean function, structure that should carry some ecological relevance if it is to be recognized.
210
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
Additionally, this method does not admit a well-defined stochastic structure that can serve as a basis for inference (Affleck, 2006). In recent years, there has been considerable interest in models for count data, such as manufacturing defects (Lambert, 1992), patent applications (Crepon and Duguet, 1997), species abundance (Welsh et al., 1996; Hall, 2000; Cunniningham and Lindenmayer, 2005), medical consultations (Gurmu, 1997), use of recreational facilities (Gurmu and Trivedi, 1996; Shonkwiler and Shaw, 1996). There are also a few studies so far have addressed this issue in forestry (e.g., Keefe, 2004; Rathbun and Fei, 2006; Eskelson et al., 2009; Flores et al., 2009; Li et al., 2011). Fortin and DeBlois (2007) applied zero-inflated Poisson model (ZIP) and zero-inflated discrete Weibull (ZIdiW) to model tree recruitment of hardwood stands. But the base condition for the positive counts in the ZIP model is that the variance equals the expectation, and to ZidiW model, it is difficult for estimating an additional vector of parameters. Negative binomial model (NB) is another basic count-data model. But, it is not suitable to deal with the data with excess of zero counts (Zeng et al., 2009). Zero-inflated negative binomial model (ZINB) and Hurdle-negative binomial model (Hurdle-NB) were common methods to analyze this type of data (Affleck, 2006; Zeileis et al., 2008). The objective of this study was to develop and compare the ZIP model, Hurdle-Poisson model, NB model, ZINB model, and HurdleNB model for predicting tree recruitment over a period of time. The modelings were restricted to include available independent-variables directly from forest inventories for forest management planning.
count data, and the probability mass function (PMF) for Poisson is characterized by as
PðY ¼ yÞ ¼
ek ky y!
ð1Þ
where y refers to random variable of recruitment count, y = 0, 1, 2. . ., and k > 0. A Poisson regression model is obtained by relating the mean k to a vector of independent variables X, by k ¼ ExpðXbÞ; where b is a vector of regression coefficients to be estimated. A character of the Poisson probability function (Eq. 1) is that the mean and the variance are equal, that is Var½YjX ¼ E½YjX ¼ k. When data are not suitable for the Poisson distribution, it is typically because of overdispersion, meaning the model’s variance exceeds the mean value. The Negative binomial regression (NB) can be considered a generalization of Poisson regression and assumes that the conditional mean k of Y is not only determined by X but also a heterogeneity component unrelated to X, which addresses the issue of overdispersion by including a dispersion parameter to accommodate the unobserved heterogeneity in the count data (MacNeil et al., 2009). The PMF for the negative binomial is given by
!h1
Cðy þ h1 Þ h1 PðY ¼ yÞ ¼ 1 1 Cðy þ 1ÞCðh Þ h þ k
k
h1 þ k
y ð2Þ
where h represents the dispersion parameter. The mean is k, as the Poisson model, but the variance is k þ hk2 , thus allowing the variance to exceed k. NB regression model is obtained by relating the mean k to a vector of independent variables, k ¼ ExpðXbÞ. The pre^ corresponds to the expectation k of the distribution. diction y
2. Data 3.2. Zero-inflated models The data, collected by the Inventory Institute of Beijing Forestry, consisted of a systematic sample of permanent plots with a 5-year re-measurement interval. The plots (squared plots), 0.067 ha each, were in Chinese pine (Pinus tabulaeformis Carr.) plantations situated on upland sites of Beijing (Fig. 1). The data consisted of measurements from 132 plots obtained between 1986 and 2001. Summary statistics of stand variables are presented in Table 1.
3. Method Tree recruitment (a threshold size of diameter at breast height for tree recruitment in China is DBH = 5 cm) was assumed to relate with location variables, site conditions and stand characteristics. The models should be biologically realistic to ensure logical behavior when applied outside the range of the modeling data. Thus, an appropriate biological interpretation was the main criterion with respect to variable selection (Lexerød, 2005). Independent variables, with an appropriate biological interpretation, were selected among those describing location (altitude Alt, slope Slo), site conditions (dominant height H), and stand characteristics (stand density N, stand arithmetic mean diameter Dm, stand basal area B, relative spacing, etc.). In this study, these independent variables above would be used to predict tree recruitment. Moreover, we analyzed the multicollinearity between independent variables through variance inflation factor (VIF) test. A common rule of thumb is that if VIF > 5 then multicollinearity exists. The VIF test was performed through DAAG package in R software.
Although NB model is better than Poisson model for dealing with dispersed data (Liu and Cela, 2008), there are some situations where a major source of overdispersion is a preponderance of zero counts. This suggests consideration of more flexible count data models such as zero-inflated models. They are mixture models combining a count distribution with a point mass at zero. In zero-inflated models there are two sources of zeros: zeros may come from both the point mass and from the count component (Lambert, 1992; Eskelson et al., 2009). Actually, two regression equations are created in zero-inflated models: one predicting whether the count occurs (logistic function) and the other predicting differences on the occurrence of the count (Poisson or NB) (Karazsia and van Dulmen, 2008). Supposing a discrete variable Y (recruitment counts) obeys to zero-inflated distribution, the PMF for zero-inflated distribution is given by:
PðY ¼ yÞ ¼
p þ ð1 pÞf ð0Þ
y¼0
ð1 pÞf ðyÞ
y>0
ð3Þ
Where p is the probability of belonging to the point mass at zero and (1 p) is the probability of belonging to the count component. Logistic regression was commonly used to fit point mass compop nent: logitðpÞ ¼ Logð1p Þ ¼ Xd, where d is a vector of parameters to be estimated. In the count component, it was the same to Poisson or NB model k ¼ ExpðXbÞ.
3.1. Negative binomial model
3.2.1. Zero-inflated Poisson model (ZIP) In count component, if Y obeys to the Poisson distribution with parameter k in Eq. (3), then we get the ZIP distribution. The PMF for ZIP is given by:
Count data models are a subset of discrete-response regression models and aim to explain the number of occurrences or counts of an event. Poisson regression is the simplest regression model for
8 < p þ ð1 pÞExpðkÞ PðY ¼ yÞ ¼ ð1 pÞExpðkÞky : y!
y¼0 y>0
ð4Þ
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
211
Fig. 1. The spatial location of Chinese pine in Beijing. Lines in the figure are the boundaries of administrative region in Beijing.
Table 1 Summary statistics of stand-level variables, including tree recruitment counts in a stand Tr, stand age A, stand dominant height H, stand density N, stand basal area B, stand arithmetic mean diameter Dm, altitude Al, slope Slo and the relative spacing Rs. The relative spacing is a function of dominant height and stand density: Rs = (10,000/ N)0.5/H. Variables
Min
Max
Mean
SD
Tr (trees/plot) A (years) H (m) N (trees/ha) B (m2/ha) Dm (cm) Rs (trees/ham) Alt (m) Slo (°)
0 12 2.5 132 0.3580 5.8 0.1712 140 5
59 55 17.4 2164.1790 28.3105 16.1260 3.4512 1400 40
7.9242 28.3561 6.6402 899.2537 8.0792 9.9097 0.6823 443.2500 25.4318
12.0061 8.5132 2.9006 453.5848 5.6711 2.3951 0.4675 266.9939 8.1082
SD, standard deviation.
3.2.2. Zero-inflated negative binomial model (ZINB) If the NB distribution is used as count component, the mixture distribution is zero-inflated negative binomial, and the PMF for ZINB is given by:
PðY ¼ yÞ ¼
8 > < p þ ð1 pÞ > : ð1 pÞ
h1 h1 þk 1
1=h
Cðyþh Þ Cðyþ1ÞCðh1 Þ
h
y¼0 1
h
1
1=h
þk
y
k h þk 1
ð5Þ y>0
When h ! 0, the ZINB is changed into ZIP model. 3.3. Hurdle models Hurdle models, originally proposed by Mullahy (1986) in the econometrics literature (e.g., Cameron and Trivedi, 1998, 2005) are another model class capable of dealing with excess zero counts. They are also called two-part models (Heilbron, 1994), thus similar to zero-inflated models. Hurdle models consider count outcome generated by two systematically different statistical processes, a binomial distribution determining if a count outcome is zero or nonzero and a truncated-at-zero distribution for count data governing all positive counts conditional on count component. However, they are slightly different from zero-inflated model with all zero counts from two different sources and assume that zero counts might come from a single statistical process (Liu and Cela, 2008). Similar to the zero-inflated models, logistic function often was used to model point mass at zero. About count component, Poisson model or negative models were used in this case. 3.3.1. Hurdle-Poisson model (HP) One of the most popular Hurdle model is called HP model, which is the combination of a logit regression modeling point mass at zero and a truncated Poisson regression modeling positive counts conditional on count component. Its PMF is given as:
212
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
8
y¼0 ð6Þ
y>0
3.3.2. Hurdle-negative binomial model (HNB) Similar to ZINB model, HNB model is another important Hurdle model, which is the combination of a logit regression modeling point mass at zero and a truncated negative binomial regression modeling count component. Its PMF is given as:
( PðY ¼ yÞ ¼
p
h i1 y 1 ðyþh1 Þ k ð1 pÞ CðhC1 ð1 þ hkÞh 1 kþh1 ÞCðyþ1Þ
y¼0 y>0 ð7Þ
^ yi y ffiffiffiffiffi i pr ¼ p v^ i
ð11Þ
^i is the estimated mean and v^ i is the variance of y ^i of the ith where y observation. To evaluate the adequacy of the models for predicting the overall counts of tree recruitments, an approximate v2 goodness-of-fit was calculated using observed and predicted counts (Eskelson et al., 2009):
v2 ¼
P 2 m X ½#ðyi ¼ jÞ Pðyi ¼ jÞ P Pðyi ¼ jÞ j¼1
ð12Þ
where # represents the frequency of observations yi in count class j, m represents the number of count classes, and P(y = j) is the predicted probability that an observation belongs to count class j.
3.4. Parameter estimation
5. Results and discussion
All vectors of parameters b, d, and h can be estimated through the optimization of the likelihood function, i.e., the maximum likelihood method (ML). Standard errors of model parameters were estimated with two methods. The first one, based on standard ML theory from the inverse of the k k (k is the number of parameters estimated in the model) observed Fisher information matrix (Affleck, 2006):
According to the VIF test, we found that there was no multicollinearity among stand age, stand dominated height, relative space, altitude, slope, and any two of the combinations of stand density, stand arithmetic mean diameter, and stand basal area (VIF < 5). Table 2 lists the parameter estimates and their standard deviation errors for NB model, ZIP model, ZINB model, HP model, and HNB model. In the positive count component of the models, the majority of variables were significant at the 0.05 level; and for ZINB model, parameters of dominant height and slope were significant at the 0.1 level. Dominant height reflects the site situation. Vanclay (1989) found that the number of recruits became larger on good sites. In this study, the relation between dominant height and tree recruitment was positive (i.e. the greater H, the greater the number of tree recruitment). Stand basal area is a measure of stand competition, reflecting both size, and density of trees in a stand, and as it increased, tree recruitments decreased. This was consistent with the results of previous studies (Vanclay, 1989, 1992; Lexerød, 2005). Similarly, the effect of stand arithmetic mean diameter (size index) was negative (i.e. the greater Dm, the smaller number of tree recruitment). This meant that tree recruitment was more likely in forests with many small trees compared to forests with larger trees (Lexerød, 2005). In contrast, as the altitude became larger, tree recruitments increased. And the parameter estimates of altitude for these five models were nearly same. On the other hand, as the slope increased, tree recruitments decreased. In the zero component of the models, the parameters of ZINB model were all significant at the 0.1 level (Table 2).Two model fit statistics are gi^ ; yÞ for ZINB model was 737.5109, as ven in Table 2. The 2Lðu compared to 1241.5210, 1240.2000, 744.4256, 748.6767 for ZIP model, HP model, NB model, and HNB model, respectively; AIC was 761.5109, versus 1257.5210, 1264.1680, 762.4734, 768.6767. ^ ; yÞ and AIC of the ZINB These numbers showed that the 2Lðu model were much smaller than ZIP model and HP model, and slightly smaller than NB model and HNB model. It meant that the ZINB regression provided the ‘‘best’’ model. Pearson residuals are shown in Fig. 2. All five regression models overestimated for small recruitment counts and underestimated for large counts. The residuals ranged from 2 to 6 for ZIP model and HP model, while ranged from 1 to 3 for NB model, ZINB model, and HNB model. The v2 value of ZINB model was 2.9944, as compared to 3.8499, 3.9821, 44.2936, 46.1286 for NB model, HNB model, ZIP model, and HP model, respectively, and p-value >0.05 for negative binomial mixture models, p-value <0.0001 for Poisson mixture models (For the v2 -test, classes were grouped together to obtain expected frequencies of 5 or greater). The lower the v2 statistic, the better the model fit. ZIP model and HP model accommodated the large zero fraction, but their Poisson component fit the positive count poorly. Fig. 3 showed the empirical
" ^ m ¼ va^r ½u
@ 2 Lð/; yÞ ju^ @/@/0
#1
"
n X @ 2 Lð/; yi Þ ¼ ju^ @/@/0 i¼1
#1 ð8Þ
^ yÞ is the log-likelihood of the developed model, u ^ is the where Lð/; ^ includes the vector of estik 1 vector of estimated parameters (u mated parameter b, d, and h depending on the model fitted), n is number of observations, and y is a column vector whose elements are the responses yi. The variance estimator is consistent if the model is correctly specified. But it is biased if the data are overdispersed. Thus, standard errors were also calculated from the matrix product.
^ r ¼ va^r½u ^ m va^r ½u
" # n X @Lð/; yi Þ @Lð/; yi Þ ^ m r ½u j ^ va^ u @/ @/0 i¼1
ð9Þ
The standard errors obtained with Eq. (9) are consistent in model specification and require only the unbiased ML estimating equations (Cameron and Trivedi, 1998). Estimation of parameters was implemented with pscl package in R software (Jackman, 2011). 4. Model selection and goodness-of-fit Negative binomial model, zero-inflated models and Hurdle models calibrated with the same data set can be compared through ^ ; yÞ) as well as Akaike information crilog-likelihood values (2Lðu ^ ; yÞ indicate that modterion (AIC). Smaller values of AIC and 2Lðu el is better.
^ ; yÞ þ 2k AIC ¼ 2Lðu
ð10Þ
Vuong test is considered as a better method to compare two non-nested models for count data (Vuong, 1989). The Vuong tests were performed in R software. Goodness-of-fit was also evaluated from predictive and residual diagnostics. Residuals measure the departure of predicted values from observed values of the dependent variable. The most commonly used residual is the raw resid^. For the classical linear model with normally ual: r ¼ y y ^Þ N½0; r2 , so that in large distributed homoscedastic error ðy y samples the raw residual has the desirable properties of being symmetrically distributed around zero with constant variance. However, for count data, the raw residual is heteroscedastic and asymmetric. The obvious correction for heteroscedasticity is the Pearson residual (Cameron and Trivedi, 1998):
213
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215 Table 2 Parameter estimations and fit statistics of ZIP model, HP model, NB model, ZINB model, and HNB model. Parameter
ZIP model
HP model
Estimation
Sig.
Positive count component of the model **** Intercept 5.3011 (0.1861) **** H 0.0721 (0.0184) **** B 0.0384 (0.0102) **** Dm 0.3592 (0.0282) **** Alt 0.0014 (0.0001) **** Slo 0.0233 (0.0038) Zero component of the model Intercept – – **
Estimation
Sig.
Estimation
Sig.
Estimation
Sig.
5.2921 (0.1881) 0.0681 (0.0186) 0.0372 (0.0103) 0.3562 (0.0286) 0.0014 (0.0001) 0.0231 (0.0038)
****
4.8353 (0.6504) 0.1220 (0.0545) 0.0559 (0.0210) 0.3515 (0.0808) 0.0014 (0.0004) 0.0235 (0.0137)
****
4.8954 (0.6197) 0.1321 (0.0507) 0.0567 (0.0317) 0.3631 (0.0792) 0.0014 (0.0004) 0.0228 (0.0136)
****
5.3638 (0.6455) 0.0875 (0.0525) –
****
–
–
–
–
*
–
–
**
–
–
*
–
–
*
–
–
0.7300 (0.1090) 744.4256 762.4734
****
247.7297 (140.0640) 8.6282 (5.1010) 3.9387 (2.0332) 284.6551 (164.2425) 0.0350 (0.0205) 0.7944 (0.1515) 737.5109 761.5109
N h
–
^ ; yÞ 2Lðu AIC
1241.5210 1257.5210
1240.2000 1264.1680
Rs
– *
HNB model
Sig.
1.1600 (0.4901) –
B
0.1042 (0.0502) –
ZINB model
Estimation
0.1487 (0.0761) 0.1327 (0.0598) 0.7956 (0.4708) 0.0009 (0.0005) –
H
NB model
****
****
****
****
****
**
**
****
***
**
***
*
****
****
*
–
–
*
0.1487 (0.0761) 0.1327 (0.0598) 0.7956 (0.4708) 0.0009 (0.0005) 0.7641 (0.2877) 748.6767 768.6767
*
***
*
*
*
*
****
20
30
40
50
60
0
30
40
HP
20
30
40
50
60
Observed recruitments
50
60
50
60
-2 0 2
4 6
ZINB 1 2 3 4
10
20
Observed recruitments
-1 0
10
Observed recruitments
Pearson residual
10
-2 0 2 4 6
Pearson residual
1 2 3 4
ZIP
-1
Pearson residual Pearson residual
0
0
10
20
30
40
Observed recruitments
1 2 3 -1
Pearson residual
HNB
0
10
20
30
40
50
–
0.4210 (0.0745) 0.0013 (0.0004) 0.0224 (0.0134)
Values in parentheses are standard deviation errors, and Sig. is significance. * P < 0.1. ** P < 0.05. *** P < 0.01. **** P < 0.001.
NB
*
60
Observed recruitments Fig. 2. Pearson residual plots of tree recruitments for ZIP model, HP model, NB model, ZINB model and HNB model.
***
*
**
*
*
***
214
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
Fig. 3. Comparison between the empirical distribution of observed counts (black solid line) and the empirical distribution of predicted counts (according to the ZIP model, HP model, NB model, ZINB model, and HNB model).
distribution of the observations and the empirical distribution of the predictions of the five models. It was clear that NB model fit the observed count outcomes as well as ZINB model and HNB model, which was much better than ZIP model and HP model. A ZIP or HP will reflect data accurately when overdispersion is caused by a preponderance of zeros. If overdispersion is attributed to factors beyond the inflation of zeros, a ZINB or HNB model is more appropriate (Long and Freese, 2006). And the ZIP model or HP model accommodated the excess of zero component, but its Poisson component fit the positive component poorly (Fig. 3). Poisson is highly restrictive as the variance of the outcome is assumed to equal its mean. Count data often exhibits overdispersion. NB offers a dispersion parameter which explains the overdispersion of positive component (MacNeil et al., 2009). Therefore, NB mixture models (NB, ZINB, and HNB) performed better than Poisson mixture models (ZIP, and HP). However, we should know that the dispersion parameter h of the NB, ZINB, or HNB model is scale dependent when dealing with real data (Pielou, 1957; Eberhardt, 1967). For the Hurdle models, the zero hurdle component describes the probability of a zero count while for the zero-inflated models, the zero-inflation component predicts the probability of observing a zero count from the point mass component and count component (Zeileis et al., 2008). Overall, both models led to the same qualitative results and very similar model fits (Table 2; Figs. 2 and 3). NB model provides flexible structures for describing empirical count data. The model showed closer agreement with the tree recruitment data than the ZIP model and HP model on the basis of AIC and in terms of their ability to reproduce the observed count distribution (Table 2; Fig. 3). The improvement fit of the NB model could be ascribed to unobserved heterogeneity in tree recruitment across and with plots over time. The NB model could also be attributed to contagion, wherein the occurrence of one event augments the probability of successive events (Boswell and Patil, 1970). However it has only one mode, like the Poisson, does not deal with the extra-zero problem even though it allows for extra dispersion (Cunniningham and Lindenmayer, 2005). And ZINB model or HNB has three advantages. Firstly, by allowing separate models for the two components, tree recruitment can be modeled com-
pared with simply assuming the data is approximately Poisson or NB. Secondly, the zero component of the model and positive count component of the model of the two models are directly interpretable in terms of the empirical features of the data. Thirdly, the two models will better model the variability seen in this data. This will ensure that the information in the data is used appropriately (Barry and Welsh, 2002). Additionally, in this study, Vuong test-statistic V was 4.9617, pvalue <0.0001 for comparing HNB model with HP model, and V was 5.0648, p-value <0.0001 for comparing ZINB model with ZIP model. This showed that ZINB model or HNB model was much better than ZIP model and HP model. In the same way, for comparing the ZINB model with NB model, V was 1.6963, p-value was 0.0449, which indicated that ZINB model performed better than NB model at the significant level 0.05. We also compared ZINB model with HNB model, V was 1.4883, p-value was 0.0683. It showed that ZINB model was better than HNB model at significant level 0.1. Overall, ZINB model provided more opportunities for exploring when the probability of tree recruitment and the number of tree recruitment may be driven by different factors. Although the ZINB model is proved to have a good fit for recruitment data with an excess of zero counts, the extra zeros are not the only feature of the recruitment patterns. This typically occurs when considering a rare species, a mast species, counting recruitment of each different tree species in mixed forests, recruitment on small plots, or recruitment in short census interval for slow-growing tree species. With increasing census interval, the probability of zero occurrence gradually decreases (Yao et al., 2001).
6. Conclusion In this paper we have reviewed several count data models for analyzing tree recruitment. When the data is overdispersed, the negative binomial model is a good choice. The NB model is more flexible as it allows for the variance to be greater than mean and considers the observed heterogeneity in the model (Yaacob et al., 2010). However, the phenomenon of tree recruitment was always characterized by large number of zero counts. In this case, zero
X. Zhang et al. / Forest Ecology and Management 270 (2012) 209–215
negative binomial model and Hurdle negative binomial model might give a more satisfactory fit to the data. Based on the model comparisons, ZINB model performed better than other four models for modeling tree recruitments. These techniques can also be extended to other applications, such as stand mortality, snag abundance, species abundance, tree generation, etc. (e.g., Barry and Welsh, 2002; Affleck, 2006; Eerujäinen et al., 2007; Eskelson et al., 2009). Acknowledgments The authors express their appreciation to State Forestry Administration (SFA), the Ministry of Science and Technology (MOST) and National Natural Science Foundation of China (NSFC) for fiscal support (Project Research Grants No. 201204510, No. 2005DIB5JI42 and No. 31170588) and the Inventory Institute of Beijing Forestry for its data. The authors also thank the two anonymous reviewers for improving the scientific quality of this manuscript and the editors for their careful work. References Adams, D.M., Ek, A.R., 1974. Optimizing the management of uneven-aged forest stands. Can. J. For. Res. 4, 274–287. Affleck, D.L.R., 2006. Poisson mixture models for regression analysis of stand-level mortality. Can. J. For. Res. 36, 2994–3006. Álvarez-González, J.G., Dorado, F.G., Gonzalez, A.D.R., López-Sánchez, C.A., Gadow, K.V., 2004. A two-step mortality model for even-aged stands of Pinus radiate D.Don in Galicia (Northwestern Spain). Ann. For. Sci. 61, 439–448. Andreassen, K., 1994. Development and yield in selection forest. Medd. Skogforsk 47, 1–37. Barry, S.C., Welsh, A.H., 2002. Generalized additive modeling and zero inflated count data. Ecol. Mod. 157, 179–188. Boswell, M.T., Patil, G.P., 1970. Chance mechanisms generating the negative binomial distributions. In: Patil, G.P. (Ed.), Random Counts in Scientific Work. The Pennsylvania State University Press, University Park, Penn, pp. 3–22. Cameron, A.C., Trivedi, P.K., 1998. Regression Analysis of Count Data. Cambridge University Press, Cambridge. Cameron, A.C., Trivedi, P.K., 2005. Microeconometrics: Methods and Applications. Cambridge University Press, Cambridge. Cao, Q.V., 2006. Predictions of individual-tree and whole-stand attributes for loblolly pine plantations. For. Ecol. Manage. 236, 342–347. Crepon, B., Duguet, E., 1997. Research and development, competition and innovation-pseudo-maximum likelihood and simulated maximum likelihood methods applied to count data models with heterogeneity. J. Econometrics 79, 355–378. Cunniningham, R.B., Lindenmayer, D.B., 2005. Modeling count data of rare species: some statistical issues. Ecology 86, 1135–1142. Eberhardt, L.L., 1967. Some developments in distance sampling. Biometrics 23, 207– 216. Eerujäinen, K., Miina, J., Valkonen, S., 2007. Models for the regeneration establishment and the development of established seedlings in uneven-aged, Norway spruce dominated forest stands of southern Finland. For. Ecol. Manage. 242, 444–461. Eskelson, B.N.I., Temesgen, H., Barrett, T.M., 2009. Estimating cavity tree and snag abundance using negative binomial regression models and nearest neighbor imputation methods. Can. J. For. Res. 39, 1749–1765. Ferguson, D.E., Stage, A.R. and Boyd, R.J., 1986. Predicting regeneration in the grand fir-cedar-hemlock ecosystem of the northern rocky mountains. For. Sci. Monog. 26, 41 p. Flores, O., Rossi, V., Mortier, F., 2009. Autocorrelation offsets zero-inflation in models of tropical saplings density. Ecol. Model. 220, 1797–1809. Fridman, J., Ståhl, G., 2001. Models for prediction of basal area mean diameter and number of trees forest stands in southeastern Norway. Scand. J. For. Res. 16, 455–466. Fortin, M., DeBlois, J., 2007. Modeling tree recruitment with zero-inflated models: the example of hardwood stands in Southern Quebec. Canada. For. Sci. 53, 529– 539. Gurmu, S., 1997. Semi-parametric estimation of hurdle regression models with an application to Medicaid utilization. J. Appl. Econom. 12, 225–242. Gurmu, S., Trivedi, P., 1996. Excess zeros in count models for recreationaltrips. J. Bus. Econ. Stat. 14, 469–477.
215
Hall, D.B., 2000. Zero-inflated Poisson and binomial regression with random effects: a case study. Biometrics 56, 1030–1039. Heilbron, D., 1994. Zero-altered and other regression models for count data with added zeros. Biometrical J. 36, 531–547. Jackman S. 2011. pscl: Classes and Methods for R Developed in the Political Science Computational Laboratory, Stanford University. Department of Political Science, Stanford University, Stanford, California. R package version 1.04.1. Available from: