Study area and dendrometric data

The study was conducted in forests on communal land of Ixtlán de Juárez, Oaxaca, México (17° 23´ 0.50´´, 17° 23´ 0.58´´ N, 96° 28´ 45´´, 96° 28´ 53´´ W, mean altitude 2780 m). The database was obtained from stem analysis of 34 P. ayacahuite individuals from cross sections extracted immediately after the stump, from diameter at breast height (DBH, 1.3 m above-ground) and after DBH, at each commercial length (2.62 m, 2.54 m and 1.25 m) up to the tree apex. The dendrometric information was age-height in each cross section and age-total height of the tree. With these data, the complete profile of each individual was reconstructed.

The data cover a wide interval of dominant heights and ages (Table 1). The statistical summary included mean, median, standard deviation and variable interval for fitting dominant height (HD) growth functions.

Biometric models used and derivation of self-referencing models

To forecast dynamic and biological behavior, we used five base models (guide curve) with non-linear structure (Table 2), classic in forest biometry.

In modeling HD growth and IS of P. ayacahuite, each model was expressed in its formulations under the guide curve and ADA (Table 3). For the case of the Hossfeld IV and Chapman-Richards models, two GADA formulations were posed (Table 4). These approaches permitted construction of different IS curves (anamorphic, polymorphic and asymptotic polymorphic).

Deriving the models with a guide curve involved fitting the base models to obtain the average curve in the entire interval of observations of HD - E, from which proportional curves were generated above and below the central curve in the following order: 1) each growth equation (in its integral form) was fitted to generate the guide curve of a base age (EB) from the predicted HD value, 2) the model was formulated in its site index expression, 3) the parameter of the growth hypothesis of interest was selected and solved from the IS expression, and 4) the solution of the parameter in the original equation (integral) was substituted.

The ADA method consisted of developing the equation to be fitted in its differentiated form, by expressing future dominant height (HD₂) as a function of future age (E₂), from the initial age (E₁) and initial dominant height (HD₁). The model of algebraic difference that defines the family of IS curves is expressed as HD₂ = f (HD₁, E₂, E₁, β), where HD₂ is the value of dominant height at age E₂ (age of the projection), HD₁ is the dominant height at age E₁ (initial age) and β is the vector of regression parameters (^{Clutter et al., 1983}; ^{Diéguez et al., 2005}; ^{Magaña et al., 2008}; ^{Santiago et al., 2013}).

In developing the GADA models, the following were considered: 1) selection of a growth model and identification of dependent parameters of site quality, 2) expression of the selected parameters as functions of the site quality defined by the variable X (not observable, independent variable that describes site productivity) and the new parameters, 3) the bidimensional base equation selected was expanded to an explicit site index tridimensional equation, and 4) the value of X was solved from the initial conditions of the site; that is, we parted from initial values of age and height, in such a way that the model was implicit and applicable (^{Cieszewski and Bailey, 2000}; ^{Cieszewski, 2002}; ^{Vargas et al., 2010}; ^{Quiñonez et al., 2015}).

Model fit

To fit the guide curve equations, 468 pairs of dominant height-age data were used. In ADA and GADA, fit was performed simultaneously to minimize errors of the equation system; HD₁ = f (E₁, b), HD₂ = f (HD₁, E₂, E₁, b) with 434 pairs of non-overlapping data. These fits were done using non-linear least squares with the procedure MODEL of SAS/ETS^® v. 9.0 (SAS Institute Inc., 2002). Graphic analysis was done with R Project 3.2.3 (^{R Development Core Team, 2015}).

Statistical indicators

The statistical indicators used to evaluate goodness of fit of each model tested are presented in Table 5. The graphic analysis of the predictive capacity and significance of the estimators of the parameters of each model were the basis of election of the best models. Significance was evaluated with the Student t test, under H₀: b _i = 0 and a critical value α = 0.05 (^{Infante and Zarate, 2012}).

