Study area

The study was conducted in mixed forests of the San Diego de Tezains and Salto de Camellones ejidos which belong to the Forestry Management Unit (UMAFOR) 1005 "Santiago Papasquiaro y Anexos". The study area is located in the northwest region of the state of Durango, Mexico, between 24°30'-25°27' N and 105°01'-106°24' W. Height above sea level fluctuates between 800 and 3 103 m. The prevailing climate is temperate with rain in summer. The average annual rainfall varies from 1 000 to 1 200 mm. The average annual temperature fluctuates from 5 to 18 °C, the coldest temperature occurs in January (-6 °C) and the hottest in May (28 °C) (^{García, 1981}; ^{Quiñonez et al., 2018}). The most important forest communities are made up of species of Pinus, Quercus, Juniperus, Cupressus, Pseudotsuga, Arbutus and Alnus. Approximately 90 % of the study area consists of mixed stands with irregular structures (^{Quiñonez et al., 2018}).

Source of the data

The mensuration information used to adjust the functional relationship h-d was obtained from 44 units or permanent sampling sites and monitoring of growth and production, established in 2008, covering all types of vegetation, site qualities, densities and ages in the stands in timber production. The dimensions of the sites are 50 × 50 m and they were established under a systematic sampling design. For each tree with a normal diameter equal to or greater than 7.5 cm, information was recorded on species, normal diameter (d, cm), total height of the tree (h, m), crown diameter (m) and its position with respect to the center of the plot (azimuth, degrees, distance, m).

A database with a total of 24 species was processed. The estimated composition was 84 % for the genus Pinus (P. arizonica Engelm., P. ayacahuite Ehrenb., P. durangensis Martínez, P. herrerae Martínez, P. lumholtzii B.L. Rob. & Fernald, P. teocote Schiede ex Schltdl. & Cham., P. douglasiana Martínez), 8 % for the genus Quercus (Q. sideroxyla Humb. & Bonpl., Q. arizonica Sarg., Q. mcvaughii Spellenb, Q. durifolia Seemen ex Loes, Q. crassifolia Bonpl., Q. jonesii Trel., Q. rugosa Neé, Q. laeta Liebm.) and 8 % for species of the genus Juniperus, Cupresus, Arbutus and Alnus (J. deppeanna Steud., J. durangensis Martínez, C. lusitanica Mill., A. arizonica (A. Gray) Sarg., A. bicolor M. González & PD Sørensen, A. madrensis M. González, A. tesselata PD Sørensen, A. xalapensis Kunth., and Alnus firmifolia Mill). Only the h-d data of Pinus species were used. The distribution of the pairs of data was examined through graphs to detect outliers, which were eliminated from the database used for equation adjustment, the main descriptive statistics of the variables analyzed are shown in Table 1.

The following indicators were estimated for each sampling unit (site): the number of trees per hectare (N, trees ha^-1), basal area per hectare (G, m² ha^-1), mean square diameter (dg, cm), dominant height estimated as the average of the 25 trees with the largest diameter in the plot (H ₀ , m) (^{Assmann, 1970}), dominant diameter estimated as the average of the 25 trees with the largest diameter in the plot (D ₀ , cm) (^{Assmann, 1970}) and Hart Index (%) estimated by the expression:

Models

Ten local equations and 10 generalized equations were analyzed and compared to describe the h-d relationship (Table 2). The adjustments were individual, by species and in group (all species of the genus Pinus) by ordinary non-linear least squares (ONLS) with the R program (^{R Core Team, 2016}) by means of the "nlsLM" function of the "minpack.lm" package (^{Elzhov et al., 2012}).

Adjustment capacity evaluation was based on residual graphs, numerical comparisons of the root mean square error (RMSE), values of the coefficient of determination (R ² ), the average bias of the residuals (Em) and Schwarz Baysian Information Criterion (BIC) (^{Schwarz, 1978}) which compares equations in terms of their parsimony (simplicity). The formulation of these statistics is expressed as follows:

Where:

yi, yi^ y y- = Observed height, estimated height and average height

n = Number of observations

k = Number of parameters of the equation

ln = Natural logarithm

