Global Analysis of Protected Area Placement

All datasets are global in scale, in raster (grid) format, and projected into Albers Equal Area projection at a one km² resolution. We used ArcGIS 9.1 to harmonize projections, cell size, and extent and used Python 2.4 in order to remove all marine areas and to create individual text files for each variable for each country. All further analyses were done in R 2.7.1 [40] .
The analytical framework we used was a general linear model (package “glm” in R 2.7.1) with a probit link because the regressions we run here involve binary outcomes. In the first set of regressions (for Table 1a), we are explaining whether or not a location is in the nation's protected area network. In the second set of regressions (for Table 1b), we are examining only the locations that are in the network and explaining whether or not a location is in a protected area that is accorded higher status. One important point is that the coefficients for a given variable not only need not be the same in those regressions and could even be different in sign (for instance, protection may be biased towards steep slopes but, within the protected network, higher status could be biased towards flat areas). Thus Table 1a and 1B do not have similar results by construction. These probabilities were generated using the “predict” command, along with the coefficients from the original regression models, in the “stats” package in R 2.7.1.
Information on PA location came from the 2007 World Database on Protected Areas (WDPA) [41] . Only PAs classified by the International Union for Conservation of Nature (IUCN) in categories 1 through 6, and only countries with PA networks of 100 km² or more, were included. When two PAs overlapped, we assigned that area the highest IUCN classification of the two. Due to high potential error rates [3] , PAs without polygon boundaries (i.e., point representation only) were not included.
For comparisons across different management categories, we included only countries with 100 km² or more of categories I – II and 100 km² or more of categories III – VI. To analyze protection over time (Figure 1), we used information included in the WDPA on date of PA creation. When analyzing over time, PAs with no date were included in each temporal step, distributing the error uniformly over the analysis.
We obtained elevation data from the Shuttle Radar Topography Mission (SRTM) [18] . The source for this data set was the Global Land Cover Facility (www.landcover.org). The SRTM gathered elevation data on a near-global scale, generating a very complete high-resolution elevation database. We calculated slope values from the SRTM elevation dataset. All slope values are degrees from horizontal. Distance to roads was calculated from a vector road network extracted from the VMAP Level0 dataset [19] . While the quality of this data is variable it is the only freely available global road dataset to characterize the global road network. Distance to urban areas was calculated using the Gridded Rural Urban Population dataset (GRUMP) [20] , which provides a gridded and global extent of urban population. Agricultural suitability was taken from a dataset provided by the International Institute for Applied Systems Analysis [21] . The dataset (plate 28) incorporates climate, soil type, land cover, and slope of terrain to measure agricultural suitability, ranking each grid cell from 0 (no constraints) to 9 (severe constraints). We used the World Wide Fund for Nature (WWF) Ecoregions product to determine ecoregion type [23] . The Ecoregions product delineates 8 biogeographic realms, 14 biomes, and 867 ecoregions. Included in the WWF Ecoregions project are data on terrestrial vertebrate species richness. These data encompass more than 26,000 terrestrial vertebrates (amphibians, reptiles, birds, and mammals) and were gathered from literature, expert opinion, and online datasets. Richness is at the resolution of ecoregion.