Development of Lichen 

Response Indexes Using a 

Regional Gradient Modeling 

Approach for Large-Scale 

Monitoring of Forests

Susan Will-Wolf and Peter Neitlich

United States 

Department of 

Agriculture

Forest Service

Pacific Northwest 

Research Station

General Technical 

Report

PNW-GTR-807

September 2010

FIA lichen model region

Pacific Northwest West

Pacific Northwest East

Northern Rockies

Northwestern California

Sierra–Southern Cascades–Modoc

Great Basin–Semidesert

Greater Central Valley

Southern California Coast

Southwest

Middle and Southern Rockies

North Central

Northeast

Ohio Valley

Mid Atlantic

South Central

Southeast

Southern Florida

Boundaries

States

Ecoregion sections

Alaska Arctic and

Interior

Southeast

and South

Central

Southwest

and Aleutian

Hawaii
The Forest Service of the U.S. Department of Agriculture is dedicated to the principle of 

multiple use management of the Nation’s forest resources for sustained yields of wood, 

water, forage, wildlife, and recreation. Through forestry research, cooperation with the 

States and private forest owners, and management of the National Forests and National 

Grasslands, it strives—as directed by Congress—to provide increasingly greater service to a 

growing Nation.

The U.S. Department of Agriculture (USDA) prohibits discrimination in all its programs and 

activities on the basis of race, color, national origin, age, disability, and where applicable, sex, 

marital status, familial status, parental status, religion, sexual orientation, genetic information, 

political beliefs, reprisal, or because all or part of an individual’s income is derived from any 

public assistance program. (Not all prohibited bases apply to all programs.) Persons with 

disabilities who require alternative means for communication of program information (Braille, 

large print, audiotape, etc.) should contact USDA’s TARGET Center at (202) 720-2600 (voice 

and TDD). To file a complaint of discrimination, write USDA, Director, Office of Civil Rights, 

1400 Independence Avenue, SW, Washington, DC 20250-9410 or call (800) 795-3272 (voice) 

or (202) 720-6382 (TDD). USDA is an equal opportunity provider and employer.

Authors

Susan Will-Wolf is a research ecologist, Department of Botany, University 

of Wisconsin-Madison, 430 Lincoln Drive, Madison, WI 53706-1381 (Fax: 

608-262-7509, Phone: 608-262-2754, E-mail: swwolf@wisc.edu). Peter Neitlich

is a research ecologist, National Park Service, 41A Wandling Road, Winthrop, 

WA 98862 (Phone: 509-996-8031, E-mails: Peter_Neitlich@ccmail.itd.nps.gov, 

pneitlich@methownet.com).

This work was performed under cooperative agreements SRS 06-CA-11330145-079 

and SRS 09-CA-11330145-101 with the Forest Inventory and Analysis Program for 

the senior author. The second author is a contractor under the first agreement.
Abstract

Will-Wolf, Susan; Neitlich, Peter. 2010. Development of lichen response indexes 

using a regional gradient modeling approach for large-scale monitoring of 

forests. Gen. Tech. Rep. PNW-GTR-807. Portland, OR: U.S. Department of 

Agriculture, Forest Service, Pacific Northwest Research Station. 65 p.

Development of a regional lichen gradient model from community data is a 

powerful tool to derive lichen indexes of response to environmental factors for 

large-scale and long-term monitoring of forest ecosystems. The Forest Inventory 

and Analysis (FIA) Program of the U.S. Department of Agriculture Forest Service 

includes lichens in its national inventory of forests of the United States, to help 

monitor the status of forested ecosystems. Development of a model for a specific 

region to calculate lichen response indexes that are correlated with air quality and 

major climate factors, and are also independent of each other, is a critical step in 

achieving program goals. These indexes are the primary lichen bioindicators used 

in FIA for assessing regional patterns and monitoring trends of lichen response to 

environment over time. This general approach is also applicable to other monitoring 

efforts. A first step in the modeling process is to identify an appropriate geographic 

region for a model. Unconstrained ordination alone, or combined with indicator 

species analysis followed by regression analysis, are two approaches borrowed from 

plant ecology that have been shown to generate successful regional lichen gradient 

models. Calculation of lichen response indexes for new plots not part of the original 

model is necessary to support long-term monitoring. We explain the rationale for 

recommended approaches, describe in detail the recommended steps in the modeldevelopment process, and explain how to document and evaluate results, all to 

support successful application of a model for monitoring. A template is included for 

documenting a model and archiving all products necessary to understand and apply 

it, as is required for each FIA model.

Keywords: Air pollution, air quality, biomonitor, climate, environmental 

response index, forest health, lichen, community, ordination.
Summary

Development of a regional lichen gradient model from community data is a 

powerful tool to derive lichen indexes of response to environmental factors for 

large-scale and long-term monitoring of forest ecosystems. The Forest Inventory 

and Analysis (FIA) Program of the U.S. Department of Agriculture Forest Service 

includes lichens in its national inventory of forests of the United States for largescale monitoring of the status of forested ecosystems with tight budget constraints. 

The purpose of the FIA Program is to inventory forests of the Nation to facilitate 

responsible use, conservation, and maintenance of function for forested ecosystems. 

A critical step to achieve program goals for the lichen indicator is development of 

lichen response indexes for a region that are correlated with air quality and major 

climate factors, and are also independent of each other. These indexes are the primary lichen bioindicators used in FIA for assessing regional patterns and monitoring trends over time. This general approach is also applicable to other monitoring 

efforts. 

A first step in the modeling process is to identify an appropriate geographic 

region for a model. A region must be large to support models that will remain useful for many years and response indexes that are sensitive primarily to large-scale 

environmental factors. Upper size is constrained by the need to limit the variation 

of the sampled lichen communities across a region to support a stable model. A 

region must be compact to allow linkage of change in lichen response indexes to 

lichen dispersal. 

Unconstrained ordination alone, or combined with indicator species analysis 

followed by regression analysis, are two approaches borrowed from plant ecology 

that have been shown to generate successful regional lichen gradient models. Other 

analysis approaches that generate lichen response indexes are also acceptable if 

they satisfy the requirements for a successful gradient model. These requirements 

are explained in detail both to set benchmarks during model development and to 

support appropriate evaluation of a completed model. Calculation of lichen response 

indexes for new plots not part of the original model is necessary to support longterm monitoring. One should always first try to develop a single ordination model. 

Such a model has the possibility to generate all needed lichen response indexes 

from one analysis procedure. A hybrid model requiring several separate analyses 

has been successful for several regions where independent lichen air quality and 

climate response indexes could not be developed from a single ordination model.

We explain the rationale for the recommended approaches, describe in detail 

steps in the model development process for the two recommended approaches, and 

explain how to document and evaluate results, all to support successful application 
of a model for monitoring. Many decisions must be made about appropriate data, 

data treatment, and data analysis at every step in the development of a model. 

These decisions differ depending on the characteristics of both the region to which 

the model will apply and the data to be used. Each decision has an impact on the 

quality of the final model, and no one set of decisions can be recommended in 

general. Evaluation of a model is required to establish the accuracy and precision 

of the lichen response indexes derived from the model; this precision in turn sets 

the limits on how large changes must be before they can be detected by monitoring 

with the indexes. A template is included for documenting a model and archiving all 

products necessary to understand and apply it, as is required for each FIA model 

and is critical for application of a model in any other situation as well.

We have provided in this document both a general overview of the usefulness of 

this approach for monitoring forest health with lichens, and detailed instructions for 

each step in developing and applying such a model in a monitoring program. The 

detailed explanations and justifications provided at each step also support application of the principles to achieve similar goals with new approaches or analytical 

techniques and in novel situations.
Contents

1 Introduction

3 Background and History of the Lichen Indicator

5 Assessment Questions

9 Lichen Community Gradient Models

9 Rationale for Approach

11 Program Requirements

13 Lichen Model Regions 

19 Lichen Gradient Model Budgets

20 Development of a Lichen Gradient Model

20 Data Needs

23 Methods and Tools

29 Steps in Building a Lichen Gradient Model

37 Evaluating a Lichen Gradient Model During Development

38 Evaluating, Comparing, and Calibrating Completed Lichen 

 Gradient Models

38 Evaluating a Completed Lichen Gradient Model

40 Comparing and Calibrating Lichen Gradient Models

40 Summary Discussion

41 Acknowledgments

42 English Equivalents

42 References

49 Appendix: Template for a Report Section Formally Documenting 

 a Lichen Gradient Model and Its Implementation

62 Glossary
1

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Introduction

The study of lichen communities in forest ecosystems allows researchers to address 

several key questions concerning natural resources: biodiversity, response to 

climate and air pollution, ability to provide ecosystem services, and sustainability 

of timber production (McCune 2000). Lichens are strong indicators of air pollution 

and are also indicators for forest response to other environmental factors such as 

climate (Bates and Farmer 1992, Nash 2008). Lichens have little economic value, so 

their responses are not confounded by deliberate human management. Monitoring 

the quality and sustainability of ecosystems with lichens has occurred worldwide 

for many different purposes, using many different methods, with many different 

constraints, and at many different spatial scales (Nimis et al. 2002). 

The Forest Inventory and Analysis (FIA) Program of the U.S. Department of 

Agriculture Forest Service (USDA FS) monitors forests to facilitate responsible use, 

conservation, and maintenance of function for forested ecosystems. The FIA Program includes lichens in its surveys of forests on permanent plots across the Nation 

to inventory status (one-time assessment) and monitor trends (repeated assessment) 

in response of forests to environmental factors over time within quite strict budgetary constraints. It is the most extensive lichen monitoring program in the world, and 

is an excellent model for a large-scale program. As the only biological indicator in 

the FIA Program that is not a vascular plant, the lichen indicator is to some extent 

a potential surrogate for the many nonrepresented members of forest communities 

(for instance birds, insects, mycorrhizal fungi) that are also not vascular plants. 

Field data collection and analysis protocols for the FIA Program have been designed 

to address program purposes within program constraints (McCune 2000, Will-Wolf 

2010). 

Three major categories of assessment questions have been identified so far for 

lichens in the FIA Program: air quality, climate, and biodiversity. Questions about 

biodiversity can be addressed with relatively simple summaries of field data, but 

more complicated analyses are needed to evaluate response to environmental factors. Development of a multivariate lichen response model from community data is 

a powerful tool to indicate forest response to multiple environmental factors in such 

situations. Lichen response indexes derived from such a model are single numbers 

that indicate the response of lichens at a single site to particular defined environmental factors, as compared with response at other sites in the same region. 

A critical step in the implementation of the lichen indicator in FIA is the 

development of such lichen response indexes from an “FIA lichen gradient model” 

for a geographically defined “FIA lichen model region.” The two phrases in quotes 

are used throughout this document with very specific meanings (see glossary). The 
2

GENERAL TECHNICAL REPORT PNW-GTR-807

goal of such a lichen gradient model is to develop lichen indexes that represent the 

response of the lichen community to the wide range of values for each important 

regional environmental factor (a response “gradient” for each factor) of interest to 

FIA for forest resource assessment (also see Jovan 2008 for an introduction to the 

model development process). Each lichen response gradient needs to be independent of all other lichen response gradients defined for that model. This powerful 

approach to developing biotic response indexes uses all lichen community data, so 

it is not as limited by idiosyncratic species distributions as are approaches that rely 

on response of a few species. These lichen response indexes are the primary lichen 

bioindicators used in FIA for assessing regional patterns and monitoring trends 

over time.

The primary target audience for this document is researchers interested in 

developing a lichen gradient model for the FIA Program. Researchers and data 

analysts, both inside and outside the FIA Program, who wish to improve their 

understanding of the rationale for FIA lichen response indexes should find this 

document helpful as well. Other researchers who want to assess patterns and 

monitor long-term trends for any biological community will also find this approach 

useful, especially in situations where monitoring of large and diverse areas must be 

accomplished within tight budget constraints. The rationale and methods for each 

step in developing a model for the FIA Program are described in detail to illustrate 

the importance of each decision along the way for ensuring the scientific validity 

and usefulness of the final product.

Authors assume at least basic understanding by the reader of general plant 

community ecology theory and practice, in particular the use of multivariate 

techniques for data analysis and analysis of plant community response to environmental gradients (see glossary). Reference to this document is not necessary to use 

lichen response indexes for the assessment of forest ecosystem patterns and trends. 

An appendix gives a template for model documentation that includes all necessary information for calculation of lichen response indexes from field data using 

a completed lichen gradient model. A glossary includes definitions of the most 

important words and phrases used in this document. Other USDA FS publications 

focus on assessment of the quality of FIA lichen data (Patterson et al., in press) and 

guidelines for standard analysis of FIA lichen data (Will-Wolf 2010).

Development of a 

lichen gradient model 

is a critical step in 

implementation.
3

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Background and History of the Lichen Indicator

A lichen is a close association (symbiosis) of a fungus with photosynthesizing algae 

or cyanobacteria; this association is named for the fungus and functions as a discrete 

individual (an honorary “plant”). A lichen community is an assemblage of these “species” of lichens (Nash 2008). The sensitivity of lichens to air pollution, the response 

of lichens and lichen communities to climate, and the use of lichens as indicators are 

well-documented in scientific literature (McCune 2000, Muir and McCune 1988, 

Nimis et al. 2002, Richardson 1988). The protocols adopted for lichen sampling and 

data management in the FIA Program are designed and standardized to enhance 

ability to accurately assess patterns and monitor trends regionwide and at larger geographic scales. Data are less suitable to assess precisely the response of lichens to local 

and within-plot factors such as stand age, tree species composition, and recent disturbance history (Will-Wolf 2010). The precision of FIA lichen indicator data is constrained by requirements for protocols to be both affordable and consistently applied 

for monitoring thousands of plots across the entire Nation into the indefinite future 

(McCune 2000; Smith et al. 1993; Will-Wolf 1988; Will-Wolf et al. 2002b, 2004). 

For the purposes of the FIA Program, the lichen sample population is restricted 

to macrolichens (leafy, tufted, or hanging lichens) found on living or standing dead 

woody substrates, including both trunks and branches of trees and shrubs. Only 

standing and recently fallen woody substrates are included, thereby standardizing the 

measurements to a class of substrates that can be found at all forested sites. Lichens 

are collected and assigned abundance scores by a trained (Will-Wolf 2007, Will-Wolf 

and Neitlich 2007) nonspecialist crew person in a time-constrained (up to 2 hours) 

search of an FIA plot following standard protocols (USDA FS 2004, 2010b), then are 

identified by a lichen specialist also following standard protocols (Will-Wolf 2009). 

These plot data collection protocols and low plot density together set limits on the 

precision achievable with lichen indicator data. 

The lichen community indicator is implemented in the FIA Program in two phases 

(fig. 1). 

1. In the calibration phase, a lichen gradient model is developed for lichen 

communities of a particular region and is calibrated to isolate and describe 

air quality, climate, and other gradients important for resource assessment. 

Lichen indexes of response to the described gradients are the products.

2. In the application phase, the model is used to calculate, for newly sampled 

plots, values for the lichen response indexes developed in the calibration 

phase. These indexes are the primary tools used to answer resource assessment questions concerning the spatial patterns and trends in the condition of 

our forest resources, as described above. 
4

GENERAL TECHNICAL REPORT PNW-GTR-807

The first lichen gradient model was published in 1997 (McCune et al. 1997a, 

1997b). Thus far, four more models have been completed and four are in progress. 

Only after such a model has been completed and implemented can a region move 

into the application phase (fig. 1).

Figure 1—The two phases of implementation of the Forest Inventory and Analysis (FIA) lichen indicator–During the calibration phase, standard FIA plot data and other data are collected and a lichen 

gradient model is developed. During the application phase, lichen response indexes are calculated 

for new plots from the lichen gradient model based on lichen data collected for each plot. These plot 

indexes are used to analyze patterns across space and trends over time related to the defined lichen 

response gradients. Other abbreviations: CVS = U.S. Department of Agriculture Forest Service 

national forest current vegetation survey, FHM = Forest Health Monitoring Program, QA = quality 

assessment and quality assurance. Adapted from McCune et al. (1994). 

LICHEN DATA

• Standard plots (FHM and FIA)

• Supplemental urban/industrial

 and other reference plots

• CVS plots, special studies

• QA plots (hot audits and blind checks)

ENVIRONMENTAL DATA

• Mensuration (tree basal area, 

 percentage conifers, other)

• Climate (precipitation, temperature, other)

• Lichen tissue content (percentage

 sulfur or nitrogen, other)

• Physical (elevation, location, other) 

Interpretation and 

isolation of lichen 

response gradients 

Final lichen gradient model 

Final lichen gradient model 

Resource assessment: spatial pattern and trend

CALIBRATION PHASE

Extraction of lichen 

community gradients 

from analysis of plotby-lichen species table, 

other techniques

APPLICATION PHASE

Calculate response index values for

lichen air quality and climate gradients

for each plot using calibration data set

and designated algorithms 

Repeat using revisited

standard plots
5

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Several lichen parameters can be calculated from FIA lichen data while the 

geographic area is still in the calibration phase (Will-Wolf 2010). The lichen species richness index is the total number of species recorded on an FIA plot. The 

total number of species in the data set for a specific geographic region represents 

regional diversity. The rate of species replacement, or turnover, across a specific 

geographic region can be calculated as regional diversity divided by average plot 

species richness index. All of these parameters are based on the count of lichen 

species at an FIA plot. 

The lichen indicator was developed in 1990–1992 for the Environmental Monitoring and Assessment Program, was implemented in the USDA FS Forest Health 

Monitoring Program 1993–1999, and has been part of the FIA Program since 2000. 

Field and laboratory protocols for the lichen indicator have not changed; data from 

all years are fully compatible. The basic guidelines for analysis of lichen indicator 

data also remain unchanged since the beginning. Data from the first round of plot 

sampling are used to assess the initial condition and patterns of lichen communities 

in a region. Data from repeat samples of the same plots are used over time to monitor change in condition and in pattern of response to air quality, climate, and other 

environmental factors.

Assessment Questions

The general rationale behind key FIA lichen indicator assessment questions is 

addressed in more detail elsewhere (Geiser and Neitlich 2007; Jovan 2008; McCune 

2000; McCune et al. 1997b; Will-Wolf 2010). In this document we discuss the major 

