Hostname: page-component-848d4c4894-pjpqr Total loading time: 0 Render date: 2024-06-16T14:37:41.642Z Has data issue: false hasContentIssue false

Factors affecting poaching risk to Vulnerable Andean bears Tremarctos ornatus in the Cordillera de Mérida, Venezuela: space, parks and people

Published online by Cambridge University Press:  10 July 2008

Ada Sánchez-Mercado*
Affiliation:
Laboratorio de Ecología y Genética de Poblaciones, Centro de Ecología, Instituto Venezolano de Investigaciones Científicas, Apartado 20632, Caracas 1020-A, Venezuela.
José R. Ferrer-Paris
Affiliation:
Laboratorio de Biología de Organismos, Centro de Ecología, IVIC, Caracas, Venezuela.
Edgard Yerena
Affiliation:
Departamento de Estudios Ambientales, Universidad Simón Bolívar, Caracas, Venezuela.
Shaenandhoa García-Rangel
Affiliation:
Wildlife Research Group, The Anatomy School, University of Cambridge, Cambridge, CB2 3DY, UK.
Kathryn M. Rodríguez-Clark*
Affiliation:
Laboratorio de Ecología y Genética de Poblaciones, Centro de Ecología, Instituto Venezolano de Investigaciones Científicas, Apartado 20632, Caracas 1020-A, Venezuela.
*
*Laboratorio de Ecología y Genética de Poblaciones, Centro de Ecología, Instituto Venezolano de Investigaciones Científicas, Apartado 20632, Caracas 1020-A, Venezuela. E-mail kmrc@ivic.ve or kmrodriguezclark@gmail.com
Rights & Permissions [Opens in a new window]

Abstract

Worldwide, many large mammals are threatened by poaching. However, understanding the causes of poaching is difficult when both hunter and hunted are elusive. One alternative is to apply regression models to opportunistically-collected data but doing so without accounting for inherent biases may result in misleading conclusions. To demonstrate a straightforward method to account for such biases, and to guide further research on an elusive Vulnerable species, we visualized spatio-temporal poaching patterns in 844 Andean bear Tremarctos ornatus presence reports from the Cordillera de Mérida, Venezuela. To create maps of poaching risk we fitted two logistic regression models to a subset of 287 precisely georeferenced reports, one ignoring and one including spatial autocorrelation. Whereas the variance explained by both models was low, the second had better fit and predictive ability, and indicated that protected status had a significant positive effect on reducing poaching risk. Poaching risk increased at lower altitudes, where all indicators of human disturbance increased, although there was scant evidence that human-bear conflicts are a major direct trigger of poaching events. Because highest-risk areas were different from areas with most bear reports, we speculate that hunting may be driven by opportunistic encounters, rather than by purposeful searches in high-quality bear habitat. Further research comparing risk maps with bear abundance models and data on poaching behaviour will be invaluable for clarifying poaching causes and for identifying management strategies.

Type
Papers
Copyright
Copyright © Fauna & Flora International 2008

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:

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:

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.

References

