Introduction

The planning, implementation and monitoring of sustainable forest management require research to support decision-making and the assessment of the established goals. The estimation of the timber stock and productivity of the stands is a main objective in forest management systems; therefore, it is essential to know the growth of commercial species (^{Aguirre-Calderón, 2015}; ^{Salas et al., 2016}).

Research on the estimation of volume, growth and increment is a key tool for understanding the dynamics of ecosystems involved in forest management; therefore, these approaches continue to be necessary for the planning and implementation of forest activities. Taper and dominant height growth equations have been widely explored topics (^{Castillo et al., 2013}; ^{Paulo et al., 2015}; ^{Rodríguez-Carrillo et al., 2015}; ^{Krisnawati, 2016}; ^{Corral-Rivas et al., 2017}; ^{Fierros-Mateo et al., 2017}; ^{Tamarit et al., 2017}) because they are biometric tools that allow characterizing the profile of the trees, as well as estimating the total volume and the commercial volume, and the forest productivity, based on the site index approach (^{Crecente-Campo et al., 2013}; ^{Quiñonez-Barraza et al., 2015}; ^{Özçelik and Crecente-Campo, 2016}).

Since the datasets used for fitting taper and height growth equations are time series obtained from measured variables on the same tree and resulting from taper and tree stem analyses, it is reasonable to assume that the observations in each tree are correlated and, also the residuals of the adjusted equations (^{Arias-Rodil et al., 2015}; ^{Quiñonez-Barraza et al., 2015}; ^{Corral-Rivas et al., 2017}); furthermore, their predictions might show a variation in the levels of the predictor variables, which is known as heteroscedasticity (^{Gujarati and Porter, 2011}).

The term autocorrelation refers to the correlation between the residuals of a regression model when series of observations arranged in time are used, e.g. time-series data, or in space, and with cross-sectional data. The linear and nonlinear regression models are based on the theoretical assumption that the residues have the same variance and are, therefore, homoscedastic. The presence of autocorrelation and heteroscedasticity leads to estimations of non-minimum variance parameters and to somewhat unreliable prediction intervals, especially in modelling taper and volume equations (^{Fortin et al., 2007}; ^{Xu et al., 2014}; ^{Tang et al., 2016}). Consequently, the usual t or F tests are not valid. Therefore, the use of generalized least squares (GNLS) approach with variance functions and correlation structures is an alternative to generate the best unbiased linear parameter estimates (^{Gujarati and Porter, 2011}).

In studies of taper and dominant height growth models, correlation structures and power functions to correct the autocorrelation and heteroscedasticity of the residuals, respectively, are frequently used (^{Quiñonez-Barraza et al., 2014}; ^{Quiñonez-Barraza et al., 2015}; ^{Sharma et al., 2015}; ^{Özçelik and Crecente-Campo, 2016}; ^{Corral-Rivas et al., 2017}; ^{Tamarit et al., 2017}); continuous-time autoregressive structures of errors with one, two or three power functions with known variance are highlighted.

Because it is important to improve the predictive capacity and the interpretation of the statistic properties in equations fitting, the objective of this study was to evaluate the combination of variance functions with correlation structures, in order to model the heteroscedasticity and the errors dependence in taper and dominant height growth equations of Pinus teocote Schiede ex Schltdl. & Cham. in Durango, Mexico.

Materials and Methods

Results and Discussion

The combinations of variance functions and correlation structures generated 36 models for taper, and 36 for height growth, based on equations 1 and 2, respectively. The fit statistics and the ranking score (RS) showed the goodness-of-fit of the equations (Table 2, for taper, and Table 3, for height growth). Also the DwM statistics and probability of Bartlett’s test of homogeneity of variances (P-value).

The ranking score exhibited the combinations of variance functions with correlation structures by hierarchical order (^{Sakici et al., 2008}). The lowest RS value was the statistically best combination, and the highest RS corresponded to the worst combination, based on the sum of the ranks for each fitting statistic (^{Tamarit et al., 2014}).

