Prediction System for Diameter Distribution and Wood Production of Eucalyptus

The aim of the present study was to propose a prediction system to estimate the diameter distribution and wood production in unthinned eucalyptus clone plantations. Data was obtained from permanent sample units, with ages from 28 to 78 months. The Weibull distribution was used to estimate the frequency of the number of trees per diameter class. Models were used to relate the Weibull distribution coefficients with forest stand attributes. The forest stand variables most correlated with the Weibull distribution parameters were the minimum, medium, maximum and quadratic diameters and the dominant height. The projected frequency of tree numbers and the production by diameter class obtained by the system of equations did not differ statistically from the observed values evaluated by the t-test at 95% probability.


INTRODUCTION
Diameter distribution is an indicator of stock growth structure and allows us to draw effective conclusions regarding forest structure (Loetsch et al., 1973), which constitutes basic knowledge for forest management and helps to monitor forest resource dynamics (Hosokawa et al., 2008). According to Binoti et al. (2014), the description of the frequency observations by diameter class is performed using probability density functions (pdf).
Among the most common probability density functions in forest management, the Weibull stands out due to its relative simplicity and flexibility when representing different distribution trends (Bailey & Dell, 1973;Duan et al., 2013). The main advantage of using it is that its parameters are easily correlated with stand variables, which makes it possible to calculate them for future ages, in addition to being useful in studies on forest growth and production (Nokoe & Okojie, 1984;Soares et al., 2010).
A pdf is the main tool to model forest stand structure in uneven-aged and heterogeneous forests (Nascimento et al., 2012). There are many methods that can be used to estimate pdf parameters (Piqué-Nicolau et al., 2011), but they depend on the required accuracy and, according to Binoti et al. (2011), the fitted distribution. The maximum likelihood estimation method is the most suitable to adjust the Weibull distribution since the computational procedures optimize the parameter estimation process (Araújo et al., 2010).
Diameter distribution estimates can be realized using prediction or projection models (Leite et al., 2013). In the first case, the parameters of the selected statistical distribution are correlated with stand parameters (Clutter et al., 1983). This method is known as the parameter prediction method (PPM) (Kangas & Maltamo, 2000). In the second case, a method known as parameter recovery method (PRM) (Kangas & Maltamo, 2000), the statistical distribution parameters are estimated based on the same parameters obtained at a previous age .
Regardless of the category, these models were mostly used to predict the diametric distribution in stands with species of the genus Pinus, managed for long cycles and for multiple-use woods (assortments). By contrast, these models have rarely been applied in eucalyptus stands. Thus, the present study was developed with the hypothesis that the proposed prediction system might generate accurate estimates for a certain number of trees and volume per hectare for clonal Eucalyptus sp. stands. Therefore, the objective of the present study was to propose a prediction system to estimate the diameter distribution and wood production in unthinned Eucalyptus sp. stands.

