Elephant carcass encounter data

African elephant carcass data were collected as part of the Monitoring the Illegal Killing of Elephants (MIKE) programme, which was instituted by the Convention on International Trade in Endangered Species (CITES) in 2002 and since then has worked with wildlife authorities across Africa to implement a ranger-based monitoring programme. The programme collates annual carcass counts from 53 sites (mostly protected areas, but often extending into neighbouring unprotected zones) in 29 countries across sub-saharan Africa. Full details of the monitoring methods are described elsewhere12, but, in essence, rangers on regular patrols record the location of any elephant carcass encountered and identify whether death was the result of natural mortality, management or illegal killing (almost always poaching for ivory, but very occasionally the result of retaliation in human-elephant conflict). Between 2002 and 2017, the programme has recorded 18,007 carcasses in Africa, of which 8860 were identified as illegal killings, providing 607 observations from 53 sites in 16 years (includes all records received by February 2018). Several sites did not report carcasses every year, or joined the programme later than 2002.

It should be noted that these carcass encounter data, collated by the MIKE programme, show a few potential limitations12: (a) variation in background mortality (i.e. carcasses resulting from natural mortality or management) is unknown, but influences PIKE as it is assumed to be constant across years and sites. Background mortality (here natural mortality) is increased during droughts and periods of low rainfall49,50, so we aimed to account for variable natural mortality by estimating the effect of site-specific annual precipitation on PIKE and setting this effect to zero for the predictions of the model. (b) Calculating PIKE across sites and years is based on the assumption that detection probabilities are the same for all carcasses, resulting from illegal activities, management or natural reasons. This might be an unlikely assumption, as the data are collected by anti-poaching patrols with the objective to deter illegal activities. It seems plausible to assume, however, that this bias is rather constant across space and time leading to an accurate estimation of trends and association with covariates. (c) Based on data from 53 sites across Africa, the prediction of poaching rates might not cover the full uncertainty of continental estimates, yet the surveyed area covers 25% of the area where African savannah elephants are extant residents51 and about 50% of Africa’s savannah elephant population6,15.

Covariates

The choice of covariates (Supplementary Table 1), considered as potential drivers of poaching intensity, was guided by previous studies7,12 and expert knowledge52. We included covariates that we considered to relate to demand or supply of elephant ivory, including factors that vary at temporal and spatial levels and two separate indicators of poverty: infant mortality rate and poverty density. Poverty is a complex, multidimensional problem that cannot easily be measured in a single variable53, but the negative impact of poverty on illegal wildlife activities has been highlighted before18 so it is important to consider multiple aspects of poverty. Not all covariates were available at the highest site-by-year resolution. Below, we present them in the following order: site-by-year, country-by-year, site-level, annual. Before the analysis, all covariates were centred and standardised to have a mean of 0 and a standard deviation of 1. We also tested for collinearity among predictors. All combinations showed Spearman’s ρ2 estimates <0.5, which we considered a non-problematic correlation (see Supplementary Fig. 1).

Infant mortality rate: The infant mortality rate (IMR) measures the number of deaths of children under one year of age per 1000 live births and is a crude indicator of development and socio-economic status levels in a community54. Note that IMR is included entirely as a proxy for one axis of poverty55: if IMR is strongly predictive of elephant poaching rates we would not interpret this as suggesting that healthcare interventions alone would be expected to impact poaching rates.

IMR estimates were available at site level for the year 2000, produced by the Centre for International Earth Science Information Network (CIESIN56). Annual IMR estimates by country were made available by the United Nations (UN) inter-agency group for child mortality estimation57. As both spatial and temporal variability are high, we combined the two datasets to obtain IMR estimates for every site in every year. In reality, improvement rates in rural and urban areas may differ, but national changes likely reflect greater improvements in the rural areas where elephant populations are most common and IMR is higher, rather than smaller changes from urbanisation58.

