Multiscale Spatially Explicit Generalized Additive Modeling of Human-Tiger Conflict Risk

doi:10.21203/rs.3.rs-1313541/v1

Download PDF

Research Article

Multiscale Spatially Explicit Generalized Additive Modeling of Human-Tiger Conflict Risk

https://doi.org/10.21203/rs.3.rs-1313541/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Context:

Spatial modeling of human-carnivore conflict has recently gained traction and predictive maps have become a great tool to understand the distribution of present and future conflict risk. However, very few such studies consider scale and use appropriate spatial modeling tools.

Objectives

This study aimed at understanding the ecological predictors of human-tiger conflict and predicting livestock predation risk by reintroduced tigers in Panna Tiger Reserve, Central India, by modelling livestock kill as a function of various tiger relevant ecological variables, at multiple scales employing spatially explicit statistical tools.

Methods

Geostatistical modeling was used to create raster layers of covariates (prey, cover, human activities) following which univariate scaling was done. Livestock loss by tiger was then modelled using geoGAM and spatially explicit conflict risk predicted.

Results

It was found that prey and shrub cover both selected at fine scale, were key ecological determinants of human-tiger conflict. Prey showed an inverse relationship with livestock predation and shrub non-linear, livestock predation increasing with increase in shrub cover, but decreasing beyond a certain point. Thus, in habitats where optimum ambush cover is available but prey presence is low at fine-scale, carnivores are more likely to depredate domestic livestock since livestock have lost most of their anti-predator behaviours.

Conclusions

Livestock kill by tiger is a culmination of predator choice and foraging tactics, and, prey vulnerability and defence mechanism. Thus, the spatially explicit predation risk map can provide a basis to guide adequate mitigation measures, as per the risk probability.

human-carnivore conflict

ecological predictors

prey

cover

domestic livestock

Panna

Spatial modelling of conflict or predictive risk modelling has become one of the major tools to understand human-carnivore conflict (HCC) (Treves et al. 2004; Kaartinen et al. 2009; Kissling et al. 2009; Marucco and McIntire 2010; Edge et al. 2011; Zarco-gonzalez et al. 2013; Soh et al. 2014; Rostro-Garcia et al. 2016; Amirkhiz et al. 2018; Mbiba et al. 2018). Statistical modelling is used to isolate the factors that are related to depredation events, and predict its distribution by extrapolating to similar areas and predicting future conflict risk (Treves et al. 2004; Kaartinen et al. 2009; Treves et al. 2011; Behdarvand et al. 2014; Rostro-Garcia et al. 2016; Broekhuis et al. 2017). The predictive/risk maps so produced help managers identify vulnerable habitats, communities, and species (Kaartinen et al. 2009; Treves et al. 2011; Mateo-tomas et al. 2012; Behdarvand et al. 2014; Davie et al. 2014; Soh et al. 2014; Broekhuis et al. 2017; Amirkhiz et al. 2018). However, since all the factors associated with HCC are neither linearly related to kill occurrence nor come into play at the same scale, scale must be considered when modelling habitat correlates of HCC.

Most ecological relationships are complex and involve several factors. And because these factors range from macro to microhabitat/environmental covariates, all of them cannot be expected to operate at and influence the relationship at the same scale. Thus, ecological relationships are scale-dependent, such that when examined at different spatio-temporal scales the relationship and its interpretation are subject to change (Weins et al. 1989; McGarigal et al. 2016). Therefore, when identifying the factors related to a process or phenomenon they need to be examined at multiple scales, to identify the meaningful scale and make ecologically sound inferences (McGarigal et al. 2016). In the absence of such a multi-scale approach, misleading conclusions may be drawn. Even if we are able to successfully identify the causal factors of a problem, we would not know at which level to intervene, without understanding the scale at which these causal factors are influencing the problem. In which case, selection of the scale at which the variables are meaningfully correlated with the issue becomes as important as, the selection of variables themselves (Mateo Sánchez et al. 2013). Thus, coupled with variable selection, scale optimization should be the first step to predictive modelling.

Multiscale models have been proven to perform better than single-scale models at identifying and predicting relationships between environmental variables and the phenomenon/process under study (Mateo Sánchez et al. 2013; Timm et al. 2016). Thus, multiscale modelling has become an important tool for studying a myriad of ecological and biological processes/problems/groupings including community ecology (Dray et al. 2012), ecological niche modelling/niche/resource partitioning (Hearn et al. 2018; Khosravi et al. 2019), habitat selection (Mateo Sánchez et al. 2013)/ habitat suitability modelling (Store and Jokimäki 2003; Kittle et al. 2018; Khosravi et al. 2019; Rather et al. 2020), predicting indicator species hotspots (Grand et al. 2004) and predicting carnivore dispersal (Krishnamurthy et al. 2016). Even though HCC often involves different variables and complex interactions, very few studies have tried to examine the factors determining HCC at multiple scales (Wilson et al. 2005; Soh et al. 2014; Miller et al. 2015; Rostro-Garcia et al. 2016; Broekhuis et al. 2017). However, these studies have found that scale influences livestock predation risk (Davie et al. 2004), with certain habitat factors influencing livestock depredation at broad scale and others at fine scale (Miller et al. 2015; Rostro-Garcia et al. 2016; Broekhuis et al. 2017). Upon comparison of multiscale model with a single-scale model, studies have concluded that scale optimization improves modelling results of livestock predation risk by large carnivores like tiger (Rostro-Garcia et al. 2016).

However, most studies employ aspatial models to predict predation risk using spatial correlates (Soh et al. 2013; Miller et al. 2015). Since most ecological variables exhibit a certain degree of spatial autocorrelation, therefore, it is important to account for the spatial nature of the data (Griffith 1992; Legendre 1993), when modelling predation risk by carnivores. In absence of which, ‘independence of data points’, a common assumption across most statistical models, is violated, leading to unreliable model outcomes (Legendre 1993; Dale and Fortin 2002). As there is a dearth of studies examining HCC at multiple scales employing appropriate spatial statistical models; in this study we have evaluated the ecological determinants of human-tiger conflict (HTC) and predicted livestock predation risk by tigers in and around Panna Tiger Reserve, Central India. Livestock kill was modelled as a function of various tiger relevant ecological variables, viz. prey, cover, water, and anthropogenic disturbance, at multiple scales employing spatially explicit Generalized Additive Model (GAM) .

Study Area

Panna Tiger Reserve (24°16ʹN to 24°42ʹN and 79°29ʹE to 80°16ʹE) covering an area of 1598.10 km² is situated in the state of Madhya Pradesh in central India. The Critical Tiger Habitat (CTH) or core of the reserve is composed of Panna National Park and Gangau Wildlife Sanctuary covering an area of 576.13 km² (Madhya Pradesh Forest Department 2007), and the buffer covers an area of about 1,021.97 km² (Madhya Pradesh Forest Department 2012). The reserve lies in Vindhyan hills, its altitude ranging between 330 and 540 m a.s.l. (Chawdhry 1996; Rodgers et al. 2002). There are four villages within the national park area, however there are seven villages in sanctuary area and 49 villages in the buffer of the reserve. The current study was focussed on the CTH and two kilometres natural buffer around it (Figure 1).

The original tiger population in Panna was lost due to poaching and as a population recovery measure tigers were reintroduced from neighbouring reserves in 2009 (Sarkar et al. 2016). Since then, the tiger population in Panna has increased exponentially from seven founders to around 60 individuals in 2020.

Multivariate multiscale predictive modelling

GAM while retaining the interpretability of generalized linear models have the flexibility of machine learning algorithms, because it does not assume a linear relationship between dependent and independent variables (Hastie and Tibshirani 1990; Larsen 2015). GAM is, as the name suggests, a generalization of the linear model, in which the linear function of the covariate is replaced with a smooth function (Hastie and Tibshirani 1990). Because of their semiparametric nature, GAMs are much more sensitive to unique data distribution than GLM, allowing for modelling of nonlinear relationships by deriving predictor functions during model building (Härdle and Turlach 1992; Larsen 2015). At the same time to avoid overfitting, one can control the smoothness or ‘wiggliness’ of the predictor function (Larsen 2015). Despite their versatility, GAM has not been explored as much as linear or machine learning models to understand the relationship between livestock depredation by carnivores and the environmental factors (Kaartinen et al. 2009; Miller et al. 2015; Rostro-Garcia et al. 2016; Broekhuis et al. 2017; Amirkhiz et al. 2018). Given that the data for this study does not fit the assumptions of linear models, and the primary aim is to not merely get accurate prediction but identify the drivers of HTC, GAM was selected, to model livestock kill locations as a function of various tiger relevant habitat factors (after determining appropriate scale) (Wood 2017).

Livestock kill quantification and kill site location

In India, domestic livestock depredation by wild carnivores is compensated by the forest department of the state. These livestock compensation records were collected from tiger reserve management to understand the intensity of HTC. But the compensation data was not geo tagged, hence the livestock kill data (that has GPS locations) which is collected by tiger monitoring teams in the reserve was also obtained, for the period 2009-2016. Since in case of Panna, there are several feral livestock, and no distinction is made between feral and domestic livestock in the kill records; kill and compensation data were matched, to obtain reliable locations for domestic livestock kill. Thus, 156 locations were matched for the period 2011-2016. For the purpose of modelling, livestock kill data was treated as presence, while random absence data was generated in GIS. Livestock kill data for the period 2011-2015 was used for training the model and 2016 data was used for testing the model.

Measuring ecological variables

A total of 27 variables were considered for ecological driver modelling, which can be grouped under the following heads

I. Prey

Distance sampling technique was used for tiger prey species population estimation (Buckland et al. 2001). A total of 41 line transects each up to 2 km in length were surveyed in the reserve, during winter 2012-13 and 2013-14. All the line transects were walked in a replicate of three. Data for both the years were combined and prey encounter rate for each transect, calculated. Further analysis was done using ‘conventional distance sampling engine’ in the software Distance 6.2 (Thomas et al. 2010).

II. Vegetation cover

Vegetation indices, viz. canopy cover and shrub abundance were quantified in 15m circular plots laid at every 400 m on the line transects (Jhala et al. 2009). A total of 234 circular plots were laid. Within the plots, ocular estimation of canopy cover was made and shrub cover was scored according to its abundance (0-4).

For further understanding the vegetation cover of the tiger reserve, LANDSAT 8 (OLI/TIRS) scenes for the reserve were downloaded from USGS website for the month of April 2013 (LANDSAT SCENE ID = LC81440432013119LGN01; Download date = 20 April 2015) and November 2013 (LANDSAT SCENE ID = LC81440432013327LGN01; Download date = 15 December 2020). Normalized Difference Vegetation Index (NDVI) (Rouse et al. 1974) was calculated based on these scenes, using Raster calculator in ArcGIS 10.1 (ESRI, 2012) by the formula: NDVI = (Near Infrared - Red)/ (Near Infrared + Red).

III. Water

Normalized Difference Water Index (NDWI) (McFeeters 1996) was calculated with raster calculator in ArcGIS 10.1 (ESRI 2012) using the same scenes as used for calculating NDVI, by the formula: NDWI = (Green-Near Infrared)/ (Green+Near Infrared). Additionally, drainage and water source data were obtained from the forest department and Euclidean distance raster created using Euclidean Distance tool in spatial analyst toolbox in ArcMap 10.4 (ESRI 2016a). Euclidean distance rasters were also created for Ken River and its tributaries, and water sources tagged perennial.

III. Topography

ASTER Global Digital Elevation Model (DEM) data was downloaded from USGS Global Visualization Viewer website. Slope tool in spatial analyst was used to calculate slope from DEM layers in ArcGIS 10.1 (ESRI 2012). Additionally, topographic ruggedness index or terrain ruggedness index (TRI) that measures elevation difference between a cell and mean of its eight neighbouring cells (Riley et al. 1999) was calculated using raster calculator in ArcGIS 10.1 (ESRI 2012), by the formula (Cooley 2016): TRI = SquareRoot (Abs((Square(“3x3max”)-Square (“3x3min”)))).