MATERIAL AND METHODS
The database was taken from permanent sample units of a continuous forest inventory from 2007 to 2010, in stands of hybrid Eucalyptus grandis W. Hill ex Maiden x Eucalyptus urophylla S. T. Blake. The stands are located in the northeastern region of the state of Bahia, Brazil, with an initial average density of 1,111 trees per hectare. Circular plots were used, with 12.25 m radius, randomly distributed, with an average sampling intensity of one plot per 12.5 hectares. The diameters at breast height with bark of all trees (d at 1.3 m from the ground), the total height (TH) of the first 15 trees and the total height of the five dominant ones (hdom) were measured in the plot. From seventy-one sample units, 48 were used for fitting and 23 for validation of the models. The ages ranged from 28 to 78 months. The synthesis of data used is presented in Table 1.
The hypsometric model was adjusted for data pairs of total height and diameter at breast height collected in the plots, as well as to data pairs of scaled trees using the Smalian's method to fit the volume equations. In total, 48 trees were scaled. The general equation (Equation 1) of total height per tree, as used by Ribeiro et al. (2010), resulted in an adjusted coefficient of determination -R 2 adjusted = 0.94 and standard error of the estimate in percentage -S yx = 4.9%.
The Richards biological model, adjusted by Rodríguez-Carrillo et al. (2015), was used to study the dominant height and age relationship (Syx = 13.5%). Based on this model, the expression for a definition of site classes (S) per plot was obtained (Equation 4), adopting 60 months as the age index (Ii). Three productivity site classes were defined: high (class 1: hdom > 25.9), medium (class 2: 21.3 < hdom ≤ 25.9 m) and low (class 3: hdom ≤ 21.3 m). Where: I = age (months); dc = tree diameter of the central diameter class; h, hdom, S, Ii, hc, dg and d = as previously defined, ln = Napierian logarithm; e = exponential. For all these equations, the coefficients were significant (p-value < 0.05).
To project the diameter distribution, the parameter prediction method (PPM) was used (parameter estimation). In this method, the coefficients of the selected probability density function (pdf) were correlated with parameters of the forest stand, generating equations to estimate them at a future age. The pdf used was the three-parameter Weibull distribution (Wendling et al., 2011) described in Equation 5, where "a" is the location parameter, "b" the scale parameter, "c" the form parameter, and xi is the diameter at breast height (d), in which: a ≤ xi ≤ ∞, a ≥ 0, b > 0 and c ( ) The location parameter (a) was independently obtained. The solver adds Microsoft  Excel  was used to optimize the parameter "a", whose selected value was that which minimized the sum of squares of the difference between the observed and estimated frequency. The amplitude of the diameter classes used was 2 cm. Kolmogorov-Smirnov's test (KS test) was applied to analyze the adherence of diameter distributions in each permanent sample plot at all ages (Téo et al., 2012).
Survival estimates in forest stands were evaluated using data from the permanent sample plots. For the adjustment of the models, failures in the stands were disregarded. The models used (6, 7, 8 and 9) are presented in Table 2 and can be found in Retslaff et al. (2012).
For the selection of the variables to be used when modeling the coefficients of the Weibull distribution, a preliminary analysis was performed on the correlation matrix. The variables most correlated with the coefficients were subjected to the stepwise process, characterizing the finally adjusted models (12, 14 and 15). Furthermore, models from the literature were tested (10, 11 and 13) ( Table 3).
Equations for projection of stand attributes were selected, which are the independent variables of the models for parameters "b" and "c" of the Weibull distribution. The equations were defined from linear relationships of variables based on the simple linear correlation matrix (Equations 16,17,18,20,21,22,24,25,26,28 and 29), in addition to the application of some models found in the literature (Equations 19, 23, 27, 30 and 31) ( Table 4).
For evaluation and selection of models, the adjusted coefficient of determination (R 2 adjusted ) and the standard error of the estimate in percentage (S yx %) were used Table 2. Survival models tested for Eucalyptus sp. stands in northeastern Bahia, Brazil.

Equation Model
Model Formulation Where: N 2 and N 1 = number of trees at future and current ages, respectively (N.sample-unit -1 ); I 2 and I 1 = future and current age, respectively (months); S 2 = site index at future age; β i = coefficients to be estimated; ln = Napierian logarithm. Where: b 2 = Weibull distribution scale parameter at future age; c 2 = Weibull distribution form parameter at future age; d 2 = average diameter at future age (cm); dg 2 = quadratic diameter at future age (cm); d min2 = minimum diameter at future age (cm); βi = coefficients to be estimated; ln = Napierian logarithm.

5/12
Prediction System for Diameter Distribution… Floresta e Ambiente 2018; 25(3): e20160548 (Tonini & Borges, 2015). The parameter estimates had their significance evaluated by p-value (p-value < 0.05). In cases where the coefficient was not significant, the variable associated with it was disregarded and the model readjusted. The selected equation system was applied for the validation of the database. Weibull distribution coefficients and the diameter distribution were used to evaluate the effectiveness when estimating the forest stand attributes for different ages and sites. Distributions were evaluated by KS test (95% probability). According to Lei (2008), it is important that the probability density function estimates are compared to the observed values, because the estimated parameters play an important role in the stand yield model, employing the parameter prediction method (PPM).
Using the hypsometric and volumetric equations, volumes by diameter class and unit area were estimated. The statistical comparison between the projections and the corresponding values from the forest inventory was done using the t-test (95% probability).