It is important to note that spatial differences in average IMR might represent differences between sites in poverty better than annual IMR measures. Annual IMR declines over time as successful medical and public health measures have improved healthcare much faster than other factors associated with poverty have improved, potentially weakening the value of annual IMR as a proxy for overall poverty. We therefore tested the correlation of site-level IMR56 with PIKE in a separate model, in which we neglected the temporal changes in IMR entirely. The results of this supplementary analysis supports the assumption that spatial variation in IMR is a better indicator of poverty than temporal variation (see Supplementary Table 5).

Precipitation: Annual precipitation per site was derived from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS59). In the analysis, we took the natural logarithm of this variable because of its left skewed distribution. This climate variable was included to allow for changes in natural elephant mortality. Variation might arise from two processes. Sites with higher precipitation may identify denser habitats, where finding carcasses due to natural mortality is more difficult, and hence PIKE may be higher due to underestimated natural mortality. Secondly, lower precipitation (within or among sites) may increase natural mortality49,50 and thus lead to underestimated poaching rates because of lower PIKE values.

Corruption perceptions index: Corruption perceptions index (CPI) was derived by Transparency International60 for every country in every year. It represents the perceived level of public sector corruption of a country according to experts and businesspeople. The index uses a scale of 0 to 100, where 0 is ‘highly corrupt’ and 100 is ‘very clean’. We included CPI as a proxy for public sector and political corruption, which has been shown to affect the presence of illegal wildlife activities27.

Poverty density: Poverty density defines the number of people per km2 earning less than US$ 1.25 per day. It represents a measure of relative poverty and thus another proxy of the multidimensional poverty problem. These site level data were provided for the year 2005 by HarvestChoice61.

Site area: Surface area of MIKE sites62 in km2. In the analysis, we took the natural logarithm of this variable because of its left skewed distribution. The expected effect of the site area on poaching intensity is somewhat ambivalent. On the one hand, larger protected areas might exhibit less of the negative edge effect, on the other hand, smaller ecosystems might be easier to patrol.