IV. LULC and forest contiguity

LULC prepared by Forest Survey of India for the entire country at 98m resolution for the year 2009 was procured (FSI 2009). Landscape characteristics and patterns, particularly habitat connectivity, were studied using multiple indices in program FRAGSTATS (ver. 4.2). FRAGSTAT analyses was run using the FSI LULC, three class-level metrics were calculated: PD, LPI (percentage of total landscape area comprised by the largest patch), and CLUMPY (a measure of fragmentation). LPI is a simple measure of dominance. And CLUMPY that ranges from -1 (patch type is maximally disaggregated) to 1 (patch type is maximally clumped) provides an index of fragmentation of the focal class that is not affected by changes in class area (McGarigal 2015).

V. Disturbance

Anthropogenic disturbance indices were quantified in the circular plots (as discussed earlier under section II). In each plot, all the lopped (only branches were cut) and cut (cut to stump) trees were counted. Camera traps were deployed in 109 locations in 2x2 km grids within national park, accounting for 7459 trap nights. The number of tigers, livestock, humans, and vehicles captured in each camera trap was then manually counted and encounter rates (total no. of captures/total trap nights) calculated.

Village and road location data were also obtained from the forest department and Euclidean distance raster, calculated. Additionally, human footprint data was downloaded from Socio Economic Data and Application Centre (SEDAC) website for the year 2009 (Sanderson et al. 2002; Venter et al. 2018). Population census data for 2011 was also downloaded from SEDAC (Balk et al., 2020).

Geostatistical modelling to create rasters: Canopy cover, prey, human, livestock, and vehicle encounter rate, cutting and lopping intensity rates, were interpolated to create rasters, using the Geostatistical wizard in ArcMap 10.4 (Cressie 2015; ESRI 2016a). Four interpolation tools were considered for the purpose: Inverse Distance Weighing, Simple Kriging (SK), Ordinary Kriging (OK), and Empirical Bayesian Kriging (EBK). Statistical measures of correctness (mean prediction error, root-mean-square error, standardized root-mean-square error, average standard error) from the cross-validation were used to compare the kriging algorithms. The best model was considered the one that has the smallest root-mean-squared prediction error (RMSE), standardized mean nearest to zero, the average standard error nearest the root-mean-squared prediction error, and the standardized root-mean-squared prediction error nearest to 1 (ESRI 2016b) (Online Resource 1).

Scale and variable selection: To construct a multiscale model, each of these variables (except human footprint, LULC and FRAGSTATS variable) were resampled at five scales: 30m (the highest resolution available), 50m (average kill drag distance for tiger, as indicated by literature (Karanth and Sunquist, 2000; Miller et al., 2015)), 100m (midpoint between fine and coarse resolution), 350m (maximum kill drag distance from literature), 1200m (coarsest resolution used by us and at which most global environmental data is available). Thus, in total there were 118 variables. All the variables were then extracted for each of the presence/absence points using the Spatial Analyst toolbox of ArcMap 10.4.

For scale and feature selection, firstly, univariate logistic regressions were run, after performing Box-Tidwell procedure to test for logistic regression linearity assumption (Shin and Ying 1994). Secondly, univariate GAMs were run to understand how much r square/deviance was explained by each of the predictors at each of the scales. The scale at which the variable was best explaining the response was selected. Additionally, Information Value and Weight of evidence were employed for feature selection (Good and Osteyee 1974). Results of univariate logistic regression, univariate GAM, and Information value were studied to select the explanatory variables at appropriate scales for constructing the predictive model. The data was then checked for multicollinearity and spatial autocorrelation. To check for multicollinearity among the selected variables, we ran appropriate tests of association (for continuous vs continuous and continuous vs ordinal variables, Kendall’s tau b; for categorical vs continuous variables, logistic regression; and for categorical vs categorical variables, Cramer’s V). Among the correlated variables, selection for inclusion in the model was made based on which variable was better explaining the response. All prey (wild prey plus livestock) was highly correlated to wild prey encounter rate (r=0.812). Since all prey encounter rate was explaining higher deviance of dependent variable, it was selected. Similarly, NDVI and NDWI were moderately correlated (r=0.636), since NDVI was explaining higher deviance, it was selected.

So, after variable selection, the global model consisted of the following 12 variables at these scales: All prey encounter rate (50m), slope (50m), elevation (1200m), slope deviation (100m), NDVI (100m), shrub abundance (50m), canopy cover (30m), distance to water (1200m), livestock encounter rate (ct) (50m), human encounter rate (1200m), distance to road (100m), distance to village (100m) (Figure 2).

Spatial autocorrelation was checked among these variables using Moran’s I statistic and it was found that many variables were spatially autocorrelated (Getis 2007; Moran 1950).

Spatial GAM

To account for the spatial autocorrelation in the data, spatial GAM model was constructed using geoGAM package in R ver. 3.6.3. (Nussbaum and Papritz 2017). It’s a procedure to build parsimonious model based on gradient boosting, smoothing splines and a smooth spatial surface to account for the spatial structure. The GAM for spatial data or geoaddtive model in its full generality is represented by

Where, f_s(s) is a smooth function of spatial coordinates, which accounts for residual autocorrelation (Nussbaum and Papritz 2017).

Since the response variables in this case is binary, Bernoulli distribution is assumed and logit link used

where,

For building parsimonious model geoGAM automatically selects factors, covariates and spatial effects using componentwise gradient boosting, following which model is further reduced using cross validation (Nussbaum and Papritz 2017).

The final model so selected was run on test data, to get model performance measures, area under the curve (AUC) and true skill statistic.

Risk map

A raster with 42m cell size was created and masked to the reserve boundary. It was then converted into points and for each of these points, values for all the explanatory variables were extracted. The final selected model was then run on this data to predict predation risk probability for each of the points. The points were then converted back to raster. Finally, risk predictions were assigned as the value of the raster, to create HTC risk map.

Summary of ecological variables