RESULTS AND DISCUSSION
The results of the adjustments of the best survival models, by site class, are presented in Table 5. In general, the best statistical evaluations were obtained for sites 1 and 2, where the Clutter-Jones model presented superior performance. For site 3, the number of tree estimations at a future age was obtained with the Pienaar and Shiver model. In all models, trends in the estimates were not detected.  Where: I 2 and I 1 = future and current ages, respectively (months); d min2 and d min1 = minimum diameter at future and current age, respectively (cm); dg 2 and dg 1 = quadratic diameter at future and current age, respectively (cm); d max2 and d max1 = maximum diameter at future and current age, respectively (cm); S 2 = site index at future age; d 2 and d 1 = average diameter at future and current age, respectively (cm); β i = coefficients to be estimated. The linear correlation matrix between forest stand variables and Weibull distribution parameters are presented in Table 6. The highest correlations were observed in the variables d, d max , dg and hdom.
The Weibull distribution location parameter (a) showed a weak correlation with the forest stand variables. The minimum diameter (d min ) showed a negative correlation, being the only statistically significant case for this parameter. This result was also found by Bailey & Dell (1973), who reported that the Weibull distribution location parameter presents a significant relationship with the smallest diameter of the forest stand. The scale (b) and form (c) parameters showed a statistically significant correlation with all forest stand variables ranging from low to moderate values, except between density (N) and the form parameter (c).
The smaller value of the sum of squares of differences between the observed and estimated frequencies was obtained when the location parameter (a) was equal to or very close to zero. Therefore, this parameter was considered equal to zero, resulting in the two-parameter Weibull distribution. This distribution was used by Gove (2003), Zasada (2013) and Sanquetta et al. (2014). Figura (2010) tested different percentages of minimum diameter to represent the location parameter of Weibull distribution in stands of Eucalyptus grandis, in southern Brazil. The best results were obtained using a percentage of minimum diameter equal to zero, similar to the results obtained in this study.
The other Weibull distribution parameters were obtained using selected equations and for the future scale parameter (b 2 ), the adjustments of the accuracy statistics were satisfactory, reaching R 2 adjusted values higher than 0.95 for all adjusted equations. The S yx values were very low, less than 1% for Equation 12 (Stepwise), which was selected to estimate the scale parameter for the three site classes ( Table 7). Estimates of form parameters (c 2 ) at future ages resulted in R 2 adjusted and S yx % values lower than those obtained for the scale parameter, however, they remained adequate results. The models selected to estimate c 2 are the Equation 13 (Retslaff, 2010), Equation 14 and Equation 15 (Stepwise), for site classes 1, 2 and 3, respectively. The estimates  showed a good performance for all models, resulting in unbiased estimates.
The models were selected to estimate the parameter of the Weibull distribution, at a future age, applying the variables: d min , d, d max and dg. These express the forest stand attributes to be estimated and to constitute the projection equation system. Table 8 presents the coefficients and statistics for evaluation of selected equations to estimate these stand attributes.
Estimates of forest stand attributes at a future age were unbiased: d 2 , dg 2, and d max2 . However, a slight trend was observed, moving from an overestimation to an underestimation for the d min2 , even though it did not impose restrictions on the equations selected.
When considering d 2 , Equation 31 (Retslaff, 2010) showed better performance for its estimation in all site classes, in which R 2 adjusted was equal to or greater than 0.95, and the S yx % was close to zero. These results are due to the high correlation between d and dg variables, observed in the previously presented correlation matrix. For dg 2 , Equation 24 (Stepwise) showed the best performance for its estimation in site class 1, while for site classes 2 and 3, Equation 27 (Retslaff, 2010 -modified) was superior.
The estimate of d max2 presented a high R 2 adjusted and low S yx %. The Equation 23 (Retslaff, 2010) resulted in better statistics to represent site class 3. For site class 2, the best result was obtained with the Equation 22 (Stepwise). For the estimation of d min2 , Where: * = significant at 95% probability; ns = not significant at 95% probability; -: Coefficient not present in the model; R² adjustment = adjusted coefficient of determination; S yx % = relative standard error of estimate; β i = adjusted coefficients. Equation 19 (Retslaff, 2010) also resulted in better performance for site classes 2 and 3.
In general, the results were considered satisfactory. Therefore, the models can be used for future estimates of forest stand attributes with precision and, consequently, to estimate Weibull distribution parameters at a future age.
Once the equations were selected to estimate the Weibull distribution parameters, in addition to the forest stand attributes at a future age, they were also applied to data of 23 surplus plots for validation purposes. Table 9 presents the statistics for the selected models to estimate the forest stand attributes, in which the R 2 adjusted ranged from 0.81 to 0.99, and S yx from 0.96 to 5.95%. The results were obtained for the equation to estimate d 2 , indicating a high correlation of this variable with dg 2 , depending on the model selected to estimate this variable.
The Weibull distribution parameters were estimated from selected equations and estimated frequencies for the number of trees in each sample plot were obtained. The observed and estimated frequencies were subjected to evaluation by Kolmogorov-Smirnov's test (95%). Of the 23 sample plots selected for validation, 20 showed statistic "D" less than the respective table value, which indicates the adherence of the distributions. The remainder did not present adherence. In one sample plot, positive asymmetry occurred and in the others, a large number of trees appeared in the initial diameter classes, which provoked difficulties when adjusting the Weibull distribution in such circumstances.
The number of trees per hectare and diameter class was obtained for each site class and year of measurement to evaluate the accuracy of these estimates (Figure 1).  Where: N 2 = number of trees at future age; d 2 = diameter at breast height at future age (cm); dg 2 = quadratic diameter at future age (cm); d max2 = maximum diameter at future age (cm); d min2 = minimum diameter at future age (cm); R² adjustment = adjusted coefficient of determination; S yx % = relative standard error of estimate.