identified assessment questions (air quality, climate, and biodiversity) specifically 

as they relate to the development of lichen gradient models. Originally (McCune 

2000; McCune et al. 1994, 1997b) the FIA lichen indicator’s primary usefulness to 

the program was seen as estimating response to sulfur (S)- and nitrogen (N)-based 

air pollution. The purpose of a lichen gradient model was to develop indexes of 

response to air pollution statistically independent of response to climate. More 

recently, it has become apparent that lichen climate and biodiversity indexes are as 

important as lichen air quality indexes for assessing and monitoring forest health 

in the FIA Program (Jovan 2008, Will-Wolf et al. 2002a). We suggest guidelines 

for relating additional categories of potentially interesting assessment questions to 

completed lichen gradient models. 

Lichen gradient models 

address assessment 

questions.
6

GENERAL TECHNICAL REPORT PNW-GTR-807

Air quality—

Does regional air quality affect our forests? Have the effects changed? If so, 

is air quality improving or deteriorating? In what areas is it changing?

Addressing these air quality assessment questions in the FIA Program with lichen 

data reflects lichen response primarily to S and N compounds in air pollution. 

A lichen air quality index requires the development of a regional lichen gradient 

model to quantify lichen response. Historically research has focused on lichen 

response to airborne acidic N and S compounds (often abbreviated in literature 

as NOx

-

 and SOx

-

) and heavy metals from industry, traffic, and urban centers 

(Hawksworth and Rose 1976, Richardson 1998, Smith et al. 1993, van Dobben 

1993). More recently, lichen communities have been shown to respond in different 

ways to ammonia (NH4

+

) and other usually alkaline N compounds that come 

primarily from agriculture and animal husbandry, and also from incomplete 

combustion of gasoline engines (Fenn et al. 2003a; Jovan and McCune 2005, 

2006; Sillett and Neitlich 1996; van Dobben and ter Braak 1998). Related research 

(summarized in Jovan 2008) has identified two possible air pollution effects from 

alkaline N air pollution: increased pH that affects both lichens and substrate 

chemistry (only indirectly related to the chemical elements involved), and direct 

fertilizing effects of N compounds. Thus air pollution effects on lichens can have 

multiple sources with effects sometimes difficult to distinguish (Jovan 2008). This 

has made interpretation of lichen response to air pollution more complicated. The 

interpretation by McCune et al. (1997b) that their lichen air quality index for the 

Southeast Lichen Gradient Model represents response primarily to acidic urban/

industrial pollution is supported by regional air quality data (NADP 2009) and 

more recent studies in a partially overlapping region (Will-Wolf et al. 2006). Geiser 

and Neitlich (2007) concluded that although all three potential types of pollution 

were present, acidifying pollutants and fertilizing N pollutants (Jovan 2009) were 

probably the most important factors affecting their lichen air quality index. Jovan 

and McCune (2005, 2006) concluded that in each of two regions, their lichen air 

pollution response index was most strongly linked to the combined effects of 

alkaline and fertilizing N compounds. For each FIA lichen gradient model, in the 

future, it will be important to consider all of these potential pollution sources and 

their effects on lichen response, and distinguish them to the extent possible as 

research progresses.
7

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Climate—

Are changes in climate affecting our forests? If so, how do we characterize 

the effects? In what areas is it changing?

Changing climate is an ever-more-likely prospect whose impact on lichens and 

forests needs to be assessed (Bates and Farmer 1992, Ellis et al. 2007, Søchting 

2004). The process for development of an FIA lichen gradient model was designed 

from the beginning to produce independent indexes for lichen response to climate 

and to air pollution (McCune et al. 1994, 1997b). Current climate response indexes 

for different lichen model regions differ in number and in the climate factors to 

which they are related; this is expected to be the case for future models as well. 

The McCune et al. (1997b) climate response index was correlated with temperature and elevation; the Geiser and Neitlich (2007) climate response index was 

correlated with temperature, moisture, and elevation. Jovan and McCune (2004) 

developed two independent indexes of lichen response to climate: one correlated 

with temperature and elevation, the other correlated with moisture. Will-Wolf et al. 

(in press) found that the FIA lichen species richness index was variously correlated 

with different climate factors in different regions of the coterminous United States, 

illustrating the variety of responses to climate exhibited by lichen communities of 

different areas. 

Biodiversity—

Is the lichen component of biodiversity changing through time? What 

environmental drivers seem to be most strongly linked to patterns and 

changes in lichen biodiversity? 

Lichens contribute a substantial proportion of the visible species of many forests 

and they participate in ecosystem processes such as nutrient cycling and species 

interactions in ecosystems (Nash 2008, Sharnoff and Rosentreter 2008). The FIA 

lichen species richness index tracks one aspect of lichen diversity patterns, is 

currently the only defined lichen biodiversity index, and does not require a lichen 

gradient model for its calculation. Species composition and abundance in lichen 

communities are not used directly to assess lichen diversity patterns; they are 

included in lichen air quality and climate indexes that require a lichen gradient 

model for their calculation. The model-derived lichen response indexes thus indirectly provide insights into potential causes for variation in lichen biodiversity and 

changes over time. Additional components of lichen biodiversity may emerge as 

products of the development of a lichen gradient model. 
8

GENERAL TECHNICAL REPORT PNW-GTR-807

Other categories of assessment questions—

Other categories of assessment questions applicable to the lichen indicator are 

expected to arise in the future. For example one category of questions currently 

being investigated concerns the effect of nearby land cover on the FIA lichen 

species richness index across the conterminous United States. Preliminary findings 

suggest that for some large geographic areas (equal to or larger than FIA lichen 

model regions), nearby land use is correlated with lichen species richness index 

independent of air quality and climate factors, an interesting and unexpected finding (Will-Wolf, unpublished data). Such a finding, if confirmed, would suggest the 

need to consider landscape pattern factors (for example, the percentage of nearby 

land in forest) along with air quality and climate factors in the development of a 

lichen gradient model for some regions. Some correlation between a lichen air quality index and a tree health index for part of the Eastern Broadleaf Forest ecoregion 

has been found in one pilot evaluation (Will-Wolf and Jovan 2008). Possible linkage of a lichen land use response index with a tree health index would be of equal 

interest to FIA for forest health assessment. Linkage between a lichen response 

index and landscape pattern at large geographic scales is more likely to be found in 

the Eastern United States, where climate variation is relatively low and much of the 

landscape has been strongly modified by human activity. 

For another example, Geiser and Neitlich (2007) identified one independent 

gradient in lichen species composition correlated with the relative abundance of 

N-fixing lichens in the community. This lichen response is probably related to pollution in some way, but it was uncorrelated with any of the authors’ environmental 

or forest structure variables, so it remains undefined. A lichen biodiversity response 

index could be calculated for that model, but its usefulness for biomonitoring will 

be limited until it can be linked to one or more potential causal factors. 

Most ecological studies are able to identify no more than three to four mutually independent environmental factors correlated with variation in community 

composition. This constraint applies to the lichen gradient modeling approach as 

well. Modeling approaches such as partial correlation or multiple regression can be 

used to explain the relative correlation of several additional environmental factors with a particular lichen response index at the scale of an entire lichen model 

region. Supplemental intensive studies designed to isolate the effects of these other 

environmental factors would be needed to assess the relative contribution of such 

additional variables to explaining lichen response indexes.
9

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Lichen Community Gradient Models

The goal of each lichen gradient model is to provide indexes for lichen community 

response to air quality, to climate, and possibly to other environmental factors that 

vary at the scale of large regions. We use the term “model” in this document to 

mean a template for defining and calculating indexes of response to environmental 

factors of interest from the lichen species composition of a field plot. It is useful for 

each index developed from such a model to be as insensitive as possible to variation 

in local site factors such as stand structure and disturbance history, to avoid interfering with the desired goal of tracking large-scale patterns.

Rationale for Approach 

The gradient model approach is particularly powerful for large-scale assessment of 

forest health with lichens and is required in the FIA Program for definitive assignment to plots of indexes representing responses by lichens to important regional 

environmental factors. The lichen species richness index (that does not require a 

model) usually relates to many factors that affect lichen communities, so interpreting how a change in this index relates to particular assessment questions is difficult 

without a gradient model (Will-Wolf 2010, Will-Wolf et al. 2006). For example, 

McCune et al. (1997b) found that the lichen species richness index for the Southeast 

Lichen Model Region was equally correlated with air quality and climate. Derived 

indexes representing response of lichen community composition that are correlated 

with individual environmental factors, on the other hand, are much more reliable 

indicators of the influences of those factors. One requirement for indexes of lichen 

response defined from a gradient model for the FIA Program is that they be independent of one another, so interpretation of response to one factor is not confounded 

by response to a different factor. Variations in these indexes across a region reflect 

the initial state of the resource with respect to assessment questions, and changes in 

these indexes over time reflect trends related to assessment questions. 

The modeling approaches recommended here involve use of several common 

and well-tested analysis techniques applied to lichen species composition at a 

plot. Indexes based on most to all species at a plot (rather than a few preselected 

individual species) are preferred to ensure they can be calculated for a majority 

of FIA plots. Given the nature of FIA lichen data (Will-Wolf 2010) and variation 

across FIA lichen model regions, even common lichen species are usually found 

at no more than 50 to 70 percent of FIA plots in a lichen model region (Geiser and 

Neitlich 2007; Jovan and McCune 2005, 2006; McCune et al. 1997b; Will-Wolf 

et al. 2006). 
10

GENERAL TECHNICAL REPORT PNW-GTR-807

A multivariate lichen community analysis method should be the first modelbuilding approach tried for two reasons: (1) it is possible to develop all necessary 

independent gradients and lichen response indexes from a single analysis, and (2) 

this modeling approach uses information for all the lichen species found at a plot. 

Unconstrained ordination is, for this purpose and for FIA Program needs, the most 

appropriate of several multivariate analysis approaches designed or adapted and 

currently in common use by ecologists to extract important information from large 

species-by-plot data tables (Legendre and Legendre 1998, McCune and Grace 2002 

[both these references use the terms “ordination,” “simple ordination,” or “indirect 

gradient analysis” rather than “unconstrained ordination,” contrasting them with 

“canonical” analysis or “constrained ordination”]; Økland 1996). In unconstrained 

ordination, the investigator identifies major gradients (ranges) of variation in 

macrolichen community composition, then compares these community gradients a 

posteriori with individual environmental factors. Unconstrained ordination is preferred for exploratory analysis and the development of a model. It avoids distortion 

of community gradients from a priori selection of environmental factors (required 

for constrained ordination), avoids bias from inclusion of correlated environmental 

factors, and allows for valid statistical evaluation of correlations between community gradients and environmental factors (McCune 1997b, Økland 1996). These 

characteristics are critical to unbiased development of lichen response indexes for 

factors of interest, to ensure indexes are scaled appropriately to response of lichens 

rather than to convenient measurement scales for environmental factors, and to 

reliably assign to new plots response index values for the developed gradients. 

Assigning lichen response indexes to plots not in the original model data set is 

a critical function to support monitoring, although use of an existing ordination 

model to do this is not a typical application for multivariate community analysis in 

plant ecology. 

Nonmetric multidimensional scaling (NMS, or NMDS in some recent literature) ordination (Kruskal 1964) is an unconstrained ordination technique that 

has been preferred by FIA lichen gradient model developers historically because 

it generates minimal distortion of original patterns compared to other popular 

unconstrained ordination techniques (Minchin 1987). Other classes of multivariate analysis methods including classification, constrained ordination (including 

Canonical Correspondence Analysis and Redundancy Analysis), and direct gradient 

analysis (evaluating variation in species composition along preselected ranges of 

environmental factors), fail on one or more of the above criteria. 

A hybrid modeling approach incorporating different analysis methods for 

response to different environmental factors may be needed when unconstrained 

Multivariate analysis 

is recommended 

to develop a lichen 

gradient model.
11

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

ordination alone does not generate a satisfactory model. For some regions, an 

ordination model has not successfully extracted a lichen response index for air 

quality independent of a stronger climate response index. In this case, an alternate 

approach using indicator species and regression modeling has proved successful to define an independent air pollution response index. Objective, quantitative 

identification of sets of at least 6 to 10 indicator lichen species (indicators of poor 

air quality or of some other environmental factor, depending on what is needed) is 

the first step. The more indicator species that are identified the better, to enhance 

the liklihood that the final index reflects response to the factor much more than just 

chance variation in species distributions. Indicator Species Analysis (Dufrêne and 

Legendre 1997) has thus far been the most successful technique used to quantitatively identify indicators of polluted areas. This can be followed by construction of 

a regression model using the proportion of indicator species at a plot and the residuals from regression against the primary lichen response factor, to model response to 

the secondary factor independent from the primary factor. 

Program Requirements

The two subsections below present the standards that must be met and the products 

that must be provided for a lichen gradient model useful in the FIA Program. 

A model developed using the approaches and techniques described later in this 

document can meet these standards and support provision of these products. New 

analysis approaches and techniques that meet these criteria would also be suitable 

for the FIA Program. The current FIA lichen indicator advisor has the responsibility to determine whether a new lichen gradient model meets all required criteria 

and is acceptable for use in the program. Many of these criteria would be appropriate for other monitoring programs as well.

Required performance standards—

To meet the goals for the FIA Program, a lichen gradient model must meet the 

following performance standards:

1. A model must define at least one lichen air quality response index and at 

least one lichen climate response index independent of each other and of 

other major environmental factors.

2. There must be an objective and unbiased method for estimating what proportion of variability in lichen community composition is explained by 

each lichen response index, and for estimating the strength and statistical 

reliability of the correlation of that index with the relevant environmental 

factor(s). Correlation of the relevant environmental factor(s) with that index 

must be strong enough to support the interpretation that the lichen index 
12

GENERAL TECHNICAL REPORT PNW-GTR-807

represents lichen community response to that environmental factor. The 

set of indexes (minimum two, No. 1 above) developed from a single NMS 

lichen gradient model (or other acceptable single model) should account 

for at least 50 percent of the information in the full model lichen data set, 

as calculated using standard approaches. For a hybrid model, a reasonable 

standard is that each index from a hybrid model should account for at least 

25 percent of total information in the data set used, as calculated using 

standard approaches.

3. The numerical values of the derived indexes should be scaled appropriately 

to the response of lichens (as opposed to convenient measurement scales 

for environmental factors). For instance, the same magnitude of differences 

between two index values at the high, middle, and low end of the numerical 

range of the index should represent approximately equivalent differences 

in lichen species composition. This simplifies analysis of trends over time. 

Such equivalent scaling of lichen responses happens automatically with 

unconstrained ordination; this must be evaluated for other model development techniques.

4. There must be an objective, repeatable, and reliable method to assign lichen 

response indexes to new plots, and to repeat samples of plots, from an existing model. 

5. The variability of response indexes from repeated samples of the same plot 

in the same year must be calculated for each defined lichen response index 

to allow statistical evaluation of patterns and trends. This repeat sample 

variability for all lichen response indexes (applicable to FIA lichen data 

meeting standard quality criteria) must be small enough to support realistic 

trend analysis (see section “Evaluating a Lichen Gradient Model During 

Development,” p. 37, for description of standards).

6. Calculation of all lichen indexes and estimates of variability must use data 

on lichen species composition at a plot. Only plot environmental factor(s) 

routinely archived by FIA (such as plot elevation or location) should be 

required for the calculation of a lichen response index, and then only when 

a lichen response index based solely on lichen species composition has 

proved unsatisfactory. No nonstandard field sampling by FIA crews should 

be required to calculate a lichen response index for FIA plots.

Required products—

Standardized documentation for the model is a required component of the final 

report to ensure that each lichen gradient model is usable in the FIA Program. We 

provide in an appendix a template for this standardized documentation. Copies of 

Models for FIA must 

meet these standards, 

provide these 

products, and cover 

regions that meet 

these criteria.
13

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

instructions and data files necessary for assignment of lichen indexes to sample 

plots are to be permanently archived and made available to government officials 

and the general public to calculate lichen response indexes from the model. Key 

required products include:

1. Clear description and rationale for the geographic region for which the 

model applies.

2. Clear description of procedures used to develop the model.

3. Clear description of all lichen response indexes developed from the model.

4. Clear description of accuracy, precision, and bias of all lichen indexes.

5. Detailed instructions for assigning lichen response indexes to new plots 

based on lichen sample data.

6. Original data sets, all other data sets, formulas, and information required 

for assignment of lichen response indexes are to be archived.

Requirements 3 and 5 are met with the formal documentation of a model using 

the template presented in the appendix. Requirements 1, 2, and 4 are usually met 

in the text of a lichen gradient model final report and peer-reviewed publication. 

Requirement 6 is met through archiving files with FIA regions and central information management offices, or another program that uses the model. Summaries of or 

information about requirements 2 and 6 are also included in the formal documentation of a model.

Lichen Model Regions 

Because composition of lichen communities varies widely across the United States, 

separate gradient models are needed for different biogeographic regions. In this section, we describe proposed lichen model regions for the FIA Program and explain 

our rationale for the proposed boundaries. A reader interested in other applications 

for our modeling approach should review our explanation of boundary delineation 

(second subsection below). Selection of an appropriate geographic region for which 

the model will apply is a critical first step in developing the model that can affect 

the model’s long-term success.

Proposed lichen model regions and funding priority—

We propose 17 potential FIA lichen model regions in the conterminous United 

States, with four additional regions for Alaska and Hawaii (fig. 2). Boundaries of 

regions where models are not currently being developed are tentative; they are 

expected to change somewhat as more information is accumulated. Table 1 presents 

our recommendations for funding priorities for these proposed lichen gradient 

models. These recommendations are based on current availability of FIA data 
14

GENERAL TECHNICAL REPORT PNW-GTR-807

Figure 2—Boundaries of existing and proposed Forest Inventory and Analysis Program lichen model regions for the United States. 

Twenty-one regions (including Alaska and Hawaii) have been outlined as of 2010. Lichen model regions with completed models (see 

table 1) have had boundaries tested and found to be ecologically supported. Those with models in progress or not yet begun have 

proposed boundaries that may be revised as a lichen gradient model is completed. Many boundaries follow Bailey et al. (1994) and 