Prey: All prey encounter rate (including livestock) was 2.98 ± 1.90 n/km (Figure 2). Wild prey encounter rate was 2.60 ± 1.71 n/km, and livestock encounter rate was 0.38 ± 0.61 n/km. Prey density estimate (including livestock) was 52.80 ± 7.10 n/km²for both winters combined. Wild ungulate prey density was 27.78 ± 2.61 n/km² while livestock density was 39.80 ± 10.90 n/km²within the CTH.

Cover: Mean canopy cover was 0.33 and shrub abundance 1.81 (Figure 2). The mean value of NDVI was 0.15 ± 0.03 during summer and 0.22 ± 0.06 during winter (Figure 2).

Water: NDWI during summer was -0.09±0.04, and during winter 0.01±0.07. Euclidean distance to water source including ponds and water falls but excluding river at 1200m resolution was 1378.21m±906.25 (Figure 2). Mean distance to river (at 30m resolution) was 4339.78m±3655.48.

Topography and contiguity: Mean elevation and slope were 350.69±83.84 and 9.58±9.93, respectively (Figure 2). The TRI value was 63.68±32.82. As for FRAGSTATS class-level metrics, PD was highest for open forest. LPI was highest for moderately dense forest. Surprisingly clumpiness was highest for non-forest class and lowest for scrub.

Disturbance: Mean lopping was calculated to be 0.60 ± 1.50, and mean woodcutting was 1.39 ± 3.28. Mean livestock encounter rate was 3.19 ± 9.6 (Figure 2). Vehicle encounter rate was 4.57 ± 7.72. As for human encounter rate, it was 1.02 ± 1.20 (Figure 2). Mean Euclidean distance to road was 4064.82m ± 2524.38 (100m) (Figure 2). Mean distance to human settlement was 2252.61m ± 1279.18 (100m) (Figure 2). Human footprint gave an average disturbance score of 37.14±6.95.

Ecological predictors of HTC

The variables selected to be included in the final geoGAM model were, all prey encounter rate (50m), elevation (1200m), NDVI (100m), shrub abundance (50m), and human encounter rate (1200m). Among these smooth terms, all prey encounter rate and shrub abundance were significant at α=0.05 level, and elevation was significant at α=0.1 level (Table 1). The effective degree of freedom (edf) is higher than 3 for most of the smooth terms, indicating that the wiggliness is high, and relationships, nonlinear (Table 1). Even more is revealed by examining the partial effect plots of smooth terms also called rug plots. A partial effect plot shows the effect of an explanatory variable on the response variable after accounting for the effects of all the other variables included in the model. Upon examination of the partial effect plot for all prey encounter rate, it is revealed that it has an inverse relationship with log odd of livestock kill i.e., the odds of livestock kill by tiger are higher when prey is low (Figure 3 a). In case of shrub abundance, a unique trend was observed, log odds of livestock kill increases with shrub abundance but only till it reaches a certain mark, after which increase in shrub abundance seems to reduce odds of livestock kill (Figure 3 b). In case of NDVI, human encounter rate and elevation, as also indicated by their chi square p values, do not seem to have a significant relationship with odds of livestock kill (Figure 3 c, d, e).

The deviance explained by the model was 44.4% and AUC of the model was 0.91. When run on test dataset the model accuracy was calculated to be 0.65 (Table 2) and AUC was found to be 0.70, indicating that the model had fair amount of prediction capability.

Spatial modelling of HCC has enabled conservationists to visualize where the risk of conflict is high and requires mitigation (Kaartinen et al. 2009; Treves et al. 2011; Zarco-González et al. 2013; Amirkhiz et al. 2018; Broekhuis et al. 2017). Potential habitat/environmental factors identified as drivers of conflict risk can help reduce conflict potential and design targeted mitigation measures (Behdarvand et al. 2014).

Habitat factors that structure the carnivore use of an area are likely to dictate livestock kill by the carnivore and, thereby, HCC. Preferred habitat parameters for tiger have been identified mainly as high prey density, forest contiguity, thick understory, proximity to water, and low human disturbance (Miquelle et al. 1999; Karanth and Sunquist 2000; Sunarto et al. 2012). Among these, past studies have linked HTC with tree cover, elevation/altitude, slope, aspect, proximity to reserve forest, proximity to water, distance to village, distance to road, and density of livestock, settlements, and roads (Li et al. 2009; Ahmed et al. 2012; Soh et al. 2014; Miller et al. 2015; Rostro-Garcia et al. 2016; Struebig et al. 2018; Ramesh et al. 2020).

