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

Authors: Petras Rupšys
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 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),Y, as a stochastic process with respect to the predictor variable (diameter at breast height − diameter in the sequel), x. 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, y(x), are described using the ordinary differential equation
(1)dy(x)dx=αβγeβx-γyx,
where α, β, and γ are unknown fixed-effects parameters and verify d≥0>ln⁡γ/β, β>0, α>1. The formula describing the Bertalanffy trajectory follows the form of a sigmoidal function
(2)yx=y01-γe-βx1-γe-βx0α,x∈x0;X0,
where y0=y(x0). Equation (1) can be seen as a generalization of the Malthusian growth model with nonconstant fertility depending on the predictor variable (diameter) by r(x)=αβγ/(eβx-γ). 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
(3)rx=αβγeβx-γ+σεx,
where σ is the diffusion coefficient and ε(x) 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:
(4)α=α0+ϕi,
where parameters ϕi of the ith stand (i=1,2,…,M) are stand-specific random effects. It is assumed that the random effects ϕi(i=1,2,…,M) are independent and normally distributed with a mean of 0 and a constant variance (ϕi~N0;σϕ2). Therefore, the response variable, Yix, x≥0, i=1,2,…,M, is described using a stochastic differential equation of the Bertalanffy type:
(5)dYix=α0+ϕiβγeβx-γYixdx+σYixdWix,PYix0=y0=1,x∈x0;X0,
where Wix, x≥0, are the independent standard Brownian Motions, Wi(x) and ϕj are assumed to be mutually independent for all 1≤i, j≤M, and M 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 g(x) that introduces the dynamics of each exogenous stand variable into the model. The function g(x) is considered as a linear combination of the exogenous stand variables and is given by
(6)gx=α0+∑k=1qαkgkx,
where gk(x) for k=1,2,…,q are termed exogenous stand variables and constitute a diameter-continuous function in x∈x0;X0, α0,αk are diameter-independent real fixed-effects parameters (to be estimated), and q is the number of exogenous stand variables. One candidate set for the basis functions gk(x) 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:
(7)αx=α0+ϕi+∑k=1qαkgkx.
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:
(8)dYix=βγα0+ϕi+∑k=1qαkgkxeβx-γYixdx+σYixdWix,PYix0=y0=1,x∈x0;X0.
By Itô’s [16] lemma, the logarithmic transformation Z≡ln⁡(Y) transforms (8) to the following form:
(9)dZix=βγα0+ϕi+∑k=1qαkgk(x)eβx-γ-σ22dx+σdWix,PZix0=ln⁡(y0)=1.
Then, by integration and Y≡exp⁡(Z), we have
(10)Yix=y0exp⁡∫x0xβγα0+ϕi+∑k=1qαkgk(u)eβu-γdu-σ22x-x0+σ∫x0xdWiu.
Taking into account that the random variable ∫x0xdWiu has a normal distribution with mean 0 and variance x-x0 we can deduce that the process Yi(x) corresponds to a lognormal distribution, Λ1μix∣y0,x0;v(x∣x0). In this case, the conditional probability density function for the considered stochastic process Yi(x) is
(11)fiy,x=1y2πv(x∣x0)×exp⁡-12v(x∣x0)ln⁡y-μix∣y0,x02,
where
(12)μix∣y0,x0=ln⁡y0βγ1-γe-βx1-γe-βx0α0+ϕi+βγ∑k=1qαk∫x0xgk(u)eβu-γdu-σ22x-x0,(13)vx∣x0=σ2x-x0.
In this paper, we approximate the functions gk(d), k=1,2,…,q, by polygonal functions. We have real observations z0ki,zjki; j=1,2,…,ni (ni is the number of observations possessed in the ith stand) for exogenous stand variables, zk, in the values of predictor variables x0 and xjki. Hence, for k=1,2,…,q we have
(14)gkx=z0ki+zjki-z0kixjkix-x0,0≤x≤xjki,j=1,2,…,ni.
More often than not each exogenous stand variable in a particular stand has the same observed value, zjki=zki; j=1,2,…,ni. From (14), we deduce that analytical expression of the integral to the right-side equation (12) is
(15)∫x0xgk(u)eβu-γdu=1βγxGkx,
where
(16)Gkx=z0kx-x0zk-z0kln⁡1-γe-βx1-γe-βx0-zj-z02βx2-x02+xln⁡1-eβxγ-x0ln⁡1-eβx0γ-2Li2eβxγ-Li2eβx0γ,
where Li2x is the dilogarithm function defined by Li2x=∫0x(ln⁡1-u/u)du. The conditional mean, mode, median, quantile, variance, and coefficient of variation functions, respectively, m(·), mo(·), me(·), mq(·)(0<q<1), and w(·), of the stochastic crown width process are given by
(17)mix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡⁡αkGkx,(18)moix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡⁡αkGkx·exp⁡-32σ2x-x0,(19)meix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡⁡αkGkx·exp⁡-12σ2x-x0,(20)mqix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡(αkGk(x))·exp⁡-12σ2x-x0+σx-x0z1-q,(21)wix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡(αkGk(x))2·exp⁡σ2x-x0-1,(22)cvx=exp⁡σ2x-x0-1.