For local equations which showed the greatest advantage in goodness-of-fit-statistics, the relationship between the parameters and the site variables (covariables) was analyzed, in order to better explain the variability between total height and the characteristics of the variables. For the inclusion of covariates in the base model, different combinations of stand variables and their transformations (square, logarithm, root, inverse, etc.) were also tested in order to improve the adjustments. The relationship was analyzed through the "cor" function of R (^{R Core Team, 2016}) and graphical analysis.

Mixed-effects model

While the h-d data pairs used in fitting the equations were obtained from different sampling units, a high correlation was detected between the measurements taken within the sampling unit, which infringes one of the fundamental assumptions of the regression model. In order to deal with the lack of independence between observations, the mixed-effects model was used, in which the variability between sampling units is explained by the introduction of a random parameter that is estimated together with the fixed parameters (^{Pinheiro and Bates, 1998}; ^{Calama and Montero, 2004}; ^{Castedo et al., 2006}).

The parameters vector in a non-linear mixed model can be defined, according to ^{Pinheiro and Bates, (1998)} in the following way:

Where:

yij = The j^th total tree height taken from the i^th sample unit

xij = The j^th measurement of the predictor variable (d and site variables) taken from the i^th sample unit

ϕi = Parameter vector rx1 (where r is the total number of parameters in the model) specific for the j^th sample

f = Non-linear function

εij = Random error in vector form

Vector ϕi can be divided in fixed and random components. Fixed components are shared by the entire population, while random components vary from plot to plot, in the following way:

Where:

λ = px1 vector for the fixed parameters (p being the number of fixed parameters)

b_i = qx1 vector for the random parameters, related to the i^th unit (q is the number of random parameters)

A_i and B_i = rxp and rxq size matrixes for the fixed and random effects, respectively, corresponding to the i^th unit

The basic theoretical assumptions of the non-linear mixed models assume that the residual vector (ε_i ) and the random effect vector (bi) have a normal distribution, with mean equal to zero and a variance-covariance matrix (R_i and D) which represents the variability between the different sampling units (^{Littell et al., 2006}).

Random Parameters Estimation

The random parameters were estimated using the maximum likelihood marginal approximation of the empirical best linear unbiased predictor (EBLUP) (^{Vonesh and Chinchilli, 1997}) from the expansion of the Taylor serie, using the algorithm of ^{Lindstrom and Bates (1990)}, implemented in the "nlme" function of R (^{R Core Team, 2016}).

The construction of the generalized equation with mixed effects inclusion was based on the selection of the local equation which showed the best goodness-of-fit results of all of the species under study. Thus, different combinations of fixed and random parameters were tested in each equation, and the RMSE, R² , Em and BIC statistics were compared in order to determine which fixed parameters should be considered as mixed. The h-d equation generated in this way was compared with the best generalized equation.

Random parameters prediction

The application of the mixed effects model requires previous information on the response variable (total tree height) that is later used to calculate the specific random parameters of that sample unit. The measured heights of one or more trees at each site can be used to predict the specific random parameters for that sampling unit and adjust the fixed part of the mixed effects model (average population response) to the stand conditions. This is known as location or calibration of the mixed effects model (^{Sharma et al., 2016}).

When a sub-sample of m_j tree heights is available, it can be used to predict the random parameters of vector b_j, with the following equation (^{Vonesh and Chinchilli, 1997}):

Where:

D^ = q × q variance-covariance matrix related to the random parameters (q being the number of random parameters included in the model), which is common to all units and is estimated from the general model adjustment

R^i = m_j × m_j variance-covariace matrix of the model error component

ε^i = m × 1 residuals vector, whose components result from the difference between the observed hight of each tree and the predicted value from the model, considering only fixed parameters

Zi = m × q matrix of the partial derivatives of the valued random parameters from b^i=0

To determine the optimal size of the sub-sample of trees to be measured within each sample unit and find the response of the calibration (RC), the following options were analyzed:

The six alternatives described were evaluated in terms of RMSE, R² and Em and compared with the statistics obtained by the equations adjusted with the ONLS and NLME procedure (i.e. when all the trees in the unit are included).

Heights of trees predictions in the different calibration options evaluated were calculated with the method used to estimate the random parameters through the following expression:

Where:

h^ij = Total estimated height (m)

ϕ^i=Aiλ^ 

λ^ = Matrixes for the specific fixed and random effects of the i^th sample unit, respectively

xij = Predictor variable vector

Local and generalized equations

