Study area

The study area was the La Victoria ejido, located in the southwest of the state of Durango, Mexico, between 105°25'39.465" W and 23°43'53.022" N. The altitude above sea level fluctuates between 2 400 and 2 850 m. The predominant climate is semi-cold temperate, with summer rains. The average annual rainfall varies from 900 to 1 200 mm, and the average annual temperature varies from 5 to 18 °C (^{García, 2004}). The forests are made up of species of the Pinus, Quercus, Juniperus, Cupressus, Pseudotsuga, Arbutus and Alnus genera, which form mixed and irregular structures.

Data

Data were collected in 2017, during the timber forest management inventory, from 2 252 temporary sampling sites with a circular shape and a surface area of 0.10 ha each, located in the La Victoria ejido, in the municipality of Pueblo Nuevo, Durango, using a stratified systematic sampling design. A surface area of 10 876 ha divided into 518 management units (stands), was covered. The database included records of eight species of conifers: Pinus cooperi C. E. Blanco, P. durangensis Martínez, P. leiophylla Schiede ex Schltdl. & Cham., P. teocote Schied. ex Schltdl. & Cham., P. engelmannii Carrière, P. strobiformis Engelm., P. herrerae Martínez y Juniperus deppeana Steud.; and eight species of broadleaves: Alnus acuminata Kunth, Arbutus bicolor S. González, M. González & P. D. Sørensen, Quercus durifolia Seemen ex Loes., Q. sideroxyla Bonpl., Q. obtusata Bonpl., Q. coccolobifolia Trel., Q. viminea Trel. y Q. candicans Née.

For each tree with a normal diameter equal to or larger than 7.5 cm, the normal diameter was measured at 1.3 m (d, cm) with model 122450 Ben Meadows^® diameter tapes; the total height of the tree (h, m) and of the stem without branches (h, m) was measured with the model Pm5/360pc SUUNTO^® clinometer. These data were used to estimate the following variables: number of trees per hectare (N, trees ha^-1), basimetric area (G, m² ha^-1), root mean square diameter (d _g , cm), dominant height and diameter, estimated as the average of the 100 largest trees per hectare (H ₀ , m and D₀, cm, respectively) (^{Assmann 1970}). Table 1 shows the main descriptive statistics of the database utilized in the adjustment of the PDFs.

Models

Evaluation of PDFs

The analysis of the PDFs' ability to fit was based on numerical and graphical comparisons of the errors (residuals). The statistics used for evaluating goodness-of-fit were as follows: (i) average bias (AB), whose optimal value is zero, implying that the deviations between the actual and predicted values are compensated for the total of the analyzed data, and (ii) Root Mean Square Error (RMSE), which is an indicator of accuracy: the lower the value, the smaller the difference between actual and predicted values (^{Gorgoso-Varela et al., 2020}). The corresponding mathematical expressions are as follows:

Where:

Fxi and F^xi = Number of observed (actual value) and estimated trees, respectively, in diameter class i

n = Number of observations

p = Number of parameters of the PDF

AB and RMSE values were estimated for each PDF as the average of the relative frequencies of each diameter class (tabulated in 2 cm classes) for each sampling site and stand. The Kolmogorov-Smirnov (KS) test was used to validate whether the empirical and theoretical cumulative frequencies are similar (^{Sokal and Rohlf, 2012}). The most notable maximum difference in absolute value between the two frequencies is given by the value of D _n in the following expression:

Where:

F(x _i ) and F ^* (x _i ) = Cumulative empirical and theoretical distribution, respectively

The value of D _n is compared to the one obtained from a table according to the number of data (n) and the chosen significance level (α=20 %).

Modeling of PDF parameters