Law enforcement adequacy: Estimates of the adequacy of law enforcement provision. For each site, MIKE specialists return a form after receiving training from (https://cites.org/eng/prog/mike/tools_training_materials/leca) the MIKE programme team, estimating the adequacy of law enforcement provision. We expected ecosystems with higher law enforcement adequacy to show lower PIKE values.

Large-scale ivory seizures: Annual weight of large scale ivory seizures (≥500 kg)63,64. In cases, where worked ivory was part of the consignment, the values were converted to raw ivory equivalent, factoring in a 30% loss during processing.

Ivory price: Annual mammoth ivory prices in the main Chinese markets (China, Hong Kong and Macao) were derived from the UN Comtrade database65. This covariate was included as a proxy for demand for elephant ivory, as we assume that mammoth ivory prices are correlated with black market elephant ivory prices (for which data do not exist). Yet, it is worth noting that price for ivory is likely not only affected by market demand, but also more general conditions of the economy. To correct the obtained trade values for varying inflation rates, we used World Bank consumer price indices66. The corrected trade values were averaged by taking the market specific net weight into account. Note that Macao only reported mammoth ivory prices for the years 2006–09 and 2014.

Ivory price and impacts of seizures on supply and demand may influence poaching rates over a variety of time-scales. While poachers in Africa may be aware of international trends, it is possible information about markets flows slowly to the field. Consequently, we repeated all our analyses with lags of up to two years in these two variables, as is common within econometric analyses67. In the main results we present the zero-lag model.

Statistical analysis

Inferring elephant poaching intensity from carcass encounter data is difficult when, as here, sampling effort is unknown. Estimating the proportion of illegally killed elephants (PIKE), a relative measure, somewhat reduces this issue, assuming sampling effort is invariant for carcasses of natural and illegal causes in a particular year and site.

To estimate PIKE for each observation i, we assumed the number of carcasses identified as illegal killings (n illegal ) to be a binomial random variable given the total number of elephant carcasses (n total ) and probability p, such that

$$n_{{{\mathrm{illegal,}}i}}\sim{\mathrm{Binomial}}( {p_i,n_{{\mathrm{total}},i}}),$$ (1)

where probability p i (=estimated PIKE) is a function of a set of a priori chosen environmental and socio-economic covariates and year-, country- and site-level normally distributed \(({\cal{N}})\) random intercepts with level-specific means (μ) and standard deviations (σ), transformed using the canonical logit link:

$$\begin{array}{*{20}{l}} {{\mathrm{logit}}\left( {p_i} \right)} \hfill & = \hfill & {\beta _0 + \beta _1\,{\mathrm{ln}}\left( {{\mathrm{Precip}}_i} \right) + \beta _2{\mathrm{IMR}}_i + \beta _3\,{\mathrm{CPI}}_{{\mathrm{country}}[i]}} \hfill \\ {} \hfill & {} \hfill & { + {\cal{N}}\left( {\mu _{{\mathrm{site}}[i]},\sigma _{{\mathrm{site}}[i]}} \right) + {\cal{N}}\left( {\mu _{{\mathrm{year}}[i]},\sigma _{{\mathrm{year}}[i]}} \right) + {\cal{N}}\left( {0,\sigma _{{\mathrm{country}}[i]}} \right).} \hfill \end{array}$$ (2)

To account for the spatial and temporal structure of the data, the hierarchical level means for site (μ site ) and year (μ year ) were modelled in detail such that

$$\mu _{{\mathrm{site}},s} = \beta _4{\mathrm{PovDens}}_s + \beta _5\,{\mathrm{LawEnf}}_s + \beta _6\,{\mathrm{ln}}\left( {{\mathrm{Area}}} \right)_s,$$ (3)

$$\mu _{{\mathrm{year}},y} = \beta _7\,{\mathrm{Seizures}}_y + \beta _8\,{\mathrm{IvoryPrice}}_y.$$ (4)

β n represent the regression coefficients, CPI is the annual country-level corruption perceptions index, PovDens (poverty density), Area (site area) and LawEnf (law enforcement adequacy) are site-level covariates, Seizures (large scale ivory seizures) and IvoryPrice (ivory price) are annual-level covariates, and Precip (precipitation) and IMR (infant mortality rate) are annual site-level covariates.

The model was fitted via Markov Chain Monte Carlo (MCMC) sampling using the software JAGS68, accessed through the R69 package R2jags70. The parameter posterior estimates were derived in three independent chains, each of 100,000 iterations, a burn-in phase of 50,000 iterations and thinned to every 50th sample. Based on estimated \(\hat R\)71 and effective sample sizes the applied MCMC algorithm fully converged (see Supplementary Table 2).

With the objective to build an interpretable model with high predictive capacity of PIKE, we regularised the model using the Bayesian lasso72 instead of applying subset selection. By imposing a penalty proportional to the absolute values of the regression coefficients (L 1 -norm penalty), the lasso73 automates variable selection using continuous shrinkage and leads to a sparse model representation. In Bayesian inference, we achieve this using Laplace priors for the regression parameters β n , such that

$$\beta _n\sim {\mathrm{Laplace}}(\mu = 0,b = \lambda ^{ - 1}),$$ (5)

where the regularisation parameter, λ, represents the inverse of the scale parameter in the Laplace distribution (or the rate in an exponential distribution), resulting in stronger shrinkage with increasing λ. We allowed the model to estimate λ from the data by setting it as a hyperparameter. For its implicit estimation, we imposed a diffuse gamma hyperprior on λ2 to maintain conjugacy:

$$\lambda ^2\sim \frac{{\delta ^r}}{{\Gamma \left( r \right)}}\left( {\lambda ^2} \right)^{r - 1}e^{\delta \lambda ^2},$$ (6)

with shape r = 1 and rate δ = 1, which resulted in a median posterior estimate of λ = 1.64 (90% CI: 1.00–2.42). We also used gamma priors with r = 1 and δ = 1 on the standard deviations of the year-, country- and site-level random effects. We tested the sensitivity of the choice of prior distributions on λ, σ site , σ year and σ country . The regression coefficients showed little difference when imposing uniform instead of gamma prior distributions (compare Supplementary Tables 2 and 3).

For an assessment of the models predictive capacity accounting for potential temporal dependencies31, we first sliced the data into temporal blocks of training and test sets. Training data comprises all records in the period 2002 to 2013 (n training = 447, i.e. ~75%). Test data are all observations between 2014 and 2017 (n test = 160, i.e. ~25%). To validate the model, we estimated PIKE for the period 2014–17 from 3000 MCMC draws of the model fitted to the training data. These estimates were compared against the respective PIKE observations in the test set (Fig. 2b). As a measure of predictive power we calculated R2 weighted by n total

Estimating annual poaching rates

While the proportion of illegally killed elephants (PIKE) overcomes the problem of unknown sampling effort, the rate of illegal killing (i.e. the proportion of illegally killed elephants of the total population) is more intuitive. Burnham52 proposed a simple conversion from PIKE to poaching rate (m p ), given an pre-defined natural mortality rate (m n ):

$$m_{\mathrm{p}} = \frac{{{\mathrm{PIKE}}\,m_{\mathrm{n}}}}{{1 - {\mathrm{PIKE}}}}$$ (7)

As such, the derived poaching rate retains a perfect 1:1 relationship with PIKE. Based on estimates collated by Wittemyer et al.7, we assumed a constant natural mortality of 3% (m n = 0.03), but compared the results to estimates with 2% (m n = 0.02) and 4% (m n = 0.04) natural mortality (see Supplementary Fig. 3). It is worth noting that as PIKE tends to 1, the estimated poaching rate increases exponentially, which might lead to implausibly high poaching rates. Therefore, when estimating continental annual poaching rates, we depicted the median across site-specific annual poaching rates. Site-specific assessment (see Supplementary Results) was based on the estimated PIKE, because in some sites PIKE values were estimated close to 1. Note we did not impose a cap on estimated PIKE, as poaching rates even in large elephant populations can be extremely high74.

To predict the annual continental poaching rate (grey lines in Fig. 2a), we drew 3000 samples from the posterior distribution to estimate PIKE for all surveyed sites and years, translated these into site-by-year poaching rates and took the annual median value. For the observed annual continental poaching rate (black crosses in Fig. 2a), we first summed up all observed carcasses across sites, derived annual continental PIKE values and turned them into annual continental poaching rates. Note that the latter might be biased downwards, because sites that report more carcasses (e.g. due to better resourced ranger patrols), and thus dominate the continentally aggregated PIKE observations, tend to have lower poaching rates than sites with fewer observations.

Identifying conservation targets

To identify potential conservation targets, we estimated the sensitivity of the estimated poaching rate to improvements in the socio-economic factors considered. We used 3000 MCMC samples from the fitted model to predict continental annual poaching rates (or region- and site-specific PIKE; see Supplementary Results) with the predictor values consecutively set to the best (i.e. most elephant friendly) observed value within all 53 sites and 15 years. These were: IMR = 17.73 deaths/1000 infants (Tarangire and Manyara National Parks, Tanzania 2016); CPI = 65 (Botswana 2012); poverty density = 4.85 people km−2 at < US$ 1.25 per day (Lopé National Park, Gabon); site area = 81,046 km2 (Selous and Mikumi National Parks, Tanzania); law enforcement adequacy = 0.83 (Etosha National Park, Namibia); large-scale ivory seizures = 790 kg (2002); mammoth ivory price = US$ 23.72 kg−1 (2002). Thus, the differences among sites and countries (see Supplementary Results) are simply a consequence of the current situation in a site or country relative to the best situation in any site or country between 2002 and 2017, and do not represent different effect sizes among sites and countries.

Spatial and temporal residual autocorrelation

We checked for spatial and temporal residual autocorrelation using the Sncf function in the ncf R package75, which allows for a spatio-temporally structured model to be analysed. Residuals were calculated as the difference of estimated and observed PIKE. We did not consider any of this further as the residuals showed neither consistent spatial nor temporal autocorrelation (Supplementary Fig. 2).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.