Study area

The site is located in the Transverse Neovolcanic Axis between latitudes N 19° 21’ - 19° 34’' and longitudes W 102° 08’ - 102° 17’ (Figure 1), with an average altitude of 2 550 m, a humid temperate climate, average annual temperature of 18 °C and rainfall of 1 600 mm. CINSJP has a total area of 18,318 ha, of which 9 914 ha are timber production, the latter comprising 8 927.5 ha of natural forest and 986.5 ha of commercial plantations. Management is carried out under the Silvicultural Development Method (MDS) with regeneration of seed trees and the Mexican Management Method for Uneven-aged Forests (MMOBI) with regeneration by selection. There are 10 production annuals, nine of which are managed under MDS and one under MMOBI (^{Dirección Técnica Forestal, 2017}).

Figure 1 Geographical location of the study area in the Indigenous Community of Nuevo San Juan Parangaricutiro (CINSJP), Michoacán, Mexico. 

Data

Size-density observations come from 9 559 temporary circular plots of 1 000 m² established in 2016 in natural forest (Table 1), in which we estimated basal area (G), density in number of trees per hectare (N), volume and quadratic mean diameter QMD=40,000πGN ).

The forest inventory identified 26 species that were organized into four groups, based on genus and shade tolerance: Pinus, Quercus, other conifers and broadleaf. Nine Pinus species were located in 98.1 % of the plots: Pinus pseudostrobus Lindl., Pinus montezumae Lamb., Pinus leiophylla Schltdl. & Cham., Pinus devoniana Lindley, Pinus douglasiana Martínez, Pinus teocote Schiede ex Schltdl., Pinus lawsonii Roezl., Pinus ayacahuite Ehren., Pinus patula Schltdl. & Cham.; seven of Quercus (72.8 % de las parcelas): Quercus rugosa Neé, Quercus laurina Bonpl., Quercus candicans Neé, Quercus castanea Neé, Quercus dysophylla Benth., Quercus obtusata Bonpl., Quercus magnoliifolia Neé; two grouped in ‘other conifers’ (17.4 %): Abies religiosa (Kunth) Schltdl. & Cham. and Cupressus lusitanica Klotzsch ex Endl.; and eight broadleaves (34.6 %): Alnus jorullensis Kunth, Alnus acuminata H. B. K., Ternstroemia pringlei (Rose) Standl., Clethra mexicana DC., Arbutus xalapensis Kunth, Carpinus caroliniana Walter, Ilex pringlei Standl. and Tilia mexicana Shltdl.

In a preliminary analysis, a selection of plots with maximum density was carried out. The method consists of calculating the Reineke’s SDI for each site and then defining a percentage to select the plots with the highest occupancy by SDI (i.e. 10 %) (^{Solomon & Zhang, 2002}). In this study, 956 sites were selected using the above technique, of which 91 % are under management with MDS and 9 % with MMOBI. The plots covered a variety of species mixtures combinations. According to the proportion of basal area in each plot, Pinus, broadleaf, Quercus and other conifer groups dominated 84.73 %, 9.21 %, 3.24 % and 2.82 % of the total, respectively. On average there were three to four species in the selected plots, with a range from one to eight species.

Estimation of the upper limit of the maximum size-density relationship

We used two approaches, the first one took into account stand composition with global parameter estimation, and the second one incorporated composition as a proportion per species group. With the first approach, MSDR was estimated with the linearized Reineke’s equation or potential model:

ln⁡N=β0+β1∙ln⁡QMD+ϵ

where, N is the number of trees, QMD is the quadratic mean diameter, ln is the natural logarithm, β₀ = ln (β₀) and ε is the error.

With the same approach, the exponential model in its linearized form proposed by ^{Quiñonez-Barraza and Ramírez-Maldonado (2019}) was also used:

ln⁡N=β0+β1∙QMD+ϵ

The second approach was based on the potential model with species proportion (SP_PM) of ^{Torres-Rojo and Velázquez-Martínez (2000}) with two assumptions:

(i) The intercept changes according to the proportion per species group and shows the same slope (SP_PM1):

(ii) ln⁡N=α1∙PS1+α2∙PS2+α3∙PS3+α4∙PS4+β∙ln⁡QMDw+ϵ

(iii) where, α _i is the intercept weighted by PS _i , i is the i-th group (1 = Pinus, 2 = Quercus, 3 = other conifers and 4 = broadleaves), PS _i is the ratio of basal area per group to the total basal area of the plot, and QMD _w is the quadratic mean diameter weighted by the PS _i ratio of each species QMDw=QMD∙PS1+QMD∙PS2+QMD∙PS3+QMD∙PS4=DCM.

