Abstract: This study presents a new characterization on tree diameter distribution problem. The main purpose of this study is to investigate the relationship between the stochastic Gompertz shape diameter growth model and diameter distribution law using stochastic differential equation methodology. The probabilistic characteristics of diameter growth model, such as the univariate transition probability density of tree diameter, the mean and variance of tree diameter is established. A generalized form of mean volume, volume increment and diameter increment of a tree is introduced which implicitly incorporates an age-height-dependent transition probability density function of tree diameter. To model the tree diameter distribution, as an illustrative experience, is used a real data set from repeated measurements on permanent sample plots of pine stands in Dubrava district at Lithuania. The results are implemented in the symbolic computational language MAPLE 11.
INTRODUCTION
Nowadays, in commercial boreal forest, forest management needs to utilize fundamental principles of forest dynamics. Forest growth, like many other natural processes, is subjected to various disturbances. Thus, the predictions computed from deterministic models limit our understanding of other possible outcomes of stand diameter and height. One way of quantifying stand growth under random perturbations is with diffusion process models.
The processes of growth play an important role in various applied areas, such as biology, medicine, biochemical industry. The environment of any real system is in general not constant but shows random fluctuations. Despite this, the growth model historically has crystallized as the deterministic logistic type process (Kar and Matsuda, 2007; Sakanoue, 2007). The forest growth is usually modeled with a logistic model (Garcia, 2005). The parameters of logistic diameter models are not directly measurable but they are estimated from the observed data set.
Stand as a community of trees is the main component of the forest. Stand consists of trees with different diameters and heights. Those differences depend on a lot of unsearchable genetic and environmental factors, therefore it leads to consideration that diameter of a tree is a random variable which depends on the age and height. Diameter dynamic is affected by many processes and varies among stands (Temesgen and Gadow, 2004). Stochastic diameter growth models allow us to reduce the unexplained variability of a diameter and to implement the randomness phenomenon, which makes a stochastic influence on diameter growth process useful in practical applications. Over the years an extensive amount of research has been devoted to the randomness of stand growth since the pioneer work of Suzuki (1971) and the successive works of Tanaka (1986) and Rupsys (2007). There are two types of approaches for this purpose. The first approach is based on ‘environment’ stochasticity, introducing a diffusion term in the ordinary differential equation of diameter dynamic (Suzuki, 1971; Tanaka, 1986; Rupsys, 2007; Rupsys et al., 2007; Rupsys and Petrauskas, 2009, 2010). The second approach is based on demographic stochasticity in which the tree size X is a random variable (Boungiorno, 2001; Lohmander and Mohammadi, 2007). In this study we follow the first approach.
The main purpose of this study is to develop the age-height-dependent probability density function on diameter size using measurements of tree variables such as age, height and diameter. This study not only provides useful stochastic models for the diameter growth modeling, but shows that it is possible to relate the diameter growth model and the diameter distribution model. The distributions of tree diameter size in stands describe forest structure and can be used for the assessment of stand volume and biomass, forest biodiversity and density management. Knowledge of the predicted age-height-dependent distribution function enables a more differentiated prediction of the assortment for a stand. This is not possible with commonly used distribution functions or yield-tables.
In even-aged stands various distribution functions, such as negative exponential, Pearson, gamma, lognormal, beta, Weibull, Johnson, Gram-Charlier, have been used in describing the diameter distributions (Mehtatalo, 2005). In uneven-aged stands have been used bivariate distributions and density mixtures (Wang and Rennolls, 2007; Wang et al., 2008). In this work we motivate the use of stochastic differential equations in forestry. The methodology is to consider a univariate distribution as arising from univariate diameter growth stochastic dynamical system. The system fluctuations, generally infiltrated from outside, are defined by a one-dimensional standard Wiener process.
In this study, a univariate age-dependent stochastic differential equation methodology of tree diameter distribution is expanded into age-height-dependent distribution function methodology. The Gompertz homogeneous and nonhomogeneous growth models are applied to analyze the trend of tree diameter, taking the height as an exogeneous variable that affects the diameter of a tree. The choice of the height among other possible exogenous factors is justified by the significant correlation with the height and age (Skovsgaard and Vanclay, 2008; Garcia, 2009). This approach is rather different from the univariate diffusion models (Suzuki, 1971; Tanaka, 1986; Rupsys and Petrauskas, 2009, 2010), since these distributions are not related with the dynamic of height.
MATERIALS AND METHODS
Growth model: Let study the dynamic behavior of tree diameter (diameter at breast height) and its relationship with diameter distribution law. For determination of diameter growth we suppose that dynamic of tree diameter is expressed in terms of the Gompertz shape stochastic differential equation with multiplicative noise. The Gompertz deterministic model is a classical continuous model useful in describing population dynamic. It was introduced by Gompertz (1825) to analyze population dynamic and to determinate life contingencies. We consider a univariate Gompertz diameter growth process facing stochastic fluctuations in the following Ito (1942) stochastic differential equation:
| (1) |
where α0, α1, β, σ>0 are unknown real parameters to be estimated, D(t) is a breast height diameter (in the sequel-diameter) at the age t, d0≥0, g(t) is an exogenous factor which is expressed by a time continuous known function, {W(t); tε[t0; T]} is a one-dimensional Wiener process and the differential dD(t) is to be understood in the sense of Ito (1942). In the sequel, the density p(d, t) of D(t) at t D(t0) = d0 at t = t0 is denominated as transition probability density function or conditional probability density function.
The height-age models can be used to detect trends in the exogenous factor g(t) because they indirectly show whether the growing conditions are changing over time (Skovsgaard and Vanclay, 2008; Garcia, 2009). In the sequel, we relate the exogenous factor g(t) as the height-age trajectory h(t) of a tree. Chronologies of height increments are a good tool to quantify the exogenous factor, thus is g(t) ≡ h(t). It is well-known that tree growth is sigmoidal (Garcia, 2005) and several sigmoidal growth models with biologically interpretable parameters have been proposed, such as Verhulst, Gomperz, Mitcherlich and Bertalanffy models. Recently, the majority of the newly developed age-height models are derived by the procedure named the algebraic difference approach (Cieszewski and Bella, 1989). In this study, we focus on two types of height-age models. The first height-age model is defined by the Gomperz shape growth model:
The formula describing the Gomperz height-age trajectory is defined in the following form:
| (2) |
where, γ, K are unknown real parameters to be estimated from the observations of a realization of a tree height for tε[t0; T], γ is the intrinsic growth rate of height, K is called the carrying capacity of the environment and commonly represents the maximum height that can be supported by the resources of the environment and h1(t0) = h0≥0. The second height-age model is defined by the Mitcherlich shape growth model:
The formula describing the Mitcherlich height-age trajectory is defined in the following form:
| (3) |
where, γ1, K1 are unknown real parameters to be estimated from the observations of a realization of a tree height for tε[t0, T], γ1 is the intrinsic growth rate of height, K1/γ1 is called the carrying capacity of tree height and h2 (t0), = h0≥0.
According to Eq. 1 the parameter β affects the term D(t) In (D(t)) and acts to slow down the drift term (α0+α1g(t))D(t). When β = 0, the diameter stochastic growth process D(t), tε[t0; T] described by Eq. 1 contains the univariate nonhomogeneous lognormal diffusion process with the exogenous factor g(t).
Using Itô’s formula for the age dependent transformation X(t) = eβt we obtain the explicit solution of original stochastic differential Eq. 1 in the following form:
| (4) |
which is a continuous nonhomogeneous Markov (Gaussian) process with transition probability density function (Gutierrez et al., 2008):
| (5) |
| (6) |
| (7) |
Noticing that the parameter α1 of the homogeneous case is equal to 0, the functions μ(t) and λ2(t) of the homogeneous case take the following forms:
Therefore, the random variable D(t)/D(t0) = d0 has one-dimensional lognormal distribution Λ(μ(t), λ2 (t)).
The mean and variance of the stochastic process D(t) defined by Eq. 1 take the following forms:
| (8) |
| (9) |
Next we address the approach of estimating the unknown parameters of stochastic
differential Eq. 1 from the following discrete sampling D(ti)
= di, h(ti) = hi, ti), i = 0,1,...,
n, assuming that ti-ti-l = 1. A natural estimation procedure
is maximum likelihood because it is possible to write the likelihood function
explicitly. Explicit knowledge of the transition probability density function
of diameter dynamic allows us to construct the likelihood function L(α0,
α1, β, σ). The transition probability density function
| (10) |
Derivation of the maximum likelihood function from Eq. 10 and the maximum likelihood estimators are given in the Appendix.
Goodness-of-fit tests allow us to verify the correspondence between the estimated theoretical model and real data set. The quantitative analysis of tree diameter distribution is usually based on the tests, such as, the Chi-squared, the Kolmogorov-Smirnov, the Anderson-Darling, the Cramer-von Mises (Thode, 2002). Most of these tests are very sensitive to the presence of outliers in the observed data. In forestry various measures for the deviation of an actual (empirical) distribution from its estimated theoretical distribution are commonly used, such as, the Reynolds error index, the absolute discrepancy, the stand stability index, the bias and standard error of estimate and many more (Reynolds et al., 1988; Cao, 2004). These measures of the goodness-of-fit can be used for comparisons between observed data sets and distribution models.
Statistical testing is often based on distributional assumption of normality. A useful technique for evaluating the normality of small and moderate size samples is the Shapiro-Wilk test statistic W (Shapiro and Wilk, 1965). In this study, we test the normality of the pseudo-residuals defined by Zucchini and MacDonald (1999). The pseudo-residuals, ri, corresponding to the observation (di, hi, ti) are defined in the following form:
| (11) |
where, Φ denotes the distribution function of the standard normal distribution,
(di, hi, ti) is the i-th observation of diameter,
height and age. Let
and Φ denotes the distribution function of a standard normal distribution. In case of residuals owning a standard normal distribution the value of statistics W tends to be close 1 and on the contrary tends to be small if residuals are from non-normal distribution. A normal probability plot of pseudo-residuals (11) is constructed by plotting r(i) against qi. The normal probability plot of pseudo-residuals enables us to evaluate visually the fit of the estimated theoretical diameter distribution to the observations.
In order to rank the performance of each transition probability density function we utilize Reynolds’ error index measure. The error index is calculated in 5 cm diameter classes for stem numbers. Thus, a relative error index (%) is defined by a sum of the absolute differences between the actual and predicted stem numbers of the diameter classes divided by the total stem number N:
| (12) |
where,
Growth data: The diameter analysis is based on measurements in pine (Pinus sylvestris) stands at Lithuania. The data were provided by the Lithuanian National Forest Inventory. We included full calliperings of permanent sample plots. Over 20 years period (1976-1996) in the even-aged uncut stand sample plots were re-measured at the most five times. The following variables were measured: age (t), number of trees per hectare, diameter at breast height (d), trees position (coordinates x, y), height (h) and descriptive variables such as alive or dead trees were also recorded. Approximately 20% of the sample trees were randomly selected for the height measurement. The measurements have been conducted in 30 occasions of permanent treatment plots and the initial planting densities are unknown. The age of stands ranges from 12 to 103 years. The diameter at breast height varies from 2.2 to 51.5 cm. Height was measured to the nearest 0.1 m with a digital height meter. Diameter was measured to the nearest 0.1 cm. For model estimation observations on 900 pines were used. The observed data sets of study plots are shown in Fig. 1 and 2. Figure 1 and 2 show the variation of diameter and height subject to age.
RESULTS AND DISCUSSION
Deterministic height-age models: Using the observed data set presented
in Fig. 1 and 2, were calculated the parameter
estimations of the exogenous curves h1(t), h2(t) defined
by Eq. 2-3. Notice that the original observed
data set was arranged by averaging the values observed in equidistant times.
Estimation of models 2-3 is achieved using Nonlinear Weighted Least Square method.
The values of the weighted least squares estimators (standard errors) are
Visual examination of the residuals versus predicted heights provided a random pattern around zero with approximately constant variance both for the Mitcherlich model and the Gompertz model (Fig. 3a, b). With the exception of some possible outliers the Mitcherlich and Gompertz models provide a good representation of the height data.
The distribution of a normal probability plot that is nearly linear suggests normal distribution of the standardized residuals. Figure 4a and b do not indicate any serious violation of the assumption of normality for standardized residuals. Typically, normal probability plots are not perfect straight. For the Gompertz and Mitcherlich height-age models the p-values of the Shapiro-Wilk (1965) statistic, W, are 0.0553 and 0.0772, respectively.
The exogenous height-age curves h1(t), h2(t), are presented in Fig. 5. For comparison, estimates for precision of the models were carried out based on the coefficient of determination (R2) and the relative error in prediction (RE%). The expressions of these statistics are defined by:
| (13) |
| (14) |
where,
Stochastic age-diameter models: Using the observed data set presented
in Fig. 1 and 2 were calculated the parameter
estimations of the stochastic diameter growth model defined by Eq.
1. We shall assume that the stochastic diameter growth process is observed
without error at a given collection of time instances
| Fig. 1: | The scatter graph of tree diameter at breast height against tree age for total set of sample trees n = 900 used for parameterization of growth equations |
| Fig. 2: | The scatter graph of tree height against tree age for total set of sample trees n = 900 used for parameterization of growth equations |
| Fig. 3: | Scatter plots of standardized residuals vs. predicted values: (a) the Gompertz exogenous model and (b) the Mitcherlich exogenous model |
| Fig. 4: | Normal probability plots of standardized residuals: (a) the Gompertz height-age model and (b) the Mitcherlich height-age model |
| Fig. 5: | Plot of the exogenous curves h1(t), h2(t) |
Figure 6a-c show the mean and standard deviation trajectories of the stochastic process D(t), tε[t0; T] of tree diameter. These functions are obtained by replacing the parameters in Eq. 8 and 9 by their estimators given in Table 1. All curves monotonically evolve to the steady state values. As we can see in Fig. 6a and b, the mean and standard deviation curves of tree diameter are very similar for both nonhomogeneous (Mitcherlich, Gompertz) models. The mean and standard deviation curves for the homogeneous model (Fig. 6c) describe a similar shape for trees less than 30 years age and subsequently get enlarged values than the nonhomogeneous ones.
Figure 7a -c show the estimated univariate transition probability density function (EDF) of tree diameter defined by Eq. 5. These density functions indicate that the EDF of tree diameter is steeper for the young stands and less steep for the mature stands. Figure 7a-c don’t fix marked difference between the EDFs using nonhomogeneous and homogeneous models.
| Table 1: | Parameter estimations |
| Fig. 6: | Plot of the mean and standard deviation dynamic of tree diameter with the parameterization data sets: (a) Mitcherlich exogenous factor; (b) Gompertzian exogenous factor and (c); homogeneous model; mean (continuous curve), mean±SD (non-continuous curve) |
| Fig. 7: | Plot of the estimated univariate transition probability density function: (a) for the Mitcherlich exogenous factor; (b) for the Gompertzian exogenous factor and (c) no exogenous; estimated surface of the age-dependent density (Eq. 5) of tree diameter (left); the density of tree diameter at the age t = 20, 60, 100 years (right) |
For the evaluation of goodness-of-fit of our presented lognormal shape univariate age-dependent transition probability density function (5) we use the Shapiro-Wilk statistic and normal probability plot. The normal probability plots of the pseudo-residuals using the estimates of parameters presented in Table 1 are shown in Fig. 8a-c. From Fig. 8a-b, it is possible to conclude that both EDFs for Mithcherlich and Gompertz exogenous factors fit not too bad. For the EDF with the Mitcherlich shape exogenous factor computed the Shapiro-Wilk statistic W yield a value 0.9792 (p-value 0.0010). For the EDF with the Gompertz shape exogenous factor computed the Shapiro-Wilk statistic W yield a value 0.9795 (p-value 0.0013). Finally, for the EDF of homogeneous model computed the Shapiro-Wilk statistic W yield a value 0.9792 (p-value 0.0008). These results lead us to a conclusion that the observed data set is compatible with the EDF (5) in all cases. It is worth remarking that the Shapiro-Wilk statistic provides a generally superior omnibus measure of non-normality. Moreover, the fitting data set was sufficiently large n = 900.
Finally, the relative error index was used in the comparisons as a measure of goodness of fit of the EDFs for the nonhomogeneous and homogeneous models. The values of the REI% measure Eq. 12 calculated for each EDF of stochastic nonhomogeneous (Mitcherlich, Gompertz) and homogeneous models were 23.00, 23.00, 21.44%, respectively. If we look at the relative error index from the age, the relative error index varies from 7 to 73% (Fig. 9). The relative error index is at its minimum at the age of 55 years. Taking into account that most of the stands covered in this study were within 12-80 years, the relative error index is a peaking function for ages greater than 80 years.
Application: The development of simple and accurate stand-specific volume
model based on easily obtainable tree and stand characteristics is a main problem
of forest mensuration. Traditionally, the mean tree volume
| (15) |
where, V(d, h) is an individual tree volume equation on diameter and height. Much greater accurateness is obtained by substituting (smoothing) a density function p(d) of tree diameter and integrating by all diameters d>0. If the tree volume regression function V(d, h) additionally depends on age (V(d, h, t)) and the density p(d) function additionally depends on age and height (p(d, h, t)), then Eq. 15 can be rewritten as follows:
| (16) |
The integral form Eq. 16 describes the mean tree volume
as an explicit function of height and age and can provide additional information
about volume dynamic. The commonly used functional dependence for volume (V(d,
h, t)) calculation takes the form of the power function
| Fig. 8: | Normal probability plots of pseudo-residuals: (a) Mitcherlich exogenous factor; (b) Gompertzian exogenous factor and (c) homogeneous model |
| Fig. 9: | Relative error index: Mitcherlich exogenous factor; Gompertzian exogenous factor and homogeneous model |
To derive the age-height-dependent density function p(d, h, t) of tree diameter for the nonhomogeneous stochastic model (1), we define Eq. 6 in the following form:
| (17) |
Substituting Eq. 17 into Eq. 5 gives the lognormal density function which depends on height and age as follows:
| (18) |
Figure 10a and b show the estimated lognormal density function of tree diameter which depends on height and age, defined by Eq. 18. These density functions are represented at heights 15 and 30 m. The higher is the height for a fixed age, the more skewed toward large diameters and the wider are densities curves. The ability to model age and height relationship in the diameter density function is often useful, especially in modeling diameter dynamic.
| Fig. 10: | Plot of the estimated lognormal density function Eq. 18 at height 15 m (a) and at height 30 m (b); estimated surface of the age-dependent density of tree diameter (left); the density of tree diameter at the age t = 20, 60, 100 years (right) |
The purpose of this contribution was to present a specific modeling approach based on the methodology of stochastic differential equations. To deal with the growth models in a numerical fashion, probabilistic means were adopted to give an understanding of the problems of the modeling of mean diameter, standard deviation of diameter and mean volume. Equations for predicting mean diameter (d) and standard deviation (s) of diameter are expressed in the general forms:
| (19) |
| (20) |
Thus the mean diameter Eq. 19, standard deviation of diameter Eq. 20 and mean volume Eq. 16 are modeled as a density-dependent set of curves. Figure 11a-c show the mean diameter of a tree, the standard deviation of tree diameter and the mean volume of a tree subject to height and age. These graphics demonstrate that the growth of mean diameter, standard deviation of diameter and mean volume is a peaking function over a single inflection point at diameter, increasing with age and height.
Using probabilistic mean diameter and volume growth models Eq.
19, 16 and height-age growth curves h = h(t) = h1(t)
or h = h(t) = h2(t) Eq. 2 and 3
we can define the current (the mean) annual diameter and volume increments
| (21) |
| (22) |
Relationships between the current annual diameter Eq. 21 and volume Eq. 22 increments against the height of a tree are illustrated in Fig. 12a and b. As we see in Fig. 12a and b, the height exerts a strong influence on current annual diameter and volume increments. The effect of a height on current annual diameter and volume increments becomes negligible above 100 year age.
Figure 13a-d shows the current and mean annual diameter and volume increments against the age and the mean diameter of a tree using the Mitcherlich height-age growth curve h = h(t) = h2(t) (Eq. 3). From Fig. 13a and c we see that the culmination of volume increment is reached even later than that of diameter increment. The peak in current and mean annual diameter increments occurred at 21 and 81 years of age (Fig. 13a-c), respectively and current and mean annual volume increments peaked at 26 and 125 years of age (Fig. 13a-c), respectively. If an objective of forest management is to maximize the produced stem volume, the trees should be retained until they attain their maximum mean annual volume increment at the age 125 years (Fig. 13c) or at the mean diameter 42.2 cm (Fig. 13d). The mean annual diameter increment is greatest at the mean diameter 13.3 cm (Fig. 13b).
| Fig. 11: | Plot of the mean diameter, standard deviation of diameter and mean volume of a tree: the height-age-dependent response surface of mean diameter (a); standard deviation of diameter (b); mean volume of a tree (c); the age-dependent mean diameter, standard deviation of diameter and mean volume of a tree (left); at heights 15 and 30 m (right) |
| Fig. 12: | Relationship between the current annual diameter and volume increments of an average tree against the height of a tree: (a) current diameter increment and (b) current volume increment at the height 20, 10 and 5 m |
| Fig. 13: | Annual diameter and volume increments: (a) annual diameter increments against the age; (b) annual diameter increments against the mean diameter; (c) annual volume increments against the age and (d) annual volume increments against the mean diameter |
The interrelations of the current annual volume increment and mean annual volume increment curves of a stand and the position of their point of intersection are of particular interest to forest management. Enlarged understanding and statistical inference in stand current and mean annual volume increment models require an adequate representation of the prediction of tree mortality (survival).
CONCLUSION
Given the importance of stochastic analysis in modern forestry, we consider the case where the governing tree diameter dynamic is defined by an elementary stochastic differential equation. A theoretical prerequisite of our presented approach was the stochastic Gompertz diameter growth law driven by one-dimensional standard Wiener process. The results obtained here have shown that it is possible to relate nonlinear stochastic diameter growth law and diameter distribution law. For a realistic representation of diameter and height growth, was used Gompertzian and Mitcherlich growth models.
Thus, the proposed method could be continued in terms of properly modifying the drift and diffusion functions of the stochastic diameter growth process and choosing exogenous factors.
The accuracy of the age-height-dependent diameter distribution (Eq. 5) depends on the amount of information available from the stand. Our methodology extends some way to inclusion of the basal-area or/and density of a stand as an exogenous factor or as an independent variable.
APPENDIX
The maximum likelihood estimates: Here, section we collect some results which were used in order to estimate the model parameters. To write the maximum likelihood function explicitly is possible because the transition probability density function of the diameter stochastic process D(t), tε[t0; T] is explicitly solved by Eq. 5.
The resulting maximum likelihood function is defined by:
| (A1) |
Where:
| (A2) |
| (A3) |
| (A4) |
| (A5) |
Uβ is 2xn matrix defined by
| (A6) |
Thanks to quadratic form in α the maximum likelihood estimators of α and σ2 are given by (Gutierrez et al., 2008):
| (A7) |
| (A8) |
Where:
| (A9) |
While the likelihood estimators of α and σ2 for stochastic Gompertz process (1) is well established by Eq. A7, A8, estimating parameter β is not straightforward. In this paper, the estimation approach we follow is to first estimate by likelihood estimation procedure the parameter β of the ordinary differential Eq. 1 (σ = 0), which represents the deterministic part of the stochastic Gompertz process (1). The maximum likelihood estimator of the parameter β is defined in the following form:
| (A10) |
where, yi = ln di.