Socioeconomic factors predict population changes of large carnivores better than climate change or habitat loss

Land-use and climate change have been linked to changes in wildlife populations, but the role of socioeconomic factors in driving declines, and promoting population recoveries, remains relatively unexplored. Here, we evaluate potential drivers of population changes observed in 50 species of some of the world’s most charismatic and functionally important fauna—large mammalian carnivores. Our results reveal that human socioeconomic development is more associated with carnivore population declines than habitat loss or climate change. Rapid increases in socioeconomic development are linked to sharp population declines, but, importantly, once development slows, carnivore populations have the potential to recover. The context- and threshold-dependent links between human development and wildlife population health are challenges to the achievement of the UN Sustainable development goals.

density) estimates. We modelled these time-series with log-linear regressions, where abundance (the response) was loge transformed, and year of abundance estimates was selected as the predictor. We added 1% of the maximum abundance value to all abundances in each time-series to handle cases where abundance equaled zero and could not be log transformed. We included a continuous Ornstein-Uhlenbeck autoregressive process to control for temporal autocorrelation in these models, where abundances are more similar when near in time. We extracted the slope coefficient which represents the annual instantaneous rate of change, sometimes called the population growth rate (rt). There are also other formats of quantitative trends in CaPTrends which fall into three broad data types, all of which we converted into an annual instantaneous rate of change (rt): 1) Finite rate of change = ( ) Where λ represents the mean annual finite rate of change.
2) Estimates of relative abundance change between two points in time (e.g. percentage or fold change in the past 10 years) = (1 + ( /100)) Where P represents the additive percentage change (e.g. a population doubling in size = 100%), and N is the difference in time (in years) between the two estimates of abundance. For fold changes, we first converted the fold change into an additive percentage change.
3) time-series of population change estimates, reported as either population lambdas or percentage changes e.g. in year 1 the population doubled (λ = 2) and in year 2 it halved (λ = 0.5). We back-converted the change estimates into abundance estimates against a constant value of 100. We then fitted log-linear regressions with abundance and year, as in the abundance time-series.
We converted all annual instantaneous rates of change into an annual rate of change percentage to improve interpretability. These rates of change ranged from -75% to 68%, but the majority of values fell within -10% to 10% ( Figure S1a). Alongside the quantitative records, 138 populations in the CaPTrends dataset were only described qualitatively with categories: Increase, Stable, and Decrease.
Population trends were assigned into these categories based on their description in the primary literature i.e. if a population was solely described as stable by the authors in the primary literature, without any form of quantification, we recorded the population as Stable. As such, we don't know what levels of abundance change would be needed to meet each category as each primary source had their own subjective definition which was not made explicit in writing.These records were more common for populations located in traditionally poorer-sampled countries (e.g. with lower human development), so whilst they are less informative (only describing the direction and not the magnitude), we deem them important to reduce known biases 1 . As a result, we used a combination of annual rate of change (%) and qualitative categories as our response in our inference model -see below. Despite including both types of trend, our data still contained many of the same spatial and temporal biases that plague macro scale biodiversity change studies ( Figure 1). Figure S1. Distribution of quantitative (a; N = 985) and qualitative (b; N = 138) trends for the 50 large carnivore species with available data.

Covariates
We extracted a variety of covariates that we considered potential drivers of population change in large carnivores (Table S1). The data used to derive these covariates is described in Table S2. Our covariates fall into four categories: land-use, climate, governance, and traits ( Figure 1). One of the challenges in identifying how covariates impact population trends is matching the spatial scale of the covariate with the population i.e. how much of the population is affected by the covariate. To tackle this problem, we used data on the area of extent of each population to generate a circular distribution zone around the population's coordinate centroid. We refer to this as the 'population area' hereafter.
In populations without a reported extent (N = 347), we searched the locality and location description online to identify the approximate size of the population area. For example, the top result in a Google search for 'Serengeti area' described the location as 30,000km 2 . In some cases, we were able to find a location on a map, but its area was not described in any of the Google search results. To handle cases like this, we assigned the population to one of the following categories: small locations (e.g. towns and counties): 1,000km 2 [N = 123], medium locations (e.g. regions and states): 10,000km 2 [N = 151], and large locations (e.g. countries): 100,000km 2 [N = 73]. For example, one trend was located in an area called 'Thekwane, Kruger National Park, South Africa', which from a Google search appeared to be a small section of Kruger National Park. As the area of Kruger is estimated at 20,000km 2 , we opted to label Thekwane as a small location (1000km 2 ). These categories were arbitrarily chosen with the purpose of trying to compile covariate data at an appropriate spatial scale. Table S1. Proposed direction of effect and rationale for inclusion for each of our 16 main effects and 7 interactive effects. For 'Effect', plus and minus signs indicate that we expect populations to increase or decrease, respectively, as the parameter value increases. An 'Effect' value of zero would indicate that we don't expect a simple main effect, and the variable is only included to allow testing of interactive effects, or to control for features in the data (e.g. the spatial size of population). The rationale is a brief description of why we expect a given effect.

