Statistical modelling of a new global potential vegetation distribution

The 'potential natural vegetation' (PNV) can be seen as the natural vegetation, i.e., in equilibrium with climate, that would exist at a given location non-impacted by human activities. A global PNV distribution is required for many purposes in environmental sciences. Examples include estimating historical changes of land-use (Ramankutty and Foley 1999), the impact of atmospheric CO₂ concentration on vegetation (Cha 1997, Notaro et al 2005), the response of vegetation to climate changes (Ni et al 2006, Notaro 2008) or palaeovegetation distributions (Crucifix et al 2005, Woillez et al 2011).

Most of these studies refer to the global PNV map derived by Ramankutty and Foley (1999) (hereafter 'RF99') from remotely-sensed observations and corrected in human-impacted regions using the vegetation model BIOME3 (Haxeltine and Prentice 1996). Other applications (e.g., palaeovegetation modelling) directly use vegetation models to simulate a PNV distribution. These two methods include biases and uncertainties from vegetation models.

In this context, Levavasseur et al (2012) (hereafter 'L12') described a new methodology to statistically model a high-resolution PNV distribution over Europe, entirely based on vegetation and climatological data. Their approach consists in using a multinomial logistic regression (MLR). MLR builds statistical relationships between vegetation data and climatological variables independently from models and without any subjective corrections. These relationships allow one to model the occurrence probabilities of each PNV type. L12 show good results over Europe in comparison with the map from RF99 or the PNV simulated by the vegetation model BIOME4 (Kaplan et al 2003). Moreover, the use of occurrence probabilities provides useful additional information as PNV fractions or uncertainty index.

In the present work we extend the L12 method to model a high-resolution gridded world-wide PNV distribution. We summarize the L12 framework and highlight the adaptations and updates that we performed in section 2. Then, we compare the global PNV modelled by MLR to the RF99 map (section 3.3) and to the PNV distribution simulated by the vegetation model BIOME4 over the globe (section 3.2). Discussions and conclusions follow in section 4.

2.1. The multinomial logistic regression (MLR)

To predict PNV types distribution, L12 used a multinomial logistic regression (MLR, Hosmer and Lemeshow 2000, Hilbe 2009). MLR builds statistical relationships between a nominal explained variable (called the predictand, the vegetation type in our case) and continuous explanatory variables (called the predictors, the climatic variables). Those relationships allow one to estimate the occurrence probability of the nominal explained variable (Y, the PNV types in our case), taking into account p numerical explanatory variables (X_k, the climatological variables):

$$\begin{eqnarray} &&\log \left (\frac{P({Y}_{i}=j)}{P({Y}_{i}=r)}\right )={\beta }_{0,j}+\sum _{k=1}^{p}{\beta }_{k,j}{X}_{i,k}^{\ast },\tqs \forall ~j\not = r, \end{eqnarray} \tag{ 1 }$$

where ${X}_{i,k}^{\ast }=\frac{{X}_{i,k}-{\mu }_{k}}{{\sigma }_{k}}$ are the standardized predictors, with μ_k the mean of the kth predictor and σ_k its standard deviation. P(Y_i = j) is the occurrence probability of the jth PNV type and i is the grid-cell. β_0,j is the intercept for the jth PNV type and β_k,j is the regression coefficient for the kth predictor and the jth PNV type. p is the number of predictors.

The predictors are standardized to obtain comparable regression coefficients (β_k,j) without units. According to their weights, the predictors will be ranked and selected in section 2.4.

Equation (1) is based on a reference category r: the desert vegetation in our case (defined in section 2.2.2). We obtain j − 1 relationships and the occurrence probabilities of the reference PNV type can be deduced in each grid-cell i with ${\mathop{\sum }\nolimits }_{j=1}^{m}P({Y}_{i}=j)=1$ (considering m PNV types including r). MLR is performed with the R package 'VGAM' (Yee and Wild 1996, Yee 2010a, 2010b) and parameters are estimated through likelihood maximization.

2.2. Vegetation data

2.2.1. The BIOME 6000 database.

As predictand for MLR, we use the BIOME 6000 database (Prentice and Jolly 2000, Harrison et al 2001, Bigelow et al 2003, Pickett et al 2004) from the Global Palaeovegetation Mapping Project³ for the modern period (i.e., 0 ka). BIOME 6000 compiles pollen and plant macrofossil data that are heterogeneously distributed over the world. Most data points are concentrated in the northern hemisphere between 30°N and 70°N, and no data point covers South America or India. Our region of interest is the globe without Antarctica, from 180°W to 180°E and from 60°S to 90°N, as shown in figure 1(a).

Zoom In Zoom Out Reset image size

BIOME 6000 data points are expressed in biomes following a 'biomization' method described in Prentice et al (1996). A biome includes characteristic vegetation types deduced from pollen under similar climatic conditions. BIOME 6000 can be classified into eight 'megabiomes': boreal forest, desert, tundra, grassland and dry shrubland, savanna and dry woodland, temperate forest, tropical forest and warm-temperate forest.

According to the authors of the database, some BIOME 6000 data points appear inconsistent in mountain areas due to pollen transport (Guiot 2012). For instance, the Alps or Pyrenees are essentially dominated by temperate forest pollens, even at high altitude. To correct for such discrepancies, we reclassify some BIOME 6000 data points using the growing degree day at 5 °C (GDD5) limits from Prentice et al (1992) (see figure 1(a)). Annual GDD5 corresponds to the sum of daily temperatures above 5 °C during a year. As no gridded dataset of daily temperature covers the world with a fine enough spatial resolution, we first compute a global GDD5 climatology from the NCEP/NCAR⁴ air temperature at the surface reanalysis daily time-series (between 1961 and 1990) on a 2.5° grid (Kalnay et al 1996). Then, we apply a statistical downscaling method based on the use of a Generalized Additive Model (GAM, Vrac et al 2007) to obtain a global GDD5 (cf appendix A) at a resolution of 10' (i.e., 1/6° in longitude and latitude, the final resolution of our map). This procedure reclassified 221 BIOME 6000 data points (i.e., 3.6% of the 6091 points) essentially localized on mountains (Alps, Himalaya, Rocky Mountains) and in boreal regions (Alaska, the coast of Labrador and Siberia).

2.2.2. 'True' deserts.

We define true deserts as sterile land areas without vegetation producing pollens. Consequently, no BIOME 6000 data point covers desert regions such as the Sahara or the Greenland ice-sheet. The desert megabiome from BIOME 6000 refers more specifically to desert vegetation dominated by sparse steppe forb and grass (Prentice et al 1996) instead of 'true' desert. Accordingly, we rename the desert megabiome from BIOME 6000 into a separated desert vegetation megabiome.

In order to represent real deserts, we build two new 'pseudo-megabiomes' (warm and cold deserts) from the IGBP-DIS land cover map⁵ (Loveland and Belward 1997). This dataset derived 17 land cover types from remotely-sensed observations between 1992 and 1993. Assuming a limited human-induced desertification, we manually add training data points as follows (see figure 1(b)):

• To reflect cold deserts, 100 points have been randomly evenly distributed over the Greenland ice-sheet and where Loveland and Belward (1997) shows polar, rock or ice deserts. 100 points are enough to homogeneously cover cold desert regions.
• Over the globe, the total warm desert area is approximately five times larger than the total cold desert area. To be consistent and reflect warm deserts, 500 data points have been distributed in the same way as where Loveland and Belward (1997) shows warm deserts.

2.3. The explanatory variables

2.4. Model selection

3.1. Comparison MLR versus BIOME 6000

3.2. Comparison MLR versus BIOME4

To ascertain our method, we directly confront the modelled PNV by MLR with the simulated vegetation from a vegetation model. The BIOME4 model (Haxeltine and Prentice 1996, Kaplan et al 2003) is driven by temperature, sunshine and precipitation monthly climatologies from the CRU database (described in section 2.3). To be consistent with the period represented by CRU climatologies (around 1980), the atmospheric CO₂ concentration is set for 360 ppm (Lüthi et al 2008). BIOME 4 has a biome scale easily translatable into our megabiomes following Harrison and Prentice (2003).

Figures 2(a) and (b) respectively show the modelled PNV by MLR in each grid-cell of our map (at 10' resolution) and the simulated PNV distribution by BIOME4. Both maps show large similarities especially concerning the distribution of tundra, temperate and boreal forests at high latitudes. Note that our BIOME 4 simulation does not show warm-temperate forests around Mediterranean Sea, in southeastern China and USA and in eastern Australia. Modelling warm-temperate forests by MLR in these regions is in agreement with BIOME 6000 database (see figure 1(a) and Levavasseur et al 2012) and with older published BIOME simulations (Prentice et al 1992, 1996, Harrison and Prentice 2003, Tang et al 2009).

Zoom In Zoom Out Reset image size

Nevertheless, some mountain areas are not well defined (e.g., tundra and boreal forests of Tian Mountains are replaced by desert vegetation) or even disappear (e.g., the Andes or US Rockies) with MLR. Differences appears in the western US, where MLR models a drier vegetation than BIOME4. Moreover, MLR models larger warm deserts than BIOME4.

Although both methods are based on CRU climatologies, BIOME4 computes some mechanistic processes (i.e., physiology, competitiveness or productivity) which may induce some of the difference.

3.3. Comparison MLR versus RF99

We compared in this paper several potential natural vegetation (PNV) distributions over the globe based on different methods:

• The PNV map from Ramankutty and Foley (1999) (RF99) built from remotely-sensed data (Loveland et al 2000) and from the vegetation model BIOME3 (Haxeltine and Prentice 1996).
• PNV simulated by the vegetation model BIOME4 (Kaplan et al 2003) driven by the CRU climatologies.
• A new high-resolution global PNV map built from multinomial logistic models.

Obvious similarities appear between reconstructions, especially with the establishment of tundra, temperate and boreal forests at high latitudes.

The vegetation model BIOME4 is partly calibrated to represent BIOME 6000 data (Kaplan et al 2003). The κ coefficient computed between BIOME4 and BIOME 6000 (appendix B) is of 0.43 over the world, while MLR obtains a κ of 0.64 (section 3.1). Despite including no physical or mechanistic processes, MLR obtains a realistic PNV distribution closer to BIOME 6000 data. A κ coefficient computed between RF99 and BIOME 6000 is not possible given the different biomes vegetation types. The RF99 map appears more heterogeneous than MLR. Some of these local details are disputable and correspond more to current observed vegetation rather than potential vegetation. For example, RF99 sees the Landes forests (in the region of Landes in France) which have been mainly planted by humans. MLR sees also a dominance of forests, but with large uncertainty.

Indeed, MLR does not only provide a vegetation distribution, since we obtain an occurrence probability by megabiome. Occurrence probabilities allows us to estimate the uncertainty of the modelled PNV distribution, taking into account the megabiome with second highest occurrence probability. This second dominant megabiome often appears in agreement in regions where the first dominant megabiome is different from other databases. The percentage of agreement with BIOME 6000 increases from 69.5% (see section 3.1) to 89.9%, taking into account the second dominant megabiome. In agreement with RF99 or BIOME4, the map of the second dominant megabiome modelled by MLR (not shown) shows tundra and boreal forests in northeastern Europe (section 3.3), tropical forest in equatorial region and cold desert in Andes capturing the effect of local-scale topography through the high-resolution CRU climatologies. An uncertainty index UI can also be computed from the difference between the two highest occurrence probabilities (figure 3):

$$\begin{eqnarray} &&\mathrm{UI}=\frac{1-({p}_{1}-{p}_{2})}{{p}_{1}} \end{eqnarray} \tag{ 2 }$$

where p_x is the occurrence probability and x the rank of the probability ranging from 1 (the highest probability) to m (the lowest probability) with m the number of megabiomes.

Zoom In Zoom Out Reset image size

A first dominant megabiome with a probability close to the probability of the second dominant megabiome has an uncertainty close to 1 and vice versa. This index appears very useful in bringing some confidence in the PNV modelled by MLR and pointing out the limits of our method:

• This uncertainty index allows us to target the regions needing some data to improve the PNV distribution, such as in South America. Indeed, the main limit of the MLR method lies in the training data (BIOME 6000 in this case). The modelled PNV by MLR highly depends on the abundance and geographical distribution of data points. If a megabiome is absent or over/under-represented, this will have a significant impact on the modelled PNV by MLR. Nevertheless, a calibration of MLR over the globe provides a geographical robustness to the statistical model. The PNV predicted by MLR in regions with no or less BIOME 6000 data appears consistent with climatic patterns (e.g., MLR shows similarities with modern biome reconstructions from Marchant et al 2009 in several regions of South America).
• This index highlights the regions where the modelled vegetation is to be taken with caution (as in western US) because the climatic signal alone is not sufficient to distinguish a dominant megabiome. The role of fires, herbivory or soil properties has been clearly identified for vegetation (e.g., Sankaran et al 2005 for savannas) and could be included as predictors in our statistical model.

Like any database, the disadvantages of our statistical approach should be discussed to better constrain its application. The 'vegetation–climate' relationship estimated by MLR from BIOME 6000 modern data is implicitly constrained by an atmospheric CO₂ concentration of about 360 ppm. As a prospect, exporting this relationship in different climatic conditions leads to a distribution ignoring the crucial effect of CO₂ on vegetation (Harrison and Prentice 2003, Woillez et al 2011). Moreover, MLR do not simulate soil–vegetation–atmosphere interactions such as photosynthesis, growth and competitiveness of plants, which may be more constant at the biome level. Vegetation models allow us to provide characteristics of vegetation as leaf area index (LAI) or net primary productivity (NPP). Statistical modelling of vegetation is an interesting and complementary alternative to process-based vegetation models.

Finally, the PNV modelled by MLR cannot claim to be fully independent of human influences. MLR is mainly based on climatological data between 1961 and 1990, impacted by human activities through climate change. Moreover, BIOME 6000 data includes modern data referring to samples dated within the past thousand years (most of pollen samples falls within the past 500 years (Bigelow et al 2003)). Man has intensively used lands for thousands of years (for example, in ancient Greece, during the Western agricultural revolution in the Middle Ages or more recently with the Green Revolution between 1960 and 1980). To warrant the 'potential' feature of the modelled vegetation by MLR, it could be relevant to calibrate MLR on BIOME 6000 data from the Holocene (−6 ka). At this period, the land-use was limited to a few scattered subtropical farm households (e.g., in China or South America).

For details about these last remarks, the method and the used data, we refer the reader to Levavasseur et al (2012). All final data (megabiomes and occurrence probabilities) are in an attached supplementary NetCDF file (available at stacks.iop.org/ERL/7/044019/mmedia).

To conclude, for the modern period, BIOME 6000 can be confidently considered as reference data collected in areas with less possible human activity, although ensuring data not impacted by humans is difficult. Accounting all our observations and statistical indices, MLR models the most realistic PNV on the regions covered by BIOME 6000. Over the rest of the world, MLR models a vegetation distribution consistent with climatic signal. The MLR method is a fast and robust alternative in vegetation modelling with several advantages. The modelled PNV map is (i) directly and only based on vegetation (BIOME 6000) and climatological (CRU) data; (ii) not subjective and independent of any vegetation model; (iii) easily updatable as soon as additional data is made available.

We acknowledge the E-OBS dataset from the EU-FP6 project ENSEMBLES (http://ensembles-eu.metoffice.com) and the data providers in the ECA&D project (http://eca.knmi.nl). G Levavasseur is supported by UVSQ, D M Roche and M Vrac by INSU/CNRS.

As GDD5 is built from temperatures, we choose to apply the statistical downscaling method developed by Vrac et al (2007) for temperature downscaling and based on the use of a Generalized Additive Model (GAM, Hastie and Tibshirani 1990). GAM models statistical relationships between local-scale observations over Europe: the high-resolution GDD5 climatology built from the ECA&D times-series (Haylock et al 2008); and global variables: the topography from ETOPO2 (Amante and Eakins 2009), the temperature from CRU (New et al 2002), and the GDD5 built from NCEP/NCAR (see section 2.2.1). GAM represents the expectation of an explained variable Y (the predictand, i.e., the GDD5 from ECA&D in our case) by a sum of nonlinear functions f_k, conditionally on explanatory variables (the predictors, i.e., the topography, the temperature and the GDD5 from NCEP/NCAR) X_k (Hastie and Tibshirani 1990):

$$\begin{eqnarray} &&E({Y}_{i}\vert {X}_{k,k=1\cdots n})=\sum _{k=1}^{n}{f}_{k}({X}_{i,k})+\epsilon , \end{eqnarray} \tag{ A.1 }$$

where is the residual or error, n is the number of predictors and i is the grid-cell. To use GAM, we need to define the distribution family of the explained variable. For simplicity, Vrac et al (2007) assumed that temperature has a Gaussian distribution so we assume that GDD5 too, which implies a zero-mean Gaussian error (Hastie and Tibshirani 1990). Then, we define the nonlinear functions as cubic regression splines (piecewise by third-degree polynomials). Finally, any SDM needs a calibration/projection procedure. The calibration is the fitting process of the splines over Europe in our case. Afterwards, we project over the world to predict a high-resolution global GDD5 climatology.

Instead of a simple bilinear interpolation of the NCEP/NCAR GDD5, we use GAM to geographically extrapolate the characteristics of the ECA&D GDD5 over Europe to the world. For more details we refer the reader to Vrac et al (2007) and Martin et al (2012). We perform this analysis within the statistical programming environment R (R Development Core Team 2011) and its 'mgcv' package (Wood 2006).

The Bayesian Information Criterion (BIC). The BIC (equation (B.1)) is a particular form of the Akaike Information Criterion (Sakamoto et al 1986), developed by Schwarz (1978) and defined by:

$$\begin{eqnarray} &&\mathrm{BIC}=-2\times \mathrm{LL}+P\times \log (n) \end{eqnarray} \tag{ B.1 }$$

where n corresponds to the number of BIOME 6000 data points (n = 6091), P is the number of parameters in the fitted model (P = n × (m − 1)) and LL is the log-likelihood of the fitted model. This criterion measures the goodness-of-fit between the statistical model and the data, balancing the risk of over-fitting. The BIC includes a penalty term depending on the sample size (n) and on the dimension of the model (P). The smaller the BIC, the better the model.

Pseudo-R². The R² is a classical statistical index in ordinary least squares regression that is often used as a goodness-of-fit measure. In logistic regression, an equivalent statistic to R² does not exist. However, to evaluate the goodness-of-fit of logistic models, several 'pseudo-R²' (ranging from 0 to 1) have been proposed. Among the different approaches, the McFadden's pseudo-R² is often used for its simplicity of calculation (equation (B.2)) and interpretation. It is defined by (Menard 2000, Shtatland et al 2002):

$$\begin{eqnarray} &&{R}^{2}=1-\frac{\mathrm{LL}}{{\mathrm{LL}}_{\mathrm{null}}} \end{eqnarray} \tag{ B.2 }$$

where LL is the log-likelihood of the selected model (i.e., with selected predictors) and LL_null the log-likelihood for the null-model (i.e., with intercept only). The ratio of log-likelihoods suggests the level of improvement over the null-model offered by the involved predictors. A small ratio of likelihoods indicates that the full model is far better than the null-model. In terms of pseudo-R², the closer the R² is to 1, the better the agreement with data is.

The kappa statistic. The κ coefficient measures the quality of the agreement (Cohen 1960, Fleiss et al 1969) between the modelled PNV by MLR in each BIOME 6000 location (figure 1(b)) and the BIOME 6000 data (figure 1(a)). This index can take values between 0 and 1 and is based on a simple counting of matching and non-matching points in a matrix used to represent errors in assigning classes (see appendix A of Levavasseur et al (2011)). The closer the κ coefficient is to 1, the better the agreement with data is. The kappa statistic is often used for spatial comparison of categorical variables, such as vegetation (Monserud and Leemans 1992).

The Brier score. The Brier score was developed by Brier 1950 to assess the accuracy of probabilistic forecasts. As MLR provides probabilities of occurrence of different megabiomes, this score is well adapted here. It measures the average squared deviation between predicted probabilities for a set of events and their binary outcomes (0 if the event does not happen and 1 if it happens). For a multinomial variable, the Brier score is defined by:

$$\begin{eqnarray} &&\mathrm{BS}=\frac{1}{n}\sum _{j=1}^{m}\sum _{i=1}^{n}({p}_{i,j}-{o}_{i,j})^{2}\in ~[0;~m] \end{eqnarray} \tag{ B.3 }$$

where n is the number of BIOME 6000 data points and m is the number of megabiomes. p_i,j corresponds to the predict probability of the jth megabiomes at the ith point/location and o_i,j is the corresponding binary outcome for this point. The Brier score can take values between 0 and m. A lower score represents higher accuracy of the prediction. The Brier score can also be reduced in two other ways:

• Taking into account all m megabiomes by grid-cell/location, we obtain a map of Brier scores:
  $$\begin{eqnarray} &&\mathrm{BS}(i)=\sum _{j=1}^{m}({p}_{j}-{o}_{j})^{2}\in ~[0;~m]. \end{eqnarray} \tag{ B.4 }$$
• Taking into account all n grid-cells/locations by megabiome, we obtain m Brier scores (i.e., for each j megabiome):
  $$\begin{eqnarray} &&\mathrm{BS}(j)=\frac{1}{n}\sum _{i=1}^{n}({p}_{i}-{o}_{i})^{2}\in ~[0;~1]. \end{eqnarray} \tag{ B.5 }$$