### 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, ϕi, i=1,…,M, are assumed to be equal to its mean value Eϕi=0). Second, the maximum log-likelihood function is derived for generalized mixed-effects model. The fixed-effects parameters α0,α1,…,αq,β,γ,σ,σϕ and random-effect parameters ϕi, i=1,…,M, 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 (l1i,l2i,…,lnii) at the diameters (d1i,d2i,…,dnii), where ni is the number of observed trees of the ith plot, i=1,2,…,M. Under the initial condition PL0=1.0=1, the associate log-likelihood function for generalized fixed-effects parameters model defined by (8) can be obtained by the following expression:
(23)LL1α0,α1,…,αq,β,γ,σ=∑i=1M∑j=1niln⁡flji,dji,
where ni is the number of observed trees of the ith plot i=1,2,…,M, conditional density function fl,d takes the form defined by (11), and random effects ϕi≡0; i=1,2,…,M. The maximum log-likelihood function for generalized mixed-effects parameters model defined by (8) takes the following form:
(24)LL2(α0,α1,…,αq,β,γ,σ,σϕ)=∑i=1M∫R∑j=1niln⁡fi(lji,dji)+ln⁡p(ϕi∣σϕ)·dϕi,
where α0,α1,…,αq,β,γ,σ,σϕ are fixed-effects parameters (the same for all stands) and ϕi are random-effects parameters (stand-specific), which are assumed to follow a normal distribution with mean 0 and standard deviation σϕ, p(ϕi∣σϕ) 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 g:R→R by
(25)giϕi∣α0,α1,…,αq,β,γ,σ,σϕ=∑j=1niln⁡fi(lji,dji)+ln⁡p(ϕi∣σϕ),i=1,2,…,M.
The integral ∫R2egiϕi∣α0,α1,…,αq,β,γ,σ,σϕ·dϕi, i=1,2,…,M, by a second-order Taylor series expansion can be approximated as the Laplace approximation [18]:
(26)ln⁡∫R2egiϕi∣α0,α1,…,αq,β,γ,σ,σϕ·dϕi≈giϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ+12ln⁡2π-12ln⁡-Hϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ,
where
(27)ϕ^i=argmaxϕigiϕi∣α0,α1,…,αq,β,γ,σ,σϕ,(28)Hi(ϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ)=∂2giϕi∣α0,α1,…,αq,β,γ,σ,σϕ∂2ϕiϕi=ϕ^i.
The log-likelihood function for a generalized mixed-effects SDE model defined by (8) is approximately given by
(29)LL2(α0,α1,…,αq,β,γ,σ,σϕ)≈∑i=1Mgiϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ+12ln⁡2π-12ln⁡-Hiϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ.
The maximization of LL2(α0,α1,…,αq,β,γ,σ,σϕ) is a nested optimization problem. The internal optimization step estimates the ϕ^i for every stand i=1,2,…,M by (27). The external optimization step maximizes LL2(α0,α1,…,αq,β,γ,σ,σϕ) after plugging the ϕ^i 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,d, and crown width, l, 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 r trees with crown width li and diameter di, i=1,2,…,r, is taken from a new stand, the random-effects parameter ϕ for the new stand can be predicted in the following form:
(30)ϕ^=1r∑i=1rln⁡1-γ^exp⁡-β^di1-γ^exp⁡-β^d0-1ln⁡li-ln⁡l0×ln⁡1-γ^exp⁡-β^di1-γ^exp⁡-β^d0-1-α^0,
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 L(d). The stochastic predictions of ϕ and L(d) can be defined in the following form, respectively:
(31)ϕstoch=ϕ^+ΦU-10;σ^ϕ2,lstoch=LNU-1μi^(d∣l0,d0);v^(d∣d0),
where ϕ^ is the estimation value of random effects obtained by (30), σ^ϕ is the estimation value of the standard deviation of random effects, ΦU-10;σ^ϕ2 is the inverse of the normal distribution with a mean of 0 and a constant variance σ^ϕ2 for a uniform random variable, U, in the interval 0.05;0.95, μi^(d∣l0,d0) and v^(d∣d0) are the estimated trend of the mean and variance (calculated by (12) and (13)) of the lognormal density of the crown width, respectively, and LNU-1μi^(d∣l0,d0);v^(d∣d0) is the inverse of the lognormal distribution with a mean of μi^(d∣l0,d0) and a variance of v^(d∣d0) for a uniform random variable, U, in the interval 0.05;0.95.


