Generalized Height-Diameter Models for Pinus montezumae Lamb. and Pinus pseudostrobus Lindl. Plantations in Michoacan, Mexico

Jonathan Hernández-Ramos^1,2, Valentín José Reyes-Hernández², Héctor Manuel De los Santos-Posadas², Aurelio Manuel Fierros-González², Enrique Buendía-Rodríguez^3*, Gerónimo Quiñonez-Barraza⁴

¹Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias (INIFAP)-Campo Experimental Bajío, Celaya, Guanajuato, México.
²Colegio de Postgraduados (COLPOS), Campus Montecillo, México.
³INIFAP, Campo Experimental Valle de México, Texcoco, México.
⁴INIFAP, Campo Experimental Valle del Guardiana, Durango, México.

DOI: 10.4236/ojf.2024.143014                

Abstract

Tree height (H) in a natural stand or forest plantation is a fundamental variable in management, and the use of mathematical expressions that estimate H as a function of diameter at breast height (d) or variables at the stand level is a valuable support tool in forest inventories. The objective was to fit and propose a generalized H-d model for Pinus montezumae and Pinus pseudostrobus established in forest plantations of Nuevo San Juan Parangaricutiro, Michoacan, Mexico. Using nonlinear least squares (NLS), 10 generalized H-d models were fitted to 883 and 1226 pairs of H-d data from Pinus montezumae and Pinus pseudostrobus, respectively. The best model was refitted with the maximum likelihood mixed effects model (MEM) approach by including the site as a classification variable and a known variance structure. The Wang and Tang equation was selected as the best model with NLS; the MEM with an additive effect on two of its parameters and an exponential variance function improved the fit statistics for Pinus montezumae and Pinus pseudostrobus, respectively. The model validation showed equality of means among the estimates for both species and an independent subsample. The calibration of the MEM at the plot level was efficient and might increase the applicability of these results. The inclusion of dominant height in the MEM approach helped to reduce bias in the estimates and also to better explain the variability among plots.

Keywords

Random Covariate, Random Effects, Variance Structure, Forest Inventories, Forest Management, Mixed Models

1. Introduction

The height of trees (H; m) in forest stands or plantations (FP) is an attribute that is useful to calculate timber estimates, to evaluate site productivity at the stand or plot level, to project growth and yield scenarios, and to create distribution tables of usable products (Crecente-Campo et al., 2010; Santiago-García et al., 2017; Santiago-García et al., 2020) . In a forest inventory, measuring the H of all the trees within a sampling unit involves an important demand of resources, time, and effort; therefore, a subsample of tree heights is usually obtained (Mehtätalo, 2005; Bronisz & Mehtätalo, 2020) .

A reliable and low-cost alternative for error reduction in tree height estimation and to reduce the costs of measuring this variable in the field, is using local or generalized height—diameter (d; cm) models (Zambrano, Suarez, & Jerez, 2001; Barrio, Álvarez, & Díaz-Maroto, 2004; Diéguez-Aranda, Barrio, Castedo-Dorado, & Álvarez, 2005) . Besides, acquiring the data necessary to make calculations at the population level does not imply additional costs, since these are mostly office work to further estimate the dimensions of the remaining trees in each plot, and the generalized expressions consider characteristics at the stand level (Diéguez-Aranda et al., 2005; Crecente-Campo et al., 2010; Gómez-García et al., 2015) .

The generalized H-d models also enable the evaluation of tree height response to different growth conditions, forest management, or forest stand characteristics (Sharma & Parton, 2007; Crecente-Campo et al., 2010; Bronisz & Mehtätalo, 2020) . Statistical fit of these models has been performed with nonlinear least squares (NLS) (Misir, 2010; García-Cuevas et al., 2013) , mixed effects models (MEM) (Vargas-Larreta, Castedo-Dorado, Álvarez-González, Barrio-Anta, & Cruz-Cobos, 2009; Corral-Rivas, Álvarez-González, Crecente-Campo, & Corral-Rivas, 2014; Corral, Silva, & Quiñonez, 2019) , and neural networks (Correia et al., 2018) .

The mathematical structure of MEMs is characterized by fixed and random parameters; the former is related to an average response of the data used, and the second represents the difference of the average plot with respect to each unit or classification variable (Sharma & Parton, 2007; Gómez-García et al., 2015) . The application of MEM is a statistically improved alternative to NLS, which allows the inclusion of covariates that decrease the error and reduce the specific variability by classification level (Baty et al., 2015; Corral et al., 2019; Flores-Ayala et al., 2023) . Consequently, the variability of the response (i.e. height) to different growth conditions can be efficiently explained and estimated with greater certainty (Corral et al., 2019; Bronisz & Mehtätalo, 2020) . In addition, they have been an effective alternative solution to the lack of independence of measurements within a sampling unit (Calama & Montero, 2004; Corral et al., 2019) and for the reduction of the estimation error through the use of NLS in the optimization of parameter estimates (Bronisz & Mehtätalo, 2020) .

In Mexico, owing to the practical utility of generalized H-d models in forest management, this type of structure has been fitted to different coniferous and broadleaf species in natural temperate forests in Durango (Vargas-Larreta et al., 2009; Corral-Rivas et al., 2014; Corral, Silva, & Quiñonez, 2019) , Michoacán (García-Cuevas et al., 2013; Hernández et al., 2018) , and Hidalgo (Hernández et al., 2015) , as well as to even-aged forest stands in Ixtlán de Juarez, Oaxaca (Santiago-García et al., 2020) and in Ixhuacan de los Reyes, Veracruz (Hernández et al., 2020) . However, to our knowledge, there is no research work that had been performed in pine plantations, in particular for the two species that are analyzed in this research, and it is anticipated that these results would be better than those reported for natural forests owing to the homogeneity of forest plantations.

