Study area

The study was carried out in the stands under forest management at Umafor 2001, in the northeast of the state of Oaxaca, Mexico, which cover an area of 195 397 ha. The study area is located at the coordinates 17°22’58” N and 96°22’29” W (Figure 1). The altitude of the sampling areas varies from 1 800 to 3 361 m. The annual temperature ranges between -3 °C and 28 °C, with a mean annual precipitation of 3 797 mm (^{Granados-Sánchez, 2009}). The predominant vegetation type corresponds to pine-oak forests, with Pinus oaxacana, P. douglasiana, P. patula, P. pseudostrobus, P. leiophylla Schiede ex Schltdl. & Cham., P. pringlei Shaw., P. rudis Endl., P. teocote Schiede ex Schltdl. & Cham., Quercus crassifolia Bonpl., Q. laurina Bonpl. and Q. rugosa Née as the main dominant species.

Figure 1 Geographic location of the study area. 

Data

The dynamic SI and dominant height equations were fitted using a database of 206 tree stem analyses of 44 specimens of Pinus oaxacana, 45 of P. douglasiana, 58 of P. patula, and 59 of P. pseudostrobus. The individuals were felled at a stump height that ranged from 0.1 m to 0.25 m, so that the first section corresponded to the stump height, the second to a height of 1.3 m, and subsequently, every 2.54 m, up to the tip of the tree. The specimens were selected from mixed stands, and an attempt was made to cover the various site qualities, as well as the diameter categories present in the areas under forest management. Although the concept of site index is designed for coetaneous and pure stands (^{Daniel et al., 1979}), the dominant character for a tree implies that these individuals have little sensitivity to the effect of the density, the competition and the mix of species, as its height remains constant within a broad interval of densities and degrees of mingling, which reflects the productive potential of a particular site throughout its life, and therefore the concept has been extended to mixed pine stands (^{Corral et al., 2004}; ^{Vargas-Larreta, 2010}; ^{Castillo et al., 2013}).

Table 1 shows the descriptive statistics of the database; as it may be observed, the sample covers a broad interval of age-dominant height combinations, and therefore it is regarded as adequate for the fitting of the dynamic growth equations.

The true heights of the sections were estimated using ^{Carmean’s method (1972)} modified by ^{Newberry (1991)}, as it has exhibited good results in other researches (^{Dyer and Bailey, 1987}; ^{Fabbio et al., 1994}; ^{Castillo et al., 2013}). The method is based on two assumptions: i) between two sections, the tree grows at a constant rate, and ii) in average, the tree is cut at the growth core, at a one-year height (^{Dyer and Bailey, 1987}; ^{Amaral et al., 2010}). The equations used to estimate the true height vary according to the tree section (equations 1-3).

Where:

H ₁ and H ₂ = Heights of the lower and upper sections of the log

N ₁ and N ₂ = Number of rings of the lower and upper sections of the log

N ₀ = Age of the tree, i.e. number of rings of the stump (in stump N ₀ = N ₁ )

T ₀ = Age of the tree when it reached the height H_1; i.e. N ₀ - N ₁

T = Whole number from 1 to N ₁ - N ₂

Dynamic site index equations

The GADA method was used to fit three equations, which have been documented in various studies in order to describe the growth in dominant height and site index: Korf (4), Hossfeld (5) and Bertalanffy-Richards (6) (^{Cieszewski, 2002}; ^{Castillo et al., 2013}; ^{Quiñonez-Barraza et al., 2015}; ^{González et al., 2016}). Table 2 shows the expression of the base model, as well as their GADA formulation.

Description of the GADA method

The development of any dynamic equation using the GADA formulation considers: i) selecting a base equation and identifying the number of parameters that will be dependent on the productivity of the site; ii) the selected parameters are expressed as functions of the site quality defined by variable X ₀ (unobservable, independent variable describing the productivity of the site, as a result of the sum of the ecological and climatic factors, management regimes, soil conditions and the new parameters; iii) the selected bi-dimensional base equation (H=f(t)) is expanded to a tri-dimensional site index equation (H=f(t, X ₀ )), and iv) the value of X ₀ is solved in terms of the baseline conditions of the site, i.e. of the baseline dominant height values - age (t ₀ , H ₀ ), so that the model can be implicitly defined and applied in practice (H=f(t, t ₀ , H ₀ )) (^{Cieszewski and Bailey, 2000}; ^{Cieszewski, 2001}; ^{Cieszewski, 2002}; ^{Cieszewski, 2003}). In the course of the process, redundant parameters are often eliminated, and a model with an equal or lower number of parameters to that of the original base equation is obtained.