### 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 (R2), Akaike’s information criteria [19], and the Bayesian information criterion [20]:
(32)AIC=-2·LLθ^+2p,BIC=-2·LLθ^+ln⁡(n)·p,
where LL(θ^) is the maximum log-likelihood function associated with model, θ^ is the maximum likelihood estimate of the parameters of the model, and p 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(33)T=-2·ln⁡LLθ^2LLθ^1,
where the test statistic T is asymptotically distributed as a chi-squared random variable χp2-p12 with degrees freedom p2-p1, where p1 and p2 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].


## 2.1. Stochastic Differential Equations Model

The focus of the present work is on the dynamics of the response variable (crown width),Y, as a stochastic process with respect to the predictor variable (diameter at breast height − diameter in the sequel), x. 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, y(x), are described using the ordinary differential equation
(1)dy(x)dx=αβγeβx-γyx,
where α, β, and γ are unknown fixed-effects parameters and verify d≥0>ln⁡γ/β, β>0, α>1. The formula describing the Bertalanffy trajectory follows the form of a sigmoidal function
(2)yx=y01-γe-βx1-γe-βx0α,x∈x0;X0,
where y0=y(x0). Equation (1) can be seen as a generalization of the Malthusian growth model with nonconstant fertility depending on the predictor variable (diameter) by r(x)=αβγ/(eβx-γ). 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
(3)rx=αβγeβx-γ+σεx,
where σ is the diffusion coefficient and ε(x) 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:
(4)α=α0+ϕi,
where parameters ϕi of the ith stand (i=1,2,…,M) are stand-specific random effects. It is assumed that the random effects ϕi(i=1,2,…,M) are independent and normally distributed with a mean of 0 and a constant variance (ϕi~N0;σϕ2). Therefore, the response variable, Yix, x≥0, i=1,2,…,M, is described using a stochastic differential equation of the Bertalanffy type:
(5)dYix=α0+ϕiβγeβx-γYixdx+σYixdWix,PYix0=y0=1,x∈x0;X0,
where Wix, x≥0, are the independent standard Brownian Motions, Wi(x) and ϕj are assumed to be mutually independent for all 1≤i, j≤M, and M 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 g(x) that introduces the dynamics of each exogenous stand variable into the model. The function g(x) is considered as a linear combination of the exogenous stand variables and is given by
(6)gx=α0+∑k=1qαkgkx,
where gk(x) for k=1,2,…,q are termed exogenous stand variables and constitute a diameter-continuous function in x∈x0;X0, α0,αk are diameter-independent real fixed-effects parameters (to be estimated), and q is the number of exogenous stand variables. One candidate set for the basis functions gk(x) 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:
(7)αx=α0+ϕi+∑k=1qαkgkx.
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:
(8)dYix=βγα0+ϕi+∑k=1qαkgkxeβx-γYixdx+σYixdWix,PYix0=y0=1,x∈x0;X0.
By Itô’s [16] lemma, the logarithmic transformation Z≡ln⁡(Y) transforms (8) to the following form:
(9)dZix=βγα0+ϕi+∑k=1qαkgk(x)eβx-γ-σ22dx+σdWix,PZix0=ln⁡(y0)=1.
Then, by integration and Y≡exp⁡(Z), we have
(10)Yix=y0exp⁡∫x0xβγα0+ϕi+∑k=1qαkgk(u)eβu-γdu-σ22x-x0+σ∫x0xdWiu.
Taking into account that the random variable ∫x0xdWiu has a normal distribution with mean 0 and variance x-x0 we can deduce that the process Yi(x) corresponds to a lognormal distribution, Λ1μix∣y0,x0;v(x∣x0). In this case, the conditional probability density function for the considered stochastic process Yi(x) is
(11)fiy,x=1y2πv(x∣x0)×exp⁡-12v(x∣x0)ln⁡y-μix∣y0,x02,
where
(12)μix∣y0,x0=ln⁡y0βγ1-γe-βx1-γe-βx0α0+ϕi+βγ∑k=1qαk∫x0xgk(u)eβu-γdu-σ22x-x0,(13)vx∣x0=σ2x-x0.
In this paper, we approximate the functions gk(d), k=1,2,…,q, by polygonal functions. We have real observations z0ki,zjki; j=1,2,…,ni (ni is the number of observations possessed in the ith stand) for exogenous stand variables, zk, in the values of predictor variables x0 and xjki. Hence, for k=1,2,…,q we have
(14)gkx=z0ki+zjki-z0kixjkix-x0,0≤x≤xjki,j=1,2,…,ni.
More often than not each exogenous stand variable in a particular stand has the same observed value, zjki=zki; j=1,2,…,ni. From (14), we deduce that analytical expression of the integral to the right-side equation (12) is
(15)∫x0xgk(u)eβu-γdu=1βγxGkx,
where
(16)Gkx=z0kx-x0zk-z0kln⁡1-γe-βx1-γe-βx0-zj-z02βx2-x02+xln⁡1-eβxγ-x0ln⁡1-eβx0γ-2Li2eβxγ-Li2eβx0γ,
where Li2x is the dilogarithm function defined by Li2x=∫0x(ln⁡1-u/u)du. The conditional mean, mode, median, quantile, variance, and coefficient of variation functions, respectively, m(·), mo(·), me(·), mq(·)(0<q<1), and w(·), of the stochastic crown width process are given by
(17)mix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡⁡αkGkx,(18)moix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡⁡αkGkx·exp⁡-32σ2x-x0,(19)meix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡⁡αkGkx·exp⁡-12σ2x-x0,(20)mqix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡(αkGk(x))·exp⁡-12σ2x-x0+σx-x0z1-q,(21)wix=y01-γe-βx1-γe-βx0α+ϕi·∏k=1qexp⁡(αkGk(x))2·exp⁡σ2x-x0-1,(22)cvx=exp⁡σ2x-x0-1.


