New Insights into Tree Height Distribution Based on Mixed Effects Univariate Diffusion Processes

The aim of this paper is twofold: to introduce the mathematics of stochastic differential equations (SDEs) for forest dynamics modeling and to describe how such a model can be applied to aid our understanding of tree height distribution corresponding to a given diameter using the large dataset provided by the Lithuanian National Forest Inventory (LNFI). Tree height-diameter dynamics was examined with Ornstein-Uhlenbeck family mixed effects SDEs. Dynamics of a tree height, volume and their coefficients of variation, quantile regression curves of the tree height, and height-diameter ratio were demonstrated using newly developed tree height distributions for a given diameter. The parameters were estimated by considering a discrete sample of the diameter and height and by using an approximated maximum likelihood procedure. All models were evaluated using a validation dataset. The dataset provided by the LNFI (2006–2010) of Scots pine trees is used in this study to estimate parameters and validate our modeling technique. The verification indicated that the newly developed models are able to accurately capture the behavior of tree height distribution corresponding to a given diameter. All of the results were implemented in a MAPLE symbolic algebra system.


Introduction
Understanding the key forces that shape tree heights distribution patterns and their dynamics through average breast height diameter within a forest stand (in the sequel-diameter) is a fundamental goal of forestry [1]. Stand volume, one of the most important variables in forest management, is heavily dependent on tree diameter and height distribution. The literature on forestry reports that tree height distribution varies across different stands and/or species. The tree height distribution is of prime importance from the point of view of the quality and quantity of a stand and its future growth. The importance of using the tree height rather than the tree diameter as a predictor of forest demographics arises from the former's high potential for predicting the properties of forest productivity as pointed out by Kempes et al. [2]. Traditional methods quantify the tree size distribution in an even-age forest stand [3]. Unfortunately, height and diameter distributions cannot be combined if they are estimated independently a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 using datasets from different stands. Much research has been conducted on the utilization of various theoretical functions for height and diameter distribution modeling techniques for improving stand volume prediction, such as Johnson's [4][5], beta distribution [6], and powernormal [3]. Recently, there were also a few results published on the use of the copula approach for the modeling of tree height and diameter distribution in stands [7], [8].
The relationship between height and diameter varies for the same tree in different forest stands, such that there is a distribution of tree heights for any given tree diameter based on environmental conditions, or a random site effect. The different height-diameter relationships affect growth predictions and stand trajectories. The new developed stochastic differential equation (SDEs) based modeling approach for complex stands uses stochastic height-diameter relationships at the individual-tree level representing tree growth and neighborhood interactions that are then aggregated to predict the stand height structure. In this study, to project height distribution for a given diameter, a one-dimensional SDE with mixed effects was employed. The main feature of mixed effects models is that they allow parameter vectors to vary from plot to plot by splitting regression coefficients into a fixed part, common to the population, and random components, specific to each plot [9]. Mixed effects models allow fixed and random parameters to be estimated simultaneously and evaluate the value of the random parameters for a location not present in the original estimation dataset. This approach is known as calibration and can be applied if a sub-sample of trees measured for the total height and breast height diameter are available [10]. Fixed effects parameter SDEs are used in a wide range of applications in environmental, engineering, and biological modeling [11][12][13][14]. Discrete stationary stochastic models defined by Markov chains have been used to describe size-structure predictions [15].
The essential features of developed height distributions for a given diameter may be explained as follows. Heights are measured at different diameters in a number of sample plots. The diameters and number of measurements differ among plots and the measurements of the diameters are not evenly spaced. The diameter based dependent tree height distribution models are assumed to have some fixed effects parameters that are common to all plots and random effects that are specific to each plot. Two sources of variation were simultaneously included for modeling tree height distribution: variability between plots using a random effects approach and variability in the individual tree height using system noise, which reflects the random fluctuations around the corresponding theoretical height-diameter model. New developed conditional probability density functions of a tree height at a given diameter based on diffusion processes can be used for calculating the mean value of growth and yield attributes and its coefficient of variation as a function of tree height at any specified diameter. The random effects SDEs height-diameter relationships allow taking into account the effect of multiple causal relations in the model, the influence of unknown covariates affecting the height growth and they allow for developing height distribution accounting for spatial variability in large-scale modeling.
In this study, the evolution of a random variable (height), H(d), for a given diameter, d, is modeled using mixed effects SDEs from the Ornstein-Uhlenbeck family [16], for example, the Vasicek, the Gomperz (3-parameters and 4-parameters), the Bertalanffy, and the Gamma. We focused on mixed effects SDEs with a deterministic term depending on random effects and a stochastic term without random effects.
The aim of this study is to present the advantages of SDEs with mixed effects in analyzing tree height distributions for a given diameter and their application to describe the evolution of the height-diameter ratio, quantile curves, mean tree height, mean stem volume and the coefficients of variation for the mean tree height and stem volume. We also discuss how a conditional height's probability density function for a given diameter can be used to construct maximum likelihood estimators using large collections of datasets provided by the Lithuanian National Forest Inventory (LNFI). A MAPLE program was used to carry out the calculations.

Materials and Methods
A typical diffusion process is modeled as a differential equation involving a deterministic (drift) term and a stochastic (diffusion) term, the latter represented by Brownian motion [14]. Traditionally used ordinary differential equation models are the Malthus, Mitscherlich, Gompertz, and Bertalanffy types [11][12][13] (see Appendix A).
There are alternative ways of introducing stochasticity in an ordinary differential equation. In this work, the tree height randomness was approximated as a standard Brownian motion [11][12][13][14]. Therefore, the complete deterministic models defined by Eqs A.1, A.3, A.4, A.7 and A.9 were converted, into stochastic models assuming that the deterministic parameter, α, varies randomly around the mean: aðdÞ ¼ ( a þ sεðdÞ; for Eqs: C:1; C:3; C:4; abg e bd À g þ sεðdÞ; for Eq: C:7; a d À b þ sεðdÞ; for Eq: C:9: where σ (σ>0) is the diffusion coefficient, which reflects random fluctuations around the corresponding theoretical height-diameter curve, and ε(d) is a Gaussian white noise process. If the magnitude of the parameter capturing system noise, σ, is zero, the entire system noise term will vanish, and the remaining part of the SDEs will simply be differential forms, the solutions to which are Eqs A.2, A.5, A.6, A.8 and A.10, respectively. The relationships between total tree height and diameter are altered by environmental conditions. Among other plot-specific characteristics such as soil type, nutrient status, resistance of trees to windthrow, competition for light, and elevation cause the parameters to differ across plots. In the case of between-plot variations, the fixed effects parameters α, β, and σ vary from plot to plot and, hence, account for these variations. For the construction of the mixed effects parameters models, the first step is to determine which parameters should be considered mixed effects and which should be considered purely fixed effects. The parameters with high variability could be considered mixed effects. The parameter α has high variation between plots for all used SDEs models, so it can be altered by adding plot-specific random effects to the fixed effects parameter to produce a plot-specific parameter in the following form: where ϕ i (i = 1, 2, . . ., M)-plot-specific random effects, M is the number of plots. It is assumed that the random effects, ϕ i , i = 1, 2, . . ., M, are independent and normally distributed with 0 mean and constant variance s 2 0 (0 i $ Nð0; s 2 0 Þ). In order to derive mixed effects SDEs height-diameter models, it is sufficient to substitute Eqs 1 and 2 into Eqs A.1, A.3, A.4, A.7 and A.9. In this study, the tree height, H i (d), i = 1,2,. . ., M, evolving in M different experimental plots randomly chosen from a theoretical population was described by the Itô [17] sense SDE of the Vasicek type: the 3-parameters Gompertz type:  An approximated maximum likelihood procedure (see Appendix C) was used for the estimation of the fixed effects parameters and random effects by assuming that tree height and diameter observations are without measurement noise.

Data
The data used for developing the models were obtained from the Lithuanian National Forest Inventory (LNFI) (2006)(2007)(2008)(2009)(2010). The NFI plots are systematically distributed using a grid of 4x4 km squares with a random starting point. The sample plots are arranged into triangle distributed clusters with a distance between angles of 2 km. Each cluster has 4 sample plots. They are situated on each 250 m length side of square 25 m from its angles [19]. At plot establishment, the following data were recorded for every sample tree: the species, the diameter over bark at 1.30 m high and measured to the nearest millimeter and the total height to the nearest quarter meter. The tree diameters were measured with outside calipers in two perpendicular directions. A total of 3,455 plots (500 m 2 ) of Scots pine trees were chosen from the LNFI 2006-2010 database. A random sample of 1,999 plots (7,343 trees) was selected for model estimation, and the remaining dataset of 1,456 plots (5,413 trees) was utilized for model validation. Only measurements from live trees without top damage were included in the statistical analysis. Summary statistics for the diameter at breast height (d), the total height (h) and the age (A) for all of the trees used in model estimation and validation datasets are presented in Table 1. Table 2 presents the distribution of the number of trees per plot measurements from both datasets. It should be noted that data on the number of plots with greater than 10 measured trees were very limited.

Results
To examine the impact of fixed and random effects on the prediction of the height distribution, the maximum likelihood estimators (Eqs A.2 and A.10) were calculated using the NLPSolve procedure in MAPLE 11 [20]. The models with fixed effects and mixed effects were evaluated based on Akaike's [21] information criterion (AIC), which was defined as follows: where LL s K is the log-likelihood function and p is the number of parameters in the model. The nested models with the smallest AIC value are considered to be the best. Using the estimation dataset presented in Table 1, the parameter estimates of the fixed effects and mixed effects SDEs height-diameter models, defined by Eqs 3-7, are summarized in Table 3. The standard errors of the parameter estimates were calculated by Eq C12. All of the parameters of the fixed effects and mixed effects SDEs height-diameter models are highly significant (p < 0.001). The AIC values for the fixed effects SDEs height-diameter models were more than for mixed effects models, indicating that random effects are needed in the height-diameter SDEs.

Height distributions
Tree height structure is a basic modeling component of many complex forest yield models relating individual tree characteristics with stand variables. The distribution of the tree height, as a diameter dependent variable, can be approximated by classifying diameter and applying the desired transformation to the mean tree of the class [22]. A more convenient way to derive tree height distributions for a given diameter is the use of SDEs. This paper described research aimed at deriving tree height probability densities for a given diameter (Eqs B.1, B.4, B.7, B.10 and B.13) by directly fitting the SDEs (Eqs 3-7) to the diameter and height observations. Several empirical methods are available for comparing conditional probability densities, as has been illustrated by [23]. In the present paper, a well-known measure of distributional accuracy named by Kullback-Leibler Information Criterion (KLIC) [24] was utilized. We are interested in comparing two conditional probability density functions f A h; dj b y s and Under appropriate conditions, the KLIC has limiting distribution under the null, and is consistent against all possible fixed alternatives. The expression for KLIC in Eq 9 depends on the unknown expectation E(Á). We consider estimating KLIC by a discrete height sample  Left-fixed effects models; right-mixed effects models; first plot-solid line and height's dataset-cross; second plot-dash line and height's datasetdiamond; third plot-dot line and height's dataset-box; mean diameter within a plot is recorded in the graph.
where n i is the number of observed trees of the ith plot, i = 1,2,. . .,M.
Analysis of paired comparison of the five conditional probability densities, described in Section 2 by Eqs B.1, B.4, B.7, B.10 and B.13 was performed by KLIC calculated using the estimation dataset. The results of comparisons are presented in Table 4. As we see in Table 4, the Vasicek type conditional probability density function of the tree height with mixed effects (see Table 4, values above diagonal) and fixed effects (see Table 4, values below diagonal) are superior to the other densities and the worst conditional probability density function is the Gompertz (3-parameters) type with mixed and fixed effects. All mixed effects density functions are superior to corresponding fixed effects parameters density functions (see bold values in diagonal Table 4).

Height-diameter models
Many comparisons between the different models or ecoregions have been carried out to identify the appropriate height-diameter relationships within stands. The height dynamics defined by Eqs B.16, B.18, B.20, B.22 and B.24 are affected by many processes and vary among stands. Fig 2 illustrates the influence of the plot within Lithuanian pine forests on the mean and standard deviation of height-diameter dynamics using the Vasicek and 4-parameters Gompertz diffusion processes and random-effects parameter, ϕ, for the 3 randomly selected plots from the estimation dataset. The parameter estimates for each plot are calculated by adding the fixed effect parameter and random effect. Therefore, considering the asymptotic maximum height parameter, α+ϕ, for the Vasicek type model, the values varied from plot to plot. , fixed effects models, mixed effects models and mixed effects models with the random effects set to zero scenarios were used to predict tree height in both the estimation and validation datasets. The performance statistics of new developed tree height's equations included three statistical indices: prediction accuracy, δ, which combines the mean bias, B, and the variation, ξ, of the biases, enabling improved assessment of model accuracy; an adjusted coefficient of determination, " R 2 , which reflects the part of the total variance explained by the equation; and Akaike's information criterion, AICC, which measure the quality of the heightdiameter equation for a given dataset. The expressions for these statistics are as follows: where n is the total number of observations used to estimate the height-diameter model, p is the number of model parameters, and y i , b y i , and " y are the measured, predicted, and average values of the dependent variable (total tree height), respectively.  Table 5 presents the performance statistics for the tree height's models for all three scenarios; these include the fixed effects model, the mixed effects model and the mixed effects model with random effects set to zero, illustrating the extent to which the inclusion of the random effects improved the performance statistics for both the estimation and validation datasets. The random effects for the validation dataset were calibrated using Eqs D.1-D.5, respectively. The results of this study show (see Table 5, Akaike's information criterion) that the SDEs Vasicek and Gompertz 4-parameters type tree height's models with mixed effects are significantly superior at predicting tree height compared to the other newly developed models. Compared to the basic fixed effects models, the mixed effects models show better performance with lower bias and prediction accuracy, and with a higher adjusted coefficient of determination evaluated over the entire dataset. The mixed effects models with random effects set to zero show the worst performance, with greater bias and prediction accuracy, and with a lower adjusted coefficient of determination evaluated over the entire dataset. The fixed effects models, the mixed effects models, and the mixed effects models with random effects set to zero have very similar fit statistics for both the estimation and validation datasets.
The plots of the residuals versus predicted heights and the lowess line [25], estimated for the validation dataset, in the fixed effects and mixed effects scenarios (random effects for the validation dataset were calibrated by Eqs D.1-D.5) are presented in Fig 3. Fig 3 shows that the residuals that were calculated using the mixed effects scenario are distributed more symmetrically around zero, with approximately constant variance, compared with the fixed effects scenario. A non-parametric smoothing line, called a lowess line, shows a clear trend in the middle range of predicted height; however, what happens at the extremes is dictated by relatively little data.

Quantile regression
The conditional tree height's mean, defined by Eqs B. 16 and for the Gompertz 4-parameters type model:

Slenderness ratio
Tree height to diameter ratio (slenderness ratio) is regarded as an index of the resistance of trees to windthrow and competition for light, and its mean value may be useful in determining stand stability. The slenderness ratio is calculated by dividing the tree height to its diameter at breast height. For the fixed effects and mixed effects SDEs height-diameter models, defined by Eqs 3-7, the slenderness functions are defined as follows: The slenderness ratio dynamics for three randomly selected plots are presented in Fig 5. The findings of our investigation generally support that the height's conditional probability densities driven by diameter correctly predict slenderness ratio (see Fig 5). Furthermore, all new developed tree height's distribution models show a decrease of slenderness ratio with increasing diameter.

Mean stem volume
The fixed effects and mixed effects height's conditional probability density functions allow us to revise mean stem volume calculation in the following form: Here V(d,h) is the stem volume regression function of power form [27], where parameters β 1 , β 2 , and β 3 are to be estimated. The selection of stem volume model was  The direct effects of stand variables such as site index and management practices and thinning could be included in the new developed models (see [28], [29]); however, their indirect effect via mixed effects (see right side Fig 6) has been included in mixed effects tree mean volume models.

Coefficient of variation for height and volume
The coefficient of variation is typically used to indicate the precision of the dispersion of datasets and is also often used to compare numerical distributions measured at different scales. Tree height based and tree volume based quantifications of the stand structural diversity can Tree Height Distribution Based on Univariate Diffusion Processes be performed using the coefficient of variation. The coefficient of variation reaches its maximum with two-storied stands. The coefficient of variation of tree height (tree volume) measures the variability of tree height (tree volume) relative to its mean and relates the mean and standard deviation by expressing the standard deviation as a percentage of the mean. To further discuss the results of this study, the coefficient of variation, which may help examine dispersion in tree heights occurring at diameter d, is defined by:

Discussion
The models commonly used of height distribution fitting in a forest stand are supplemented by tree height's measurements. However, in the Lithuanian National Forest Inventory foresters measure no more than 15 heights (see Table 2) of pine trees per stand. For estimating the parameters by traditionally used maximum likelihood technique such sample sizes are too small [5], [30]. New developed height distribution models based on mixed effects parameters diffusion processes overcome such weakness. The pioneer of the SDE approach in forest growth modeling is Suzuki [31]. In this paper for height-diameter evolution were used linear and non-linear SDEs from the Ornstein-Uhlenbeck family by incorporating random effects into deterministic (drift) term. This extended model describes the within-stand variation in data through the system noise reflecting the random fluctuations around the corresponding theoretical height-diameter curve and the between-stand variation in data through the random effects. The maximum likelihood estimation procedure converged for all five diffusion processes using the estimated dataset from LNFI.
In order to predict the parameters of the tree size probability density function for a new stand, traditionally were carried out regression models from the different stand variables [32]. If the diameter and height of a sub-sample of trees are known, then for new developed height distributions based on univariate diffusion processes the random effects can be calibrated by Eqs D.1-D.5.
Quantifying variability in tree height at a given diameter by a distribution law has both theoretical and practical value. First, knowledge of the tree height distribution in forest stands is important for understanding of competition and self-thinning which must studied not in the mean size of trees but in the size structure of trees in a forest stand. Second, understanding tree height distribution at a given diameter is important for improving estimates of stand biomass and carbon storage. To describe how tree height distribution vary across regional scales, we developed new empirical distributions of tree height at a given diameter across the Scots pine trees in Lithuania. In this paper our specific objectives were to test (1) what new developed probability density function based on stochastic differential equation height-diameter evolution provides the best fit, (2) how new developed models explain observed variation in the probability density functions, the mean height-diameters, quantile height-diameter, mean slenderness ratio and mean stem volume relationships across the Scots pine trees in Lithuania, and (3) how to describe the mean height-diameter and mean stem volume relationships fit in terms of the relative sizes defined by coefficient of variation.
The Kullback-Leibler Information Criterion [24] was used to compare all new developed conditional probability density functions using the estimation dataset. The conditional probability density function derived from the Vasicek type height-diameter univariate diffusion process showed better results than the other used stochastic processes (Table 4). All mixed effects parameters probability density functions are superior to corresponding fixed effects parameters density functions (see bold values in diagonal, Table 4). Theoretical validating that a height dataset observed at discrete diameters follows univariate probability density functions defined by Eqs B.1, B.4, B.7, B.10 and B.13 is not easy and there is no simple statistical test. The goodness of fit of the estimated univariate density functions (Eqs B.1, B.4, B.7, B.10 and B.13) graphically were illustrated in Fig 1 using fixed effects and mixed effects parameters scenarios and three randomly selected plots from an estimation dataset by plotting the estimated probability density functions and height's measurements. Fig 1 showed that the mixed effects and fixed effects parameters estimated probability density functions well capture the main features of the data from three randomly selected plots. The height-diameter evolution can be written using a wide range of mathematical relationships from linearized fixed effects regression equations to nonlinear mixed effects generalized relationships. Mathematical technique of a system of uniform diameter and height regional functions is the approach known as the generalized model. The mixed effects regression models are able to achieve the same results than the generalized model [10,33]. In this study new developed mixed effects parameters height-diameter relationships demonstrated similar statistical indexes as in the nonlinear generalized heightdiameter regression models presented by Petrauskas et al. [34].
In addition, one of the advantages of using diffusion processes for quantifying tree height distribution is that it allows to derive the first two moments about height's and volume's evolutions through diameter and to calculate the relative standard deviation (coefficient of variation) for the height and volume. Fig 8 shows the variation of the coefficient of variation in pine trees forest stands from the estimation dataset from LNFI as a function of mean plot diameter using the mixed effects Vasicek type diffusion process. There is an exponential increase of the coefficient of variation as the mean diameter per plot decreases; the coefficients of variation for tree height varies from 6.94% to 24.72% and for stem volume varies from 4.77% to 17.05%.

Conclusions
This study demonstrated the use of SDEs to quantify tree height distribution at a given diameter in a forest stand using the Lithuanian National Forest Inventory dataset. The results indicated that it is possible to measure mean tree height and volume evolution with an acceptable accuracy over a broad area of Lithuania. Overall, the models explained over 90% of the variation in height predictions observed in the LNFI (2006-2010) dataset. The remaining variation was likely to have resulted from stand variables. Better performance can be expected by introducing stand variables [29]. The diffusion processes based SDE models described here implicitly model spatial effects. The technique we described can be used for developing a new generation of forest growth models.
A system of bivariate stochastic differential equations with mixed-effects parameters could be used to develop tree diameter and height at a given age (or trivariate: diameter, height and stand density at a given age) distribution model. This extension to multivariate SDEs come with an increased computational burden.
Results for both tree height and volume predictions using the mixed effects SDE Vasicek type height-diameter model indicate that the coefficient of variation over all plots for the tree height and volume (at the mean diameter of a plot) takes values from the interval 6.9%-24.8% and 1.7%-16.0%, respectively, and evolves to a stationary value from the interval 6.6%-19.8% and 1.7%-13.0%, respectively.
The field of SDEs is a large and growing area of applied mathematics that is being increasingly used to model biological systems. In this paper, new mixed effects height's probability density functions for a given diameter were developed using an Ornstein-Uhlenbeck SDE family. Unfortunately, measurements from at least one tree in a stand, or their measure of central tendency (mean, median, mode of diameter and height) are required for the practical calibration of the random effects for a new stand. The use of the mixed effects model enables us to develop a simple model structure without including additional predictor stand variables.
The results showed that the mixed effects Vasicek type tree heights distribution models are superior to the other new developed models.
The variance functions developed here can be applied to generate weights in every linear and nonlinear least squares regression height-diameter model by the weighted least squares form.

Appendix A Deterministic models
The mathematical representation of Mitscherlich growth [35] is derived from physical chemistry, where it describes a first order irreversible chemical reaction. The deterministic heightdiameter model used to describe the individual growth of a tree in terms of its size (height), h (d), at instant (diameter), d, can be written in the form of an autonomous differential equation given by the following: where D 0 is the upper limit on the diameter at the breast height. Height dynamics are irreversible, and the growth rate is proportional to the difference between the asymptotic maximum height, α, and the already formed tree height, h(d), β is the proportionality constant (β>0). The formula describing a Mitscherlich type height-diameter trajectory takes the form: hðdÞ ¼ a þ ð1:3 À aÞexpðÀ bdÞ; d 2 ½0; D 0 : ðA:2Þ The changes in tree height, h(d), using deterministic ordinary differential equations, developed by Gompertz [36], for 2-parameters and 3-parameters models, respectively, are described The formulas describing a Gompertz type height-diameter trajectory for 2-parameters and 3-parameters models, respectively, are as follows: where α is the intrinsic growth rate of the height, β is the growth deceleration factor, γ is a threshold parameter, and exp a b represents the largest height size that the tree can sustain.
Von Bertalanffy (for a review, see example in Román-Román et al. [13]) 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, the tree height, h(d), are described using an ordinary differential equation: hðdÞ; d 2 ½0; D 0 : ðA:7Þ where α, β, and γ are unknown fixed effects parameters. The formula describing the Bertalanffy trajectory follows the form of a sigmoidal function: The changes in tree height, h(d), using the well-known regulated Malthusian growth process [37], are described in the following form: The conditional mean, m(d), and variance, v(d), functions of the tree height, H(d), for all the models (Eqs 14-17) are given by the following expressions [28], [29], [38]: specific), which are assumed to follow a normal distribution with 0 mean and constant variance s 2 0 , and p(ϕ i |σ ϕ ) is the normal density of the random effects. The integral in Eq C.4 does not have a closed form solution. Because analytic expression for the integrand in Eq C.4 is known, the Laplace method [39], [40] may be used. Let us define a function g:R!R as follows: where b 0 i is the global max of g K ð0 i jy 2 K Þ and the root of: . . . ; M; K 2 fV; G3; G4; B; Gg: ðC:9Þ The log-likelihood function for the mixed-effects SDEs height-diameter models is approximately given by: The maximization of LL 2 K ðy 2 K Þ is a two-step optimization problem. The internal optimization step estimates the b 0 i for every plot i = 1,2,. . .,M with Eq C.9. The external optimization step maximizes LL 2 K ðy 2 K Þ after plugging the b 0 i into Eq C.10.
To assess the asymptotic standard errors of the maximum likelihood estimators for the stochastic height-diameter models, a study of the Fisher [41] information matrix was performed. The approximate asymptotic variance of the approximated maximum likelihood estimators (Eq C. 10 The approximate asymptotic standard errors of the fixed effects parameters are defined by the diagonal elements of the matrix½I e ð b y s K Þ À 1 , s = 1,2, K 2 {V,G3,G4,B,G} by: ðC:12Þ

Appendix D Calibration and stochastic prediction
In the literature on forestry, calibration means that random effects are predicted using a supplementary sample of observations taken from a sampling unit. The tree heights for new stand can be predicted either by using random effects set to zero, or by adding random effects that were predicted from prior observations. When the diameter and height of a sub-sample of trees are known, the predicted random effects are added to the fixed effects parameters to obtain localized parameters for this sub-sample plot. Let us assume that a sub-sample of m trees with height h j and diameter d j , j = 1, 2, . . ., m, is taken from a new plot. Using height-diameter models defined by Eqs 3-7, the random effect, ϕ, for a new stand can be approximately calibrated as follows: ln h j À Á À ln 1:3 ð Þ ln h j À Á À ln 1:3 ð Þ þ b b d j À 0:001 À Á ln d j À Á À ln 0:001 ð Þ À b a 0 @ 1 A ; ðD:5Þ where a b ; b b ; b gare the parameter estimates calculated using the approximated maximum likelihood procedure (Eq C.10). The height of another tree from the same plot can be predicted by adding the random effect calibrated by Eqs D.1-D.5 to the fixed effects parameter a b , respectively. The random effects height distribution models explain much more variability than the fixed effects models and provide better height-diameter model fitting. The calibrated height distribution models allow accurate results to be obtained with a very small sampling effort, making this approach highly effective and useful. Mixed effects models incorporate the variability between plots using the expression of the model's parameters in terms of both fixed and random effects. Random effects are conceptually random variables; they can be simulated as such, in terms of utilizing their distribution. To address this, we can also add a random component to the height prediction. This stochastic prediction approach uses distribution functions of random variable, H(d), and their confidence intervals. The stochastic predictions, h stoch,K , K 2 {V,G3,G4,B,G}, of a tree height can be defined as follows: