Use of the Weibull function to model maximum probability of abundance of tree species in northwest Mexico

Abstract

Key message

Plant abundance is controlled by environmental factors and, theoretically, a suitable function can be used to predict site preferences of particular species. This study represents a first attempt to use the Weibull function to define the bioclimatic niche of some trees, which may help in understanding particular adaptations of individual species to the environment.

Context

The plant–environment relationships studied in this work, and the definition of points of maximum probability of abundance, constitute an important basis for understanding species-specific adaptation to certain environmental conditions. The method presented here can be applied to any location.

Aims

The aim of this study was to test the efficacy of the Weibull function (with two and three parameters) to model the relationship between environmental variables and abundance of tree species in the temperate forests of Durango (Mexico).

Methods

Six environmental variables and 16 tree species were selected for study in 1,804 sample plots located in the Sierra Madre Occidental, Durango, Mexico.

Results

Good models fits were obtained for the abundance of 11 out of the 16 species in relation to one or more of the six predictors; this enabled points of maximum probability of abundance to be defined for each species. Significant fits were obtained for 31 models: 15 with the 3-parameter Weibull function, six with the 2-parameter approach, and 10 with both equations.

Conclusion

The Weibull function was generally flexible enough to model the relationship between environmental variables and tree species abundance.

1 Introduction

The abundance of forest species is determined by the action of multiple factors, mainly abiotic and biotic factors, geographical accessibility and the capacity of organisms to adapt to site conditions. These factors exert effect with different intensities and at different levels in all living organisms in natural systems (Soberón and Peterson 2005; Saupe et al. 2012).

Determination of the bioclimatic preferences of forest tree species is important for understanding the ecology and genetic evolution of the species and for predicting their capacity to adapt to new environments; it is also very useful in the study of global climate change and how this affects vegetation (Hegel et al. 2010; Lenoir et al. 2010). Models that help visualize different scenarios and biogeographic boundaries are usually used to define bioclimatic preferences (Sáenz-Romero et al. 2010).

Several efforts have been made to define the climatic preferences of forest tree species, as well as to study the relationship between the distribution of tree species and environmental variables (Mellert et al. 2011) and to determine the response to climatic gradients of oak species (Uğurlu and Oldeland 2012).

Knowledge about the bioclimatic niche of different species is important for the above-mentioned purposes and also for good planning and decision making in the integrated and sustainable management of forest resources. Nevertheless, accurate determination of the bioclimatic niche of a species is not an easy task, because of the complex interactions between the components of natural systems. However, if the specific values of the main climatic elements in a site where a species grows and develops are calculated, this would be an important step in defining its bioclimatic niche and thus its potential distribution.

It is possible to use the Weibull probability density function (PDF) to characterize the behaviour of a population in relation to certain phenomena. The advantage of using this mathematical tool lies in its flexibility for representing various forms in the curve, based on the values of the parameters (mainly the shape parameter) (Murthy et al. 2004). The Weibull function enables density curves to be constructed for a variable of interest in relation to one or more independent variable (Moritz et al. 2009). Several authors have used the Weibull distribution function in forestry, particularly for modelling the diameter distribution of forest stands (Mehtatalo 2004; Nord-Larsen and Cao 2006; Castedo-Dorado et al. 2007; Vallejos-Barra et al. 2009). However, we did not find any studies in which this model had been used to represent the abundance of tree species in relation to climatic variables.

The aims of this study were as follows: (1) to test whether two- and three-parameter Weibull functions can accurately represent the abundance of 16 tree species of Durango’s temperate forest, in relation to six climate variables, and (2) to define the climatic preferences of selected tree species.

The results of this study will be useful for identifying the range of distribution and the potential area of maximum probability of abundance of species of interest on the basis of environmental factors. This information could be used as a reliable basis for constructing ecograms and eco-maps representing bioclimatic niches, which is important for forestry management planning and for developing other types of studies.

2 Materials and methods

2.1 Study area

The study area is located in the Sierra Madre Occidental mountain range, within the State of Durango, Mexico, between geographical coordinates 26° 50′ and 22° 17′ north and 107° 09 and 102° 30′ west and covering an area of approximately 6.33 million ha (Fig. 1). The main vegetation in the region consists of forests of Pinus, Pinus-Quercus and Quercus-Pinus, small portions of temperate mesophytic forest, and some tropical deciduous and semi deciduous forests in the western flanks (González-Elizondo et al. 2012).