## 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, ϕi, i=1,…,M, are assumed to be equal to its mean value Eϕi=0). Second, the maximum log-likelihood function is derived for generalized mixed-effects model. The fixed-effects parameters α0,α1,…,αq,β,γ,σ,σϕ and random-effect parameters ϕi, i=1,…,M, 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 (l1i,l2i,…,lnii) at the diameters (d1i,d2i,…,dnii), where ni is the number of observed trees of the ith plot, i=1,2,…,M. Under the initial condition PL0=1.0=1, the associate log-likelihood function for generalized fixed-effects parameters model defined by (8) can be obtained by the following expression:
(23)LL1α0,α1,…,αq,β,γ,σ=∑i=1M∑j=1niln⁡flji,dji,
where ni is the number of observed trees of the ith plot i=1,2,…,M, conditional density function fl,d takes the form defined by (11), and random effects ϕi≡0; i=1,2,…,M. The maximum log-likelihood function for generalized mixed-effects parameters model defined by (8) takes the following form:
(24)LL2(α0,α1,…,αq,β,γ,σ,σϕ)=∑i=1M∫R∑j=1niln⁡fi(lji,dji)+ln⁡p(ϕi∣σϕ)·dϕi,
where α0,α1,…,αq,β,γ,σ,σϕ are fixed-effects parameters (the same for all stands) and ϕi are random-effects parameters (stand-specific), which are assumed to follow a normal distribution with mean 0 and standard deviation σϕ, p(ϕi∣σϕ) 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 g:R→R by
(25)giϕi∣α0,α1,…,αq,β,γ,σ,σϕ=∑j=1niln⁡fi(lji,dji)+ln⁡p(ϕi∣σϕ),i=1,2,…,M.
The integral ∫R2egiϕi∣α0,α1,…,αq,β,γ,σ,σϕ·dϕi, i=1,2,…,M, by a second-order Taylor series expansion can be approximated as the Laplace approximation [18]:
(26)ln⁡∫R2egiϕi∣α0,α1,…,αq,β,γ,σ,σϕ·dϕi≈giϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ+12ln⁡2π-12ln⁡-Hϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ,
where
(27)ϕ^i=argmaxϕigiϕi∣α0,α1,…,αq,β,γ,σ,σϕ,(28)Hi(ϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ)=∂2giϕi∣α0,α1,…,αq,β,γ,σ,σϕ∂2ϕiϕi=ϕ^i.
The log-likelihood function for a generalized mixed-effects SDE model defined by (8) is approximately given by
(29)LL2(α0,α1,…,αq,β,γ,σ,σϕ)≈∑i=1Mgiϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ+12ln⁡2π-12ln⁡-Hiϕ^i∣α0,α1,…,αq,β,γ,σ,σϕ.
The maximization of LL2(α0,α1,…,αq,β,γ,σ,σϕ) is a nested optimization problem. The internal optimization step estimates the ϕ^i for every stand i=1,2,…,M by (27). The external optimization step maximizes LL2(α0,α1,…,αq,β,γ,σ,σϕ) after plugging the ϕ^i 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,d, and crown width, l, 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 r trees with crown width li and diameter di, i=1,2,…,r, is taken from a new stand, the random-effects parameter ϕ for the new stand can be predicted in the following form:
(30)ϕ^=1r∑i=1rln⁡1-γ^exp⁡-β^di1-γ^exp⁡-β^d0-1ln⁡li-ln⁡l0×ln⁡1-γ^exp⁡-β^di1-γ^exp⁡-β^d0-1-α^0,
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 L(d). The stochastic predictions of ϕ and L(d) can be defined in the following form, respectively:
(31)ϕstoch=ϕ^+ΦU-10;σ^ϕ2,lstoch=LNU-1μi^(d∣l0,d0);v^(d∣d0),
where ϕ^ is the estimation value of random effects obtained by (30), σ^ϕ is the estimation value of the standard deviation of random effects, ΦU-10;σ^ϕ2 is the inverse of the normal distribution with a mean of 0 and a constant variance σ^ϕ2 for a uniform random variable, U, in the interval 0.05;0.95, μi^(d∣l0,d0) and v^(d∣d0) are the estimated trend of the mean and variance (calculated by (12) and (13)) of the lognormal density of the crown width, respectively, and LNU-1μi^(d∣l0,d0);v^(d∣d0) is the inverse of the lognormal distribution with a mean of μi^(d∣l0,d0) and a variance of v^(d∣d0) for a uniform random variable, U, in the interval 0.05;0.95.