where donde yi, y^i and y- are the observed and predicted values and the mean of the observed values, p is the number of parameters of the model and n the number of observations.

Guide curve

All the non-linear models generated adequate fit. The highest R²-adj values were found with the Levakovic II, Hossfeld IV and Chapman-Richards, explaining nearly 93.6 % of the total variance in dominant height. Moreover, the three models generated the lowest values of SCE and, consequently, of RECM. Each had significant parameters, with P value below 1 % (Table 6).

The graphic analysis indicated that the Schumacher model underestimated growth in early ages. In contrast, the Gompertz model overestimated growth up to 1.5 m at 2 years of age. According to the estimated parameters (Table 6), the asymptotic values of the Schumacher, Levakovic II and Hossfeld IV models were higher than those the species can reach (up to 40 m) (^{Perry, 1991}).

The Chapman-Richards model, which generated a more real asymptotic value, was selected for construction of the site index curves, and growth estimations were compared with the ADA and GADA approaches.

Dynamic equations

The ADA projection models in their anamorphic and polymorphic expressions had excellent fit (Table 7), and their parameters were significant (p ≤ 0.0001) (Table 8). The exception was the Schumacher model. The polymorphic models derived from Chapman-Richards (b ₂ ), Levakovic II (b ₂ ) and Hossfeld IV (b ₁ ) were outstanding for their higher value of the adjusted coefficient of determination (R²-adj) and for explaining approximately 99 % of the total variability. These models had lower SCE and CME values. The Chapman-Richard model b ₂ and Hossfeld IV model b ₁ had lower standard errors, as well as the most real asymptotic values observed. Thus, in this phase of the numerical analysis, the last two models were better than the others

According to the Student t test, the parameters of all the fitted models were statistically significant, with p ≤ 0.0001, lower than the critical value (α = 0.05). Thus, any value lower than this limit represented a rejection of H₀. The graphic analysis verified the results of the fit by showing the dynamic trend for P. ayacahuite (Figure 1).

Figure 1 Families of site index curves for Pinus ayacahuite at 40 years base age. A) and B): by the guide curve method, anamorphic expression (b ₀ ) and Chapman-Richards polymorphic b ₂ . C) and D): ADA-type dynamic growth curves, polymorphic Chapman-Richards b ₂ model and Hossfeld IV polymorphic b ₁ . E) and F): GADA-type dynamic growth curves of the Chapman-Richards and Hossfeld models, respectively. 

For the derivation of the GADA structures, the Chapman-Richards and Hossfeld IV models (Table 9), which have been widely used for modeling forest growth (^{Diéguez et al., 2006}; ^{Vargas et al., 2013}; ^{Tamarit et al., 2014}; ^{Quiñonez et al., 2015}), were considered.

R²-adj of the two models explained more than 98 % of the total variance. The values of SCE and CME were low, like those of the selected ADA models, and RECM had plausible values for both equations. All the estimations of the parameters were significant (p ( 0.01), and there were no major differences between the two adjusted models. The curves were compared to evaluate the performance of the models with GADA structure (Figure 1).

Families of site index curves for Pinus ayacahuite

From 12, 15, 18, 21 and 24 m to the base age of 40 years of the best models, the IS curves constructed from the models of similar fit followed different growth patterns (Figure 1A - F). Different models can have the same goodness of fit statistics or of comparison but different response (some underestimate in the early ages and overestimate with advanced ages and vice versa) (^{Diéguez et al. 2006}).

The superimposed trajectories of the IS data observed in a family of curves, in this case of the anamorphic type of the Chapman-Richards model, showed that the dominant height growth pattern of P. ayacahuite follows different trends as a reflection of polymorphism (Figure 1A). In this respect, ^{Corral et al. (2004)} point out that anamorphic functions are generally not the most adequate for representing dominant height growth because the curve shapes vary among sites. Therefore, dominant height growth is actually polymorphic. For this reason, anamorphic structure was ruled out for describing dominant height growth.

Moreover, larger discrepancies are observed at early ages. In this respect, ^{Huang (1999)}) suggested that erratic height growth at young ages is due to young trees’ reaction to abnormal climatic extremes that can affect their total height, which is frequently determined by factors different from site quality.