In the case of taper, the varPower function combined with corARMA3-2 performed best, while the worst combination was with corCAR1. In all cases of correlation structures combined with varPower, values ranges from 1.634 to 1.862 were displayed for the DwM statistic. However, for the test of homogeneity of variances, all the combinations were heteroscedastic. In the variance function varExp, all the combinations with correlation structures, except for corCompSymm, were homoscedastic (P>0.01). And the corARMA1-1 structure performed best according to the RS, with a DwM value of 1.553.

The varConstPower and varComb functions combined with the correlation structures have consistent statistics and independent residuals (DwM values from 1.722 to 1.973); however, they exhibit heteroscedasticity in Bartlett’s test. Structures corARMA3-1 and corARMA2-1 were the best for varConstPower and varComb, respectively.

In combinations of variance functions and correlation structures, the height growth generated DwM statistics of approximately 1.12. However, for most equations, they exhibited homogeneous variances. The varPower function is likewise statistically combined with the corARMA2-1, corARMA2-2 and corARMA3-2 structures, since the RS was 22 for all three cases. The corCAR1 structure had the lowest fitting, with unequal variances. Table 3 shows the behavior of the variance functions varExp, varConstPower and varComb with the combinations of correlation structures.

In all the fitted equations, the estimated parameters were statistically different from zero at a significance level of 5 % (P<0.05). The estimated parameters for the best combinations of each variance function with the correlation structures are summarized in Table 4, for both taper (T) and height growth (HG). In the case of taper, the combinations varExp&corARMA3-2, varExp&corARMA2-2 and varConstPower&corARMA3-1 represented 14 parameter estimates, and the test of homogeneity of variances displayed unequal variances; while the combination varExp&corARMA1-1, with 12 parameter estimates, exhibited constant variances.

The estimated parameters defining the changes of the dendrometric stem shapes of the segmented taper model were found for the change from neiloid into paraboloid, as 4.0 % to 7.2 %; while the change from paraboloid into cone occurs in average at 70 % of the total height in the tree stem profiles. Similar results have been documented by researchers for different Pinus species (^{Uranga-Valencia et al., 2015}; ^{Özçelik and Crecente-Campo, 2016}; ^{Corral-Rivas et al., 2017}; ^{Tamarit et al., 2017}), with values range from 4 % to 7 % for the first inflection point, and range from 55 % to 75 % for the corresponding second inflection point. Furthermore, continuous-time autoregressive structures of orders 1 and 2 and power functions were used, and a known variance was assumed. Nevertheless, the studies do not include the corresponding test of homogeneity of variances.

The tendency of the residuals by relative (h/H; %) of the four best combinations of variance functions with correlation structures for taper are shown as box and whisker plots (Figure 1). Although the tendencies are very similar, only the varExp&corARMA1-1 (H2-A5) combination exhibited a constant variance, and the dispersion of the extreme values are shown closer to the values of the corresponding quantiles.

H1 = varPower; H2 = varExp; H3 = varConstPower; H4 = varComb; A5 = corARMA1-1; A6 = corARMA2-2; A8 = corARMA3-1; A9 = corARMA3-2.

Figure 1 Box and whisker plots for the distribution of the taper residuals by relative height for combinations of varFunc and corStruct.  

The combinations in the height growth equations with the largest number of parameter estimates (9) were varExp&corARMA3-2 and varConstPower&corARMA2-2; the varExp&corARMA2-1 and varComb&corARMA2-0 combinations had seven estimated parameters.

The use of the equations with multiple parameters can cause an overparameterization, and the predictions resulting from them could not be the most efficient (^{Gregoire and Schabenberger, 1996}). In the four cases, the variances were constant, which guarantees that the parameters are unbiased and efficient and have a minimum variance (^{Gujarati and Porter, 2011}; ^{Tang et al., 2016}).

The selected combinations of the height growth equations generated a dispersion of residuals close to zero line, which appear in the form of boxes and whisker plots by relative height (h/H, %) in Figure 2. Although there is an overestimation from 1 % to 10 % relative height classes, and an underestimation from 10 % to 25 % relative height classes, the variances were constant for all four cases, according to Bartlett’s test of homogeneity of variances (^{Bartlett, 1954}; ^{Hidding et al., 2013}).

