General Height-Diameter Equation Depending on the Stand Variables, for Eucalyptus benthamii

1Universidade do Oeste de Santa Catarina, Xanxere, SC, Brasil Abstract The aim of this study was to develop a general applicability equation to represent the hypsometric relationship of Eucalyptus benthamii Maiden et Cambage stands in Irani, Santa Catarina. The data came from a continuous forest inventory, with a simple random sampling design and the Bitterlich method. The functions to represent the hypsometric relationship were adjusted with their parameters linear function of stands variables, based on the Gauss-Newton algorithm. The models were evaluated according to the criteria of adjustment and precision, accuracy and validation. In general, the Harrison model, adjusted with the coefficient “c” linear function of age, dominant diameter and basal area, presented the best statistical performance to represent the hypsometric relationship of Eucalyptus benthamii.


INTRODUCTION AND OBJECTIVES
In Brazil, the first Eucalyptus benthamii Maiden et. Cambage stands were planted in 1988 by EMBRAPA Florestas, in the municipality of Colombo, state of Paraná. The forest species aroused interest soon after its introduction in Brazil, due to its tolerance to low temperatures, good development, high growth rate and healthy appearance of the plants. (Alves et al., 2011).
Thus, in recent years, a considerable expansion in the areas of E. benthamii production has been verified in southern Brazil, especially in regions of frost occurrence. However, in spite of all the results presented by the species and the massive increase of plantation area, the main work carried out, until now, are basically related to its potential use in the pulp and paper industry (Müller et al., 2014).
Nevertheless, it is important to know the allometric relations of the species, among them the height-diameter relationship, which can provide technical directions to increase the productive potential of the species. It is extremely important to study the dendrometric variables of a forest, since they are used as the basis for appropriate silvicultural treatments in a forest stand, supporting forest management and planning.
The height measurement is a costly activity, generally making it impossible to obtain the height of all trees within a forest inventory plot. Based on the aforementioned, it is common to measure the diameter of all trees in the plot and the height of some, during forest inventory. Then, through height-diameter observations, mathematical functions must be adjusted and tested to estimate the heights of the non-measured trees. The use of these height-diameter equations reduces the cost of forest inventories (Ribeiro et al., 2010). Azevedo et al. (2011), reported that the use of height-diameter equations is of great practical significance and an important aspect to be considered in forest inventories.
The height-diameter relationship is susceptible to several factors, such as site quality (Retslaff et al., 2015;Vendruscolo et al., 2015), age (Machado et al., 2011;Paulo et al., 2011;Soares & Tomé, 2002;Téo et al., 2017), stand density (Araújo et al., 2012;Silva et al., 2016;Téo et al., 2017), crown size (Machado et al., 2015), species Donadoni et al., 2010;Hess et al., 2014;Sanquetta et al., 2013) and the sociological position of the trees (Costa et al., 2014;Leduc & Goelz, 2009;Martins et al., 2016). Retslaff et al. (2015) reported that the use of the variables dominant height and site index produced height-diameter equations with more homogeneous and narrow residues distribution. The effect of stand density on the height-diameter relationship was evaluated by Silva et al. (2016), who found an accuracy reduction of the height estimates of Tectona grandis L. F., for wider spacings; similar results were found by Bartoszeck et al. (2004). Concerning the effect of the variables sociological position and crown length, Costa et al. (2014) and Machado et al. (2015) found height-diameter equations with better statistical performance fitted into sociological position and crown length classes, respectively.
In this context, it becomes imperative to select the most appropriate equation to estimate height in each forest stand, considering the large number of variables that effect the height-diameter relationship (Oliveira et al., 2011).
The height-diameter equations may have local or general applicability. Commonly, the local height-diameter equations are only dependent on tree diameter at breast height (d) and can only be applied to the stand where the data were gathered. On the other hand, the general or regional height-diameter equations are dependent on diameter at breast height and stand variables, such as age (t), dominant height (hdom) and stand density, developed for general application to a species or region (Soares & Tomé, 2002).
In conducting the development of height-diameter equations of a Eucalyptus globulus Labill stand in Portugal, Tomé et al. (2007) recommended an equation with its parameters expressed as linear combination of dominant height (hdom), dominant diameter (ddom) and quadratic mean diameter (dg) of the forest stand as independent variables.
In Brazil, it is common to use empirical models to represent height-diameter relationship (Oliveira et al., 2011;Rufino et al., 2010). However, it is more appropriate to use theoretical models, which have an implicit hypothesis associated with the characteristics of the height-diameter relationship and considering the causes and effects of different variables on the height-diameter relationship (Téo et al., 2017).
Thus, this research aimed to develop general applicability height-diameter equations, with their parameters expressed as function of forest stand variables for Eucalyptus benthamii Maiden et Cambage, in Irani, state of Santa Catarina, in order to indicate the one that best describes the height-diameter relationship.