The parameters of the future distribution were estimated using the moments method of ^{Shifley and Lentz (1985)}, so that the parameters of the PDF are directly related to the moments of first and second order with respect to the origin; i.e. the average diameter d¯, and the variance (s2). In turn, the latter is related to the root mean square diameter of the stand (d _g ). The d _g value of the stand at any given age (years) can be estimated based on the values of density (N, trees ha^-1) and basal area (G, m² ha^-1) [dg=4∙Gπ∙N], obtained from the models to be built for the species present in the stand. However, since the value of d¯ cannot be obtained directly from the stand variables, a relationship between d _g and stand variables was adjusted.

Given that the d _g is always equal to or greater than d¯, a mathematical model was used to reduce the bias generated by the use of the originally estimated d _g [d¯=dg-exp⁡(β∙X)], where X is a matrix of stand variables that characterize the state of the stand at any age and can be calculated using a static or dynamic stand model; β is the vector of parameters to be estimated. Once the value of the stand’s d¯ and dg has been determined, the variance (s2) (second order moment of the distribution with respect to the mean) is calculated by means of the mathematical relation: s2=dg2- d¯2 (^{Frazier, 1981}).

In the case of the S_B Johnson and Beta PDFs, the parameter recovery using the moments method required fitting another equation that related the maximum diameter of the distribution to the stand variables estimated from the growth models [dmax=⁡(β∙X)]. With these equations, it will be possible to use the relationships to estimate λ, γ, and δ, assuming that the location parameter ε corresponds to the minimum inventoriable diameter, which in this case is zero (^{Gorgoso et al., 2012}).

Adjustment and evaluation of the PDFs

Table 2 presents the statistics used to evaluate the quality of the fit, as well as the number and percentage of sampling sites that passed the KS test for the evaluated PDFs.

The statistics reveal a significant difference between the ability of the Weibull PDF to fit the Beta and S_B Johnson PDFs. This was demonstrated by the lower AB and RMSE values in all diameter classes and the higher percentage of sample sites passing the KS test at a significance level of 20 %. This gives a clear idea of the validity of using this distribution to adequately describe the wide range of shapes shown by the diameter distributions for this type of forest stands. ^{Pece et al. (2000)} suggest that the 20 % significance level of the KS test is the most demanding, as it reduces the minimum deviations for non-rejection of concordance.

The results reflect the flexibility and parsimony of the Weibull PDF because it only requires the estimation of two parameters, compared to the four that have to be estimated for the Beta and S_B Johnson PDFs. Several authors also cite the flexibility of the Weibull PDF to describe a wide range of unimodal distributions (j-inverted, exponential and normal), furthermore, they highlight its ability to be expressed in its closed form, as a cumulative distribution function (mathematical simplicity) (^{Bailey and Dell, 1973}; ^{Reynolds et al., 1998}; ^{Maldonado and Návar, 2002}; ^{Liu et al., 2004}; ^{Palahí et al., 2006}; ^{Lei, 2008}; ^{Pogoda et al., 2019}).

When analyzing the Weibull PDF adjustment methodologies, a marginal gain was observed in the number of plots that passed the KS test by using the MM (70.8 %), in contrast to the ML method (67.0 %). In this regard, papers such as the one by ^{Corral-Rivas et al. (2015)} confirm the suitability of the Weibull PFD adjustment according to the MM, compared to four alternative methodologies evaluated for the mixed and patchy forests of northwestern Durango. For their part, ^{Sun et al. (2019)} recommend estimating the parameters of this distribution with the method of moments and nonlinear regression for mixed and irregular masses in China.

For the detection of some kind of systematic pattern, the behavior of the PDFs was analyzed based on the evolution of the AB and RMSE in the prediction of the number of trees per diameter class (Figure 1).

Figure 1 Evolution of the average bias and Root mean squared error (RMSE) in the prediction of trees by diameter class of the PDFs evaluated at the adjustment stage. 