Fitting method and statistics for the comparison between models

The global and site-specific parameters were estimated using the nested iterative procedure (^{Tait et al., 1988}; ^{Cieszewski, 2003}; ^{Krumland and Eng, 2005}), which is utilized with databases derived from tree stem analyses or permanent sampling plots. The models were fitted with a continuous-time autoregressive error structure (CAR_i), in order to correct potential autocorrelation problems of the term of error (^{Zimmerman et al., 2001}; ^{Nord-Larsen, 2006}; ^{Crecente-Campo et al., 2009}) and obtain unbiased, efficient parameter estimators (^{Parresol and Vissage, 1998}).

Where:

Hij = Prediction of the height i using H _j (height j), t _i (age i), and t _j (age j≠i) as predictive variables

β = Vector of parameters to be estimated

eij = Error

d1=1 for j> 1 and zero when j=1

ρ1 = Autoregressive parameter of order 1, whose value is to be estimated

tij-tij-1 = Temporal difference separating the jth observation from observation jth-1 in each tree stem analysis, t_ij>t_ij-1

The goodness-of-fit of the models was determined based on the results of quantitative analyses, such as the root mean square error (RMSE), the adjusted coefficient of determination (R ² _adj ), the mean bias (e) and the residuals graphs for comparing the following aspects: i) the biological behavior of the predictions of the various fitted models; ii) the representation of the residuals versus the values predicted by the model, and iii) the representation of the residuals versus residuals with various lags, in order to verify the correction of the autocorrelation of the errors by modeling the error structure (^{Goelz and Burk, 1992}; ^{Sharma et al., 2011}). The presence of autocorrelation in the utilized data was calculated with the Durbin-Watson statistic (^{Durbin and Watson, 1951}). The expressions of these statistics are:

Where:

H,H^ and H- = Observed, predicted and mean values in pairs of height dominant data - age

n = Number of observations

p = Number of parameters of the model

e = Residual value of the fitted model

The simultaneous fitting of growth and error structure equations (CAR2) was carried out with the MODEL procedure of the SAS/ETS^TM statistical package (^{SAS, 2004}), which allows a dynamic updating of the residuals in the fitting procedure.

GADA equations generated without considering the error structure

The estimated parameters, the standard errors and the fitting statistics of the dynamic equations analyzed for Pinus oaxacana, P. douglasiana, P. patula and P. pseudostrobus were obtained for all the pairs of dominant height - age, with the fit of equations 4, 5 and 6, without considering the modeling of the error structure (Table 3). The results evidence autocorrelation problems of the residues, due to height measurements in the same tree (Figure 4).

Altura dominante = Dominant height; Edad(años) = Age(years).

Figure 2 Overlapping of SI curves generated using the Korf (discontinuous line), Hosffeld (dotted line) and Bertalanffy-Richards (continuous line) models.  

Altura dominante = Dominant height; Edad(años) = Age(years).

The curves were generated with the dynamic model of Bertalanffy-Richards, with (continuous line) and without (dotted line) correction of the autocorrelation using a CAR2 autoregressive structure.

Figure 3 Site index curves at the base age of 40 years for Pinus oaxacana Mirov (Po), P. douglasiana Martínez (Pd), P. patula Schiede ex Schltdl. & Cham. (Pp) and P. pseudostrobus Lindl. (Pps).  

Autocorrelación = Autocorelation; Dos errores estándar = Two standard errors.

Figure 4 Graphs of the Partial Autocorrelation Function of the Bertalanffy-Richards equation (Equation 6): a) fitting without correction of the autocorrelation, b) considering a CAR1 model, and c) considering a CAR2 model. 

The coefficients of determination (R ² _adj ) adjusted by the number of parameters of the models of Korf, Hossfeld and Bertalanffy-Richards for the different species show that it is possible to account for more than 96 % of the variance observed in the dominant height growth, in terms of the of the age, with estimations of significant parameters for most combinations of equations and species (p<0.01). The mean errors in the three models range between 1.18 m (Pd) and 1.85 m (Pp) (Table 3).

