Introduction
Illegal hunting is a major threat to many large-bodied animals, and understanding the factors that promote or deter poaching is fundamental to effective management (Rowcliffe et al., Reference Rowcliffe, de Merode and Cowlishaw2004; Shivik, Reference Shivik2006). However, detailed data are scarce. Many poaching reports are collected incidentally, and may be indirect (Madhusudan & Karanth, Reference Madhusudan and Karanth2002). Such reports may be difficult to analyse because of lack of a sampling design and uncertain locations (Reutter et al., Reference Reutter, Helfer, Hirzel and Vogel2003). Nonetheless, if data are abundant and their limitations are taken into account, heterogeneous reports can offer important insights, particularly in species that are low-density, shy, or otherwise difficult to study (Milner-Gulland et al., Reference Milner-Gulland, Kholodova, Bekenov, Bukreeva, Grachev and Amgalan2001).
The Andean bear Tremarctos ornatus is a prime example of a large-bodied mammal threatened by poaching throughout its distribution (Servheen et al., Reference Servheen, Herrero and Peyton1998). Along with habitat loss, poaching has caused the species to be categorized as Vulnerable on the IUCN Red List (IUCN, 2007) and Endangered in Venezuela (Rodríguez & Rojas-Suárez, Reference Rodríguez and Rojas-Suárez1999). Although bear hunting is illegal in Venezuela (República de Venezuela, 1992, 1996), an estimated 2.5 bears per year have been hunted for the past 80 years across the Cordillera de Mérida region of the Andes alone (Yerena, Reference Yerena1994).
In Venezuela poaching has been attributed to a variety of socio-economic and cultural causes but few data have been available to determine their relative importance and thereby propose effective management. One possible solution is to examine spatial patterns of risk, because different causes are expected to create different patterns. For example, if subsistence hunting, religious practices, rites of passage, or commercial traffic are the most important causes, then poaching risk will be highest where bears are concentrated, in areas of high-quality bear habitat (the habitat hypothesis). In this case, increasing the extent and effectiveness of protected areas may be most effective (Herrera et al., Reference Herrera, Nassar, Michelangeli, Rodríguez and Torres1994; Goldstein et al., Reference Goldstein, Paisley, Wallace, Jorgenson, Cuesta and Castellanos2006). Alternatively, if poaching risk is highest in areas where human populations are highest, one of several encounter hypotheses may be more valid: poaching may be due to opportunistic hunting (the encounter-kill hypothesis) or conflict over cattle predation and crop raiding (the encounter-conflict hypothesis; Goldstein, Reference Goldstein1991; Vineyard & Torres, Reference Vineyard and Torres2004). Targeted environmental education may be the best approach in the first case but electric pasture fences, loss compensation schemes, or changes in agricultural practices would be more appropriate for the second (Goldstein et al., Reference Goldstein, Paisley, Wallace, Jorgenson, Cuesta and Castellanos2006; Shivik, Reference Shivik2006). Finally, factors such as protection status may also modulate poaching risk (Yerena, Reference Yerena1994), producing a spatial pattern similar to the other encounter hypotheses (the encounter-protection hypothesis) but implying management actions similar to those appropriate to the habitat hypothesis.
Given the difficulty of studying poaching directly (Goldstein et al., Reference Goldstein, Paisley, Wallace, Jorgenson, Cuesta and Castellanos2006), our purpose here is to present a framework to employ data often excluded from conventional analyses to understand spatial patterns of threat to a threatened species. Our specific aim was to test hypotheses about the factors driving T. ornatus poaching in the Venezuelan Andes by fitting logistic regression models to examine the distribution of poaching probability in the Cordillera de Mérida.
Study area
Tremarctos ornatus is found only in the Andes, which in Venezuela comprise the Cordillera de Mérida and the Sierra de Perijá (Fig. 1). Given the scarcity of bear observations in the latter (13 in total; Viloria et al., Reference Viloria, Mondolfi, Yerena and Herrera1995; EY unpubl. data), we focused only on the former. This 26,702 km2 area spans altitudes of 500–5,000 m, and is covered by a complex mosaic of cloud and dry forests, semi-arid formations, and natural páramo savannah. Because no relative abundance models yet exist for T. ornatus in Venezuela, we had to assume that bears were uniformly distributed in this area. Protected areas presently cover 50% of the region, and include 11 national parks and one natural monument, which are thought to contain high-quality bear habitat (Yerena, Reference Yerena, Servheen, Herrero and Peyton1998). Although population centres, roads, and agricultural development are concentrated at lower altitudes, human population density is > 250 people per km2 and is increasing rapidly (Vivas, Reference Vivas1992).
Fig. 1 Distribution of reports of Andean bear in the Cordillera de Mérida, Venezuela: (a) Political boundaries and bear distribution, with rectangle indicating study area, (b) density of all reports per 20 km2, and (c) density of poaching reports per 20 km2.
Methods
Presence data
We compiled 844 bear presence reports from the literature (Herrera et al., Reference Herrera, Nassar, Michelangeli, Rodríguez and Torres1994; Yerena, Reference Yerena1994; Goldstein, Reference Goldstein2002; Viloria et al., Reference Viloria, Mondolfi, Yerena and Herrera1995, and references therein), from the unpublished field and interview data of SG-R, and from interviews with D. Torres (Fundación AndígenA, Mérida) and H. Zambrano (EcoVida, Táchira). All statistical analyses were conducted using R (R Development Core Team, 2005). We first examined the effect of geopolitical location (state) on the density of poaching reports (resolution 20 km2), using a χ2 test and an analysis of standardized residuals, and then mapped poaching report density using the kriging & simulation module in Idrisi (Eastman, Reference Eastman2001).
We then selected the 287 direct observations with detailed location information and assigned UTM coordinates using global positioning system field data, or by locating place names on 1 : 100,000 maps (ICNSB, 1969). Reports were assigned additional attributes with Idrisi, using various source maps (rescaled to a 500 m raster resolution; Table 1). Vegetation index (veg) and distance to the nearest protected area (dprot) were assigned by record year but this was impossible for the distance to the nearest road (road), or for population density (hpop), because no temporal data exist. This could diminish the apparent influence of road and hpop on poaching, because, for example, bear reports from 1964 had to be associated with present human population densities. Nevertheless, this bias was minor for population densities because the relative sizes of population centres have not changed greatly (Vivas, Reference Vivas1992). The potential bias associated with roads is less clear, because although officially they have not expanded, many formerly dirt roads have been paved recently in some areas, and new dirt roads have been added (SG-R, unpubl. data).
Table 1 Geographical and anthropogenic characteristics for georeferenced reports of T. ornatus in Venezuela.
In addition to these explanatory variables, we assigned each record a response variable, poach, which identified it as being an event of illegal hunting or not. In some cases both record types shared identical geographical coordinates. We therefore assigned a weight to each location, calculated as the number of poaching reports divided by the total, and used these in model fitting (see below; R Development Core Team, 2005).
Data transformation and evaluation of spatial autocorrelation
Three of our explanatory variables had highly skewed distributions, which can destabilize regression coefficient estimates (Quinn & Keough, Reference Quinn and Keough2002). We therefore rescaled them by applying a square-root transformation to road and dprot, and a log transformation to hpop. We also searched for possible collinearity among explanatory variables using Fox & Monette's (Reference Fox and Monette1992) variance inflation factor (VIF) but found it only in latitude (x) and longitude (y; VIF = 12.11 and 12.59, respectively). Because eliminating one of these would have reduced geographical resolution, we retained both in initial analyses.
Poaching reports may be spatially autocorrelated because the probability of additional events may increase in areas where hunting has been previously reported (Segurado & Araújo, Reference Segurado and Araújo2004). Ignoring autocorrelation can lead to poor model selection, highly biased factor weightings, and misleading predictions (Carroll & Pearson, Reference Carroll and Pearson2000; Quinn & Keough, Reference Quinn and Keough2002). Therefore, we next examined the raw data for autocorrelation using spatial correlograms, which plot the covariance between each pair of observations against the distance separating them (Legendre, Reference Legendre1993). Because poach was weighted and binomial rather than continuous and normal, we used the modifications suggested by Bjørnstad & Falck (Reference Bjørnstad and Falck2001) to fit a covariance function with splines, and to generate confidence intervals with 100 bootstrap replicates (implemented with the spline.correlog function of the R ecolsupl package).
Autocorrelation can be generated by several causes: in a given area there may be more presence reports because more research effort was invested, or because it contains more bears or more hunters. We therefore generated separate correlograms for all data, and non-poaching and poaching reports, to distinguish these causes as far as possible. If autocorrelation is found one possible approach is to discard observations until independence is achieved. However, a more data-efficient approach, employed below, is to fit statistical models that directly correct for autocorrelation (Legendre, Reference Legendre1993).
Model 1 Assuming spatial independence
The first set of models presumed spatial independence among observations with a binomial error distribution (McCullagh & Nelder, Reference McCullagh and Nelder1983). Logistic regression captured the categorical response of the dependent variable (poach) in the full model as a function of five continuous independent variables: altitude (z), vegetation index (veg), transformed human population density (hpop), distance to the nearest road (road) and to the nearest protected area border (dprot). We did not include interaction terms because no plot of the residuals from a model without interactions versus the product of any pairs of variables had a slope different from zero (Quinn & Keough, Reference Quinn and Keough2002). Thus, the full Model 1 was:
![$$\eqalignno{{\rm logit}(poach) &= \log \left( {{{w(poach)} \over {1 - w(poach)}}} \right) \cr&= \beta _0 + \beta _1 z + \beta _2 veg \cr&\quad+ \beta _3 \log [hpop] + \beta _4 \sqrt {road} \cr&\quad+ \beta _5 \sqrt {prot}$$](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160921020601289-0415:S0030605308006996:S0030605308006996_equ1.gif?pub-status=live)
in which w is the weighting factor described above and βx is the relevant regression coefficient.
After fitting the full Model 1 we extracted the minimum subset of explanatory variables that maximized predictive power to arrive at a final Model 1, using mixed stepwise selection, which selects the smallest model among those with the minimum Akaike Information Criterion (AIC; R Development Core Team, 2005). We evaluated the importance of each variable using Wald's test (the standardized regression coefficient divided by its standard error; Quinn & Keough, Reference Quinn and Keough2002). Finally, we evaluated the autocorrelation in the residuals of the final Model 1 using ecolsupl as described above.
Model 2 Assuming spatial autocorrelation
Spatial autocorrelation may be incorporated in regression models in several ways, including incorporating geographical coordinates as covariates in classic auto-regressive models (Legendre, Reference Legendre1993) or using geostatistical methods to define a spatial covariance function (Hoeting et al., Reference Hoeting, Davis, Merton and Thompson2006). Because both approaches yielded similar results (data not shown), we present the simpler classic method here. The full Model 2 thus included both latitude (x) and longitude (y) in addition to the variables described for the full Model 1 above. Although we found Spearman correlations > 0.28 between x and all other variables, plotting their products versus residuals indicated no reason to include interaction terms (as above), leaving the full Model 2 as:
![$$\eqalignno{{\rm logit}(poach) &= \log \left( {{{w(poach)} \over {1 - w(poach)}}} \right) \cr&= \beta _0 + \beta _1 z + \beta _2 veg \cr&\quad+ \beta _3 \log [hpop] + \beta _4 \sqrt {road} \cr&\quad+ \beta _5 \sqrt {prot} + \beta _6 x + \beta _7 y$$](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160921020601289-0415:S0030605308006996:S0030605308006996_equ2.gif?pub-status=live)
in which all variables are as described in Model 1. The methods described above were then used to select the final Model 2, to evaluate the importance of selected variables, and to construct a spatial correlogram of residuals.
Model performance and risk map construction
We compared the performance of both final models by evaluating model fit and prediction capacity (Segurado & Araújo, Reference Segurado and Araújo2004). For model fit, we used deviance in a χ2 test (Quinn & Keough, Reference Quinn and Keough2002). For predictive capacity, we examined both sensitivity and precision (Liu et al., Reference Liu, Berry, Dawson and Pearson2005). Sensitivity measured the proportion of observed poaching events that were correctly assigned, whereas precision measured the proportion of predicted poaching events that truly were poaching events. To estimate both, observed reports must be classified as correctly or incorrectly predicted. Because logistic models produce probabilities, a threshold is needed for classification. Because the optimal threshold is expected to vary with the prevalence of positives (Segurado & Araújo, Reference Segurado and Araújo2004), we fixed an objective threshold that minimized erroneous classifications (Liu et al., Reference Liu, Berry, Dawson and Pearson2005). Lastly, we selected the models with the best overall performance, and used their predictions to generate poaching risk maps at a resolution of 2.5 km2 (predict.glm function; R Development Core Team, 2005).
Results
Temporal and geographical patterns in Andean bear reports
We found a total of 844 bear reports for 1940–2004. Although 31% had no date, the temporal distribution of the remainder indicated a constant poaching prevalence for the first 50 years, after which it declined greatly, in spite of a simultaneous increase in sampling effort (Fig. 2). Nearly all dated poaching reports associated with cattle predation or crop raiding occurred in the 1990s, and comprised just 0.8% of the total. Of the total, 557 reports had precise locations, of which 21% were poaching reports; the rest were direct sightings (57%), physical remains of natural deaths (1%), and secondhand reports of sightings (21%, excluded from further analyses due to unreliability). We also excluded all reports below 500 m, for the same reason. This left 287 georeferenced reports for subsequent analyses, in which poaching frequency was 29% (84/287), and the odds of poaching was 0.414.
Fig. 2 Temporal distribution of accumulated reports of Andean bear presence and poaching reports in the Cordillera de Mérida, Venezuela.
Táchira State had the lowest density of total reports, and Lara and central Mérida had the highest (Fig. 1a,b). Poaching reports were distributed throughout the study area; however, they were proportionally more dense in Mérida than in any other State (χ2 test, P = 0.0081, standardized residual = 1.989; Fig. 1c).
Model selection and spatial autocorrelation
The final Model 1 had only altitude (z) as a significant effect (Table 2; Fig. 3a), whereas the final Model 2 also included longitude (x) and the square root distance to the closest protected area border (dprot). In Model 2, the highest risk occurred below 1,300 m (Table 2; Fig. 3b), and between 71° 25′ W 8° 22′ N and 72° 39′ W 7° 39′ N (Fig. 3c). dprot had a smaller but still significant effect on poaching, with lowest risk inside protected areas, or within c. 2.5 km beyond park borders (Fig. 3d).
Fig. 3 Predictions of logistic regression models. Plots of predicted poaching versus significant explanatory variables in Model 1 (a, spatial independence), and Model 2 (b-d, spatial autocorrelation). See text for details of models. Points above and below the classification threshold (poaching odds) indicate high and low poaching, respectively.
Table 2 Regression coefficients (±SE) estimated for both final models (assuming spatial independence and autocorrelation), with standardized regression coefficient (SRC) indicating the relative importance of each variable in explaining poaching probability, and the Akaike Information Criterion (AIC).
*P ≤ 0.05; ***P ≤ 0.001; ****P ≤ 0.0001 (Student's t test)
Deviance analysis indicated that Model 2 had a significantly better fit than Model 1 (P << 0.0001), highlighting the importance of spatial autocorrelation. However, even Model 2 explained only 12% of the variance in the data (Table 2) indicating, perhaps unsurprisingly, that poaching probability was also influenced by factors not included in this model. Correlograms confirmed that reports were autocorrelated, and were spatially independent only when pairs of reports were considered at distances > 40.3 km (Fig. 4a1). Autocorrelation was due mainly to factors operating in non-poaching reports: poaching reports considered alone had no significant autocorrelation (Figs 4a2,3). Model 1 residuals had a covariance and amplitude similar to the raw data (Fig. 4b), whereas Model 2 reduced both to an important degree (Fig. 4c), suggesting that it better accounted for spatial structure.
Fig. 4 Correlograms of (a1-3) raw presence data, (b) Model 1 residuals (presuming spatial independence), and (c) Model 2 residuals (including autocorrelation), with 95% confidence intervals. Raw data are divided into: (a1) all reports, (a2) non-poaching reports, and (a3) poaching reports. a = amplitude in km, the minimum record distance for correlation to be zero; cov = maximum spatial covariance.
Model performance and risk maps
Although neither model explained a high proportion of the variance, Model 2 was superior, and had predictions with higher precision (0.64 vs 0.49, respectively) and slightly higher sensitivity (0.44 vs 0.36). A poaching probability threshold of 26% minimized erroneous classifications of reports, and was thus used in estimates of sensitivity and precision; however, Model 2 outperformed Model 1 regardless of where this threshold was set.
Finally, Model 2 highlighted spatial patterns that were not evident in Model 1 (Fig. 5). Although both models predicted higher poaching risk at the margins of the study area (Fig. 5a,b), this was much more marked in Model 2, which also predicted concentrated poaching risk in the south. Similarly, although both models predicted lower risk in the centre of the study area, it was only in Model 2 that this lower risk was clearly associated with protected areas. Errors were low across much of the study area for both models (Fig. 5c,d), and tended to increase at the limits of the species distribution, although more so in the case of Model 2.
Fig. 5 Predicted poaching risk to Andean bears in the Cordillera de Mérida according to (a) Model 1 and (b) Model 2, with associated standard errors, (c) and (d), respectively. The rectangle of the figure is that in Fig. 1a.
Discussion
Poaching patterns and the importance of spatial autocorrelation
Within the Cordillera de Mérida in Venezuela it is clear that poaching has been a long-standing and widespread threat to Andean bears: poaching reports spanned 66 years and were distributed throughout the study area. However, reports as a whole were highly spatially autocorrelated, being most likely to occur within 40 km of each other (Fig. 4a1). This does not appear to be caused by autocorrelation in hunting behaviour, because poaching reports considered alone had a covariance of essentially zero (Fig. 4b3). Instead, factors operating in non-poaching reports, which had a similar covariance and amplitude as the overall dataset (Fig. 4b2), may have been more important. At 44.2 km, the average study area radius of the five researchers responsible for the majority of reports was similar to the amplitude of the autocorrelation (Fig. 4a1), suggesting that research effort was strongly biased. However, ecological factors determining T. ornatus presence may have also played a role, for example, if remaining habitat patches or bear home ranges have a radius of c. 40 km. Because these factors remain unknown, further research on T. ornatus habitat use will be vital to distinguish causes of autocorrelation.
Whatever its origin, spatial autocorrelation was clearly important to take into account. Although values for both were low, Model 2 was superior to Model 1 in overall fit, in the proportion of variation explained, and in the sensitivity and precision of its predictions. Also, Model 2 revealed more factors as significant than did Model 1 (Table 2). These results support the conclusion that models ignoring spatial autocorrelation tend to underestimate the significance of important variables and to include only those with the highest autocorrelation (Carroll & Pearson, Reference Carroll and Pearson2000).
The habitat and encounter-protection hypotheses
Our results also provide direction for future research to clarify the factors influencing Andean bear poaching in the Cordillera de Mérida. The habitat hypothesis posits that poaching risk should be highest simply where bears are most frequent. However, lower poaching risk was significantly associated with protected areas (Table 2, Fig. 5b). These areas are presumed to encompass most of the highest-quality Andean bear habitat in Venezuela (Yerena, Reference Yerena1994; Servheen et al., Reference Servheen, Herrero and Peyton1998). Thus, in the last 20 years (from which the majority of our data come; Fig. 2), poaching does not appear to have been highest where bears have been presumably most frequent, suggesting that hunters may not specifically seek out bears, and undermining the habitat hypothesis. Clearly, a detailed habitat model for this species in Venezuela is needed to refute this hypothesis further.
Instead, poaching risk significantly increased at low- and medium altitudes (< 1,300 m; Table 2, Fig. 3b), where all human activities, including agricultural and urban expansion, tend to overlap with bear habitat (Vivas, Reference Vivas1992; Yerena, Reference Yerena1994). This suggests that poaching events may be more related to the probability of encounter with humans, supporting the encounter hypotheses. One possibility is that the significant relationship between protected status and lowered poaching is causal (the encounter-protection hypothesis). Supporting this possibility is the fact that seven of today's 12 protected areas were established between 1988 and 1993 (56% of the 10,500 km2 presently protected), when poaching prevalence decreased greatly (Fig. 2).
However, that protected status could have a causal impact on poaching is surprising, given that structural adjustment policies adopted in the midst of these park designations subsequently crippled the institutions charged with their management. Between 1988 and 1993 both the Venezuelan Ministry of Environment and the National Parks Institute had budgets and staff cut by half (Markandya et al., Reference Markandya, Yero, Cariola, Lacabana, Velasco, Caraballo and Reed1996). Because this institutional dismantling has continued, our results may be among those supporting the controversial notion that even so-called paper parks can play a positive role in conservation (Naughton-Treves et al., Reference Naughton-Treves, Holland and Brandon2005). However, other events in this period make simple interpretations impossible. For example, the 1992 Penal Code of the Environment changed bear poaching from a fined infraction to a felony requiring jail time (República de Venezuela, 1992), and may be partly responsible for reduced poaching (although data from other countries suggest that increased penalties may simply increase incentives to avoid detection rather than actually decreasing poaching; Rowcliffe et al., Reference Rowcliffe, de Merode and Cowlishaw2004).
The encounter-conflict and encounter-kill hypotheses
An alternative explanation for increased risk where humans are encountered is that agriculture/habitat overlaps cause human-bear conflicts that directly encourage poaching events (Nielsen et al., Reference Nielsen, Herrero, Boyce, Mace, Benn, Gibeau and Jevons2004; Goldstein et al., Reference Goldstein, Paisley, Wallace, Jorgenson, Cuesta and Castellanos2006). Given the presence of inhabited areas near park borders it has been feared that increased human-bear contact in these areas could in fact increase the risk of poaching (Yerena, Reference Yerena1994; Nielsen et al., Reference Nielsen, Herrero, Boyce, Mace, Benn, Gibeau and Jevons2004). However, this does not seem to be the case: parks had a significant and positive effect, which extends a few kilometres beyond their borders (Fig. 3d). Furthermore, we found no evidence that such conflicts were a widespread factor in poaching: < 1% of poaching reports were associated with cattle predation or crop-raiding events, and all occurred in a short time span. Nevertheless, it may be that such conflict is a general motivation for hunting rather than a direct trigger, and future research on hunting behaviour should focus on this point.
A final explanation for increased poaching risk at low and medium altitudes is that opportunities are simply increased for encounters between bears and hunters in pursuit of other animals (the encounter-kill hypothesis), because altitude was also moderately correlated with human population density and roads. If most poaching is opportunistic rather than conflict-driven, then management interventions aimed at minimizing conflicts with humans and bears (Gunther et al., Reference Gunther, Haroldson, Frey, Cain, Copeland and Schwartz2004; Vineyard & Torres, Reference Vineyard and Torres2004; Goldstein et al., Reference Goldstein, Paisley, Wallace, Jorgenson, Cuesta and Castellanos2006) will be less useful at reducing poaching than other actions, such as educational campaigns to reduce negative perceptions and increased law enforcement (Herrera et al., Reference Herrera, Nassar, Michelangeli, Rodríguez and Torres1994).
To distinguish definitively among the various encounter hypotheses, information about poacher attitudes, behaviours and needs are clearly required. Such data could be obtained in interviews using modifications of public-health survey techniques developed for sensitive topics (Lensvelt-Mulders et al., Reference Lensvelt-Mulders, Hox, van der Heuden and Maas2005). Interviews could also provide valuable data about the relative importance of increased penalties and expanded protected areas.
High-risk areas
Our final important result is that the areas of most intense poaching risk, in the south and along the borders of the species distribution in the Cordillera de Mérida (Fig. 5b), were different from the areas with the highest number of poaching reports (Fig. 1c). This underscores the importance of employing unified multivariate analyses on properly-filtered data, because future research should clearly focus on areas of true risk rather than simply on areas that have been researched extensively in the past. As such, all of the following areas are of interest for future research because they are areas of possible high risk but low sampling effort: the southern slope and border of the Sierra Nevada National Park, the north slope and border of La Culata National Park, and Tapo-Caparo, Páramos de Batallón y la Negra and El Tamá National Parks (Fig. 5c).
Our model may predict increased risk in these southern and border regions for several reasons. Firstly, increased poaching may be precisely why the edge of the species distribution is presently where it is, with risk so high beyond this edge that the species cannot persist there (Nielsen et al., Reference Nielsen, Herrero, Boyce, Mace, Benn, Gibeau and Jevons2004). Secondly, the south is more amenable to agriculture, ranching, and urban development because it has gentler slopes and easier access. Finally, three of the four protected areas in this region were among the last to be declared, and may have little historical tradition to bolster their efficacy.
Although tempting, it would be premature to conclude that these predicted high-risk areas also represent high priority areas for conservation action. This is partly because conservation priority-setting should incorporate societal and cultural values in addition to risk assessments (e.g. Mace & Lande, Reference Mace and Lande1991). More importantly, however, even our best regression model correctly predicted only 64% of poaching reports, and explained only 12% of total variation in risk: effects were significant but small nonetheless. Clearly, major factors influencing poaching risk, such as temporal patterns in human population and road construction, could not be included in our models. For this reason, further data on species presence and abundance, current habitat availability, and poaching behaviour are clearly needed to distinguish among poaching hypotheses and for the proposal of conservation actions.
Conclusions
In spite of these limitations, we conclude that analysing and visualizing the spatial distribution of reports and of inferred risk are useful for examining hypotheses about effective interventions for conserving threatened species. Using highly dispersed and heterogeneous presence-only data for a species that is difficult to study directly, we showed that spatial analysis can effectively account for autocorrelation, and found convincing evidence that protected areas have a positive effect on a known threat to that species. This bolsters the case for initiatives currently underway to create new protected areas and strengthen existing ones (Yerena et al., Reference Yerena, Padrón, Vera, Martínez and Bigio2003). We also found scant evidence that human-bear conflicts are direct triggers of poaching, suggesting that environmental education initiatives aimed at influencing hunters who accidentally encounter bears may be more important than conflict-reduction approaches (AndígenA, Reference AndígenA2006). We expect that as additional reports and more factors are added, including information on bear distribution and poacher behaviour, the resulting risk maps will become increasingly useful as conservation planning tools. We also expect that similar initiatives for other elusive, large-bodied mammals of this area (such as Panthera onca and Tapirus terrestris) will produce useful perspectives on poorly-understood threats, direct future research, and provide a baseline against which to judge future actions.
Acknowledgements
Funds for this study were provided by the Instituto Venezolano de Investigaciones Científicas (IVIC), the Iniciativa de Especies Amenazadas of Provita, the Venezuelan Fondo Nacional de Ciencia, Tecnología, e Investigación, Idea Wild, the International Bear Association, the Scott Neotropical Fund of the Cleveland Metroparks Zoo, the Rufford Small Grants for Nature Conservation, the Denver Zoological Foundation, the Pittsburgh Zoo and PGG Aquarium, Fundación Gran Mariscal de Ayacucho, the National Parks Institute, and Fundación para la Defensa de la Naturaleza (FUDENA). We are grateful for assistance in the field from Alfredo Freitez Sánchez, Francisco Daza Russo, Eugenio Gusman Pérez, Henrry Sánchez, Marcial Quiroga, Anna Veit, Dorangel Nuñez, and the Pacheco, Mateus-Leal, Freitez-Sánchez, Montilla and Saavedra families, as well as from Frecopal, the Junta Ambientalista de Ospino, the Guardiantes de Terepaima, the Frente Ecológico de Cubiro, the Parque Zoológico y Botánico Bararida, and the Guardaparques Universitarios. We also thank the National Parks Institute for guidance during field work, R. Lazo for help constructing our geographical information system, G. Barreto, D. Chivers and D. Augeri for sage advice, the EcoSIG of the Centro de Ecología of IVIC, and the Laboratorio de Manejo y Conservación de Fauna Silvestre of Universidad Simón Bolívar for use of their facilities and software, Fundación AndígenA and EcoVida for their presence data, I. Herrera, G. Canale, and two anonymous reviewers for useful comments on the manuscript, and FUDENA and Fundación AndígenA for their support.
Biographical sketches
Ada Sánchez-Mercado is interested in habitat modelling and population viability analysis and, with EY, SG-R and others, belongs to the organizing committee for updating the Andean Bear Conservation Action Plan in Venezuela. José R. Ferrer-Paris works with distributional database management and modelling biodiversity patterns. He is currently developing NeoMapas, a project generating distribution and abundance maps for a variety of Neotropical taxa. Edgard Yerena researches planning and policy for protected areas and biodiversity conservation. Shaenandhoa García-Rangel is interested in modelling Andean bear habitat- and landscape-use with GIS databases, and in wildlife management and ursid ecology. She collaborates with FUDENA and Fundación AndígenA in Andean bear conservation. Kathryn M. Rodríguez-Clark has broad interests in statistical, modelling and genetic approaches to conservation problems. She also collaborates with Venezuelan NGO Provita in environmental education.
