A methodology for assessing departure of current plant communities from historical conditions over large landscapes

Abstract

Fire frequency and severity, and vegetation composition and structure have been altered across much of North America during the past century because of fire exclusion and other land management practices. The cumulative results are now recognized to be partly responsible for dramatic increases in wildland fire severity and declines in ecosystem health. In response, the Departments of Agriculture and Interior are developing a strategy for reducing fire hazard and restoring forest health. A key component of this strategy is the assessment of the departure of current plant communities from historical conditions. Assessing departure is difficult because of limited spatial coverage of data on successional class distribution prior to European settlement (historical) conditions. This article discusses a two-step approach to the problem that first generates data on historical vegetation composition and structure by simulating succession and disturbance processes under historical conditions, and then measures departure by comparing these simulated data to an observation on current conditions. The data are observations on high-dimensional multinomial categorical variables, and pose several problems that limit the usefulness of conventional statistical methods. We propose a method that constructs a linear approximation of the current observation vector from the simulated historic data, and measures the departure of the current observation vector by the length of the residual error vector. A simulation study indicates that this departure measure is substantially more sensitive than conventional outlier detection methods. Our methodology is demonstrated using a pilot region encompassing lands in Utah.

Introduction

Many forests and rangelands of North America experienced frequent, low severity fires prior to European settlement (circa 1900) Wright and Heinselman, 1973, Heinselman, 1981, DeBano et al., 1998. Since European settlement, the frequency and severity of fire, vegetation composition and structure, and fuel loading have changed substantially (Ferry et al., 1995, Frost, 1998, Kolb et al., 1998, Arno et al., 2000, Keane et al., 2002b). Fire exclusion, livestock grazing, and logging, and generally, land management policies are widely attributed to these changes. A probable consequence of these changes are dramatic increases in the number, size, and intensity of wildfires in recent years (GAO, 2002). In response, Congress and the Executive Branch mandated the development of a National Fire Plan (www.fireplan.gov) for reducing fire hazard and restoring forest health. Through this vehicle, the Departments of Agriculture and Interior have been directed by Congress to develop a cohesive strategy for implementing the Plan (Laverty and Williams, 2000). These efforts, along with the passage of the Healthy Forests Restoration Act (www.healthyforests.gov), provide support for reducing fire hazard using silvicultural thinning and prescribed burning.

The Cohesive Strategy uses fire regime condition class as a key measure for implementing the National Fire Plan. Fire regime condition class measures the departure of current conditions from historical conditions with respect to fuel, fire regime, and vegetation (Laverty and Williams, 2000), and is used in part to direct funding and resources to those lands with the greatest need for restoration. However, assessing departure is difficult because data on historical conditions are extremely sparse (Keane et al., 2002a). Several strategies have been identified for quantifying local historical conditions (Landres et al., 1999), but the only practical strategy for assessing departure on a national scale is to generate historical data by spatially simulating plant community succession as disturbances across landscapes Keane et al., 2002a, Landres et al., 1999. The use of simulated data in ecology and resource management has become commonplace for large scale problems (Waring and Running, 1998). Examples include the analysis of net primary production (Jiang et al., 1999), forest planning (Bettinger et al., 2004), and nutrient cycling (Nour et al., 2006). In particular, succession and fire modeling has been an active area of ecological research for 3 decades (see Keane et al., 2004; Mladendoff and Baker, 1999). These models are now mature in their sophistication and accuracy, and limitations arising from model assumptions are widely understood Botkin and Schenk, 1996, Keane et al., 2001. The current generation of successional models are spatially explicit, and model fire events as a function of topography, fuel, weather and prevailing wind direction (Keane et al., 2004). The subject of this article is the comparison of simulated historical and current observations, and departure assessment. The simulated data pose several statistical challenges common to data generated by automated collection devices and other recent technological innovations (Rao, 2004). These challenges include high-dimensional observations, massive data sets, and a lack of appropriate distributional models.

This article presents a statistical method of assessing the departure of current vegetation conditions from simulated historical conditions for each map unit in a landscape. The essence of the method is the comparison of an observation on the current distribution of successional classes to a simulated set of observations on the historical distribution of successional classes. The population of interest consists of the landscape units within geographic region, and an observation vector consists of the areal covers of each possible successional class within each observed potential vegetation class. Areal cover is quantified by the number of pixels within a landscape unit belonging to a particular potential vegetation class and successional class combination. As most of these class combinations are rare, or not present at a particular point in time, the data are observations on high-dimensional multinomial categorical variables dominated by zero counts. The tendency for many classes to be dominated by small counts implies that many conventional distribution-based statistical methods [e.g., the Chi-square test for homogeneity of proportions (Ott and Longnecker, 2001, Chapter 10)] are inappropriate for these data. Moreover, the historical data sets contain multiple linear dependencies, and so the data matrices are less-than-full rank. Most multivariate methods are based on using full rank data matrices and must be modified to accommodate a less-than-full rank data matrix.

These considerations motivate our approach. We construct a linear approximation of the current data vector using the historical data, and measure the departure of current from historical conditions by the approximation error. This is accomplished by using the singular value decomposition of the historic data matrix $\text{X}$ to construct a projection matrix for the space spanned by the rows of $\text{X}$. Then, the departure of an observation vector $X_0$ is defined to be the length of the error, or residual vector associated with the projection of $X_0$ onto the space spanned by the rows of $\text{X}$. To determine the observed significance level, an empirical distribution of departure values is constructed by removing each simulated historical observation from $\text{X}$ and computing its departure from the remaining data.

Many federal land management agencies use a two-tier hierarchical plant community classification system Anderson et al., 1998, Grossman et al., 1998 for describing large scale landscapes (for example, see Hessburg et al., 1999, Hann et al., 1998). The first tier, called the potential vegetation type, describes the biophysical conditions at a site by identifying the plant community that would exist if succession were to progress to the stable or climax community Daubenmire, 1966, Pfister and Arno, 1980. The second tier of the classification system identifies the existing successional class at a site. Succession is the progression of plant communities from colonization to climax community. Usually, succession is described as a series of classes, and transition between classes is viewed as a deterministic process governed by species-specific environmental requirements, inter-species competition, and environmental conditions. In the plant community classification system, successional class is nested within potential vegetation type, and identifies the vertical structure and dominant overstory species. Fig. 1 shows an example of the successional pathways for a subalpine forest potential vegetation type found in Utah.

Disturbance events such as wildland fire interrupt succession and often modify the subsequent successional progression, or pathway. The resulting successional class and pathway depends on the disturbance severity. For example, a low intensity fire in a mid-successional forest class may remove understory trees leaving a single-strata overstory tree canopy that will progress to a two-tiered canopy, whereas a high intensity fire may produce a grass and shrub community that will progress towards a single-strata canopy. Low intensity fires tend to perpetuate some non-climax successional classes. For example, frequent fires in a grass and shrub community will prevent the establishment of an overstory, thereby maintaining the grass and shrub community. Mid-successional forests are also maintained by low intensity fires that eliminate climax tree species from the understory. On the other hand, early-successional forests comprised of saplings and poles may be set back to grass and shrub communities by low intensity fire. Recognition of the prevalence of multiple successional pathways implies in a natural ecosystem, the distribution of successional classes is within the historical range of variation Keane et al., 2002a, Landres et al., 1999. Before proceeding further, it is convenient to establish notation and terminology.

Let $X_i$ denote a vector those elements are the frequency of occurrence for p successional classes for some landscape unit at time i. There are n such vectors sampled from a simulation of historical conditions, and these are indexed by $i=1\text{,}\dots \text{,}n$. The vector $X_0$ is an observation on the current successional class distribution. In this article, a landscape unit is comprised of 900, $30\text{m}\times 30$  m pixels, each of which belongs to a single successional class at time i. Consequently, the successional class vector $X_i$ consists of the number of pixels in each of the successional classes, subject to the constraint ${\sum }_j{X}_{ij}=900$. The data, historical and current, are collected as an $n+1\times p$ matrix

$$\text{X}=\left(\begin{array}{c}\hfill X_0^{\text{T}}\hfill \\ \hfill X_1^{\text{T}}\hfill \\ \hfill ⋮\hfill \\ \hfill X_n^{\text{T}}\hfill \end{array}\right)\text{,}$$

where $X_i^{\text{T}}$ denotes the transpose of $X_i$. The i, jth element of $\text{X}$ is denoted by $x_{ij}$. The jth column mean is denoted by ${\overline{x}}_{\cdot j}=n^{-1}{\sum }_{i=1}^n{x}_{ij}$, and the p-vector of column means is $\overline{X}=({\overline{x}}_{\cdot 1}\text{,}\dots \text{,}{\overline{x}}_{\cdot p})$. The $L^2$ norm of the vector X is $\Vert X{\Vert }_2=\sqrt{X^{\text{T}}X}$. The $p\times p$ sample variance matrix is $\text{D}$, and its spectral decomposition is $\text{D}=\text{V}\Lambda {\text{V}}^{\text{T}}$, where $\Lambda =\text{diag}({\lambda }_1\text{,}{\lambda }_2\text{,}\dots \text{,}{\lambda }_p)$ is a diagonal matrix comprised of the eigenvalues of $\text{D}$ with the eigenvalues arranged in descending order, i.e., ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_p$ (Healy, 1986, Chapter 7). The matrix $\text{V}=(V_1\text{,}{V}_2\text{,}\dots \text{,}{V}_p)$ is comprised of the normalized eigenvectors of $\text{D}$ arranged in concordance with the ordering on the eigenvalues. For the problem at hand, $\text{X}$ is not full rank and $\text{D}$ is singular because there are multiple constraints on the columns of $\text{X}$. Suppose that the rank of $\text{X}$ is $r\le \min (n+1\text{,}p)$. Let ${\Lambda }_r$ denote the square submatrix of $\Lambda$ containing the r positive eigenvalues and ${\text{V}}_r=(V_1\text{,}\dots \text{,}{V}_r)$ denote the submatrix of $\text{V}$ corresponding to the positive eigenvalues. The Moore–Penrose generalized inverse of $\text{D}$ is

$${\text{D}}^{-}={\text{V}}_r{\Lambda }_r^{-1}{\text{V}}_r^{\text{T}}$$

(Healy, 1986, Chapter 5; Schott, 1997, Chapter 5). Because ${\lambda }_i>0$ for $i=1\text{,}\dots \text{,}r$, ${\text{D}}^{-}$ is positive definite and can be used to define a vector norm $\Vert X_0{\Vert }_{{\text{D}}^{-}}=\sqrt{X_0^{\text{T}}{\text{D}}^{-}{X}_0}=\Vert {\Lambda }_r^{-1/2}{\text{V}}_r^{\text{T}}{X}_0{\Vert }_2$, where ${\Lambda }_r^{-1/2}$ is the diagonal matrix with diagonal elements ${\lambda }_i^{-1/2}$, and a distance function $\rho (X\text{,}Y)=\Vert X-Y{\Vert }_{{\text{D}}^{-}}$. Finally, the matrix constructed by removing the ith row, or observation, from x is denoted by ${\text{X}}_i$.