MATERIALS AND METHODS
The present study was undertaken in a forest stand of Eucalyptus benthamii, located in the municipality of Irani, state of Santa Catarina, Brazil (Figure 1). According to Alvares et al. (2013), the study area presents an altitude of 887 m a.s.l., Cfb climate type based on the Köppen classification, that is, humid subtropical zone, oceanic climate, without a dry season and with temperate summers. The average annual temperature of the region is 16,6 °C, the average temperature of the hottest month is 20,7 °C and the coldest month is 12,2 °C. The average annual relative humidity ranges from 76 and 78%, the total annual precipitation is 1836 mm, evenly distributed all year long. The main soil types in the region are Haplic Cambisols, Humic Cambisols, Litholic Neosol, Red Nitisol and Haplic Nitisol (EMBRAPA, 2011).
The forest stand of Eucalyptus benthamii has an area of 13,292 ha, initial spacing of 2,5 x 2,5 m, or 1600 plants per hectare. The forest stand was sampled by a continuous forest inventory with two phases, the first in June 2015 at 7,5 years of age and the second in June 2017 at 9,5 years of age. Twelve sampling units were measured in each occasion, distributed according to the unrestricted random design and using the Bitterlih sampling method. The sampling method was operationalized using the Criterion RD 1000, with the basal-area factor equal to 1.
For each Eucalyptus benthamii tree counted by Bitterlih sampling method, were measures the following variables: i) circumference at breast height (c), directly with a metric tape, which were subsequently transformed in diameter and ii) total height (h) in meters, measured with TruPulse 200B hypsometer, at least ten trees per sampling unit.
After data collection, four models were adjusted and tested to represent the hypsometric relationship of Eucalyptus benthamii (Table 1).