Intercept and slope change according to the proportion or grouping (SP_PM2); the estimators are calculated with the following equation:

ln⁡N=α1∙PS1+α2∙PS2+α3∙PS3+α4∙PS4+(β1∙PS1+β2∙PS2+β3∙PS3+β4∙PS4)∙ln⁡QMDw+ϵ

Fitting techniques

Base models and variants of the SP_PM were fitted with two techniques. One of these is stochastic frontier regression (SRF) (^{Aigner et al., 1977}; ^{Meeusen & van den Broeck, 1977}) which incorporates a compound error (ε) with two random terms v∼N(0,σv2) and u for which the semi-normal SFR-SN ui∼N+(0,σu2) and normal-truncated SFR-NT ui∼N+0,σu2 (^{Meeusen & van den Broeck, 1977}). Values of σ2=σv2+σu2 and γ=σv2σ2 were obtained with the parameterization of ^{Battese and Corra (1977}). With a likelihood ratio test with X ² and under H ₀ : µ = 0 we evaluated whether the SN distribution was more appropriate than NT for µ (^{Bi, 2004}).

The other technique is quantile regression (QR), which uses quantiles τ∈(0,1) to describe non-central positions of a distribution (^{Koenker & Bassett Jr., 1978}). Parameter estimation solves an optimization problem by minimizing an asymmetric function with eq. β^τ=minβ∈R⁡∑i=1nρτ∙N-N^, where ρτ∙N-N^=N-N^∙τ-1∙N-N^<0 is a checking function, and thus:

ρτ∙N-N^=N-N^∙(τ-1),  si N-N^<0N-N^∙(τ),  si N-N^≥0;  τ=∈(0.95, 0.975, 0.99)

The slope estimated with the three quantiles used in the CR was compared with a second likelihood ratio test.

Fittings were carried out with the R® statistical package (^{R Core Team, 2020}), evaluating the significance of the parameters with α = 0.05. For SFR the ‘frontier’ package was used with the maximum likelihood technique, and for QR ‘quantreg’ with the ‘simplex’ algorithm. Standard errors and confidence intervals for QR were calculated using the Hall-Sheather method (^{Koenker, 2019}). The adjustment with SFR was evaluated with the Akaike information criterion (AIC) and the log-likelihood (LogLik) (^{Quiñonez-Barraza et al., 2018}), and QR was evaluated with the Pseudo R², the AIC and the LogLik (^{Condés et al., 2017}). In addition, the dispersion of N vs. QMD was visually evaluated relative to the respective fit line.

The best model of the first approach was used to construct two density guides that related density (N) and average size (QMD) in logarithmic scale. From the maximum density SDI value (SDI_max) of the potential (SDI_R) and exponential (SDI_Q) models, three isolines were plotted to identify three stages of stand development: self-thinning or imminent mortality, constant growth and free growth. The isolines were projected to 70 %, 40 % and 20 % of the SDI_max with the following equations for the potential and exponential models, respectively (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}):

SDIR=NQMDRQMDβ^1

SDIQ=Ne-β^1QMD-QMDR

where, SDI is the stand density index, QMD_R is the reference quadratic mean diameter (i. e. QMD_R = 25 cm), e is the exponential and β^1 is the slope parameter of the potential or exponential model as appropriate.

Upper limit of the maximum size-density relationship excluding stand composition

Table 2 indicates that the estimated parameters for the potential and exponential models with the two regression techniques were significant (P < 0.0001), although the slope of the potential model was not different from -1.605 according to the t-student test (α = 0.05). The intercept and slope, estimated with SFR, were similar for the semi-normal (SN) and normal-truncated (NT) forms. Values of / σv2 were 0.015 potential) and 0.018 (exponential), with σu2 close to zero in both models; SFR-SN was better than SFR-NT for µ according to the likelihood ratio test (P > 0.05). Estimates of µ with SFR-NT were not significant (P > 0.05) and close to zero, so their fit is not appropriate (^{Weiskittel et al., 2009}).

In both the potential and exponential models, the intercepts and slopes were greater in QR than in SFR; as τ increased in the QR fit, the parameters also increased in both models (Table 2). The likelihood ratio test detected no significant differences in slope consistent with the quantiles used (P > 0.05).