To select the local and generalized equation with the best individual goodness-of-fit result for each species, firstly, all equations in Table 2 were adjusted by species. During this stage, convergence problems and non-significant parameters were identified for some of them (Pinus ayacahuite, P. herrerae, P. lumholtzii and P. douglasiana) with equations (10), (13) and (17); after that, it was decided to carry out the subsequent analyzes with the joint database for all Pinus genera.

When comparing the values of goodness-of-fit statistics for the group of local equations (Table 3), Equations (3), (8), (9) and (10) were the most accurate and with the least bias. Although these four equations provided very similar values in the comparison statistics used, equation (3) presented non-significant parameters in the fit, even though the BIC value was slightly lower than the rest of the equations. Finally, the graphical bias analysis in height estimation by diametric classes indicated that equation (8) showed the best behavior in the whole data interval, so it was selected as the base model to study the inclusion of variables of stand in its structure (Figure 1a).

Figure 1 Bias behavior by diametric class of local and generalized equations with the best goodness-of-fit results. 

^{Huang et al. (2000)} and ^{Peng et al. (2001)} point out that Equation (8) is perhaps the most flexible and versatile to model the h-d relationship in forest stands. In this regard, ^{Uzoh (2017)} showed that the Equation (8) of Chapman-Richards is statistically the most robust (i.e. least biased) out of 12 different options to estimate the height of Pinus ponderosa Douglas ex Lawson in a mixed-effects model in the west of the United States of America.

While comparing the goodness-of-fit statistics values of the generalized equations group (Table 3), equations (11) and (13) stand out as the most accurate, considering R² and RMSE statistics. However, Equation (13) exhibited non-significant parameters, even though it has the lowest value in BIC results. Therefore, the generalized equation of ^{Sharma and Parton’s (2007)} (11) was that of the best advantages in precision and parsimony (BIC) of all evaluated equations. In addition to this, the graphical analysis of bias while estimating heights by diameter classes had an unbiased behavior (Figure 1b). This equation is quite flexible to describe the h-d relationship and has been widely used in similar studies (^{Vargas et al., 2009}; ^{Crecente et al., 2010}; ^{Corral et al., 2014}).

The best goodness-of-fit statistics in terms of precision and parsimony are derived from the generalized equations (Table 3), since they include stand characteristics as predictor variables, which decrease variability within each sample unit, and in general, are widely recommended to describe the h-d relationship (^{López et al., 2003}; ^{Barrio et al., 2004}; ^{Canga et al., 2007}; ^{Vargas et al., 2009}; ^{Crecente et al., 2010}; ^{Corral et al., 2014}; ^{Uzoh, 2017}).

The correlation analysis between stand variables and the parameters in Equation (8) (Table 4) show that only three variables (H₀ , N and G) were significantly correlated to improve their accuracy, resulting in the generalized equation with the following expression:

Where:

b₀ - b₄ = Parameters of the equation; the rest of the previously defined variables

Where:

H₀ = Dominant height (m)

D₀ = Dominant diameter (cm)

N = Density (number of trees ha^-1)

G = Total basal area (m² ha^-1)

d_g = Mean square diameter (cm)

IH = Hart Index (%)

b_i = Estimated parameters

With the inclusion of dominant height, basal area and stand density, in addition to normal diameter in Equation (28), RMSE was reduced by 18 % with respect to the Chapman-Richards base model (Equation 8) (from 2.95 m to 2.43 m) and, consequently, R² was increased by 12 % (from 0.731 to 0.818). BIC was reduced by 8 % (from 20 200 to 18 651) since it includes the sum total of the residual squares.

On the other hand, while comparing the precision and parsimony of Equation (28) with respect to Equation (11), developed by ^{Sharma and Parton (2007)}, a marginal improvement was also found since RMSE decreased by 0.4 % and the value of R² increased by 0.16 %, while BIC had a reduction of 1.5 %.

Based on goodness-of-fit statistics and on the graphical behavior of bias by diametric class (Figure 1b), Equation (28) was taken as base model to adjust it using the mixed effects technique in the h-d ratio of the seven species under study. This equation includes dominant height as an indicator of productivity, basal area as a simple and objective measure of trees occupancy and that in general is better than the number of trees per hectare, since it combines the number of trees with its size (^{Huang et al., 2009}), and the number of trees generates predictions subject to instantaneous changes according to the thinnings or silvicultural interventions to which the masses are subjected (^{Crecente et al., 2010}). These variables have been used in modeling the h-d relationship by different authors for a wide variety of species (^{López et al., 2003}; ^{Calama and Montero, 2004};^{Canga et al., 2007}; ^{Crecente et al., 2010}; ^{Corral et al., 2014}; ^{Sharma and Breidenbach, 2015}).