9/12
Prediction System for Diameter Distribution… Floresta e Ambiente 2018; 25(3): e20160548 Palahí et al. (2006) obtained similar results in modeling the diameter distribution in Pinus spp. stands using the Weibull distribution. There, the parameter optimization method, when minimizing the differences between observed and estimated accumulated frequencies, provided the best estimates in relation to regression methods.
The hypsometric and volumetric equations were applied to the central value of each diameter class, in accordance with the frequency estimated by Weibull distribution. The product of the number of trees in each class by the volume of the central tree made it possible to obtain the volume per class. The sum of the volume of each class resulted in the volume of the sample plot and, by extrapolation, the results per hectare were obtained. The results of observed and estimated production values are presented in Table 10.
The production estimates followed the same trend in the various site classes. Site 1 was the only one that presented production values in the upper diameter classes. For sites with high and medium yield, there was a tendency to underestimate the production values,    Where: Vobs and Vest = observed and estimated volumetric production, respectively (m 3 .ha -1 ). while in lower productivity sites, the estimated values were above the respective values observed for all years of measurement. For all cases, there has been progress in the diameter class production with increasing age, as expected.
The "F" test was carried out to verify the homogeneity of the variances between the total observed and estimated outputs, by site class and year of measurement. The results for this test were: 1.86 (p-value = 0.31), 1.97 (p-value = 0.30) and 1.55 (p-value = 0.36), for site 1, 2 and 3, respectively, indicating homoscedasticity. Consequently, the t-test was applied for two independent samples, assuming equal variances. The "t" value was 0.90 (p-value = 0.40), 0.26 (p-value = 0.81) and 0.69 (p-value = 0.52), for sites 1, 2 and 3, respectively, indicating that the total outputs observed and estimated by site class did not differ statistically at 95% probability.
All systems provided consistent results, however, differences in the yield estimates at some ages were observed. This was also observed by Abreu et al. (2002) for Eucalyptus grandis, in the state of São Paulo, Brazil, and Nogueira et al. (2005) for thinned Eucalyptus sp. stands, in northeastern Bahia, Brazil.
Based on these results, the system was considered suitable for projection of diameter distribution and production in Eucalyptus sp. stands in northeastern Bahia, Brazil. The synthesis of the equations system stratified by site classes is as follows:

CONCLUSIONS
According to the results, the proposed system of equations is efficient to project the frequencies of the number of trees as well as the wood production by diameter class for Eucalyptus sp. stands, and may be used as a tool to aid in planning forest production.