Introduction

Site quality can be defined as the productive capacity of a particular site for the growth of trees of a given species in response to the totality of existing environmental conditions; it commonly refers to the volume of wood a stand produces in a rotation forest (^{Diéguez et al., 2009}; ^{Skovsgaard & Vanclay, 2008}). Site quality assessment methods are divided into direct and indirect approaches. Direct methods are based on observations of dominant or co-dominant height; this variable is unaffected by stand density and is a fundamental component in growth and yield models (^{Pyo, 2017}; ^{Seki & Sakici, 2017}). On the other hand, indirect ones use the relationships between canopy species, the characteristics of undergrowth vegetation, topographic, climatic and edaphic factors, and the chemical composition of foliage (^{Chen, Klinka, & Kabzems, 1998}; ^{Wang, 1998}).

To characterize the forest productivity of the stands, the dominant height expected at a certain age is used and is called the site index (SI) (^{Martín-Benito, Gea-Izquierdo, Del Río, & Cañellas, 2008}). Modeling dominant height and SI is an important tool for classifying the productivity of forest lands and defining management strategies in the application of silvicultural treatments (^{López-Sánchez, Álvarez-González, Diéguez-Aranda, & Rodríguez-Soalleiro, 2015}; ^{Quiñonez-Barraza et al., 2015}). One of the most commonly used methods for generating SI curves is the guide curve, which characterizes the height status to determine the current average condition; in contrast, dynamic equations estimate dominant height as a function of current and future age and height (^{Tamarit-Urias et al., 2014}). Dynamic equations are based on the algebraic difference approach (ADA) and generalized algebraic difference approach (GADA). In the ADA method, a parameter of the equation is considered site-specific. Depending on the site-dependent parameter, anamorphic or polymorphic curves are obtained (^{Bailey & Clutter, 1974}). In the GADA equations, two parameters are site-specific and generate families of polymorphic curves with multiple asymptotes (^{Cieszewski & Bailey, 2000}).

Dummy variable (DV) and mixed-effects model (MEM) approaches have been applied in ADA and GADA equations to model dominant height and SI (^{De los Santos-Posadas, Montero-Mata, & Kanninen, 2006}; ^{Nigh, 2015}; ^{Tamarit-Urias et al., 2014}; ^{Wang, Borders, & Zhao, 2007}); results have shown ability to predict and project height as a function of age. In the DV approach, global and local parameters are generated, and in the MEM approach, parameters are fixed and random for each tree or plot. These approaches fit global (common) and specific (local) parameters simultaneously, but differ in the form of estimation (^{Wang, Borders, & Zhao, 2008}). The former is characterized by the addition of dummy variables to the parameter that explains the effect of the tree or plot; it considers local parameters as fixed, but different for each tree or plot. On the other hand, the MEM approach provides a mean response if only the fixed parameters and a specific response per sampling unit, which together form the mixed-effects parameters, are considered. These are estimated by defining a variance-covariance matrix in the model structure (^{Calama & Montero, 2004}; ^{Fang & Bailey, 2001}).

In the state of Michoacan, Mexico, Pinus pseudostrobus Lindley is the most important forest species in economic terms, which is why commercial forest plantations have been established that represent a component in timber production (^{López-Upton, 2002}; ^{Sáenz-Romero et al., 2012}). Therefore, it is necessary to develop accurate silvicultural techniques for the planning and execution of sustainable forest management activities. In this context, the objective of the research was to fit and compare dynamic dominant height and SI equations with the DV and MEM approaches for P. pseudostrobus in commercial forest plantations in Nuevo San Juan Parangaricutiro, Michoacan, Mexico.

Materials and methods

Results and discussion

Table 2 shows the estimators and standard errors of the global parameters for the DV approach and the fixed ones for MEM. All parameters were different from zero (P < 0.0001). The means of the local and mixed parameters were similar in each fitted equation; the values ranged from 21.71 to 23.76 and the variance from 4.21 to 9.59, respectively (Table 3).

According to the comparison of the fitting of the modeling approaches with ADA equations, the ADA2 equation (equation 2) with DV obtained the maximum value in LL (-1 014.74) and the lowest in AIC (2 115.47). On the other hand, the GADA equation (equation 4) with DV presented values of -1 015.04 and 2 118.07, in the LL and AIC, respectively, which were similar to those obtained for ADA2 (Table 3). In the equations with DV, the local parameters were totally independent in comparison to MEM, where the random parameters were added to the fixed parameter. This explained the favorable results for DV in terms of LL and AIC, which can be considered as a modified version of LL to take into account the effect of the number of model parameters, even if the same results are obtained in terms of goodness-of-fit (^{Wang et al., 2008}).

In general, the ADA equations with DV were relatively higher than those fitted with MEM. In the ADA equations, the R² value was higher (98.58 %) and the RMSE one lower (1.17 m) in the ADA2 equation (equation 2) with DV, although with MEM it had the lowest bias. On the other hand, the MEM-fitted GADA equation (equation 4) had the highest R² value (98.66 %) and the lowest RMSE one (1.14 m); however, the bias was lower with the DV approach (Table 3).

Figure 1 shows that the distribution of local and mixed parameters resembles a normal one and that it was similar for DV and MEM. The random MEM parameters for the ADA group showed a normal distribution with a different degree of kurtosis, and the mean values were close to 0 (−7.32×10⁻¹¹, −9.09×10⁻¹⁶, 7.31×10⁻¹¹ for ADA1, ADA2 and ADA3), while GADA had a value of -0.0030 with the Gauss optimization method (Figure 2). The condition of obtaining a normal distribution of the random parameters with zero mean and known variance is an aspect that must be considered; if this does not happen, it could be a reason to use DV, since in this approach there are no considerations about local parameters, as they are independent of each other and there are no restrictions on the values they can take (Nigh, 2015; Verbeke & Molenberghs 2000; ^{Wang et al., 2008}).

Figure 1 Normal distribution function of the local (Hs _i ) and mixed ( Hs+Hsi ) parameters of the fitted equations (algebraic difference [ADA] and generalized algebraic difference [GADA]) with the dummy variable (DV) and mixed-effects model (MEM) approaches. 

Figure 2 Normal distribution function of the random parameters (As _i ) of algebraic difference (ADA) and generalized algebraic difference (GADA) equations with mixed-effects models (MEM).  

Figure 3 shows the families of dominant height growth curves with the modeling approaches and the ADA and GADA equations. The SI categories were symmetrical: 18, 22, 26 and 30 m at a base age of 20 years.

Figure 3 Families of dominant height growth curves of Pinus pseudostrobus for the site index (SI) categories 18, 22, 26 and 30 m, from the fitted algebraic difference (ADA) and generalized algebraic difference (GADA) equations with the dummy variable (DV) and mixed-effects model (MEM) approaches. 

At ages from one to 10 years, an underestimation of the dominant height was observed with the curves obtained from the ADA1 equation (equation 1), while for ADA3 (equation 3) an overestimation was observed. The ADA2 (equation 2) and ADA3 equations generate polymorphic curves, and GADA (equation 4), polymorphic with different asymptotes. The use of polymorphic models has been suggested in the literature to adequately represent the age-height relationship (^{Álvarez, Ruiz, Rodríguez, & Barrio, 2005}; ^{Kahriman, Sönmez, & Gadow, 2018}; ^{Tewari, Álvarez-González, & García, 2014}); however, ADA3 on the 30-m SI curve showed a disarticulation on the basis of the data trend. The trajectories of the GADA-generated curve families better described dominant height growth behavior for the SIs, at the base age of 20 years. The SI values estimated with GADA with DV for the communities of Pario, Huerekutini, Tejamanil I and Tejamanil II were 34, 32, 37 and 36 m, respectively. The commercial forest plantations of Pario and Huerekutini are located in the northeast of Nuevo San Juan Parangaricutiro and have similar topographic characteristics with slopes of 11 to 20 %, while Tejamanil I and Tejamanil II are located in the southwest and slopes vary from 12 to 23 %.

The dynamic equations based on the Chapman-Richards model had results similar to those reported by Pacheco, Santiago, Martinez, and Ortiz (2016), who indicated that the ADA2 and GADA equations were more accurate; however, the latter better described the data and covered the dominant height of P. montezumae with greater amplitude; ^{Quiñonez-Barraza et al. (2015)} generated curve families with the GADA equation used in this study, which adequately described the growth of the dominant height of Pinus species in mixed-species stands of Durango. Similarly, González, Cruz, Quiñonez, Vargas, and Nájera (2016) selected the a GADA equation because it presented better goodness-of-fit to predict height growth of P. pseudostrobus.

Figure 4 shows that the trend of the average bias by age category for the equations was similar between modeling approaches and that the highest values occurred in the two-year category. The distribution of bias was higher in the six- to 14-year categories with ADA3 (equation 3), while the values were homogeneous for ADA1 (equation 1), ADA2 (equation 2) and GADA (equation 4).

Figure 4 Average bias by age category of the fitted algebraic difference (ADA) and generalized algebraic difference (GADA) equations for Pinus pseudostrobus with the dummy variable (DV) and mixed-effects model (MEM) approaches. 

The GADA equation with DV and MEM better described the growth in dominant height; in addition, it presented the greatest efficiency in the fitting statistics. In forest modeling, these approaches have allowed the development of accurate equations with parameters unique to a tree or plot (^{Jerez-Rico et al., 2011}; ^{Pyo, 2017}; ^{Sharma, Subedi, Ter-Mikaelian, & Parton, 2015}; ^{Tamarit-Urias et al., 2014}), which coincides with the results of this study. Although the results of the modeling approaches were similar, the GADA equation with the DV approach had the greatest accuracy; therefore, it is considered that it can be used in the description of the dominant height growth of P. pseudostrobus for the commercial forest plantations studied. The results were similar to those found by ^{Wang et al. (2007}, ²⁰⁰⁸), who concluded that DV and MEM had an almost equivalent yield for dominant height and SI for Pinus tadea L. According to Nigh (2015), a GADA model presented better results with the DV approach for Picea engelmannii Parry ex Engelm. On the other hand, in biomass studies for Pinus massoniana Lamb., both approaches were statistically accurate; however, the use of MEM was more efficient (^{Fu, Zeng, Tang, Sharma, & Li, 2012}).

Conclusions

The dummy variable and mixed-effects modeling approaches were statistically accurate for the algebraic difference (ADA) and generalized algebraic difference (GADA) equations and allowed modeling dominant height growth and generating site index families for P. pseudostrobus. Results suggest that the ADA2 polymorphic equation adequately describes the growth pattern and projects maximum height growth in a rotation forest; however, the GADA equation satisfactorily covered the dominant height growth trajectories and had the highest accuracy with the dummy variable approach. Therefore, the use of the GADA equation with the dummy variable approach allows classifying the productive potential of commercial forest plantations in Nuevo San Juan Parangaricutiro, Michoacan, which have higher site index values in the Tejamanil I and Tejamanil II communities, and lower values in Pario and Huerekutini.