Spatial separation of prey from livestock facilitates coexistence of a specialized large carnivore with human land use

# Spatial separation of prey from livestock facilitates coexistence of a specialized large carnivore with human land use [https://doi.org/10.5061/dryad.n5tb2rbv6](https://doi.org/10.5061/dryad.n5tb2rbv6) This dataset contains all data used in the data analysis procedure of this study, including SL\_watershed0.csv to construct the GLMM for modelling snow leopard habitat use, and BSden\_unit0.csv to construct the GLMM for modelling bharal density. The other three files, BS.occupancy.csv, BS.dis.csv, BS.sea.csv, BS.site.csv were used to construct occupancy model for bharal spatial distribution. As authors, we confirm the data didn't contain any location data for any of the species involved, thus won't affect the conservation of snow leopard. ## Description of the data and file structure Abbreviations in SL\_watershed0.csv: ID (valley ID), SL (snow leopard detection number), CD (camera working day), BS (bharal density/km^2), LS (livestock grazing intensity in sheep unit/km^2),EVI (mean grassland EVI), WP (mean winter precipitation),WT (mean winter temperature), DEM (mean elevation), VRM (mean ruggedness), Slope (mean slope), House (house density/km^2), Rock (rock coverage), Site (study site name) For bharal occupancy model: BS.occupancy.csv: Bharal detection status in each 2 by 2 km grid across 8 surveys. BS.site.csv: occupancy covariates. Grid (grid ID), LS (livestock presence), Rock (rock coverage), DEM (elevation), EVI (grassland EVI), Slope, VRM (ruggedness), Site (study site name) BS.sea.csv and BS.dis.csv are two detection covariates used in occupancy model, one indicating the survey seasona and one indicating the distance from the transects to bharal herds. ## Code used for plotting the figures: library(plyr) library(reshape2) library(ggplot2) library(hrbrthemes) #preparing data setwd("D:/phD/DATA\_ANALYSIS/Raw\_data") bs\_ls<-read.csv("BS\_LS\_loc.csv", sep=",") head(bs\_ls) bs\_ls$Species<-as.factor(bs\_ls$Species) bs\_ls$Season<-as.factor(bs\_ls$Season) bs\_ls$Groups<-as.factor(bs\_ls$Groups) attach(bs\_ls) library(car) library(rgl) install.packages("barplot3d") library(barplot3d) scatter3d\(Elevation\~EVI\+Slope\|Species\, data=bs\_ls\, axis\_ticks=T\, surface=F\) melted <- melt(bs\_ls, id.vars=c("Species")) com<-ddply(melted, c("Species", "variable"), summarise, mean = mean(value), sd = sd(value), sem = sd(value)/sqrt(length(value))) com$upper<-com$mean+1.96\*com$sem com$lower<-com$mean-1.96\*com$sem EVI<-subset(com,variable=="EVI") Slope<-subset(com,variable=="Slope") Elevation<-subset(com,variable=="Elevation") \# This part below was just to play around with the data to see how it looks \# scatter plot with confidence intervals p1<-ggplot(EVI, aes(Species, mean))+xlab("Species")+ylab("EVI")+ geom\_pointrange(aes(ymin = lower, ymax = upper))+ theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) p2<-ggplot(Slope, aes(Species, mean))+xlab("Species")+ylab("Slope")+ geom\_pointrange(aes(ymin = lower, ymax = upper)) + theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) p3<-ggplot(Elevation, aes(Species, mean))+xlab("Species")+ylab("Elevation")+ geom\_pointrange(aes(ymin = lower, ymax = upper))+ theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) require(gridExtra) grid.arrange(p1, p2, p3, ncol=3) \# density plot p1 <- ggplot(data=bs\_ls, aes(x=EVI, group=Species, fill=Species)) + geom\_density(adjust=1.5, alpha=.4) + theme\_ipsum(axis\_title\_size = 12, axis\_title\_face = "plain", grid = F, axis\_col = "#cccccc", axis = T, ticks = T) p2 <- ggplot(data=bs\_ls, aes(x=Slope, group=Species, fill=Species)) + geom\_density(adjust=1.5, alpha=.4) + theme\_ipsum(axis\_title\_size = 12, axis\_title\_face = "plain", grid = F, axis\_col = "#cccccc", axis = T, ticks = T) p3 <- ggplot(data=bs\_ls, aes(x=Elevation, group=Species, fill=Species)) + geom\_density(adjust=1.5, alpha=.4) + theme\_ipsum(axis\_title\_size = 12, axis\_title\_face = "plain", grid = F, axis\_col = "#cccccc", axis = T, ticks = T) grid.arrange(p1, p2, p3, nrow=3) \# Script for figure 2 setwd("D:/Manuscripts/Land sharing in mountainous refuge") c<-read.csv("Coefficients\_no GI.csv", sep=",") head(c) c$Variable <- factor(c$Variable, levels = rev(c("Grazing intensity","Bharal density","Winter temperature","Winter precipitation","Elevation","Vision coverage"))) p1<-ggplot(c, aes(Coefficients, Variable, colour=Model, order=ID))+labs(fill="Models")+ scale\_colour\_discrete(limits=c("Snow leopard habitat use", "Bharal density"))+ xlab("Standardized coefficient with 95% CI")+ylab("Covariates")+ geom\_pointrange(aes(xmin = Lower, xmax = Upper),position=position\_dodge(width=0.5))+ geom\_vline(xintercept = 0, linetype = "dashed")+ xlim(-1.2, 2)+ theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) p1 bs<-read.csv("BSoccu.csv", sep=",") head(bs) bs$X <- factor(bs$X, levels = rev(c("p (Survey season)","p (Distance to transect)","psi (Livestock presense)","psi (Elevation)","psi (Slope)","psi (EVI)"))) p2<-ggplot(bs, aes(Coefficient, X, colour=Model)) + scale\_colour\_manual(values=c("goldenrod4", "goldenrod2"))+xlab("Standardized coefficient with 95% CI")+ylab("Covariates")+ geom\_pointrange(aes(xmin = Lower, xmax = Upper),position=position\_dodge(width=0.5))+ geom\_vline(xintercept = 0, linetype = "dashed")+ theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) p2 require(gridExtra) grid.arrange(p1, p2, nrow=2) #PCA analysis for niche partitioning visualization \# Figure 3b bs\_ls<-read.csv("ALL\_BS\_LS\_PCA2.csv", sep=",") head(bs\_ls) library(ggfortify) df <- bs\_ls[1:8] pca\_res <- prcomp(df, scale. = TRUE) autoplot(pca\_res) p1 <- autoplot(pca\_res, data = bs\_ls, colour = 'Species',shape= 'Season',loadings = TRUE, loadings.colour = 'blue', loadings.label = TRUE, loadings.label.size = 4) + stat\_ellipse(geom="polygon", aes(fill = Species\_Season), alpha = 0.2, show.legend = FALSE, level = 0.95)+ theme\_minimal() + theme(panel.grid = element\_blank(), panel.border = element\_rect(fill= "transparent")) p1 #PCA with FactoMineR library(FactoMineR) res.pca = PCA(bs\_ls, quali.sup=9:11, graph=T) summary(res.pca) plot.PCA(res.pca, label = c("ind.sup","quali","var"), habillage="Species") require(factoextra) require(ggplot2) require(tidyr) require(dplyr) require(MASS) require(reshape2) require(cowplot) \# extract pc scores for first two component and add to dataframe bs\_ls$pc1 <- res.pca$ind$coord[, 1] # indexing the first column bs\_ls$pc2 <- res.pca$ind$coord[, 2] # indexing the second column bs\_ls$pc3 <- res.pca$ind$coord[, 3] # indexing the third column write.csv(bs\_ls[,11:14], "PC\_BS\_LS.csv") \# At last we didn't use the figures produced here p <- ggplot(data = bs\_ls, aes(x = pc1, y = pc2, color = Species, shape = Season)) + geom\_hline(yintercept = 0, lty = 2) + geom\_vline(xintercept = 0, lty = 2) + geom\_point(alpha = 0.8) p <- p + stat\_ellipse(geom="polygon", aes(fill = Species), alpha = 0.2, show.legend = FALSE, level = 0.95) + xlab("PC 1 (27.67%)") + ylab("PC 2 (17.95%)") + theme\_minimal() + theme(panel.grid = element\_blank(), panel.border = element\_rect(fill= "transparent")) p \# Wilcoxon test on PC1 and PC2 wilcox.test(pc1\~Species, data=bs\_ls) wilcox.test(pc2\~Species, data=bs\_ls) #nicheRover for niche partitioning visualization (Figure 3c, d) library(nicheROVER) bsls<-read.csv("PC\_BS\_LS.csv", sep=",") head(bsls) aggregate(bsls[2:3], bsls[1], mean) \# generate parameter draws from the "default" posteriors of each fish nsamples <- 1e3 system.time({ bsls.par <- tapply(1:nrow(bsls), bsls$Species, function(ii) niw.post(nsamples = nsamples, X = bsls[ii,2:3])) }) \# 2\-d projections of 10 niche regions clrs <- c("black", "red", "blue", "orange") # colors for each species nsamples <- 10 bsls.par <- tapply(1:nrow(bsls), bsls$Species, function(ii) niw.post(nsamples = nsamples, X = bsls[ii,2:3])) \# format data for plotting function bsls.data <- tapply(1:nrow(bsls), bsls$Species, function(ii) X = bsls[ii,2:3]) p2 <- niche.plot(niche.par = bsls.par, niche.data = bsls.data, pfrac = .05, iso.names = expression(PC1, PC2), col = clrs, xlab = expression("scores")) \# niche overlap plots for 95% niche region sizes nsamples <- 1000 bsls.par <- tapply(1:nrow(bsls), bsls$Species, function(ii) niw.post(nsamples = nsamples, X = bsls[ii,2:3])) \# Overlap calculation\. use nsamples = nprob = 10000 \(1e4\) for higher accuracy\. \# the variable over\.stat can be supplied directly to the overlap\.plot function over.stat <- overlap(bsls.par, nreps = nsamples, nprob = 1e3, alpha = c(.95, 0.99)) #The mean overlap metrics calculated across iteratations for both niche #region sizes (alpha = .95 and alpha = .99) can be calculated and displayed in an array. over.mean <- apply(over.stat, c(1:2,4), mean)\*100 round(over.mean, 2) over.cred <- apply(over.stat\*100, c(1:2, 4), quantile, prob = c(.025, .975), na.rm = TRUE) round(over.cred[,,,1]) # display alpha = .95 niche region \# Overlap plot\. Before you run this\, make sure that you have chosen your alpha level\. clrs <- c("black", "red", "blue", "orange") # colors for each species over.stat <- overlap(bsls.par, nreps = nsamples, nprob = 1e3, alpha = .95) p3 <- overlap.plot(over.stat, col = clrs, mean.cred.col = "turquoise", equal.axis = TRUE, xlab = "Overlap Probability (%) -- Niche Region Size: 95%") \# posterior distribution of \(mu\, Sigma\) for each species nsamples <- 1000 fish.par <- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) \# posterior distribution of niche size by species bsls.size <- sapply(bsls.par, function(spec) { apply(spec$Sigma, 3, niche.size, alpha = .95) }) \# point estimate and standard error rbind(est = colMeans(bsls.size), se = apply(bsls.size, 2, sd)) \# boxplots clrs <- c("black", "red", "blue", "orange") # colors for each species p4 <- boxplot(bsls.size, col = clrs, pch = 16, cex = .5, ylab = "Niche Size", xlab = "Species\_Season") require(gridExtra) grid.arrange(p1, p2, p3, nrow=2, ncol=2)