AndígenA, (2006) Proyecto Oso Andino. Http://www.andigena.org/proyecto_oso_andino/ [accessed 24 August 2006].Google Scholar
Bjørnstad, O.N. & Falck, W. (2001) Nonparametric spatial covariance functions: estimation and testing. Environmental and Ecological Statistics, 8, 5370.CrossRefGoogle Scholar
Carroll, S.S. & Pearson, D.L. (2000) Detecting and modeling spatial and temporal dependence in conservation biology. Conservation Biology, 14, 18931897.CrossRefGoogle ScholarPubMed
Dobson, J.E., Bright, E.A., Coleman, P.R., Durfee, R.C. & Worley, B.A. (2000) LandScan: a global population database for estimating populations at risk. Photogrammetric Engineering & Remote Sensing, 66, 849857.Google Scholar
Eastman, J.R. (2001) Idrisi 32. Guide to GIS and Image Processing. Clark University, Worcester, USA.Google Scholar
Fox, J. & Monette, G. (1992) Generalized collinearity diagnostics. Journal of the American Statistical Association, 87, 176183.CrossRefGoogle Scholar
Global Forest Watch (2002) Venezuelan Roads, 1st edition. US National Imagery and Mapping Agency, Washington, DC, USA.Google Scholar
Goldstein, I. (1991) Spectacled bear predation and feeding behavior on livestock in Venezuela. Studies on Neotropical Fauna and Environment, 26, 231235.CrossRefGoogle Scholar
Goldstein, I. (2002) Andean bear-cattle interactions and tree nest use in Bolivia and Venezuela. Ursus, 13, 369372.Google Scholar
Goldstein, I., Paisley, S., Wallace, R., Jorgenson, J.P., Cuesta, F. & Castellanos, A. (2006) Andean bear–livestock conflicts: a review. Ursus, 17, 815.CrossRefGoogle Scholar
Gunther, K.A., Haroldson, M.A., Frey, K., Cain, S.L., Copeland, J. & Schwartz, C.C. (2004) Grizzly bear–human conflicts in the Greater Yellowstone ecosystem, 1992–2000. Ursus, 15, 1022.2.0.CO;2>CrossRefGoogle Scholar
Herrera, A.-M., Nassar, J.M., Michelangeli, F., Rodríguez, J.P. & Torres, D. (1994) The spectacled bear in the Sierra Nevada National Park of Venezuela. International Conference of Bear Research and Management, 9, 149156.Google Scholar
Hoeting, J.A., Davis, R.A., Merton, A.A. & Thompson, S.E. (2006) Model selection for geostatistical models. Ecological Applications, 16, 8798.CrossRefGoogle ScholarPubMed
ICNSB (1969) Proyección Provisional de Sur América (PSAD_1956) UTM Zona 19N. (Sheets 5737-40, 5837-41, 5939-42, 6040-3, 6142-5, 6244-6, 6345-6, 5543-4, and 5643–8), Ministerio de Obras Públicas, Dirección de Cartografía Nacional, Caracas, Venezuela.Google Scholar
IUCN (2007) 2007 IUCN Red List of Threatened Species. IUCN, Gland, Switzerland. Http://www.iucnredlist.org [accessed 29 January 2008].Google Scholar
Legendre, P. (1993) Spatial autocorrelation: trouble or new paradigm? Ecology, 74, 16591673.CrossRefGoogle Scholar
Lensvelt-Mulders, G.J.L.M., Hox, J.J., van der Heuden, P.G.M. & Maas, C.J.M. (2005) Meta-analysis of randomized response research: thirty-five years of validation. Sociological Methods & Research, 33, 319348.CrossRefGoogle Scholar
Liu, C., Berry, P.M., Dawson, T.P. & Pearson, R.G. (2005) Selecting thresholds of occurrence in the prediction of species distributions. Ecography, 28, 385393.CrossRefGoogle Scholar
Mace, G.M. & Lande, R. (1991) Assessing extinction threats: toward a reevaluation of IUCN threatened species categories. Conservation Biology, 5, 148157.CrossRefGoogle Scholar
Madhusudan, M.D. & Karanth, K.U. (2002) Local hunting and the conservation of large mammals in India. Ambio, 31, 4954.CrossRefGoogle ScholarPubMed
Markandya, A., Yero, L., Cariola, C., Lacabana, M., Velasco, F.J., Caraballo, A. et al. (1996) Case study for Venezuela. In Structural Adjustment, the Environment, and Sustainable Development (ed. Reed, D.), pp. 201224. Earthscan Publications, London, UK.Google Scholar
McCullagh, P. & Nelder, J.A. (1983) Generalized Linear Models. Chapman and Hall, London, UK.CrossRefGoogle Scholar
Milner-Gulland, E.J., Kholodova, M.V., Bekenov, A., Bukreeva, O.M., Grachev, Iu.A., Amgalan, L. et al. (2001) Dramatic declines in saiga antelope populations. Oryx, 35, 340345.CrossRefGoogle Scholar
Naughton-Treves, L., Holland, M.B. & Brandon, K. (2005) The role of protected areas in conserving biodiversity and sustaining local livelihoods. Annual Review of Environment and Resources, 30, 219252.CrossRefGoogle Scholar
Nielsen, S.E., Herrero, S., Boyce, M.S., Mace, R.D., Benn, B., Gibeau, M.L. & Jevons, S. (2004) Modelling the spatial distribution of human-caused grizzly bear mortalities in the Central Rockies ecosystem of Canada. Biological Conservation, 120, 101113.CrossRefGoogle Scholar
Quinn, G.P. & Keough, M.J. (2002) Experimental Design and Data Analysis for Biologists. Cambridge University Press, Cambridge, UK.CrossRefGoogle Scholar
R development core team (2005) R. R Foundation for Statistical Computing, Vienna, Austria.Google Scholar
República de Venezuela (1992) Ley Penal del Ambiente de la República de Venezuela. Gaceta Oficial Nº 4.358, January 03. Caracas, Venezuela.Google Scholar
República de Venezuela (1996) Decreto 1485: Animales Vedados para la Caza de las Especies Incluidas o No en la Lista de Animales de Caza. Gaceta Oficial Nº 36.059, October 07. Caracas, Venezuela.Google Scholar
Reutter, B.A., Helfer, V., Hirzel, A.H. & Vogel, P. (2003) Modelling habitat-suitability using museum collections: an example with three sympatric Apodemus species from the Alps. Journal of Biogeography, 30, 581590.CrossRefGoogle Scholar
Rodríguez, J.P., Lazo, R., Solórzano, L.A. & Rojas-Suárez, F. (2005) Cartografía digital básica de las áreas naturales protegidas de Venezuela, CIET, IVIC, CI, UNESCO and MARN, Caracas, Venezuela.Google Scholar
Rodríguez, J.P. & Rojas-Suárez, F. (1999) El Libro Rojo de la Fauna Venezolana, 2nd edition. Provita-Fundación Polar, Caracas, Venezuela.Google Scholar
Rowcliffe, J.M., de Merode, E. & Cowlishaw, G. (2004) Do wildlife laws work? Species protection and the application of a prey choice model to poaching decisions. Proceedings of the Royal Society of London B: Biological Sciences, 27, 26312636.CrossRefGoogle Scholar
Segurado, P. & Araújo, M.B. (2004) An evaluation of methods for modelling species distributions. Journal of Biogeography, 31, 15551568.CrossRefGoogle Scholar
Servheen, C., Herrero, S. & Peyton, B. (eds) (1998) Bears: Status Survey and Conservation Action Plan. IUCN/Species Survival Commission Bear and Polar Bear Specialist Groups, Gland, Switzerland.Google Scholar
Shivik, J.A. (2006) Tools for the edge: what's new for conserving carnivores. BioScience, 56, 253259.CrossRefGoogle Scholar
Tucker, C.J., Pinzon, J.E. & Brown, M.E. (2004) Satellite Drift-corrected and NOAA-16 Incorporated Normalized Difference Vegetation Index (NDVI), Monthly 1981. Global Inventory Modelling and Mapping Studies, University of Maryland, College Park, USA.Google Scholar
USGS (2006) The GTOPO30 Global Digital Elevation Model. Http://edc.usgs.gov/products/elevation/gtopo30/gtopo30.html [accessed 5 September 2005].Google Scholar
Viloria, A.L., Mondolfi, E., Yerena, E. & Herrera, F. (1995) Nuevos registros de osos de anteojos o frontino (Tremarctos ornatus, F. Cuvier) en la Sierra de Perijá. Memoria de Fundacion La Salle de Ciencias Naturales, 55, 313.Google Scholar
Vineyard, T. & Torres, D.A. (2004) Solar electric fencing to prevent Andean bear-human conflicts in the Venezuelan Andes. International Bear News, 13, 26.Google Scholar
Vivas, L. (1992) Los Andes Venezolanos. Academia Nacional de la Historia, Caracas, Venezuela.Google Scholar
Yerena, E. (1994) Parques Nacionales y Conservación Ambiental: Corredores Ecológicos en los Andes de Venezuela, Fundación Polar, Caracas, Venezuela.Google Scholar
Yerena, E. (1998) Spectacled bear conservation action plan (Tremarctos ornatus): status and management of the spectacled bear in Venezuela. In Bears: Status Survey and Conservation Action Plan (eds Servheen, C., Herrero, S. & Peyton, B.), pp. 309. IUCN/Species Survival Commission Bear and Polar Bear Specialist Groups, Gland, Switzerland.Google Scholar
Yerena, E., Padrón, J., Vera, R., Martínez, Z. & Bigio, D. (2003) Building consensus on biological corridors in the Venezuelan Andes. Mountain Research and Development, 23, 215218.CrossRefGoogle Scholar
Figure 0

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.

Figure 1

Table 1 Geographical and anthropogenic characteristics for georeferenced reports of T. ornatus in Venezuela.

Figure 2

Fig. 2 Temporal distribution of accumulated reports of Andean bear presence and poaching reports in the Cordillera de Mérida, Venezuela.

Figure 3

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.

Figure 4

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).

Figure 5

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.

Figure 6

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.