Our spatial modelling revealed that the potential ecological drivers of livestock depredation by tigers in Panna Tiger Reserve, were prey, and shrub, at fine scale. The model suggests that when prey encounter is low at fine scale (50m), i.e., tiger encounters less prey, it is more likely to predate upon domestic livestock. Low availability of prey has been linked with livestock depredation by carnivores including tiger (Fritts et al. 2003; Bhattarai and Fischer 2014; Burgas et al. 2014; Khorozyan et al. 2015). Moreover, vulnerability of prey influences predator choice (Greene 1986; Onkonburi, and Formanowicz' Jr 1997; Provost et al. 2006; Cresswell et al. 2010), predators are known to select a kill which is easier to catch (Mueller 1977; Lang and Gsödl 2001; Weise et al. 2020). Livestock having lost most of their anti-predator behaviour during the domestication process, are vulnerable to becoming easy prey for predators in the absence of human herders (Linnell et al. 1999; Laporte et al. 2010; Flörcke, and Grandin 2013; Weise et al. 2020). Thus, in predator-occupied habitats where there is low availability of wild prey, if the optimal foraging theory (large prey, high in abundance, easy to catch) is applied (Emlen 1966; MacArthur and Pianka 1966; Werner and Hall 1974), the tiger kills what it can with least effort i.e., livestock.

The second explanatory variable is shrub abundance, which is also selected at the same scale (50m) as prey. Shrub abundance seems to have a unique relationship with livestock kill, resulting in increase of livestock kill up to a certain point after which increase in shrub abundance is actually decreasing the odds of livestock kill. Although it is difficult to explain such a complex relationship, but it can be examined in the light of predation technique of tigers. Tiger is an ambush predator, therefore in areas where there is very low cover, it might be very difficult to make a kill, but in areas where cover is high the chances of success may improve (Greene 1986; Murray et al. 1995; Karanth and Sunquist 2000; Sunquist 2010). Studies done on tiger and other carnivores have also found that livestock predation risk was higher in habitats with high shrub density, because it provides cover for these predators (Davie et al. 2014; Miller et al. 2015). However, if the cover is too dense it is also less likely for grazers like livestock to venture into such patches because they would be devoid of grasses. Thus, making the relationship curve between livestock kill and shrub cover, bell-shaped. Livestock kill by tiger is thus a culmination of predator choice and foraging tactics, and, prey vulnerability and defense mechanism.

Studying ecological drivers of HTC in Panna Tiger Reserve, therefore revealed that in predator-occupied habitat if prey availability is low at fine scale, domestic livestock availability is high and ambush cover is available, the odds of predator depredating livestock become high.

Predation risk map

The risk map produced using spatial modelling shows that the risk is higher in the south eastern part of the reserve encompassing Panna Range and parts of Gehrighat Range (Figure 4).

Summary and recommendation

Ecological drivers of human-carnivore conflict are complex and scale dependent, with the likelihood of conflict being high in large carnivore habitat that has low prey encounter and influx of domestic cattle. In case of Panna, it is suggested that mitigation efforts should be focussed on the administrative units flagged as high risk by our study. Furthermore, detailed study should be conducted to understand lower availability of wild prey and higher availability of livestock in certain parts of the reserve, based on which, where required and deemed feasible, prey augmentation should be considered.

Funding: This study was funded by National Tiger Conservation Authority (NTCA), India and Madhya Pradesh Forest Department, India

Competing interests: It is stated that there is no conflict of interest among authors, funding agency, or with any other party. The authors have no relevant financial or non-financial interests to disclose.

Authors' contributions: All authors contributed to the study conception and design. Funding for the study was secured by Ramesh Krishnamurthy. Material preparation, data collection, and analysis were performed by Manjari Malviya. The first draft of the manuscript was written by Manjari Malviya and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Availability of data and material: The datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request keeping in view sensitivity and data use policy of Government of India.

Ethics approval (include appropriate approvals or waivers): Not applicable.

Consent to participate (include appropriate statements): Not applicable.

Consent for publication (include appropriate statements): Not applicable.

Code availability (software application or custom code): We used ArcGIS 10.4 and R ver. 3.6.3. All R codes are available from the corresponding author on reasonable request.

Ahmed RA, Prusty K, Jena J, Dave C., Vihar SC (2012) Prevailing Human Carnivore Conflict in Kanha-Achanakmar Corridor, Central India. J Zool 7(2): 158–164. https://doi.org/10.5829/idosi.wjz.2012.7.2.6335
Amirkhiz RG, Frey JK, Cain JW, Breck SW, Bergman DL (2018) Predicting spatial factors associated with cattle depredations by the Mexican wolf (Canis lupus baileyi) with recommendations for depredation risk modeling. Biol Conserv 224: 327-335. https://doi.org/10.1016/j.biocon.2018.06.013
Balk D, Montgomery MR, Engin H, Lin N, Major E, Jones B (2020) Spatial Data from the 2011 India Census. Palisades, NY: NASA Socioeconomic Data and Applications Center (SEDAC). https://doi.org/10.7927/gya1-wp91. Accessed 25 August 2020.
Behdarvand N, Kaboli M, Ahmadi M, Nourani E, Mahini AS, Aghbolaghi MA (2014) Spatial risk model and mitigation implications for wolf–human conflict in a highly modified agroecosystem in western Iran. Biol Conserv 177:156-164. https://doi.org/10.1016/j.biocon.2014.06.024
Bhattarai BR, Fischer K (2014) Human–tiger Panthera tigris conflict and its perception in Bardia National Park, Nepal. Oryx 48(4): 522-528. https://doi.org/10.1017/S0030605313000483
Broekhuis F, Cushman SA, Elliot NB (2017) Identification of human–carnivore conflict hotspots to prioritize mitigation efforts. Ecol Evol 7(24): 10630-10639. https://doi.org/10.1002/ece3.3565
Buckland ST, Anderson DR, Burnham KP, Laake JL, Borchers DL, Thomas L (2001) Introduction to Distance sampling: Estimating abundance of biological populations. Oxford University Press, Oxford.
Burgas A, Amit R, Lopez BC (2014) Do attacks by jaguars Panthera onca and pumas Puma concolor (Carnivora: Felidae) on livestock correlate with species richness and relative abundance of wild prey?. Rev. Biol. Trop. 62(4):1459-67.
Chawdhry PK (1996) Management Plan for Panna National Park. Madhya Pradesh Forest Department, Government of Madhya Pradesh.
Chundawat RS, Gogate N, Johnsingh AJT (1999) Tigers in Panna: preliminary results from an Indian tropical dry forest. In: Seidensticker J, Christie S, Jackson P (eds) Riding the tiger: tiger conservation in human-dominated landscapes. Cambridge University Press, Cambridge, pp 123-129.
Cooley SW (2016) Terrain Roughness – 13 ways. http://www.gis4geomorphology.com/ roughness-topographic-position. Accessed January 2017.
Cressie N (2015) Statistics for spatial data. John Wiley and Sons, New Jersey USA.
Cresswell W, Lind J, Quinn JL (2010) Predator‐hunting success and prey vulnerability: Quantifying the spatial scale over which lethal and non‐lethal effects of predation occur. J Anim Ecol 79(3): 556-562. https://doi.org/10.1111/j.1365-2656.2010.01671.x
Dale MR, Fortin MJ (2002) Spatial autocorrelation and statistical tests in ecology. Ecoscience 9(2): 162-167. https://doi.org/10.1080/11956860.2002.11682702.
Davie HS, Murdoch JD, Lhagvasuren A, Reading RP (2014) Measuring and mapping the influence of landscape factors on livestock predation by wolves in Mongolia. J Arid Environ 103: 85-91. https://doi.org/10.1016/j.jaridenv.2014.01.008
Dray S, Pélissier R, Couteron P, Fortin MJ, Legendre P, Peres-Neto PR, Bellier E, Bivand R, Blanchet FG, De Cáceres M, Dufour AB, Heegaard E, Jombart T, Munoz F, Oksanen J, Thioulouse J, Wagner HH (2012) Community ecology in the age of multivariate multiscale spatial analysis. Ecol Monogr 82(3): 257-275. https://doi.org/10.1890/11-1183.1
Edge JL, Beyer DE, Belant JL, Jordan MJ, Roell BJ (2011) Adapting a predictive spatial model for wolf Canis spp. predation on livestock in the Upper Peninsula, Michigan, USA. Wildlife Biol 17(1): 1-10. https://doi.org/10.2981/10-043
Emlen JM (1966) The role of time and energy in food preference. Am Nat 100(916): 611-617.
ESRI (Environmental Systems Resource Institute) (2016a) ArcMap 10.4. ESRI, Redlands, California.
ESRI (Environmental Systems Resource Institute) (2016b) Comparing models. Evaluating interpolation results. ArcGIS Desktop. https://desktop.arcgis.com/en/arcmap/10.4/extensions/geostatistical-analyst/comparing-models.htm
ESRI (Environmental Systems Resource Institute) (2012) ArcMap 10.1. ESRI, Redlands, California.
Flörcke C, Grandin T (2013) Loss of anti-predator behaviors in cattle and the increased predation losses by wolves in the Northern Rocky Mountains. Open J Anim Sci 3(3): 248. http://dx.doi.org/10.4236/ojas.2013.33037
Fritts SH, Stephensen RO, Hayes RD, Boitani L (2003) Wolves and humans. In Mech LD, Boitani, L. (eds) Wolves: Behavior, Ecology and Conservation. University of Chicago Press, Illinois, pp 289-316.
FSI (2009) Indian State of Forest Report (ISFR) 2009. Forest Survey of India. MoEF, GOI, Dehradun, pp 199.
Getis A (2007) Reflections on spatial autocorrelation. Reg Sci Urban Econ 37(4): 491-496.
Good IJ, Osteyee DB (1974) Information, weight of evidence. The singularity between probability measures and signal detection (Vol. 376). Springer-Verlag, Berlin.
Grand J, Buonaccorsi J, Cushman SA, Griffin CR, Neel MC (2004) A multiscale landscape approach to predicting bird and moth rarity hotspots in a threatened pitch pine–scrub oak community. Conserv Biol 18(4): 1063-1077. https://doi.org/10.1111/j.1523-1739.2004.00555.x
Greene CH (1986) Patterns of prey selection: implications of predator foraging tactics. Am Nat 128(6): 824-839.
Griffith DA (1992) What is spatial autocorrelation? Reflections on the past 25 years of spatial statistics. L'Espace géographique 21(3): 265-280.
Härdle WK, Turlach BA (1992) Nonparametric Approaches to Generalized Linear Models. In Fahrmeir L, Francis B, Gilchrist R, Tutz G (eds) Advances in GLIM and Statistical Modelling. Lecture Notes in Statistics, vol 78. Springer, New York, pp 213-225. https://doi.org/10.1007/978-1-4612-2952-0_33
Hastie TJ, Tibshirani RJ (1990) Generalized additive models (Vol. 43). Chapman and Hall/CRC press, London.
Hearn AJ, Cushman SA, Ross J, Goossens B, Hunter LT, Macdonald DW (2018) Spatio-temporal ecology of sympatric felids on Borneo. Evidence for resource partitioning?. PloS one 13(7): e0200828. https://doi.org/10.1371/journal.pone.0200828
Jhala YV, Qureshi Q, Gopal R, Amin R (2009) Field Guide: Monitoring tigers, co-predators, prey and their habitats, 3^rd ed. Wildlife Institute of India and National Tiger Conservation Authority, Dehradun. https://projecttiger.nic.in/WriteReadData/PublicationFile/FIELD%20GUIDE%202013.pdf
Kaartinen S, Luoto M, Kojola I (2009) Carnivore-livestock conflicts: determinants of wolf (Canis lupus) depredation on sheep farms in Finland. Biodivers. Conserv. 18(13): 3503–3517. https://doi.org/10.1007/s10531-009-9657-8.
Karanth KU, Sunquist ME (2000) Behavioural correlates of predation by tiger (Panthera tigris), leopard (Panthera pardus) and dhole (Cuon alpinus) in Nagarahole, India. J Zool 250(2): 255-265.
Khorozyan I, Ghoddousi A, Soofi M, Waltert M (2015) Big cats kill more livestock when wild prey reaches a minimum threshold. Biodivers. Conserv. 192: 268-275. https://doi.org/10.1016/j.biocon.2015.09.031
Khosravi R, Hemami MR, Cushman SA (2019) Multi-scale niche modeling of three sympatric felids of conservation importance in central Iran. Landsc. Ecol. 34(10), 2451-2467. https://doi.org/10.1007/s10980-019-00900-0
Kissling DW, Fernández N, Paruelo JM (2009) Spatial risk assessment of livestock exposure to pumas in Patagonia, Argentina. Ecography 32(5): 807-817. https://doi.org/10.1111/j.1600-0587.2009.05781.x
Kittle AM, Watson AC, Cushman SA, Macdonald DW (2018) Forest cover and level of protection influence the island-wide distribution of an apex carnivore and umbrella species, the Sri Lankan leopard (Panthera pardus kotiya). Biodivers. Conserv. 27(1): 235-263. https://doi.org/10.1007/s10531-017-1431-8
Krishnamurthy R, Cushman SA, Sarkar MS, Malviya M, Naveen M, Johnson JA, Sen S (2016) Multi-scale prediction of landscape resistance for tiger dispersal in central India. Landsc. Ecol. 31(6): 1355-1368. https://doi.org/10.1007/s10980-016-0363-0
Laporte I, Muhly TB, Pitt JA, Alexander M, Musiani M (2010) Effects of wolves on elk and cattle behaviors: implications for livestock production and wolf conservation. PloS one 5(8): e11954. https://doi.org/10.1371/journal.pone.0011954
Larsen K (2015) GAM: the predictive modeling silver bullet. Multithreaded. Stitch Fix 30: 196-223.
Legendre P (1993) Spatial autocorrelation: trouble or new paradigm?. Ecology 74(6): 1659-1673.
Li B, Zhang E, Liu Z (2009) Livestock depredation by Amur tigers in Hunchun Nature Reserve, Jilin, China. Acta Theriol 29: 231-238.
Linnell JDC, Odden J, Smith ME, Aanes R, Swenson JE (1999) Large Carnivores that kill livestock: do “problem individuals” really exist? Wildlife Soc B 27(3): 698–705.
MacArthur RH, Pianka ER (1966) On optimal use of a patchy environment. Am Nat 100(916): 603-609.
Madhya Pradesh Forest Department (August 9, 2012) The Gazette of Madhya Pradesh Extraordinary No.15-21-2010-X-2. Government of Madhya Pradesh, Bhopal.
Madhya Pradesh Forest Department (December 24, 2007) The Gazette of Madhya Pradesh Extraordinary No. 15-31-2007-X2. Government of Madhya Pradesh, Bhopal.
Marucco F, McIntire EJB (2010) Predicting spatio‐temporal recolonization of large carnivore populations and livestock depredation risk: wolves in the Italian Alps. J Appl Ecol 47(4): 789-798. https://doi.org/10.1111/j.1365-2664.2010.01831.x
Mateo Sánchez MC, Cushman SA, Saura S (2013) Scale dependence in habitat selection: the case of the endangered brown bear (Ursus arctos) in the Cantabrian Range (NW Spain). Int J Geogr Inf Sci 28(8): 1531-1546. https://doi.org/10.1080/13658816.2013.776684
Mbiba M, Muvengwi J, Ndaimani H (2018) Environmental correlates of livestock depredation by spotted hyaenas and livestock herding practices in a semi-arid communal landscape. Afr J Ecol 56(4): 984-992. https://doi.org/10.1111/aje.12529
McFeeters SK (1996) The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. Int J Remote Sens 17(7): 1425-1432.
McGarigal K (2015) FRAGSTATS help. University of Massachusetts, Amherst, MA, USA, 182.
McGarigal K., Wan HY, Zeller KA, Timm BC, Cushman SA (2016) Multi-scale habitat selection modeling: a review and outlook. Landsc. Ecol. 31(6): 1161-1175. https://doi.org/10.1007/s10980-016-0374-x
Miller JR., Jhala YV, Jena J, Schmitz OJ (2015) Landscape‐scale accessibility of livestock to tigers: implications of spatial grain for modeling predation risk to mitigate human–carnivore conflict. Ecol Evol 5(6): 1354-1367. https://doi.org/10.1002/ece3.1440
Miquelle DG, Smirnov EN, Merrill TW, Myslenkov AE, Quigley HB, Hornocker MG, Schleyer B (1999) Hierarchical spatial analysis of Amur tiger relationships to habitat and prey. In: Seidensticker J, Christie S, Jackson P (eds) Riding the Tiger: tiger conservation in human-dominated landscapes. Cambridge University Press, Cambridge, pp 71-99.
Moran PA (1950) Notes on continuous stochastic phenomena. Biometrika 37(1/2): 17-23.
Mueller HC (1977) Prey selection in the American kestrel: experiments with two species of prey. Am Nat 111(977): 25-29.
Murray DL, Boutin S, O'Donoghue M, Nams VO (1995) Hunting behaviour of a sympatric felid and canid in relation to vegetative cover. Anim Behav 50(5): 1203-1210.
Nussbaum M, Papritz A (2017) geoGAM: Select Sparse Geoadditive Models for Spatial Prediction. R package version 0.1-2.
Onkonburi J, Formanowicz' Jr DR (1997) Prey choice by predators: effect of prey vulnerability. Ethol Ecol Evol 9(1): 19-25.
Provost C, Lucas E, Coderre D (2006) Prey preference of Hyaliodes vitripennis as an intraguild predator: Active predator choice or passive selection?. Biol Control 37(2): 148-154. https://doi.org/10.1016/j.biocontrol.2006.01.005
Ramesh T, Kalle R, Milda D, Gayathri V, Thanikodi M, Ashish K, Giordano AJ (2020) Patterns of livestock predation risk by large carnivores in India’s Eastern and Western Ghats. Glob. Ecol. Conserv. 24: e01366. https://doi.org/10.1016/j.gecco.2020.e01366
Rather TA, Kumar S, Khan JA (2020) Multi-scale habitat modelling and predicting change in the distribution of tiger and leopard using random forest algorithm. Sci Rep-UK 10(1): 1-19. https://doi.org/10.1038/s41598-020-68167-z
Riley SJ, DeGloria SD, Elliot R (1999) A terrain ruggedness index that quantifies topographic heterogeneity. Intermountain. J. Sci. 5(1-4): 23-27.
Rodgers WA, HS Panwar, VB Mathur (2002) Wildlife Protected Area Network in India: A Review. Wildlife Institute of India publication, Dehradun, India, 44pp.
Rostro-García S, Tharchen L, Abade L, Astaras C, Cushman SA, Macdonald DW (2016) Scale dependence of felid predation risk: identifying predictors of livestock kills by tiger and leopard in Bhutan. Landsc. Ecol. 31(6): 1277-1298. https://doi.org/10.1007/s10980-015-0335-9
Rouse J Jr, Haas RH, Schell JA, Deering DW (1974) Monitoring vegetation systems in the Great Plains with ERTS. NASA Spec. Publ. 351: 309.
Sanderson EW, Jaiteh M, Levy MA, Redford KH, Wannebo AV, Woolmer G (2002) The Human Footprint and the Last of the Wild: The human footprint is a global map of human influence on the land surface, which suggests that human beings are stewards of nature, whether we like it or not. BioScience 52(10): 891–904. https://doi.org/10.1641/0006-3568(2002)052[0891:THFATL]2.0.CO;2
Sarkar MS, Ramesh K, Johnson JA, Sen S, Nigam P, Gupta SK, Murthy RS, Saha GK (2016) Movement and home range characteristics of reintroduced tiger (Panthera tigris) population in Panna Tiger Reserve, Central India. Eur J Wildlife Res 62(5): 537-547. https://doi.org/10.1007/s10344-016-1026-9
Shin RT, Ying JS (1994) Imposing linear homogeneity on Box–Tidwell flexible functional forms. J Bus Econ Stat 12(2): 261-265.
Soh YH, Carrasco LR, Miquelle DG, Jiang J, Yang J, Stokes EJ, Tang J, Kang A, Liu P, Rao M (2014) Spatial correlates of livestock depredation by Amur tigers in Hunchun, China: relevance of prey density and implications for protected area management. Biol Conserv 169: 117-127. https://doi.org/10.1016/j.biocon.2013.10.011
Store R, Jokimäki J (2003) A GIS-based multi-scale approach to habitat suitability modeling. Ecol Model 169(1): 1-15. https://doi.org/10.1016/S0304-3800(03)00203-5
Struebig MJ, Linkie M, Deere NJ, Martyr DJ, Millyanawati B, Faulkner SC, Le Comber SC, Mangunjaya FM, Leader-Williams N, McKay JE, John FAS (2018) Addressing human-tiger conflict using socio-ecological information on tolerance and risk. Nat Commun 9(1): 3455. https://doi.org/10.1038/s41467-018-05983-y
Sunarto S, Kelly MJ, Parakkasi K, Klenzendorf S, Septayuda E, Kurniawan H (2012) Tigers need cover: Multi-scale occupancy study of the big cat in Sumatran forest and plantation landscapes. Plos One 7(1): e30859. https://doi.org/10.1371/journal.pone.0030859
Sunquist M (2010) What is a tiger? Ecology and behavior. In: Tilson R, Nyhus PJ (eds) Tigers of the world. Academic Press, London, pp 19-33.
Thomas L, Buckland ST, Rexstad EA, Laake JL, Strindberg S, Hedley SL, Bishop JRB, Marques TA, Burnham KP (2010) Distance software: design and analysis of distance sampling surveys for estimating population size. J Appl Ecol 47(1): 5–14. https://doi.org/10.1111/j.1365-2664.2009.01737.x.
Timm BC, McGarigal K, Cushman SA, Ganey JL (2016) Multi-scale Mexican spotted owl (Strix occidentalis lucida) nest/roost habitat selection in Arizona and a comparison with single-scale modeling results. Landsc. Ecol. 31(6): 1209-1225. https://doi.org/10.1007/s10980-016-0371-0
Treves A, Martin KA, Wydeven AP, Wiedenhoeft JE (2011) Forecasting environmental hazards and the application of risk maps to predator attacks on livestock. BioScience 61(6): 451-458. https://doi.org/10.1525/bio.2011.61.6.7
Treves A, Naughton‐treves L., Harper, E. K., Mladenoff, D. J., Rose, R. A., Sickley, T. A., and Wydeven, A. P. (2004) Predicting human‐carnivore conflict: a spatial model derived from 25 years of data on wolf predation on livestock. Conserv Biol, 18(1), 114-125. https://doi.org/10.1111/j.1523-1739.2004.00189.x
Venter O, Sanderson EW, Magrach A, Allan JR, Beher J, Jones KR, Watson JE (2018) Last of the wild project, version 3 (LWP-3): 2009 human footprint, 2018 release. NASA Socioeconomic Data and Applications Center (SEDAC), 10, H46T40JQ44, Palisades, NY.
Wiens JA (1989) Spatial scaling in ecology. Functional ecology 3(4):385-397.
Weise FJ, Tomeletso M, Stein AB, Somers MJ, Hayward MW (2020) Lions Panthera leo prefer killing certain cattle Bos taurus types. Animals 10(4): 692. https://doi.org/10.3390/ani10040692
Werner EE, Hall DJ (1974) Optimal foraging and the size selection of prey by the bluegill sunfish (Lepomis macrochirus). Ecology 55(5): 1042-1052.
Wilson SM, Madel MJ, Mattson DJ, Graham JM, Burchfield JA, Belsky JM (2005) Natural landscape features, human-related attractants, and conflict hotspots: a spatial analysis of human-grizzly bear conflicts. Ursus 16(1): 117-129. https://doi.org/10.2192/1537-6176(2005)016[0117:NLFHAA]2.0.CO;2
Wood SN (2017) Generalized additive models: an introduction with R. CRC press, Florida.
Zarco-González MM, Monroy-Vilchis O, Alaníz J (2013) Spatial model of livestock predation by jaguar and puma in Mexico: conservation planning. Biol Conserv 159: 80-87. https://doi.org/10.1016/j.biocon.2012.11.007

Table 1

Approximate significance of smooth terms of final geoGAM model predicting livestock predation by tiger
Smooth terms	Effective degree of freedom (edf)	Ref. df	Chi. sq	p-value
s(lt_all_50) (prey encounter rate at 50 m)	3.17	3.85	16.65	0.002
s(sa_ebk_50) (shrub abundance at 50 m)	3.23	3.97	10.95	0.025
s(hum_1200) (human encounter rate at 1200 m)	3.26	4.02	6.46	0.170
s(ndvi_100) (NDVI at 100 m)	3.18	3.95	5.72	0.190
s(dem_1200) (elevation at 1200m)	3.41	4.13	9.37	0.089

Table 2

True skill statistic of final geoGAM model predicting livestock predation by tiger
Accuracy	0.65
95% CI	0.52-0.77
Kappa	0.30
Sensitivity	0.61
Specificity	0.69
Positive Predicted Value	0.66
Negative Predicted Value	0.65
Prevalence	0.49
Detection Rate	0.30
Detection Prevalence	0.46
Balanced Accuracy	0.65

ESM1.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Multiscale Spatially Explicit Generalized Additive Modeling of Human-Tiger Conflict Risk

Status:

Version 1

Abstract

Context:

Objectives

Methods

Results

Conclusions

Figures

Introduction

Study Area

Methodology

I. Prey

II. Vegetation cover

III. Water

III. Topography

IV. LULC and forest contiguity

V. Disturbance

Results

Discussion

Declarations

References

Tables

Supplementary Files

Status:

Version 1