Fig. 1

Location of the study area

The sampling units used in this research were established by the National Forestry Commission (CONAFOR) for the National Inventory of Forests and Soils (INFyS 2004–2009). The sampling design was systematic and stratified by clusters. The Land Use and Vegetation Series III chart (scale, 1:250,000), produced by INEGI (National Institute of Statistics, Geography and Informatics), was used to project the sampling points (CONAFOR 2009). Each main sample was referred to as a “conglomerate” and comprised four sub units (called “sites”) each of 400 m². Thus, each conglomerate comprised an area of 1,600 m². A total of 1,804 conglomerates were distributed at 5 km from each other in the study area (CONAFOR 2009).

Sixteen of the most abundant and representative species of the Sierra Madre Occidental (González-Elizondo et al. 2012) were chosen for this analysis, since they are the most important both economically and ecologically. The selection included eight species of the genus Pinus and five of the genus Quercus: P. arizonica Engelm, P. strobiformis Engelm, P. cembroides Zucc, P. cooperi D. Don, P. durangensis Martínez, P. engelmannii Carr, P. leiophylla Schl. & Cham, P. teocote Schl. & Cham, Q. arizonica Sarg, Q. crassifolia Humb & Bonpl, Q. grisea Liebm, Q. magnoliifolia Née and Q. sideroxyla Humb. & Bonpl. The other three cases considered were Alnus acuminata Kunth, Juniperus deppeana Steud and a group of Arbutus species, including Arbutus bicolor González, Arbutus arizonica Gray, Arbutus madrensis González, Arbutus tessellata Sørensen, Arbutus occidentalis McVaugh & Rosatti, Arbutus mollis HBK and Arbutus xalapensis HBK.

In the Forest National Inventory, all trees of the genus Arbutus were recorded as Arbutus xalapensis Kunth. However, it has been specified recently that, in the study area, these trees correspond to the species mentioned above, which have very similar characteristics and are often confused with each other (González-Elizondo et al. 2012). Therefore, in this work we considered these trees together and refer to them only as Arbutus spp.

In this study, density was used as an indicator of abundance (number of trees per species per conglomerate, with diameter greater than or equal to 7.5 cm at a height of 1.3 m above ground). For tree species that were not recorded within the conglomerate, the density value was equal to zero. The density term used here is that defined by Jonsson et al. (1992).

For each species, the mean values and ranges of the following six predictor variables analyzed are shown below in Table 2: mean cumulative precipitation between April and September, which is referred to as the growing season precipitation (GSP, mm); mean temperature in the coldest month (MTCM, °C); mean temperature in the warmest month (MTWM, °C); degree days above 5 °C (DD5, °C); Julian date of last spring frost (SDAY, day) and elevation above sea level (EASL, m). According to a previous study (Martínez-Antúnez et al. 2013), these six variables correlated closely with species density; these variables have also been used and recommended by Nadezda et al. (2006), Rehfeldt et al. (2006) and Sáenz-Romero et al. (2010) for studies related to the presence and distribution of forest species. The correlation matrix for these climatic variables and their principal component analysis (PCA) were reported by Silva-Flores et al. (2014).

The data on climate variables were obtained for each conglomerate by accessing the server of the Forest Service Department of Agriculture of the United States (http://forest.moscowfsl.wsu.edu/climate, which models climatic data between 1961 and 1990, and considering data from more than 6,000 weather stations in Mexico, South America, Guatemala, Belize and Cuba, of which 183 are in the state of Durango (Crookston et al. 2008; Sáenz-Romero et al. 2010). The EASL was determined in the field by using a global positioning system (GPS).

2.2 Statistical analysis

The density of each species was modelled by using six predictor variables in the two- and three-parameter Weibull probability density functions (PDFs) (W2p and W3p) in order to evaluate the relative performance of both equations. This density curve was used to calculate the maximum probability of abundance of each species (the highest point of the density curve).

The Weibull function was chosen because it generates robust and flexible models, which allow analytical expression of the value of the integral by means of the cumulative distribution functions (Murthy et al. 2004; Torres-Rojo 2005).

The density curves generated were characterized by their parameters of shape and scale (W2p) and also by the location parameter (W3p). For easy interpretation, the results were plotted graphically in terms of relative frequency.

The equation for the PDF and cumulative PDF of W2p is

$$ f\left(x\Big|c,b\right)=\frac{c}{b^c}{x}^{c-1}{e}^{-{\left(x/b\right)}^c} $$

(1)

whose cumulative function is

$$ F\left(x\Big|c,b\right)=1-{e}^{-{\left(x/b\right)}^c} $$

(2)

where c > 0 is the shape parameter responsible for the skew of the distribution, and b > 0 is the scale parameter.

The equation of PDF for W3p is

$$ f\left(x\Big|c,b,a\right)=\frac{c}{b}{\left(\frac{x-a}{b}\right)}^{c-1}{e}^{\left(-{\left(\frac{x-a}{b}\right)}^c\right)} $$

(3)

and its cumulative function is

$$ F\left(x\Big|c,b,a\right)=1-{e}^{\left(-{\left(\frac{x-a}{b}\right)}^c\right)} $$

(4)

where c >0 and b >0 mean the same as in W2p, while a >0 is the location parameter.

The parameters, the probability density curves, and the Kolmogorov-Smirnov test (K-S) for the models W3p and W2p were estimated using R software, version 2.13.1 (Marsaglia et al. 2003). The parameters were estimated by the method of maximum likelihood (MLE).

Although the initial values were estimated using the method proposed by Broyden, Fletcher, Goldfarb and Shanno (BFGS) and described by Nocedal and Wright (2006), the MLE method was finally chosen because it provides consistent and asymptotically efficient and unbiased estimators (Zarnoch and Dell 1985; Seguro and Lambert 2000). This yielded the most probable values of the parameters of the distribution for a data set by maximizing the value of what is known as the “likelihood function”, which is based on the PDF for a given distribution.

The K-S test is defined by

$$ {D}_n= \sup x\left|ECD(X)-OC{D}_n(x)\right| $$

(5)

where D _n is the absolute value of the maximum difference between the estimated and the observed cumulative probability distributions.

ECD(X):

= estimated cumulative probability distribution of the variable x.

OCD _n (X):

= observed cumulative probability distribution of the variable x.

sup x :

= supremum or maximum D _n value at a given x value.

The null hypothesis (H _o) states that the differences between the observed and estimated distributions are not statistically significant, considering a significance level of α = 0.20.

The null hypothesis is rejected when the calculated value of D _n is greater than the tabulated value of D _α. For large samples, D _α can also be calculated using the following formula (Marsaglia et al. 2003):

$$ {D}_{\alpha }=\sqrt{\frac{- \ln \kern0.1em \left(\frac{1}{2}\alpha \right)}{2n}} $$

(6)

where D _α is an estimation of the tabulated value of D, with a level of significance of α = 0.20 and n = number of tree species under study.

3 Results

Visually, the density of most tree species studied was represented adequately by W2p and W3p. Comparison between Figs. 2 and 3 suggest that there may be a difference between the models. In Fig. 2, the difference between the two curves cannot be distinguished at first glance. Alternatively, Fig. 4 illustrates a perceptible difference between the curves of the two models. This type of result was common for several species and variables.

Fig. 2

Distribution of observed and projected probability density of Pinus arizonica in relation to elevation above sea level (EASL) for the W2p and W3p models

Fig. 3

Distribution of observed and projected probability density of Quercus grisea in relation to mean temperature in the coldest month (MTCM) for the W2p and W3p models

Fig. 4

Distribution of observed and projected probability density of Pinus engelmannii in relation to Julian date of last spring frost (SDAY) for W2p and W3p models

Figures 4 and 5 show that both models represent the real data quite well, except for the curve generated by W2p in Fig. 5; however, considering the K-S test, the curve generated by W3p in Fig. 5 is not significant either, despite the apparently good graphical output.

Fig. 5

Distribution of observed and projected probability density of Pinus cembroides in relation to mean temperature in the warmest month (MTWM) for W2p and W3p models

When using the EASL and SDAY as predictors, both equations provide significant fits (Dn < Dα) for the observed density of 6 out of the 16 species studied (Table 1). However, in some cases the data fitted significantly to only one of the equations.

Table 1 Parameters estimated and absolute maximum difference (D _n) of the Kolmogorov-Smirnov test for Weibull functions W2p and W3p. EASL Elevation above sea level, MTCM mean temperature in the coldest month, MTWM mean temperature in the warmest month, DD5 degree-days above 5 °C, GSP precipitation in the growing season, SDAY Julian date of the last frost spring

In total, 31 models provided significant fits, of which 15 were with the W3p function, 10 with both functions and 6 with W2p (Table 1).

In summary, for the following 11 out of the 16 selected species a significant fit was obtained for the density by using at least one of six predictors: Alnus acuminata, Arbutus spp., Juniperus deppeana, Pinus arizonica, P. strobiformis, P. cembroides, P. cooperi, P. durangensis, Quercus arizonica, Q. grisea and Q. magnolifolia.

The condition in which the maximum probability of abundance of a tree species occurs was represented by the maximum value of the Weibull function, i.e. the highest point of the probability density curve (see Figs. 2–5), thus yielding the values given in Table 2.

Table 2 Climatic values for the maximum probability of abundance (MaxProb) and the respective extreme values (MinXtrm and MaxXtrm) where each species was found

The maximum probability of abundance values for each tree species and variable may differ relative to those for other species, except in very rare cases, even though all of the species grow in the same region. In some species, the change is only in a few units, as in the case of Juniperus deppeana, Pinus durangensis and Arbutus spp. where the difference considering EASL is less than or equal to 52 m, between each maximum point (Table 2).

In some tree species, the difference between the maximum values of the predictor variables are evident, e.g. in Quercus magnolifolia and Pinus arizonica, which show a difference in EASL of 1,237 m, or Alnus acuminata and Q. magnolifolia, which show a difference of 3,727 °C of DD5, or between Quercus grisea and P. strobiformis, with a difference of 47 days in SDAY (Table 2).

4 Discussion

We propose a robust approach based on the Weibull distribution for describing the relationship between environmental variables and tree species abundance in temperate forests. The two- and three-parameter Weibull functions both provided significant fits for the density of various tree species of temperate climate. The density curves helped determine points of maximum probability of abundance for these species.

According to the study results, W3p performed better than W2p in describing the density of most of the selected species by using the predictor variables, which may be explained by the inclusion of the location parameter in the former (Table 1). Theoretically, the location parameter contributes significantly to defining the distribution curve, because it has a definite integral (Borders et al. 1987). In the W3p approach, the location parameter indicates the point of departure of the tree density distribution along the abscissa of the respective environmental variable, and its position seems to be influenced by an uneven distribution of the climatic records that exist in the study area, which provides insights into the ecological significance of this parameter (Martínez-Antúnez et al. 2013).

In this study, the W2p provided the best description of the relationships for six of the tree species. Pece et al. (2000) showed that the W2p provided better fits than W3p for modelling the diameter distribution of Melia azedarach.

In relation to the location parameter, Knowe et al. (1997) indicated that this parameter does not make a significant contribution to the model if it is a negative value, and therefore it can be ignored or equated to zero to yield a two-parameter model. In the present study, the negative value of the location parameter did not appear to affect the fit of the model as some of the Dn values were very low, as in the following cases: Arbutus spp. for the EASL variable, with a Dn of 0.01; Quercus grisea for the SDAY variable, with a Dn = 0.05; Juniperus deppeana and P. strobiformis for EASL, with a Dn = 0.05.

There is evidence that negative values of the location parameter are compensated by high positive values of the other two parameters (i.e. presence of correlation between parameters), especially for values of the scale parameter. In this regard, Vallejos-Barra et al. (2009) observed a highly significant linear relationship between the location parameter and the other two parameters, particularly the scale parameter.

For several of the other tree species, a positive location parameter was obtained with most of the variables analysed. This suggests the need to estimate this parameter individually for each tree species in relation to each variable of interest and not to use pre-established values (Vallejos-Barra et al. 2009).

Unlike the location or the scale parameters, the shape parameter has a notable effect on the model because its value depends on the symmetry and kurtosis of the distribution of the function (Borders et al. 1987).

The findings of the present study indicate that there is not enough evidence to demonstrate that one of the functions is always better than the other for modelling the density of the species under study (Table 1). Thus, there are two main points to consider when deciding whether to use the two- or three-parameter Weibull distribution. If the data follow a curve and a value of α other than zero can be justified, then the three-parameter Weibull distribution is recommended. The high correlation between parameter estimates observed in other studies (Vallejos-Barra et al. 2009) suggests that it may often be preferable to use the two-parameter form of the model. However, if the three-parameter model is used, and high correlations are observed, the parameters estimates will obviously be more reliable.

It is difficult to determine clearly and accurately the bioclimatic niche of tree species, due to the complex interaction that occurs in the natural system (Soberón and Peterson 2005). However, determining the main predictor of species density and the values in which the maximum probability of abundance of a tree species occurs may represent an important advance, and this information can be used in planning and decision-making for the sustainable use and management of forest resources.

Most of the maximum values of tree species abundance modelled here with the Weibull density curve are within the ranges cited in previous studies for temperate forest communities (García-Arévalo and González-Elizondo 2003; González-Elizondo et al. 2012). According to model W3p, values of the EASL at which Pinus cooperi, Arbutus spp. and Juniperus deppeana occur at maximum probability of abundance are very similar. Also in MTWM, the maximum probability of abundance values coincided for Arbutus spp. and Pinus durangensis. The SDAY and GSP variables showed similar values for the same species.

In the study area, large amounts of coniferous trees are used for timber production, which may have affected the findings, especially for Pinus durangensis, Pinus cooperi, Pinus arizonica and Pinus engelmannii, which are of great commercial interest in the area (García-Arévalo and González-Elizondo 2003).

In relation to the species for which density was not modelled well by any of the Weibull functions tested, several factors may have had an effect:

• (1) The sample size may have been insufficient for some species (Kolb et al. 2006).

• (2) More than one maximum point of abundance will probably occur for some tree species; i.e. the data may have a multimodal distribution that cannot be described by the Weibull approaches tested. Some of the histograms indicated a decrease in the abundance of a particular species or absence of that species in relation to some of the variables analyzed (Figs. 2, 4 and 5): Alnus acuminata for EASL, MTCM and MTWM; Quercus magnolifolia for EASL and MTCM; Pinus cooperi for GSP; Juniperus deppeana for DD5; and Pinus engelmannii for SDAY.

• (3) The parameter estimation method may also have had an effect (Lockhart and Stephens 1994).

• (4) The effect of the variability on climatic data is also important (Saupe et al. 2012). The climatic data used here were obtained from a model output, and they may be imprecise due to the limited number of weather stations (183) scattered across the study area (6.33 million ha).

• (5) Previous studies detected correlation between the predictor variables studied (Martínez-Antúnez et al. 2014).

• (vi) In addition, some other factors were not included in this study; for example, the soil physical and chemical characteristics, the geographic accessibility and the adaptive capacities of the trees (Kolb et al. 2006).

The two- and three-parameter functions based on the Weibull distribution both yielded significant fits for the density of various tree species growing in a temperate climate. The density curves can be used to determine the points of maximum probability of abundance of these species. The three-parameter Weibull equation yielded better fits than the two-parameter equation in most cases, although, on a few occasions both were equally good and could have been used indistinctively.

In addition to the numerical results, the graphical results should also be analyzed to assess the behaviour of the model and to choose the best function. When the distribution curves of these functions are very similar, as shown in Figs. 2 and 3, the W2p function is preferred, since the parameters are relatively easy to estimate.

The relationships studied here and the definition of points of maximum probability of abundance, constitute an important basis for understanding species-specific adaptation to certain environmental conditions. The method presented here can be applied to any location, and in future studies it would be interesting to compare the W2p and W3p distributions with other models including more parameters, which may even yield multimodal frequency curves.

The models conventionally used to classify plant species, such as the generalized linear model (GLM), generalized additive models (GAMS) and canonical correspondence analysis (CCA) (Guisan and Zimmermann 2000), are important for studying populations and their environments. However, Martínez-Antúnez et al. (2014) reported that use of multivariate linear models to evaluate density as a function of climate variables yielded poor results.

We also suggest exploring alternatives to complement conventional methods (Phillips et al. 2006; Soberón and Nakamura 2009), e.g. use of the Gaussian mixture model (Rasmussen 1999) to represent optimal and sub-optimal ranges of abundance of forest species, using 0.75 and 0.25 quantiles as classification limits.

The maximum probability of abundance considered here may serve as a reference for constructing ecograms or ecomaps to provide a multidimensional perspective of the fundamental requirements for each species.