Next Article in Journal
Optimal Auxiliary Functions Method for a Pendulum Wrapping on Two Cylinders
Next Article in Special Issue
A Multivariate Hybrid Stochastic Differential Equation Model for Whole-Stand Dynamics
Previous Article in Journal
Rational Design of a Genetic Finite State Machine: Combining Biology, Engineering, and Mathematics for Bio-Computer Research
Previous Article in Special Issue
Non-Gaussian Quadrature Integral Transform Solution of Parabolic Models with a Finite Degree of Randomness
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Construction of Reducible Stochastic Differential Equation Systems for Tree Height–Diameter Connections

1
Agriculture Academy, Vytautas Magnus University, 53361 Kaunas, Lithuania
2
Department of Mathematics and Statistics, Vytautas Magnus University, 53361 Kaunas, Lithuania
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(8), 1363; https://doi.org/10.3390/math8081363
Submission received: 21 July 2020 / Revised: 6 August 2020 / Accepted: 10 August 2020 / Published: 14 August 2020
(This article belongs to the Special Issue Stochastic Differential Equations and Their Applications 2020)

Abstract

:
This study proposes a general bivariate stochastic differential equation model of population growth which includes random forces governing the dynamics of the bivariate distribution of size variables. The dynamics of the bivariate probability density function of the size variables in a population are described by the mixed-effect parameters Vasicek, Gompertz, Bertalanffy, and the gamma-type bivariate stochastic differential equations (SDEs). The newly derived bivariate probability density function and its marginal univariate, as well as the conditional univariate function, can be applied for the modeling of population attributes such as the mean value, quantiles, and much more. The models presented here are the basis for further developments toward the tree diameter–height and height–diameter relationships for general purpose in forest management. The present study experimentally confirms the effectiveness of using bivariate SDEs to reconstruct diameter–height and height–diameter relationships by using measurements obtained from mountain pine tree (Pinus mugo Turra) species dataset in Lithuania.

1. Introduction

The structural class hierarchy in a forest stand has interested researchers for more than a century. This study shows that the framework of the stochastic differential equations (SDEs) can easily be generalized to incorporate symmetric or non-symmetric size component distributions to explain empirical datasets in a theoretically consistent way.
Tree height (h) and tree diameter at breast height (d) (in the sequel—diameter) are traditionally used variables for tree volume calculation, aboveground carbon-stock estimation, and modeling of tree growth and yield. Unfortunately, it is difficult to measure tree height in the field, and errors tend to be observed. Tree height and diameter are traditionally formalized by using their regression relationship [1], the artificial neural network (ANN) [2,3], or stochastic differential equations [4,5]. Height–diameter regression equations have been developed for the local (stand) level and the generalized (ecoregional) level, by introducing additional stand variables, as well as random parameters [6,7]. The application of stand-level models to a wider region would probably lead to bias in predictions, as the deterministic height–diameter equations are related to the growth conditions and stand characteristics. The bias of tree-height predictions may affect the accuracy of the estimation of the stand volume and aboveground carbon stocks. To achieve reliable predictions of tree height, European forest statisticians have applied different techniques, such as generalized height–diameter models that include additional stand variables [8]. A generalized height equation evaluates the contribution of stand attributes, such as the stand density per hectare, the basal area per hectare, the quadratic mean diameter, and others in height–diameter models [9,10]. The inclusion of stand variables improves the predictive capacity of the selected height–diameter equation, but this technique requires additional sampling effort and reduces its practical application. The second technique that is used to introduce stand attributes to the height–diameter equation uses stand-specific effects defined for each stand by a normally distributed random variable, named a random effect [7], which allows models to adapt to diverse growth environments.
In evenly aged stands, the use of the previously listed approaches is appropriate for explaining and modeling tree development. On the other hand, the height–diameter relationship is not static; instead, it varies with stand age within the same stand. Despite the progress in linear and nonlinear regression modeling techniques, these earlier-proposed height–diameter models do not address fundamental factors of tree height and diameter development, such as tree age and covariance between size variables. In seeking a more comprehensive understanding of the dynamics of the tree-size variables, diffusion processes defined by SDEs have been proposed to express the shape of the tree height or crown width distribution versus the diameter [4,11]. SDE individual-tree growth models are a fundamental framework to link tree-size variables with stand age and to warrant mean stand volume predictions at a particular stand age. Most of the published SDE height–diameter models are applied to medium- and large-sized trees. However, the difficult challenge of forest growth and yield modeling is in tree species that differ from the idealized shape, size (diameter and height), age, and density. Previous SDE models have focused on bivariate (height and diameter) [12], trivariate (quadratic diameter, mean height, and density per hectare) [13], four-variate (height, diameter, crown base height, and crown width) [14,15,16], and five-variate (height, diameter, crown base height, crown width, and density per hectare) cases [17] for Scots pine trees in Lithuania. The main reason for developing SDE models is to gain the capacity to model highly nonlinear biological dynamics and their abnormalities [18].
The goal of this study was to develop SDE-type individual height–diameter models that are specific to the mountain pine tree (Pinus mugo Turra) species in Lithuania and to offer new insights into stand dynamics. The specific objectives were (1) to compare the newly developed SDE-type height–diameter equations with traditionally used regression equations; (2) to evaluate individual-tree and stand size variables such as the mean, median, mode, quantiles, and much more; (3) to evaluate and present the long-term behavior of predictive size variables; and (4) to account for between-stand variability by introducing random effects (random variables). All results were determined by using the symbolic algebra system Maple.

2. Materials and Methods

2.1. Height–Diameter Regression Models