Mixed effects model

While adjusting Equation (28), it became clear that the parameters with the greatest variability were those that define the asymptote and the shape of the h-d curve (Table 5). Therefore, in a first step it was adjusted with the parameters that define the asymptote (b₀ , b₁ ) and the form (b₃ , b₄ ) as random; however, since no convergence was found, several combinations were tested with one and two random parameters associated with the fixed parameters until convergence was achieved.

The best combination was found by associating the shape parameter (b₃ ) with a random parameter u_i in an additive way (b₃ + u_i ). Similar results have been reported by ^{Sharma and Parton (2007)}, ^{Vargas et al. (2009)} and ^{Corral et al. (2014)}, who considered only one random parameter within the mixed model.

The generalized final model with mixed effects is expressed as follows:

Where:

b₀ - b₄ = Fixed parameters (common to all sample units)

uj~N0,σu = Random parameter (specific for each sample unit) with a normal distribution with mean equal to zero and variance due to random effect (σu2)

hij and εij = Height and error estimated by the model for the i^th tree observation in the j^th sample unit

The mixed effects model (Equation 29) showed an improvement in precision of 6.6 % in RMSE value with respect to the base model (Equation 28); in addition to this, with the inclusion of the random parameter among the sample units, a 2.8 % increase on the value of R² was achieved, explaining the observed variability. On the other hand, the value of the average bias in both equations was similar, but BIC, which evaluates the simplicity of the model, was reduced by 2.9 % (Table 5).

Figure 2a shows residual dispersion against heights predicted by (Equation 29); a homogeneous variance is observed throughout the data interval, and there is no systematic pattern in the variability of residuals. Therefore, the random effect helped eliminate heterogeneity in the variance (^{Fang and Bailey, 2001}). On the other hand, the graph of residuals against first delay (Figure 1b) does not indicate an autocorrelation tendency, which means that it was removed when the random effect among sample units was modeled by expanding parameter b₃.

Figure 2 Residuals against the predicted heights (a), residuals versus the first delay (b) estimated with Equation 29. 

Model Calibration

The results of calibration alternatives to improve predictions are shown in Table 6. The highest values of RMSE were obtained by including a subsample of one to three trees with extreme diameter categories (minimum and maximum). In contrast, with the inclusion of trees belonging to central diametric categories, a slight reduction in RMSE was identified. Although there is no methodology for defining the most appropriate sample size in the calibration of a mixed-effects model, it is noted that the RC6 alternative was found to be more accurate while measuring total height in a subsample of three randomly selected trees with diameters close to the second quartile of the site’s diametric structure, given that it presented a 4.6 % reduction in RMSE in regard to the estimated value during the adjustment of Equation (28) by ONLS. This value was very close to the one generated with the mixed effects model adjustment (measurement of the total height of all trees) of (Equation 29) by NLME procedure.

The mixed effects model registered an increase in accuracy, because it includes the variability that exists in the total height of trees within and between plots (^{Calama and Montero, 2004}). Given that mixed forests have a greater variability in total tree height compared to those that come from stands with one or two species, it is estimated that the predictions of the developed mixed effects model is acceptable, since it explains 82 % of the data variability. In experiences with data from regular masses, a close to 90 % prediction capacity has been achieved (^{Calama and Montero, 2004}; ^{Diéguez et al., 2005}; ^{López et al., 2012}).

The calibration of mixed-effects models is an advantage for forest managers, since it is possible to determine the height of all the trees in the stand from only a few trees. ^{Calama and Montero (2004)} propose four randomly selected trees; ^{Corral et al. (2014)} define that three trees close to 10 % of the average diameter are sufficient, while ^{Castedo et al. (2006)} and ^{Crecente et al. (2010)} propose the three smallest trees of the stand's diametric structure. As the number of trees selected to calibrate a mixed-effect model increases, RMSE is generally reduced and the predictive capacity of the model is improved. However, this is not always justified since the cost of sampling increases, so calibration with a small sample of trees is the main advantage of the mixed effects models when applied in the forest field (^{Castedo et al., 2006}).