## 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 (R2), Akaike’s information criteria [19], and the Bayesian information criterion [20]:
(32)AIC=-2·LLθ^+2p,BIC=-2·LLθ^+ln⁡(n)·p,
where LL(θ^) is the maximum log-likelihood function associated with model, θ^ is the maximum likelihood estimate of the parameters of the model, and p 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(33)T=-2·ln⁡LLθ^2LLθ^1,
where the test statistic T is asymptotically distributed as a chi-squared random variable χp2-p12 with degrees freedom p2-p1, where p1 and p2 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].


## 3. Results and Discussion

### 3.1. Data

We focus on the modeling of Scots pine (Pinus SylvestrisL.) tree dataset. Scots pine tree stands dominate Lithuanian forests, grow onArenosolsandPodzols forest sites,and cover 725,500 ha. At plot establishment, the following data were recorded for every sample tree: crown width, l, diameter over bark at 1.30 m, d, and crown height, hc. 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), d, crown height, hc, and crown width, l, of all the trees used for model estimation and validation are presented in Table 1.Table 1
Characterization of datasets.
Data
Number of trees
Min
Max
Mean
St. Dev.
Number of trees
Min
Max
Mean
St. Dev.
Estimation
Validation
l (m)
1630
0.65
7.20
3.26
1.22
699
0.55
7.30
3.26
1.24
d (cm)
1630
1.7
51.1
24.29
8.33
699
6.4
49.3
24.04
8.63
h
c (m)
1630
1.8
27.8
16.16
4.47
699
0.9
27.6
15.38
3.68