One of the most important links between tree-size variables is the relationship between tree height and diameter. Estimating tree height from tree diameter at breast height is a common practice for forest inventory, model simulation, and forest management. Height–diameter equations enable a better understanding of the rates of changes in tree-size variables and the development of stand volume, forest biomass, and carbon stocks [19,20]. Previous studies have summarized a great number of available height–diameter linear and nonlinear relationships and compared their performance for different tree species [1,21]. It should be noted that the main limitation of the height–diameter equations is that they can produce very different results if applied to stands that differ from fitted stands. Therefore, mixed-effects modeling techniques were developed to account for stand variability [8,22]. The use of mixed-effects modeling techniques has made it possible to calibrate random effects to local stands by using a subsample of heights and diameters [22].
This study compared three nonlinear regression models named the power, Richards–Chapman, and q-exponential models [23,24,25]. The stand-specific random effects, φi (i = 1,…,M), were included in each examined model for only one of parameters that possessed higher variability. The general equations of the mixed-effects power, Richards–Chapman, and q-exponential models are as follows:
Model 1.Power function
h = ( β 1 + φ i ) d β 2 .
Model 2. Richards function
h = 1.3 + ( β 1 + φ i ) ( 1 e x p ( β 2 d ) ) β 3 .
Model 3. Q-exponential function
h = 1.3 + ( β 1 + φ i ) [ ( 1 β 2 ( 1 β 3 ) d ) ] + 1 1 β 3 .
where β1β3 are the unknown fixed-effect parameters to be estimated, [ a ] + = { a ,   i f   a 0 , 0 ,   i f   a < 0 and φi (i = 1,…,M) are normally distributed random variables with a mean of 0 and constant variance σ φ 2 , and M is the number of fitted stands.

2.2. SDE Models

More forest statisticians promote the use of bivariate distribution (diameter and height) for the derivation of the static height–diameter equation [26,27]. The mathematical description of the diameter or height distribution provides more detailed information about the stand [28]. It is worth noting that, by using measurements conducted by remote-sensing technologies, the height–diameter problem might be reversed, generating the need for a diameter–height relationship. A bivariate distribution is a framework used for symmetric interchange of diameter and height [29]. A bivariate SDE diameter and height stochastic process would be more appropriate for describing the development of the diameter and height structure in a stand over time (age). Therefore, we can define a dynamic form of the probability density function representing the diameter and height distribution, which changes in shape over time. The height–diameter or diameter–height equation can be derived by using conditional probability density functions. This study focused on four types of nonhomogeneous bivariate SDE models, X i ( t ) = ( D i ( t ) , H i ( t ) ) T , of the tree diameter, D i ( t ) , and height, H i ( t ) , dynamics, as follows [30]:
Model 4.Vasicek-type stochastic process
d X i ( t ) = A i ( X i ( t ) ) d t + B 1 2 · d W i ( t ) ,   P ( X i ( t 0 ) = x 0 ) = 1 ,   i = 1 , , M ,
A i ( x ) = ( β d ( α d + φ d i d ) , β h ( α h + φ h i h ) ) T ,   B = ( σ d d σ d h σ d h σ h h ) .
Model 5.Gompertz-type stochastic process
d X i ( t ) = A i ( X i ( t ) ) d t + D ( X i ( t ) ) B 1 2 · d W i ( t ) ,   P ( X i ( t 0 ) = x 0 ) = 1 ,   i = 1 , , M ,
A i ( x ) = ( ( ( α d + φ d i ) β d l n ( d ) ) d , ( ( α h + φ h i ) β h l n ( h ) ) h ) T ,   D ( x ) = ( d 0 0 h ) ,
G ( x ) = ( D ( x ) B 1 2 ) ( D ( x ) B 1 2 ) T = ( d 0 0 h ) ( σ d d σ d h σ d h σ h h ) ( d 0 0 h ) .
Model 6.Von Bertalanffy–type stochastic process
d X i ( t ) = A i ( X i ( t ) ) d t + D ( X i ( t ) ) B 1 2 · d W i ( t ) ,   P ( X i ( t 0 ) = x 0 ) = 1 ,   i = 1 , , M ,
A i ( x ) = ( ( α d + φ d i ) β d e β d ( t t 0 ) 1 d , ( α h + φ h i ) β h e β h ( t t 0 ) 1 h ) T ,   D ( x ) = ( d 0 0 h ) ,
G ( x ) = ( D ( x ) B 1 2 ) ( D ( x ) B 1 2 ) T = ( d 0 0 h ) ( σ d d σ d h σ d h σ h h ) ( d 0 0 h ) .
Model 7.Gamma-type stochastic process
d X i ( t ) = A i ( X i ( t ) ) d t + D ( X i ( t ) ) B 1 2 · d W i ( t ) ,   P ( X i ( t 0 ) = x 0 ) = 1 ,   i = 1 , , M ,
A i ( x ) = ( ( α d + φ d i t β d ) d , ( α h + φ h i t β h ) h ) T ,   D ( x ) = ( d 0 0 h ) ,
G ( x ) = ( D ( x ) B 1 2 ) ( D ( x ) B 1 2 ) T = ( d 0 0 h ) ( σ d d σ d h σ d h σ h h ) ( d 0 0 h ) .
Here, W i ( t ) = ( W 1 i ( t ) , W 2 i ( t ) ) T represents independent bivariate Brownian motions; t [ t 0 ; T ] is a finite horizon, T < ∞; φ d i and φ h i , i = 1, …, M are independent and normally distributed random variables with zero mean and constant variances ( φ d i ~ N ( 0 ; σ d 2 ) and φ h i ~ N ( 0 ; σ h 2 ) ); and { α d , α h , β d , β h , σ d d , σ d h , σ h h , σ d , σ h } are fixed-effect parameters to be estimated. In the sequel, the initial condition takes the following form: if t = t0, then X(t0) = x0 = (d0, h0).
The solution of the Vasicek-type SDE (4) has a bivariate normal distribution, N 2 ( μ i ( t ) ; Σ ( t ) ) , and the Gompertz-type, von Bertalanffy–type, and gamma-type SDEs (5)–(7) have a bivariate lognormal distribution, L N 2 ( μ i ( t ) ; Σ ( t ) ) , i = 1, … M, with the mean vectors μ i ( t ) = ( | μ d i ( t ) | μ h i ( t ) ) and the variance–covariance matrices Σ ( t ) = ( v d ( t ) v d h ( t ) v d h ( t ) v h ( t ) ) listed in Table 1.
The marginal distributions of D i ( t ) | D i ( t 0 ) = d 0 and H i ( t ) | H i ( t 0 ) = h 0 are also normal N 1 ( μ g i ( t ) ; v g ( t ) ) for the Vasicek type diffusion and lognormal L N 1 ( μ g i ( t ) ; v g ( t ) ) for the other diffusions ( g { d , h } ) [30]. The marginal mean, median, mode, p-quantile (0 < p < 1), and variance trajectories m g i ( t ) , m e g i ( t ) , m o g i ( t ) , m q g i ( t , p ) , and w g ( t ) of the tree diameter and height, g { d , h } , are listed in Table 2.
The conditional distributions D i ( t ) | D i ( t 0 ) = d 0 at a given ( H i ( t ) = h ) and H i ( t ) | H i ( t 0 ) = h 0 at a given ( D i ( t ) = d ) are univariate normal N 1 ( η g i ( t , g ¯ ) ; λ g ( t ) ) for the Vasicek-type diffusion and univariate lognormal L N 1 ( η g i ( t , g ¯ ) ; λ g ( t ) ) , g { d , h } for the other diffusions, with the mean and variance given in Table 3. Here, if g = d , then g ¯ = h , and if g = h , then g ¯ = d .
The conditional mean, median, mode, p-quantile (0 < p < 1), and variance trajectories m g i ( t , g ¯ ) , m e g i ( t , g ¯ ) , m o g i ( t , g ¯ ) , m q g i ( t , p , g ¯ ) , and w g i ( t , g ¯ ) of the tree diameter and height, g { d , h } , are listed in Table 4.

2.3. Maximum Likelihood Procedure

2.3.1. Regression Equation Models

For a sample { ( d 1 i , h 1 i ) , ( d 2 i , h 2 i ) ,   , ( d n i i , h n i i ) } of n ( n = i = i m n i ) observations on the dependent and independent variables, the maximum log-likelihood function takes the following form [31]:
L L 1 ( θ 1 ) = i = 1 M j = 1 n i l n ( 1 2 π σ e x p ( 1 2 σ ( g j i f ( g ¯ j i | β 1 , β 2 , β 3 , 0 ) ) 2 ) ) ,   θ 1 = { β 1 , β 2 , β 3 , σ } .
Here, the fixed-effect parameter nonlinear general regression model is defined as g = f ( g ¯ | β 1 + φ , β 2 , β 3 ) + ε = f ( g ¯ | β 1 , β 2 , β 3 , φ ) + ε , with the assumption that the error (ɛ) is a normally distributed random variable with zero mean and constant variance, σ2, and random effects, φ i = E ( φ i ) = 0 , i = 1 , , M .
For the mixed-effect parameter nonlinear general regression model g = f ( g ¯ | β 1 + φ , β 2 , β 3 ) + ε = f ( g ¯ | β 1 , β 2 , β 3 , φ ) + ε , with the assumption that the normal random error is ε ~ N 1 ( 0 ; σ 2 ) and the normal random effect is φ ~ N 1 ( 0 ; σ φ 2 ) , the two step approximated maximum log-likelihood function takes the following form [32]:
L L 2 ( θ 2 | φ i ) i = 1 M ( g ( φ i | θ 2 ) + 1 2 l n ( 2 π ) 1 2 l n ( d e t ( [ 2 g ( φ i | θ 2 ) ( φ i ) 2 ] ) | φ i = φ i ) ) ,  
where i = 1 , , M , φ i = a r g m a x φ i g ( φ i | θ 2 ) , θ 2 = { β 1 , β 2 , β 3 , σ , σ φ } ,
g ( φ i | θ 2 ) = j = 1 n i l n ( 1 2 π σ e x p ( 1 2 σ ( h j i f ( d j i | θ 1 , φ i ) ) 2 ) ) + l n ( ϕ ( φ i | σ φ ) ) .

2.3.2. SDE Models

The maximum likelihood parameter estimator of the SDE seeks to make the conditional (transition) probability density function the most likely fit for the observed diameter and height estimation dataset { ( d 1 i , h 1 i ) , ( d 2 i , h 2 i ) , , ( d n i i , h n i i ) } at discrete times { t 1 i , t 2 i , , t n i i } , starting at the initial time, t0, diameter, and height X ( t 0 ) = x 0 = ( 0.1 , 1.3 ) T , from which a given set of the SDEs (4)–(7) evolves. The parameter estimators for the fixed-effect model θ 3 = { α d , α h , β d , β h , σ d d , σ d h , σ h h } and for the mixed-effect model θ 4 = { α d , α h , β d , β h , σ d d , σ d h , σ h h , σ d , σ h } can be calculated by the maximum-likelihood and approximated-maximum-likelihood techniques. For the fixed-effect parameter scenario, the maximum log-likelihood function takes the following form:
L L 3 ( θ 3 ) = i = 1 M j = 1 n i l n ( p ( d j i , h j i | θ 3 , 0 ) ) ,
where p ( d , h | θ 3 , 0 ) is a normal (the Vasicek type) or lognormal (the others) probability density function, as defined in Table 1, with random effects φ i = E ( φ i ) = 0 .
For the mixed-effect parameter scenario models, the two-step approximated maximum log-likelihood procedure takes the following form:
L L 4 ( θ 4 | φ i ) i = 1 M ( g ( φ i | θ 4 ) + l n ( 2 π ) 1 2 l n ( d e t ( [ 2 g ( φ i | θ 4 ) ( φ i ) 2 ] ) | φ i = φ i ) ) ,
where φ i = ( φ d i , φ h i ) ,
φ i = a r g m a x φ i g ( φ i | θ 4 ) , i = 1 , , M ,
g ( φ i | θ 4 ) = j = 1 n i l n ( p ( d j i , h j i , t j i | θ 3 , φ d i , φ h i ) ) + l n ( ϕ ( φ d i | σ d ) + l n ( ϕ ( φ h i | σ h ) ) .

2.4. Standard Errors of Parameter Estimates

To assess the accuracy of the parameter estimates obtained through the maximum-likelihood and approximated-maximum-likelihood procedures, we considered the standard errors associated with the estimates. Estimates of these standard errors can be obtained as the inverse of the observed Fisher [33] information matrix (see, for instance, Reference [30]). The approximate asymptotic standard errors of the fixed-effects parameters are defined by the diagonal elements of the observed Fisher information matrix [ I ˜ ( θ s ) ] 1 , s = 1,2,3,4
S E ( θ s ) = [ I i i ˜ ( θ s ) ] 1 ,
where the matrix is I ˜ ( θ s ) = [ 2 L L s ( θ s ) θ i s θ j s ] T | θ s = θ s .

2.5. Random Effects Calibration

To calibrate the random effects, φ d , φ h , of the bivariate SDE models defined by Equations (4)–(7) for a new subsample, we used a full dataset of observations taken from the validation sample unit. Using new subsample observations { ( d 1 , h 1 ) , ( d 2 , h 2 ) , , ( d m , h m ) } at discrete times { t 1 , t 2 , , t m } , the random effects of the SDE models φ = ( φ d , φ h ) for diameter and height can be calibrated by the following:
φ = a r g m a x ( φ d , φ h ) j = 1 m l n ( p ( d j , h j , t j | θ 3 , φ d , φ h ) ) + l n ( ϕ ( φ d | σ d ) ) + l n ( ϕ ( φ h | σ h ) ) .
The random effect φg of the regression models can be calibrated by the following:
φ g = a r g m a x φ g j = 1 m l n ( 1 2 π σ e x p ( 1 2 σ ( g j f ( g ¯ j | θ 1 , φ ) ) 2 ) ) + l n ( ϕ ( φ g | σ φ ) ) ,   g { d , h } .  

2.6. Data

The study site was located in western Lithuania (Kuršių Nerija). A mountain pine tree (Pinus mugo Turra) species was planted to stop sand dunes of Curonian Spit from drifting in the nineteenth to twentieth centuries. The data used for modeling consisted of measurements collected from 31 mountain pine plots located in Western Lithuania (Kuršių Nerija). The ages of the sampled plots varied from 53 to 123 years. The size of the round sample plot was 150 m2. In all plots, the diameter at breast height of all trees larger than 1 mm was measured (7005 trees). The total diameter at breast height and tree height was measured for 702 trees. At plot establishment, the following data were recorded for every sample tree: diameter over bark at 1.30 m, tree height, and age. Diameter was measured to an accuracy of 1 mm. Height was measured by using clinometer, with precision to the nearest 0.1 m. Tree age was identified with stand age, which provides a general timeframe for when stands were first established. The dataset was randomly divided into estimation and validation datasets. A random sample of 23 plots was selected for model estimation, and the remaining dataset of 8 plots was used for model validation. The summary statistics for the diameter at breast height (d), the total height (h), and age (t) for all of the trees used in the model estimation and validation datasets are presented in Table 5.

3. Results and Discussion

3.1. Estimates of Parameters

All used fixed and mixed-effect parameter regression models described above, and the newly developed bivariate fixed and mixed-effect parameter SDE models, are nonlinear. In order to establish the mixed-effect parameter height–diameter and diameter-height models, the sample plot-level effects were considered by including the random effects for a particular parameter. This study estimated the parameters of all used regression and bivariate SDE models by maximizing the maximum likelihood functions defined by Equations (8)–(12). Standard errors of parameter estimates were calculated by using the observed Fisher information matrix, defined by Equation (15). All results were determined by using Maple mathematical software. The results of the parameter estimation and the standard errors for all used models are listed in Table 6 and Table 7. Parameter estimates of all fitted models were significant (p < 0.05).

3.2. Bivariate Diameter and Height Distributions

Newly developed bivariate probability density functions obtained from the bivariate Vasicek-, Gompertz-, von Bertalanffy–, and gamma-type SDEs defined by Equations (4)–(7) can be visualized in graphs for different stand and parameter scenarios. Additionally, to demonstrate that data from the validation dataset follow the bivariate probability density functions derived from SDEs (4)–(7), we used a simple graphical technique represented by contour plots. The random effects, φ d ,   φ h , of the diameter and height for a randomly selected stand from the validation dataset were calibrated by Equation (16). The estimated bivariate mixed-effect parameter probability density functions and contour plots for a particular randomly selected stand from the validation dataset and the observed values are given in Figure 1. We can see that our developed lognormal probability density functions for the Gompertz, von Bertalanffy, and gamma diffusions inherit circular structures but also have unequal tails and show signs of skewness.
For mountain pine trees, the height and diameter probability density function presented steeper and more symmetrical surfaces for the Vasicek-type diffusions (see Figure 1); however, with increasing age, the surfaces became flatter.

3.3. Comparison of Diameter, Height, Height–Diameter, and Diameter–Height Models

Tree size dimensions vary continuously through natural processes and specific silvicultural activities. Accurate tree-height relationships are necessary for most forest inventories, as knowledge of tree heights is needed to calculate tree volumes, biomass, carbon content, and economic value [34]. The newly developed fixed- and mixed-effect marginal and conditional univariate normal and lognormal probability density functions, which are listed in Table 1 and Table 2, respectively, enable the results of the diameter and height behavior to be formalized from an evolutionary (time) perspective, with an emphasis on models of the mean diameter and mean height in a forest stand. By using marginal and conditional univariate normal and lognormal probability density functions, trajectories of the mean tree diameter and height in a forest stand were derived (see Table 2 and Table 4). Fixed-effect parameter height–diameter SDE models are suitable for local application, and more generalized mixed-effect parameter SDE models, which do not use additional stand variables, can be applied at the regional level.
The principal objective of the newly framed SDE height–diameter and diameter–height models, as presented in Table 4, is to show their precedence over regression models. The statistical measures commonly used for the comparison of all used height–diameter and diameter–height models are shown in Table 8 and Table 9. The SDE height–diameter and diameter–height models were found to be better than the regression models. The mixed-effect parameter SDE models showed better statistical measures than those of fixed-effect parameter models, as, by definition, they do not take into account the variability between stands and therefore do not provide good predictions for a dataset sampled from a new stand. Statistical measures for all fitted height–diameter and diameter–height SDE models gave very similar values. For the validation dataset, the SDE height–diameter and diameter–height models appeared to work better than the regression models. With reference to all used statistical measures and to all used height–diameter and diameter–height models, the bivariate Gompertz-type SDE, in general terms, provided the highest accuracy.
Figure 2 shows plots of the residuals against the predicted values resulting from the validation dataset. The p-values (entered up in the plots) of the Shapiro–Wilk test for the residuals demonstrate that the residuals follow a normal distribution. The visual examination of the residual plots presented in Figure 2 does not show the presence of systematic errors (bias). Based on the statistical measures presented in Table 8 and Table 9, the residual analysis of the mixed-effect parameter height–diameter and diameter–height SDE models, the impossibility of the negative values, and the asymmetry of the tree size component distribution, it was concluded that the Gompertz-type diffusion is superior to other diffusions. On the other side, the four mean curves from the SDE models did not show a notable difference.
For the mixed-effect parameter Gompertz-type SDE using the marginal diameter and height distributions presented in Table 1, the mean, mode, median, 5% quantile, 95% quantile, and lower and upper quartile trajectories (see Table 2) for two randomly selected plots from the validation dataset are presented in Figure 3 (random effects were calibrated by Equation (16)). Figure 3 demonstrates that the stochastic bivariate Gompertz-type differential equations with a sigmoidal form solution are more reliable for future diameter and height predictions and allow better interpretability of model parameters than nonlinear regression models. An essential research finding pertaining to the tree-diameter and -height distributions is the fact that they are positively skewed (see Figure 3). Figure 3 illustrates that the tree diameter distribution is more asymmetric than the tree-height distribution. Theoretical studies on the growth of the tree-size components proved that the tree-size-component frequency distribution is characterized by many small trees and few larger trees [35].
For the bivariate Gompertz-type diffusion process, it is possible to derive the stationary univariate marginal and conditional distributions if parameters βd and βh are positive [36]. As time, t, goes to infinity for the Gompertz-type diffusion, the diameter and height marginal distributions (see Table 1) and conditional distributions (see Table 3) converge to a lognormal stationary distribution with the means and variances listed in Table 10 ( g { d , h } ). The asymptotic mean, median, mode, p-quantile (0 < p < 1), and variance trajectories of the tree diameter and height marginal and conditional processes are listed in Table 11. The stationary conditional lognormal type distribution derived from the Gompertz SDE (see Table 10) enables us to estimate the full distribution of the dependent variable (diameter or height) when the value of the independent variable (height or diameter) is known. This technique is a flexible rule that can be used to describe the development of the diameter or height distribution against the height or diameter. One advantage of the stationary conditional distribution is that the conditional quantile functions, and different measures of central tendency and variability can be obtained, providing a more comprehensive analysis of the relationship between diameter and height. Until now, the quantile regression has been widely conducted in tree-variable modeling [10,16,37].
By using the stationary conditional diameter and height lognormal distributions presented in Table 10, we present the mean, mode, median, 5% quantile, 95% quantile, and lower and upper quartile trajectories (see Table 11) for two stands randomly selected from the validation dataset in Figure 4 (random effects were calibrated by Equation (16).

3.4. Mean Tree Basal Area Models

The mean tree basal area, BA, is one of the central attributes of a forest stand. It is simply the average of the cross-sectional area at breast height of all trees in a forest stand. In order to better understand how the mean tree basal area is changing over time, we used the Gompertz diffusion process framework developed in previous subsections. The mean tree basal area can be used to estimate the mean tree volume and to generate management procedures. The mean tree basal area dynamics describe tree development with age and are defined by the following:
B A ( t ) = π 400 0 + x 2 p ( x , t | α d , β d , σ d d , σ d , φ d ) d x ,
where p ( · , t | α d , β d , σ d d , σ d , φ d ) is the marginal probability distribution density of the tree diameter at a given age, t. The probability density functions (pdf) can be divided into four scenarios:
p ( x , t ) = { | p d f   o f   L N 1 ( μ d i ( t ) ;   v d d ( t ) ) ,   t r a n s i t i o n   d e n s i t y | p d f   o f   L N 1 ( η d i ( t , h ) ;   λ d ( t ) ) ,   t r a n s i t i o n   d e n s i t y | p d f   o f   L N 1 ( μ d i ;   v d d ) ,   s t a t i o n a r y   d e n s i t y | p d f   o f   L N 1 ( η d i ( h ) ;   λ d ) ,   s t a t i o n a r y   d e n s i t y ,
where the moments μ d i ( t ) , v d d ( t ) , η d i ( t , h ) , λ d ( t ) , μ d i , v d d , and η d i ( h ) , λ d are listed in Table 1, Table 3, Table 10 and Table 11, respectively.
An important contribution to the development of tree basal area prediction methodology is compatibility between its growth and increment equations, which can be derived from the growth rate equation by differentiation. Traditionally, the majority of authors have used a system of permanent experimental plots to provide partial time-series data covering different stand densities and thinning strategies [38]. Using the mean tree basal area growth, Equation (18), the current and mean annual increments of the tree basal area (CAIB and MAIB, respectively) can be defined in the following forms:
C A I B ( t ) = d d t B ( t ) ,
M A I B ( t ) = 1 t B ( t ) .
According to the mean diameter relationships listed in Table 2, the current and mean annual increments of the diameter (CAID and MAID, respectively) have the following forms:
C A I D ( t ) = d d t 0 + x · p ( x , t | α d , β d , σ d d , σ d , φ d ) · d x ,
M A I D ( t ) = 1 t 0 + x · p ( x , t | α d , β d , σ d d , σ d , φ d ) · d x .
Figure 5 illustrates the mean basal area dynamics, over time, with the observed values for two randomly selected stands from the validation dataset, using the mixed-effect scenario (the random effects were calibrated by Equation (16) and for all stands from the validation dataset, using the fixed-effect scenario). Figure 6 provides graphs of the basal area and diameter increments (acceleration) over time. We should note that all increments illustrated in Figure 6 are positive; therefore, they correspond to mathematical principles of the derivative, because trees increase in basal area and diameter for as long as they live. Figure 6 demonstrates that maximum basal area growth occurs later than the maximum diameter growth. Our conclusion is not related to the selection of a specific function to model the basal area or diameter increments. This is the fundamental quality of the used modeling technique that deals with probability density functions of the tree diameter and height over time. The maximum values of CAI and MAI, defined by Equations (20)–(23), occur at different times, according to the quality of the stand.
Table 12 shows the goodness-of-fit statistics for tree basal area predictions. From Table 12, it is clear that the mixed-effect scenario model is more accurate for predicting the tree basal area than the fixed-effect scenario model. Statistical indexes for the estimation and validation datasets produced similar values. The best statistical index values were obtained with the mixed-effect parameter model, which includes the tree height as an additional predictor. The mixed-effect parameter stationary models produced statistical index values very close to the non-stationary process.

3.5. Mean Tree Volume Models

Prediction of growth and yield in trees are important for decision-making in sustainable forest management. In order to estimate the mean tree volume, the diameter at breast height and total height variables are traditionally used as predictors. The newly developed bivariate lognormal probability density function of the diameter and height distributions (see Table 1) and an additional stem volume regression function allow us to assess the mean tree volume and its dynamics in the forward and backward directions, as follows:
V ¯ ( t ) = 0 + 0 + V ( x , y ) · p ( x , y , t | θ 4 , φ d , φ h ) · d x · d y ,
where, in this study, V ( d , h ) is the individual tree volume regression function of the q-exponential form [25]:
V = β 1 + β 2 h β 3 [ β 4 ( 1 e x p ( ( 1 β 5 ) d ) ) ] + 1 1 β 5 .
The estimates of the parameters β1, β2, and β3 were calculated by the weighted least-squares technique. The estimators were β 1 0.0007 , β 2 0.7181 , β 3 0.0013 , β 4 0.2658 , and β 5 - 314.9569 . The observed mean stem volumes for all stands were computed by regression, using Equation (25).
The SDE framework allows the construction of dynamic volume growth models, which offer a better description of the development of a stand over time than static models and can be considered superior for prediction purposes. Additionally, the dynamic growth model provides forest decision-makers with more detailed information about the volumes of merchantable products. The variations in the mean tree volume with age and among stands are shown in Figure 7. Insight from diffusion processes allowed us to formalize a unique way to study the long-term dynamics of the mean tree volume at the stand level. Figure 8 provides a graph of the mean tree volume increments over time. Figure 6 and Figure 8 demonstrate that the maximum individual tree volume growth occurs at a later age than the maximum basal area and diameter growth. Table 13 shows the goodness-of-fit statistics for mean tree volume predictions. From Table 13, it is clear that the mixed-effect scenario bivariate SDE model can accurately predict the mean tree volume. Statistical indexes for both estimation and validation datasets produced similar values.

4. Conclusions

SDEs were developed at the beginning of the twentieth century, to quantify all aspects of stochastic processes. This study evaluated the use of SDEs to model the tree diameter and height at any given age in mountain pine tree (Pinus mugo Turra) species of Lithuania. We developed a few new models for diameter and height evolution by using well-defined diffusion processes, such as the symmetric Vasicek diffusion process and asymmetric Bertalanffy-, Gompertz- and gamma-diffusion processes. In both scenarios, fixed- and mixed-effect parameters were applied to model the evolution of the diameter and height. We performed fixed- and mixed-effect parameters estimation via maximum likelihood procedure and using a real-world dataset of mountain pine trees in Lithuania. These models were compared with traditionally used regression models by using performance statistics and residual analysis. Overall, the best goodness-of-fit statistics were produced by the Gompertz-type diffusion model. On the other side, the superiority of the Gompertz-type diffusion can been confirmed by the existence of the steady state variance, the impossibility of the negative values, and the asymmetry of the tree-size component distribution. All results were determined by using the Maple computer algebra system.

Author Contributions

Conceptualization, M.N. and P.R.; methodology, P.R.; software, P.R.; validation, M.N., E.P., and P.R.; data curation, M.N. and E.P.; writing—original draft preparation, M.N., E.P., and P.R.; writing—review and editing, M.N., E.P., and P.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zea-Camaño, J.D.; Soto, J.R.; Arce, J.E.; Pelissari, A.L.; Behling, A.; Orso, G.A.; Guachambala, M.S.; Eisfeld, R.L. Improving the Modeling of the Height–Diameter Relationship of Tree Species with High Growth Variability: Robust Regression Analysis of Ochroma pyramidale (Balsa-Tree). Forests 2020, 11, 313. [Google Scholar] [CrossRef] [Green Version]
  2. Shen, J.; Hu, Z.; Sharma, R.P.; Wang, G.; Meng, X.; Wang, M.; Wang, Q.; Fu, L. Modeling Height–Diameter Relationship for Poplar Plantations Using Combined-Optimization Multiple Hidden Layer Back Propagation Neural Network. Forests 2020, 11, 442. [Google Scholar] [CrossRef] [Green Version]
  3. Tavares Júnior, I.S.; Rocha, J.E.C.; Ebling, Â.A.; Chaves, A.S.; Zanuncio, J.C.; Farias, A.A.; Leite, H.G. Artificial Neural Networks and Linear Regression Reduce Sample Intensity to Predict the Commercial Volume of Eucalyptus Clones. Forests 2019, 10, 268. [Google Scholar] [CrossRef] [Green Version]
  4. Rupšys, P. New insights into tree height distribution based on mixed-effects univariate diffusion processes. PLoS ONE 2016, 11, e0168507. [Google Scholar] [CrossRef]
  5. Rupšys, P.; Petrauskas, E. A New Paradigm in Modelling the Evolution of a Stand via the Distribution of Tree Sizes. Sci. Rep. 2017, 7, 15875. [Google Scholar] [CrossRef] [Green Version]
  6. Bronisz, K.; Mehtätalo, L. Mixed-Effects Generalized Height–Diameter Model for Young Silver Birch Stands on Post-Agricultural Lands. For. Ecol. Manag. 2020, 460, 117901. [Google Scholar] [CrossRef]
  7. Xie, L.; Widagdo, F.R.A.; Dong, L.; Li, F. Modeling Height–Diameter Relationships for Mixed-Species Plantations of Fraxinus mandshurica Rupr. and Larix olgensis Henry in Northeastern China. Forests 2020, 11, 610. [Google Scholar] [CrossRef]
  8. Stankova, T.V.; Diéguez-Aranda, U. Height–diameter relationships for Scots pine plantation in Bulgaria: Optimal combination of model type and application. Ann. For. Res. 2013, 56, 149–163. [Google Scholar]
  9. Temesgen, H.; Gadow, K.V. Generalized height–diameter models–an application for major tree species in complex stands of interior British Columbia. Eur. J. For. Res. 2004, 123, 45–51. [Google Scholar] [CrossRef]
  10. Zhang, B.; Sajjad, S.; Chen, K.; Zhou, L.; Zhang, Y.; Yong, K.K.; Sun, Y. Predicting Tree Height–Diameter Relationship from Relative Competition Levels Using Quantile Regression Models for Chinese Fir (Cunninghamia lanceolata) in Fujian Province, China. Forests 2020, 11, 183. [Google Scholar] [CrossRef] [Green Version]
  11. Rupšys, P. Stochastic Mixed-Effects Parameters Bertalanffy Process, with Applications to Tree Crown Width Modeling. Math. Probl. Eng. 2015, 2015, 375270. [Google Scholar] [CrossRef] [Green Version]
  12. Rupšys, P.; Petrauskas, E. Evolution of Bivariate Tree Diameter and Height Distribution via Stand Age: Von Bertalanffy Bivariate Diffusion Process Approach. J. Forest Res.-Jap. 2019, 24, 16–26. [Google Scholar] [CrossRef]
  13. Rupšys, P. Modeling Dynamics of Structural Components of Forest Stands Based on Trivariate Stochastic Differential Equation. Forests 2019, 10, 506. [Google Scholar] [CrossRef] [Green Version]
  14. Rupšys, P. The Use of Copulas to Practical Estimation of Multivariate Stochastic Differential Equation Mixed-Effects Models. AIP Conf. Proc. 2015, 1684, 080011. [Google Scholar]
  15. Rupšys, P.; Petrauskas, E. A Linkage among Tree Diameter, Height, Crown Base Height, and Crown Width 4-variate Distribution and Their Growth Models: A 4-variate Diffusion Process Approach. Forests 2017, 8, 479. [Google Scholar] [CrossRef] [Green Version]
  16. Rupšys, P. Understanding the Evolution of Tree Size Diversity within the Multivariate Nonsymmetrical Diffusion Process and Information Measures. Mathematics 2019, 7, 761. [Google Scholar] [CrossRef] [Green Version]
  17. Rupšys, P. Modeling Perspectives of Forest Growth and Yield: Framework of Multivariate Diffusion Process. AIP Conf. Proc. 2019, 2164, 060017. [Google Scholar]
  18. Petrauskas, E.; Rupšys, P.; Narmontas, M.; Aleinikovas, M.; Beniušienė, L.; Šilinskas, B. Stochastic Models to Qualify Stem Tapers. Algorithms 2020, 13, 94. [Google Scholar] [CrossRef] [Green Version]
  19. Deng, C.; Zhang, S.; Lu, Y.; Froese, R.E.; Ming, A.; Li, Q. Thinning Effects on the Tree Height–Diameter Allometry of Masson Pine (Pinus massoniana Lamb.). Forests 2019, 10, 1129. [Google Scholar] [CrossRef] [Green Version]
  20. Liu, G.; Wang, J.; Dong, P.; Chen, Y.; Liu, Z. Estimating Individual Tree Height and Diameter at Breast Height (DBH) from Terrestrial Laser Scanning (TLS) Data at Plot Level. Forests 2018, 9, 398. [Google Scholar] [CrossRef] [Green Version]
  21. Mensah, S.; Pienaar, O.L.; Kunneke, A.; Dutoit, B.; Seydack, A.; Uhl, E.; Pretzsch, H.; Seifert, T. Height Diameter allometry in South Africa’s indigenous high forests: Assessing generic models performance and function forms. For. Ecol. Manag. 2018, 410, 1–11. [Google Scholar] [CrossRef]
  22. Dean, W.C.; Lee, Y.J. A Mixed-Effects Height–Diameter Model for Individual Loblolly and Slash Pine Trees in East Texas. South. J. Appl. For. 2011, 35, 12–17. [Google Scholar]
  23. Huxley, A. Problems of Relative Growth; The Dial Press: New York, NY, USA, 1932. [Google Scholar]
  24. Richards, F.J. A flexible growth function for empirical use. J. Exp. Bot. 1959, 10, 290–300. [Google Scholar] [CrossRef]
  25. Petrauskas, E.; Rupšys, P.; Memgaudas, R. Q-Exponential Variable-form of a Steam Taper and Volume Model for Scots Pine (Pinus sylvesteris L.) in Lithuania. Baltic For. 2011, 17, 118–127. [Google Scholar]
  26. Xingji, J.; Fengri, L.; Weiwei, J.; Lianjun, Z. Modeling and Predicting Bivariate Distributions of Tree Diameter and Height. Sci. Silvae Sin. 2013, 49, 74–82. [Google Scholar]
  27. Ogana, F.N.; Osho, J.S.A.; Gorgoso-Varela, J.J. An approach to modeling the joint distribution of tree diameter and height data. J. Sustain. Forest. 2018, 37, 475–488. [Google Scholar] [CrossRef]
  28. Pogoda, P.; Ochał, W.; Orzeł, S. Performance of Kernel Estimator and Johnson SB Function for Modeling Diameter Distribution of Black Alder (Alnus glutinosa (L.) Gaertn.) Stands. Forests 2020, 11, 634. [Google Scholar] [CrossRef]
  29. Mønness, E. The bivariate power-normal distribution and the bivariate Johnson system bounded distribution in forestry, including height curves. Can. J. For. Res. 2015, 45, 307–313. [Google Scholar] [CrossRef]
  30. Rupšys, P. Univariate and Bivariate Diffusion Models: Computational Aspects and Applications to Forestry. In Stochastic Differential Equations: Basics and Applications; Tony, G., Deangelo, T.G., Eds.; Nova Science: New York, NY, USA, 2018; pp. 1–77. [Google Scholar]
  31. Davidian, M.; Giltinan, D.M. Some simple methods for estimating intra-individual variability in nonlinear mixed-effects model. Biometrics 1993, 49, 59–73. [Google Scholar] [CrossRef]
  32. Joe, H. Accuracy of Laplace Approximation for Discrete Response Mixed Models. Comput. Stat. Data An. 2008, 52, 5066–5074. [Google Scholar] [CrossRef]
  33. Fisher, R.A. On the Mathematical Foundations of Theoretical Statistics. Philos. T. Roy. Soc. A 1922, 222, 309–368. [Google Scholar]
  34. MacPhee, C.; Kershaw, J.A.; Weiskittel, A.R.; Golding, J.; Lavigne, M.B. Comparison of Approaches for Estimating Individual Tree Height–Diameter Relationships in the Acadian Forest Region. Forestry 2018, 91, 132–146. [Google Scholar] [CrossRef]
  35. Ishihara, M.I.; Konno, Y.; Umeki, K.; Ohno, Y.; Kikuzawa, K. A New Model for Size-Dependent Tree Growth in Forests. PLoS ONE 2016, 11, e0152219. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Ramos-Ábalos, E.M.; Gutiérrez-Sánchez, R.; Nafidi, A. Powers of the Stochastic Gompertz and Lognormal Diffusion Processes, Statistical Inference and Simulation. Mathematics 2020, 8, 588. [Google Scholar] [CrossRef] [Green Version]
  37. Narmontas, M.; Rupšys, P.; Petrauskas, E. Models for Tree Taper Form: The Gompertz and Vasicek Diffusion Processes Framework. Symmetry 2020, 12, 80. [Google Scholar] [CrossRef] [Green Version]
  38. Burkhart, H.E.; Tomé, M. Modelling Forest Trees and Stands; Springer Sience + Business Media: Dordrecht, The Netherlands, 2012. [Google Scholar]
Figure 1. Bivariate probability density functions and their contour plots for a particular randomly selected stand from the validation dataset: (a1,b1,c1,d1) mixed-effect bivariate densities, (a2,b2,c2,d2) Contour plots (levels: 5−2, 5−3, 5−4), (a1,a2) Vasicek type, (b1,b2) Gompertz type, (c1,c2) Bertalanffy type, and (d1,d2) Gamma type. Values observed from the validation dataset are represented by circles.
Figure 1. Bivariate probability density functions and their contour plots for a particular randomly selected stand from the validation dataset: (a1,b1,c1,d1) mixed-effect bivariate densities, (a2,b2,c2,d2) Contour plots (levels: 5−2, 5−3, 5−4), (a1,a2) Vasicek type, (b1,b2) Gompertz type, (c1,c2) Bertalanffy type, and (d1,d2) Gamma type. Values observed from the validation dataset are represented by circles.
Mathematics 08 01363 g001
Figure 2. Residuals of the diameter–height (left) and height–diameter (right) models. The solid black line shows the loess regression, and observed values from the validation dataset are represented by boxes.
Figure 2. Residuals of the diameter–height (left) and height–diameter (right) models. The solid black line shows the loess regression, and observed values from the validation dataset are represented by boxes.
Mathematics 08 01363 g002
Figure 3. Dynamics of the mean, mode, median, 5% quantile, 95% quantile, and lower- and upper-quartile trajectories for two stands within the validation dataset: top—the first stand; bottom—the second stand; mean—solid line in black; median—solid line in red; mode—solid line in blue; 5% and 95% quantiles—dotted lines; lower and upper quartiles—dashed lines; values are shown as circles.
Figure 3. Dynamics of the mean, mode, median, 5% quantile, 95% quantile, and lower- and upper-quartile trajectories for two stands within the validation dataset: top—the first stand; bottom—the second stand; mean—solid line in black; median—solid line in red; mode—solid line in blue; 5% and 95% quantiles—dotted lines; lower and upper quartiles—dashed lines; values are shown as circles.
Mathematics 08 01363 g003
Figure 4. Dynamics of the stationary conditional mean, mode, median, 5% quantile, 95% quantile, and lower- and upper-quartile trajectories for two stands within the validation dataset: left panel—diameter–height models; right panel—height–diameter models; (a1,a2)—the first stand mixed-effect scenario; (b1,b2)—the second stand mixed-effect scenario; (c1,c2)—the fixed-effect scenario for all stands from the validation dataset; mean—solid line in black; median—solid line in red; mode—solid line in blue; 5% and 95% quantiles—dotted lines; lower and upper quartiles—dashed lines; values are shown by circles.
Figure 4. Dynamics of the stationary conditional mean, mode, median, 5% quantile, 95% quantile, and lower- and upper-quartile trajectories for two stands within the validation dataset: left panel—diameter–height models; right panel—height–diameter models; (a1,a2)—the first stand mixed-effect scenario; (b1,b2)—the second stand mixed-effect scenario; (c1,c2)—the fixed-effect scenario for all stands from the validation dataset; mean—solid line in black; median—solid line in red; mode—solid line in blue; 5% and 95% quantiles—dotted lines; lower and upper quartiles—dashed lines; values are shown by circles.
Mathematics 08 01363 g004
Figure 5. Dynamics of the mean tree basal area at different ages. (M) Mixed-effects scenario for two stands; (F) fixed-effects scenario for all stands; solid line—mean tree basal area curve determined by using the marginal probability density function; dotted line—mean tree basal area determined by using the conditional probability density function where the height is equal to the mean height of the stand; black in (M)—the first stand; red in (M and F)—the second stand.
Figure 5. Dynamics of the mean tree basal area at different ages. (M) Mixed-effects scenario for two stands; (F) fixed-effects scenario for all stands; solid line—mean tree basal area curve determined by using the marginal probability density function; dotted line—mean tree basal area determined by using the conditional probability density function where the height is equal to the mean height of the stand; black in (M)—the first stand; red in (M and F)—the second stand.
Mathematics 08 01363 g005
Figure 6. Dynamics of the current and mean annual tree basal area increments of the basal area and diameter, over time, for two randomly selected stands from the validation dataset and mixed-effects scenario: (a) basal area increments; (b) diameter increments; black color represents the first stand; red color represents the second stand; solid line—the current increment; dotted line—the mean annual increment.
Figure 6. Dynamics of the current and mean annual tree basal area increments of the basal area and diameter, over time, for two randomly selected stands from the validation dataset and mixed-effects scenario: (a) basal area increments; (b) diameter increments; black color represents the first stand; red color represents the second stand; solid line—the current increment; dotted line—the mean annual increment.
Mathematics 08 01363 g006
Figure 7. Dynamics of the mean tree volume, over time, for two randomly selected stands from the validation dataset and mixed-effect scenario. Black color—the first stand; red color—the second stand; circles—the observed values.
Figure 7. Dynamics of the mean tree volume, over time, for two randomly selected stands from the validation dataset and mixed-effect scenario. Black color—the first stand; red color—the second stand; circles—the observed values.
Mathematics 08 01363 g007
Figure 8. Dynamics of the current and mean annual tree volume increments versus age for two randomly selected stands from the validation dataset and mixed-effect scenario. Black color—the first stand; red color—the second stand.
Figure 8. Dynamics of the current and mean annual tree volume increments versus age for two randomly selected stands from the validation dataset and mixed-effect scenario. Black color—the first stand; red color—the second stand.
Mathematics 08 01363 g008
Table 1. Mean vectors and variance–covariance matrixes of bivariate probability density functions.
Table 1. Mean vectors and variance–covariance matrixes of bivariate probability density functions.
SDE 1 Mean   Vector   μ ( t ) Variance Covariance   Matrix   Σ ( t )
Vasicek ( | d 0 e β d ( t t 0 ) + α d + φ d i β d ( 1 e β d ( t t 0 ) ) | h 0 e β h ( t t 0 ) + α h + φ h i β h ( 1 e β h ( t t 0 ) ) ) ( 1 e 2 β d ( t t 0 ) 2 β d σ d d 1 e ( β d + β h ) ( t t 0 ) β d + β h σ d h 1 e ( β d + β h ) ( t t 0 ) β d + β h σ d h 1 e 2 β h ( t t 0 ) 2 β h σ h h )
Gompertz ( | l n ( d 0 ) e β d ( t t 0 ) + α d + φ d i β d ( 1 e β d ( t t 0 ) ) | l n ( h 0 ) e β h ( t t 0 ) + α h + φ h i β h ( 1 e β h ( t t 0 ) ) )
Bertalanffy ( | l n ( d 0 ) + ( α d + φ d i ) l n ( 1 e x p ( β d t ) 1 e x p ( β d t 0 ) ) σ d d 2 ( t t 0 ) | l n ( h 0 ) + ( α h + φ h i ) l n ( 1 e x p ( β h t ) 1 e x p ( β h t 0 ) ) σ h h 2 ( t t 0 ) ) ( σ d d ( t t 0 ) σ d h ( t t 0 ) σ d h ( t t 0 ) σ h h ( t t 0 ) )
Gamma ( | l n ( d 0 ) + ( α d + φ d i ) l n ( t t 0 ) ( β d + σ d d 2 ) ( t t 0 ) | l n ( h 0 ) + ( α h + φ h i ) l n ( t t 0 ) ( β h + σ h h 2 ) ( t t 0 ) )
1: SDE = stochastic differential equation.
Table 2. Marginal mean, median, mode, p-quantile (0 < p < 1), and variance trajectories.
Table 2. Marginal mean, median, mode, p-quantile (0 < p < 1), and variance trajectories.
SDETrajectory TypeEquation
VasicekMean, median and mode μ g i ( t )
Quantile (0 < p < 1) Φ p 1 ( μ g i ( t ) ; v g g ( t ) ) 1
Variance v g g ( t )
Gompertz
Bertalanffy
Gamma
Mean e x p ( μ g i ( t ) + 1 2 v g g ( t ) )
Median e x p ( μ g i ( t ) )
Mode e x p ( μ g i ( t ) v g g ( t ) )
Quantile (0 < p < 1) e x p ( μ g i ( t ) + v g g ( t ) Φ p 1 ( 0 ; 1 ) )  
Variance e x p ( 2 μ g i ( t ) + v g g ( t ) ) · ( e x p ( v g g ( t ) ) 1 )
1: Φ p 1 ( · ; · ) is the inverse of the standard normal distribution function.
Table 3. Mean and variance of conditional probability density functions.
Table 3. Mean and variance of conditional probability density functions.
SDE Mean   η g i ( t , g ¯ ) Variance   λ g ( t )
Vasicek μ g i ( t ) + v g g ¯ ( t ) v g ¯ g ¯ ( t ) ( h μ g ¯ i ( t ) ) v g g ( t ) ( v g g ¯ ( t ) ) 2 v g ¯ g ¯ ( t )
Gompertz
Bertalanffy
Gamma
μ g i ( t ) + v g g ¯ ( t ) v g ¯ g ¯ ( t ) ( l n ( g ¯ ) μ g ¯ i ( t ) )
Table 4. Marginal mean, median, mode, p-quantile (0 < p < 1), and variance trajectories.
Table 4. Marginal mean, median, mode, p-quantile (0 < p < 1), and variance trajectories.
SDETrajectory TypeEquation
VasicekMean, median, and mode η g i ( t , g ¯ )
Quantile (0 < p < 1) Φ p 1 ( η g i ( t , g ¯ ) ; λ g ( t ) )
Variance λ g ( t )
Gompertz
Bertalanffy
Gamma
Mean e x p ( η g i ( t , g ¯ ) + 1 2 λ g ( t ) )
Median e x p ( η g i ( t , g ¯ ) )
Mode e x p ( η g i ( t , g ¯ ) λ g ( t ) )
Quantile (0 < p < 1) e x p ( η g i ( t , g ¯ ) + λ g ( t ) Φ p 1 ( 0 ; 1 ) )
Variance e x p ( 2 η g i ( t , g ¯ ) + λ g ( t ) ( t ) ) · ( e x p ( λ g ( t ) ) 1 )
Table 5. Summary statistics.
Table 5. Summary statistics.
DataNumber of Plots MinMaxMeanSDNumber of PlotsMinMaxMeanSD
EstimationValidation
d (cm)231.2013.403.961.6181.808.004.151.25
h (m)231.507.903.551.1182.076.563.840.91
t (year)2353.00123.0092.2120.82853.00103.0080.6017.50
Table 6. Estimates of parameters and their standard errors (in parenthesis) for all fitted bivariate SDE models.
Table 6. Estimates of parameters and their standard errors (in parenthesis) for all fitted bivariate SDE models.
ModelParameters of SDE Models
αdβdσddαhβhσhhσdhσdσh
4 fixed
effect
4.0760
(0.0037)
0.0501
(0.0003)
0.2661
(0.0016)
3.5518
(0.0021)
0.1861
(0.0056)
0.4596
(0.0139)
0.3400
(0.0081)
--
4 mixed
effect
4.2783
(0.0023)
0.0841
(0.0009)
0.2165
(0.0024)
3.7127
(0.0011)
0.2348
(0.0640)
0.1453
(0.0396)
0.01285
(0.0258)
1.5029
(0.0097)
1.1790
(0.0076)
5 fixed
effect
0.1136
(0.0004)
0.0772
(0.0003)
0.0224
(0.0001)
0.3599
(0.0442)
0.2754
(0.0338)
0.0457
(0.0056)
0.0303
(0.0029)
--
5 mixed
effect
0.1414
(0.0003)
0.0946
(0.0002)
0.0156
(0.0001)
0.3639
(0.0004)
0.2749
(0.0003)
0.0147
(4.6 × 10−5)
0.0112
(4.4 × 10−5)
0.0274
(0.0002)
0.0750
(0.0005)
6 fixed
effect
5.2039
(0.0173)
0.0668
(0.0002)
0.0021
(5.8 × 10−6)
75.0842
(37.1990)
0.4296
(0.0492)
0.0013
(3.4 × 10−6)
0.0013
(4.0 × 10−6)
--
6 mixed
effect
7.2647
(0.0166)
0.0907
(0.0002)
0.0012
(3.3 × 10-6)
67.3215
(0.0788)
0.4159
(0.0001)
0.0004
(1.1 × 10−6)
0.0004
(1.6 × 10−6)
0.5511
(0.0035)
17.1946
(0.1104)
7 fixed
effect
3.1269
(0.0024)
0.0395
(6.3 × 10−5)
0.0023
(6.1 × 10−6)
1.0109
(0.0018)
0.0153
(4.8 × 10−5)
0.0013
(3.4 × 10−6)
0.0014
(4.1 × 10−6)
--
7 mixed
effect
3.1690
(0.0017)
0.0389
(4.6 × 10−5)
0.0012
(3.3 × 10−6)
1.0551
(0.0010)
0.0155
(2.6 × 10−5)
0.0004
(1.1 × 10−6)
0.0004
1.6 × 10−6
0.1480
(0.0010)
0.1368
(0.0009)
Table 7. Estimates of parameters and their standard errors (in parenthesis) for all fitted regression models.
Table 7. Estimates of parameters and their standard errors (in parenthesis) for all fitted regression models.
ModelParameters of Regression Models
β1β2β3Σσϕ
DiameterHeight
1 fixed
effect
1.0772
(0.0022)
1.0248
(0.0014)
-0.9898
(0.0027)
-1.4949
(0.0021)
0.6405
(0.0009)
-0.4830
(0.0013)
-
1 mixed
effect
0.9214
(0.0017)
1.1529
(0.0013)
-0.7511
(0.0021)
0.1130
(0.0007)
2.2444
(0.0012)
0.3685
(0.0004)
-0.1724
(0.0005)
0.4405
(0.0028)
2 fixed
effect
198.4369
(21.2221)
0.0157
(0.0013)
1.4884
(0.0050)
0.9846
(0.0026)
-23.7733
(1.2223)
0.0277
(0.0018)
1.0383
(0.0055)
0.4759
(0.0013)
-
2 mixed
effect
73.9614
(0.2618)
0.0499
(0.0002)
1.8394
(0.0029)
0.7466
(0.0020)
12.3094
(0.0802)
25.1334
(0.0386)
0.0037
(2.3 × 10−5)
0.5584
(0.0010)
0.1745
(0.0005)
7.3270
(0.0471)
3 fixed
effect
0.0018
(1.9 × 10−5)
−61.6512
(0.4348)
0.3074
(0.0018)
0.9843
(0.0078)
-0.0186
(0.0048)
−34.9128
(9.6510)
−0.0357
(0.0022)
0.4764
(0.0013)
-
3 mixed
effect
0.1861
(0.0006)
−1.3312
(0.0036)
0.5858
(0.0008)
0.7405
(0.0020)
0.0314
(0.0002)
0.2888
(0.0004)
−6.1644
(0.0472)
−0.7894
(0.0034)
0.1750
(0.0005)
0.0834
(0.0005)
Table 8. Statistical measures for all fitted mixed-effect parameter diameter, D, and height, H, models.
Table 8. Statistical measures for all fitted mixed-effect parameter diameter, D, and height, H, models.
Model
(Predictors)
Estimation Dataset 1Validation Dataset
B
(%B)
AB
(%AB)
RMSE
(%RMSE)
R2B
(%B)
AB
(%AB)
RMSE
(%RMSE)
R2
Diameter
Vasicek
(t)
0.0068
(−8.47)
0.8278
(23.70)
1.1114
(28.10)
0.5213−2.5 × 10−10
(−6.91)
0.8675
(22.65)
1.0674
(25.74)
0.2768
Vasicek
(t,h)
0.0025
(−4.60)
0.6462
(17.72)
0.8548
(21.62)
0.7168−1.9 × 10−10
(−4.56)
0.7243
(18.62)
0.9406
(22.69)
0.4384
Gompertz
(t)
0.0076
(−8.49)
0.8302
(23.61)
1.1207
(28.36)
0.51320.1339
(−3.43)
0.8615
(21.70)
1.0687
(25.77)
0.2750
Gompertz
(t,h)
0.0030
(−4.91)
0.6418
(17.55)
0.8567
(21.66)
0.71560.0951
(−2.30)
0.7177
(17.96)
0.9266
(22.35)
0.4549
Bertalanffy
(t)
−0.0226
(−9.38)
0.8348
(23.90)
1.1256
(28.46)
0.50900.1442
(−3.45)
0.8621
(21.76)
1.0698
(27.86)
0.2741
Bertalanffy
(t,h)
−0.0127
(−5.33)
0.6448
(17.68)
0.8587
(21.71)
0.71420.0950
(−2.28)
0.7184
(17.98)
0.9272
(22.36)
0.4543
Gamma
(t)
−0.0216
(−9.28)
0.8334
(23.85)
1.1222
(28.37)
0.51210.1338
(−3.44)
0.8623
(21.75)
1.0689
(25.78)
0.2751
Gamma
(t,h)
−0.0135
(−5.30)
0.6442
(17.67)
0.8574
(21.68)
0.71510.0950
(−2.28)
0.7185
(17.99)
0.9273
(22.37)
0.4543
Power
(h)
0.0064
(−4.50)
0.6398
(17.51)
0.8507
(21.51)
0.71950.0026
(−4.59)
0.7199
(18.57)
0.9269
(22.35)
0.4546
Richards
(h)
0.0006
(−5.07)
0.6397
(17.74)
0.8485
(21.45)
0.7209
0.0056
(−4.63)
0.7167
(18.50)
0.9269
(22.36)
0.4546
Q-exp
(h)
−0.0004
(−5.50)
0.6387
(17.97)
0.8482
(21.45)
0.72110.0147
(−4.53)
0.7183
(18.55)
0.9326
(22.49)
0.4479
Height
Vasicek
(t)
0.0026
(−2.69)
0.4049
(12.45)
0.5448
(15.35)
0.7609−4.0 × 10−11
(−1.92)
0.3912
(10.80)
0.5160
(13.44)
0.6822
Vasicek
(t,h)
0.0034
(−1.55)
0.3219
(9.75)
0.4191
(11.81)
0.85850.0106
(−1.01)
0.3363
(9.05)
0.4566
(11.89)
0.7512
Gompertz
(t)
0.0049
(−2.69)
0.4063
(12.44)
0.5476
(15.42)
0.75850.0352
(−0.99)
0.3941
(10.76)
0.5162
(13.44)
0.6821
Gompertz
(t,h)
−0.0004
(−1.54)
0.3171
(9.55)
0.4084
(11.50)
0.8657−0.0189
(−1.68)
0.3435
(9.22)
0.4657
(12.13)
0.7412
Bertalanffy
(t)
−0.0048
(−2.98)
0.4069
(12.48)
0.5484
(15.45)
0.75780.0352
(−1.01)
0.3943
(10.77)
0.5164
(13.48)
0.6819
Bertalanffy
(t,h)
−0.0060
(−1.69)
0.3174
(9.57)
0.4087
(11.51)
0.8654−0.0196
(−1.71)
0.3443
(9.25)
0.4668
(12.16)
0.7400
Gamma
(t)
−0.0031
(−2.92)
0.4065
(12.47)
0.5477
(15.43)
0.75830.0352
(−1.03)
0.3941
(10.77)
0.5162
(13.44)
0.6821
Gamma
(t,h)
−0.0050
(−1.66)
0.3175
(9.57)
0.4087
(11.51)
0.8654−0.0195
(−1.71)
0.3444
(9.25)
0.4669
(12.17)
0.7398
Power
(h)
0.0010
(−1.51)
0.3161
(9.56)
0.4066
(11.45)
0.8668
0.0118
(−0.93)
0.3422
(9.15)
0.4652
(12.11)
0.7418
Richards
(h)
−0.0002
(−1.77)
0.3197
(9.73)
0.4091
(11.52)
0.86520.0169
(−0.86)
0.3450
(9.20)
0.4709
(12.26)
0.7354
Q-exp
(h)
−0.0032
(−1.82)
0.3202
(9.76)
0.4112
(11.58)
0.86380.0156
(−0.90)
0.3442
(9.18)
0.4699
(12.24)
0.7365
1: The mean prediction bias ( B = 1 n i = 1 n ( y i y i ) ), the percentage mean prediction bias ( % B = 1 n i = 1 n y i y i y i ), the absolute mean (AB) prediction bias ( A B = 1 n i = 1 n | y i y i | ), the percentage mean absolute prediction bias ( % B = 1 n i = 1 n | y i y i y i | ), the root mean square error (RMSE) ( R M S E = 1 n 1 i = 1 n ( y i y i ) 2 ), the percentage root mean square error ( % R M S E = 1 n 1 i = 1 n ( y i y i y i ) 2 ), and the coefficient of determination ( R 2 = 1 i = 1 n ( y i y i ) 2 i = 1 n ( y i y ¯ ) 2 ). Here, n = i = 1 M n i is the total number of observations used to fit the model; M is the number of stands; ni is the number of measured trees in the ith stand; and y i , y i , and y ¯ are the measured, estimated, and average values of the dependent variable.
Table 9. Statistical measures for all fitted fixed-effect parameter diameter, D, and height, H, models.
Table 9. Statistical measures for all fitted fixed-effect parameter diameter, D, and height, H, models.
Model
(Predictors)
Estimation DatasetValidation Dataset
B
(%B)
AB
(%AB)
RMSE
(%RMSE)
R2B
(%B)
AB
(%AB)
RMSE
(%RMSE)
R2
Diameter
Vasicek
(t)
−0.0025
(−15.31)
1.1752
(33.58)
1.6318
(41.26)
00.2470
(−2.78)
1.0200
(25.08)
1.2912
(31.14)
0
Vasicek
(t,h)
−0.0006
(−5.60)
0.7237
(19.53)
0.9833
(24.86)
0.6253−0.0888
(−6.90)
0.7815
(20.54)
0.9960
(24.02)
0.3702
Gompertz
(t)
0.2425
(−2.94)
1.0273
(25.29)
1.2999
(31.35)
00.2425
(−2.94)
1.0273
(25.29)
1.2999
(31.35)
0
Gompertz
(t,h)
0.0018
(−5.78)
0.7172
(19.43)
0.9757
(24.67)
0.6310−0.0732
(−6.67)
0.7714
(20.27)
0.9771
(23.56)
0.3940
Bertalanffy
(t)
−0,1293
(−19.04)
1,2091
(35.47)
1.6420
(41.52)
00.1323
(−5.87)
1.0373
(26.22)
1.3092
(31.58)
0
Bertalanffy
(t,h)
−0.0069
(−6.24)
0.7185
(19.55)
0.9785
(24.74)
0.6290−0.0662
(−6.61)
0.7710
(20.24)
0.9751
(23.52)
0.3964
Gamma
(t)
−0,0982
(−18.76)
1.2670
(37.10)
1.7045
(43.10)
0−0.0485
(−10.95)
1.1219
(29.41)
1.3827
(33.25)
0
Gamma
(t,h)
−0.0207
(−6.49)
0.7323
(19.97)
0.9885
(24.99)
0.6213−0.1528
(−8.83)
0.7935
(21.17)
1.0064
(24.27)
0.3571
Power
(h)
0.0043
(−5.57)
0.7357
(19.78)
0.9977
(25.23)
0.6142−0.1330
(−7.97)
0.7862
(20.78)
1.0147
(24.47)
0.3465
Richards
(h)
−0.0017
(−6.21)
0.7423
(20.11)
0.9932
(25.11)
0.6177−0.0961
(−7.15)
0.7847
(20.58)
1.0127
(24.42)
0.3490
Q-exp
(h)
−0.0024
(−6.26)
0.7409
(20.15)
0.9930
(25.11)
0.6178−0.0926
(−7.08)
0.7843
(20.56)
1.0125
(24.42)
0.3493
Height
Vasicek
(t)
−0.0016
(−8.73)
0.8132
(23.77)
1.1142
(31.39)
00.2881
(2.68)
0.6839
(16.67)
0.9155
(23.84)
0
Vasicek
(t,h)
−0.0010
(−3.33)
0.4980
(14.42)
0.6708
(18.89)
0.63760.1526
(1.69)
0.5297
(13.27)
0.6963
(18.13)
0.4214
Gompertz
(t)
0.0062
(−8.49)
0.8114
(23.67)
1.1142
(31.38)
00.2957
(2.89)
0.6848
(16.66)
0.9154
(23.84)
0
Gompertz
(t,h)
0.0045
(−3.28)
0.5006
(14.50)
0.6711
(18.90)
0.63720.1340
(1.15)
0.5311
(13.35)
0.6923
(18.03)
0.4281
Bertalanffy
(t)
−0.0921
(−11.50)
0.8369
(25.05)
1.1142
(31.38)
00.1975
(0.19)
0.6787
(16.95)
0.9154
(23.84)
0
Bertalanffy
(t,h)
−0.0385
(−4.34)
0.5029
(14.75)
0.6686
(18.83)
0.63990.0886
(0.04)
0.5369
(13.64)
0.6963
(18.13)
0.4214
Gamma
(t)
−0.0389
(−9.95)
0.8328
(24.84)
1.1207
(31.57)
00.1214
(−1.95)
0.6948
(17.74)
0.9296
(24.21)
0
Gamma
(t,h)
−0.0144
(−3.57)
0.4940
(14.46)
0.6547
(18.44)
0.65470.0960
(0.37)
0.5331
(13.58)
0.6901
(17.97)
0.4316
Power
(h)
0.0037
(−3.07)
0.5162
(14.87)
0.6970
(19.63)
0.60870.1598
(1.81)
0.5534
(13.82)
0.7325
(19.08)
0.3596
Richards
(h)
0.0014
(−3.31)
0.5095
(14.66)
0.6905
(19.45)
0.61590.1783
(2.26)
0.5439
(13.51)
0.7224
(18.81)
0.3772
Q-exp
(h)
−0.0006
(−3.44)
0.5096
(14.67)
0.6908
(19.46)
061560.1826
(2.33)
0.5395
(13.38)
0.7196
(18.74)
0.3820
Table 10. Means and variances of stationary Gompertz-type marginal and conditional distributions.
Table 10. Means and variances of stationary Gompertz-type marginal and conditional distributions.
TypeMeanVariance
Marginal μ g = α g + φ g i β g v g g = σ g g 2 β g
Conditional η g i ( g ¯ ) = μ g i + v g g ¯ v g ¯ g ¯ ( l n ( g ¯ ) μ g ¯ i ) λ g = v g g ( v g g ¯ ) 2 v g ¯ g ¯
Table 11. Gompertz-type stationary marginal and conditional mean, median, mode, p-quantile (0 < p < 1), and variance trajectories.
Table 11. Gompertz-type stationary marginal and conditional mean, median, mode, p-quantile (0 < p < 1), and variance trajectories.
TypeTrajectory TypeEquation
MarginalMean e x p ( μ g i + 1 2 v g g )
Median e x p ( μ g i )
Mode e x p ( μ g i v g g )
Quantile (0 < p < 1) e x p ( μ g i + v g g Φ p 1 ( 0 ; 1 ) ) 1
Variance e x p ( 2 μ g i + v g g ) · ( e x p ( v g g ) 1 )
ConditionalMean e x p ( η g i ( g ¯ ) + 1 2 λ g )
Median e x p ( η g i ( g ¯ ) )
Mode e x p ( η g i ( g ¯ ) λ g )
Quantile (0 < p <1) e x p ( η g i ( g ¯ ) + λ g Φ p 1 ( 0 ; 1 ) )
Variance e x p ( 2 η g i ( g ¯ ) + λ g ) · ( e x p ( λ g ) 1 )
1: Φ p 1 ( · ; · ) is the inverse of the normal distribution function.
Table 12. Statistical measures for all fitted basal, BA, area (m2) models.
Table 12. Statistical measures for all fitted basal, BA, area (m2) models.
PdfEstimation DatasetValidation Dataset
B
(%B)
AB
(%AB)
RMSE (%RMSE)R2B
(%B)
AB
(%AB)
RMSE
(%RMSE)
R2
Mixed-effect scenario
L N 1 ( μ d i ( t ) ; v d d ( t ) ) 0.0078
(0.23)
0.0187
(8.69)
0.0300
(15.94)
0.96020.0066
(4.08)
0.0080
(5.08)
0.0079
(4.84)
0.9886
L N 1 ( η d i ( t , h ) ; λ d ( t ) ) 0.0040
(0.27)
0.104
(5.23)
0.0151
(8.02)
0.98990.0062
(3.72)
0.0072
(4.45)
0.0075
(4.56)
0.9899
L N 1 ( μ d i ; v d d ) 0.0038
(−1.67)
0.0171
(8.51)
0.0286
(15.19)
0.96390.0066
(4.08)
0.0080
(5.09)
0.0079
(4.84)
0.9886
L N 1 ( η d i ( h ) ; λ d ) −0.0001
(−1.69)
0.0090
(4.99)
0.0135
(7.19)
0.99190.0062
(3.72)
0.0072
(4.45)
0.0075
(4.56)
0.9899
Fixed-effect scenario
L N 1 ( η d i ( h ) ; λ d ) 0.0119
(−0.37)
0.0394
(18.67)
0.0684
(36.35)
0.7930−0.0066
(−3.00)
0.0224
(14.27)
0.0260
(15.80)
0.8787
Table 13. Statistical indexes for the mean tree volume, V ¯ , (m3) model, using the mixed-effect scenario.
Table 13. Statistical indexes for the mean tree volume, V ¯ , (m3) model, using the mixed-effect scenario.
Estimation DatasetValidation Dataset
B
(%B)
AB
(%AB)
RMSE
(%RMSE)
R2B
(%B)
AB
(%AB)
RMSE
(%RMSE)
R2
0.0003
(−0.35)
0.0009
(10.27)
0.0017
(24.64)
0.9500−3.5 ×10−5
(1.41)
0.0003
(5.07)
0.0005
(9.34)
0.9765

Share and Cite

MDPI and ACS Style

Narmontas, M.; Rupšys, P.; Petrauskas, E. Construction of Reducible Stochastic Differential Equation Systems for Tree Height–Diameter Connections. Mathematics 2020, 8, 1363. https://doi.org/10.3390/math8081363

AMA Style

Narmontas M, Rupšys P, Petrauskas E. Construction of Reducible Stochastic Differential Equation Systems for Tree Height–Diameter Connections. Mathematics. 2020; 8(8):1363. https://doi.org/10.3390/math8081363

Chicago/Turabian Style

Narmontas, Martynas, Petras Rupšys, and Edmundas Petrauskas. 2020. "Construction of Reducible Stochastic Differential Equation Systems for Tree Height–Diameter Connections" Mathematics 8, no. 8: 1363. https://doi.org/10.3390/math8081363

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop