Abstract
A stochastic modeling approach based on the Bertalanffy law gained interest due to its ability to produce more accurate results than the deterministic approaches. We examine tree crown width dynamic with the Bertalanffy type stochastic differential equation (SDE) and mixed-effects parameters. In this study, we demonstrate how this simple model can be used to calculate predictions of crown width. We propose a parameter estimation method and computational guidelines. The primary goal of the study was to estimate the parameters by considering discrete sampling of the diameter at breast height and crown width and by using maximum likelihood procedure. Performance statistics for the crown width equation include statistical indexes and analysis of residuals. We use data provided by the Lithuanian National Forest Inventory from Scots pine trees to illustrate issues of our modeling technique. Comparison of the predicted crown width values of mixed-effects parameters model with those obtained using fixed-effects parameters model demonstrates the predictive power of the stochastic differential equations model with mixed-effects parameters. All results were implemented in a symbolic algebra system MAPLE.
1. Introduction
Ordinary differential equation is a powerful tool for analyzing scientific data by model-driven approach. Many real-life phenomena can be described by a nonlinear ordinary differential equation [1, 2]. Usually, the nonlinear model parameters are easier to interpret compared to those from linear regression models because the parameters in nonlinear models generally have a natural physical interpretation. The dynamics of biological systems are largely driven by stochastic processes and are subject to random external perturbations. A stochastic differential equation (SDE) is a differential equation in which some of the terms evolve according to Brownian Motion. This approach assumes that the dynamic is partly driven by noise. SDEs are often used to model the stochastic dynamics of biological systems [3].
Tree growth dynamical models project the growth and development of forest ecosystems, known to be stochastic in nature [4]. In order to understand this complex dynamics behavior computational modeling is inevitable. As an alternative to model-based approach, data-based methods have been developed during the past several decades. The recent surveys by Yin et al. [5, 6] provide the excellent review of the current development of the advanced and sophisticated data-based methods, schemes, and their applications. However, traditionally used multivariate statistical methods [7] are not suitable to handling the process dynamics of forest ecosystems and operation condition changes.
Tree crown width models can be classified into two basic types according to their independent variables [8]. The first type (local model) assumes that tree crown width is completely dependent on the breast height diameter. The second type (generalized model) includes the breast height diameter and other individual tree variables (exogenous variables) such as tree crown height, total tree height. The first model type requires only low sampling effort and is usually locally applied, whereas the second model type demands high sampling effort and is often applied regionally. All the methods traditionally used in nonlinear crown regression models ignore the variation of residuals that assumes independence and constant variance.
In this study, we use crown width as a test case and examine how the behavior of a small-size tree is altered when stochasticity is introduced into the deterministic Bertalanffy type equation. To detect an explanation about the behavior of the crown width by the diffusion a set of external factors (exogenous variables) is introduced in the model. The breast height diameter and exogenous stand variables behaviors are assumed to be known and they must contribute to the description of the evolution of the stochastic crown width process. Several SDE models with fixed-effects parameters [4, 9–12] and mixed-effects parameters [13] have previously been developed for modeling the tree diameter and height. In this paper, more sophisticated Bertalanffy type stochastic process of a response variable is proposed including mixed-effects parameters and nonhomogeneous drift coefficient which depends on deterministic functions that describe the dynamics of certain stand variables named as exogenous variables. The importance of the mixed-effects parameters SDE Bertalanffy type model lies in its ability to split the total variation into within- and between-stands components. The new developed model describes the within-stand variation in data through two sources of noise: the measurement noise representing the uncorrelated part of the residual variability associated with assay or sampling errors and the system noise reflecting the random fluctuations around the corresponding theoretical growth model. If the magnitude of the parameter of the system noise is zero, then the entire system noise term will vanish and the remaining part of the SDE will simply be the ordinary differential Bertalanffy form whose solution is the regression term of the mixed-effects model. Maximum likelihood estimations of the fixed- and mixed-effects parameters are calculated using Laplace approximations of the maximum likelihood. The maximum likelihood methodology developed allows us to statistically fit the model to the sample dataset considered.
The objective of this paper is to develop a novel SDE Bertalanffy type model for growth modeling and to put forward the advantages of using diffusion processes, random parameters, and exogenous stand variables in the analysis of crown width dynamics. We also discuss how conditional density function can be used to construct maximum likelihood estimators of the fixed- and mixed-effects parameters. A MAPLE program was implemented to carry out the calculations required for this study.
2. Material and Methods
2.1. Stochastic Differential Equations Model
The focus of the present work is on the dynamics of the response variable (crown width), , as a stochastic process with respect to the predictor variable (diameter at breast height − diameter in the sequel), . In this study, we selected using the generalization of the Bertalanffy type ordinary differential equation. The von Bertalanffy [14] hypothesized that the growth of an organism could be represented as the difference between the synthesis and degradation of its building materials. There are few theoretical equations formulated specifically for biology applications. In this paper, we use deterministic model developed by Román-Román et al. [15] as the basis of the newly developed stochastic model. The changes in response variable, , are described using the ordinary differential equation where , , and are unknown fixed-effects parameters and verify , , . The formula describing the Bertalanffy trajectory follows the form of a sigmoidal function where . Equation (1) can be seen as a generalization of the Malthusian growth model with nonconstant fertility depending on the predictor variable (diameter) by . There are alternative ways of introducing stochasticity into the behavior of the response variable. In this work, the randomness in the operation of response variable was defined by standard Brownian Motion. Therefore, the complete deterministic model defined by (1) was converted into a stochastic model assuming that the fertility varies randomly around the mean where is the diffusion coefficient and is a Gaussian white noise process. The relationship between the response and predictor variables is altered by environmental conditions. Stand-specific characteristics, such as soil type, nutrient status, and elevation, cause parameters to differ between different stands. In the case of between-stand variation, the parameters , , and vary from stand to stand and hence account for this variation. The interest of this study is the development of growth models for a large geographic region rather than localized areas. In order to compensate for the randomness found in stands we can add a random term to our model. Thus, specific stands may have what are generally termed “random parameters” in mixed-effects model terminology. For the construction of a mixed-effects model, the model first needs to determine which parameters should be considered mixed and which should be considered purely fixed. The parameter with high variability could be considered as mixed-effects parameter. The parameter exhibits high variation between stands and thus can be altered by adding stand-specific random effects to the fixed-effects parameter to produce a stand-specific parameter in the following form: where parameters of the th stand are stand-specific random effects. It is assumed that the random effects are independent and normally distributed with a mean of 0 and a constant variance . Therefore, the response variable, , , , is described using a stochastic differential equation of the Bertalanffy type: where , , are the independent standard Brownian Motions, and are assumed to be mutually independent for all , , and is the total number of stands used for model fitting.
Moreover, to account for the variation between different growth conditions, the parameter can be modified by including function that introduces the dynamics of each exogenous stand variable into the model. The function is considered as a linear combination of the exogenous stand variables and is given by where for are termed exogenous stand variables and constitute a diameter-continuous function in , are diameter-independent real fixed-effects parameters (to be estimated), and is the number of exogenous stand variables. One candidate set for the basis functions that is widely used is polygonal functions. Therefore, the multifunctional nature of the parameter is established in terms of a random-effect parameter and a linear combination of exogenous stand variables functions, which are deterministic diameter functions, in the following form: Equation (7) can be put into the above formulae (5) to determine the effects of the covariates (exogenous stand variables) and random effects on response variable. In this situation, the nonhomogeneous SDE of the response variable is expressed as follows: By Itô’s [16] lemma, the logarithmic transformation transforms (8) to the following form: Then, by integration and , we have Taking into account that the random variable has a normal distribution with mean 0 and variance we can deduce that the process corresponds to a lognormal distribution, . In this case, the conditional probability density function for the considered stochastic process is where In this paper, we approximate the functions , , by polygonal functions. We have real observations ; ( is the number of observations possessed in the th stand) for exogenous stand variables, , in the values of predictor variables and . Hence, for we have More often than not each exogenous stand variable in a particular stand has the same observed value, ; . From (14), we deduce that analytical expression of the integral to the right-side equation (12) is where where is the dilogarithm function defined by . The conditional mean, mode, median, quantile, variance, and coefficient of variation functions, respectively, , , , , and , of the stochastic crown width process are given by
2.2. Computation of the Parameter Estimators
We consider the generalized SDE crown width (response variable) model defined by (8) in two aspects. First, the maximum log-likelihood function is derived for generalized fixed-effects model (in this case the parameters of random effects, , , are assumed to be equal to its mean value ). Second, the maximum log-likelihood function is derived for generalized mixed-effects model. The fixed-effects parameters and random-effect parameters , , are estimated by means of the maximum likelihood procedure using discrete sampling and conditional probability density function defined by (11). We assume that all observations are independent (no repeated measurements are used in the dataset for model fitting). Let us consider a discrete sample of the crown width process () at the diameters (), where is the number of observed trees of the th plot, . Under the initial condition , the associate log-likelihood function for generalized fixed-effects parameters model defined by (8) can be obtained by the following expression: where is the number of observed trees of the th plot , conditional density function takes the form defined by (11), and random effects ; . The maximum log-likelihood function for generalized mixed-effects parameters model defined by (8) takes the following form: where are fixed-effects parameters (the same for all stands) and are random-effects parameters (stand-specific), which are assumed to follow a normal distribution with mean 0 and standard deviation , is the normal density of the random effects. Unfortunately, the integral in (24) has no closed form solution. Since analytic expression for the integrand in (24) is known, the Laplace method may be used [17]. Let us define a function by The integral , , by a second-order Taylor series expansion can be approximated as the Laplace approximation [18]: where The log-likelihood function for a generalized mixed-effects SDE model defined by (8) is approximately given by The maximization of is a nested optimization problem. The internal optimization step estimates the for every stand by (27). The external optimization step maximizes after plugging the into (29).
2.3. Calibration
In forestry literature, calibration requires the prediction of the random-effects parameter using a supplementary sample of observations collected at the same sampling unit. The crown width of trees in new stand can be predicted either by setting the random effects to zero or by adding random parameter predicted from prior observations. When the diameter, , and crown width, , of a subsample of trees are known, the predicted random parameter is added to the fixed parameter to obtain localized parameter for the corresponding stand. If a subsample (without exogenous stand variables) of trees with crown width and diameter , , is taken from a new stand, the random-effects parameter for the new stand can be predicted in the following form: where are estimations of the parameters calculated by the maximum likelihood procedure for mixed-effects parameters model. The crown width of another tree from the same stand can be estimated by adding the random effects predicted by (30) to parameter . Mixed-effects parameters SDE model incorporates the variability between stands through the expression of the model’s parameters and in terms of both fixed and random effects. Random effects are conceptually random variables and can be simulated as such in terms of their distribution. To address this, a random component to the random-effects parameter prediction, , and the crown width predictions, , can additionally be added. This stochastic approach uses the distribution functions and confidence intervals of random variables, , and . The stochastic predictions of and can be defined in the following form, respectively: where is the estimation value of random effects obtained by (30), is the estimation value of the standard deviation of random effects, is the inverse of the normal distribution with a mean of 0 and a constant variance for a uniform random variable, , in the interval , and are the estimated trend of the mean and variance (calculated by (12) and (13)) of the lognormal density of the crown width, respectively, and is the inverse of the lognormal distribution with a mean of and a variance of for a uniform random variable, , in the interval .
2.4. Statistical Analysis
Numerical and graphical analyses of the residuals were used as criteria for comparing the generalized fixed-effects, mixed-effects SDE crown width models. The performance statistics of the SDE crown width models included seven fit statistics: the mean prediction bias (B), the percentage mean prediction bias (PB), the root mean squared error (RMSE), the percentage root mean squared error (PRMSE), an adjusted coefficient of determination (), Akaike’s information criteria [19], and the Bayesian information criterion [20]: where is the maximum log-likelihood function associated with model, is the maximum likelihood estimate of the parameters of the model, and is the number of parameters in the model.
In practice, many stand exogenous variables are introduced to reduce the possible model deviation. Generally speaking, the larger the number of introduced exogenous stand variables is, the less the crown width biases are. However, many of them may be insignificant and should be excluded from the final model to increase the accuracy of prediction. So it is necessary to study the exogenous stand variable selection procedure. Additionally, we used a statistical test defined by the ratio of the maximum log-likelihood of the model to determine whether exogenous stand variables are significant. The test statistic of the ratio of the maximum log-likelihoods of the models is defined by where the test statistic is asymptotically distributed as a chi-squared random variable with degrees freedom , where and are the number of free parameters and , are the maximum likelihood parameter estimates for models 1 and 2, respectively.
To assess the standard errors of the maximum log-likelihood estimators for the SDE Bertalanffy type models, a study of the Fisher [21] information matrix was performed. The asymptotic variance of the maximum likelihood estimator is given by the inverse of the Fisher’ information matrix [21].
3. Results and Discussion
3.1. Data
We focus on the modeling of Scots pine (Pinus Sylvestris L.) tree dataset. Scots pine tree stands dominate Lithuanian forests, grow on Arenosols and Podzols forest sites,and cover 725,500 ha. At plot establishment, the following data were recorded for every sample tree: crown width, , diameter over bark at 1.30 m, , and crown height, . A random sample of 12 plots (1630 trees) was selected for model estimation, and the remaining dataset of 5 plots (698 trees) was utilized for model validation. Summary statistics for the breast height diameter over bark (diameter), , crown height, , and crown width, , of all the trees used for model estimation and validation are presented in Table 1.
3.2. Performance of SDE Models
To examine the effect of fixed, random parameters and exogenous stand variables on crown width predictions, the parameters were estimated by the maximum likelihood method using the NLPSolve procedure in a mathematical software package MAPLE [22]. A large estimation dataset allowed us to obtain precise estimates of all fixed-effects and random-effects parameters as well as their standard errors. Estimation results are presented in Table 2. All parameters are highly significant ().
Table 3 lists the fit statistics for all the new developed fixed- and mixed-effects parameters and without and with exogenous stand variable SDE crown width models for both estimation and validation datasets. On the whole, for fixed- and mixed-effects SDE crown width models the mean prediction bias (the percentage mean prediction bias) proves to be in intervals −0.0504 m–−0.0477 m (−6.70%–−4.74%) and 0.0350 m–0.0374 m (−3.41%–−2.68%), respectively, for the estimation dataset. For the estimation dataset the percent of variation explained attains high levels 65.84% and 73.66%, respectively. For the validation dataset the crown width’s estimate for the fixed- and mixed-effects SDE models proves satisfactory, with the mean prediction biases (the percentage mean prediction bias) −0.0236 m–0.1117 m (−6.21%–−5.18%) and 0.0344 m–0.1115 m (−2.01%–−0.70%), respectively. For the validation dataset the percent of variation explained attains high levels too, 70.66% for the fixed-effects methodology and 73.19% for the mixed-effects methodology.
The fixed- and mixed-effects parameters SDE Bertalanffy type models with exogenous stand variable (crown height) showed no difference of the fit statistics for the validation dataset (see Table 3). The statistical significance of the difference between two models can be assessed by the simple test defined by (33). For the estimation dataset the values calculated under the appropriate chi-squared distribution showed that the generalized fixed- and mixed-effects SDEs Bertalanffy type models including crown height, as exogenous stand variable, fit the estimation dataset significantly better than the SDEs models without exogenous stand variable, and we infer that this additional stand variable is biologically meaningful.
Another way to evaluate and compare the SDE crown width models is to look at the graphics of the residuals with lowess regression for the estimation dataset. The residuals are the differences between measured and predicted crown widths. Positive residuals mean underestimation and negative residuals mean overestimation. The plots of the residuals do not indicate any serious tendency towards overestimation or underestimation of predicted crown width values (see Figure 1) for the fixed- and mixed-effects parameter models. Both mixed-effects parameter models also had an approximately homogeneous variance over the full range of the predicted crown width values, as well as no systematic pattern in the variation of the residuals (see Figure 1). In the estimation dataset, the greatest prediction errors were obtained for larger diameter classes. It was also found that the overestimation (negative prediction error) increases slightly with crown widths, in particular, for values higher than 5 m. This can be due to the very few observations of a crown width of more than 5 m in our estimation dataset.
(a)
(b)
In order to illustrate the behavior of the new developed mixed-effects SDE Bertalanffy type models, crown width versus diameter was plotted in Figure 2 for all the validation plots and without exogenous (stand) variable. Crown width was calculated using the mixed-effects model defined by (8) without exogenous variable. The calibrated value of the random parameter for a stand present in the validation dataset was calculated by (30). The crown width-diameter relationships developed in this study are constrained to pass through the point (; ).
The models traditionally used in nonlinear mixed-effects mode of regression analysis are supplemented by a model for the between-stand variation in the model parameters and a model for the variation of the residuals that assumes independence and constant variance. However, the variance over the full range of the predicted values is not homogeneous, and it may lead to inexact estimates. New developed nonlinear mixed-effects models based on SDE extend the usual nonlinear mixed-effects regression models by including system noise as an additional source of variation in the first-stage model. This extended model describes the within-stand variation in data through two sources of noise: the measurement noise representing the uncorrelated part of the residual variability associated with assay or sampling errors and the system noise reflecting the random fluctuations around the corresponding theoretical crown width model. If the magnitude of the parameter of the system noise is zero, then the entire system noise term will vanish and the remaining part of the SDE will simply be the ordinary differential form whose solution is the regression term of the mixed-effects model. Regarding stochastic differential equations, as far as I know, in tree crown width modeling the generalized mixed-effects parameters model methodology has not previously been applied. Unfortunately, mixed-effects parameters SDE models can be implemented using a symbolic calculus program for computing the analytic expressions for the Hessians.
4. Conclusions
The new Bertalanffy type crown width models were developed using stochastic differential equation. The SDE method provides more mathematical sophisticated and narrow examination analysis tools compared to regression approaches. Comparison of the predicted crown width calculated using SDE models with the values obtained using the existing regression models revealed a comparable predictive power of the SDE crown width model with mixed-effects parameters.
The developed SDE crown width models may be recommended for both their ease of fitting procedure and the biological interpretation of the relevant parameters.
The SDE approach allows us to incorporate new tree and stand variables and new forms of stochastic dynamic.
The variance functions developed here can be applied to generate weights in every linear and nonlinear least squares regression crown width model by the weighted least squares form.
Finally, SDE methodology may be of interest in diverse areas of research that are far beyond the modeling of tree crown width.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors would like to express thanks to the unknown referees for their careful reading and helpful comments.