-10
Prodan modified by Tomé (1988)  The candidate functions were adjusted with their parameters expressed as a linear function of different categories of stand variables. When the Harrison et al. (1986) model was adjusted with its parameter "a" expressed as a linear function of stand variables, it was denominated Harrison_a, when the adjustment was proceeded with its parameter "b" expressed as linear function of stand variables, it was denominated Harrison_b, and finally, Harrison_c (Equation 1). The models by Michailoff, Prodan and Stoffels & Van Soest have only one parameter "a" expressed as linear combination of stand variables, so they were denominated Michailoff_a, Prodan_a and Stoffels & Van Soest_a.
Firstly, the stand variables which the parameters of the height-diameter functions were expressed as linear combination (Equation 1), were organized into different categories: I) age; II) sample plot; III) stand density.
The first category was represented by age (t) in years. The second category was represented by: mean diameter (d mean ), in cm; dominant diameter (d dom ), in cm; coefficient of variation of the diameters (CV d ), in %; quadratic mean diameter (d g ), in cm, calculated according to Machado & Figueiredo Filho (2014); maximum diameter of the sample plot (d max ), in cm; and dominant height (h dom ), in m. The variables of the category stand density were: basal area (G), in m 2 .ha -1 ; number of trees (N), trees.ha -1 ; the transformations of the number of trees N-1 and 100*N-1; and Wilson Factor (WF) (Equation 2) (Wilson, 1951).
For each category of independent variables was applied the multiple linear regression, where total height (h) of the Eucalyptus benthamii trees was the dependent variable and the stand variables of each category was the regressor variables. The selection of the best regressor variable within each category was made with PROC REG procedure in the software program SAS ® University Edition, based on the rsquare variable selection method.
The linear combination of stand variables, presented in Equation 1 for each parameter, was based on linear multiple regression with the best stand variable of each category as regressor variables and again, the total height (h) as independent variable. For the final selection of the stand variables, we analyzed the significance of parameters (Student's value) and the variance inflation factor (VIF), to avoid collinearity, according to Myers (1986).
The candidate functions (Table 1) to represent the heightdiameter relationship were adjusted for 289 values of total height (h) and diameter at breast height. The parameter estimation of the functions was made with PROC NLIN procedure in the software program SAS ® University Edition, based on the Gauss-Newton iterative method. The studies of Harrison et al. (1986) and Téo et al. (2017) provided the basis to obtain the initial values of the parameters of the heightdiameter functions. The final estimates of the parameters of the height-diameter equations should exclude zero in its 95% confidence band, indicating that there are only nonzero values for the parameters and then, they are always significant (Téo et al., 2017).
In the initial selection phase, the criteria for comparing the performance of the height-diameter equations were:  th HAT matrix element.
The graphical analysis of the studentized residuals was performed on the estimated dependent variable Y i  , which was aimed at the uniform distribution of residuals with absence of patterns. The verification of the heteroscedasticity was made by visual analysis of the distribution of studentized residuals, seeking to identify the funnel effect indicates that as the estimated dependent variable Y i  gets large, the deviations of the residuals from zero become greater (Myers, 1986;Draper & Smith, 1998;Montgomery et al., 2006). Once heteroscedasticity of the distribution of studentized residuals of the equation resulting from the adjustment by the ordinary least squares method was detected, the weighted least squares method was adjusted. The weights tested included various transformations and combinations of stand density variables (basal area, in m².ha -1 ; number of trees, in trees. ha -1 ) and sample plot variables (mean diameter, quadratic mean diameter, maximum diameter, dominant diameter, in cm; coefficient of variation of the diameters, in %), based on the protocol described by Parresol (1993), which assumes that the variance of the residuals is an exponential function of multiple explanatory variables of the model, as well as their transformations and combinations. To verify whether the studentized residuals had normal distribution or not, a graphical representation of these increasingly ordered values against the theoretical quantiles of the normal distribution was constructed. Once the non-normal distribution of the studentized residuals was verified, the weighted least squares method was adjusted, with weights assigned to the studentized residuals exceeding the interval ± 2, based on the method described by Huber (1964) and recommended by Myers (1986).
In addition to the statistics calculated by the residuals (R² ad , syx%, MD e MAD), the height-diameter equations were submitted to validation using PRESS residuals (e i , -i ) (Equation 8). The PRESS residual of the observation i is calculated by subtracting from the observed dependent variable (Yi) the value of the dependent variable estimated by the equation adjusted for the dataset without the observation i (Y i i  , − ). Consider a set of data in which the first observation was left aside from the sample, and the remaining "n -1" observations were used to estimate the coefficients of the height-diameter model and obtain Y i i  , − for the first observation. The first observation is then replaced, and the second observation withheld with the coefficients estimated again and obtain Y i i  , − for the second observation. This procedure continued until all observations were removed one by one, and thus the model was adjusted "n" times (Myers, 1986).

RESULTS AND DISCUSSION
All the models tested to represent the hypsometric relationship of Eucalyptus benthamii presented significant estimates of the parameters (α = 0,05), in the way they are in Table 1.
The Harrison et al. (1986) model adjusted with its parameters expressed as a linear function of stand variables showed the best statistical performance to represent the height-diameter relationship of Eucalyptus benthamii (Table 2). Among the different versions of the Harrison et al. (1986) model, Harrison_b presented the best adjustment (R2ad), the smaller error (syx%) and the value of MAD closer to zero, while Harrison_a and Harrison_c presented equal values of R 2 ad and syx%. The Harrison et al. (1986) with its parameter "c" expressed as a linear function of stand variables presented the bias (MD) closer to zero.
Similar to the results found in this study, Téo et al. (2017) reported the best statistical performance for Harrison et al. (1986) model adjusted with its parameters expressed as a linear function of stand variables, to represent the heightdiameter relationship of Pinus taeda L. in the Midwest region of Santa Catarina State, Brazil. The different versions of the Harrison et al. (1986) model, adjusted by Téo et al. (2017), also presented negative bias (MD).
The models by Michailoff, Prodan and Stoffels & Van Soest, modified by Tomé (1988), showed similar statistical performance when compared to each other. However, when they were compared with the Harrison et al. (1986) model, they presented poorer adjustment (R 2 ad ), greater error (syx%) and bias (MD) ( Table 2). Where: R 2 ad = adjusted coefficient of determination; syx% = standard error of estimate; MD = mean of differences; MAD = mean of absolute differences.
The Harrison et al. (1986) model adjusted with its parameters expressed as a linear function of stand variables, by ordinary least squares, produced distribution of studentized residuals without heteroscedasticity and without evidences of non-normality, then weighted least squares regression with Parresol (1993) and Huber (1964) methods were not necessary (Figure 2). Téo et al. (2017) reported studentized residuals without heteroscedasticity, but with evidences of non-normality of the residuals, especially with large residuals at both extremes of the distribution of errors for the Harrison et al. (1986) model.
The models by Michailoff, Prodan and Stoffels & Van Soest adjusted with its parameters expressed as a linear function of stand variables, by ordinary least squares, produced distribution of studentized residuals without heteroscedasticity, however with characteristics of non-normality of the residuals, mainly at the negative values of the distribution of errors. The Michailoff, Prodan and Stoffels & Van Soest equations tended to underestimate the predicted values of total height smaller than 20 m. (Figure 3).
Similarly to what happened with the selection statistics, the Harrison et al. (1986) model adjusted with its parameters expressed as a linear function of stand variables, showed the best performance for validation statistics to height-diameter of Eucalyptus benthamii (Table 3). Among the different versions of the Harrison et al. (1986)    According to the selection statistics (Table 2), analysis of studentized residuals (Figures 2 and 3) and validation statistics (Table 3), the most appropriate equation to represent the height-diameter relationship of Eucalyptus benthamii stands was Harrison_c (Equation 12). As occurred in this study, Téo et al. (2017) selected the Harrison et al. (1986) model adjusted with its parameter "c" expressed as a linear function of stand variables to represent the height-diameter relationship of Pinus taeda stands in the West region of Santa Catarina State. Soares & Tomé (2002) also selected the Harrison et al. (1986)  where: c = 0,620953 + 0,088471t + 0,011543d dom -0,012183G; h = total height (m); d = diameter at breast height (cm); h dom = dominant height of the sample plot (m); t = age of the forest stand (years); d dom = dominant diameter (cm); G = basal area of the sample plot (m².ha-¹); e = Euler's constant (2,718281829...).
The parameter "c" is also known as rate parameter of the Harrison et al. (1986) model, which guarantees that the height increase rate is smaller for the greatest diameter at breast height (d). This parameter "c" was expressed as a linear combination of age (t), dominant diameter (ddom) and basal area of the sample plot (G). These formulation and adjustment model attributes allow portraying the differences in height-diameter relationship of Eucalyptus benthamii, under the effect of these several stand variables, as well as, ensure more general applicability to the equation.
The hypsometric relationship of Eucalyptus benthamii represented by the Equation 12, with increasing age, moves to the right and rises of level following the diameter and height growth (Figure 4). The same pattern was found by Araújo et al. (2012), for the hypsometric relationship of Eremanthus erythropappus Mac Leish stands, with ages ranging from 4 to 8 years.
Concerning the stand density, the lower the densities, the higher the position of the height-diameter curves of Eucalyptus benthamii predicted by the Equation 12 (Figure 4). Similar results were found by Bartoszeck et al. (2004) when evaluating Mimosa scabrella Benth. stands, in municipality of Curitiba -PR, the height-diameter curves are located at a higher level for the lower density stands, at age of 7,5 years. However, Bartoszeck et al. (2004) reported for the others ages they evaluated, that the height-diameter curves were at the same level and with similar slope, so they concluded that the stand density did not affect the hypsometric relationship of Mimosa scabrella.

CONCLUSIONS
The Harrison et al. (1986) function, adjusted with the parameter "c" expressed as a linear combination of age (t), dominant diameter (d dom ) and basal area of the sample plot (G) was the most appropriate to represent the height-diameter curves of Eucalyptus benthamii stands, in the municipality of Irani -SC. The formulation and adjustment characteristics of the modeling procedure used in this study allow the representation of the height-diameter relationship of stands with different ages and stand densities.
The height-diameter curves move to the right and level up with increasing age, because of the diameter and height growth. The lower the stand density, the higher the level of the height-diameter curves, suggesting that trees with equal diameters are taller in the lower density forest stands.