Stochastic Mixed-Effects Parameters Bertalanffy Process , with Applications to Tree Crown Width Modeling

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 crownwidth 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 fixedeffects 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.


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 databased 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 2 Mathematical Problems in Engineering 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][10][11][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 mixedeffects 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 mixedeffects 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.

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  ≥ 0 > ln()/,  > 0,  > 1.The formula describing the Bertalanffy trajectory follows the form of a sigmoidal function where  0 = ( 0 ).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 fixedeffects parameter to produce a stand-specific parameter in the following form: where parameters   of the th stand ( = 1, 2, . . ., ) are stand-specific random effects.It is assumed that the random effects   ( = 1, 2, . . ., ) are independent and normally distributed with a mean of 0 and a constant variance (  ∼ (0;  2  )).Therefore, the response variable,   (),  ≥ 0,  = 1, 2, . . ., , is described using a stochastic differential equation of the Bertalanffy type: where   (),  ≥ 0, are the independent standard Brownian Motions,   () and   are assumed to be mutually independent for all 1 ≤ ,  ≤ , 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  = 1, 2, . . .,  are termed exogenous stand variables and constitute a diameter-continuous function in  ∈ [ 0 ;  0 ],  0 ,   are diameter-independent real fixedeffects 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  ≡ ln() transforms ( 8) to the following form: Then, by integration and  ≡ exp(), we have Taking into account that the random variable ∫   0 () has a normal distribution with mean 0 and variance  −  0 we can deduce that the process   () corresponds to a lognormal distribution, Λ 1 (  ( |  0 ,  0 ); V( |  0 )).In this case, the conditional probability density function for the considered stochastic process   () is where In this paper, we approximate the functions   (),  = 1, 2, . . ., , by polygonal functions.We have real observations   0 ,    ;  = 1, 2, . . .,   (  is the number of observations possessed in the th stand) for exogenous stand variables,   , in the values of predictor variables  0 and    .Hence, for  = 1, 2, . . .,  we have More often than not each exogenous stand variable in a particular stand has the same observed value,    =    ;  = 1, 2, . . .,   .From ( 14), we deduce that analytical expression of the integral to the right-side equation ( 12) is where where  2 () is the dilogarithm function defined by  2 () = ∫  0 (ln(1 − )/).The conditional mean, mode, median, quantile, variance, and coefficient of variation functions, respectively, (⋅), (⋅), (⋅), (⋅) (0 <  < 1), and (⋅), of the stochastic crown width process are given by The integral ∫  2    (  | 0 , 1 ,...,  ,,,,  ) ⋅   ,  = 1, 2, . . ., , by a second-order Taylor series expansion can be approximated as the Laplace approximation [18]: after plugging the φ into (29).

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   ,  = 1, 2, . . ., , is taken from a new stand, the randomeffects parameter  for the new stand can be predicted in the following form: where α0 , β, γ 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 α0 .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, l, 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, Φ −1  (0; σ2  ) is the inverse of the normal distribution with a mean of 0 and a constant variance σ2  for a uniform random variable, , in the interval [0.05; 0.95], μ ( |  0 ,  0 ) and V ( |  0 ) are the estimated trend of the mean and variance (calculated by ( 12) and ( 13)) of the lognormal density of the crown width, respectively, and  −1  ( μ ( |  0 ,  0 ); V ( |  0 )) is the inverse of the lognormal distribution with a mean of μ ( |  0 ,  0 ) and a variance of V( |  0 ) for a uniform random variable, , in the interval [0.05; 0.95].

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 ( 2 ), Akaike's information criteria [19], and the Bayesian information criterion [20]: where LL( θ) 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 loglikelihoods of the models is defined by where the test statistic  is asymptotically distributed as a chisquared random variable  2  2 − 1 with degrees freedom  2 − 1 , where  1 and  2 are the number of free parameters and θ1 , θ2 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].

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.

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 ( < 0.05).
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 fixedand mixed-effects SDE crown width models the mean prediction bias (the percentage mean prediction bias) proves  to be in intervals −0.0504 m-−0.0477m (−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.
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 ( = 1.0;  = 0.0).
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 firststage 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 mixedeffects 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.

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.

Figure 1 :
Figure 1: Residuals against the predicted crown width and lowess regression curve for all fitted crown width SDE Bertalanffy type models: fixed-effects models (a) and mixed-effects models (b).

Table 2 :
Estimated parameters (standard errors) of all the used fixed-effects and mixed-effects models applied to the estimation dataset.

Table 3 :
Estimation goodness of fit statistics for all the used fixed-effects models applied to the estimation and validation datasets * .The best values of the performance statistics for all scenarios of stand variables are in bold, * * the mean prediction bias  = (1/)∑  =1 (  −   − )2), and  is the number of parameters in the model,  = ∑  =1   is the total number of observations used to fit the model,  is the number of stands,   is the number of trees in th stand, and   , ∧   , and  are the measured, estimated, and average values of the dependent variable (crown width). *