According to Table 3, the potential model was the best for analyzing MSDR with higher Loglik and Pseudo R² values in the adjustments with SFR and QR, respectively. SFR-SN was better than SFR-NT based on Loglik and AIC, while, in QR, the statistics with τ = 0.95 were better than the other two quantiles in both the potential and exponential models.

MSDR trend was similar in SFR-SN and SFR-NT; however, QR fit lines were above the previous ones and, as the value of τ was higher, these approached the upper limit of the data (Figure 2). QR fit with τ = 0.99 was adequate for both the potential and exponential models, by its closer approach to the MSDR limit of the selected sites with maximum density (dark dots) and of the total plots (gray dots); nine and four plots exceeded the MSDR limit estimated with the potential and exponential models, respectively.

Figure 2 Trajectory of the maximum size-density relationship obtained with the potential and exponential models adjusted with quantile regression (QR) and regression with stochastic semi-normal (SFR-SN) and normal-truncated (SFR-NT) frontier. QMD: quadratic mean diameter. 

Upper limit of the maximum size-density relationship according to species ratio

The estimators for SP_PM1 and SP_PM2 were different from zero (P < 0.0001) except for β^2 for SP_PM2 using QR with τ = 0.99 (Table 4). The estimators found with the two variants of SFR were similar, although different from those of QR. The slopes β^'s ) of SP_PM1 were steeper in QR compared to SFR. Initial occupancy α^i in terms of N, when QMD approaches zero, was lower in the Pinus and Quercus groups than in broadleaves and other conifers, for SFR and QR with τ = 0.95 and τ = 0.975; however, in QR with τ = 0.99 the opposite was observed. SP_PM2 intercept estimators were similar to those of SP_PM1, although with variable mortality rates per species group. SP_PM2 slopes were steeper for broadleaves and other conifers relative to Pinus and Quercus.

SP_PM2 estimation was better than SP_PM1 with lower AIC with QR, higher Loglik at QR and SFR, and higher Pseudo R2 at all three τ of QR (Table 5). Also, SFR-SN was better than SFR-NT, and QR with τ = 0.95 showed the best goodness-of-fit statistics compared to the other two quantiles for SP_PM1 and SP_PM2. Consistent with the above results, QR with τ = 0.95 of SP_PM2 is higher for defining MSDR with species proportion.

Upper limit of the maximum size-density relationship excluding stand composition

In this analysis, QR adequately estimated MSDR for the potential and exponential models. The benefits of QR in relation to SFR center on the assumptions of the error distribution, as it allows inferences to be made correctly with other quantiles (^{Salas-Eljatib & Weiskittel, 2018}). However, QR can be subjective because a τ needs to be previously selected for adjustment (^{Condés et al., 2017}; ^{Tian et al., 2021}). Despite the above, in this study we sought an upper limit that covered most of the plots with maximum density, so QR is considered appropriate to define MSDR. QR has been used with the potential model in other studies and satisfactory results have been found estimating MSDR, when τ values close to one are used (^{Aguirre et al., 2018}; ^{Andrews et al., 2018}; ^{Condés et al., 2017}).

The use of the potential and exponential models in mixed-species forests has been analyzed and compared in few studies (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}). Due to differences in the mathematical structures of both models, the value of the slope varies significantly, and it is useless to make a statistical comparison. Nevertheless, at the biological level, the models show distinct forms defining the allometric dynamics between average tree size and growing space (^{Condés et al., 2017}; ^{Pretzsch & Biber, 2005}), as well as the representation of interspecific competition (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}; ^{Weiskittel et al., 2009}).

The results of the potential model showed constant mortality rates during stand development with mixed species. The MSDR obtained with the potential model is essential in identifying the maximum average size that stands reach before self-thinning (^{Zeide, 1985}), regardless of the initial density (^{VanderSchaaf, 2010}). Therefore, the maximum density line of the potential model is only representative for a specific interval of stand dynamics (VanderSchaaf, 2010; ^{Zeide, 2010}), mainly when complete crown closure occurs (^{Quiñonez-Barraza et al., 2018}; ^{Zeide, 2005}). In the present study, the MSDR estimated with the potential model covered an interval from 25 cm to 50 cm of average tree size (Figure 2), showing that mixed stands with high densities reach the maximum density line faster (^{Ningre et al., 2016}).

The average MSDR estimated with the exponential model showed a different pattern from that of the potential model, and its concave trajectory covered most of the plots under study. The value of the intercept estimated with the exponential model showed a more realistic value than the potential model and better described the initial stages of the stands for the study area. The maximum density line of the exponential model can be considered as an average representation that characterizes the various competition events (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}), according to the combination of size-density observations evaluated.