Cleland et al. (2005).

FIA lichen model region

Pacific Northwest West

Pacific Northwest East

Northern Rockies

Northwestern California

Sierra–Southern Cascades–Modoc

Great Basin–Semidesert

Greater Central Valley

Southern California Coast

Southwest

Middle and Southern Rockies

North Central

Northeast

Ohio Valley

Mid Atlantic

South Central

Southeast

Southern Florida

Boundaries

States

Ecoregion sections

Alaska Arctic and

Interior

Southeast

and South

Central

Southwest

and Aleutian

Hawaii

within the region, pollution severity (based on National Air Deposition Program 

ionic deposition maps, NADP 2009), and known regional forest health issues. 

Because lichen communities may well change over a 10-year period, data sets 

for development of a lichen gradient model should include data collected over no 

more than 5 to 8 years. Priorities reflect the need to develop models for some areas 

soon so older data do not become unusable. Indexes for forested plots not currently 

located in a lichen model region will be calculated from the model for the nearest 

adjacent lichen model region, then tested for their usefulness.

As of spring 2010, the Southeast Lichen Model Region (McCune et al. 1997b), 

Colorado (McCune et al. 1998), most of California (Jovan and McCune 2004, 2005, 

2006), and Washington and Oregon west of the Cascade crest (Geiser and Neitlich 
15

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Table 1–Forest Inventory and Analysis Program (FIA) lichen model regions and progress toward completion 

of models 

Lichen Funding

model region Status priority Regional issues Publications

New England In progress Funded High urban/industrial pollution in 

 southern and coastal areas. Strong 

 latitudinal gradient in regional pollution.

Mid Atlantic In progress Funded Urban/industrial pollution high in local 

 areas; regional pollution affects coastal 

 and northern sections.

Southeast Completed Funded Moderate localized pollution. McCune et al. 

 1997a, 1997b

North Central In future Medium Moderate/high pollution in some areas. 

 1994–1999 plot data available. Gradient 

 model will require recent data.

Ohio Valley In future Low Pollution high. Plot data sparse as of 2009.

South Central In future Low Moderate pollution. No plot data as of 2009.

Northern Rockies In progress Funded Moderate localized pollution. Plot data 

 available.

Middle and (Colorado (Funded for Moderate localized pollution. Test Colorado Colorado: McCune

Southern completed) Colorado) model for applicability to the broader et al. 1998

Rockies High region. If funded soon, an expanded model 

(Colorado) can include available older plot data.

Southwest In future High Moderate pollution. Plot data available.

Pacific Northwest Completed Funded Urban/industrial pollution moderate/high Geiser and Neitlich 

West in some areas; high lichen diversity in 2007

 most areas.

Pacific Northwest In progress Funded Moderate local pollution; high lichen

East/Northern diversity in most areas.

Great Basin

Greater Central Completed Funded Agricultural and urban/industrial pollution Jovan and McCune

Valley (of severe. 2004, 2005

California)

California Sierras Completed Funded Agricultural and urban/industrial pollution Jovan and McCune

 moderate/high in places. 2004, 2006

Central Great Basin In future Medium Moderate localized pollution. Plot data 

 available; construct gradient model soon 

 to link with older data.

Alaska—Southeast In progress Funded Coastal point sources of pollution; most 

and South Central of region high lichen diversity. Plot 

 data available.

Alaska—Southwest In future Low Pollution extremely low. No plot data

and Aleutian available as of 2009.

Alaska—Interior In future Low Arctic haze, Eurasian pollution moderate.

and Arctic No plot data available as of 2009.

Hawaii In future Low Pollution moderate. No plot data available 

 as of 2009.

Notes: Twenty-one lichen model regions have been identified as of 2009. Three regions on figure 3–South Florida, Northwestern California, and 

Southern California Coast–are not listed in this table. Their small size probably precludes use of standard FIA lichen gradient model development 

techniques; indexes of lichen response to environmental drivers of interest may be developed in other ways. Also see table 2. 
16

GENERAL TECHNICAL REPORT PNW-GTR-807

2007) have completed lichen gradient models in place and are in the application 

phase (fig. 1). Lichen gradient model development is in progress (table 1) in several 

other regions. 

The number of proposed lichen model regions is a compromise between competing constraints on development of useful models. If a lichen model region is too 

large, representation of regional patterns with lichen response indexes becomes less 

reliable; repeat sample variability of indexes may be too high to support realistic 

trend analysis. Conversely, delineation of more and smaller model regions would 

inflate the cost of completing all models beyond the capacity of FIA Program funding and might not adequately represent large-scale patterns. The average lifetime 

of a lichen gradient model also affects cost to the program of the modeling effort. 

Climate change might reduce the effectiveness of trend detection from a more 

precise model for a smaller region faster than for a more general model for a larger 

region, and thus shorten the lifetime of a model for a too-small region. To balance 

all these factors, we have outlined the largest lichen model regions expected to 

support acceptable models. 

Five currently delineated FIA lichen model regions might be studied for possible consolidation into other adjacent regions (table 2): three because they might 

be too small for model development using standard approaches, and two because 

gradient models for adjacent regions might well be successfully applied to calculate 

lichen response indexes for plots. The extent to which the existing model for the 

state of Colorado may apply to the broader Middle and Southern Rockies region 

(fig. 1, table 1) should be explored, as should the application of the recently funded 

model for the Northern Rockies region to some northern Colorado plots in areas 

ecologically similar to Wyoming and Montana mountains. Virginia plots, shown 

in figure 2 as part of the Mid Atlantic Lichen Model Region, are also treated as 

part of the Southeast Region. This assignment to two model regions will facilitate 

evaluation of the boundary between the two regions in the future. The feasibility of 

such region consolidations and boundary adjustments should be investigated after 

Table 2—Separate lichen model regions (see fig. 2, table 1) worthy of study to group with 

other regions

Lichen model region Alternate grouping

Southern Florida Small; possibly combine with Southeast.

Ohio Valley Possibly combine with adjacent regions, especially South Central.

Northwestern California Small; possibly combine with Pacific Northwest West.

Southern California Coast Small; possibly combine with Southwest.

Alaska—Southwest and Aleutian Possibly combine with Alaska—Arctic and Interior.
17

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

models for adjacent areas have been established and more data are available. The 

FIA data archiving protocols have been designed to accommodate plots assigned 

lichen response indexes from multiple gradient models; this feature will facilitate 

future reevaluation of model region boundaries without massive reanalysis.

Model region boundary delineation—

We have based our designation of lichen model regions here on a combination of 

ecoregion boundaries (Bailey et al. 1994, Cleland et al. 2005, Omernik 1987, US 

EPA 2008), annual precipitation (Daly and Taylor 2000) for western regions, political boundaries, and estimated lichen species composition change across regions. 

Considering both usefulness and economical concerns for lichen gradient models 

(see above), we combined as many ecoregion provinces (or sections) as possible 

using these two moderately subjective but reasonable criteria:

1. Lichen community variation cannot be too high within a lichen model 

region. The collected wisdom of past and current gradient model developers is that if estimated change across a region in lichen species composition 

is too high, myriad problems arise in multivariate community analysis. For 

the FIA Program an estimate for such change is calculated as total lichen 

species number in the region divided by average number of species per 

plot (Will-Wolf 2010); if this value is more than about 10, problems arise 

(see McCune and Grace 2002). This species change estimate is related to 

the terms “beta diversity” and “species turnover rate” used in ecological 

literature. In practice, very high estimates of species change across regions 

are typically attained by combining adjacent areas that have very different floras. For instance, a combined data set for the Pacific Northwest West 

and East model regions has 284 lichen species at 846 plots, an average 

of 17.9 species per plot, and a species change estimate of 284/17.9 or 15.9 

(Neitlich, unpublished data). A combined data set for the New England and 

Mid Atlantic model regions has 184 lichen species at 271 plots, an average of 11.5 species per plot, and a species change estimate of 184/11.5 or 

16.0 (Will-Wolf, unpublished data). Both these species change estimates 

are much higher than desired, so the original lichen model regions should 

stay separate. Species change estimates are less than 10 for each of the four 

original model regions. Species change estimates were not examined on 

a quantitative basis countrywide because FIA lichen data are either very 

sparse or are absent for parts of several proposed model regions. We made 

judgment-based estimates for regions as needed. 
18

GENERAL TECHNICAL REPORT PNW-GTR-807

2. Lichen model regions must be geographically compact. A requirement 

for geographically compact lichen model regions is included because geographic cohesion makes it reasonable to expect that lichens can disperse/

migrate throughout the region. Thus a possible mechanism for change 

(lichen migration in response to changing climate or air pollution, for 

instance) can be linked with observation of changes in the related lichen 

response indexes over time. In the West, boundaries of compact regions 

were qualitatively aligned with precipitation range from a high-resolution 

annual precipitation data model (Daly and Taylor 2000). In the East, with 

more gradual change in climate across regions, boundaries were loosely 

related to the correspondence of ecoregion section boundaries with state 

boundaries.

The delineation process was therefore a guided, subjective, and iterative process 

based on the expert opinion of the authors and other gradient model developers, 

with the expectation that boundaries will be modified as gradient models are developed and tested. Figure 2 results from at least the third iteration of the delineation 

process following the two above criteria. Forested FIA plots are found in the white 

areas of figure 2, but plot density in those areas is estimated to be too low (probably 

fewer than one tenth of FIA forest health plots in the area are forested) to support 

the development of separate lichen gradient models. After these areas are sampled 

for lichens, indexes will be calculated for forested plots using the nearest adjacent 

lichen gradient model. If the resulting lichen response indexes are not adequate for 

pattern estimation and trend analysis there, alternate analytical approaches will be 

considered.

A quantitative, completely objective delineation process could have been 

attempted, based on current lichen data and environmental data, but was rejected at 

this time for several reasons:

1. There are currently several large geographic gaps in lichen data (no data for 

states in the Great Plains and South Central United States, Florida, Hawaii, 

and most of Alaska) and several areas of very sparse data (e.g., New 

Mexico, North Central States). 

2. Defining thresholds of acceptable variation within and between regions 

based on lichen community variation, precipitation, annual temperature, 

and other environmental factors would still be ultimately subjective (like 

any attempt to divide a continuous variable into discrete groups).

3. In the Eastern United States, temperature appears to correlate more 

strongly with lichen community variation than does precipitation, whereas 
19

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

in the West the reverse often happens (Geiser and Neitlich 2007; Jovan and 

McCune 2004; McCune et al. 1997b; Will-Wolf et al. 2006, in press). Such 

differences make a uniform national choice of climate-based factors for 

defining lichen model regions difficult.

4. The large amount of time and funding support that would be needed to 

quantitatively delineate lichen model region boundaries before developing the models was judged to be better spent completing individual lichen 

gradient models. Once models for several adjacent lichen model regions 

are completed, the subject of boundary delineation between these adjacent 

regions should be revisited.

Lichen Gradient Model Budgets

The FIA Program lichen gradient models are estimated to cost a minimum of 

$40,000 to $45,000 each (2002 dollars). Models completed with substantial 

involvement of FIA lichen indicator advisors and/or other FIA staff whose salary originates from other sources (i.e., other funding agreements or staff salary) 

might require less model-specific funding, whereas model development completely 

outsourced to contractors is likely to cost more than a “typical” estimate (table 3). 

A project to develop lichen response indexes from a lichen gradient model independent of an established data collection program would need to include much more 

funding for data collection and processing than estimated here (table 3). 

Table 3—Sample budget for a basic Forest Inventory and Analysis 

Program lichen gradient model contracted to a university researcher

Item Amount

2002 dollars

Field work (30 lichen specialist crew-days) 9,000 

Field work (30 GS-07 biotechnician crew-days) 5,700 

Sample identification from 60 supplemental plots 3,600 

Travel (6 people + expenses for 60 crew days) 9,240 

Data analysis 8,960 

Reporting 7,000 

 Subtotal direct costs: 43,500 

Institutional indirect costs (0 to >50%)a

 Total (includes indirect costs): At least 43,500

Notes: This sample budget assumes lichen data for many (100 or more) plots are available from 

an existing field inventory program to be used in the development of a model. If they are not, 

substantial additional costs for field data collection and sample processing must be planned.

a

Depends on kind of institution, negotiated indirect cost rate 
20

GENERAL TECHNICAL REPORT PNW-GTR-807

The Pacific Northwest West (completed) and Pacific Northwest East (in progress) models cost at least $60,000 per model, involving substantial supplemental 

funding and in-kind contribution through FIA field staff participation from the 

Pacific Northwest FIA region plus an equal partnership with National Forest System Region 6 via staff time and plot data. This is an unusual situation; other efforts 

have had or are expected to have funding and resources more similar to that shown 

in table 3. Two California lichen gradient models were completed via a contracted 

supervised graduate thesis at a cost of $75,000 for one complete model and field 

work for the second model, with about $16,000 additional funding from the Pacific 

Northwest FIA region to complete the second model (approximately $48,000 per 

completed model in 2003-2006). Funding gradient models gradually, no more than 

one or two per year, allows for greater oversight of model development and thus for 

higher quality models.

Lichen tissue element analysis is often included to facilitate better interpretation 

of air quality gradients; approximately $50 per plot (2002 dollars) must be added 

to typical cost estimates (table 3) to cover field, laboratory, and data management 

costs for this procedure. A report is planned to present the comparative value and 

costs of a variety of element analysis options for FIA lichen gradient models and 

other FIA research applications. 

Development of a Lichen Gradient Model

In this section we describe the data needed for a lichen gradient model in the FIA 

Program, compare and recommend particular methods and tools for modeling, 

outline steps in developing the model, and propose methods to evaluate a model. 

All FIA lichen gradient models completed thus far have been developed using NMS 

ordination alone or in combination with Indicator Species Analysis (Dufrêne and 

Legendre 1997) followed by regression analysis. 

Data Needs

An FIA lichen gradient model developer uses a combination of standard FIA lichen 

and other plot data, additional environmental data available from other public 

sources, and supplemental field data collected specifically for the model development project. Those planning independent projects can use these guidelines to 

organize their data collection efforts.

1. Data collected using standard FIA training and quality assessment protocols 

(Will-Wolf 2007), field protocols (USDA FS 2004, 2010a), and laboratory data 

management protocols (Will-Wolf 2009) for abundance of lichen species at 
21

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

standard FIA plots distributed across an entire lichen model region for at least 

one year, preferably from multiple years. Data for at least 100 plots are needed 

for a single model region.

2. Data for abundance of lichen species at supplemental plots collected using 

standard FIA field and laboratory protocols in areas with known/measured 

poor air quality (urban/industrial areas and perhaps agricultural areas), spread 

across major climate and topographic gradients in the project region. A minimum of 20 plots is needed; double that number would be desirable.

3. Data from supplemental plots as needed from areas with documented very 

good air quality, spread across major climate and topographic gradients in 

the project region. Many standard FIA plots are in areas with relatively good 

air quality; a minimum of 20 plots total from known very good air quality 

areas must be present in the data set for model development. It is desirable to 

include in the model data set plots from “clean air” areas with field data collected by lichen specialists.

4. A minimum of 20 plots spread across all other major gradients of interest in 

the lichen model region. Some supplemental plots may be required to meet 

this objective. For instance, if regional gradients of young to old-growth forests or hardwoods (deciduous trees) to softwoods (conifer trees) are thought to 

be important, data from supplemental plots near endpoints of these gradients 

may need to be collected. 

The final data set of plots and lichen 

species (see 1 through 4 above) should be 

relatively free of geographic bias such as 

massive over- or underrepresentation of 

geographic areas or ecological macrohabitats. If plot density is sufficient, balancing 

the model data set across the demonstrated major environmental gradients is 

recommended. A hypothetical example of 

adequately balanced representation of 100 

plots across three subdivisions of each of 

two major environmental gradients 

is illustrated in table 4. 

Table 4—Hypothetical allocation of 

plots in an exercise to balance 100 

plots across three categories each 

for two important environmental 

factors for a lichen gradient model 

data set

 Number

Elevation Pollution of plots

Low Low 10

Med Low 13

High Low 10

Low Med 10

Med Med 12

High Med 14

Low High 12

Med High 10

High High 9

Data to support model 

development must meet 

several criteria.
22

GENERAL TECHNICAL REPORT PNW-GTR-807

5. Forest mensuration data for all plots. At minimum, total plot live basal area 

and percentage of live basal area in hardwoods or softwoods are needed. 

Such data for standard FIA plots are available from program databases; 

equivalent data must be collected for all supplemental plots.

6. Climate and topographic data for all plots. This is best obtained with public databases and geographic information system coverages. Location and 

elevation are available for standard FIA plots from program databases.

7. Air quality estimates for all plots. This can be as simple as a binary value; 

plot has polluted air (1) or not (0). It is preferable to have quantitative estimates of SOx, NOx, and NH4 for each plot. Quantitative emissions data for 

point sources near “poor air quality” supplemental plots and interpolated 

estimates of regional pollution can be obtained from public sources such 

as state databases, the National Air Deposition Program (NADP 2009), 

the National Acid Precipitation Assessment Program (NAPAP 1991), the 

Community Multiscale Air Quality System (CMAQ 2009), or other sources 

such as Fenn et al. (2003b). Available estimates should be carefully evaluated before being included in model data sets. For instance, McCune et al. 

(1997b) evaluated the available modeled estimates for air pollutants in the 

Southeast Lichen Model Region and rejected them in favor of a simple polluted/not polluted plot assignment.

8. If lichen tissue element content is also available, it is a tremendous asset to 

the set of data for environmental factors. This provides accurate and consistent data measured at plots to represent concentration of air pollutants 

(Geiser and Neitlich 2007, Jovan 2008, Nimis et al. 2002) for comparison 

with variation in lichen species composition and with modeled or estimated 

air quality data from off-plot monitoring networks. Data for S, N, and 

heavy metals are most useful. Tissue element data from all plots are desirable, though data from a subset of plots can be useful. Tissue element data 

are most reliable when all data are from the same lichen species. Finding 

one species at all plots is difficult for a large-scale project, so cross-species 

calibration is often necessary. There are a number of constraints limiting 

selection of lichen species for tissue element analysis; additional field sampling cost, time, and training need to be considered as well. Until an FIA 