Figure 1 shows that the S_B Johnson PDF exhibited the highest AB and RMSE values in the smaller diameter classes, indicating that it underestimates the frequencies in the smaller diameter classes and overestimates in the diameter classes above 45 cm. The Beta PDF also underestimates the lower frequencies and overestimates the intermediate frequencies (30-40 cm), improving its accuracy from the 45 cm diameter class upwards. Finally, the estimation of the frequencies with better precision over the whole interval of the observed diameter classes to have a very similar trend was achieved by adjusting the Weibull PFD according to the MM and the ML method.

When comparing the MM and ML adjustment methods of the Weibull PDF, we observed that there are no significant differences between the values of the goodness-of-fit statistics or their behavior in the prediction of the frequencies in the whole interval of the diameter classes. In their studies in northwestern Durango, Mexico and northwestern Spain, respectively, ^{Gorgoso et al. (2012)} and ^{Quiñonez et al. (2015)} obtained the best results with the MM and ML adjustment methods. Therefore, in order to select the Weibull PDF fitting methodology, the number of excluded plots was analyzed with the KS test, it was determined that the MM was slightly larger (70.8 %), therefore, the rest of the analyses are limited to this function with the parameters estimated with the MM for the modeling of the stand diameter distributions at any time interval.

Modeling of the diameter distribution

The estimation of the parameters (scale and shape) of the Weibull PDF using the parameter recovery methodology with the MM, exhibited simplicity and yielded desirable results according to the values of the evaluated goodness-of-fit statistics. Of the total number of sampling sites, 62 % passed the KS test at a significance level of 20 % (1 394 sampling sites). This result validates the Weibull PDF as suitable for adequately describing the wide range of shapes observed in the diametric distributions of forests under forest management at any stage of their development. Similar results were obtained by ^{Siipilehto and Mehtätalo (2013)}, who compare methodologies for parameter recovery and prediction of the two-parameter Weibull PDF in Pinus sylvestris L. stands in Finland; they demonstrated that parameter recovery provides better compatibility with stand characteristics, especially as it estimates the volume per diameter class with greater precision. Furthermore, ^{Quiñonez et al. (2015)} show that in the case of incoetaneous masses with a mixture of species, the use of the MM for the recovery of the parameters of the Weibull PDF yields more accurate results.

Since the mean diameter of a stand is always less than or equal to its root mean square diameter, the adjustment relationship must meet this compatibility condition. For this purpose, it was determined that the mathematical model that meets this constraint has the following structure:

Where:

d _g = Root mean square diameter (cm)

N = Number of trees per hectare

G = Basimetric area (m² ha^-1)

RS = Relative spacing index (%)

β _i = Parameters to be estimated

The fit of Equation 19 was good, given that the values of the coefficient of determination (R ² ) and of the RMSE were 0.933 and 0.722 cm, respectively. Table 3 shows the estimated values of the parameters with Equation 19.

Figure 2 shows the validation process derived from the representation of the evolution of the average bias and the Root Mean Square Error by diameter classes obtained in the adjustment and modeling stage of the Weibull PDF.

Figure 2 Evolution of the average bias and the root mean square error by diameter class with the Weibull PDF in the fitting and modelling phase. 

Figure 2 compares the error in terms of the number of predicted trees per diameter class, showing that the largest deviations overestimate and are located in the diameter classes below 15 cm. The reason may be the result of the influence of silvicultural interventions and the effect of competition on young stands with different species mixtures. This same behavior was observed by ^{Corral-Rivas et al. (2015)} for conifer and broadleaf stands in northeastern Durango, Mexico. Therefore, since these are stands with the highest number of trees per unit area, any change in their structure increases their variation. In this regard, ^{Sandoval et al. (2012)} conclude that the accuracy in the prediction of diameter structures with the Weibull PDF increases with the age of the plantation of three dendroenergetic species in Chile.

The results of the present study are very similar to those documented by ^{Sun et al. (2019)} for mixed and patchy stands in China, where the intermediate clear-cutting of Pinus tabulaeformis Carrière resulted in stands with a bimodal distribution.