As for the residuals, although the DwM test generated values of approximately 1.2, they were considered to be independent within the measurements of height and age carried out in the same tree (^{Crecente-Campo et al., 2013}; ^{Quiñonez-Barraza et al., 2015}), and therefore the estimated parameters are unbiased, consistent, and efficient (^{Williams et al., 2013}).

In dominant height growth and site index equations with an algebraic difference approach (ADA) or generalized (GADA) equations, the power or exponential functions have been used for the correlation of the heteroscedasticity, in which unequal variances are assumed to exist in the generalized least squares adjustment process (^{Castillo et al., 2013}; ^{Quiñonez-Barraza et al., 2015}; ^{Rodríguez-Carrillo et al., 2015}; ^{González et al., 2016}), in addition to continuous-time autoregressive structures of the errors, in order to correct the autocorrelation with tree stem analysis of various species in natural forests and commercial forest plantations, such as Pinus durangensis Martínez, P. arizonica Engelm., Juniperus deppeana Steud, P. teocote Schiede ex Schltdl. & Cham. and P. ayacahuite Ehrenb. ex Schltdl. The approach developed in this paper considers the correction of the heteroscedasticity and of the autocorrelation as combinations of functions (^{Rodríguez et al., 2013}; ^{Williams et al., 2013}).

H1 = varPower; H2 = varExp; H3 = varConstPower; H4 = varComb; A4 = corARMA2-0; A6 = corARMA2-2; A7 = corARMA2-2; A9 = corARMA3-2.

Figure 2 Box and whisker plots for the distribution of the residuals of height growth by relative height for combinations of varFunc and corStruct. 

Figure 3 contrasts the predictions of taper and height growth equations with the selected combinations of power variance functions and correlation structures; it shows the observed tendency of the profile of a tree, the fitted equation without heteroscedasticity and autocorrelation correction (NC) and the combinations H1-A9, H2-A5, H3-A8 and H4A6. Only the combination H2-A5 exhibited constant variances. With regard to the height growth equation, the comparison between three SI curves for the fitted equation without correction (IS10, IS14 and IS18) and also the equations with corresponding combinations H1-A6, H2-A9, H3-A7 and H4-A4 are graphically represented. Furthermore, this figure shows the growth tendencies of the trees in dataset; in this case, all the combinations displayed constant variances, whereby the desirable properties in the estimated parameters are guaranteed (^{Beale et al., 2010}; ^{Vogelsang, 2012}).

H1 = varPower; H2 = varExp; H3 = varConstPower; H4 = varComb; A4 = corARMA2-0; A5 = corARMA1-1; A6 = corARMA2-1; A7 = corARMA2-2; A8 = corARMA3-1; A9 = corARMA3-2.

Figure 3 Taper prediction charts for a tree profile, and height growth curves by site index (SIs 10, 14 and 18 m at the baseline age of 60 years) with the combinations of variance functions and autocorrelation structures, and the equations without variance or autocorrelation (IS10, IS14, IS18) structures. 

Conclusions

The variance functions in combination with the correlation structures corrected the heteroscedasticity assumptions of variances and error autocorrelation in the taper and dominant height growth equations, using a generalized nonlinear least square approach; as a result, unbiased parameters with a minimum variance were obtained.

The predictions of the selected taper equations are more efficient, with consistent fitting statistics; the combination of the exponential variance function with an autoregressive-moving average correlation structure (corARMA1-1) produces constant variances by relative height categories of the stem profiles. The dominant height and site index model, at a base age of 60 years, exhibits realistic predictions. The selected combinations (H1-A6, H2-A9, H3-A7 and H4-A4) exhibit constant variances by relative height categories at specific ages. In both, the residuals are independent, and the properties of the tests of hypothesis of the estimated parameters are guaranteed.

The use of compatible taper and commercial volume equations, as well as of height growth and index site equations, is defined by the use of the intrinsic parameters in each equation; therefore, the parameters of the variance functions and correlation structures are only statistical indicators for rendering the fitting more efficient.