In general, the three fitted equations show satisfactory fitting statistics for predicting the dominant height growth in terms of age, without expanding the terms of error in order to correct the autocorrelation. On the other hand, the graphic analysis shows a realistic biological behavior (similar predictions to the observed data). However, in Pinus oaxacana the parameter b ₁ was not statiscally significant in the Korf (4) and Hossfeld (5) models; therefore, the null hypothesis that the parameter estimations are equal to zero is accepted; this scenario also occurred for the Bertalanffy-Rihards model (6) in parameter b ₂ . Parameter b ₂ was not significant in the Hossfeld model (5) for Pinus patula Scheide ex. Schltdl. & Cham.; this situation occurred for the same parameter in the Bertalanffy-Richards model (6) for Pinus pseudostrobus, all at a significance level of 5 %.

The results show that the differences between the statistics of goodness-of-fit were minimal between the evaluated models; therefore, the selection of the most appropriate model was based on the graphic analysis of families of dominant height growth curves superimposed to the growth trajectories of the observed data; SI categories of 15 m to 35 m were considered for Pinus oaxacana, of 10 m to 30 m for Pinus douglasiana, with 5 m intervals, and of 14 m to 38 m for Pinus patula and Pinus pseudostrobus, with 6 m intervals, at a reference base age for all cases (Figure 2). This results from the fact that different models can exhibit the same goodness-of-fit or comparison statistics, but a different response to the trajectories of the observed data; thus, some underestimate these at the earlier ages and overestimate them at the later ages, or vice versa (^{Diéguez-Aranda et al., 2006}).

The Betalanffy-Richards model (6) was selected as the most adequate for describing the growth trajectories of the trees included in the database, and in order to compare the growth estimations versus the same model, based on the modeling of the error structure, with CAR1 and CAR2 type continuous-time autoregressive error structures, as this model provides a slightly better description of the individual tendencies in dominant height growth, mainly at ages above 50 years (Figure 2).

The Bertalanffy-Richards equation, which considers a CAR2 structure, was able to model the growth tendencies of the experimental data for the studied species and describe variation patterns in the growth of the developed curves, as well as different site potentialities. This result is consistent with other studies which also recommend this model for the construction of site index curves for other species of commercial interest (^{Kiviste et al., 2002}; ^{Rodríguez-Carrillo et al., 2015}; ^{Castillo et al., 2013}).

Error structure modeling

Table 4 shows the fitting statistics of the Bertalanffy-Richards model; this incorporates a continuous-time autoregressive model of order 2 (CAR2), which adequately corrected the autocorrelation of the database errors (DW value near 2).

The average RMSE errors of model 6, with the error structure, ranged between 0.79 m (Pd) and 1.18 m (Pp) (Table 4), and accounted for more than 98 % of the variance in the dominant height growth. The continuous-time autoregressive model of order 2 (CAR2) ensured the independence of the residuals in the fitting process.

Figure 3 shows the SI statistics calculated with the Bertalanffy-Richards equation, without taking into account the error structure (dotted line), as well as those that do consider the CAR2 error structure (continuous line) at a base age of 40 years, overlapped on the experimental data utilized for fitting it. The SI curves cover various ages and describe the relationship of dominant height growth and age with a logical and biological behavior (complex polymorphism) at ages above 50 years (^{Cieszewski, 2003}; ^{Álvarez-González et al., 2010}). This figure shows, graphically, that the studied species exhibit various dominant height growth patterns, indicating that stands with dominance of one of the four pine species require to be managed with different forestry schemes (^{Attis et al., 2015}).

In the graphs below, the partial autocorrelation function (PACF) resulting from the fitting of the Bertalanffy-Richards equation with 10 lags of the residuals for each tree (Figure 4), is observed to change from a correlogram with a positive autocorrelation of order 2 to a correlogram without evidence of autocorrelation, estimated through the fit of the dynamic equation with a CAR2 type continuous-time autoregressive error structure; the results are similar to those registered by ^{Vargas-Larreta et al. (2010)}, ^{Vargas-Larreta (2013)} and ^{Quiñonez-Barraza et al. (2015)} for forest tree species in various regions of Mexico.

When correcting the autocorrelation, the random pattern of the residuals is defined around the zero line, with a homogeneous variance and without detecting a particular tendency (Figure 5). Also shown are the evolution of the bias and the RSME, as well as a considerable increase in the predicted dispersion and heights above the age of 50 years, due mainly to the lack of height-age data pairs for these age classes.

Residual = Residual; Altura predicha = Predicted height; Sesgo = Bias; Clase de edad (años) = Age category (years); REMC) RSME.

Figure 5 Residuals vs predicted heights and evolution of the bias and root mean square error (RSME) by age category in the Bertalanffy-Richards model, using the CAR2 autoregressive model.