Program document with recommendations is produced, researchers should 

consult the current FIA lichen indicator advisor and the most recent scientific literature when incorporating tissue analysis into a model development 

project.
23

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Methods and Tools

In this section, we evaluate in detail two methods for building gradient models and 

briefly review available tools. Unconstrained ordination has been used successfully 

to develop at least part of all models completed to date. Identification of indicator 

lichen species coupled with regression modeling has been used to successfully 

define air quality gradients independent of climate for three of the models developed to date, where unconstrained ordination alone was not adequate. Researchers 

are encouraged to create and/or evaluate other methods and tools for possible use 

in developing lichen gradient models for the FIA Program. Any additional method 

or tool must be shown to meet the criteria stated in subsections “Required performance standards,” p. 11, and “Required products,” p. 12, and explained in detail in 

the next two subsections below before it is used to develop a final lichen gradient 

model. We include detailed discussion of how criteria are met for each step, to make 

transparent how each decision satisfies the purpose, rationale, and specific criteria 

for model development. “Scores” are the numerical values for a plot generated 

during analysis. These become lichen response “indexes” as their interpretation is 

defined by the analysis procedure.

Unconstrained ordination—

Currently the most widely used, most robust, and most effective unconstrained 

ordination technique for multivariate data reduction and identification of important 

gradients of species composition, when faced with heterogeneous community data 

sets, is NMS (or NMDS: Kruskal 1964, Legendre and Legendre 1998, McCune and 

Grace 2002, Minchin 1987). An NMS ordination uses a species-by-plot data set and 

is well-suited to species data that are non-normal or are on arbitrary or discontinuous scales (Mather 1976). The NMS can be used both as an ordination method and 

as a technique for assessing the dimensionality of a data set. 

Other unconstrained ordination techniques and other multivariate modelbuilding techniques may be acceptable if those techniques meet all requirements 

for the program and are shown in published literature to be robust and effective for 

FIA Program purposes. For instance, use of unconstrained ordination to develop 

a model for lichen response to important environmental gradients followed by 

constrained ordination (Jongman et al. 1995, Roberts 2008) to refine that model 

might be acceptable if all the requirements described in subsection “Required 

performance standards,” p. 11, can be met. 

A critical component of any ordination technique is the choice of a distance (or 

dissimilarity) measure; an algorithm that calculates a single number representing 

the difference between any pair of plots in the data set. The set of distances between 
24

GENERAL TECHNICAL REPORT PNW-GTR-807

all pairs of plots is the internal data set used by an unconstrained ordination algorithm to generate a multivariate model. The currently most robust and popular 

distance measure for use in NMS with a heterogeneous species-by-plot data set (one 

with many zeros for species in plots) is the Sørenson distance measure applied to 

species abundance data, also known as Bray-Curtis distance (Beals 1984, Faith et 

al. 1987). Other mathematically similar distance measures have been reported in the 

literature and would also be appropriate for such use. Critical characteristics for a 

robust distance measure in these circumstances are that (1) both rare and common 

species contribute substantially to the calculated distance (i.e., species abundance 

terms have no exponents in the distance algorithm), (2) the abundances of species 

common to both plots are standardized by the abundances of all species at either 

plot, and (3) joint absences from a pair of plots are not factored into the distance 

algorithm (Legendre and Legendre 1998, McCune and Grace 2002, Roberts 2008).

Ordination techniques that incorporate Reciprocal Averaging/Correspondence 

Analysis algorithms do not allow the user to select a distance method; the 

algorithms used imply a chi-squared distance measure (Jongman et al. 1995). 

This distance has been shown to be less robust and more subject to bias than the 

Sørenson-abundance/Bray-Curtis distance (Faith et al. 1987) for heterogeneous 

community data. So this entire class of ordination techniques (including also 

the popular canonical correspondence analysis and redundancy analysis) is less 

appropriate for lichen gradient modeling than are techniques that either allow the 

user to choose a more appropriate distance measure or incorporate a distance 

measure related to those recommended above.

Species abundance data may be relativized so that all plots have the same 

total abundance, or not, depending on the variation between plots in total species 

abundance. If variation in total lichen abundance at plots is moderate to high (coefficient of variation of total plot abundance is at least 40 to 50 percent), relativization 

removes the substantial noise from this source and enhances the signal from differences between plots based on variation in species composition (McCune and Grace 

2002), which is desirable for developing a model useful for FIA Program purposes. 

Several available proprietary or shareware (R packages, for instance) statistical 

or multivariate analysis software packages provide algorithms for NMS ordination 

and other ordination techniques. Many allow selection of the Sørenson-abundance/

Bray-Curtis distance or a closely-related distance algorithm, and most allow data 

modification choices consistent with our recommendations. A tested software package used for published analyses and reviewed in recent published literature should 

be selected to reduce the chances that unrecognized software errors will affect the 

final lichen gradient model. 
25

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

For the first lichen gradient model, developed for the Southeast Lichen Model 

Region, two lichen response indexes were developed using only NMS ordination, 

the Sørenson-abundance/Bray-Curtis distance, a data set of 203 plots, and species abundance data not relativized by plot (the coefficient of variation for plot 

abundance was much less than 50); that model accounted for 58 percent of original 

variation in distances between lichen plots (McCune et al. 1997b). The Pacific 

Northwest West model (Geiser and Neitlich 2007) was also successfully developed 

using NMS alone. A model for the California Sierras is described as a hybrid NMS 

model because two separate NMS ordinations using notably different data sets were 

required. The lichen air quality index applies to the California Sierras only (Jovan 

and McCune 2006), but the two lichen climate indexes apply to most of California 

(Jovan and McCune 2004). All of these models were developed using the PC-ORD1

software package (McCune and Mefford 2006). 

We describe the Pacific Northwest West lichen gradient model (Geiser and 

Neitlich 2007) in greater detail to illustrate both general and unique characteristics 

of one successful NMS lichen gradient model. From a data archive of more than 

1,500 sampled plots, authors randomly selected 293 plots to balance the analysis 

data set for polluted versus clean air plots stratified across two climate/topographic 

gradients. The availability of so many plots with lichen data is likely to be unique 

to the Pacific Northwest region. Lichen species abundance data were not relativized 

by plot because variation among plots in lichen abundance was relatively small. The 

amount of such variation differs by region, and recommended practice is related to 

the amount of variation. Tissue element content from a single lichen species gave 

the best quantitative estimators for air pollution impact at a plot. Such data are not 

universally available for FIA lichen plots, nor is it likely that a single lichen species suitable for tissue element analysis would be found at each plot in most model 

regions. The final model illustrated in figure 3 explained about 83 percent of the 

original variation among plots in lichen species composition. Three axes (gradients) were extracted and two lichen response indexes were developed from the 

model. Axis 1 on the ordination joint plot (fig. 3), representing a lichen air quality 

response gradient, and axis 2, representing a lichen climate response gradient, each 

explained about 37 percent of original variation among plots. This is much better 

than the minimum expectation for a successful lichen gradient model as stated in 

“Required performance standards,” No. 2, p. 11. The third and undefined lichen 

gradient from this model, represented by axis 3 in figure 3, might be considered 

1

 The use of trade or firm names in this publication is for reader information and does 

imply endorsement by the U.S. Department of Agriculture of any product or service.

The PNW West model 

illustrates features of 

a successful lichen 

gradient model derived 

using ordination.
26

GENERAL TECHNICAL REPORT PNW-GTR-807

Figure 3—Nonmetric multidimensional scaling (NMS) ordination diagrams for the Pacific Northwest West lichen gradient model 

(PNWW). These three two-dimensional (two axes) graphs illustrate plot distribution on the three NMS axes (lichen response gradients) 

developed for the PNWW model and their interpretation. Axis 1 represents an air pollution response gradient, axis 2 represents a climate 

response gradient, and axis 3 represents an undefined lichen community response gradient. Each axis label includes the percentage 

of total information about differences between plots that is represented on that axis. Lines superimposed on each graph represent the 

correlation of the named factor with the two axes. In such an ordination joint plot (see glossary) from an NMS model, the information 

portrayed by these lines is external to the ordination analysis (secondary information). Length of line from the central point indicates 

strength of correlation, and direction of line indicates both the sign of the correlation and the degree to which the factor is correlated with 

both axes. A line strictly parallel with one axis has no significant correlation with the other axis. An axis is assigned meaning based on 

its correlation with these factors. Coding of plot symbols and interpretation of names for lines are given in the legend. Values assigned 

from the model for an FIA plot represent position along axes 1 and 2, and become the lichen air quality index and the lichen climate 

index, respectively, for that plot. Modified from Geiser and Neitlich (2007). 

N Plagla

Elev Lon

Continentality

Ann Temp

Coast Dis

Rel Hum

Min Dec Temp

S Plagla

Axis 1 37 percent

Axis 2 37 percent

Elev

Lon

Ann Temp Continentality

Coast Dis

Rel Hum

Min Dec Temp

S Sum

Axis 2 37 percent

Axis 3 9 percent

N Plagla S Plagla

S

Sum

Axis 1 37 percent

Axis 3 9 percent

Variables:

 Clean air plot 

 Polluted air plot 

N Plagla: nitrogen content (percentage dry weight) 

of the lichen Platismatia glauca

S Plagla: sulfur content (percentage dry weight) 

of P. glauca

Elev: elevation

Coast Dis: distance east from the ocean 

Lon: longitude

Continentality: continentality index

Rel Hum: mean annual relative humidity (percent) 

Ann Temp: mean annual temperature

Min Dec Temp: mean minimum December 

temperature 

S: number of lichen species recorded on plot

Sum: sum of all species’ abundance codes on plot
27

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

unimportant because it explains only 9 percent of original variation. It is however 

the only gradient from that model that is correlated with both number of lichen species at a plot and total lichen abundance. If a causal factor could be defined for this 

gradient, the lichen species richness index could become a more useful bioindicator 

for this region. For all other completed lichen gradient models, the lichen species 

richness index has been strongly correlated with at least one of the defined lichen 

response indexes, so this situation for the Pacific Northwest West lichen gradient 

model is unusual.

Calculate plot indexes from an unconstrained ordination model—

There is currently only one tool widely available to calculate lichen response 

indexes for new plots or repeat samples of plots from an existing ordination model. 

This tool, submodule NMS Scores under module Ordination in software package 

PC-ORD (McCune and Mefford 2006), is an iterative procedure to assign to a set 

of plots numerical values for each ordination axis (lichen response gradient) of an 

existing model, using an algorithm based on the mathematics of NMS. This tool 

uses lichen abundance by species for each new plot and requires the choice of a 

distance measure. The recommendation is to use the same distance measure used 

to develop the model. A data set of the lichen abundance by species in the model 

for each plot in the model (the model data set) and a data set of numerical values 

(scores) of each model plot for each model axis are also required. The software also 

calculates estimates of fit of each new plot to the original model, allows the user to 

specify criteria for poor fit, and provides flags for three kinds of poor fit. The score 

for each plot on each model axis (gradient) then becomes its lichen response index 

for that gradient (formally related to a named environmental factor). This tool was 

developed specifically to calculate scores for new plots along axes (defined environmental response gradients) of an existing NMS model. It could possibly be used to 

calculate scores for new plots along axes (gradients) from any kind of unconstrained 

ordination model, and possibly even along axes (gradients) from a constrained 

ordination model, but this should be tested thoroughly before being implemented 

for a new lichen gradient model. Other software to accomplish this purpose can 

be recommended by gradient model developers, as long as it has been thoroughly 

tested. A possible future option is development of FIA Program–specific software 

(by program staff or through contracting) to assign lichen response indexes to new 

plots using any kind of ordination model that meets program specifications.

Regression model using indicator species—

Regression modeling after selection of indicator species has been a successful tool 

used in a hybrid lichen gradient model when a single ordination model does not 
28

GENERAL TECHNICAL REPORT PNW-GTR-807

generate independent lichen community response gradients that can be linked to 

climate and air pollution. So far, the problem has been to derive a second gradient 

representing lichen response to air quality that is independent of a strong primary 

lichen climate response gradient. The following procedure has proved successful. 

First a set of lichen species that indicate poor air quality has been developed using 

Indicator Species Analysis (Dufrêne and Legendre 1997) to identify species that are 

significantly associated with plots having poor air quality. Other means to identify 

significant indicator species for one extreme of a particular gradient can also be 

used for this approach, subject to meeting requirements and guidelines in subsection “Program Requirements,” p. 11. To be useful across an entire lichen model 

region, probably at least 6 to 10 lichen indicator species must be identified. Next a 

raw lichen air quality (or other factor) response score for each plot is calculated. A 

further step, if necessary, is to calculate an independent adjusted lichen air quality 

index based on a regression or other statistical model that includes both the raw 

lichen air quality score and the confounding primary climate-related factor. A 

linear regression model to express the adjusted lichen air quality index for a plot as 

the residual after regressing the raw air quality score on the primary and confounding climate factor would be a simple example for this further analysis step. Such 

a simple example is illustrated in detail in section “NMS ordination plus indicator 

species plus regression,” p. 33. All major statistical software packages (including 

shareware R packages) provide appropriate tools for developing various kinds of 

regression models. 

For the construction of a lichen gradient model for the state of Colorado 

(McCune et al. 1998), an NMS two-dimensional ordination of 126 lichen plots 

clearly showed that climate as represented by elevation was strongly correlated with 

the most important variation in lichen composition, generating a lichen climate 

index. Lichen composition was also strongly correlated with an air quality factor, 

but this factor was correlated as well with the lichen climate/elevation gradient. 

This meant an independent lichen air quality index could not be developed using 

NMS alone. A list of indicator species was extracted using module “Indicator 

Species Analysis” in PC-ORD (McCune and Mefford 2006), then a lichen air 

quality response index independent of climate was calculated from a simple linear 

regression model using a proprietary statistical software package. For this hybrid 

model, the lichen climate index explained 37 percent of original variation among 

all plots in the model, and the regression model for the lichen air quality index 

explained 35 percent of variation (had an adjusted r2

 of 0.35) among a subset of the 

plots. This is better than the minimum expectation for a successful hybrid lichen 

gradient model as stated in “Required performance standards,” No. 2, p. 11.
29

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

The need for hybrid models may occur in several West lichen model regions 

where air pollution is worse at lower elevations. For the California Sierras hybrid 

model (Jovan and McCune 2006), an independent lichen air quality index was 

developed from an indicator species/regression model approach, and two independent climate indexes were developed from an NMS model for most of California 

(Jovan and McCune 2004). Similar problems may also occur in parts of the East 

where air pollution is strongly correlated with latitude or longitude. 

Steps in Building a Lichen Gradient Model

One should always first attempt to develop a lichen gradient model using NMS 

ordination or another acceptable multivariate analysis procedure, then proceed to 

a hybrid model only if the multivariate analysis model is not completely successful 

(see discussion in subsection “Rationale for Approach,” p. 9. The steps below are 

described for NMS, indicator species analysis, and regression modeling, techniques 

that have been successful for completed FIA lichen gradient models. A lichen gradient model developed using other techniques should have modifications of these 

steps documented as part of the final report.

Multivariate NMS ordination model—

1. Compile a lichen species-by-plot model data set including both standard 

and supplemental plot data that represent the range of very good to very 

poor air quality conditions across the full range of other important environmental conditions (elevation, temperature, precipitation, etc.) in a region 

(see subsection “Data Needs” above, p. 20). For FIA data, modify the lichen 

species names in the final data set to match usage in the earliest inventory 

year included in the data set, following instructions in the official FIA reference lichen species comments table. This standardizes taxonomic usage 

across all years the data were collected, and is necessary only if taxonomic 

usage changed during those years. The latest version of the official FIA 

table must be obtained from the internal program database; contact the FIA 

lichen indicator advisor or an FIA information management office. The version available on the Web (USDA FS 2010a) is not as current, but may be 

adequate depending on the inventory years for the lichen data being used. 

Second, further modify the lichen species names as needed for analysis, for 

instance by combining several related and ecologically similar uncommon 

species into a single species group. Perform this second kind of modification of lichen species categories as little as possible, to facilitate calculation 

of indexes for new plots from the completed model as easily as possible. 

Any modifications to use of lichen species categories applied to the model 

Details matter for 

developing a model.
30

GENERAL TECHNICAL REPORT PNW-GTR-807

data set must also be applied to all new plots to be assigned lichen response 

indexes from the model in future. Examine the original lichen data set for 

outlier plots, plots having very few species, and species found at very few 

plots; remove plots or species as desired to reduce noise and generate a 

more stable solution. This is standard practice in development of ordination 

models (Legendre and Legendre 1998, McCune and Grace 2002).

2. Compile a second matrix of environmental and other variables associated 

with each plot in the model data set. Variable classes represented in the 

second matrix should include climate, pollution, tissue element content if 

available, forest structure measurements, and lichen community summary 

values such as lichen species count and summed lichen abundance for each 

plot. Consider including as separate variables logarithmic or other transformations of original environmental data; response of biotic communities to 

such factors is often not linear with respect to the scale in which scientists 

measure an environmental factor (Legendre and Legendre 1998, McCune 

and Grace 2002). Identifying such instances is important to accurate interpretation of lichen response. 

3. Perform an NMS ordination using many independent runs (up to 50 is common) at each of several dimension settings. Choose an appropriate distance 

measure such as the Sørenson-abundance/Bray-Curtis distance measure 

(see subsection “Unconstrained ordination,” p. 23). Determine the dimensionality of the data set by comparing solutions at different dimensions, and 

choosing the dimensionality that on average deviates from random most 

strongly. Test this by comparing with multiple runs of randomized data. 

Next, perform multiple NMS runs (often up to 100) with the chosen number 

of dimensions to ensure that the chosen final solution is stable and represents a configuration with the lowest possible deviation from the matrix 

of original pairwise distances between plots (usually called lowest stress) 

as reported by the software. All of these tasks can be performed in several 

available multivariate analysis program packages (including shareware R 

packages). 

4. Rotate the chosen solution to maximize the correlation of ordination 

axes with the second-matrix environmental factors of greatest interest 

to address assessment questions (a “climate” axis and an “air quality” 

axis, for instance). These axes then represent lichen response to those 

environmental factors in the model, and plot scores on each axis become 

the indexes of lichen response to those factors. The ordination axis most 
31

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

strongly correlated with air quality should be rotated so that plots reflecting 

poorer air quality (more pollution) have higher values on the axis. This 

convention facilitates interpretation of lichen air quality indexes from 

this model as estimating potential risk of injury from air pollution. It is 

useful for presentation but not necessary for analysis to rotate the solution 

so that axis 1 expresses the most variation among plots attributed to any 

axis. The final solution must meet all the requirements stated below in 

subsection “Evaluating a Lichen Gradient Model During Development,” 

p. 37, including both strength and independence of correlations with 

environmental factors of interest. Many, but not all, available multivariate 

analysis program packages (including shareware R packages) allow such 

rotations. If the ability to perform such a rigid rotation of a two-, three-, or 

four-dimensional point cluster is not available in the multivariate analysis 

package used, a graphics software package can be used to perform the 

rotation. Model developers must then statistically evaluate the proportion 

of original variation among plots on each final rotated axis in a separate 

operation.

5. Archive the final model data set and the set of coordinates of all plots on 

the final set of rotated axes chosen from step 4. Together these two data sets 

compose the final ordination gradient model. Complete the template for 

documentation of a gradient model, including all information and instructions needed to assign to new plots or to resampled plots response indexes 

for each defined gradient from the model.

6. Investigate correlations of second-matrix variables and species abundance 

values with the axes in the final gradient model (step 5) to support full 

interpretation of the final gradient model and variation in lichen community composition in the region. Calculate correlations of abundance of each 

individual lichen species with final axes to interpret a species’ response 

along climate and air quality gradients. Report the proportion of original variation among plots on each final model axis. Lists of lichen species 

characteristic of different climate subregions and of species characteristic 

of low and high air quality are also useful products. Archive an ordination 

joint plot showing second matrix variables plotted over the basic ordination 

diagram as lines representing strength of correlation of second matrix variables with ordination axes (see glossary).

7. Calculate lichen response indexes for new FIA plots as well as for plots 

included in the model data set using the original plot data (not the final 
32

GENERAL TECHNICAL REPORT PNW-GTR-807

model data set) as modified following instructions for species categories. 

Such calculations require the archived model data set and the archived set 

of lichen model axis scores. Follow all instructions in the model documentation for preparation of the data set of new plots and selection of options 

for assigning scores to plots. Calculations can be done using submodule 

NMS Scores in PC-ORD (McCune and Mefford 2006), which derives 

indexes for a set of plots by using the lichen gradient model and finding 

the position of best fit for each of these plots on each model axis without 

shifting the positions (scores) of the archived set of model plots. Additional 

software may be available in the future for this purpose. In PC-ORD/

NMS Scores, if more than one lichen response index is to be derived from 

an NMS model, choose option “simultaneous axes” to optimize overall 

fit of plots to all ordination axes in the model (McCune and Grace 2002). 

If only one response index is to be derived from the NMS model, choose 

option “one axis at a time.” NOTE: the archived scores for model plots are 

based on a data set modified for model development. Calculation of lichen 

response indexes for these same plots should be based on the original plot 

data rather than data modified for the model, so these indexes are equivalent to indexes calculated for new plots and repeat samples of plots (full plot 

data not included in the archived model). 

8. Estimate how well the model applies to each plot. Submodule NMS 

Scores in PC-ORD reports two kinds of user-specified poor fit of a plot 

to a model: too far beyond an extrapolation limit for either end of the 

model gradient, and exceeding a final stress criterion (calculated using 

the same algorithms as for stress in NMS ordination (McCune and Grace 

2002). Recommendations for FIA lichen gradient models are to set the 

extrapolation limit criteria to 10 percent of original model axis (gradient) 

length beyond either end of the original gradient. Any plot whose score 

is beyond this range (+ or –) is flagged for poor fit. For stress-based fit 

calculated in NMS Scores, it is recommended that the final stress criterion 

be selected to be two standard deviations above the average plot stress 

parameter as calculated for a population of at least 100 standard FIA plots 

from the same model region. Recommended extrapolation limits and the 

selected final stress criterion are then archived and included with model 

documentation, and are entered into NMS Scores before calculation of 

scores for new plots. Any plot whose calculated fit is outside the entered 

criteria is flagged for poor fit. For plots fitted to an NMS lichen gradient 
33

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

model using other software, or for plots fitted to a lichen gradient 

model developed using a different multivariate analysis approach, the 

extrapolation limit criteria above can easily be applied. A new softwarespecific criterion for defining stress-related poor fit to the model would 

need to be established and documented as part of model development.

Hybrid model: NMS ordination plus indicator species plus regression—

1. Perform steps 1 through 3 as for an NMS ordination lichen gradient model 

above. If the best NMS solution fails to generate a model with a lichen air 

quality gradient that is independent of the major lichen climate gradient(s), 

then consider building a hybrid model. The steps below are written for 

deriving an independent lichen air quality index in this situation using “pollution” indicator species and a simple linear regression model, the simplest 

combination of methods published so far. Either choice might be different for an actual lichen model region. Specific examples are given from 

the Colorado model (McCune et al. 1998) and the California Sierras model 

(Jovan and McCune 2006). If instead there is a need to derive an independent lichen response index of some other kind, follow the steps below as 

modified for the relevant environmental factors. 

2. Examine the lichen climate response gradient (primary gradient) with 

which the air quality factor (or other secondary environmental factor) is 

most closely confounded. Create a “subset” lichen species-by-plot data set 

by removing all nonpolluted plots that have values for the climate factor most strongly correlated with the primary gradient that are outside the 

range of values for that same factor in the set of “polluted air” plots (or the 

range of another secondary environmental factor). 

3. Develop a list of lichen indicator species using the subset of lichen data and 

an appropriate software package (see subsection “Regression model using 

indicator species,” p. 27) as well as other approaches to identify lichen species that have strong indicator value for identifying plots in that region that 

have poor air quality. Each species included must have a quantitative estimate of its strength as an indicator in that region. NOTE: just as a different 

NMS lichen gradient model must be developed for each different bioclimatic region, a different set of “pollution” (or another environmental factor) indicator lichen species must be derived for each different bioclimatic 

region. 
34

GENERAL TECHNICAL REPORT PNW-GTR-807

4. Develop a formula to calculate a raw lichen air quality score for each plot 

based on the entire lichen species composition at a plot. In this example, it 

is accomplished by calculating the relative abundance of pollution indicator 

lichen species at the plot:

Raw lichen air quality score = 100 × (SUMpoll/SUMall)

where SUMpoll is the sum of the abundance scores of all lichen species 

that indicate polluted conditions and SUMall is the sum of the abundance 

scores of all lichen species recorded for the plot. This raw lichen air quality 

score has a maximum of 100 when all lichen species are pollution-tolerant 

indicators. The minimum possible value is zero when none of the lichen 

species present are pollution indicators. Note this raw score is calculated so 

that plots with more polluted air have high scores, to facilitate interpretation of the final lichen response index as estimating potential risk of injury 

from air pollution. This is the recommended strategy. A similar approach 

was used to calculate raw lichen air quality scores for the California 

Sierras model (Jovan and McCune 2006), whereas for the Colorado model 

(McCune et al. 1998) the original raw lichen air quality score was highest 

for clean air plots. 

5. Include the raw lichen air quality score (calculated for all plots in the original data set) as a variable in the environmental data set. Correlate with the 

best NMS model (step 1 in this section) to see if air quality as represented 

by the raw score is now decoupled from the confounding climate variable. 

If it is, this becomes the final lichen air quality index. If it is not, which is 

to be expected and as was the case with both the Colorado and California 

Sierras lichen gradient models, proceed to the next step.

6. Raw lichen air quality scores can be adjusted to reduce the influence of the 

confounding climate variable by using a regression approach. In the next 

step for the example introduced in step 4, it is done by calculating residuals 

from linear regression of raw lichen air quality scores on the confounding 

climate-related factor:

Adjusted lichen air quality index =

 [(Raw air quality score – (a × Climate factor + b)]/SD

where a is the slope of the regression equation, b is the intercept, and SD 

is the standard deviation of the residuals from a regression model with raw 

air quality score as the dependent variable and the climate factor as the 

independent variable. The adjusted index expresses the lichen response 
35

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

to air quality of a plot as the number of standard deviations away from 

expectation for a given value of the climate factor. If more than one independent variable is entered into a regression model, forced simultaneous 

entry rather than any form of stepwise entry should be used for the final 

regression model to confirm there has been no bias from order of entry. For 

the Colorado model, a linear regression model based on plot elevation was 

adequate, but for the California Sierras model, a more complicated regression model (also based on plot elevation) was required. NOTE: the choice of 

climate-related factor to include in a regression model is constrained for use 

in the FIA Program by the need to have values for that factor for all sample 

plots to which the model will be applied (see “Required performance standards” No. 6, p. 12). In practice, this means that latitude, longitude, elevation, and combinations of these factors are the primary choices, as location 

and elevation are the only climate-related factors routinely available for FIA 

plots. 

7. Include this adjusted lichen air quality index (calculated for all plots in 

the original data set from either all lichen species or only those in the 

model data set, as directed in model documentation) as a variable in the 

environmental data set, and compare with the best NMS model (step 3 

above) to see if adjusted lichen air quality index is now decoupled from 

the confounding climate index. It is expected to be, as was the case with 

the Colorado and California Sierras Lichen Gradient Models. This then 

becomes the final lichen air quality response index. The final hybrid model 

must meet all the requirements stated below in subsection “Evaluating 

a Lichen Gradient Model During Development,” p. 37, including both 

strength and independence of correlations with environmental factors of 

interest. 

8. Archive the final NMS model data set, the final NMS model axis scores, 

the list of air quality indicator lichen species for the region, the list of all 

plots and lichen species included in the subset lichen air quality gradient 

data set, and the formula to calculate the adjusted lichen air quality index 

for new plots and repeat samples of plots including all needed numbers such 

as standard deviation of regression residuals, all of which together compose 

the hybrid lichen gradient model for that region.

9. Assign lichen response indexes for the climate factor(s) to all plots included 

in the original lichen data set plus any additional standard plots based on 

full original plot lichen data. Use recommended software and all necessary 
36

GENERAL TECHNICAL REPORT PNW-GTR-807

final model data sets and products listed in instructions for applying the 

model. Follow all instructions in the documentation for preparation of the 

data set of plots to be assigned lichen indexes and selection of options in the 

software to be used (see step 7 in subsection “Multivariate NMS ordination model,” p. 31). If using routine “NMS Scores” in software package 

PC-ORD (McCune and Mefford 2006), choose option “simultaneous axes” 

to calculate two or more lichen response indexes from an NMS model, or 

choose option “one axis at a time “ to calculate a single lichen response 

index for each plot (McCune and Grace 2002).

10. Assign lichen air quality indexes (or other indicator species/regressionbased lichen response indexes) to all plots included in the original lichen 

data set plus any additional standard plots based on full original plot lichen 

data by calculating relative abundance of “pollution” indicator lichen species at a plot (steps 3 and 4), then calculating adjusted lichen air quality 

index (step 6) for each plot. Follow all instructions in the documentation for 

preparation of the data set of plots to be assigned lichen response indexes 

and selection of options for assigning indexes to plots.

11. Estimate how well both parts of a hybrid model apply to each new plot. 

Criteria for assessing how accurately the ordination part of a hybrid model 

generates a lichen response index for each plot are explained in step 8 of 

section “Multivariate NMS ordination model,” p. 32. Extrapolation limit 

criteria can be established in the same way for an indicator species/regression model. There is no standard approach for calculating a stress criterion 

for applying an indicator species/regression model to new plots. One possible estimate is the proportion of lichen species at a plot that are not found 

in the regression model data set. We have no established guidelines for 

interpreting fit to a regression model using this criterion to define stressrelated poor fit. Analysis of statistics for fit of plots to such models should 

be a priority in the near future to aid in developing such guidelines.

12. As an alternative to using pollution indicator lichen species in step 3, the 

above sequence may be attempted by identifying clean air indicator lichen 

species. Such an approach may be helpful in a lichen model region that has 

generally relatively clean air or for special studies (e.g., evaluation monitoring) where few or no pollution indicator lichen species may be found over 

broad parts of the area of interest. This approach is likely to generate useful 

lichen gradient response indexes primarily for areas with relatively low ecological variation, as clean air lichen indicator species must be found in most 
37

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

habitats in the region to be reliable indicators. The use of pollution or clean 

air indicator lichen species will probably result in less ability to discriminate on the end of the lichen air quality response gradient not represented 

by indicator lichen species.

Evaluating a Lichen Gradient Model During Development

An effective NMS lichen gradient model (including all gradients defined) should 

explain at least 50 to 60 percent of the variation among distances between pairs of 

plots in the model data set. This criterion should be applied to the NMS ordination 

component of a hybrid gradient model as well. Such a statistic is almost always 

reported by ordination software; its exact calculation differs by ordination technique and by software. Report the basis for this calculation in the software used 

to develop the model. If the model does not achieve this goal, then two remedies 

should be explored. Both remedies require iterative alterations and model testing. 

a. The data set can be enlarged, by adding more years of standard plots, by 

collecting more supplemental plot data, or both. 

b. The boundaries of the model region can be adjusted. 

An effective lichen gradient model should have correlation r2

 values of 0.3 to 

0.4 or higher between defining variables (climate or air quality factors) and major 

axes of the ordination. The factors used to define and align the ordination axes 

that will generate lichen response indexes must be relatively uncorrelated with one 

another (correlations not significant or correlations closer to zero than ± 0.5). If this 

is so, the major factors will appear close to orthogonal on an ordination joint plot. 

Failure to achieve at least 75 to 80° (maximum 90°) separation (geometric degrees) 

or ≥85 percent (maximum 100 percent) orthogonality for both ordination axes and 

the major environmental factors correlated with them should trigger use of a hybrid 

model development approach. For a hybrid model, the lichen response index derived 

from indicator species selection and regression modeling must achieve the above 

level of separation from all ordination axes used to define lichen response indexes. 

If this cannot be achieved with the steps outlined above, a new procedure for a 

hybrid model must be developed.

With an effective lichen gradient model, we have found that if a field crew 

person obtains 65 percent or more of the count of species found on that plot by a 

lichen specialist, then the lichen response indexes from field crew data for that 

plot will mostly fall within 10 to 12 percent of the indexes from expert data, as 

compared with model gradient length (see glossary). With within-plot variation at 

10 to 12 percent or less, it is reasonable to expect that deviations in lichen response 

index over time of 15 percent or more of model gradient length suggest real change 
38

GENERAL TECHNICAL REPORT PNW-GTR-807

has occurred (Will-Wolf 2010). The ability of lichen data meeting the field measurement quality objective (MQO) to generate lichen response indexes that meet 

laboratory MQOs (table 5) needs to be demonstrated for each new lichen gradient 

model. Failure to achieve this relationship indicates the model is not adequate and 

an improved model must be developed.

Field crew certification exams by specialists conducted during training and data 

quality checks during the field season (Will-Wolf 2007) provide the data needed to 

test whether achievement of MQOs for field data generate lichen response indexes 

that meet laboratory MQOs. If such data (repeat samples of the same plot) from 

training or data quality assessment are not available for a region, a reference plot 

exercise should be conducted concurrent with lichen gradient model development, 

to provide appropriate data for this test.

Evaluating, Comparing, and Calibrating 

Completed Lichen Gradient Models

Evaluating a Completed Lichen Gradient Model

A lichen gradient model should be expected to have a finite lifespan, beyond which 

it is no longer valid. Reorganization of biotic communities over time in response to 

climate change or other large-scale alterations of environmental context is a distinct 

possibility in today’s world. The lifespan of a lichen gradient model is yet to be 

tested; we hope that models will remain valid for more than 50 years, or about five 

current sampling cycles. The validity of the lichen gradient model can be tested 

after each sampling cycle is completed.

Every completed lichen gradient model includes instructions for assessing how 

well plots fit with the model. When a new plot is assigned lichen indexes from 

an NMS gradient model using module NMS Scores in PC-ORD (see step 8 in 

subsection “Multivariate NMS ordination model,” p. 32), plots are flagged for two 

different types of poor fit: too far outside the range of the original model gradient 

at either end (fit flags 2 and 3), plus unspecified poor fit (fit flag 1). For plots fitted 

to an NMS lichen gradient model using other software, or for plots fitted to a lichen 

gradient model developed using a different multivariate analysis approach, similar 

flags for poor fit are also reported. For a model component developed using an 

indicator species/regression model approach, methods for estimating how well a 

plot can be fitted to the model are less well established (see step 11 in subsection 

“Hybrid model: NMS ordination plus indicator species plus regression,” p. 36.

If for a particular lichen model region the proportion of plots with poor fit 

increases substantially over time, this may suggest that the lichen gradient model is 

losing its validity. A suggested guideline is that if the percentage of all lichen plots 

We hope a model will 

be useful for 50 years.
39

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Table 5—Measurement quality objectives (MQO) for Forest Inventory and Analysis Program lichen data, and 

their method of assessment

Field and laboratory data MQO Method of assessment

Minimum plot standard (field MQO) 65 percent Percentage of a lichen specialist’s count of species that a 

crew person finds on the same plot (laboratory data). This is 

the minimum score for crew to pass training certification, 

field training recertification, and for minimum acceptable 

quality of field crew species (laboratory) data as compared 

 with specialists’ species (laboratory) data. 

Completeness 90 percent Percentage of forested plots searched for lichen data 

(field data).

Accuracy of archived laboratory data 95 percent Percentage of plots with archived lichen laboratory data 

 that have been 100 percent proofread and represent exactly 

 the plot lichen species count and abundance from the lichen 

identification specialist’s report.

Quality of archived laboratory data 90 percent Percentage of crew samples from plots achieving the 

 minimum plot standard for data quality: 65 percent of a 

 lichen specialist’s species count (laboratory data) based 

 on a blind resurvey of the same plot.

Lichen response indexes MQO Method of assessment

Bias 10–12 percent Signed deviation from "true" indexes (from lichen specialist 

 laboratory data). In practice, plot laboratory data achieving 

 the minimum standards (above) should yield indexes that 

 meet this MQO.

Accuracy 10–12 percent Absolute deviation (sign ignored) from "true" indexes (from 

 lichen specialist laboratory data). In practice, plot laboratory 

 data achieving the minimum standards (above) should yield 

 indexes that meet this MQO.

Precision 10–12 percent Deviation between indexes from measurement of the 

same plot in the same year by both field crew and lichen 

 specialists. Data that meet the MQOs for bias and accuracy 

 should meet this MQO. For pattern and trend analysis, 

 deviations of indexes between subregions and/or years are 

 considered evidence of differences or change only if they 

 exceed this precision standard.

Quality of archived lichen 95 percent Percentage of archived plot values for one particular lichen

response indexes response index that achieve the target MQO for accuracy, 

 bias, or precision. 

Notes: For a successful lichen gradient model, it is demonstrated during development that if the lichen species field and laboratory MQOs are met, 

MQOs for lichen response indexes are met. All four target MQOs for lichen response indexes require that a lichen gradient model be in place for 

their calculation.
40

GENERAL TECHNICAL REPORT PNW-GTR-807

for a year with flags for poor fit increases by 10 percent from one sampling cycle to 

the next or reaches 20 percent, the lichen gradient model should be reevaluated. The 

type of flags will help diagnose the nature of the poor fit and suggest remedies.

In several instances, plots from an area not included in a particular lichen 

model region may need to be assigned lichen response indexes by using the nearby 

existing lichen gradient model. For instance, northern Florida plots will in the 

future be assigned lichen response indexes from the Southeast Lichen Gradient 

Model, which included no Florida plots. A suggested guideline is that if the proportion of plots in the new area with any flag for poor fit is 10 percent or more when 

at least 50 plots are tested, that lichen gradient model does not adequately represent 

the new area. Remedies will differ depending on the location. If an area is located 

between or near two lichen model regions with completed models, geographically 

defined groups of plots should be assigned to the lichen model region whose model 

provides the best average fit to the entire group of plots.

Comparing and Calibrating Lichen Gradient Models

Calibrating between lichen response indexes developed for different regions would 

allow one to monitor trends in lichen response across multiple lichen model regions 

in a single analysis. The first requirement for such calibration is the development 

of lichen gradient models for adjacent lichen model regions. Measured air quality 

data for a range of plots in each region can be used to compare with ranges of lichen 

air quality indexes between regions. Data on lichen tissue element content can also 

be of help in comparing lichen air quality indexes between regions. Lichen climate 

indexes between regions can also be linked this way, and the relation of individual 

climate factors to lichen community composition can be compared between regions. 

Calibration of lichen response indexes between regions will also allow quantitative comparison of differential response of individual lichen species between 

regions, and will contribute much-needed improvement in the understanding of 

geographical variation in the effectiveness of environmental indicators. Evaluation 

and modification of lichen model region boundaries can also probably be conducted 

most efficiently by comparing and calibrating lichen gradient models for adjacent 

regions.

Summary Discussion

Development of a lichen gradient model for calculation of lichen response indexes 

to assess patterns and monitor trends in response to environmental factors has been 

an effective and powerful tool for the FIA Program to implement large-scale and 
41

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

long-term monitoring within tight budget constraints. This approach should be considered for other large monitoring programs. Developing the model is an exacting 

process. Choice of an appropriate lichen model region and modeling approach are 

critical first steps, but they do not guarantee success. Appropriate choices must be 

made at each of the many steps in the development of a model. Failure to do this at 

any step along the way can seriously impair the usefulness and scientific validity of 

the entire model. It is reasonable to expect that inappropriate choices at early stages 

in model development may have multiplicative rather than just additive impacts on 

the accuracy and precision of the monitoring that can be done with the final lichen 

response indexes. 

With this in mind, we have examined the rationale and techniques for each step 

in the modeling process in great detail. We have provided a reporting template that 

supports application of a completed model to monitoring for many years. We anticipate this attention to detail will make the model development process more transparent for future FIA Program lichen gradient models. We also hope this approach 

will clarify the model development process and the usefulness of our recommended 

approaches and techniques for researchers designing their own biomonitoring 

projects and programs.

Acknowledgments

Funding was provided under cooperative agreements between Southern Research 

Station Forest Inventory and Analysis (FIA) and the University of WisconsinMadison. Author Will-Wolf has been FIA Program lichen indicator co-advisor 

1995–2001, advisor 2002–2007, co-advisor 2008–2009, and is developing two FIA 

lichen gradient models. Author Neitlich has been FIA lichen indicator co-advisor 

1995–1998, advisor 1999–2001, has developed one completed FIA lichen gradient 

model, and is developing another model. 

This report is substantially modified from Neitlich et al. (2001), but takes much 

of its organization and content from that document. All other lichen gradient model 

developers, especially FIA lichen indicator developer Bruce McCune (Oregon State 

University [OSU]) and current FIA lichen indicator advisor Sarah Jovan (OSU and 

Pacific Northwest FIA), have contributed to the content of this document. Sarah 

Jovan, Paul Patterson (Interior West FIA), and Paul Rogers (Utah State University) 

provided insightful reviews of drafts of this report. The final product is the responsibility of the authors.
42

GENERAL TECHNICAL REPORT PNW-GTR-807

English Equivalents

1 hectare (ha) = 2.47 acres (ac)

1 kilometer (km) = 0.621 mile (mi)

1 meter (m) = 3.28 feet (ft)

1 square meter per hectare (m2/ha) = 4.37 square feet per acre (ft2/ac)

References

Bailey, R.G.; Avers, P.E.; King, T.; McNab, W.H. 1994. Ecoregions and 

subregions of the United States. [1:7,500,000]. Colored map. Washington, 

DC: U.S. Geological Survey. Accompanied by a supplementary table of map 

unit scale descriptions compiled and edited by McNab, W.H. and Bailey, R.G. 

Prepared for the U.S. Department of Agriculture, Forest Service.

Bates, J.W.; Farmer, A.M., eds. 1992. Bryophytes and lichens in a changing 

environment. Oxford, United Kingdom: Clarendon Press. 404 p.

Beals, E.W. 1984. Bray-Curtis ordination: an effective strategy for analysis of 

multivariate ecological data. Advances in Ecological Research. 14: 1–55.

Cleland, D.T.; Freeouf, J.A.; Keys, J.E.; Nowacki, G.J.; Carpenter, C.A.; 

McNab, W.H. 2005. Ecological subregions: sections and subsections for the 

conterminous United States. In: Sloan, A.M., tech. ed. Washington, DC: U.S. 

Department of Agriculture, Forest Service. [1:3,500,000]. Colored map. Also on 

CD-ROM as a geographic information systems coverage in ArcINFO format.

Community Multiscale Air Quality System [CMAQ]. 2009. CMAQ model 

version 4.7. U.S. EPA Atmospheric Science Modeling Division active opensource development project: suite of programs for conducting air quality model 

simulations. http://www.cmaq-model.org/. (September 2009).

Daly, C.; Taylor, G. 2000. United States average monthly or annual precipitation, 

temperature, and relative humidity 1961–90. Arc/INFO and ArcView coverages. 

Corvallis, OR: Spatial Climate Analysis Service, Oregon State University.

Dufrêne, M.; Legendre, P. 1997. Species assemblages and indicator species: the 

need for a flexible asymmetrical approach. Ecological Monographs. 67: 345–366.

Ellis, C.J.; Coppins, B.J.; Dawson, T.P.; Seaward, M.R.D. 2007. Response of 

British lichens to climate change scenarios: trends and uncertainties in the 

projected impact for contrasting biogeographic groups. Biological Conservation. 

140: 217–235.
43

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Faith, D.P.; Minchin, P.R.; Belbin, L. 1987. Compositional dissimilarity as a 

robust measure of ecological distance. Vegetatio. 69: 57–68.

Fenn, M.E.; Baron, J.S.; Allen, E.B.; Rueth, H.M.; Nydick, K.R.; Geiser, L.; 

Bowman, W.D.; Sickman, J.O.; Meixner, T.; Johnson, D.W.; Neitlich, P. 

2003a. Ecological effects of nitrogen deposition in the Western United States. 

BioScience. 53(4): 404–420.

Fenn, M.E.; Haeuber, R.; Tonnesen, G.S.; Baron, J.S.; Grossman-Clarke, S.; 

Hope, D.; Jaffe, D.A.; Copeland, S.; Geiser, L.; Rueth, H.M.; Sickman, J.O. 

2003b. Nitrogen emissions, deposition, and monitoring in the Western United 

States. BioScience. 53(4): 391–403.

Geiser, L.H.; Neitlich, P. 2007. Air pollution and climate gradients in western 

Oregon and Washington indicated by epiphytic macrolichens. Environmental 

Pollution. 145: 203–218.

Hawksworth, D.L.; Rose, F. 1976. Lichens as pollution monitors. London: 

Edward Arnold. 60 p.

Jongman, R.H.G.; ter Braak, C.J.F.; van Tongeren, O.F.R., eds. 1995. Data 

analysis in community and landscape ecology. Cambridge, United Kingdom: 

Cambridge University Press. 299 p.

Jovan, S. 2008. Lichen bioindication of biodiversity, air quality, and climate: 

baseline results from monitoring in Washington, Oregon, and California. Gen. 

Tech. Rep. PNW-GTR-737. Portland, OR: U.S. Department of Agriculture, Forest 

Service, Pacific Northwest Research Station. 115 p. http://www.fs.fed.us/pnw/

pubs/pnw_gtr737.pdf. (March 2010).

Jovan, S. 2009. Baseline results from the Lichen Community Indicator in the 

Pacific Northwest: air quality patterns and evidence of a nitrogen problem. In: 

Ambrose, M.J.; Conkling, B.L., eds. Forest Health Monitoring 2006 national 

technical report. Gen. Tech. Rep. SRS-117. Asheville, NC: U.S. Department of 

Agriculture, Forest Service, Southern Research Station. Chap. 5.

Jovan, S.; Geiser, L.H.; Neitlich, P. 2008. Pacific Northwest West lichen model 

region gradient model documentation and specifications for calculating lichen 

response indexes for new plots. Internal FIA report. 16 p. On file with: Sarah 

Jovan, U.S. Department of Agriculture, Forest Service, Forestry Sciences 

Laboratory, 620 SW Main, Suite 400, Portland, OR 97205.
44

GENERAL TECHNICAL REPORT PNW-GTR-807

Jovan, S.; McCune, B. 2004. Regional variation in epiphytic macrolichen 

communities in northern and central California forests. Bryologist. 107: 328–339.

Jovan, S.; McCune, B. 2005. Air-quality bioindication in the greater central valley 

of California, with epiphytic macrolichen communities. Ecological Applications. 

15: 1712–1726. 

Jovan, S.; McCune, B. 2006. Using epiphytic macrolichen communities for 

biomonitoring ammonia in forests of the greater Sierra Nevada, California. 

Water, Air, and Soil Pollution. 170: 69–93.

Kruskal, J.B. 1964. Nonmetric multidimensional scaling: a numerical method. 

Psychometrika. 29: 115–129.

Legendre, P.; Legendre, L. 1998. Numerical ecology. 2nd English ed. Vol. 20. 

Developments in Environmental Modeling. New York: Elsevier. 853 p.

Mather, P.M. 1976. Computational methods of multivariate analysis in physical 

geography. London: J. Wiley and Sons. 532 p.

McCune, B. 1997. Influence of noisy environmental data on canonical 

correspondence analysis. Ecology. 78: 2617–2623.

McCune, B. 2000. Lichen communities as indicators of forest health. The 

Bryologist. 103: 353–356.

McCune, B.; Dey, J.P.; Peck, J.E.; Cassell, D.; Heiman, K.; Will-Wolf, S.; 

Neitlich, P.N. 1997a. Repeatability of community data: species richness versus 

gradient scores in large-scale lichen studies. The Bryologist. 100: 40–46.

McCune, B.; Dey, J.P.; Peck, J.E.; Heiman, K.; Will-Wolf, S. 1994. Lichen 

communities. In: Lewis, T.E.; Conkling, B.L., eds. Forest health monitoring. 

Southeast loblolly/shortleaf pine demonstration interim report. Washington, DC: 

Office of Research and Development, U.S. Environmental Protection Agency: 

8.1–8.32. Chap. 8.

McCune, B.; Dey, J.P.; Peck, J.E.; Heiman, K.; Will-Wolf, S. 1997b. Regional 

gradients in lichen communities of the Southeast United States. The Bryologist. 

100: 145–158.

McCune, B.; Grace, J.B. 2002. Analysis of ecological communities. Gleneden 

Beach, OR: MjM Software Design. 300 p.

McCune, B.; Mefford, M.J. 2006. PC-ORD. Multivariate analysis of ecological 

data. Version 5 (updates for version 4 also available). Gleneden Beach, OR: MjM 

Software Design. http://home.centurytel.net/~mjm/index.html. (March 2010).
45

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

McCune, B.; Rogers, P.; Ruchty, A.; Ryan, B. 1998. Lichen communities for 

forest health monitoring in Colorado, USA. Report to the USDA Forest Service. 

Ogden, UT: U.S. Department of Agriculture, Forest Service, Rocky Mountain 

Research Station. 98 p.

Minchin, P.R. 1987. An evaluation of the relative robustness of techniques for 

ecological ordination. Vegetatio. 69: 89–107.

Muir, P.S.; McCune, B. 1988. Lichens, tree growth, and foliar symptoms of 

air pollution: Are the stories consistent? Journal of Environmental Quality. 

17: 361–370.

Nash, T.H., III, ed. 2008. Lichen biology. 2nd ed. Cambridge, United Kingdom: 

Cambridge University Press. 486 p.

National Acidic Deposition Program [NADP]. 2009. Ionic deposition maps. 

http://www.nadp.gov. (September 2009).

National Acid Precipitation Assessment Program [NAPAP]. 1991. Acidic 

deposition: state of science and technology reports. Volumes I–IV. Washington, 

DC: U.S. Government Printing Office. 

Neitlich, P.N.; Will-Wolf, S.; McCune, B. 2001. The lichen community indicator 

in the Forest Inventory and Analysis (FIA) and the Forest Health Monitoring 

(FHM) Programs. Development of regional gradient models: methods and 

priorities. Internal report. Washington, DC: U.S. Department of Agriculture, 

Forest Service. 15 p.

Nimis, P.L.; Scheidegger, C.; Wolseley, P., eds. 2002. Monitoring with lichens–

monitoring lichens. Introduction. NATO Science Series. The Hague, The 

Netherlands: Kluwer Academic Publishers. 424 p.

Økland, R.H. 1996. Are ordination and constrained ordination alternative or 

complementary strategies in general ecological studies? Journal of Vegetation 

Science. 7: 289–292.

Omernik, J.M. 1987. Ecoregions of the conterminous United States. Map. 

[1:7,500,000]. Annals of the Association of American Geographers. 

77(1): 118–125.

Patterson, P.; Will-Wolf, S.; Trest, M.T. 2009. FIA lichen indicator—analysis of 

QA data. In: Westphal, J., ed. Forest Inventory and Analysis national assessment 

of data quality for forest health indicators. Gen. Tech. Rep. NRS-53. St. Paul, 

MN: U.S. Department of Agriculture, Forest Service, Northern Research Station: 

Chap. 3.
46

GENERAL TECHNICAL REPORT PNW-GTR-807

Richardson, D.H.S. 1988. Understanding the pollution sensitivity of lichens. 

Botanical Journal of the Linnean Society. 96(1): 31–43.

Roberts, D.V. 2008. Statistical analysis of multidimensional fuzzy set ordinations. 

Ecology. 89: 1246–1260.

Sharnoff, S.; Rosentreter, R.R. 2008. Lichen use by wildlife in North America 

[and] lichens and invertebrates: a brief review and bibliography. http://www.

lichen.com/animals.html. (September 2008).

Sillett, S.C.; Neitlich, P.N. 1996. Emerging themes in epiphyte research in 

westside forests with special reference to cyanolichens. Northwest Science. 

70: 54–60.

Smith, C.; Geiser, L.; Gough, L.; McCune, B.; Ryan, B.; Showman, R. 1993. 

Species and communities. In: Stolte, K.; Mangis, D.; Doty, R.; Tonnessen, K.; 

Huckaby, L.S., eds. Lichens as bioindicators of air quality. Gen. Tech. Rep. RM224. Fort Collins, CO: U.S. Department of Agriculture, Forest Service, Rocky 

Mountain Research Station: 41–66. Chap. 4.

Søchting, U. 2004. Flavoparmelia caperata—a probable indicator of increased 

temperatures in Denmark. Graphis Scripta. 15(1/2): 53–56.

U.S. Department of Agriculture, Forest Service [USDA FS]. 2004. Forest 

Inventory and Analysis national core field guide. Volume 2: field data collection 

procedures for phase 3 plots. Version 2.0. Internal report. On file with: U.S. 

Department of Agriculture, Forest Service, Forest Inventory and Analysis, 

Rosslyn Plaza, 1620 North Kent Street, Arlington, VA 22209.

U.S. Department of Agriculture, Forest Service [USDA FS]. 2010a. Forest 

Inventory and Analysis. FIA DataMart. http://fiatools.fs.fed.us/fiadb-downloads/

fiadb3. (March 2010).

U.S. Department of Agriculture, Forest Service [USDA FS]. 2010b. Forest 

Inventory and Analysis national core field guide. FIA field methods for phase 3 

measurements. Version 3.0. On file with: U.S. Department of Agriculture, Forest 

Service, Forest Inventory and Analysis, Rosslyn Plaza, 1620 North Kent Street, 

Arlington, VA 22209. http://fia.fs.fed.us/library/field-guides-methods-proc/. 

(March 2010).

U.S. Environmental Protection Agency [US EPA]. 2008. Level III ecoregions. 

http://www.epa.gov/wed/pages/ecoregions/level_iii.htm. (June 2008).
47

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

van Dobben, H.F. 1993. Vegetation as a monitor for deposition of nitrogen and 

acidity. The Netherlands: Utrecht University. Privately published. 214 p. Ph.D. 

dissertation.

van Dobben, H.F.; ter Braak, C.J.F. 1998. Effect of atmospheric NH3

 on 

epiphytic lichens in the Netherlands: the pitfalls of biological monitoring. 

Atmospheric Environment. 32: 551–557.

Will-Wolf, S. 1988. Quantitative approaches to air quality studies. In: Nash, 

T.H., III; Wirth, V., eds. Lichens, bryophytes and air quality. Bibliotheca 

Lichenologica 30. Berlin-Stuttgart, Germany: J. Cramer in der Gebruder 

Borntraeger Verlagsbuchhandlung: 109–140.

Will-Wolf, S. 2007. Forest Inventory and Analysis Program lichen specialist 

training and QA procedures. http://fia.fs.fed.us/lichen/methods/. (September 

2009). On file with: Sarah Jovan, FIA lichen indicator advisor, USDA Forest 

Service, Forestry Sciences Laboratory, 620 SW Main, Suite 400, Portland, OR 

97331. 

Will-Wolf, S. 2009. Forest Inventory and Analysis Program lichen specialist 

identification procedures. http://fia.fs.fed.us/lichen/methods/. (September 2009). 

On file with: Sarah Jovan, FIA lichen indicator co-advisor (West), USDA Forest 

Service, Forestry Sciences Laboratory, 620 SW Main, Suite 400, Portland, OR 

97331. 

Will-Wolf, S. 2010. Analyzing lichen indicator data in the Forest Inventory 

and Analysis Program. Gen. Tech. Rep. PNW-GTR-818. Portland, OR: U.S. 

Department of Agriculture, Forest Service, Pacific Northwest Research Station. 

61 p.

Will-Wolf, S.; Ambrose, M.J.; Morin, R.S. [In press]. Relationship of a lichen 

species diversity indicator to environmental factors across Continental United 

States. In: Conkling, B.L.; Ambrose, M.J., eds. Forest Health Monitoring 2007 

national technical report. Gen. Tech. Rep. Asheville, NC: U.S. Department of 

Agriculture, Forest Service, Southern Research Station: Chap. 3.

Will-Wolf, S.; Geiser, L.H.; Neitlich, P.; Reis, A. 2006. Comparison of 

lichen community composition with environmental variables at regional and 

subregional geographic scales. Journal of Vegetation Science. 17: 171–184.
48

GENERAL TECHNICAL REPORT PNW-GTR-807

Will-Wolf, S.; Hawksworth, D.L.; McCune, B.; Sipman, H.J.M.; Rosentreter, 

R. 2004. Assessing the biodiversity of lichenized fungi. In: Mueller, G.M.; Bills, 

G.F.; Foster, M.S., eds. Biodiversity of fungi: standard methods for inventory and 

monitoring. San Diego, CA: Elsevier: 173–195. 

Will-Wolf, S.; Jovan, S. 2008. Lichens, ozone, and forest health—exploring 

cross-indicator analyses with FIA data. In: McWilliams, W.; Moisen, 

G.; Czaplewski, R., comps. 2008 Forest Inventory and Analysis (FIA) 

symposium. Proc. RMRS-P-56CD. Fort Collins, CO: U.S. Department of 

Agriculture, Forest Service, Rocky Mountain Research Station. One CD. 

http://www.treesearch.fs.fed.us/pubs/33326. (March 2010).

Will-Wolf, S.; Neitlich, P. 2007. Training manual: the lichen community indicator 

in FIA and FHM. USDA Forest Service, FIA Program. On file with: Sarah 

Jovan, FIA lichen indicator co-advisor (West), USDA Forest Service, Forestry 

Sciences Laboratory, 620 SW Main, Suite 400, Portland, OR 97231. 100 p.

Will-Wolf, S.; Neitlich, P.; Esseen, P.-A. 2002a. Monitoring biodiversity and 

ecosystem function: forests. In: Nimis, P.L.; Scheidegger, C.; Wolseley, P., 

eds. Monitoring with lichens—monitoring lichens. NATO Science Series. The 

Hague, The Netherlands: Kluwer Academic Publishers: 203–222.

Will-Wolf, S.; Scheidegger, C.; McCune, B. 2002b. Methods for monitoring 

biodiversity and ecosystem function. In: Nimis, P.L.; Scheidegger, C.; Wolseley, 

P., eds. Monitoring with lichens—monitoring lichens. NATO Science Series. The 

Hague, The Netherlands: Kluwer Academic Publishers: 147–162. 
49

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Appendix: Template for a Report Section Formally 

Documenting a Lichen Gradient Model and Its 

Implementation

This required part of a final Forest Inventory and Analysis (FIA) Program lichen 

gradient model report should be written for ecologists and analysts both inside and 

outside the FIA Program who will calculate lichen response indexes for new plots 

from data collected within an FIA lichen model region. All required instructions 

and pointers to needed data files and software must be included; details whose 

purpose is to support understanding of the model development process should not 

be included here.

[LICHEN MODEL REGION] LICHEN GRADIENT MODEL 

DOCUMENTATION AND SPECIFICATIONS FOR CALCULATING 

LICHEN RESPONSE INDEXES FOR NEW PLOTS 

1. BASIC MODEL DOCUMENTATION

Fill in the specific information for each topic listed in bold. Useful information and 

general instructions are written for each. Labels or values that change with each 

gradient are indicated between brackets.

Creation Date:

Updated Date:

Description of Changes:

Lichen Model Region Name:

Geographic Extent/Region Boundaries (verbal description): Include a map at 

the end of the document.

Geographic Extent/Region Boundaries (operational definition): If available, 

describe in detail how plots are assigned to this region.

FIA Administrative Regions included:

States included:

Number and brief description of all lichen response indexes 1- [X] defined 

for the region: List by order of importance, if order is defined. Index numbering is cumulative across all response indexes developed for a single lichen model 
50

GENERAL TECHNICAL REPORT PNW-GTR-807

region, whether from a single ordination model or from a hybrid model. Name 

lichen response index by major interpretation, i.e., pollution or air quality response 

index, climate/moisture or climate/temperature response index, etc. Name should 

clearly distinguish it from all other defined lichen response indexes in that region, 

i.e. lichen climate/elevation response index vs. lichen climate/latitude response 

index, etc.—list an undefined lichen community response gradient if necessary, 

for instance an undefined ordination axis required for model development. Indicate 

interpretation of index values, i.e., larger lichen response indexes indicate plots with 

more polluted air, or plots in cooler conditions, etc. Briefly describe the type, and 

number if needed, of model procedure(s) from which the lichen response index was 

developed, i.e., NMS ordination (lichen response index 1, 2, etc.); indicator species 

identification plus regression (lichen response index 1, 2, etc.); etc.

[named lichen response index 1]; [model type]; [interpretation of index values] 

[named lichen response index X]; [model type]: [interpretation of index values] 

Lichen gradient model developers: List names, institutions, phone numbers, and 

emails.

2. DETAILED MODEL DOCUMENTATION AND SPECIFICATIONS FOR 

DETERMINING LICHEN RESPONSE INDEXES FOR NEW PLOTS 

Include a separate subsection for each type of model used to develop lichen 

response indexes listed in Part 1. 

List the number and type of models described in this section. Label as subsection 

2A, 2B, 2C, etc.

2A. 1. DESCRIPTION OF MODEL [X]. Explanations for a nonmetric 

multidimensional scaling (NMS) ordination type of model for lichen response 

indexes 1, 2, etc. Describe another kind of ordination model similarly.

References for full description of model: Include internal and/or public Forest 

Service final report from model developers and location stored as well as publication of the model in a peer-reviewed journal.

Software and version used to develop the model:

Amount of original information in the final model, as compared with information in the model data set: This is often reported as a percentage.
51

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Method used by the software to calculate this amount of information in the 

final model: Most methods calculate proportion of variance of distances between 

pairs of plots that is represented on the final ordination axes. The calculation 

method is specific to the ordination technique and the software.

Brief description of lichen response indexes and length of gradients developed 

from this model: Copy lines from first section and add text—some of this is 

redundant, but redundancy is useful and important here. For gradient length, report 

the range of model scores for plots (length of axis for an NMS or other ordination 

model) based on the final reduced model data set. Do not report here the range of 

the lichen response indexes for plots in the model data set that were calculated using 

complete plot data. This second range based on lichen response indexes is likely to 

be different from the first range, as the model gradient is usually developed from a 

data set with some uncommon species excluded, and lichen response indexes for the 

same plot are calculated with all uncommon species included.

Number of plots included in the model: 

Criteria for choice and/or exclusion of plots: Briefly document all criteria: e.g., 

plots chosen to balance clean and dirty plots for 12 classes of environment; plots 

with tree basal area less than a minimum excluded and reason; plots with fewer 

than x species excluded; plots missing environmental data excluded; plots determined to be outliers as defined by developers or by software were excluded; etc.

Number of species included in the model:

Taxonomic usage year for model lichens data: No matter when lichen data 

were actually collected, the “taxonomic usage year” for the model is the year 

(= variable INVYR in FIA database tables) for which the REF_LICHEN_SPECIES 

and REF_LICHEN_SPP_COMMENTS tables best match the taxonomic choices 

made for the model data set.

Criteria for modification and/or exclusion of taxa: species at < x plots excluded, 

for others excluded state reason. List by at least acronym from the REF_LICHEN_

SPECIES table all taxa combined into groups for analysis, clearly indicate which 

acronym remains in the data set, how abundances for combined taxa are determined. All species names, numerical codes, and species acronyms MUST match 

a record in the FIA REF_LICHEN_SPECIES table.
52

GENERAL TECHNICAL REPORT PNW-GTR-807

Rules for determining abundances for combined taxa at one plot:

Abundances: individual taxa Abundance: final combined taxon

1 + 1 . . . . . . . . . . . . . . . . . . 2

1 + 1 + 1 + 1 + 1 . . . . . . . . . 2

More than five 1s . . . . . . . . 3

1 + 2 . . . . . . . . . . . . . . . . . . 2

2 + 2 . . . . . . . . . . . . . . . . . . 2

1 + 1 + 2 . . . . . . . . . . . . . . . 2

1 + 1 + 1 + 2 . . . . . . . . . . . . 3

1 + 2 + 2 . . . . . . . . . . . . . . . 3

3 + any others . . . . . . . . . . . 3

4 + any others . . . . . . . . . . . 4

Example list of combined taxa from internal FIA documentation (Jovan et al. 2008) 

for the Pacific Northwest West model (Geiser and Neitlich 2007):

For the following list, the “merged into” acronym is the one retained in the final 

model data set [Statements have been ordered here alphabetically by the “merged 

into” acronym.]:

Bryfri merged into Bry 

Vulcan merged into Cetcan 

Clafur, Clanor, Clapyx, Claumb merged into Cladonia. [This action combined four 

uncommon species that would probably have each been removed from the model 

data set, into a genus group large enough to be retained in the model data set.]”

Clacon merged into Claoch 

CLASUBSQ merged into Clasqu (the former acronym is found in a Region 6 

[National Forest] lichen species master list but not in the FIA REF_LICHEN_

SPECIES table; the final archived model data set has only acronyms in the FIA list)

Clabel merged into Clatra

Usndip (FIA), USNSST, and USNPAC merged into Usn (the two acronyms 

USNSST and USNPAC are found in a Region 6 [National Forest] lichen species 

masterlist but not in the FIA REF_LICHEN_SPECIES table; the final archived 

model data set has only acronyms in the FIA list)

Usncha merged into Usnfil 

Xanulo merged into Xanfal
53

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

List of taxa included in the ordination model: Model developers can name a 

separate archived file rather than include a list here. The list should include species 

code, species acronym, genus, and species name. In a separate archived file or a list 

if included in the documentation file (or both), all species names, numerical codes, 

and species acronyms must match a record in the appropriate FIA REF_LICHEN_

SPECIES table

Abundances relativized to row totals for the model? [YES, NO]

Distance measure used: For instance, Sørenson-abundance/Bray-Curtis distance. 

Gradient model lichen data set and format: Note—because computer software 

and file formats change frequently and archived files need to be usable for many 

(up to 50) years, choice of file format for archiving necessary data sets is very 

important. Only very general file formats should be used, and it may be advisable to submit each file in two formats, for instance the format specified by the 

recommended software to be used, and also as an ASCII text file or a file in xml 

format. It may be deemed necessary to submit hard copy versions of each file for 

emergency backup of electronic files. Apply a consistent and informative naming 

convention to all electronic files for a particular lichen gradient model, to facilitate 

identifying the needed files long after they have been archived.

File name and type: For example, [Lichen model region name] lichen model 

region model data set [filetype extension] is a plot-by-species table in [specified]format, with FIA P3_PLOT numbers and/or other codes [specify the 

code(s)] designating plots and species acronyms as in field SPP_ACRONYM 

of FIA database table REF_LICHEN_SPECIES. 

The format specified should be acceptable for software used to assign 

lichen response indexes to new plots, for instance module NMS Scores of PCORD [give version number] (McCune and Mefford 2006) or other software. If 

necessary, include a crosswalk table of public plot codes in the model data set to 

private FIA plot codes for internal FIA Program use. Give a brief rationale for 

the format submitted, and the software needed to use this file to assign lichen 

response indexes to new plots.

Where archived: Give the appropriate FIA office to contact. When the model 

data set is archived as a file linked to a database or on the Web, give appropriate location information. NOTE: all such files for currently completed models 

identify plots with FIA P3(FHM)HEXID number, which is considered private 
54

GENERAL TECHNICAL REPORT PNW-GTR-807

information. This file MUST be available to any FIA regional analyst calculating plot lichen indexes for archiving. Any outside user must first negotiate a 

confidentiality agreement with FIA, or a version of this file with public plot 

identifiers can be posted on the Web.

NMS model axis (gradient) scores:

File name and type: For example, [Lichen model region name-model-axisscores] [filetype extension], in [specified file format]. A file named xxx.gph in 

ASCII text format is the format currently required to assign lichen response 

indexes to new plots using module NMS Scores of PC-ORD v5.2 (McCune 

and Mefford 2006). Give rationale for the format submitted, and the software 

needed to use this file for assigning lichen response indexes to new plots. The 

codes used to identify plots in this file must match exactly the codes used to 

identify plots in the model data set.

Where archived: Give the appropriate FIA office to contact. When the axes 

data set is archived as a file linked to a database or on the Web, give appropriate location information. NOTE: all such files for currently completed models 

identify plots with P3(FHM)HEXID number, which is considered private information. This file MUST be available to any FIA regional analyst calculating 

plot lichen response indexes for archiving. Any outside user must first negotiate 

a confidentiality agreement with FIA, or a version of this file with public plot 

identifiers must be posted on the Web.

Model environmental data set: NOTE—this file should be archived for documentation of the lichen gradient model, but it should not be necessary for assignment of 

lichen response indexes to new plots with an ordination model. 

File name and type: For example, [Lichen model region name] lichen gradient 

model environmental data set [file type extension] is a plot x environmental 

factor table in [specify] format, with FIADB P3_PLOT numbers or other codes 

designating plots. One choice is the format acceptable for plotting environmental factor vectors on the final model graph using module GRAPH of PC-ORD 

[specify software version] (McCune and Mefford 2006). The codes used to 

identify plots in this file must match exactly the codes used to identify plots in 

the model data set. Give rationale for the format submitted, and the software 

needed to use this file to re-create the model diagram. See model references for 

description and origin of environmental factors.
55

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Where archived: Give the appropriate FIA office to contact. When the 

environmental variables data set is archived as a file linked to a database or 

on the Web, give appropriate location information. NOTE: currently all such 

files identify plots with P3(FHM)HEXID number, which is considered private 

information. This file may be useful to regional analysts and outside users 

analyzing and interpreting plot indexes. Any outside user must first negotiate 

a confidentiality memorandum of understanding with FIA, or a version of this 

file with public plot identifiers can be posted on the Web.

2A. 2. SPECIFICATIONS FOR CALCULATING LICHEN RESPONSE 

INDEXES FOR NEW PLOTS FROM MODEL [X]. Explanations are for specifications for calculating lichen response indexes from an NMS ordination model.

Software used to calculate lichen response indexes: At publication of this 

report, lichen response indexes are calculated for new plots using an NMS 

model by fitting plots to model axes with module NMS Scores in PC-ORD 

software v. 4.30 or later (McCune and Mefford 2006). A cautionary note about 

running NMS Scores–our experience is this routine can be very time consuming. It is convenient to allow the routine to run overnight. Test the run time for 

fitting 10 plots, then based on the time for this run, include only the number of 

plots at one time that are likely to be completed overnight. Alternately, run this 

routine on a computer that can be set aside and allowed to run for as long as 

necessary. For PC-ORD, each plot must have a unique identifier in a single field 

of maximum eight spaces. The private FIA variable P3(FHM)HEXID meets 

this criterion. Give directions and settings for any other software used in detail 

similar to that included below.

Preparation of data: All taxa should be combined exactly as for the model 

data set. First consult the relevant REF_LICHEN_SPP_COMMENTS table 

and follow all instructions for reconciling taxonomic changes between the 

taxonomic usage year for the model data set and the inventory year(s) for the 

plots to be fitted to the model. Next, follow any instructions given in model 

documentation above for further modification or combining taxa for the model. 

Assign lichen response indexes to new plots from the lichen gradient model 

using this modified data set.

Calculate lichen response indexes using all taxa, or calculate indexes 

using only taxa included in the model [follow recommendation for the 

model]. Each model should have one of these recommendations stated. 
56

GENERAL TECHNICAL REPORT PNW-GTR-807

The default recommendation should be to use all taxa at each plot. If the 

recommendation is to use only taxa included in the model, documentation 

should include a brief justification.

Species abundances relativized by plot totals for the model?: [YES, NO] 

For assigning response indexes to new plots, follow the procedure adopted 

for the model.

Selecting setup options for fitting plots to model [X] using module 

NMS Scores in PC-ORD:

• Select [the distance measure used for the model analysis] distance 

measure.

• Select Simultaneous method for fitting. This choice is most often the 

appropriate choice for assigning lichen indexes to a new FIA plot for 

program purposes, to better maintain the status of indexes being independent of one another. In some cases, model developers may recommend and justify the choice to fit axes one at a time.

• Select user-defined unspecified fit (stress) criterion, then enter this fit 

criterion. The unspecified fit (stress) criterion for fitting to this model 

is determined by model developers. A reasonable criterion might be 

two standard deviations higher than average unspecified fit (stress) calculated from assigning lichen response indexes for the model to more 

than 100 plots in a single NMS Scores run.

• Set NMS Scores extrapolation limit to ± 10 percent.

• Be sure to DESELECT the NMS Scores option “List Intermediate 

Results.”

Saving the output from an NMS Scores run: For saving output from some 

other software, give instructions in equivalent detail.

Save run specifications. Save run specifications as text in a file that can be 

archived as metadata for the assigned plot lichen response indexes. Consult 

the FIA Program information management unit before submitting files for 

archiving, to decide appropriate file format. Add calculation date, software 

name, and version to this file. AXIS SUMMARY STATISTICS in NMS 

Scores output refer only to that particular set of plots and are not of general 

interest, so they should not be archived.

Save plot lichen response indexes and flags for poor fit. The assigned 

axis scores become the plot lichen response indexes. Information in the 
57

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

output table BEST-FIT SCORES ON EACH AXIS from all runs of NMS 

Scores for the population of plots in one inventory year should be saved 

in a separate file to be archived, and when data structures are in place, to 

be uploaded to the FIA database. Consult the FIA Program information 

management unit before submitting files for archiving, to decide appropriate file format. For simultaneous fit, only one column for Fit and Flag is 

output for all the axes on which scores are assigned. For archiving in the 

FIA database, duplicate the Fit and Flag columns to accompany values for 

each defined lichen response index (= axis), and add a column for inventory 

year. 

It is convenient to maintain in this file both a single-field plot identifier 

the analyst might choose to use and the STATECD, COUNTYCD, NIMS_

PLOT and FIADB_PLOT (FIA database field names) identifier fields in a 

master file, to facilitate relinking the results of analysis with the standard 

FIA plot identifiers for upload to the FIA database and crosswalking with 

other data. Neither NIMS_PLOT nor FIADB_PLOT is by itself a unique 

plot identifier, so state and county codes must always be linked with one of 

these plot codes.

Field names and form of file to archive: These field names should be included in 

each record for each plot. FIA P3(FHM)HEXID (or other unique plot identifier(s) 

if used), STATECD, COUNTYCD, NIMS_PLOT, FIADB_PLOT, Inventory year, 

Lichen gradient model region code, Response_index_1, Response_index_1_fit, 

Response_index_1_flag, Response_index_2_index, Response_index_2_fit, 

Response_index_2_flag, etc. 

The default value for the Response index X flag field is “0.” All plots having a 

value higher than the user-specified maximum acceptable stress-related fit criterion 

are assigned a flag value of “1” for general poor fit, entered in the “flag” field. If 

the lichen response index for a plot is more than 10 percent higher than the maximum for a plot in the model, as compared with full model gradient (= axis) length 

(see glossary), a flag value of “2” is entered; if it is more than 10 percent lower than 

the minimum for a plot in the model as compared with full gradient length, a flag 

value of “3” is entered in the “flag” field for that plot. If a plot qualifies for two fit 

flags, the higher flag value is entered. 

Each plot has one row or record. The file to be submitted should be in a format 

agreed upon after consultation with FIA Program information managers. After an 

official FIA database structure has been developed for such files, this file will need 

to conform to the designated upload format. 
58

GENERAL TECHNICAL REPORT PNW-GTR-807

Estimating fit and assigning poor-fit flags with alternate software: Identify 

plots with lichen response index values beyond the acceptable range for each 

gradient using criteria analogous to those detailed above. Estimate stress-related or 

unspecified fit of plots to the model; identify those plots whose fit value is beyond 

the criterion. Assign the same three flag values as above to indicate the different 

kinds of poor fit.

2B. 1. DESCRIPTION OF MODEL [Y] Explanations for indicator species/

regression model type, for lichen response indexes 2, 3, etc.

Model description: Indicate briefly the exact kind of model developed, for instance 

indicator species identification based on abundance of species, followed by linear 

regression.

References for full description of model: Include reference to the internal Forest 

Service final report from model developers (and location of archived copy) as well 

as publication of the model in a peer-reviewed journal.

Software and version(s) used to develop the model: Include name, date, version 

number for all software, such as “the Indicator Species Analysis module in PCORD v5.2 (McCune and Mefford 2006)” and the name and version of a statistics 

package used to develop a regression model, etc.

Brief description of lichen response indexes and length of gradients developed 

from this model: Give number and brief description of lichen response indexes 

developed from this model. Copy lines from the first section—this is redundant, but 

redundancy is useful here. For each response index, give the range of model lichen 

scores (gradient length) calculated from the final reduced model data set, not the 

range of lichen response indexes calculated for plots included in the model using 

full plot data.

Number of plots included:

Criteria for choice and/or exclusion of plots: Write “Same as for model X,” 

“Subset of plots in model X within a certain elevation range,” etc. or include full 

documentation of criteria. 

Number of species included:
59

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Taxonomic usage year for model lichens data: The “taxonomic usage year” is the 

year for which notes in an FIA REF_LICHEN_SPP_COMMENTS table match the 

taxonomic choices made for the model data set, no matter when data were actually 

collected.

Criteria for modification and/or exclusion of taxa: Write “Same as for model 

X,” “Subset of species in model X” with additional criteria [specify], etc., or give 

full documentation of criteria. It is very desirable for these criteria to be the same 

for all response indexes developed for a hybrid model, to reduce confusion and to 

simplify the task of modifying data before calculation of all lichen response indexes 

for plots. All species names, numerical codes, and species acronyms MUST match 

a record in the FIA REF_LICHEN_SPECIES table.

List of [pollution] indicator species: Note: it is not necessary that all of these species are included in the model data set. The list should include species code, species 

acronym, genus, and species name for each species.

Source for list of [pollution] indicator taxa: For instance, list software name 

and version used to determine indicator species, list published literature sources, 

or both. Note which indicator taxa were shown to be significant indicators 

quantitatively.

List of all taxa, including indicator species, included in the model data set: List 

in documentation or provide the name of an archived file that lists the taxa. The list 

and/or file should include species code, species acronym, genus, and species name. 

All species names, numerical codes, and species acronyms must match a record in 

the NIMS_REF_LICHEN_SPECIES table.

Abundances relativized to row totals for the model? [YES, NO]

Model equation: Includes all parameters needed to define the model, such as 

regression coefficients, standard deviations, etc. and precise definition including 

units for all terms in a model equation, or other needed inputs.

List plot environmental data, if any, and units required as inputs to this model: 

Include elevation in meters, latitude in decimal degrees, etc. Give FIA database 

field name and the table where this information resides, if appropriate. Geographic 

information such as elevation, latitude, or longitude related to exact plot locations 

should be used.
60

GENERAL TECHNICAL REPORT PNW-GTR-807

2B.2. SPECIFICATIONS FOR DETERMINING LICHEN RESPONSE 

INDEXES FOR NEW PLOTS WITH MODEL Y. Explanations of specifications 

for determining lichen response indexes from indicator species/regression model 

type.

Software used to calculate lichen response indexes: Indicate any specific 

software needed. For a common kind of analysis such as linear regression that 

is available in many software packages, indicate particular form of the model, 

choices, specification of error terms, etc. that should be used when calculating 

response indexes for new plots from the model.

Preparation of data: All taxa should be combined exactly as for the model. 

Consult the most recent FIA REF_LICHEN_SPP_COMMENTS table; follow 

all instructions for reconciling taxonomic changes between the taxonomic 

usage year for the model data set and the inventory year(s) for the new plots to 

be assigned lichen response indexes from the model.

Calculate lichen response indexes using all taxa, or calculate using only 

taxa included in the model [follow recommendation for the model]. Each 

model should have the recommendation stated. The default recommendation should be to use all taxa, modified according to instructions. If the 

recommendation is to use only taxa included in the model, documentation 

should include a brief justification.

Abundances relativized by plot totals for the model? [YES, NO] For 

assigning lichen indexes to new plots, follow the procedure adopted for the 

model.

Required environmental data. Obtain for each plot all environmental 

data, if any, required as inputs to this model (elevation, latitude, etc.). All 

values should be in the units used for the model. Obtain data for the model 

from the FIA database field name and table where this information resides, 

if appropriate. Information related to exact plot locations should be used.

Calculating lichen response indexes for new plots with the model: Insert a 

detailed step-by-step protocol for calculating lichen response indexes for new 

plots with the model. Emphasize units for values wherever needed. If appropriate, append a spreadsheet file with formulas embedded.
61

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

Estimating the degree of fit of a plot to the model: This explanation is a suggested method for estimating the degree of fit to an indicator species/regression 

type model. Other methods chosen should be documented in similar detail.

The relative abundance of species at a plot that are NOT included in the 

regression model data set (AbNOT) is an estimate of the degree of fit of that 

plot to the model. This value can be calculated as: 

AbNOT = 

100 × 

∑ abundance scores for all species at the plot NOT in the model data set

∑ abundance scores for all species at the plot

Assigning a flag for poor fit: The default value for this field is “0.” If for this 

model, a maximum acceptable value of AbNOT for adequate fit to the model 

has been recorded, then all plots having a value higher than this maximum 

acceptable value should be assigned a value of “1” for general poor fit recorded 

in the “flag” field. If the lichen response index for a plot is more than 10 percent higher than the maximum for a plot in the model as compared with the full 

range of lichen gradient scores for plots in the original model (gradient length), 

enter a flag value of “2” for too high; if it is more than 10 percent lower than 

the minimum for a plot in the model as compared with gradient length in the 

original model, enter a flag value of “3” for too low. If a plot qualifies for two 

fit flags, enter the higher flag value. 

If no threshold value of AbNOT for adequate fit (or no other general fit 

criterion) has been determined, assign flags for poor fit based only on criteria 

related to gradient length.

Saving the lichen response indexes: For archiving and uploading into the 

FIA database, add lichen response index, fit value, and fit flag to an existing file that includes other response indexes for the same plots in the same 

region. Recall that lichen response indexes for a single region are numbered 

consecutively even if they are derived with different techniques. For instance, 

lichen response indexes 1 and 2 might have come from an ordination model as 

described in section A of documentation, and lichen response index 3 has come 

from a regression model as described in section B of documentation. All three 

response indexes with their fit values and fit flags would be archived in the 

same file.

Field names and form of file to archive: Use the field names and format 

as described in section “Saving the output,” p. 56, for an ordination gradient 

model X. 
62

GENERAL TECHNICAL REPORT PNW-GTR-807

Glossary

application phase—The second phase in implementation of the lichen indicator 

in the Forest Inventory and Analysis (FIA) Program. In this phase, lichen data are 

collected from standard FIA plots and are assigned indexes for response to environmental factors defined from an existing regional “FIA lichen gradient model.” 

Many analyses based on these indexes are possible in this phase.

calibration phase—The first phase of implementation of the lichen indicator in the 

FIA Program. In this phase, lichen data are collected from standard FIA plots and 

limited analyses are conducted. Supplemental data are collected and an “FIA lichen 

gradient model”for a particular “FIA lichen model region” is developed. This phase 

lasts for a region until a gradient model is completed and implemented.

FHM—The Forest Health Monitoring Program in the U.S Forest Service. Current 

focus is nationwide reports and special projects.

FIA—The Forest Inventory and Analysis Program in the U.S. Forest Service. 

Responsible for nationwide collection of data describing the Nation’s forests, from 

plots located on a permanent grid.

FIA lichen gradient model—A quantitative model developed for a particular “FIA 

lichen model region.” Such a model relates lichen community composition in the 

region to major environmental factors of interest, such as climate and air quality. 

Abbreviated after first use under a topic as “lichen gradient model” or just “model” 

in appropriate context and when usage is unambiguous.

FIA lichen indicator—Collectively refers to lichen species and communities 

as used in the standard national FIA Program as biomonitors of forest health, 

including all the different indexes developed from the lichen data collected, and 

more loosely to all aspects of lichens in the program. Often abbreviated in text as 

“lichen indicator.”

FIA lichen indicator advisor—An FIA employee or contractor whose responsibility it is to advise and collaborate with the FIA Program on all aspects of 

the FIA lichen indicator from data collection through analysis and reporting to 

long-term planning. This person has a Ph.D. or equivalent professional training 

in both lichenology and community ecology. A major responsibility is to oversee 

maintenance of high standards related to lichenology expertise in all aspects of the 

program. Often abbreviated in text as “indicator advisor” or “lichen IA.”
63

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

FIA lichen model region—A geographic area defined for the FIA Program for 

which a quantitative “FIA lichen gradient model” will link lichen community 

composition to major environmental factors of interest including climate and air 

quality. After such a gradient model has been developed, FIA plots within this 

region have response indexes calculated from lichen data for all environmental factors defined from the model for this region. Abbreviated after first use under a topic 

as “lichen model region” or just “region” in appropriate context and when usage is 

unambiguous.

gradient—The term as used in plant ecology and for this document is defined as 

a range of values ordered from small to large (or vice versa) for either an environmental factor or a derived number representing some biotic factor of interest. For 

instance, an ordination axis score is a derived number that for a lichen gradient 

model represents an aspect of lichen species composition at one plot. The entire 

ordered range of scores for one ordination axis thus represents a species composition gradient. And the full ordered range of average summer high temperatures for 

a set of plots represents a gradient for summer temperature.

gradient length—The largest minus the smallest score on a published lichen model 

gradient. Specifically, gradient length is the full range of scores on one ordination 

axis or full range of values from a regression or other model, for plots in the original model based on the archived modified model data set. Gradient length is not 

the range of lichen response indexes for the same set of plots as calculated based on 

full plot data. Gradient length is reported in documentation for the lichen gradient 

model and does not change over time. Differences between resampled plots for 

lichen response indexes are converted to percentage of the published model gradient 

length for comparison with data quality objectives.

inventory—Defined for the FIA Program as a one-time assessment across a 

geographic region of the patterns of forest vegetation and its components based on 

field sampling of permanent plots.

inventory year—Defined for the FIA Program as the calendar year in which the 

majority of plots in that group were sampled: FIA database variable INVYR. A 

group of plots is identified each calendar year for sampling following standard 

program protocol, with most sampled during the spring through fall seasons of 

that year. 
64

GENERAL TECHNICAL REPORT PNW-GTR-807

lichen air quality index—A unitless numerical index that represents the position of 

an FIA plot at one particular time on an air quality gradient defined from a quantitative “FIA lichen gradient model” for a particular “FIA lichen model region.” This 

index is one of a class of derived “lichen response indexes.” The index is calculated 

based on lichen species composition from one sample of a plot. The air quality 

gradient as defined is statistically independent of all other lichen response gradients 

defined for that particular lichen gradient model. The air quality gradient may be 

related to measured or modeled quantities of particular pollutants, but the index is 

not defined in terms of those measured units.

lichen climate index—A unitless numerical index that represents the position of 

an FIA plot at one particular time on a climate gradient defined from a quantitative “FIA lichen gradient model” for a particular “FIA lichen model region.” This 

index is one of a class of derived “lichen response indexes.” The index is calculated 

based on lichen species composition from one sample of a plot. The lichen climate 

gradient as defined is statistically independent of all other lichen response gradients 

defined for that particular lichen gradient model. The lichen climate gradient may 

be related to measured or modeled values of particular climate components such as 

average summer temperature or average annual rainfall, but the index is not defined 

in terms of those measured units.

lichen index—A unitless numerical value based on FIA lichen data that represents 

some defined relative performance or response of lichens at an FIA plot. Both 

primary indexes and response indexes derived from a lichen gradient model are 

included under this phrase.

lichen response index—A unitless numerical value that represents the status of 

an FIA plot at one particular time with respect to one of the gradients defined 

for lichen response to an environmental factor from a quantitative lichen gradient model. The derived index is an indicator of lichen community response to 

the environmental factor. A response index is calculated based on lichen species 

composition from one sample of a plot. Each lichen response gradient defined for a 

model is statistically independent of all other lichen response gradients defined for 

that particular model. 

lichen species richness index—The count of macrolichen species recorded from 

an FIA plot surveyed following standard protocols. This primary index is available 

for any FIA plot surveyed for lichens; it is an indicator of relative (not absolute, 

see Will-Wolf 2010) diversity of the lichen community at that plot at that time. 

Often abbreviated after first use under a topic as “richness index” when usage is 

unambiguous.
65

Development of Lichen Response Indexes Using a Regional Gradient Modeling Approach for Large-Scale Monitoring of Forests

measurement year—The year in which a plot was actually visited in the field for 

sampling: FIA database variable MEASYR. Most plots in an “inventory year” are 

visited for lichen sampling in the spring through fall seasons of the inventory year, 

but visits in the winter of the following year for lichen sampling are possible. In the 

latter case, measurement year and inventory year would differ.

monitoring—Defined for the FIA Program as repeated assessment over multiple 

years and across a geographic region of the patterns and trends of forest vegetation 

and its components from repeated field sampling of the same population of permanent plots.

multivariate analysis—A kind of data analysis, here of ecological communities, 

that identifies in a single analysis the most important patterns of many species at 

many plots, and often also the related patterns of many environmental variables 

measured at many plots.

ordination—A particular kind of multivariate data analysis whose goal for ecological studies is to graphically portray the most important patterns among many 

sampling units (often plots) and many species (or other attributes of the sampling 

units), and the relations of those patterns to many environmental variables (either 

primary or secondary to ordination analysis), in very few (usually two to four) 

dimensions, or independent axes.

ordination joint plot—By convention, a graph that displays field plots with proximity representing similarity (or “distance”) based on the lichen species data, with 

superimposed lines that represent the direction and strength of correlation between 

environmental factors and the displayed axes of the graph (Legendre and Legendre 

1998, McCune and Grace 2002). Depending on the ordination technique used, the 

superimposed lines may represent information external or internal to the ordination 

analysis itself.

Pacific Northwest Research Station

Web site http://www.fs.fed.us/pnw

Telephone (503) 808-2592

Publication requests (503) 808-2138

FAX (503) 808-2130

E-mail pnw_pnwpubs@fs.fed.us

Mailing address Publications Distribution

Pacific Northwest Research Station

P.O. Box 3890

Portland, OR 97208-3890
U.S. Department of Agriculture

Pacific Northwest Research Station

333 SW First Avenue

P.O. Box 3890

Portland, OR 97208-3890

Official Business

Penalty for Private Use, $300