The curvilinear relationship of the exponential model showed slopes close to zero when the stand has high densities and a small average size, while the slope becomes steeper as tree size increases (^{Zeide, 2010}). This behavior suggests that older stands or stands with mature trees are not as dense as young stands due to the increasing proportion of gaps (^{Zeide, 1995}). The same trend was found when the Reineke model was modified to describe MSDR in curvilinear form for stands of Pinus elliotti var. elliotti Engelm. (^{Cao et al., 2000}), Pinus tadea (L.) (^{VanderSchaaf & Burkhart, 2008}) and Picea abies (L.) Karst. and Fagus sylvatica L. (^{Schütz & Zingg, 2010}). However, few studies have yet applied these new formulations to mixed-species stands.

With the slope estimated by both models, SDI_max values were calculated with a QMD_R equal to 25 cm (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}; ^{Reineke, 1933}). For the potential model, SDI_max was 883 and for the exponential model it was 800, values that together with of the respective model were used to delimit the following growth zones: imminent mortality (SDI = 618 and 560), constant growth (SDI = 309 and 280) and free growth (SDI = 177 and 160) for potential and exponential models, respectively (Figure 3). Differences between density guides lie in the linear and curvilinear projections of curve families of the respective SDI. For the potential model, N projections decrease linearly as the average tree size increases, meanwhile such a trend is curvilinear for the exponential model (^{Cao et al., 2000}; ^{Schütz & Zingg, 2010}).

Figure 3 Density guidelines based on potential model (a) and exponential model (b) for mixed forests of Nuevo San Juan Parangaricutiro, Michoacán, Mexico. QMD: quadratic mean diameter. 

Upper limit of the maximum size-density relationship according to species ratio

MSDR estimation in mixed stands is complex (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}; ^{Rivoire & Moguedec, 2012}), because all possible combinations of species found in a specific location have to be considered (^{Torres-Rojo & Velázquez-Martínez, 2000}). Plots meeting this condition were selected in this study; however, the lack of observations showing QMD sizes approximating 0 cm up to 15 cm may be a factor for not having a better representation of MSDR (^{VanderSchaaf & Burkhart, 2012}).

The parameters estimated with the best model using species ratio described the differences between mortality rates, tolerance, and intra- and interspecific competition (^{Aguirre et al., 2018}; ^{de Prado et al., 2020}; ^{Torres-Rojo & Velázquez-Martínez, 2000}), as well as the allometric relationship to reach the maximum potential of a combination of species (^{Weiskittel et al., 2009}).

In general, other conifer and broadleaf groups tend to have high densities at the early stand stage, but show a high mortality rate as they develop, because they demand more growing space at the mature stage (^{Pretzsch & Biber, 2005}), as well as robust canopy development (^{Zeide, 1985}). Pinus and Quercus groups showed the lowest intercept values, indicating that in the mature stage of the stand they will dominate with a lower demand for growing space (Pretzsch & Biber, 2005). These groups show high growth and survival capacity in the presence of interspecific competition (^{de Prado et al., 2020}).

β _i slopes varied considerably among species groups. The other conifer and broadleaf groups were more shade tolerant, because they showed more negative values compared to Pinus and Quercus. This suggests that each species has its own allometric size-density relationship, resource demand and growing space (^{Pretzsch & Biber, 2016}). ^{Weiskittel et al. (2009}) analyzed MSDR in plantations and natural forests of three species and found that Tsuga heterophylla (Raef.) Sarg. had a more negative slope with higher shade tolerance than Alnus rubra (Bong.) Carr. and Pseudotsuga mensiezii var. menziesii (Mirb.) Franco. Other studies also concluded that species with steep β _i slopes, in mixed-species stands, are more shade tolerant (^{Andrews et al., 2018}; ^{de Prado et al., 2020}; ^{Torres-Rojo & Velázquez-Martínez, 2000}).

Shade tolerant species may not have guaranteed survival in mature stages of the stand, as this may be influenced by the self-tolerance of the species, mainly by growth capacity and survival of competing trees in the stand (^{Zeide, 1985}). The results indicate that species less tolerant to shade may be more self-tolerant than tolerant species, a characteristic that may represent a benefit for the permanence in the stand with lower mortality rates (^{Zeide, 2005}). However, it is important to analyze other characteristics such as crown, site productivity and wood specific gravity to understand the specific differences in self-thinning ratio for mixed stands, as well as its allometric parameter (^{Woodall et al., 2005}).