### 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 theNLPSolve procedure in a mathematical software packageMAPLE [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 2
Estimated parameters (standard errors) of all the used fixed-effects and mixed-effects models applied to the estimation dataset.
Models(number of parameters)
Parameters
α
0
β
γ
α
1
σ
σ
ϕ
Fixed(4)
2.4943(0.0153)
2.83* 10−6 (2.69 * 10−7)
0.9999(1.50 * 10−5)
—
0.0546(0.0010)
—
Fixed +hc (5)
0.9532(0.1141)
−0.0159(0.0029)
1.1031(0.0074)
−0.0238(0.0018)
0.0495(0.0074)
—
Mixed(5)
1.9439(0.0964)
−0.0177(0.0015)
2.8998(0.3025)
—
0.0742(0.0008)
0.3394(0.0671)
Mixed +hc (6)
1.7490(0.00670)
7.88* 10−5 (3.04 * 10−5)
1.0010(0.0009)
0.0068(0.0009)
0.0452(0.0008)
0.0638(0.0134)Table3 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.Table 3
Estimation goodness of fit statistics for all the used fixed-effects models applied to the estimation and validation datasets*.
Models(number of parameters)
Estimation
Validation
B
*
*, m (PB, %)
RMSE***, m(PRMSE, %)
AIC(BIC)
R
2
*
*
*
*
B, m (PB, %)
RMSE, m (PRMSE, %)
AIC
R
2
Fixed(4)
−0.0504(−6.70)
0.7104(21.80)
10922.1(10943.79)
0.6584
−0.0236 (−5.18)
0.6686(20.46)
4012.63(4030.82)
0.7066
Fixed +hc (5)
−0.0477(−4.74)
0.6516(20.00)
10634.64(10661.61)
0.7123
0.1117(−6.21)
0.6909(21.14)
4059.46(4082.20)
0.6867
Mixed(5)
0.0350 (−3.41)
0.6414(19.67)
10590.63(10617.60)
0.7215
0.0344(−2.01)
0.6390 (19.55)
3950.60 (3973.74)
0.7319
Mixed +hc (6)
0.0374(−2.68)
0.6238 (19.13)
10501.09 (10533.46)
0.7366
0.1105(−0.70)
0.6716(19.98)
4020.96(4048.25)
0.7270
*The best values of the performance statistics for all scenarios of stand variables are in bold, **the mean prediction bias B=(1/n)∑i=1n(yi-yi∧) and the percentage mean prediction bias PB=(1/n)∑i=1n((yi-yi∧)/yi),  ***the root mean squared error RMSE=(1/(n-p))∑i=1n(yi-yi∧)2 and the percentage root mean squared error PRMSE=(1/(n-p))∑i=1n((yi-yi∧)/yi)2, ****an adjusted coefficient of determination R2=(1-((n-1)/(n-p)))(∑i=1n(yi-yi∧)2/∑i=1n(yi-y-)2), and p is the number of parameters in the model, n=∑i=1Mni is the total number of observations used to fit the model, M is the number of stands, ni is the number of trees in ith stand, and yi, yi∧, and y¯ are the measured, estimated, and average values of the dependent variable (crown width).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 Table3). The statistical significance of the difference between two models can be assessed by the simple test defined by (33). For the estimation dataset the P 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 Figure1) 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.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).
(a)
(b)In order to illustrate the behavior of the new developed mixed-effects SDE Bertalanffy type models, crown widthversus 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 (l=1.0; d=0.0).Figure 2
Predicted mean (equation (17)) and 5% and 95% quantiles (equation (20)) functions of the crown width: mean—solid line, quantiles—dash line, and observations—circle.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.


## 3.1. Data