The polymorphic b ₂ derivation of Chapman-Richards and b ₁ of Hossfeld IV were robust in the ADA models, according to the goodness of fit indicators and the significance of the estimated parameters (Tables 7 and 8). The graphic results (Figure 1C and D) show that both models satisfactorily group real observed data of all site qualities. Nevertheless, the model derived from Chapman-Richards grouped the data on early ages better than the Hossfeld IV model. Therefore, the Chapman-Richards polymorphic b ₂ structure was superior to the Hossfeld IV b ₁ derivation for predicting dominant height growth. The latter model was similar to that of Hossfeld IV, selected with the GADA methodology, which had adequate goodness of fit and graphic performance. In this way, the two specified ADA models were selected.

The numerical analysis of the models derived from Chapman-Richards and Hossfeld IV with the GADA formulation did not show major differences. However, the Hossfeld IV model was more realistic in the graphic result because it generated predictions with trends similar to the observed profiles. For sites of better quality, the Chapman-Richards overestimated growth; the opposite occurred in low quality sites (Figure 1E). For this reason, we concluded that to predict dominant height growth of P. ayacahuite with GADA structure, the model derived from Hossfeld IV is adequate. Consequently, to compare growth and productivity estimations of P. ayacahuite stands with the three methodologies, the selection is the following: Chapman-Richards polymorphic b ₂ model with the guide curve, Chapman-Richards polymorphic b ₂ and Hossfeld IV b ₁ in ADA and Hossfeld IV in GADA.

The Chapman-Richards and Hossfeld IV models are widely used to model forest growth by the guide curve method and their versions of algebraic difference and generalized algebraic difference approaches. As examples for the Chapman-Richards model, ^{Diéguez et al. (2006)} constructed dynamic growth curves for Pinus tadea L. in the USA, and ^{Vargas et al. (2010)} studied Pinus cooperi Blanco in Durango, Mexico, with the compared models and derivations of Korf and Hossfeld IV. In addition, ^{Quiñonez et al. (2015)} developed a GADA expression to compare it with the structure derived by ^{Kumland and Eng (2005)} in forests of Durango.

^{Gómez et al. (2009)} exemplified the Hossfeld IV model to describe dominant height in plantations of E. grandis Hill ex Maiden and E. urophylla S. T. Blake, in Oaxaca, and ^{De los Santos et al. (2006)} evaluated growth in plantations of Terminalia amazonia (Gmel.) Excell in two regions of Costa Rica. ^{Santiago et al. (2013)} investigated the construction of a compatible system of growth and yield for even-aged stands of Pinus patula Schl. et Cham.

Graphs of the residuals for the models selected with the three methodologies are presented in Figure 2. We can see that in the ADA and GADA dynamic models, the dispersion of the residuals was less (-4 to 3 m), relative to the model selected with the guide curve (-6 to 8 m). Moreover, the latter model showed a slight funnel-shape tendency at early ages; the distribution of the residuals of the other models did not exhibit heteroscedasticity and showed a random pattern of the residuals around the reference line (line zero). This is due to the nature of the data used in the ADA and GADA dynamic equations, which enabled simultaneous fit and comparison of the change rates. For this reason, the studies conducted for construction of site index curves were developed with the methodology of the equations in algebraic differences or in its generalized form (^{Castillo et al., 2013}; ^{Vargas et al., 2013}; ^{Quiñonez et al., 2015}).

Figure 2 Dispersion of residuals against predicted heights. Models: A) Chapman-Richards polymorphic b ₂ by the guide curve method, B) Chapman-Richards polymorphic b ₂ under ADA, C) Hossfeld IV polymorphic b ₁ under ADA, and D) Hossfeld IV under GADA.