In Michoacán, Mexico, pine species have been established in 16,000 ha of FPs (CONAFOR, 2018) . Pinus pseudostrobus Lindl. and Pinus montezumae Lamb. are the most frequently used species as a result of their high growth rates, their potential to successfully expand the area already established, and consequently to increase Michoacán’s timber yield (Muñoz, Sáenz, García, Hernández, & Anguiano, 2011; Muñoz, Sáenz, García, Coria, & Muñoz, 2015) . These two Pinus species play also a major role in the resin, sawmill, furniture, and handicrafts industries of Nuevo San Juan Parangaricutiro (NSJP) (González, Gasca, & Heredia, 2014) . The objective of this research was to fit and propose a generalized H-d model for Pinus montezumae and Pinus pseudostrobus in the FPs of NSJP, Michoacán, Mexico with a mixed effects modelling approach.

2. Materials and Methods

2.1. Study Area

The FPs evaluated are located in the forests of Nuevo San Juan Parangaricutiro, Michoacán, between 19˚17' and 19˚30'N and 102˚06' and 102˚17'W, at an altitude ranging from 2000 to 3100 m. The region has a rugged terrain topography, with an extrusive igneous geology and an Andosol-type dominant soil—typical of the physiographic sub-province known as sierra Neovolcánica Tarasca. The Cw-type climate has a 1,500 mm average annual precipitation and 15˚C average temperature (INEGI, 2017) . The study area is shown in Figure 1.

2.2. Data Collection Methods

In this research, fourteen 7 - 32 years old Pinus montezumae and thirty-four 7 - 37 years old Pinus pseudostrobus plantations were analyzed. One-hundred

Figure 1. Area of distribution of the pine plantations that were evaluated. (a) Location of the Pinus pseudostrobus Lindl. and Pinus montezumae Lamb plantations evaluated, (b) location the Mexico in the North America, and (c) location of study area in the Mexico.

ninety-one 400 m² (20 × 20 m) quadrangular sampling plots were systematically located in these plantations. The d and H of each tree within each plot were measured, and the mean height (mH; m) was further calculated. The average maximum dominant height of the five tallest trees (dH; m) was estimated with the information from each plot, and this was defined according to the concept of dominant height that uses the 100 thickest trees per hectare (Alder, 1980) ; also, root mean square diameter (Msd; cm), basal area per plot (Ba; m²∙plot⁻¹) and density (N; tree∙plot⁻¹) were calculated. Table 1 shows the data descriptive statistics.

2.3. Data Analysis and Statistical Processing

Ten generalized models used in similar studies were selected (Table 2). Seventy percent out of the 1261 and 1751 pairs of H-d data from Pinus montezumae and Pinus pseudostrobus, respectively, were randomly selected for statistical fitting. The remaining thirty percent of the data was used to validate the fitted models.

In a first approach, the equations in Table 2 were fitted through the NLS, using the “nls” function in the R^® software (R Core Team, 2021) , in order to select a base model into which the random effects would be subsequently included. The most appropriate model was selected evaluating the statistical fitting and assuming observations independence (Temesgen et al., 2008) . Meanwhile, the deviation in the estimates was determined based on the highest values of the coefficient of determination (R², Equation (11)), the lowest values of root-mean-square error (RMSE, Equation (129), and the Akaike Information Criterion (AIC, Equation (13)) (Crecente-Campo et al., 2010; Corral et al., 2019; Bronisz & Mehtätalo, 2020) . In addition, we verified the parameters significance (p > 0.05), and the parsimony of each structure was also considered (García-Cuevas et al., 2013; Santiago-García et al., 2020) .

Table 1. Descriptive statistics of the tree and plot variables used in the research.

d = diameter at breast height (cm). H = total height (m). md = mean diameter at breast height of the plot (cm). mH = mean height of the plot. dH = dominant height of the plot (m). Msd = mean square diameter (cm). Ba = basal area per plot (m²). N = number of trees per plot.

Table 2. Generalized height-diameter equations (H-d) selected for the fitting.

ID = model identifier; a₀, a₁ and a₂ = parameters to be estimated; e = exponential function; d = diameter at breast height (cm); H = total height (m); md = mean diameter at breast height of the plot (cm); mH = mean height of the plot; dH = dominant height of the plot (m); Msd = mean square diameter (cm); Ba = basal area per plot (m²); and N = number of trees per plot.

$R^2=\frac{{\displaystyle {\sum }_{i=1}^n{\left({\hat{y}}_i-{\bar{y}}_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^n{\left(y_i-{\bar{y}}_i\right)}^2}}$ (11)

$RMSE=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{\left(y_i-{\hat{y}}_l\right)}^2}}{n-1}}$ (12)

$AIC=2\cdot p+n\cdot \ln \left(\frac{{\displaystyle {\sum }_{i=1}^n{\left(y_i-{\hat{y}}_l\right)}^2}}{n}\right)$ (13)

where $y_i$ , ${\hat{y}}_l$ and ${\bar{y}}_i$ are the observed, estimated, and mean values, respectively; n is the total number of data used in model fitting (Table 2); p refers to the number of parameters; and ln is the natural logarithm.

Once the generalized H-d model was chosen from the models fitted, a classification variable corresponding to the plot was included additively in the model, using the MEM approach (for the MEM fit). This analysis would enable the inclusion of fixed effects and optional random effects; the former applies to the whole population and the latter are specific for each grouping level. In Addition, it reduces the random variability that NLS fails to explain, it groups the H response per classification level, and it reduces the estimation bias by eliminating the initial conditions that influence the growth dynamics within each plot (Calama & Montero, 2004; Mehtätalo, de-Miguel, & Gregoire, 2015; Ferraz et al., 2018) .

The structure of the nonlinear models was: $Y_{ij}=f\left(X_{ij},{\theta }_{ij}\right)+{\epsilon }_{ij}$ ; where $f$ is the nonlinear function (selected from Table 2), $Y_{ij}$ and $X_{ij}$ are the ith-dependent and -independent observation, taken from the ith classification unit (plot), and ${\theta }_{ij}$ is the $r\times 1$ parameter vector—where r corresponds to the number of parameters of the model—and it is specific to the j-th classification level. Furthermore, this vector can be divided for the fixed and random parameters defined as ${\theta }_{ij}=A_i\gamma +B_i{b}_i$ . $A_i$ and $B_i$ matches the $r\times p$ and $r\times q$ size matrices, for the fixed and random effects, respectively; these effects are specific to each level. Meanwhile, and $b_i$ match the $p\times 1$ and $q\times 1$ vector of the fixed and random parameters, respectively (Baty et al., 2015; Corral et al., 2019) .

The MEM was fitted using the “nlme” function of R^® (R Core Team, 2021) and the marginal approximation of maximum likelihood of the empirical best linear unbiased prediction (EBLUP) (Mehtätalo & Lappi, 2020) . In addition, three variance functions were evaluated to correct for heteroscedasticity of the residuals and stabilize the variance of the inconsistent error in the estimates: 1) power function (varPower, Equation (14)); 2) exponential function (varExp, Equation (15)); and 3) constant and power function (varConstPower, Equation (6)) (Pinheiro & Bates, 2000; Zuur, Ieno, Walker, Saveliev, & Smith, 2009; Gałecki & Burzykowski, 2013; Mehtätalo & Lappi, 2020) .

$\mathrm{var}\left({\epsilon }_{ij}\right)=v_i^{2\cdot {\delta }_1}$ (14)

$\mathrm{var}\left({\epsilon }_{ij}\right)={\exp }^{2\cdot {\delta }_1\cdot v_i}$ (15)

$\mathrm{var}\left({\epsilon }_{ij}\right)={\left({\delta }_1+{\left|v_i\right|}^{{\delta }_2}\right)}^2$ (16)

where $\mathrm{var}\left({\epsilon }_{ij}\right)$ is the variance function evaluated in the variance covariate of the predictor residuals ( $v_i$ ), while ${\delta }_1$ and ${\delta }_2$ are the coefficients of the variance function, which will be specific to each level ( ${\delta }_i$ ).

The best structure was selected through the evaluation of the statistical fit with R², RMSE, and AIC values, as well as the significance of the parameters (p > 0.05) (Corral et al., 2019; Bronisz & Mehtätalo, 2020) . Similarly, the likelihood ratio of each fit was verified through the ANOVA test included in the nlme function (Pinheiro & Bates, 2000) . In the event of a tie or equality, priority will be given to the analysis that presents the smallest deviation (RMSE) and that provides a greater explanation of the sample variability (R²), since an integral evaluation between both mentioned criteria as a system of adjustment evaluation statistical is the best alternative (Mayer & Butler, 1993; Sakici et al., 2008) . The estimation bias (Bias) Equation (17) for the best structure was calculated in order to verify the quantitative deviations resulting from the model application.

$Bias={\displaystyle {\sum }_{i=1}^n\left(\frac{y_i-{\hat{y}}_l}{n}\right)}$ (17)

In addition, a means comparison (as independent populations) was carried out with a t test at a = 0.05 (Infante & Zarate, 2012) . For this purpose, 30 % of the randomly selected sample (which was not used for model calibration) was used (i.e., 378 and 525 pairs of H-d data for Pinus montezumae and Pinus pseudostrobus, respectively). The null hypothesis (H₀) tested was that there is no difference between heights of both independent populations $\left({\mu }_1={\mu }_2\right)$ , while the alternative hypothesis (H₁) was that the real value of the population mean (H) is different from the value established by H₀ $\left({\mu }_1\ne {\mu }_2\right)$ .

Further, a subsample ( $m_i$ ) of the H-d data was used to calibrate the MEM and to estimate the value of random parameters of the vector $b_i$ , (Vonesh & Chinchilli, 1997; Corral-Rivas et al., 2014; Corral et al., 2019) ; data from Equation 9 and Equation 12 independent sampling plots were used for Pinus montezumae and Pinus pseudostrobus, respectively. The equation used was the following:

${\hat{b}}_l=\hat{D}{\hat{Z}}_l^T{\left({\hat{R}}_l\hat{D}{\hat{Z}}_l^T\right)}^{-1}{\hat{\epsilon }}_l$ (18)

where $\hat{D}$ common variance-covariance matrix of the random parameters for which $q\times q$ refers to the number of parameters with the mixed effect, $Z_i$ : $m\times q$ matrix of the values of the partial derivatives, corresponding to the dimension of the trees selected for the calibration according to the chosen criterion, and evaluated as ${\hat{b}}_l=0$ , ${\hat{R}}_l$ : variance-covariance matrix of the global error of the fitted model $m_j\times m_j$ ; ${\hat{\epsilon }}_l$ : error vector of the fixed model, which represents the residuals between the observed values minus the predicted ones $m\times 1$ , and $Z^T$ ; refers to the transpose matrix of each expression.

Four calibration options (C) or localization of the MEM were analyzed to obtain the best estimation strategy of the random parameters for each sampling unit (plot), and to adjust the response of the model fixed component to each specific growth condition in the FPs (Sharma et al., 2016) . The criteria to calibrate the specific response per sampling unit of the MEM and obtain a value with which the random parameter ( $+{\phi }_i$ ) was fit in an additive way to estimate the height of the trees based on the normal diameter, were the following:

C1: Measure the heights of the two thickest trees in normal diameter (d) within the plot.

C2: Measure the heights of the two thinnest trees in d within each plot.

C3: Measure the height of the thickest and thinnest trees in d per plot.

C4: Measure the heights of the three thickest trees in d within each plot.

Once the calibration was completed, tree heights on each plot was estimated by strategy (C1, C2, C3 and C4); the estimates were evaluated through the RMSE (Equation 12), Bias (Equation 17) and the mean absolute percentage error values (E%AM: Equation19) (Mayer & Butler, 1993) .

$E\text{\%}AM=100\cdot \left[{\displaystyle \sum \left(y_i-{\hat{y}}_l/y_i\right)}\right]/n$ (19)

3. Results

The two-parameter models (Equation 1 and Equation 2) had the lowest precision and provided the less accurate explanation. The Wang and Tang model (Equation 9) was chosen as the best owing to it goodness-of-fit and parsimony, being significant all of its parameters, which is not the case with Equations 3 and 10. Therefore, this expression was used for the analysis with MEM (Table 3).

The random effect at plot level ( ${\phi }_i$ ) was additively included within the generalized H-d structure in the Wang and Tang model (Equation 9.1, 9.2, and 9.3); however, the fit resulted in non-significant parameters when including ( ${\phi }_i$ ) in the ( $a_o$ ) parameter or jointly with ( $a_1$ ) and ( $a_2$ ); therefore, those parameters are not shown.

$H=1.3+a_{0}d{H}^{\left(\left(a_1+{\phi }_i\right){\text{e}}^{\frac{-a_2}{d}}\right)}$ (9.1)

$H=1.3+a_{0}d{H}^{\left(a_1{\text{e}}^{\frac{-\left(a_2+{\phi }_i\right)}{d}}\right)}$ (9.2)

$H=1.3+a_{0}d{H}^{\left(\left(a_1+{\phi }_i\right){\text{e}}^{\frac{-\left(a_2+{\phi }_i\right)}{d}}\right)}$ (9.3)

The inclusion of an additive random effect ( $+{\phi }_i$ ) in the parameter $a_1$ for the MEM within the Wang and Tang model, and in the combination of $a_1$ and $a_2$ provided good results. These results indicated that the inclusion with the best qualification is given in the combination of these parameters, and that a statistical gain in the fit is obtained with respect to NLS in R², RMSE, and AIC: 4.3%, 10.9%, and 3.2%, for Pinus montezumae and 4.9%, 24.8%, and 6.5% for Pinus

Table 3. Goodness-of-fit statistics for the generalized height-diameter models.

ID = model identifier; R² = Coefficient of determination; RMSE = Root-mean-square error; and AIC = Akaike Information Criterion.

pseudostrobus, respectively. The inclusion of the varExp structure in the MEMs improved the results for both species since it regulated the unequal error variance (Table 4).

According to the ANOVA, the Wang and Tang expression that included effects on the parameters $a_1$ and $a_2$ , and the varExp heteroscedasticity correction structure was the one that showed better results for Pinus montezumae (Table 5). The test indicated that reasonable fit is obtained for Pinus pseudostrobus when including the varPower and varConstPower variance structures (Table 5). However, for Pinus montezumae greater deviations of the response variable estimations, less explanatory power (Table 4), and non-significant parameters were obtained; hence, we decided to select the setting (Equation 9.2 + Equation 15: varExp) for both species, with the inclusion of random effects in the parameters $a_1$ and $a_2$ .

A bell-shaped curve (i.e. Gaussian) of the residuals was obtained for both species (Figure 2(a) and Figure 2(b)) indicating compliance of this regression assumption, which is a desirable pattern. The analysis of the homoscedasticity and comparison of the fittings made (NLS and MEM) showed a more compact distribution (close to zero) when using the MEM approach; this is due to the fact that the increasing trend of residuals as the predicted value increases is corrected, and they are held constant when the structure that models the varExp type variance is included in the fit (Figure 2(c), Figure 2(d)).

An individual deviation of −0.1110 m for Pinus montezumae and 0.1560 m for Pinus pseudostrobus was observed for the H estimates calculated with the fixed parameters model, obtained with the MEMs approach (Table 6). In addition, the evaluation by diameter class indicated deviations lower than the unity in all cases, as well as the largest bias in those categories with limited data (Figure 3).

Table 4. Statistical indicators of the goodness of fit for the generalized H-d models under the mixed effects model.

where, ID = model variant identifier; R² = coefficient of determination; RMSE = root-mean-square error; and AIC = Akaike Information Criterion.

Table 5. Comparison of height-normal diameter (H-d) generalized full models in its nested form under the mixed effects approach.

where ID = model variant identifier; BIC = Bayesian Information Criterion; logLik: log likelihood.

Figure 2. Frequency and distribution of residuals for the Wang and Tang generalized model according to its statistical fit—Non-linear Least Squares (NLS) and Mixed Effects Models (MEM)—for the two pine species analyzed.

Table 6. Fixed parameters with the mixed-effects-by-species model approach for the generalized Wang and Tang model.

Figure 3. Bias by diameter category (Cd) when applying the generalized Wang and Tang model for the two Pinus species.

Table 7 shows the variance-covariance matrix (vcov) values; the diagonal terms are the variances and the non-diagonal ones are covariances. Also, parameter values of the function included to model the residual distribution is presented (weights = varPower (value = 0.5)), the standard error given for each one by group (intervals (model, level = 0.95)), and residual values of the fit used to calibrate the random effects obtained with the nlme function used to fit the MEM.

The t test between the validation subsample data (i.e. observed) and the estimates indicated equality between populations for both species (Pinus montezumae: t = −0.8354 and p = 0.4038; Pinus pseudostrobus: t = −0.0792 and p = 0.9369). Therefore, H_a: $\left({\mu }_1={\mu }_2\right)$ is accepted. The graphical representation confirms these results (Figure 4).

The evaluation of predictions with the calibration alternatives shows that the minor global deviations for Pinus montezumae are obtained by including the heights of two trees, one with the largest and another with the smallest d in the plot, whereas for Pinus pseudostrobus the heights of three trees with the largest d in the plot should be used to calibrate the value of the random parameter (Table 8).

Table 7. Values of the variance-covariance matrix, value of the homoscedasticity function incorporated, within -group standard error and residual value of random effects.

Figure 4. Trend of height estimates as a function of diameter at breast height and mean dominant height when the Wang and Tang generalized model is used.

Table 8. Goodness-of-fit statistics for the calibration alternatives of the Wang and Tang generalized model to estimate total height as a function of normal diameter.

where, C1: Measure the heights of the two thickest trees in normal diameter (d) in the plot. C2: Measure the heights of the two thinnest trees in d in the plot. C3: Measure the height of the thickest and thinnest tree in d in the plot. C4: Measure the heights of the three thickest trees in d in the plot. RMSE: root mean square error. E%AM: mean absolute percentage error.

4. Discussion

Bronisz & Mehtätalo (2020) and Liu et al. (2017) pointed out that two-parameter models have some accuracy disadvantages with regard to three-parameters models, hence, in order to estimate tree height based on the diameter at breast height, choosing the latter is recommended. The Wang and Tang model proposed here differs from structures selected by other authors. For example, Misir (2010) proposed the Schnute model which includes dominant diameter (Dd) as predictor for Populus tremula L., whereas Bronisz & Mehtätalo (2020) projected a Schumacher-type expression with Msd and Ba as predictors for Eucalyptus globulus L.; in our case, the equations in which Msd was included (Equations 1, 6, 7, 8, and 10) showed non-significant parameters for both species. Conversely, our results concurred with the findings of Liu et al. (2017) for Metasequoia glyptostroboides Hu & Cheng; these authors concluded that including tree density as explanatory variables did not contribute significantly to models improvement after fitting 16 generalized models.

The inclusion of an additive effect ( $+a_i$ ) in the Wang and Tang model fit with MEM is similar to the procedure followed by Sharma & Parton (2007) , Vargas-Larreta et al. (2009) , Corral et al. (2014) , and Corral et al. (2019) . These authors also recorded a statistical gain in the model’s fit with respect to NLS in terms of the R², RMSE and AIC values. Our results are also similar to those reported by Crecente-Campo et al. (2010) who used the Bertalanffy-Richards model to estimate H as a function of d and dH, obtaining in average 6.4% and 18.0% gains for the R² and RMSE values, respectively. Therefore, the use of MEM with the inclusion of random effects makes it possible to model the random variability that results from different growth conditions of each site; in addition, more efficient and precise predictors are obtained, which is reflected in the goodness-of-fit statistics (Crecente-Campo et al., 2010; Ercanli, 2015) .

The R² and RMSE values were higher than those reported by Santiago-García et al. (2020) who used a power function to weight the variance for Pinus patula Schiede ex Schltdl. et Cham. (R² = 0.75, RMSE = 4.3), Pinus oaxacana Mirov (R² = 0.82, RMSE = 4.3), Pinus ayacahuite Ehren. (R² = 0.83, RMSE = 3.5), Pinus teocote Schlecht. & Cham. (R² = 0.80, RMSE = 3.6), and Pinus leiophylla Schiede & Deppe (R² = 0.83, RMSE = 3.4). They used the Sharma and Parton expression for the first three species and the Wang and Tang and Nilson for the rest. These differences might be due to the homogeneity of the plantations with respect to natural forests conditions; unfortunately, to our best knowledge, no work has been developed on forest plantations in Mexico in this subject.

In addition, our results also agree with Álvarez-González et al. (2007) findings, who concluded that including a structure to model the variance within the statistical fitting helps to correct heteroscedasticity problems in residuals distribution. Similarly, our findings are consistent with Ercanli (2015) and Corral et al. (2019) results, who reported that the inclusion of random effects in two parameters within the Schnute and Chapman-Richards generalized model improves the fit and precision indicators.

Our results of the mixed effects models obtained for both species (Equation 9.2 + Equation 15: varExp) showed an efficient correction of heteroscedasticity for the Wang and Tang generalized model, in comparison with both NLS and MEM fit without a variance structure (varExp), and this helped to reduce the deviations in heights (H) estimation (better RMSE) and to maximize the explanation of the sample variability (R²) (Bronisz & Mehtätalo, 2020; Santiago-García et al., 2020) . Furthermore, a lower estimation error was achieved with the same field sampling effort by including the dominant height as inherent to the site (Corral et al., 2019) . Likewise, the biases calculated for both species are similar to those reported in similar works; for example, Corral et al. (2014, 2019) reported bias values of −0.013 m and −0.009 m, in seven pine species in Durango, Mexico for a generalized model that considered d, dH, density, and Msd as predictors of H.

The results also showed that by taking into consideration variables inherent to the population (mean square diameter: Dq, mH, A₀, Ba, or age) to explain the growth trend of H in different stands is an alternative to improve the quality of the estimates, and increases the applicability of the expressions proposed on the H-d local equations (Misir, 2010; Liu et al., 2017) . The only differences were found in the explanatory variable(s) of the forest stand or FP, given the particular ecological conditions and requirements of each species studied.

Moreover, including dH as an explanatory variable helps to reduce bias in the estimates, and this is one of the variables with the greatest contribution to explain the sample variability and the dynamics of H in the stand or FP. The dH is directly related to site quality, it is sensitive to the environmental conditions in which the species grows, it remains constant after thinning (as long as it is not on the upper part), and it has shown better results than mH; besides, its registration does not entail a greater effort in the field (Crecente-Campo et al., 2010; Ferraz et al., 2018; Santiago-García et al., 2020) .

The model proposed can be a reliable quantitative tool for the development of forest management plans, owing to its statistical robustness and precision. The model is included in the growth and yield modeling, when it is combined with a site index equation or when it is used as a support in forest inventories (Ercanli, 2015; Santiago-García et al., 2020) . Furthermore, a dominant height by site index tag structure can be used within the structure proposed for the Wang and Tang model, which could help to expand its applicability and allow application of the results outside the research area or to other plantations of these two species, due to the direct relationship among site quality, tree age and dominant height within the site or plantation (Spurr & Barnes, 1982; Prodan et al., 1997; Torres & Magaña, 2001) .

The calibration of a MEM with a subsample of tree that does not generate additional sampling costs, allows the estimation of height of all the trees within a sampling plot, based on the d of an independent subsample of trees, and gives insight into it application in different areas. In addition, this is an alternative approach that allows for saving time and costs in the planning and execution of forest management programs (Corral et al., 2019) . Even though there are several calibration criteria for this type of fit (Calama & Montero, 2004; Crecente-Campo et al., 2010; Corral-Rivas et al., 2014) , several authors agree that the use of few trees (<4 individuals) yields specific results for each site or plot condition, reduces the deviations in the estimates and improves the predictive capacity of the equation proposed.

5. Conclusion

We obtained a significant statistical improvement with the use of mixed effects models (MEM) when fitting a generalized height-diameter (H-d) model, over nonlinear least squares (NLS). In addition, when a function to model the variance structure is included in the MEM, the residuals decreased and this helped to correct undesirable residuals trend.

The goodness of fit statistics, bias, means comparison test, and distribution of the estimates suggested that the Wang and Tang model that uses d and dominant height can accurately estimate H for Pinus montezumae and Pinus pseudostrobus forest plantations in Michoacán, Mexico.

The calibration or localization process of the MEM also allows the application of the expression proposed in plantations of these two species that were not included in the analysis, reduces the estimation error and improves the estimates of the H-d relationship. Finally, this does not represent an additional sampling effort or time invested in the inventory.

Acknowledgements