We focus on the modeling of Scots pine (Pinus SylvestrisL.) tree dataset. Scots pine tree stands dominate Lithuanian forests, grow onArenosolsandPodzols forest sites,and cover 725,500 ha. At plot establishment, the following data were recorded for every sample tree: crown width, l, diameter over bark at 1.30 m, d, and crown height, hc. 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), d, crown height, hc, and crown width, l, of all the trees used for model estimation and validation are presented in Table 1.Table 1
Characterization of datasets.
Data
Number of trees
Min
Max
Mean
St. Dev.
Number of trees
Min
Max
Mean
St. Dev.
Estimation
Validation
l (m)
1630
0.65
7.20
3.26
1.22
699
0.55
7.30
3.26
1.24
d (cm)
1630
1.7
51.1
24.29
8.33
699
6.4
49.3
24.04
8.63
h
c (m)
1630
1.8
27.8
16.16
4.47
699
0.9
27.6
15.38
3.68


## 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 theNLPSolve procedure in a mathematical software packageMAPLE [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 2
Estimated parameters (standard errors) of all the used fixed-effects and mixed-effects models applied to the estimation dataset.
Models(number of parameters)
Parameters
α
0
β
γ
α
1
σ
σ
ϕ
Fixed(4)
2.4943(0.0153)
2.83* 10−6 (2.69 * 10−7)
0.9999(1.50 * 10−5)
—
0.0546(0.0010)
—
Fixed +hc (5)
0.9532(0.1141)
−0.0159(0.0029)
1.1031(0.0074)
−0.0238(0.0018)
0.0495(0.0074)
—
Mixed(5)
1.9439(0.0964)
−0.0177(0.0015)
2.8998(0.3025)
—
0.0742(0.0008)
0.3394(0.0671)
Mixed +hc (6)
1.7490(0.00670)
7.88* 10−5 (3.04 * 10−5)
1.0010(0.0009)
0.0068(0.0009)
0.0452(0.0008)
0.0638(0.0134)Table3 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.Table 3
Estimation goodness of fit statistics for all the used fixed-effects models applied to the estimation and validation datasets*.
Models(number of parameters)
Estimation
Validation
B
*
*, m (PB, %)
RMSE***, m(PRMSE, %)
AIC(BIC)
R
2
*
*
*
*
B, m (PB, %)
RMSE, m (PRMSE, %)
AIC
R
2
Fixed(4)
−0.0504(−6.70)
0.7104(21.80)
10922.1(10943.79)
0.6584
−0.0236 (−5.18)
0.6686(20.46)
4012.63(4030.82)
0.7066
Fixed +hc (5)
−0.0477(−4.74)
0.6516(20.00)
10634.64(10661.61)
0.7123
0.1117(−6.21)
0.6909(21.14)
4059.46(4082.20)
0.6867
Mixed(5)
0.0350 (−3.41)
0.6414(19.67)
10590.63(10617.60)
0.7215
0.0344(−2.01)
0.6390 (19.55)
3950.60 (3973.74)
0.7319
Mixed +hc (6)
0.0374(−2.68)
0.6238 (19.13)
10501.09 (10533.46)
0.7366
0.1105(−0.70)
0.6716(19.98)
4020.96(4048.25)
0.7270
*The best values of the performance statistics for all scenarios of stand variables are in bold, **the mean prediction bias B=(1/n)∑i=1n(yi-yi∧) and the percentage mean prediction bias PB=(1/n)∑i=1n((yi-yi∧)/yi),  ***the root mean squared error RMSE=(1/(n-p))∑i=1n(yi-yi∧)2 and the percentage root mean squared error PRMSE=(1/(n-p))∑i=1n((yi-yi∧)/yi)2, ****an adjusted coefficient of determination R2=(1-((n-1)/(n-p)))(∑i=1n(yi-yi∧)2/∑i=1n(yi-y-)2), and p is the number of parameters in the model, n=∑i=1Mni is the total number of observations used to fit the model, M is the number of stands, ni is the number of trees in ith stand, and yi, yi∧, and y¯ are the measured, estimated, and average values of the dependent variable (crown width).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 Table3). The statistical significance of the difference between two models can be assessed by the simple test defined by (33). For the estimation dataset the P 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 Figure1) 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.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).
(a)
(b)In order to illustrate the behavior of the new developed mixed-effects SDE Bertalanffy type models, crown widthversus 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 (l=1.0; d=0.0).Figure 2
Predicted mean (equation (17)) and 5% and 95% quantiles (equation (20)) functions of the crown width: mean—solid line, quantiles—dash line, and observations—circle.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.



---

*Source: 375270-2015-02-23.xml*