Parameter
Effect Rationale

Land-use
Primary habitat loss -Habitat loss could cause declines in food availability and denning habitat, leading to population declines. Habitat loss could also indirectly cause more human interaction which may lead to population declines e.g. persecution, retaliatory killings etc.
Change in human density -High human density is likely associated with more human-carnivore interactions, which could increase persecution and retaliatory killings.
Change in natural land + Increased natural land could improve food availability and denning habitat.

Governance
Protected area coverage 0 Populations in protected areas are expected to be stable. A positive effect here could mean declines outside protected areas exceed those within. A negative effect could indicate that population growth is greater outside protected areas.
Governance + High governance is associated with carnivore population growth, as high governance countries have: 1) lower population baselines in the available data; and 2) effective rule of law, management of corruption, and species protection will promote population growth.
Change in human development -Rapid growth in development will lead to population declines as the growth results in societal changes e.g. shift to industrialisation, urbanization, reduced tolerance etc.
Note: We only study the change in human development and not governance, as governance is relatively static over the studied timeframes.
Human development + High development is associated with carnivore population growth, as high development countries have: 1) lower population baselines in the available data; and 2) improved quality of life promotes tolerance with carnivores.
War present -War results in poaching and degradation of habitat leading to population declines Table S2. Description and citation of data sources used to derive covariates.

Category Data source Description
Land-use Land-use harmonization 2 (LUH2) 3 Model based estimates of 11 land-use categories, ranging from primary vegetation to urban land-uses. Each pixel (resolution: 0.25 degrees) within the global extent describes the proportional representation of each of the 11 land-use types at a given point in time. This dataset is available for all years from 1850 -2100 with annual time steps.
This dataset is used in the creation of the 'Primary land loss' and 'Change in natural land' covariates.
Land-use Global human settlementpopulation (GHS -POP) 4 Model based estimates of human population counts across a global extent, with a pixel resolution of 0. This dataset contains six governance metrics that combined create a governance score. This score is available for the majority of countries and is measured annually, ranging from 1996 -2020.
Used in the creation of the 'Governance' covariate Governance UN human development index 9 This dataset describes the human development index scores for the majority of countries, is measured annually, and ranges from 1990 to 2020.
Used in the creation of the 'Human development' and 'Change in human development' covariates Governance World database of protected areas 10 This dataset describes the spatial distribution of protected areas from 1960-2020, with protected areas represented by a series of polygons, with annual temporal resolution.
Used in the creation of the 'Protected area coverage' covariate. Used in the production of a variety of trait covariates Given population areas regularly exceeded 10,000km 2 ( Figure S2a), it was not computationally feasible to extract covariates over the entire area; thus, we sampled from a random selection of points within each population area, sampling more frequently in larger areas (range: 13 -295 sampling points, Figure S2b). Random sampling was only used for land-use and climate covariates, as governance covariates are measured at the national level, and all traits (except for population area itself) represent species level averages. The population areas and corresponding sampling points were defined with a Mollweide equal-area projection, but we transformed areas and points back into a WGS84 projection to match all covariate rasters (see below). In all covariates, 'population monitoring period' refers to the period (start and end year) the population was monitored for.

Land-use
We extracted three land-use covariates: Primary habitat loss, Change in natural land, and Change in human density. Primary habitat loss and Change in natural land were derived from the land-use harmonization dataset (Table S2), which reports the annual proportional coverage of 11 land-use types between 1850 and 2015AD, at a 0.25° spatial resolution. To make the land-use types more biologically relevant to predators, we amalgamated a selection of the 11 types into two summarytypes: primary habitat -the sum of 'forested primary' and 'non-forested primary'; and natural landthe sum of 'potentially forested secondary', 'potentially non-forested secondary', 'managed pasture' and 'rangeland'. To estimate Primary habitat loss we found the mean primary habitat across sampling points in each population area for each year in the population monitoring period. We then estimated the rate of loss in primary habitat over time by dividing the rate of loss in each year by the previous year, and then converted this to a percentage loss. We defined the mean Primary habitat loss (%) for each population area as the average across this time-series of loss rates. We followed an identical procedure for Change in natural land. Importantly, our estimate of Change in primary habitat could only decrease or remain stable, whilst natural habitat can fluctuate up and down.
We estimated the Change in human density using the Global human settlement human population raster (Table S2), which describes the human density per km 2 for four years: 1975, 1990, 2000, 2015.
For each year, we reduced the spatial resolution to 0.1° by averaging (mean) over finer-resolution pixels. In order to estimate the Change in human density for each population area's monitoring period, we had to estimate missing human density values (years) in each population area. To do this, we first extracted the mean human density across each population area, in all four of the available years. In each population area, we then used a log-linear regression to model human density (base log transformed) against year, and then predicted human density between 1960 and 2015. We modelled year (predictor) with a cubic fit because human density was non-linear. We then extracted the backtransformed predicted values of human density for all years in each population area. As we were only working with four data points, model predictions were highly uncertain. This uncertainty was included by resampling our model with 100 bootstrap iterations. For each population monitoring period and iteration, we extracted the predicted human densities and estimated the rate of change (%) as calculated for the other land-use covariates. Finally, we calculated the mean human density rate of change (%) across all iterations, as well as the standard deviation, which was used to represent uncertainty in the values within the inference model (see below).

Climate
Our two climatic covariates, Change in extreme heat and Change in drought, describe how the number of months exceeding an extreme heat or drought threshold (respectively) changed between a preindustrial period (1901 -1920) and the population monitoring period (variable time periods between 1960 and 2015). In previous work looking at the link between climate change and population abundance change, it has been common to study the change in the average temperature or precipitation over time like in 16 , instead of extreme events. However, we opted to instead focus on extreme events as we suspected they were more likely to represent direct mortality risk (e.g. exceeding thermal maxima) for species than changes in average temperatures or precipitation.
To derive the Change in extreme heat covariate, we compiled a raster time-series of monthly maximum temperature from 1901 to 2015, at a 0.008° resolution (Table S2). From this, we calculated the mean and standard deviation of the monthly maximum temperature in the pre-industrial period To derive our Change in drought covariate, we required two raster time-series describing the mean monthly temperature and precipitation, both from 1901 to 2015. Whilst the monthly precipitation data (total rainfall in mm) was readily available (Table S2), we had to derive a proxy for the mean temperature data which required two steps. First, we used the monthly maximum and minimum temperatures from CHELSAcruts 5 to derive a raster time-series of midpoint temperature values from 1901 to 2015. Then, we used a linear regression of the mean monthly temperatures, which were available from CHELSA V1.2 6 a subset of years (1979-2013), against the midpoint monthly temperatures to determine how midpoint values could be corrected to represent a proxy of mean values. This regression was nearly a 1:1 relationship with near perfect fit (R-squared = 0.99). Using our proxy mean temperature, we calculated Thornthwaite's evapotranspiration across the raster timeseries, which uses the mean temperature, latitudinal position and number of daylight hours to estimate the evapotranspiration rate 17 . Next, we subtracted this monthly evapotranspiration estimate from the monthly precipitation estimate to produce Thornthwaite's standardised precipitationevapotranspiration index (spei), a standard metric used to describe water availability 18 . We then proceeded to estimate a spei threshold and the mean difference in months overlapping the threshold (pre-industrial vs. population monitoring period) in an identical way to how monthly maximum temperature was used to estimate the Change in extreme heat covariate.

Governance
From the literature, we identified five governance covariates that we considered important to large predator population trends, four of which were measured at the country-level: War present 19 , National governance 20 , Human development 21 , and Change in human development. We used three datasets to populate these covariates: 1) For War-present, we used the UCDP/PRIO Armed conflict dataset (Table S2), which lists conflicts (between 1946-2019) where fatalities per year exceeded 25 and at least one of the parties was governmental. We summarised this dataset into a time-series that describes whether a war was taking place in each country's territory in each year between 1960 and 2016. 2) For Governance, we extracted the world governance indicator metrics (Table S2), which presents six annual governance time-series for each country between 1996 and 2016. 3) Finally, for Human development we sourced the UN human development index (Table S2), which provides an annual time-series describing life expectancy, education level, and income per capita between 1990 and 2016 for 189 countries.
As the governance and human development indicator data only stretch back until 1996 and 1990, respectively, some of the trend data (22% and 7%, respectively) preceded the indicator values. We imputed missing values through a multiple imputation chained equations (MICE) framework 22 . We  Using the imputed datasets, we extracted the mean value across the six governance indicators in each country, year, and imputation chain. We then calculated the mean and standard deviation of this combined governance across the imputation chains to produce an annual governance time-series (and associated error) for each country. For Human development, we averaged over the 50 stored imputation chains to calculate the mean and associated standard deviation for each country and year.
We ensured imputed values fell within the natural constraints of the data (e.g. within the range of 0-1 for human development) and connected seamlessly with the observed human development and governance values ( Figure S4). For governance, there is no way of knowing if imputed values are accurate, but for human development, certain countries have estimates of 'historical index of human development (HIHD)' 23 , which whilst calculated slightly differently, still represent a reasonable proxy of human development change pre-1990. We compare our imputed estimates of human development to the HIHD estimates in Figure S5, where in the majority of cases, HIHD estimates fall within the confidence intervals of our imputed estimates, even capturing complex non-linearity like that observed in Botswana (BWA), and to some extent Rwanda (RWA). Across all countries (including the nine displayed in Figure S5   To address these different means in each country and aid direct comparison between the datasets, we adjusted the mean of HIHD dataset, by adding the difference between the two datasets to the HIHD dataset, per country. After we derived the governance and human development time-series, we began extracting the covariates. For War present, we created a binary variable that described whether war(s) had occurred in the country where the population is located, at any point during the population monitoring period.  As a result, our five traits of interest were reliant on collecting values for 12 common traits.
We sourced values for our traits from three different trait datasets: PanTHERIA, AnAge, and an unpublished large predator trait dataset (Table S2). We used multiple trait datasets to populate missing values at the species level. However, the values sometimes differed between the trait datasets, and in these cases, we created multiple records for the species to capture this uncertainty in the trait value.
As a result, many species had more than one value for a given trait. However, despite using multiple trait datasets, values were still missing for some species in some traits ( phylogeny to estimate missing values -Rphylopars is considered one of the best imputation methods 26 . In our Rphylopars model, we trialled three 27 Carnivora phylogenies to ensure the imputations did not drastically change depending on the phylogeny ( Figure S6). Once we confirmed the choice of phylogeny had little impact, we proceeded with the rest of the analysis only using the phylogeny considered 'best' by 27 . We also included all 12 traits mentioned above in the imputation model, and six other traits to attempt to account for biases in the imputation model, specifically: species area of occurrence, minimum absolute latitude species occurs at, maximum absolute latitude species occurs at, difference in maximum and minimum latitude, maximum mean monthly temperature species occurs at, and minimum mean monthly temperature species occurs at. As the phylogenies we used were not perfectly matched to the CaPTrends and Living Planet Index taxonomies, we corrected synonymous species names in the phylogeny, and where species included in the taxonomy were absent from the phylogeny, we appended the species to a close relative node (inferred from taxonomy).  Figure   S5), and also checked between-trait correlations were acceptable i.e. sufficient variance in the correlation ( Figure S7).

Cleaning data
We opted to remove a selection of the population trend and covariate data as the values were deemed unreliable or unsuitable. Specifically, we removed any population trend records beginning before

Response
The core of our model is a linear regression with fixed and random effects attempting to predict some latent state (unknown annual rates of change; Equation 1). Our observed annual rates of change are realisations of this latent state and are linked to linear model predictions through a state-space structure (Equation 2).

Eq. 2 ro ~ rp + e
Where rp represents the deterministically predicted latent annual rates of change from a linear model (simplified here) containing an intercept and random intercepts ( ), as well as fixed effect beta ( ) coefficients. ro represents the observed (or realised) annual rates of change which are stochastically linked to a normal distribution of rp. The uncertainty in the normal distribution is determined by error e. When building the model, we identified that the model residuals exhibited a heavy tailed tdistribution and so we transformed our responses into a gaussian distribution with an inversehyperbolic sine transformation 31 . Figure S8. Model structure of hierarchical linear model, describing distributions of priors and hyperpriors, as well as the process for incorporating overall error, imputation error, trend weights (see Weighted error below), and censoring within the model. We use five distributions (parameters described in brackets) within the model: normal (μ = mean, σ = standard deviation), beta (shape1, shape 2), half-t distribution (μ = mean, σ = standard deviation, df = degrees of freedom), U/uniform (minimum, maximum), and Bernoulli (B = probability).
Our observed population trends fall into two types: quantitative annual rates of change and qualitative descriptions of change. Both data types were modelled with the same normal error prior (Eq2), but to deal with the different data types, and the unknown values of the qualitative descriptions, we censored the qualitative records to indicate that the true value is unknown, but it occurs within a specified range. We specified these annual rate change ranges as -50% to 0%, -5 to 5%, and 0% to 50% within the decrease, stable and increase categories, respectively.

Equation 3 is a distributional representation of Equation 2
describing how the observed annual rate of change is constrained (or censored) to occur within specific values.
The censoring range thresholds are similar to the range of the observed rates of change (-75% to 68%). Many of the qualitative records address known data biases as they occur in less-well represented regions, species, and time-periods 1 . However, these lower quality records can be more prone to error. As a result, we conduct sensitivity analysis (see Sensitivity analysis below) to assess how including censored observations altered model fit, compared to only using quantitative, and high-quality quantitative (derived from at least three abundance observations), trends.

Eq. 3d. ro ~ Uncensored distribution of quantitative records
Qualitative decrease distribution censored to range from -50% to 0% Qualitative stable distribution censored to range from -5% to 5% Qualitative increase distribution censored to range from 0% to 50% To allow the latent rates of change to vary (e in Equation 2), we set the standard deviation of the latent states normal distribution as a half-t (or half-cauchy) hyperprior (centred at zero, with a standard deviation of 0.001 and one degree of freedom). However, as some observations were likely to be more robust than others, we altered the standard deviation of the latent state normal distributions depending on each observation's quality. Specifically, we varied the standard deviation by multiplying the half-t hyperprior by a deterministic weighting term; essentially the standard deviation was inflated for lower quality observations (see Weighted error below).

Random intercepts
We used a hierarchical model structure to account for phylogenetic and spatial non-independence in the data, including species as a random intercept nested with genus, and country as a random intercept nested within sub-regions, as defined by the United Nations (https://www.un.org/about-us/memberstates). These parameters were fit with a normal distribution centred at zero and their error terms were given a vague uniform hyper prior, with a standard deviation ranging from 1e -10 to 100.

Coefficients and covariates
With a combined 23 covariates and interactive effects, we were conscious of overparameterizing the model. As a result, we split these parameters into three groups: 1) core parameters -which included main effects that have previously been reported as influential, are expected to be influential, or control for other parameters and methodological features; 2) optional parameters -which included main effects we considered interesting but with little evidence to-date of any influence on trends; and 3) interaction parameters -which includes all interaction terms between parameters. Core parameters included: Change in human density, Primary land loss, Population area, Body mass, Change in extreme heat, Governance, and Protected area coverage. These core parameters were included in every model, but we used Kuo and Mallick variable selection 32 to identify important parameters from the optional and interaction groups, where variables were only included in an iteration if they were selected from Bernoulli priors. Our optional parameter group was assigned a Bernoulli prior, which sampled from a beta hyperprior (a = 2, β = 8), such that approximately 20% of optional effects would be included in any iteration, on average, but this could range from 0 -100%. The interaction parameter group had an identical, but separate prior setup. Crucially, this interaction prior was only activated if both main effect parameters were present in the model. For example, for the Change in extreme heat and Change in drought interaction to be selected, it would require Change in drought to be selected from the optional Bernoulli prior, and then the interaction itself would need to be selected from the interactive Bernoulli prior. As variable selection can be highly influenced by the standard deviation of the parameter slope coefficients, we specified the slope standard deviation as a vague uniform hyperprior ranging from 1e -10 to 100.

Imputation uncertainty
Six of the covariates in the model contained missing values that were filled using imputation (see coefficients, except 'War present' which is a categorical variable, we also had to rescale the associated imputation standard deviations. As standard deviations cannot be rescaled in the same way as the imputed estimates, we first converted the standard deviation into confidence intervals, we then z-transformed the intervals using the mean and standard deviation of the covariate, and then back calculated the standard deviation from these intervals.

Weighted error
When developing the model, we were conscious that all not rates of change should contribute equally to the fit. For example, whilst including the censored records could decrease taxonomic and spatial biases in the data, they may also introduce error, as these censored records are unlikely to be as accurate as the quantitative trends. As a result, we included a weight term to inflate the uncertainty in these lower quality records, where the half-t hyperprior discussed above is multiplied by a weight term defined as the inverse of the estimated error in the rate of change. This weight term was developed through simulation (see below), and these simulated error weights inflated the variance around the trend in all low-quality observations, not just the qualitative ones.
When simulating the trend weights, we considered our real trend data to be estimates of true trends with some degree of error. This error would be influenced by the certainty of the population abundance estimates, the sampling intensity (e.g. Is the population sampled every year or only in 50%  Figure S9). We sampled from these newly created abundance distributions to produce new error-prone abundance estimates, reminiscent of real uncertainty in abundance estimation. Secondly, we removed a random sample (between 0% and 90%) of the observations in each time-series, producing time-series' with varying levels of completeness.
We then re-calculated the annual rate of change (%) as in the true trend to produce the estimated trend. The distribution of the estimated trend was largely similar to the true trend obtained from the complete dataset ( Figure S10).  were unavailable, so we were unable to directly calculate the coefficient of variation for each trend.
However, we did have data describing the quality of the sampling and modelling which could act as a proxy for the accuracy of the abundance values. Specifically, we scored trends separately in three areas (Table S4), where trends could only be assigned to one category per area; we then added the score across the three areas: Sampling -how systematically was the population sampled? Modellinghow robust was the approach for modelling abundance values? Low-quality record -does the record meet any of the criteria for being considered low quality? For example, a trend with systematic population sampling (+0.04), in which sampling effort was accounted for (+0.04), and none of the low-quality criteria were met (+0), would be given a coefficient of variation score of 0.08. For an abundance value of 100, this coefficient of variation score would allow the abundance to vary between 75 and 125. Admittedly, our scoring criteria are arbitrary, simply designed to add uncertainty around trends that used less robust methods, rather than describe the true uncertainty in the trend.
However, as these arbitrary values only contribute one feature of three in the weighting system, there impact is likely minimal, and is tested in sensitivity analysis regardless (see below).  Table S4. Scoring criteria used to define a coefficient of variation (CV), uncertainty, in abundances.

Sampling
Method of population sampling is not described or is unsystematic/biased. 0.08 Method of population sampling is systematic. 0.04 All individuals in the population identified. 0.01

Modelling
Method of deriving abundance from population sampling is not described or values are just reported in their raw format.

0.08
Sampling effort accounted for in abundance estimates. 0.04 Abundance derived through complex modelling, or total abundance known. 0.01

Low quality record
Abundance values derived from genetic or harvest data; or the trend is labelled as inaccurate within the primary literature; or trend describes asymptotic instead of observed growth; or trend metric is unconventional.

0.04
After predicting the error in the real trend data using the simulated weight model, we scaled and flipped the values so that 1 indicates low error and 1e-10 indicates high error. These values had to be flipped, as the weight term in our hierarchical linear model ( Figure S8) is multiplied by the precision (i.e. uncertainty) around each trend observations, in which a precision would then be deflated (i.e. uncertainty inflated) if multiplied by a high error trend. For example, for observation A with a low error of 0.9, a precision of 10 would be deflated to 6, whilst for observation B with a high error of 0.1, a precision would be deflated to 1, so A would receive 6 times more weight than B.
To ensure our weight term benefitted the model fit, we conducted sensitivity analysis to compare the model fit under four options: 1) the simulated error weight (described above); 2) weighting by trend sample size, whereby trends derived from more abundance observation are given more weight; and 3) unweighted i.e. all observations are treated equally.

Sensitivity analysis
We conducted sensitivity analysis to test how the different weighting, censoring, and temporal lag options influenced our model results, with the aim of selecting parameters which maximised model marginal and conditional R 2 , whilst also balancing this decision against potential risks. For example, including censored observations may reduce model fit but this could still be worthwhile if it reduces taxonomic and spatial biases. For weighting, we ran models separately under each of the three options, including censored observations and a 5-year lag on all covariates in all cases. After identifying the simulated error weighting as the best option for maximising fit and minimising bias (see Supplementary results) we tested the censoring options, again holding all covariates at the 5-year lag. Including censored observations was valuable, so we included the censored observations when assessing the different temporal lag models, from which we identified that using a 10-year lag improved model fit. In each case, we ran the model through two chains, each with 10,000 iterations and discarding the first 5,000. We thinned the complete chains to store every other iteration (thinning factor of 2). We monitored convergence of key parameters within each model, specifically: standard deviation of the model intercept, standard deviation of beta coefficients, standard deviation of each random effect (regions, countries, genus and species), standard deviation of the overall model error, the optional parameter beta hyperprior, and the interactive parameter beta hyperprior. We ensured the multivariate potential scale reduction factor was less 1.1 across all models in the sensitivity analysis.

Model running
After selecting the simulated error weighting, censored observations, and a 10-year lag from the from the sensitivity analysis (see Sensitivity analysis in the supplementary results), we ran the full model through three chains, each with 120,000 iterations. The first 20,000 iterations in each chain were discarded, and we only stored every 10 th iteration along the chain (thinning factor of 10). We opted for a large chain and burn-in due to the model complexity, and to allow a broad selection of parameter combinations to be tested under variable selection. We assessed convergence of the full model on all parameters monitored in the sensitivity analysis, as well as the model intercept, all 16 main effects and all 7 interactive effects (23 slope coefficients in total). We checked the standard assumptions of a mixed effect linear model (normal residuals and heterogeneity of variance), and tested the residuals to ensure no spatial (Moran's test) or phylogenetic (Pagel's lambda) autocorrelation. We also conducted posterior predictive checks to ensure independently simulated values were broadly reminiscent of model predicted values.
After model running, we calculated how frequently, as a proportion, each of the 23 main and interactive effects occurred within the iterations. For the optional parameters, this was derived by dividing the frequency of occurrence by the total count of iterations. For the interactive parameters, whose inclusion was dependent on the frequency by which their derivative main effects were selected, we divided the frequency of occurrence by the total count of iterations where both derivative main effects were present. We also report the median slope coefficient and associated credible intervals for each of the main and interactive effects, and produce marginal effect plots for a selection of important parameters. These marginal effects hold all other covariates at zero (which is the equivalent of the mean, as covariates were z-transformed). We also display the distribution of the random intercepts e.g. for each region, country, genus and species.

Modelling -Counterfactual scenarios
To explore how observed changes in land-use, climate and human development have influenced population trends, we developed three counterfactual scenarios, where we compared observed population change to predicted population change if land-use, climate, and human development remained static. For instance, in the climate change counterfactual scenario, we predicted each observations population trend using the available covariate data (e.g. land-use, governance and trait covariates), as well as taxa and location data (to provide sensitivity to the models varying random intercepts), but set the climate change covariate data to zero (in this case, change in extreme heat and change in drought). Using this scenario data, we predicted each trend using the global model (all covariate parameters). We then subtracted these counterfactual predictions from the observed trends to define 'Difference in annual rate of change (%)', whereby a positive value indicates carnivore populations would be in better shape (fewer declines) under the counterfactual scenario, and vis-versa.
We summarise counterfactual scenarios by reporting the median Difference in annual rate of change and 95% quantiles across the observed 1,127 populations.
Modelling -Socio-economic development and non-linearity in carnivore trends

Sensitivity analysis
We assessed how trend weighting, including censored observations, and specifying a lag period on the covariates influenced model fit and inference, in part to assess if results were particularly sensitive to specific parameters, but also to help choose the parameters which optimised fit and spatio-taxonomic coverage. Using censored observations and a lag period of 5 years on covariates, model fit was greater when using the simulated error weights compared to the unweighted model and the model weighted by sample size (Table S5). Using simulated error weights and a lag period of 5 years on covariates, we then tested the impact of including censored observations which showed higher marginal and conditional R 2 when censored observations were included. Using only high quality time-series (compared to including censored observation) resulted in a higher conditional R 2 , but at the cost of excluding 19 countries and 2 species from the dataset. We considered the gain in model fit did not outweigh the added spatial and taxonomic coverage. Finally, using simulated error weights and the full dataset (including censored observations), we tested how the lag period of covariates influenced fit. All lag periods offered a similar fit (marginal and conditional R 2 ), and so we selected the lag most supported by the literature -10 years, with suggestions peak population change occurs 8 years after environmental change (specifically forest loss) in mammals 27 . In our final model, we used the simulated error weights, included censored observations, and used a 10-year covariate lag. While these decisions optimized fit and data coverage, the type of weightings, the data or time lags used had little impact on inference, as model coefficients were largely similar across all parameter types (Figure S13 -S16).

Parameter frequency and effects
Under variable selection within our model, some parameters occurred far more frequently than others ( Figure S16). Generally, parameters with an effect at the 95% credible intervals were likely to occur more frequently, which is a consequence of the model identifying and favouring the inclusion of important variables. Change in human development was the most common optional parameter, occurring in all of the selected iterations. In contrast, Change in extreme heat * Protected area coverage was the most common interaction parameter ( Figure S16). Figure S16. a) Proportion of iterations where parameter is present. Core parameters (purple) are present in all models, whilst optional and interaction parameters (pink) only occur in models when they are selected from Bernoulli distributions. b) Standardised coefficients of main effect and interactive model parameters. Parameters in purple have an effect at the 95% credible interval, whilst those in grey do not. The different credible interval thresholds are shown for each parameter, with darkest centre representing the 50% credible intervals, followed by the 80%, 90%, and 95% thresholds (the maximum and minimum point on each bar).

Model assumptions and checks
The inference model passed all standard linear mixed effect model assumptions, with residuals not showing signal of spatial or phylogenetic autocorrelation ( Figure S17). However, the inference model failed to represent the more extreme observed annual rate of change (%) values, with the predicted rates of change largely ranging from -2 to 2, whilst the observed rates of change range -5 to 5; both on an inverse hyperbolic sine scale ( Figure S18a). The qualitative-censored predictions largely agreed with the observed data, where censored-increasing values were primarily predicted to increase, and censored-stable values exhibited small increases and decreases ( Figure S18b). The only category with reasonably poor alignment was censored-decreasing, where populations were predicted to be both increasing and decreasing. Our posterior predictive checks suggest the model produced broadly plausible values, with the independently simulated values occurring within the distribution of the observed trends, albeit failing to represent extreme values ( Figure S18c). For the qualitative-censored trends, the quasi-observed values matched the simulated values almost identically ( Figure S18d), which is to be expected. The inference model had a median root mean square error of 9.2%, a median marginal R 2 of 0.15, and a median conditional R 2 of 0.4 ( Figure S19).  Spatial and temporal biases in the predictors As our data exhibit spatial and temporal biases, there is a risk that our predictors don't fairly represent past climate, land-use and governance changes. For instance, if population trends only occur in countries with high governance, our ability to assess features like the impact of human development change on population trends will be limited, as we will not capture the full predictor parameter space.
To assess if we are capturing the parameter space for our most influential predictor (human development change), we compare the observed values in our models to a random sample (N = 1000) at each of the following scales: global, in the United States, and in Tanzania  year. At the national scale, we skipped the country selection stage and extracted human development change solely from the relevant country at different randomly selected time points. We repeated this process 1000 times ( Figure S20). This analysis shows that our observed human development change data closely matches the random sample at both the global and national scales. This analysis suggests, that whilst spatial and temporal biases are present, we are at least capturing much of the parameter space of our most influential variable Quadratic modelling of human development As earlier work has found that human development has a quadratic effect on wildlife population trends 21 -a feature we do not test in our model presented in Figure S8 -here, using the quantitative records (N = 983) and a simplified model, we test if our significant human development change result persists when we allow the human development score to modelled with a quadratic term. Specifically, we regressed annual rates of change (%; inverse hyperbolic-sine transformed) against human development change and a second-order polynomial of human development (i.e. quadratic), with a random intercept of species nested in genera and countries within UN sub-regions. This simplified model supported the main conclusions from the manuscript, where human development change remains a more important factor than human development score, even after accounting for a quadratic relationship between rate of change and human development score (Table S6).