The first auhor would like to sincerely thank CONAHCYT for financially supporting this work. The authors gratefully acknowledge the COLPOS and INIFAP for their kind cooperation, and guidance during the study. We would also like to thank the anonymous reviewers for their edifying comments which have improved the quality of this manuscript.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Alder, D. (1980). Estimación del Volumen Forestal y Predicción del Rendimiento con Referencia Especial a los Trópicos (Vol. 2, 118 p.). Organización de las Naciones Unidas para la Agricultura y la Alimentación.
[2]	Álvarez-González, J. G., Rodríguez-Soalleiro, R., & Rojo-Alboreca, A. (2007). Resolución de problemas del ajuste simultaneo de sistemas de ecuaciones: Heterocedasticidad y variables dependientes con distinto número de observaciones. Cuadernos de la Sociedad Española de Ciencias Forestales, 23, 35-42.
[3]	Barrio, A. M., Álvarez, J. G., & Díaz-Maroto, H. (2004). Elaboración de una tarifa con clasificación de productos para Quercus robur L. en Galicia, basada en un modelo de volumen porcentual. Investigación agraria. Sistemas y Recursos Forestales, 13, 506-517. https://doi.org/10.5424/srf/2004133-00849
[4]	Baty F., Ritz, C., Charles, S., Brutsche, M., Flandrois, J. P., & Delignette-Muller, M. L. (2015). A Toolbox for Nonlinear Regression in R: The Package Nlstools. Journal of Statistical Software, 66, 1-21. https://doi.org/10.18637/jss.v066.i05
[5]	Bronisz, K., & Mehtätalo, L. (2020). Mixed-Effects Generalized Height-Diameter Model for Young Silver Birch Stands on Post-Agricultural Lands. Forest Ecology and Management, 460, Article ID: 117901. https://doi.org/10.1016/j.foreco.2020.117901
[6]	Calama, R., & Montero, G. (2004). Interregional Nonlinear Height-Diameter Model with Random Coefficients for Stone Pine in Spain. Canadian Journal of Forest Research, 34, 150-163. https://doi.org/10.1139/x03-199
[7]	CONAFOR (2018). Programas específicos de intervención institucional: Programa de plantaciones forestales comerciales 2014-2018. https://www.gob.mx/conafor/documentos/plantaciones-forestales-comerciales-27940
[8]	Corral, R. S., Silva, A. M., & Quiñonez, G. B. (2019). Modelo generalizado no-lineal altura-diámetro con efectos mixtos para siete especies de Pinus en Durango, México. Revista Mexicana de Ciencias Forestales, 10, 86-117. https://doi.org/10.29298/rmcf.v10i53.500
[9]	Corral-Rivas, S., Álvarez-González, J. G., Crecente-Campo, F., & Corral-Rivas, J. J. (2014). Local and Generalized Height-Diameter Models with Random Parameters for Mixed, Uneven-Aged Forests in Northwestern Durango, Mexico. Forest Ecosystems, 1, 1-9. https://doi.org/10.1186/2197-5620-1-6
[10]	Correia, V. G., Ribeiro, A. M., Fernandes, G. S., Zanetti, S. S., Marques, M. S., & Rosa, A. S. (2018). Prognoses of Diameter and Height of Trees of Eucalyptus Using Artificial Intelligence. Science of the Total Environment, 619-620, 1473-1481. https://doi.org/10.1016/j.scitotenv.2017.11.138
[11]	Crecente-Campo, C. F., Tomé, M., Soares, P., & Diéguez-Aranda, U. (2010). A Generalized Nonlinear Mixed-Effects Height-Diameter Model for Eucalyptus globulus L. in Northwestern Spain. Forest Ecology and Management, 259, 943-952. https://doi.org/10.1016/j.foreco.2009.11.036
[12]	Diéguez-Aranda, U., Barrio, A. M., Castedo-Dorado, F., & Álvarez, G. (2005). Relación altura-diámetro generalizada para masas de Pinus sylvestris L. procedentes de repoblación en el noroeste de España. Investigación agraria. Sistemas y recursos forestales, 14, 229-241. https://doi.org/10.5424/srf/2005142-00886
[13]	Ercanli, Í. (2015). Nonlinear Mixed Effect Models for Predicting Relationships between Total Height and Diameter of Oriental Beech Trees in Kestel, Turkey. Revista Chapingo Serie Ciencias Forestales y del Ambiente, 21, 185-202. https://doi.org/10.5154/r.rchscfa.2015.02.006
[14]	Ferraz, F. A. C., Mola-Yudego, B., Ribeiro, A., Scolforo, R. S. J., Loos, R. A., & Scolforo, H. F. (2018). Height-Diameter Models for Eucalyptus sp. Plantations in Brazil. CERNE, 24, 9-17. https://doi.org/10.1590/01047760201824012466
[15]	Flores-Ayala, E., Pineda-Ojeda, T., Hernández-Ramos, J., Buendía-Rodríguez, E., Flores, A., & Guerra-de la Cruz, V. (2023). Models for Determining Allometric Relationships in Juniperus deppeana Steud. in the State of Tlaxcala. Revista Mexicana de Ciencias Forestales, 14, 30-53. https://doi.org/10.29298/rmcf.v14i80.1363
[16]	Gałecki, A., & Burzykowski, T. (2013). Linear Mixed-Effects Models Using R (574 p.). Springer Science-Business Media. https://doi.org/10.1007/978-1-4614-3900-4
[17]	García-Cuevas, X., Hernández-Ramos, J., García-Magaña, J. J., Aguilar-Toral, C., & Hernández-Ramos, A. (2013). Ecuaciones generalizadas de altura-diámetro para rodales de Pinus montezumae en Nuevo San Juan Parangaricutiro, Michoacán, México. In M. P. Barrón-González, & S. Moreno-Limón (Eds.), Aporte científico ante el cambio climático y el desarrollo sostenible: Sociedad de cambio climático y desarrollo sostenible (pp. 9-24). Editorial Académica Española.
[18]	Gómez-García, E., Fonseca, T. F., Crecente-Campo, F., Almeida, L. R., Diéguez-Aranda, U., Huang, S., & Marques, C. P. (2015). Height-Diameter Models for Maritime Pine in Portugal: A Comparison of Basic, Generalized and Mixed-Effects Models. iForest, 9, 72-78. https://doi.org/10.3832/ifor1520-008
[19]	González, C. E., Gasca, E. M., & Heredia, P. D. (2014). Cultura organizacional del sistema empresarial de la Comunidad Indígena de Nuevo San Juan Parangaricutiro: Un manejo sustentable forestal. Revista Mexicana de Agronegocios, 35, 1023-1034.
[20]	Hernández, R. J., Avilés, A. C., García, J. J. M., Hernández, A. R., García, X. C., & Flores, C. L. (2020). Ecuaciones locales y generalizadas de altura-diámetro para Pinus patula Schl. et Cham. en Veracruz, México. Ecosistemas y Recursos Agropecuarios, 7, e2457.
[21]	Hernández, R. J., García, J. J. M., García, X. C., García, G. G. E., Hernández, A. R., Muñoz, H. J. F., & Martínez, M. S. (2018). Ecuaciones generalizadas altura-diámetro para bosques de Pinus pseudostrobus Lindl. en Michoacán, México. Madera y Bosques, 24, e242494. https://doi.org/10.21829/myb.2018.242494
[22]	Hernández, R. J., García, X. C., Hernández, A. R., García, J. J. M., Muñoz, H. J. F., Flores, C. L. G. G., & García, E. (2015). Ecuaciones altura-diámetro generalizadas para Pinus teocote Schlecht. & Cham. en el estado de Hidalgo. Revista Mexicana de Ciencias Forestales, 6, 8-21. https://doi.org/10.29298/rmcf.v6i31.192
[23]	INEGI (2017). Anuario estadístico y geográfico de Michoacán de Ocampo 2017 (723 p.). INEGI. Aguascalientes, México.
[24]	Infante, G. S., & Zárate, G. P. L. (2012). Métodos estadísticos: Un enfoque interdisciplinario (3ª ed., 624 p.). Colegio de Postgraduados. Montecillo, Edo. de Méx., México.
[25]	Liu, M., Feng, Z., Zhang, Z., Ma, C., Wang, M., Lian, B., Sun, R., & Zhang, L. (2017). Development and Evaluation of Height Diameter at Breast Models for Native Chinese Metasequoia. PLOS ONE, 12, e0182170. https://doi.org/10.1371/journal.pone.0182170
[26]	Mayer, D. G., & Butler, D. G. (1993). Statistical Validation. Ecological Modelling, 68, 21-32. https://doi.org/10.1016/0304-3800(93)90105-2
[27]	Mehtätalo, L. (2005). Height-Diameter Models for Scots Pine and Birch in Finland. Silva Fennica, 39, 55-66. https://doi.org/10.14214/sf.395
[28]	Mehtätalo, L., & Lappi, J. (2020). Biometry for Forestry and Environmental Data with Examples in R (411 p.). Taylor & Francis Group, LLC. https://doi.org/10.1201/9780429173462
[29]	Mehtätalo, L., de-Miguel, S., & Gregoire, T. G. (2015). Modeling Height-Diameter Curves for Prediction. Canadian Journal of Forest Research, 45, 826-837. https://doi.org/10.1139/cjfr-2015-0054
[30]	Mısır, N. (2010). Generalized Height-Diameter Models for Populus tremula L. Stands. African Journal of Biotechnology, 9, 4348-4355.
[31]	Muñoz, F. H. J., Sáenz, J. T. R., García, J. J. M., Coria, V. M. Á., & Muñoz, Y. V. (2015). Áreas potenciales para establecer plantaciones comerciales de pino en la sierra Purhépecha, Michoacán. Foresta Veracruzana, 17, 35-42.
[32]	Muñoz, F., Sáenz, J. T. R., García, J. J. S., Hernández, E. M., & Anguiano, J. C. (2011). Áreas potenciales para establecer plantaciones forestales comerciales de Pinus pseudostrobus Lindl. y Pinus greggii Engelm. en Michoacán. Revista Mexicana de Ciencias Forestales, 2, 19-44.
[33]	Pinheiro, J., & Bates, D. (2000). Mixed-Effects Models in S and S-PLUS. Statistics and Computing (528 p.). Springer-Verlag.
[34]	Prodan, M., Peters, R., Cox, F., & Real, P. (1997). Mensura Forestal. Serie Investigación y Educación de Desarrollo Sostenible (586 p.). Instituto Interamericano de Cooperación para la Agricultura (IICA)/BMZ/GTZ sobre agricultura, recursos naturales y desarrollo sostenible, San José.
[35]	R Core Team (2021). R: A Language and Environmental for Statistical Computing. R Foundation for Statistical Computing. http://www.R-project.org/
[36]	Sakici, O. E., Misir, N., Yavuz, H., & Misir, M. (2008). Stem Taper Functions for Abies nordmanniana subsp. bornmulleriana in Turkey. Scandinavian Journal of Forest Research, 23, 522-533. https://doi.org/10.1080/02827580802552453
[37]	Sánchez-González, M., Cañellas, I., & Montero, G. (2007). Generalized Height-Diameter and Crown Diameter Prediction Models for Cork Oak Forests in Spain. Investigación Agraria: Sistemas y Recursos Forestales, 16, 76-88. https://doi.org/10.5424/srf/2007161-00999
[38]	Santiago-García, W., Jacinto-Salinas, A. H., Rodríguez-Ortiz, G., Nava-Nava, A., Santiago-García, E.,Ángeles-Pérez, G., & Enríquez-del Valle, J. R. (2020). Generalized Height-Diameter Models for Five Pine Species at Southern Mexico. Forest Science and Technology, 16, 49-55. https://doi.org/10.1080/21580103.2020.1746696
[39]	Santiago-García, W., Pérez-López, E., Quiñonez-Barraza, G., Rodríguez-Ortiz, G., Santiago-García, E., Ruiz-Aquino, F., & Tamarit-Urias, J. C. (2017). A Dynamic System of Growth and Yield Equations for Pinus patula. Forests, 8, Article No. 465. https://doi.org/10.3390/f8120465
[40]	Sharma, M., & Parton, J. (2007). Height-Diameter Equations for Boreal Tree Species in Ontario Using Mixed-Effects Modeling Approach. Forest Ecology and Management, 249, 187-198. https://doi.org/10.1016/j.foreco.2007.05.006
[41]	Sharma, R. P., Vacek, Z., & Vacek, S. (2016). Nonlinear Mixed Effect Height-Diameter Model for Mixed Species Forests in the Central Part of the Czech Republic. Journal of Forest Science, 62, 470-484. https://doi.org/10.17221/41/2016-JFS
[42]	Spurr, S. H., & Barnes, B. V. (1982). Ecología Forestal (Trad. C. L. Raigorodsky Z., 690 p.). A. G. T. Editor, México, D. F.
[43]	Temesgen, H., Monleon, V. J., & Hann, D. W. (2008). Analysis and Comparison of Nonlinear Tree Height Prediction Strategies for Douglas-Fir Forests. Canadian Journal of Forest Research, 38, 553-565. https://doi.org/10.1139/X07-104
[44]	Torres, R. J. M., & Magaña, S. O. (2001). Evaluación de Plantaciones Forestales (472 p.). Limusa Noriega Editores, México, D. F.
[45]	Vargas-Larreta, B., Castedo-Dorado, F., Álvarez-González, J. G., Barrio-Anta, M., & Cruz-Cobos, F. (2009). A Generalized Height-Diameter Model with Random Coefficients for Uneven-Aged Stands in El Salto. Durango (Mexico). Forestry: An International Journal of Forest Research, 82, 445-462. https://doi.org/10.1093/forestry/cpp016
[46]	Vonesh, E. F., & Chinchilli, V. M. (1997). Linear and Nonlinear Models for the Analysis of Repeated Measurements (560 p.). Marcel Dekker Inc. https://doi.org/10.1201/9781482293272
[47]	Wang, M., & Tang, S. (1997). Research on Universal Height Diameter Curves. Forest Research, 10, 259-264.
[48]	Zambrano, T., Suarez, M., & Jerez, M. (2001). Evaluación de la efectividad del ajuste de algunos modelos de regresión utilizados para estimar la relación altura-diámetro en parcelas permanentes de rendimiento y aclareo en teca (Tectona grandis L. F.). Revista Forestal Venezolana, 45, 163-173.
[49]	Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A., & Smith, G. M. (2009). Mixed Effects Models and Extensions in Ecology with R (574 p.). Springer Science Business Media. https://doi.org/10.1007/978-0-387-87458-6