- Smithsonian National Zoo and Conservation Biology Institute, Conservation Ecology Center, 1500 Remount Rd, Front Royal, VA 22630, USA
- Institute of Crop Science, University of Hoffenheim, Stuttgart, Germany
- Lolldaiga Hills Research Programme, P.O. Box 26, Nanyuki, Kenya
- Sustainability Research Institute, School of Earth and Environment, University of Leeds, Woodhouse Lane, Leeds, LS9 9JT, UK
- Directorate of Resource Surveys and Remote Sensing, P.O. Box 47146-00100, Nairobi, Kenya
- Conservation and Ecology Group, University of Groningen, PO Box 11103, 9700 CC, Groningen, the Netherlands.
- Mpala Research Centre, P.O. Box 555, Nanyuki, Kenya
- Ecology and Evolutionary Biology, Princeton University, 106A Guyot Ln, Princeton, NJ 08544, USA
Access original publication at:
https://www.sciencedirect.com/science/article/abs/pii/S0006320719314466 link
Abstract
Private lands are critical for maintaining biodiversity beyond protected areas. Across Kenyan rangelands, wild herbivores frequently coexist with people and their livestock. Human population and livestock numbers are projected to increase dramatically over the coming decades. Therefore, a better understanding of wildlife-livestock interactions and their consequences for biodiversity conservation on private lands is needed. We used a Bayesian hierarchical, multi-species and multi-year occupancy model on aerial survey data of 15 wild-herbivore species, spanning 15 years (2001-2016) to investigate a) spatiotemporal trends in species occurrence and richness across a mosaic of properties with different land uses in Laikipia County, central Kenya; and b) the effects of distance to water, vegetation and livestock relative abundance on species occurrence and richness. Although mean herbivore species richness varied little over time, we observed high spatial variation in species occurrence across Laikipia, mainly driven by negative effects of high livestock relative abundance. As expected, ‘wildlife friendly’ properties had higher herbivore species richness than other areas. However, high variability suggests that some pastoral properties support rich herbivore communities. The area occupied by five species with global conservation concerns (reticulated giraffe, Grevy’s zebra, Beisa Oryx, Defassa waterbuck and gerenuk) and for which Laikipia County is one of the last refuges was <50% across years. We conclude that ‘wildlife friendly’ properties remain crucial for conservation, although some pastoralist areas offer suitable habitats for wild herbivores. Effective management of stocking rates is critical for maintaining ecosystems able to sustain livestock and wildlife on private lands, ensuring protection for endangered species.
Introduction
Protected areas are essential for the conservation of global biodiversity (Watson et al., 2014). Yet, the effectiveness of the global system of protected areas is recognized as largely insufficient (Ceballos et al., 2005). Few reserves are large enough to satisfy the home range requirements of many species (Caro and Paul, 2007), while fewer still incorporate the broad-scale variation in resources necessary to maintain large-scale seasonal migrations (Fynn and Bonyongo, 2010; Tack et al., 2019) and/or ensure protection of large populations that are critical for long-term persistence (Newmark, 2008; Western et al., 2009b). Consequently, significant populations of large mammals occur on lands that lack formal protection (Ceballos et al., 2005; Ogutu et al., 2016, Western et al., 2009b). Enhancing conservation actions in private lands where wildlife coexist with human activities will become increasingly important if wildlife are to persist into the future (Drescher and Brenner, 2018; Nelson, 2008).
Across grassland ecosystems globally, large herbivores play critical roles in maintaining biodiversity (Mortensen et al., 2018; Post, 2013). The spatiotemporal heterogeneity in structure, productivity, phenology and composition of plant communities dictates the diversity and population stability of large herbivores (Fynn et al., 2016; Owen-Smith, 2004). Grazing by large herbivores, at the same time, improves pasture quality and maintains heterogeneity, enabling grass-dominated ecosystems to support more herbivore biomass than other terrestrial habitats (Frank et al., 1998; Knapp et al., 1999; Huntly, 1991). Across Africa, savannahs support more diversity and abundance of large herbivore species than any other continent, with wildlife coexisting alongside pastoralists and their livestock (du Toit and Cumming, 1999; Reid et al., 2008). Ecosystem benefits has been shown to increase when livestock are kept at moderate densities (Keesing et al., 2018). The presence of livestock, for instance, can reduce ectoparasites abundance when using acaricides (Tallis et al., 2017), improve pasture quality (Young et al., 2018) and promote vegetation heterogeneity through concentration of nutrients (glades; Augustine et al., 2003). At higher stocking rates, however, livestock can negatively affect wildlife due to competition for forage, water and space (du Toit and Cumming, 1999; Georgiadis et al., 2007; Prins, 2000; Western et al., 2009a).
Laikipia County in central Kenya represents an example of successful conservation across private lands where wildlife, people and livestock coexist. This rangeland supports an abundant wildlife community, second only to the Greater Mara ecosystem in Kenya (Ogutu et al., 2016). In Laikipia, local people find economic benefits from the integration of wildlife tourism with livestock commercialization (Keesing et al., 2018). However, coexisting with wildlife can also pose serious challenges, such as economic losses due to increased human-wildlife conflict, the transfer of zoonotic diseases and competition for resources (du Toit et al. 2017). The rapid human population growth across the region, together with a complex process of sedentarization, erosion of traditional governance strategies and an increase in livestock densities and associated overgrazing (Letai and Lind, 2013), is exacerbating pressures on wildlife and the ecological stability of the ecosystem.
Across Laikipia’s rangelands, there has been a general decline in the density of large-wild herbivores within properties that do not actively protect wildlife, as well as evidence of spatial segregation of wildlife and livestock over a 21-year period (1985-2005; Georgiadis et al., 2007). However, the species studied accounted for <15% of regional large mammal richness (Kinnaird and O’Brien, 2012). Kinnaird and O’Brien (2012) assessed a larger number of species using camera trap data on eight properties across Laikipia (10% of the county) during 2008-2010 and found segregation in species occupancy and richness by land use, noting that species were concentrated in ‘wildlife friendly’ properties. To our knowledge, spatially explicit temporal trends in herbivore species richness and species-specific occupancy and their drivers have not been investigated across this system. Yet, this information is critical to provide an improved understanding of the interplay between wildlife communities and the increasing livestock abundance that can guide future conservation actions in rangeland ecosystems.
In this study, we expand upon previous research by investigating a) trends in herbivore occupancy and richness across different land uses in Laikipia County over the past two decades and b) how environmental variables and livestock relative abundance affect herbivore occupancy and species richness dynamics. For this we used the aerial survey dataset collected by Kenya’s Directorate of Resource Surveys and Remote Sensing (DRSRS) and a Bayesian hierarchical multi-species and multi-year occupancy modeling approach (Dorazio et al., 2006; Goijman et al., 2015; Royle and Dorazio, 2008; Zipkin et al., 2009). It is well recognized that animals are imperfectly detected during aerial surveys (Jachmann, 2002; Schlossberg et al., 2018), leading to many species with too few sightings for individual analysis (Royle and Dorazio, 2008; Zipkin et al., 2009). Recent advances in statistical modeling of species occurrence address these problems by accounting for detection probability through repeated sampling and improving parameter estimates by sharing information across species – important for species with few detections (Dorazio et al., 2006; Royle and Dorazio, 2008; Zipkin et al., 2009).
Materials and Methods
Load libraries and set the working directory
library(raster)
library(rgdal)
library(tidyr)
library(spatialEco)
library(dplyr)
library(GISTools)
library(geosphere)
library(jagsUI)
library(R2OpenBUGS)
library(coda)
library(spdep)
library(RColorBrewer)
library(ggplot2)
library(mgcv)
library(gridExtra)
library(Hmisc)
library(tmap)
library(tmaptools)
library(sp)
library(heatmaply)
library(classInt)
library(FSA)
library(DescTools)Study area
Laikipia County (Fig. 1), Kenya, encompasses c. 9700 km2, covered by a mosaic of semi-arid savannah and woodlands. The vegetation composition varies according to the two main soil types; clay-rich ‘black cotton’ soils (vertisols) and well-drained sandy ‘red soils’ (luvisols). In ‘black cotton’ soils, woodlands are primarily composed of Acacia drepanolobium. Grass cover is dominated by Pennisetum mezianum, P. stramineum, Bracharia lachnantha, Themeda triandra and Lintonia nutans (Young et al., 1997). In ‘red soils’, the most important woody species include A. mellifera, A. etbaica and A. brevispica. The more sparse grass cover in these areas is characterized by an abundance of Digitaria milanjiana, Cynodon dactylon, P. stramineum and and Tetrapogon (formerly Chloris) roxburghiana (Butynski and de Jong, 2014). Dry forests exist across Laikipia but are restricted to elevated areas and dominated by Juniperus procera, Olea europaea, Afrocarpus gracilior, Euclea divinorum, Acokanthera schimperi and Croton megalocarpus (Butynski and de Jong, 2014).
Mt. Kenya (5199 m.a.s.l.) in the south-east and the Aberdare highlands (3990 m.a.s.l.) in the south-west determine the rainfall gradient, ranging from 1200 mm of rainfall close to the mountain ranges to 400 mm in the northern areas of the county (Butynski and de Jong, 2014). Rainfall is weakly trimodal, with the majority of precipitation falling in March-May (long rains) and October-December (short rains). Occasional continental rains occur in June-September, particularly in the west. The main dry season is January-February (Schmocker et al., 2015). From Mt. Kenya and the Aberdare highlands flow a series of permanent and seasonal rivers across the plateau. Rainfall patterns are changing, with impacts on both wildlife and livestock movement across this landscape (Schmocker et al., 2015).
A long period of colonial and post-colonial policies has resulted in a dynamic mosaic of land uses and land tenures across the region. These properties include communally owned pastoralist areas, privately owned agricultural plots, commercial cattle ranches and conservancies (Sundaresan and Riginos, 2010). For this study, we focused on the semi-arid rangelands, excluding most of the areas dedicated to small-scale agriculture where most wildlife has been extirpated (Georgiadis et al., 2007). We classified properties into four categories (‘wildlife only’, ‘ranching and wildlife’, ‘ranching’, and ‘pastoralist’) following the County Government of Laikipia (2018) standards (Fig. A1). ‘Wildlife only’ properties are areas dedicated exclusively to wildlife conservation. In these areas, livestock are generally absent. ‘Ranching and wildlife’ properties include conservancies and ranches that are managed to protect wildlife, to conduct science and for tourism, but that also manage livestock at moderated stocking levels. Collectively, these two land uses are referred to as ‘wildlife friendly’. ‘Ranching’ properties are exclusively dedicated to commercial ranching, and in general, do not tolerate wildlife. ‘Pastoralist’ areas are private, communal or government lands in which livestock production is the primary economic activity for people, resulting in areas managed with high levels of livestock (e.g., > 25 total livestock units per 1 km2). For further information on land uses see: Sundaresan and Riginos (2010) and Kinnaird and O’Brien (2012).Laikipia County, Kenya, encompasses c. 9700 km2 of plateau area, covered by a mosaic of semi-arid savannah and woodlands. The vegetation composition varies according to the main two soil types; clay-rich ‘black cotton’ soils (vertisols) and well-drained sandy ‘red soils’ (luvisols). In ‘black cotton’ soils, woodlands are primarily composed of Acacia drepanolobium. Grass cover is dominated by Pennisetum mezianum, P. stramineum, Bracharia lachnantha, Themeda triandra and Lintonia nutans (Young et al., 1997). In ‘red soils’, the most important woody species include A. mellifera, A. etbaica and A. brevispica. The more sparse grass cover in these areas is characterized by an abundance of Digitaria milanjiana, Cynodon dactylon, P. stramineum and Cloris roxburghiana (Butynski and de Jong, 2014). Dry forests exist across Laikipia but are restricted to elevated areas and dominated by Juniperus procera, Olea europaea, Afrocarpus gracilior, Euclea divinorum, Acokanthera schimperi and Croton megalocarpus (Butynski and de Jong, 2014).
Fig. 1. The study area in Laikipia County, Kenya. The county is divided into a mosaic of properties with different land uses including small agricultural plots (Mixed farming), pastoralism, commercial cattle ranches (Ranching) and pro-wildlife ranches and conservancies (Ranching and wildlife). A few properties are exclusively used for wildlife conservation (Wildlife only). The study is focused on the semi-arid rangelands excluding most of the mixed farming properties where most wildlife has been extirpated.
## Reading layer `Property_boundary_2015' from data source `C:\Users\CregoRD\Documents\GitHub\CommunityOccupancy\Shapefiles\Property_boundary_2015.shp' using driver `ESRI Shapefile'
## Simple feature collection with 315 features and 18 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 186885.9 ymin: -32707.34 xmax: 321651.3 ymax: 96041.66
## projected CRS: WGS 84 / UTM zone 37NAerial survey dataset
To evaluate the spatiotemporal dynamics of herbivore richness and occupancy across the region, we used aerial survey data collected during 2001, 2004, 2006, 2008, 2010, 2012, 2015 and 2016 by DRSRS, in partnership with the Laikipia Wildlife Forum (LWF) and the Mpala Research Centre. Aerial surveys are generally conducted across Kenya during the end of the dry season (February-March), with the entire county of Laikipia surveyed in 5-14 days (2001: 10-23 February; 2004: 20 February-1 March; 2006: 15-22 February; 2008: 6-10 March; 2010: 16-26 February; 2012: 4-9 March; 2015: 17 Feb-12 March; 2016; 14-19 April). In addition to the pilot, the survey crew consists of one front-seat observer and two rear-seat observers.
Flights are conducted at a constant speed of approximately 200 km/h at 120 m above ground level. Rear-seat observers count all animals detected between a 141 m strip width demarcated by parallel rods on the wing struts on each side of the plane. For herds >10 animals, photographs are taken for later corroboration of group size. Parallel transects regularly spaced 2.5 km apart were flown north-to-south following the Universal Transverse Mercator (UTM) coordinate system. Transects were sub-divided at 5 km sub-unit intervals. For more details about the aerial survey, see Georgiadis et al., (2007) and Ogutu et al. (2016).
To obtain the replication needed to estimate detectability, we arranged the aerial survey data in 288, 5 × 5 km cells (primary sampling units), each composed of two transect segments (secondary sampling units) of 2.5 × 5 km (Fig. 2). Using spatial, as opposed to temporal, replicates poses the risk of violating the assumption of constant occupancy status across replicates during the sampling period, confounding non-detection with absence and inflating occupancy estimates (Guillera-Arroita, 2011; Kendal and White, 2009). In this study, we assumed that if the species was present at one spatial replicate, it was also present at the second one (i.e., constant occupancy in the cell), given that ungulates are highly mobile species during the dry season (Owen-Smith, 2014) and transects were visited sequentially (Kendal and White, 2009).
We considered 15 species of wild herbivores for which we had detection data, including nine grazers (African buffalo (Syncerus caffer), Grevy’s zebra (Equus grevyi), plains zebra (Equus quagga), hartebeest (Alcelaphus buselaphus lelwel), Defassa waterbuck (Kobus ellipsiprymnus defassa), Grant’s gazelle (Nanger granti), Thomson’s gazelle (Eudorcas thomsonii), common warthog (Phacochoerus africanus), and ostrich (Struthio camelus)), two browsers (reticulated giraffe (Giraffa reticulata) and gerenuk (Litocranius walleri)) and four mixed-feeders (savanna elephant (Loxodonta africana), eland (Taurotragus oryx), impala (Aepyceros melampus), and Beisa oryx (Oryx beisa)).
Fig. 2. Primary (cells) and secondary (segments) sampling units used to assess spatiotemporal dynamics of 15 species of the herbivore community of Laikipia, Kenya, using a Bayesian hierarchical multi-species and multi-year occupancy models.
Data pre-preparation
#Loading and orginizig aerial survey data per year
Lai2001 <- read.csv("./Data/Laikipia2001.csv", header = T)
Lai2004 <- read.csv("./Data/Laikipia2004.csv", header = T)
Lai2006 <- read.csv("./Data/Laikipia2006.csv", header = T)
Lai2008 <- read.csv("./Data/Laikipia2008.csv", header = T)
Lai2010 <- read.csv("./Data/Laikipia2010.csv", header = T)
Lai2012 <- read.csv("./Data/Laikipia2012.csv", header = T)
Lai2015 <- read.csv("./Data/Laikipia2015.csv", header = T)
Lai2016 <- read.csv("./Data/Laikipia2016.csv", header = T)
#Besides years 2001 and 2010, we need to sum the observations of both right and left observers per transect.
Lai2001[is.na(Lai2001)] <- 0 #Replace NAs for 0
Lai2004[is.na(Lai2004)] <- 0 #Replace NAs for 0
Lai2004<-as.data.frame(Lai2004 %>% group_by(UNIT) %>% summarise_each(funs(sum))) #Remove columns and add counts per observer
Lai2006[is.na(Lai2006)] <- 0 #Replace NAs for 0
Lai2006<-as.data.frame(Lai2006 %>% group_by(UNIT) %>% summarise_each(funs(sum)))
Lai2008[is.na(Lai2008)] <- 0 #Replace NAs for 0
Lai2008<-as.data.frame(Lai2008 %>% group_by(UNIT) %>% summarise_each(funs(sum)))
Lai2010[is.na(Lai2010)] <- 0 #Replace NAs for 0
Lai2012[is.na(Lai2012)] <- 0 #Replace NAs for 0
Lai2012<-as.data.frame(Lai2012 %>% group_by(UNIT) %>% summarise_each(funs(sum)))
Lai2015[is.na(Lai2015)] <- 0 #Replace NAs for 0
Lai2015<-as.data.frame(Lai2015 %>% group_by(UNIT) %>% summarise_each(funs(sum)))
Lai2016[is.na(Lai2016)] <- 0 #Replace NAs for 0
Lai2016<-as.data.frame(Lai2016%>% group_by(UNIT) %>% summarise_each(funs(sum))) Load aerial survey units shapefile to match the aerial survey data.
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\CregoRD\Documents\GitHub\CommunityOccupancy\Shapefiles", layer: "SurveyUNITS"
## with 707 features
## It has 4 fieldsMerge the survey information with the sample units. The sample unit is the unique number of each 5 km transect. This will result in polygons with the abundance data as attributes.
L2001<-sp::merge(UNITS,Lai2001, by="UNIT")
L2004<-sp::merge(UNITS,Lai2004, by="UNIT")
L2006<-sp::merge(UNITS,Lai2006, by="UNIT")
L2008<-sp::merge(UNITS,Lai2008, by="UNIT")
L2010<-sp::merge(UNITS,Lai2010, by="UNIT")
L2012<-sp::merge(UNITS,Lai2012, by="UNIT")
L2015<-sp::merge(UNITS,Lai2015, by="UNIT")
L2016<-sp::merge(UNITS,Lai2016, by="UNIT")Load the centroid points of each smapling unit (2.5 x 5 km) to extract the abundnce data.
From all the 2.5 x 5 km units from Laikipia, we selected those with low anthropogenic land cover.
Create a column with Cell ID (primary sampling units) and another column with Segment ID (secondary sampling units).
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\CregoRD\Documents\GitHub\CommunityOccupancy\Shapefiles", layer: "OccPoints2"
## with 576 features
## It has 8 fieldsExtract survey data for all species. Here we use the point.in.poly function to extract data on animal abundance from the sampling units. After, we save the attribute table of each year as new data frames.
dataset01<-point.in.poly(points, L2001)
data01<-dataset01@data
dataset04<-point.in.poly(points, L2004)
data04<-dataset04@data
dataset06<-point.in.poly(points, L2006)
data06<-dataset06@data
dataset08<-point.in.poly(points, L2008)
data08<-dataset08@data
dataset10<-point.in.poly(points, L2010)
data10<-dataset10@data
dataset12<-point.in.poly(points, L2012)
data12<-dataset12@data
dataset15<-point.in.poly(points, L2015)
data15<-dataset15@data
dataset16<-point.in.poly(points, L2016)
data16<-dataset16@dataCreate a loop to orginize the data in species, year, cell and segment.
#Create a vector with the name of the 15 species to work with
SpList<-c("EL","GF","ZB","ZG","TG","GG","KG","IM","BF","ED","OS","WH","OX","WB","GN")
##2001
det01<-data.frame()
for (i in 1:length(SpList)){
sp1<-data01[,SpList[i]]
sp2<-data.frame(Species=SpList[i], Det=sp1, Cell=data01$Cell, Rep=data01$Seg, Year=2001)
det01<-rbind(det01,sp2)
}
##2004
det04<-data.frame()
for (i in 1:length(SpList)){
sp1<-data04[,SpList[i]]
sp2<-data.frame(Species=SpList[i], Det=sp1, Cell=data04$Cell, Rep=data04$Seg, Year=2004)
det04<-rbind(det04,sp2)
}
##2006
det06<-data.frame()
for (i in 1:length(SpList)){
sp1<-data06[,SpList[i]]
sp2<-data.frame(Species=SpList[i], Det=sp1, Cell=data06$Cell, Rep=data06$Seg, Year=2006)
det06<-rbind(det06,sp2)
}
##2008
det08<-data.frame()
for (i in 1:length(SpList)){
sp1<-data08[,SpList[i]]
sp2<-data.frame(Species=SpList[i], Det=sp1, Cell=data08$Cell, Rep=data08$Seg, Year=2008)
det08<-rbind(det08,sp2)
}
##2010
det10<-data.frame()
for (i in 1:length(SpList)){
sp1<-data10[,SpList[i]]
sp2<-data.frame(Species=SpList[i], Det=sp1, Cell=data10$Cell, Rep=data10$Seg, Year=2010)
det10<-rbind(det10,sp2)
}
##2012
det12<-data.frame()
for (i in 1:length(SpList)){
sp1<-data12[,SpList[i]]
sp2<-data.frame(Species=SpList[i], Det=sp1, Cell=data12$Cell, Rep=data12$Seg, Year=2012)
det12<-rbind(det12,sp2)
}
##2015
det15<-data.frame()
for (i in 1:length(SpList)){
sp1<-data15[,SpList[i]]
sp2<-data.frame(Species=SpList[i], Det=sp1, Cell=data15$Cell, Rep=data15$Seg, Year=2015)
det15<-rbind(det15,sp2)
}
##2016
det16<-data.frame()
for (i in 1:length(SpList)){
sp1<-data16[,SpList[i]]
sp2<-data.frame(Species=SpList[i], Det=sp1, Cell=data16$Cell, Rep=data16$Seg, Year=2016)
det16<-rbind(det16,sp2)
}
#Combine datasets into one data frame.
det<-rbind(det01,det04,det06,det08,det10,det12,det15,det16)
#Extract the number of observed species.
n<-length(SpList)
n## [1] 15#Find the number of unique sampling locations and calculate the number of primary sampling units of cells.
ucells <- as.character(unique(points$Cell))
J<-length(ucells)
J## [1] 288#Find the number of unique secondary sampling units or segments.
useg <- as.character(unique(points$Seg))
I<-length(useg)
I## [1] 2## [1] 8Create a loop to transpose the segments into a data frame per year with the capture history for each cell. Check the structure and head of the new data frame.
detections<-data.frame()
for (y in 1:Y){
for (i in 1:length(SpList)){
dat0<-subset(det, Year==Years[y])
dat1<-subset(dat0, Species==SpList[i])
for (j in 1:J){
dat2<-subset(dat1, Cell==j)
dat3<-as.data.frame(t(dat2[-1]))
dat4<-data.frame(Species=SpList[i], Year=Years[y], X1=dat3[1,1],X2=dat3[1,2])
detections<-rbind(detections, dat4)
}
}
}
str(detections)## 'data.frame': 34560 obs. of 4 variables:
## $ Species: Factor w/ 15 levels "EL","GF","ZB",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Year : num 2001 2001 2001 2001 2001 ...
## $ X1 : num NA NA 0 0 0 0 0 0 0 0 ...
## $ X2 : num 0 0 0 0 0 0 0 0 0 0 ...Reshape the detection/non-detection data into a four-dimensional array DetHist, where the first dimension, J, is the cell; the second dimension, I, is the seg; the third dimension, n, is the species, and the fouth dimension, Y, is the year.
DetHist <- array(NA,dim=c(J,I,n,Y))
dimnames(DetHist)[[3]] <- unique(det$Species)
dimnames(DetHist)[[4]] <- unique(det$Year)
for (y in 1:dim(DetHist)[4]){
dat0<-subset(detections, Year==Years[y])
for (i in 1:dim(DetHist)[3]){
dat1<-subset(dat0, Species==SpList[i])
DetHist[,,i,y]<-as.matrix(dat1[,-(1:2)])
}
}Convert abundance data into presence/absence data.
DetHist <- replace(DetHist,which(DetHist>0),1)
DetHist[1:10,,1,1] #Print first 10 cells for elephants for 2001## [,1] [,2]
## [1,] NA 0
## [2,] NA 0
## [3,] 0 0
## [4,] 0 0
## [5,] 0 0
## [6,] 0 0
## [7,] 0 0
## [8,] 0 0
## [9,] 0 0
## [10,] 0 0This is the array to use for modeling.
Covariates for detectability
We hypothesized that detectability may be affected by group size and woody vegetation cover (Jachmann, 2002; Schlossberg et al., 2018). We obtained records of group size from the aerial surveys for years 2006, 2008 and 2010. We calculated the mean group size for each species to account for its potential effect on species detectability (Table 1).
Calculate group size
#2006 - Original data only had the final count.
GrpSize2006 <- read.csv("./Data/groupsize2006.csv", header = T)
GrpSize2006$UNIT<-as.factor(GrpSize2006$UNIT)
GrpSize2006$Sp<-as.factor(GrpSize2006$Sp)
#2008 - Original data has visual counts and counts on photographs (group size > 10 individuals).
GrpSize2008 <- read.csv("./Data/groupsize2008.csv", header = T)
GrpSize2008$UNIT<-as.factor(GrpSize2008$UNIT)
GrpSize2008$VCnt[is.na(GrpSize2008$VCnt)] <- 0 #Fill NA with zeros as the original csv files has no data for zeros.
GrpSize2008$PCnt[is.na(GrpSize2008$PCnt)] <- 0 #Fill NA with zeros as the original csv files has no data for zeros.
GrpSize2008<-within(GrpSize2008, CCnt<-ifelse(PCnt>0, PCnt, VCnt)) #Create a new count column that utilze the value of counts done on photographs when available.
#2010 - Original data only had the final count.
GrpSize2010 <- read.csv("./Data/groupsize2010.csv", header = T)
GrpSize2010$UNIT<-as.factor(GrpSize2010$UNIT)
#Extract only desire UNITS
SamplingUNITS<-unique(data04$UNIT)[-1] #-1 to eliminate the first value that is NA
StAreaGrpSize06<-subset(GrpSize2006, UNIT=SamplingUNITS)
StAreaGrpSize08<-subset(GrpSize2008, UNIT=SamplingUNITS)
StAreaGrpSize10<-subset(GrpSize2010, UNIT=SamplingUNITS)
StAreaGrpSize<-rbind(StAreaGrpSize06,StAreaGrpSize08,StAreaGrpSize10)
#Create a new dataframe with group size per species
GroupSize<-data.frame()
for (i in 1:length(SpList)){
x<-subset(StAreaGrpSize, Sp==SpList[i])
spmean<-mean(x$CCnt)
spmedian<-median(x$CCnt)
spSD<-sd(x$CCnt)
spMin<-min(x$CCnt)
spMax<-max(x$CCnt)
datfra<-data.frame(Species=SpList[i],MeanGrSize=spmean,MedianGS=spmedian,SDGS=spSD,MinGS=spMin,MaxGS=spMax,NumGrps=nrow(x))
GroupSize<-rbind(GroupSize,datfra)
}
GroupSize[1:2,]## Species MeanGrSize MedianGS SDGS MinGS MaxGS NumGrps
## 1 EL 6.056701 5 5.194844 1 29 194
## 2 GF 3.568627 2 3.098736 1 19 153Table 1. Group size summaries for 15 species of wild herbivores observed during aerial surveys for years 2006, 2008 and 2010 in Laikipia, Kenya.
| Species | Scientific name | Mean group size | Median group size | Group size SD | Min. group size | Max. group size |
|---|---|---|---|---|---|---|
| Elephant | Loxodonta africana | 6.06 | 5.00 | 5.19 | 1 | 29 |
| Reticulated giraffe | Giraffa reticulata | 3.57 | 2.00 | 3.10 | 1 | 19 |
| Plains zebra | Equus quagga | 11.18 | 8.00 | 11.44 | 1 | 116 |
| Grevy’s zebra | Equus grevyi | 4.33 | 3.50 | 3.50 | 1 | 16 |
| Thomson’s gazelle | Eudorcas thomsonii | 6.76 | 6.00 | 5.35 | 1 | 28 |
| Grant’s gazelle | Nanger (granti) notata | 5.11 | 4.00 | 4.04 | 1 | 24 |
| Hartebeest | Alcelaphus buselaphus lelwel | 6.53 | 5.00 | 6.82 | 1 | 34 |
| Impala | Aepyceros melampus | 9.14 | 6.00 | 10.43 | 1 | 65 |
| African buffalo | Syncerus caffer | 16.06 | 6.50 | 25.86 | 1 | 146 |
| Eland | Taurotragus oryx | 6.13 | 4.00 | 8.76 | 1 | 76 |
| Ostrich | Struthio camelus | 2.02 | 1.00 | 2.84 | 1 | 26 |
| Common warthog | Phacochoerus africanus | 3.62 | 3.00 | 2.40 | 1 | 11 |
| Beisa oryx | Oryx beisa | 7.14 | 2.00 | 10.12 | 1 | 43 |
| Defassa waterbuck | Kobus ellipsiprymnus defassa | 3.96 | 3.00 | 2.37 | 1 | 11 |
| Gerenuk | Litocranius walleri | 1.59 | 1.00 | 0.75 | 1 | 4 |
Additionally, we used data on above ground woody vegetation biomass (WVB) from 2010 to account for vegetation effects on animals’ detectability (Bouvet et al., 2018). To test the assumption that this dataset accurately represented woody vegetation cover in Laikipia, we randomly selected 100, 50 × 50 m cells across the study area. We then used Google Earth high-resolution imagery from 2009-2012 to visually estimate woody vegetation cover in each cell to the nearest 5%, achieving r = 0.82 (Fig. 3). As the exact path taken by the airplane during flights was not provided with the survey metadata, we estimated mean WVB from Bouvet et al., (2018) for the entire 2.5 × 5 km segment using zonal statistics in QGIS3.2 (QGIS Development Team, 2018). We assumed WVB to be consistent through time.
#Load raster dataset (clipped and reprojected to UTM 37N)
AGB<-raster("BushCover/AGB_50m_Laikipia_UTM37N.tif")
#Load random points (center of each 50 x 50 random cell)
randmpoints<-readOGR(dsn = "./BushCover", layer = "RandPointsForAGB")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\CregoRD\Documents\GitHub\CommunityOccupancy\BushCover", layer: "RandPointsForAGB"
## with 100 features
## It has 1 fields
## Integer64 fields read as strings: id#Load estiamated tree/bush cover
EstimatedCov<-read.csv("BushCover/RandPoints50x50ForAGB.csv")
#Extract AGB pixel value
EstimatedCov$AGB <- raster::extract(AGB, randmpoints)
#Create a new dataframe with estimated cover and AGB values
ForCorr<-data.frame(EstimatedCov$EstimatedCover09.12,EstimatedCov$AGB)
names(ForCorr)<-c("Est","Obs")
#Eliminate rows with NA values
ForCorr<-na.omit(ForCorr)
#Correlation
round(cor(ForCorr$Obs,ForCorr$Est, method="pearson"),2)## [1] 0.82##
## Call:
## lm(formula = ForCorr$Obs ~ ForCorr$Est)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.167 -5.215 -1.564 3.525 50.794
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.56435 1.80611 2.527 0.0131 *
## ForCorr$Est 0.68804 0.04943 13.919 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.6 on 96 degrees of freedom
## Multiple R-squared: 0.6687, Adjusted R-squared: 0.6652
## F-statistic: 193.7 on 1 and 96 DF, p-value: < 2.2e-16Figure. 3. Correlation between estimated bush/tree cover and woody vegetation biomass obtained from Bouvet et al. (2018).
We obtained a correlation of 0.82 (p < 0.01). We concluded that we can use AGB to estimate bush/tree cover in each segment sampled.
As the exact path taken by the airplane during flights is not provided with survey metadata, we estimated mean WVB from Bouvet et al., (2018) for the entire 2.5 x 5 km segment using zonal statistics in QGIS3.2 (QGIS Development Team, 2018). We assumed WVB to be consistent across time periods.
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\CregoRD\Documents\GitHub\CommunityOccupancy\Shapefiles", layer: "SegmForOccu"
## with 730 features
## It has 6 fieldssegcov2<-point.in.poly(points, segcov)
segcov3<-segcov2@data
segcov3<-data.frame(Cell=points$Cell,VegCov=segcov3$AGBmean)
segcov3<-round(segcov3,4)
#Arrange woody vegetation biomass in a matrix of 288 x 2
segcov4<-data.frame()
for (j in 1:J){
dat1<-subset(segcov3, Cell==j)
dat2<-as.data.frame(t(dat1[-1]))
dat3<-data.frame(X1=dat2[1,1],X2=dat2[1,2])
segcov4<-rbind(segcov4, dat3)
}
segcovar<-as.matrix(round(segcov4,4))
head(segcovar)## X1 X2
## [1,] 12.7990 4.2489
## [2,] 4.2556 7.7199
## [3,] 67.5055 71.9713
## [4,] 53.8961 66.3287
## [5,] 40.0471 33.0021
## [6,] 26.0596 26.7305Covariates for occupancy
We incorporated distance to water as a covariate to account for its known effects on herbivore distribution (Ogutu et al., 2014b, 2010; Redfern et al., 2003; Sitters et al., 2009; Tyrrell et al., 2017). For each cell, we estimated the Euclidean distance from the center of the cell to the closest water source, including permanent rivers (obtained from the World Resources Institute [datasets.wri.org] and corrected using Google Earth Imagery) and permanent artificial dams. To incorporate dams, we used a global surface water layer (Pekel et al., 2016), retaining pixels containing water for >10 months per year between 1984 and 2015. We accounted for forage availability effects on herbivore distribution by incorporating the Normalized Difference Vegetation Index (NDVI) from the 1 km MODIS/Terra data product (MYD13A2). The NDVI has been proven to be a good predictor for distribution and abundance of herbivore species (Pettorelli et al., 2011; Sitters et al., 2009; Tyrrell et al., 2017). For each surveyed year, we obtained all images between January 1 (i.e., beginning of the dry season) and the last day of the survey and calculated the median value for each pixel. For each year, we then calculated the median value of NDVI for each cell using zonal statistics.
#Load shapefile with cell covariate values
cellcov<-readOGR(dsn = "./Shapefiles", layer = "CellsForOccu")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\CregoRD\Documents\GitHub\CommunityOccupancy\Shapefiles", layer: "CellsForOccu"
## with 365 features
## It has 22 fields
## Integer64 fields read as strings: Area#Create a new point shapefile that is at the center of each cell
segppointcoord<-data.frame(points@coords)
#Estimate the mean coordinate value
newcoord<-data.frame()
for (i in 1:length(segppointcoord$coords.x1)){
dat1<-(segppointcoord[i,]+segppointcoord[i+1,])/2
newcoord<-rbind(newcoord, dat1)
}
#Remove the cells that are not needed
toDelete <- seq(2, length(segppointcoord$coords.x1), 2)
newcoord2 <- newcoord[-toDelete, 1:2]
#Create the new shapefile
coordinates(newcoord2)<-newcoord2
proj4string(newcoord2) <- CRS("+proj=utm +zone=37 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
#Add cell id
newcoord2@data$Cell<-seq(1,J,by=1)
#Extract covariates
Pointcellcov<-point.in.poly(newcoord2, cellcov)
#Estimate shortest distance from center of cell to water
#Load permanent river shapefile
PermanRivers<-readOGR(dsn = "./Shapefiles", layer = "MayorRivers&DamsUTM37N")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\CregoRD\Documents\GitHub\CommunityOccupancy\Shapefiles", layer: "MayorRivers&DamsUTM37N"
## with 210 features
## It has 5 fieldsPermanRivers<-spTransform(PermanRivers,CRS("+proj=longlat +ellps=WGS84"))
coordforriv<-spTransform(newcoord2,CRS("+proj=longlat +ellps=WGS84"))
#Calulate shortest distance to river
DistWater<-dist2Line(coordforriv, PermanRivers)
#Add covariate to a data frame
data<-Pointcellcov@data
data$DistWater<-DistWater[,1]Google Earth Engine code to process and download Modis NDVI information
var aoi = ee.Geometry.Polygon(
[[[36.1, 0.95],
[36.1, -0.5],
[37.7, -0.5],
[37.7, 0.95]]])
Map.centerObject(aoi)
var modis = ee.ImageCollection("MODIS/006/MOD13A2");
var modis01 = ee.ImageCollection(modis)
.filterDate('2001-01-01', '2001-2-23')
.select('NDVI');
var modisclip01 = modis01.median().clip(aoi);
var modisclip01 = modisclip01.multiply(0.0001);
var modis04 = ee.ImageCollection(modis)
.filterDate('2004-01-01', '2004-3-1')
.select('NDVI');
var modisclip04 = modis04.median().clip(aoi);
var modisclip04 = modisclip04.multiply(0.0001);
var modis06 = ee.ImageCollection(modis)
.filterDate('2006-01-01', '2006-2-22')
.select('NDVI');
var modisclip06 = modis06.median().clip(aoi);
var modisclip06 = modisclip06.multiply(0.0001);
var modis08 = ee.ImageCollection(modis)
.filterDate('2008-01-01', '2008-3-10')
.select('NDVI');
var modisclip08 = modis08.median().clip(aoi);
var modisclip08 = modisclip08.multiply(0.0001);
var modis10 = ee.ImageCollection(modis)
.filterDate('2010-01-01', '2010-2-26')
.select('NDVI');
var modisclip10 = modis10.median().clip(aoi);
var modisclip10 = modisclip10.multiply(0.0001);
var modis12 = ee.ImageCollection(modis)
.filterDate('2012-01-01', '2012-3-09')
.select('NDVI');
var modisclip12 = modis12.median().clip(aoi);
var modisclip12 = modisclip12.multiply(0.0001);
var modis15 = ee.ImageCollection(modis)
.filterDate('2015-01-01', '2015-3-12')
.select('NDVI');
var modisclip15 = modis15.median().clip(aoi);
var modisclip15 = modisclip15.multiply(0.0001);
var modis16 = ee.ImageCollection(modis)
.filterDate('2016-01-01', '2016-4-19')
.select('NDVI');
var modisclip16 = modis16.median().clip(aoi);
var modisclip16 = modisclip16.multiply(0.0001);
Map.addLayer(modisclip16, {min: 0, max: 1, palette: ['FFFFFF','CC9966','CC9900',
'996600', '33CC00', '009900','006600','000000']}, 'NDVI');
Export.image.toDrive({
image: modisclip01,
description: "Laikipia01",
scale: 1000,
folder:"GIS",
maxPixels: 50100100100
});
Export.image.toDrive({
image: modisclip04,
description: "Laikipia04",
scale: 1000,
folder:"GIS",
maxPixels: 50100100100
});
Export.image.toDrive({
image: modisclip06,
description: "Laikipia06",
scale: 1000,
folder:"GIS",
maxPixels: 50100100100
});
Export.image.toDrive({
image: modisclip08,
description: "Laikipia08",
scale: 1000,
folder:"GIS",
maxPixels: 50100100100
});
Export.image.toDrive({
image: modisclip10,
description: "Laikipia10",
scale: 1000,
folder:"GIS",
maxPixels: 50100100100
});
Export.image.toDrive({
image: modisclip12,
description: "Laikipia12",
scale: 1000,
folder:"GIS",
maxPixels: 50100100100
});
Export.image.toDrive({
image: modisclip15,
description: "Laikipia15",
scale: 1000,
folder:"GIS",
maxPixels: 50100100100
});
Export.image.toDrive({
image: modisclip16,
description: "Laikipia16",
scale: 1000,
folder:"GIS",
maxPixels: 50100100100
});We incorporated livestock abundance effects on wild herbivore occurrence (Georgiadis et al., 2007; Keesing et al., 2018; Ogutu et al., 2014a; Sitters et al., 2009) by calculating the total counts of cattle, sheep, goats, camels and donkeys for each cell from the aerial survey dataset as a measure of relative abundance. We did not account for detectability in livestock counts because livestock occur in large highly detectable herds. Additionally, herds >10 animals are confirmed by photographs. We used single-season occupancy models to confirm that livestock detectability was close to one and that it was not affected by WVB. We assumed that livestock counts are a good measure of livestock relative abundance across the study area and that there is no spatial variation given by errors in counting animals that could obscure effects when using livestock counts as a covariate for herbivore occurrence.
#Livestock abundance (cattle, shoats, cammels, donkeys)
SpListDom<-c("CX","SH","DN","CM")
##2001
detdom01<-data.frame()
for (i in 1:length(SpListDom)){
sp1<-data01[,SpListDom[i]]
sp2<-data.frame(Species=SpListDom[i], Det=sp1, Cell=data01$Cell, Rep=data01$Seg, Year=2001)
detdom01<-rbind(detdom01,sp2)
}
##2004
detdom04<-data.frame()
for (i in 1:length(SpListDom)){
sp1<-data04[,SpListDom[i]]
sp2<-data.frame(Species=SpListDom[i], Det=sp1, Cell=data04$Cell, Rep=data04$Seg, Year=2004)
detdom04<-rbind(detdom04,sp2)
}
##2006
detdom06<-data.frame()
for (i in 1:length(SpListDom)){
sp1<-data06[,SpListDom[i]]
sp2<-data.frame(Species=SpListDom[i], Det=sp1, Cell=data06$Cell, Rep=data06$Seg, Year=2006)
detdom06<-rbind(detdom06,sp2)
}
##2008
detdom08<-data.frame()
for (i in 1:length(SpListDom)){
sp1<-data08[,SpListDom[i]]
sp2<-data.frame(Species=SpListDom[i], Det=sp1, Cell=data08$Cell, Rep=data08$Seg, Year=2008)
detdom08<-rbind(detdom08,sp2)
}
##2010
detdom10<-data.frame()
for (i in 1:length(SpListDom)){
sp1<-data10[,SpListDom[i]]
sp2<-data.frame(Species=SpListDom[i], Det=sp1, Cell=data10$Cell, Rep=data10$Seg, Year=2010)
detdom10<-rbind(detdom10,sp2)
}
##2012
detdom12<-data.frame()
for (i in 1:length(SpListDom)){
sp1<-data12[,SpListDom[i]]
sp2<-data.frame(Species=SpListDom[i], Det=sp1, Cell=data12$Cell, Rep=data12$Seg, Year=2012)
detdom12<-rbind(detdom12,sp2)
}
##2015
detdom15<-data.frame()
for (i in 1:length(SpListDom)){
sp1<-data15[,SpListDom[i]]
sp2<-data.frame(Species=SpListDom[i], Det=sp1, Cell=data15$Cell, Rep=data15$Seg, Year=2015)
detdom15<-rbind(detdom15,sp2)
}
##2016
detdom16<-data.frame()
for (i in 1:length(SpListDom)){
sp1<-data16[,SpListDom[i]]
sp2<-data.frame(Species=SpListDom[i], Det=sp1, Cell=data16$Cell, Rep=data16$Seg, Year=2016)
detdom16<-rbind(detdom16,sp2)
}
detdom<-rbind(detdom01,detdom04,detdom06,detdom08,detdom10,detdom12,detdom15,detdom16)
## Transpose the segments into a dataframa with capture history for each cell
detectionsDom<-data.frame()
for (y in 1:Y){
for (i in 1:length(SpListDom)){
dat0<-subset(detdom, Year==Years[y])
dat1<-subset(dat0, Species==SpListDom[i])
for (j in 1:J){
dat2<-subset(dat1, Cell==j)
dat3<-as.data.frame(t(dat2[-1]))
dat4<-data.frame(Species=SpListDom[i], Year=Years[y], X1=dat3[1,1],X2=dat3[1,2])
detectionsDom<-rbind(detectionsDom, dat4)
}
}
}
#Estimate total abundance per cell and per year.
detectionsDom$sum<-apply(detectionsDom[,3:4], 1, sum, na.rm=T)
#Sum the abundance of four species per year
AbunDom <- matrix(NA, nrow = J, ncol = Y)
dat2<-matrix(NA, nrow = J, ncol = length(SpListDom))
for (y in 1:Y) {
dat0<-subset(detectionsDom, Year==Years[y])
for (i in 1:length(SpListDom)){
dat1<-subset(dat0, Species==SpListDom[i])
dat2[,i]<-dat1$sum
dat3<-apply(dat2, 1, sum)
}
AbunDom[,y]<-dat3
}
AbunDom<-as.data.frame(AbunDom)
colnames(AbunDom)<-c("2001","2004","2006","2008","2010","2012","2015","2016")
AbunDom<-replace(AbunDom,is.na(AbunDom),0)
head(AbunDom)## 2001 2004 2006 2008 2010 2012 2015 2016
## 1 94 396 38 512 252 44 307 266
## 2 0 22 223 148 151 41 274 338
## 3 96 24 149 190 0 26 263 27
## 4 44 85 370 161 0 0 77 198
## 5 340 27 151 0 0 0 216 0
## 6 129 62 22 0 50 0 86 0To test this assumption, we estimated livestock detectability using a single-season occupancy model (Royle and Dorazio, 2008). From the aerial survey, we obtained the presence or absence of livestock (cattle, sheep, goats, camels or donkeys) per segment. We used year as a random effect. We allowed detectability to vary as a function of above ground woody vegetation biomass (WVB) obtained from Bouvet et al. (2018). Because we found significant autocorrelation in the models (Moran’s I test, p-value < 0.01) and autocorrelograms based on 999 simulations showed evidence of high autocorrelation for all years (Fig. 1C), we incorporated a spatial auto-covariate structure in the model (Royle and Dorazio, 2008). We used a 4 neighborhood structure for the auto-covariate specification.
We implemented models using program OpenBUGS and JAGS, through packages R2OpenBUGS and jagsUI in R programing language (R Development Core Team, 2016). We ran 3 parallel Markov chain Monte Carlo (MCMC) chains to find 30,000 samples of the parameters of interest after a 10,000 burn-in period and a 1:10 thinning rate. Non-informative priors were used in all cases.
Estimate livestock detection probability
#Create detection history for livestock
detectionsDom<-data.frame()
for (y in 1:Y){
for (i in 1:length(SpListDom)){
dat0<-subset(detdom, Year==Years[y])
dat1<-subset(dat0, Species==SpListDom[i])
for (j in 1:J){
dat2<-subset(dat1, Cell==j)
dat3<-as.data.frame(t(dat2[-1]))
dat4<-data.frame(Species=SpListDom[i], Year=Years[y], X1=dat3[1,1],X2=dat3[1,2])
detectionsDom<-rbind(detectionsDom, dat4)
}
}
}
DomDetHist <- array(NA,dim=c(J,I,4,Y))
dimnames(DomDetHist)[[3]] <- unique(detdom$Species)
dimnames(DomDetHist)[[4]] <- unique(detdom$Year)
for (y in 1:dim(DomDetHist)[4]){
dat0<-subset(detectionsDom, Year==Years[y])
for (i in 1:dim(DomDetHist)[3]){
dat1<-subset(dat0, Species==SpListDom[i])
DomDetHist[,,i,y]<-as.matrix(dat1[,-(1:2)])
}
}
##Combine four species
domY<-array(NA, dim=c(J,I,Y))
dimnames(domY)[[3]] <- unique(detdom$Year)
for (y in 1:8){
Dom1<-DomDetHist[,,,y]
Dom2<-matrix(mapply(sum, Dom1[,,1], Dom1[,,2], Dom1[,,3], Dom1[,,4], MoreArgs=list(na.rm=T)),ncol=2)
domY[,,y]<-replace(Dom2,which(Dom2>0),1)
}
#Covariable for detectability
COV<-as.matrix(segcovar)
sd.COV<-sd(c(COV),na.rm=T)
mean.COV<-mean(c(COV),na.rm=T)
COVs <- (COV - mean.COV) / sd.COV #Standarize covariate
#############################################################
#############################################################
# Compute the neighborhood (2nd order)
#We suppose the neighbors are the quadrats immediately to the north, south,
#east and west, so a site has a maximum of 4 neighbors.
neigh <- dnearneigh(newcoord2, d1 = 0, d2 = 5000)
# Store the neighborhood in the format needed by the BUGS model
cardneigh <- card(neigh)
adjmat <- matrix(nrow = length(neigh), ncol = max(cardneigh))
for (i in 1:length(neigh)) {
adjmat[i, 1:cardneigh[i]] <- neigh[[i]]
}
#Calculte number of neighbors
nn <- replace(adjmat,which(adjmat>0),1)
numnn<-apply(nn, 1, sum,na.rm=TRUE)###################################
## Single-season occupancy model
cat("
model {
## Model
## Priors
mean.psi~ dunif(0,1)
sd.psi ~ dunif(0,8)
tau.psi <- pow(sd.psi,-2)
mean.p~ dunif(0,1)
sd.p ~ dunif(0,8)
tau.p <- pow(sd.p,-2)
for (t in 1:Y){
psi0[t] ~ dnorm(mean.psi,tau.psi)
p0[t] ~ dnorm(mean.p,tau.p)
}
pcov ~ dnorm(0, 0.01)
#State process
for(i in 1:nsite) {
for (t in 1:Y) {
z[i,t] ~ dbern(psi[i,t])
logit(psi[i,t]) <- psi0[t]
# Detection process
for(j in 1:nrep){
logit(p[i,j,t]) <- p0[t] + pcov * COV[i,j]
mu.p[i,j,t] <- p[i,j,t]*z[i,t]
y[i,j,t] ~ dbern(mu.p[i,j,t])
}
yres[i,t] <- psi[i,t]-z[i,t]
}
}
}
",file="occmodeltime.txt")
occDatat <- list(
y = domY[,,],
COV=COV,
nsite=J,
nrep=I,
Y=Y
)
OccPar <- c("mean.p","psi0","p0","pcov","yres")
zst<-array(1,dim=c(J,Y))
OccInit <- function() {
list(psi0 = runif(8), p0 = runif(8), pcov = runif(1), z = zst)
}
Occt <- jags(occDatat, OccInit, OccPar,"occmodeltime.txt", n.chains = 3, n.iter = 30000, n.burnin = 10000, n.thin=10,
parallel = T)
######Moran's I
#Compute Moran I test
#
lw <- nb2listw(neigh, style="B")
#Run the Monte Carlo simulation with 1000 simulations.
MoranI<-matrix(NA,8)
for (t in 1:8){
MI<-moran.mc(as.vector(Occt$mean$yres[,t]),lw,nsim=999,zero.policy = NULL)
p<-MI$p.value
MoranI[t]<-p
}
#P > 0.05 i.e., no evidence of spatial autocorrelation
#Plot correlograms
par(mfrow=c(4,1), mar=c(4,4,2,2))
for(i in 1:8) {
bayres<-correlogram(newcoord2, v = as.vector(Occt$mean$yres[,i]), dist =5000, ns=999, latlong = F)
}
#Create autologistic model
cat("
model{
## Priors
mean.psi~ dunif(0,1)
sd.psi ~ dunif(0,8)
tau.psi <- pow(sd.psi,-2)
mean.p~ dunif(0,1)
sd.p ~ dunif(0,8)
tau.p <- pow(sd.p,-2)
for (t in 1:Y){
psi0[t] ~ dnorm(mean.psi,tau.psi)
p0[t] ~ dnorm(mean.p,tau.p)
}
beta~dnorm(0,0.01)
pcov~dnorm(0,0.01)
## Model
# State process
for(i in 1:nsite) {
for (t in 1:Y) {
z[i,t] ~ dbern(psi[i,t])
x[i,1,t]<-0
for(k in 1:numnn[i]){
x[i,k+1,t]<-x[i,k,t]+z[adjmat[i,k],t]
}#k
autocov[i,t]<-(x[i,(numnn[i]+1),t]/numnn[i])
logit(psi[i,t])<- psi0[t] + beta*autocov[i,t]
# Detection process
for(j in 1:nrep){
logit(p[i,j,t]) <- p0[t] + pcov * COV[i,j]
p.eff[i,j,t] <- p[i,j,t]*z[i,t]
y[i,j,t] ~ dbern(p.eff[i,j,t])
}#j
yres[i,t] <- psi[i,t]-z[i,t]
}#t
}#i
}
",file="AUToccmodeltime.txt")
##Bugs Data
#For all
OccParaut <- c("mean.p","psi0","p0","pcov","beta","autocov","yres")
zst<-array(1,dim=c(J,Y))
OccInitaut <- function() {
list(psi0 = runif(8), p0 = runif(8), pcov = runif(1), beta = runif(1), z = zst)
}
OccDatataut <- list(
y = domY[,,],
COV=COV,
nsite=J,
nrep=I,
Y=Y,
numnn=numnn,
adjmat=adjmat
)
M <- bugs(OccDatataut, OccInitaut, OccParaut,"AUToccmodeltime.txt", n.chains = 3,n.iter = 30000, n.burnin = 10000, n.thin=10,
working.directory = getwd(), debug = T)
######Moran's I
#Compute Moran I test
#
lw <- nb2listw(neigh, style="B")
#Run the Monte Carlo simulation with 1000 simulations.
MoranIauc<-matrix(NA,8)
for (t in 1:8){
MI<-moran.mc(as.vector(M$mean$yres[,t]),lw,nsim=999,zero.policy = NULL)
p<-MI$p.value
MoranIauc[t]<-p
}
#P > 0.05 i.e., no evidence of spatial autocorrelation
#Plot correlograms
par(mfrow=c(4,1), mar=c(4,4,2,2))
for(i in 1:8) {
bayres<-correlogram(newcoord2, v = as.vector(M$mean$yres[,i]), dist =5000, ns=999, latlong = F)
}## [1] 0.693## [1] 0.6592374## [1] 0.7206927Model residuals did not show signs of autocorrelation after incorporating the auto-covariate structure (Fig. 4). Mean detection probability for livestock was high (0.69; [95% credible intervals: 0.66 – 0.72]) and not related to woody vegetation biomass (-0.05 [-0.13-0.04]; Fig. 5). Therefore, the probability of detecting at least one livestock in a cell with the two sampling replicates was 〖1-(1-0.69)〗^2 = 0.9039. We concluded that detectability of livestock is high regardless of woody vegetation density. We assumed that livestock counts are a good measure of livestock relative abundance across the study area and that there is no spatial variation given by errors in counting animals that could obscure effects when using livestock counts as a covariate for herbivore species occurrence.
Figure 4. Correlograms showing the spatial correlation for model residuals as a function of distance (m) for occupancy models with (left) and without (right) a spatial auto-covariate specification.
Figure 5. Probability of detection (95% Bayesian Credible Interval) as a function of above ground biomass of woody vegetation for livestock on Laikipia rangelands, Kenya.
P=seq(1,70, by=1)
Pst <- (P - mean.COV) / sd.COV
p.pred.mean <- plogis(mean(M$sims.list$mean.p) + mean(M$sims.list$pcov) * Pst)
M01mc<-cbind(M$sims.list$mean.p, M$sims.list$pcov)
p.predS_mean_dist <- matrix(NA, nrow = nrow(M01mc), ncol = 70)
for (i in 1:nrow(p.predS_mean_dist)){
p.predS_mean_dist[i,] <- plogis(M01mc[i,1] + M01mc[i,2] * Pst)
}
credible_lowerS <- apply(p.predS_mean_dist, MARGIN = 2, quantile, prob = 0.025)
credible_upperS <- apply(p.predS_mean_dist, MARGIN = 2, quantile, prob = 0.975)
##Code for fig. 5
plot(P, p.pred.mean, type="n", ylim=range(c(0.5,0.8)),
ylab="Probability of detection (95% CRI)", xlab="Woody vegetation biomass", cex.lab=0.9, cex.axis=0.8)
polygon(c(P, rev(P)), c(credible_lowerS, rev(credible_upperS)), col = "gray", border = NA)
lines(P, plogis(mean(M$sims.list$p0[,1]) + mean(M$sims.list$pcov) * Pst), lty=2, col=34)
lines(P, plogis(mean(M$sims.list$p0[,2]) + mean(M$sims.list$pcov) * Pst), lty=3, col=142)
lines(P, plogis(mean(M$sims.list$p0[,3]) + mean(M$sims.list$pcov) * Pst), lty=4, col=143)
lines(P, plogis(mean(M$sims.list$p0[,4]) + mean(M$sims.list$pcov) * Pst), lty=5, col=259)
lines(P, plogis(mean(M$sims.list$p0[,5]) + mean(M$sims.list$pcov) * Pst), lty=6, col=4)
lines(P, plogis(mean(M$sims.list$p0[,6]) + mean(M$sims.list$pcov) * Pst), lty=2, col=497)
lines(P, plogis(mean(M$sims.list$p0[,7]) + mean(M$sims.list$pcov) * Pst), lty=3, col=109)
lines(P, plogis(mean(M$sims.list$p0[,8]) + mean(M$sims.list$pcov) * Pst), lty=4, col=50)
lines(P, p.pred.mean, lwd=2)
legend(0.8,0.65, legend = c("Mean",2001,2004,2006,2008,2010,2012,2015,2016), lty = c(1,2,3,4,5,6,2,3,4),
col = c(1,34,142,143,259,4,497,109,50), lwd = c(2,1,1,1,1,1,1,1,1), cex=0.85)
Model specifications
Multispecies site-occupancy models can be formulated as a hierarchical state-space model that links two binary regression models, one model for the occupancy process of each species and a second model for the observation process conditional on occupancy (Kery and Royle, 2016; Royle and Dorazio, 2008; Zipkin et al., 2009). We considered the history of sightings obtained from the aerial survey for i = 1, 2, …, 15 species at j = 1, 2, …, 288 sites (cells), for the spatial replicates (segments) k = 1 and 2 and over t = 1, 2, …, 8 years. The occupancy status of cells can be modeled as the outcome of a Bernoulli distribution as z(j,i,t) ~ Bern (ψ(j,i,t)) where ψ(j,i,t) is the probability that species i is present at site j in year t. The state variable z(j,i,t) is then conditional to the observation process x(j,k,i,t) for species i at site j for the spatial replicate k and year t, which is also assumed to follow a Bernoulli distribution as x(j,k,i,t) ~ Bern (p(j,k,i,t)* z (j,i,t)) where p(j,k,i,t) is the detection probability of for species i and spatial replicate k for year t, if the species is present at site i. In this model, detectability is zero when the species does not occur on a specific site (i.e., z (j,i,t)=0) (Kery and Royle, 2016; Royle and Dorazio, 2008; Zipkin et al., 2009).
As several species were rarely observed, estimating all parameters would not be possible in a species-specific analysis. For this reason, we used the multispecies occupancy approach that incorporates community hierarchical components into the model (Royle and Dorazio, 2008; Zipkin et al., 2009). We did not incorporate unobserved species into the estimations, important for assessing total species richness in a study system, because it was beyond the scope of our study. In the specification of the model, species-level parameters (intercepts u and v and α and β coefficients for each covariate on occupancy and detection probability respectively) are treated as random effects. We allowed occupancy and detection probability to be influenced by covariates that were incorporated into the model using the logit-link function (Royle and Dorazio, 2008). Following Goijman et al., (2015), we incorporated random time effects on the species-specific intercepts (u and v) to control for potential variation across years due to climatic conditions that could affect species occurrence in subsequent years and different observers during counts across yearly surveys.
Across our study area, there is an association between covariates and the different land uses categories (Table 2; Fig. 6). We tested for a potential association between continuous covariates (distance to water, NDVI and livestock relative abundance) and land uses using a non-parametric analysis of variance, the Kruskal-Wallis test. The null hypothesis is that the continuous covariate and the land uses are not associated. A significant Kruskal-Wallis test indicates that there is an association. We found that all covariates are associated with the land-use categories (Table D1; Fig. D1). This is because wildlife and ranching and wildlife areas have fewer stocking rates than ranching and pastoralist areas. Vegetation productivity, measured as NDVI, is generally higher at wildlife and ranching areas than at wildlife and ranching or pastoralist areas, and the same is true for distance to water. Wildlife and ranching areas have more artificial dams than wildlife and ranching or pastoralist areas which result in areas in higher distances to water for cells in many wildlife and ranching or pastoralist sites. Moreover, several cells (~35%) in our study area include more than one type of land use for which it will be wrong to assign a specific land use type. For these reasons, we modeled occupancy probability as a function of distance to water, NDVI and their quadratic effects and livestock relative abundance, a set of covariates that may be then acting as the proximate causes of wildlife occurrence across the different land uses.
Figure 6. Boxplots showing the relationship between (a) distance to water (m), (b) NDVI, and (c) livestock abundance at four different land uses (Pastoralist, ranching, ranching and wildlife, and wildlife only) for all sites sampled in Laikipia, Kenya.
boxplotCovariates<-rbind(data.frame(LandUse=data$LU100, NDVI=data$NDVI01me_1, LA= AbunDom$`2001`, Year=2001), data.frame(LandUse=data$LU100, NDVI=data$NDVI04me_1, LA= AbunDom$`2004`, Year=2004), data.frame(LandUse=data$LU100, NDVI=data$NDVI06me_1, LA= AbunDom$`2006`, Year=2006), data.frame(LandUse=data$LU100, NDVI=data$NDVI08me_1, LA= AbunDom$`2008`, Year=2008), data.frame(LandUse=data$LU100, NDVI=data$NDVI10me_1, LA= AbunDom$`2010`, Year=2010), data.frame(LandUse=data$LU100, NDVI=data$NDVI12me_1, LA= AbunDom$`2012`, Year=2012), data.frame(LandUse=data$LU100, NDVI=data$NDVI15me_1, LA= AbunDom$`2015`, Year=2015), data.frame(LandUse=data$LU100, NDVI=data$NDVI16me_1, LA= AbunDom$`2016`, Year=2016))
boxplotCovariates<-subset(boxplotCovariates, !LandUse=="MixedFarmi")
boxplotCovariates<-subset(boxplotCovariates, !LandUse=="NA")
boxplotCovariates$LandUse<-factor(boxplotCovariates$LandUse)
boxplotCovariates$Year<-factor(boxplotCovariates$Year)
levels(boxplotCovariates$LandUse) <-c("Pastoralist", "Ranching","Ranching and wildlife","Wildlife only")
boxplot1 <- ggplot(boxplotCovariates, aes(x=LandUse, y=LA)) + geom_boxplot() + facet_wrap (~Year, ncol=4)+
labs(x = "Land use", y = "Livestock abundance \n(number of individuals)") + theme(axis.text.x = element_text(angle=45, hjust=1))+
theme(axis.text = element_text(size = 10), axis.title = element_text(size = 14), strip.text = element_text(size=12))
boxplot2 <- ggplot(boxplotCovariates, aes(x=LandUse, y=NDVI)) + geom_boxplot() + facet_wrap (~Year, ncol=4)+
labs(x = "Land use", y = "NDVI") + theme(axis.text.x = element_text(angle=45, hjust=1))+
theme(axis.text = element_text(size = 10), axis.title = element_text(size = 14), strip.text = element_text(size=12))
boxplotCovariates2<-data.frame(LandUse=data$LU100, DW = data$DistWater)
boxplotCovariates2<-subset(boxplotCovariates2, !LandUse=="MixedFarmi")
boxplotCovariates2<-subset(boxplotCovariates2, !LandUse=="NA")
boxplotCovariates2$LandUse<-factor(boxplotCovariates2$LandUse)
levels(boxplotCovariates2$LandUse) <-c("Pastoralist", "Ranching","Ranching and wildlife","Wildlife only")
boxplot3 <- ggplot(boxplotCovariates2, aes(x=LandUse, y=DW)) + geom_boxplot()+
labs(x = "Land use", y = "Distance to water") + theme(axis.text.x = element_text(angle=45, hjust=1))+
theme(axis.text = element_text(size = 10), axis.title = element_text(size=14))
#Plot figure
grid.arrange(boxplot3 ,boxplot2 , boxplot1 , ncol = 2)
## Standarize covariates
#Distance to Water
sd.DW<-sd(data$DistWater,na.rm=T)
mean.DW<-mean(data$DistWater,na.rm=T)
DW <- (data$DistWater - mean.DW) / sd.DW
#NDVI
sd.NDVI<-sd(as.matrix(data[,15:22]),na.rm=T)
mean.NDVI<-mean(as.matrix(data[,15:22]),na.rm=T)
NDVI01 <- (data$NDVI01me_1 - mean.NDVI) / sd.NDVI
NDVI04 <- (data$NDVI04me_1 - mean.NDVI) / sd.NDVI
NDVI06 <- (data$NDVI06me_1 - mean.NDVI) / sd.NDVI
NDVI08 <- (data$NDVI08me_1 - mean.NDVI) / sd.NDVI
NDVI10 <- (data$NDVI10me_1 - mean.NDVI) / sd.NDVI
NDVI12 <- (data$NDVI12me_1 - mean.NDVI) / sd.NDVI
NDVI15 <- (data$NDVI15me_1 - mean.NDVI) / sd.NDVI
NDVI16 <- (data$NDVI16me_1 - mean.NDVI) / sd.NDVI
NDVI <- as.matrix(data.frame(NDVI01,NDVI04,NDVI06,NDVI08,NDVI10,NDVI12,NDVI15,NDVI16))
#Dom abundance
sd.DA<-sd(as.matrix(AbunDom),na.rm=T)
mean.DA<-mean(as.matrix(AbunDom),na.rm=T)
DA01 <- (AbunDom$`2001` - mean.DA) / sd.DA
DA04 <- (AbunDom$`2004` - mean.DA) / sd.DA
DA06 <- (AbunDom$`2006` - mean.DA) / sd.DA
DA08 <- (AbunDom$`2008` - mean.DA) / sd.DA
DA10 <- (AbunDom$`2010` - mean.DA) / sd.DA
DA12 <- (AbunDom$`2012` - mean.DA) / sd.DA
DA15 <- (AbunDom$`2015` - mean.DA) / sd.DA
DA16 <- (AbunDom$`2016` - mean.DA) / sd.DA
DA <- as.matrix(data.frame(DA01,DA04,DA06,DA08,DA10,DA12,DA15,DA16))
#Covariable for detectability
sd.COV<-sd(c(segcovar),na.rm=T)
mean.COV<-mean(c(segcovar),na.rm=T)
COV <- (segcovar - mean.COV) / sd.COV
#Check correlation among covariates for each year
forcor01<- data.frame(DW = DW, NDVI=NDVI01, DA01)
cor(forcor01, method ="spearman")## DW NDVI DA01
## DW 1.00000000 0.06943354 0.05889265
## NDVI 0.06943354 1.00000000 0.12149073
## DA01 0.05889265 0.12149073 1.00000000## DW NDVI DA04
## DW 1.000000000 -0.04786046 -0.002685415
## NDVI -0.047860464 1.00000000 -0.199985911
## DA04 -0.002685415 -0.19998591 1.000000000## DW NDVI DA06
## DW 1.00000000 -0.086771481 -0.040799201
## NDVI -0.08677148 1.000000000 0.009371144
## DA06 -0.04079920 0.009371144 1.000000000## DW NDVI DA08
## DW 1.0000000 0.1152529 0.1030167
## NDVI 0.1152529 1.0000000 0.1243975
## DA08 0.1030167 0.1243975 1.0000000## DW NDVI DA10
## DW 1.00000000 -0.1751294 -0.07141242
## NDVI -0.17512942 1.0000000 0.10680872
## DA10 -0.07141242 0.1068087 1.00000000## DW NDVI DA12
## DW 1.000000000 -0.01353994 -0.007107288
## NDVI -0.013539936 1.00000000 -0.238334491
## DA12 -0.007107288 -0.23833449 1.000000000## DW NDVI DA15
## DW 1.00000000 -0.1243988 0.09277324
## NDVI -0.12439879 1.0000000 -0.14078572
## DA15 0.09277324 -0.1407857 1.00000000## DW NDVI DA16
## DW 1.0000000 -0.1300181 -0.0716818
## NDVI -0.1300181 1.0000000 -0.1811002
## DA16 -0.0716818 -0.1811002 1.0000000#####################
#Plot covariate for the 8 years
# Store the results in a SpatialPixelsDataFrame
covariatesWD <- data.frame(cbind(newcoord2@coords, data$DistWater))
coordinates(covariatesWD) <- newcoord2@coords
gridded(covariatesWD) <- TRUE
covariatesWD <- brick(covariatesWD[3])
names(covariatesWD) <- paste(c("DistWater"))
plot(covariatesWD, axes = FALSE)
#NDVI
covariatesNDVI <- data.frame(cbind(newcoord2@coords, data$NDVI01me_1,data$NDVI04me_1,data$NDVI06me_1,data$NDVI08me_1,data$NDVI10me_1,data$NDVI12me_1,data$NDVI15me_1,data$NDVI16me_1))
coordinates(covariatesNDVI) <- newcoord2@coords
gridded(covariatesNDVI) <- TRUE
covariatesNDVI <- brick(covariatesNDVI[3:10])
names(covariatesNDVI) <- paste("NDVI", c(2001,2004,2006,2008,2010,2012,2015,2016))
plot(covariatesNDVI, axes = FALSE)
#Livestock
covariatesDA <- data.frame(cbind(newcoord2@coords, AbunDom$`2001`, AbunDom$`2004`,AbunDom$`2006`, AbunDom$`2008`, AbunDom$`2010`, AbunDom$`2012`, AbunDom$`2015`, AbunDom$`2016`))
coordinates(covariatesDA) <- newcoord2@coords
gridded(covariatesDA) <- TRUE
covariatesDA <- brick(covariatesDA[3:10])
names(covariatesDA) <- paste("Livestock abundance", c(2001,2004,2006,2008,2010,2012,2015,2016))
plot(covariatesDA, axes = FALSE)
#Kruskal-Wallis test
forKW<-rbind(data.frame(LandUse=data$LU100, NDVI=NDVI01, LA= DA01, Year=2001),
data.frame(LandUse=data$LU100, NDVI=NDVI04, LA= DA04, Year=2004),
data.frame(LandUse=data$LU100, NDVI=NDVI06, LA= DA06, Year=2006),
data.frame(LandUse=data$LU100, NDVI=NDVI08, LA= DA08, Year=2008),
data.frame(LandUse=data$LU100, NDVI=NDVI10, LA= DA10, Year=2010),
data.frame(LandUse=data$LU100, NDVI=NDVI12, LA= DA12, Year=2012),
data.frame(LandUse=data$LU100, NDVI=NDVI15, LA= DA15, Year=2015),
data.frame(LandUse=data$LU100, NDVI=NDVI16, LA= DA16, Year=2016))
forKW<-subset(forKW, !LandUse=="MixedFarmi")
forKW<-subset(forKW, !LandUse=="NA")
forKW$LandUse<-factor(forKW$LandUse)
forKW$Year<-factor(forKW$Year)
#Table D1
Summarize(NDVI ~ LandUse, data = forKW)## LandUse n mean sd min Q1 median Q3
## 1 Pastoralis 688 -0.0040843 1.1639827 -1.692157 -0.8977277 -0.2419050 0.5060687
## 2 Ranching 64 0.0430137 0.8162524 -1.163684 -0.6517830 0.0187579 0.5046282
## 3 RanchingWi 648 -0.2274239 0.8035536 -1.497625 -0.8637238 -0.3042175 0.2235071
## 4 Wildlife 144 0.8434286 0.9011551 -1.238325 0.1442764 0.8613619 1.6088634
## max
## 1 3.642824
## 2 2.133252
## 3 2.475859
## 4 2.811990forKW01<-subset(forKW, Year==2001)
forKW04<-subset(forKW, Year==2004)
forKW06<-subset(forKW, Year==2006)
forKW08<-subset(forKW, Year==2008)
forKW10<-subset(forKW, Year==2010)
forKW12<-subset(forKW, Year==2012)
forKW15<-subset(forKW, Year==2015)
forKW16<-subset(forKW, Year==2016)
#NDVI
kruskal.test(NDVI ~ LandUse, data=forKW01)##
## Kruskal-Wallis rank sum test
##
## data: NDVI by LandUse
## Kruskal-Wallis chi-squared = 18.927, df = 3, p-value = 0.0002831##
## Kruskal-Wallis rank sum test
##
## data: NDVI by LandUse
## Kruskal-Wallis chi-squared = 10.704, df = 3, p-value = 0.01344##
## Kruskal-Wallis rank sum test
##
## data: NDVI by LandUse
## Kruskal-Wallis chi-squared = 27.1, df = 3, p-value = 5.609e-06##
## Kruskal-Wallis rank sum test
##
## data: NDVI by LandUse
## Kruskal-Wallis chi-squared = 31.504, df = 3, p-value = 6.659e-07##
## Kruskal-Wallis rank sum test
##
## data: NDVI by LandUse
## Kruskal-Wallis chi-squared = 20.433, df = 3, p-value = 0.000138##
## Kruskal-Wallis rank sum test
##
## data: NDVI by LandUse
## Kruskal-Wallis chi-squared = 18.934, df = 3, p-value = 0.0002821##
## Kruskal-Wallis rank sum test
##
## data: NDVI by LandUse
## Kruskal-Wallis chi-squared = 37.12, df = 3, p-value = 4.339e-08##
## Kruskal-Wallis rank sum test
##
## data: NDVI by LandUse
## Kruskal-Wallis chi-squared = 26.583, df = 3, p-value = 7.201e-06##
## Kruskal-Wallis rank sum test
##
## data: LA by LandUse
## Kruskal-Wallis chi-squared = 35.303, df = 3, p-value = 1.051e-07##
## Kruskal-Wallis rank sum test
##
## data: LA by LandUse
## Kruskal-Wallis chi-squared = 59.029, df = 3, p-value = 9.475e-13##
## Kruskal-Wallis rank sum test
##
## data: LA by LandUse
## Kruskal-Wallis chi-squared = 64.739, df = 3, p-value = 5.704e-14##
## Kruskal-Wallis rank sum test
##
## data: LA by LandUse
## Kruskal-Wallis chi-squared = 55.005, df = 3, p-value = 6.848e-12##
## Kruskal-Wallis rank sum test
##
## data: LA by LandUse
## Kruskal-Wallis chi-squared = 31.182, df = 3, p-value = 7.781e-07##
## Kruskal-Wallis rank sum test
##
## data: LA by LandUse
## Kruskal-Wallis chi-squared = 45.643, df = 3, p-value = 6.754e-10##
## Kruskal-Wallis rank sum test
##
## data: LA by LandUse
## Kruskal-Wallis chi-squared = 36.79, df = 3, p-value = 5.096e-08##
## Kruskal-Wallis rank sum test
##
## data: LA by LandUse
## Kruskal-Wallis chi-squared = 60.33, df = 3, p-value = 4.997e-13#DW
forKWDW<-data.frame(LandUse=data$LU100, DW = DW)
forKWDW<-subset(forKWDW, !LandUse=="MixedFarmi")
forKWDW<-subset(forKWDW, !LandUse=="NA")
forKWDW$LandUse<-factor(forKWDW$LandUse)
kruskal.test(DW ~ LandUse, data=forKWDW)##
## Kruskal-Wallis rank sum test
##
## data: DW by LandUse
## Kruskal-Wallis chi-squared = 16.661, df = 3, p-value = 0.0008297Therefore, we modeled occupancy probability as a function of a set of covariates that may be acting as the proximate causes of wildlife occurrence across the different land uses: distance to water (DW), NDVI and their quadratic effects and livestock relative abundance (LRA). None of these continuous covariates were highly correlated (Pearson p < 0.5) and all were standardized to a zero mean and unit (1) standard deviation.
The occupancy model for species i at site j and year t was:
\(logit(ψ_i,_j,_t) = u_i,_t + α1_1 * DW_j + α2_i*(DW_j)^2 + α3_i*NDVI_j,_t + α4_i*(NDVI_j,_t)^2 + α5_i*LA_j,_t\)
We modeled detectability in a similar way to occupancy, but allowed detectability to be affected by WVB, which was also standardized as above:
\(logit (p_j,_k,_i,_t )=v_i,_t+β_i*WVB_j,_k\)
Species-level parameters are treated as random effects, each governed by community-level hyperparameters. For instance, we assumed that α1i ~ N(μ_α1,σ_α1) followed a normal distribution where μα1 and σα1 are the mean and standard deviation across the herbivore community. Finally, we incorporated the mean group size effect on the hyperparameter governing species-specific detectability by using a linear regression model on the mean and a linear regression model with log link function to constrain the variance as follows (Kery and Royle, 2016):
\(μ_vi=delta0.p+delta1.p*Group size_i\)
\(log (σ_vi )=phi0.p+phi1.p*Group size_i\)
We calculated the proportion of sites occupied (PSO) per species and per year by dividing the estimated occupied sites by the total number of cells. We also calculated site-specific richness for each year by summing the number of estimated species occurrences per cell. We summarized herbivore species richness in relation to land use by extracting species richness estimates for cells that contained only one land use type to avoid issues of ‘mixed’ cells (i.e., cells containing two or more land uses). Finally, we plotted the relationship of estimated cell-specific herbivore species richness against the three covariates and fit smoothing splines for visualizing trends.
We implemented the model using program JAGS (Plummer, 2016), using the jagsUI package in the R programing language (R Development Core Team, 2016). We used independent, un-informative priors for the community-level hyper-parameters. We checked that those parameters provided strong identifiability by calculating the overlap between each prior and its posterior distribution (i.e., tau << 0.35, Gimenez et al., 2009) using the MCMCvis package. Each of the three parallel Markov chain Monte Carlo (MCMC) chains was run for 150,000 iterations, discarding the first 100,000 as burn-in. We thinned the remaining posterior samples at a rate of 1:10. We evaluated model convergence by visually inspecting chain outputs and using the Gelman-Rubin diagnostic (Gelman and Rubin, 1992). We assessed model fit by estimating the discrepancy between the deviance residuals of the observed and simulated data from the fitted model and by calculating the Bayesian p-value, where values larger than 0.95 or smaller than 0.05 indicate poor fit to the data and a value of 0.5 indicates perfect model fit (Broms et al., 2016). We obtained a p-value of 0.56, indicating that the model provided a good fit to the data. We present parameter-effect sizes in terms of the probability of positive or negative relationships, expressed as the percent of posterior draws above or below zero respectively. We used the Moran’s-I test statistic and visual inspection of correlograms to ensure that the residuals of the occupancy models were not strongly autocorrelated, determining that results were robust to potential pseudo-replication given the grid design.
Model
##The model takes 62 hours to run.
cat("
model{
#Define prior distributions for community-level model parameters
mean.mu.u ~ dunif(0, 1)
mean.u <- log(mean.mu.u)- log(1-mean.mu.u)
delta0.p ~ dnorm(0, 0.01)
delta1.p ~ dnorm(0, 0.01)
phi0.p ~ dnorm(0, 0.01)
phi1.p ~ dnorm(0, 0.01)
mua1 ~ dnorm(0, 0.01)
mua2 ~ dnorm(0, 0.01)
mua3 ~ dnorm(0, 0.01)
mua4 ~ dnorm(0, 0.01)
mua5 ~ dnorm(0, 0.01)
mub1 ~ dnorm(0, 0.01)
sd.a1~dunif(0,10)
tau.a1<-pow(sd.a1,-2)
sd.a2~dunif(0,10)
tau.a2<-pow(sd.a2,-2)
sd.a3~dunif(0,10)
tau.a3<-pow(sd.a3,-2)
sd.a4~dunif(0,10)
tau.a4<-pow(sd.a4,-2)
sd.a5~dunif(0,10)
tau.a5<-pow(sd.a5,-2)
sd.b1~dunif(0,10)
tau.b1<-pow(sd.b1,-2)
sd.mu.u ~ dunif(0,10)
tau.mu.u <- pow(sd.mu.u,-2)
for (i in 1:n) {
#Create priors for species i from the community level prior distributions
#and incorporate group size effect on detectability
mu.u[i] ~ dnorm(mean.u, tau.mu.u)
sd.u[i] ~ dunif(0,10)
tau.u[i] <- pow(sd.u[i],-2)
mu.v[i] ~ dnorm(mean.v[i], tau.mu.v[i])
mean.v[i] <- delta0.p + delta1.p * grsize[i]
log(var.p[i]) <- phi0.p + phi1.p * grsize[i]
tau.mu.v[i] <- 1/var.p[i]
sd.v[i] ~ dunif(0,10)
tau.v[i] <- pow(sd.v[i],-2)
a1[i] ~ dnorm(mua1, tau.a1)
a2[i] ~ dnorm(mua2, tau.a2)
a3[i] ~ dnorm(mua3, tau.a3)
a4[i] ~ dnorm(mua4, tau.a4)
a5[i] ~ dnorm(mua5, tau.a5)
b1[i] ~ dnorm(mub1, tau.b1)
# Incorporate time random effect on the intercepts
for (t in 1:Y){
u[i,t] ~ dnorm(mu.u[i],tau.u[i])
v[i,t] ~ dnorm(mu.v[i],tau.v[i])
}
# Ecological model for latent occurrence
for (j in 1:J) {
for (t in 1:Y) {
logit(psi[j,i,t]) <- u[i,t] + a1[i]*DW[j] + a2[i]*(DW[j]*DW[j]) + a3[i]*NDVI[j,t] + a4[i]*(NDVI[j,t]*NDVI[j,t]) + a5[i]*DA[j,t]
Z[j,i,t] ~ dbern(psi[j,i,t])
# Observation model for detection
for (k in 1:I) {
logit(p[j,k,i,t]) <- v[i,t] + b1[i]*pcov[j,k]
mu.p[j,k,i,t] <- p[j,k,i,t]*Z[j,i,t]
X[j,k,i,t] ~ dbern(mu.p[j,k,i,t])
# Calculate observed deviance
dev[j,k,i,t]<-X[j,k,i,t]*log(mu.p[j,k,i,t])+(1-X[j,k,i,t])*log(1-mu.p[j,k,i,t])
# Predict new observations and compute deviance
X.new[j,k,i,t] ~ dbern(mu.p[j,k,i,t])
dev.sim[j,k,i,t]<-X.new[j,k,i,t]*log(mu.p[j,k,i,t])+(1-X.new[j,k,i,t])*log(1-mu.p[j,k,i,t])
}#k detection
# Calculate residuals
yres[j,i,t]<-psi[j,i,t]-Z[j,i,t]
}#Y time
}#J site
}#i species
# Estimate Bayesian p-value
sum.dev<-sum(dev[,,,])
sum.dev.sim<-sum(dev.sim[,,,])
test<-step(sum.dev.sim - sum.dev)
# Create a loop to determine point level richness estimates for the whole
#community per year.
for(t in 1:Y){
for(j in 1:J){
Nsite[j,t]<- sum(Z[j,1:n,t])
}
}
# Create loop to calculate proportion of sites occupied per species and per
#year
for(i in 1:n){
for(t in 1:Y){
PSO[i,t] <- sum(Z[,i,t])/J
}
}
",file="covarmodel_time_splevel.txt")
#Load all the data
sp.data.time <- list(
n=n,
J=J,
I=I,
X=DetHist,
Y=Y,
pcov=COV,
NDVI=NDVI,
DA=DA,
DW=DW,
grsize=GroupSize$MeanGrSize
)
#Specify the parameters to be monitored
sp.params.time <- c("mean.u",'mean.v',"mu.u","mu.v","mua1","mua2","mua3","mua4","mua5","mub1","delta0.p","delta1.p","phi0.p","phi1.p",'u',"v",'a1','a2','a3',"a4","a5","b1","Nsite","yres","PAO","sum.dev","sum.dev.sim","test")
#Specify the initial values
## Initialize z to be sites where at least 1 detection
zst<-array(1,dim=c(J,n,Y))
sp.inits.time <- function() {
list(Z = zst,
delta0.p = rnorm(1),
delta1.p = rnorm(1),
phi0.p = rnorm(1),
phi1.p = rnorm(1),
mu.u = rnorm(n),
mu.v = rnorm(n),
sd.u = runif(n,0.1,10),
sd.v = runif(n,0.1,10),
a1 = rnorm(n),
a2 = rnorm(n),
a3 = rnorm(n),
a4 = rnorm(n),
a5 = rnorm(n),
b1 = rnorm(n)
)}
#Run the model and call the results ?fit?
model1 <- jags(sp.data.time, sp.inits.time, sp.params.time, "covarmodel_time_splevel.txt", codaOnly=c("Nsite","yres","PAO","sum.dev","sum.dev.sim","test"), parallel = F,
n.chains=3, n.iter=150000, n.burnin=100000, n.thin=10)#### Check for prior influence on posterior estimates
#simulate data from the prior used in your model
#number of iterations should equal the number of draws times the number of chains (although the function will adjust if the correct number of iterations is not specified)
#in JAGS: parameter ~ dnorm(0, 0.01)
nsamp<-model1$mcmc.info$n.samples
PR <- rnorm(nsamp, 0, 10)
MCMCtrace(model1, params = 'mua1', priors = PR, pdf = FALSE)
MCMCtrace(model1, params = 'mua2', priors = PR, pdf = FALSE)
MCMCtrace(model1, params = 'mua3', priors = PR, pdf = FALSE)
MCMCtrace(model1, params = 'mua4', priors = PR, pdf = FALSE)
MCMCtrace(model1, params = 'mua5', priors = PR, pdf = FALSE)
MCMCtrace(model1, params = 'mub1', priors = PR, pdf = FALSE)
MCMCtrace(model1, params = 'delta0.p', priors = PR, pdf = FALSE)
MCMCtrace(model1, params = 'delta1.p', priors = PR, pdf = FALSE)
MCMCtrace(model1, params = 'phi0.p', priors = PR, pdf = FALSE)
MCMCtrace(model1, params = 'phi1.p', priors = PR, pdf = FALSE)
# Following Gimenez et al. (2009), overlap should be < 0.35We assessed potential autocorrelation in our multi-species occupancy model using Moran’s I statistics and visual inspection of correlograms. Both methods were evaluated based on 999 simulations. We observed that 32.5% of the individual models presented highly significant autocorrelation based on the Moran’s I test (p-value < 0.01). However, visual inspection of correlograms generally show that autocorrelation was not a problem in the models (Fig. F1). We opted to not incorporate a spatial auto-covariate structure into the multispecies community-occupancy model. All impala models showed signs of high autocorrelation. Therefore, model results for these species should be interpreted with caution as parameter variation for this species is likely higher than estimated.
Table F1. P-values for the Moran I test for the residuals of each species and year model output.
load("resultsmodel.Rdata")
neigh <- dnearneigh(newcoord2, d1 = 0, d2 = 5000)
lw<-nb2listw(neigh,style="B")
MoranI<-matrix(NA, 15,8)
for (t in 1:8){
for(i in 1:15) {
MI<-moran.mc(as.vector(model1$mean$yres[,i,t]),lw,nsim=999,zero.policy = NULL)
p<-MI$p.value
MoranI[i,t]<-p
}}
rownames(MoranI)<- SpList
colnames(MoranI)<- c(2001,2004,2006,2008,2010,2012,2015,2016)
MoranI<-as.data.frame(MoranI)
MoranI## 2001 2004 2006 2008 2010 2012 2015 2016
## EL 0.006 0.030 0.381 0.089 0.012 0.001 0.006 0.015
## GF 0.001 0.121 0.118 0.057 0.013 0.001 0.005 0.001
## ZB 0.001 0.486 0.001 0.016 0.011 0.016 0.130 0.366
## ZG 0.181 0.719 0.016 0.001 0.027 0.057 0.002 0.001
## TG 0.001 0.065 0.001 0.001 0.037 0.001 0.003 0.068
## GG 0.222 0.015 0.006 0.009 0.114 0.665 0.001 0.011
## KG 0.060 0.733 0.111 0.413 0.173 0.141 0.032 0.921
## IM 0.001 0.003 0.001 0.001 0.001 0.001 0.001 0.013
## BF 0.181 0.210 0.463 0.038 0.975 0.011 0.001 0.021
## ED 0.151 0.518 0.020 0.615 0.128 0.465 0.042 0.073
## OS 0.057 0.560 0.033 0.076 0.110 0.380 0.165 0.263
## WH 0.001 0.184 0.295 0.012 0.024 0.063 0.141 0.256
## OX 0.001 0.003 0.187 0.101 0.865 0.042 0.145 0.001
## WB 0.460 0.775 0.063 0.240 0.001 0.895 0.050 0.992
## GN 0.919 0.583 0.114 0.001 0.985 0.999 0.001 0.001Figure F1. Correlograms showing the spatial autocorrelation of residuals as a function of distance (m) for each species and for each year.
par(mfrow=c(3,5), mar=c(4,4,2,2))
for(i in 1:15) {
bayres<-correlogram(newcoord2, v = as.vector(model1$mean$yres[,i,1]),dist =5000, ns = 999)
}
for(i in 1:15) {
bayres<-correlogram(newcoord2, v = as.vector(model1$mean$yres[,i,2]),dist =5000, ns = 999)
}
for(i in 1:15) {
bayres<-correlogram(newcoord2, v = as.vector(model1$mean$yres[,i,3]),dist =5000, ns = 999)
}
for(i in 1:15) {
bayres<-correlogram(newcoord2, v = as.vector(model1$mean$yres[,i,4]),dist =5000, ns = 999)
}
for(i in 1:15) {
bayres<-correlogram(newcoord2, v = as.vector(model1$mean$yres[,i,5]),dist =5000, ns = 999)
}
for(i in 1:15) {
bayres<-correlogram(newcoord2, v = as.vector(model1$mean$yres[,i,6]),dist =5000, ns = 999)
}
for(i in 1:15) {
bayres<-correlogram(newcoord2, v = as.vector(model1$mean$yres[,i,7]),dist =5000, ns = 999)
}
for(i in 1:15) {
bayres<-correlogram(newcoord2, v = as.vector(model1$mean$yres[,i,8]),dist =5000, ns = 999)
}
## [1] 0.5574Results
Spatiotemporal trends of herbivore species richness
Overall, mean herbivore species richness remained constant over time but showed high intra-annual spatial variability across Laikipia County (Fig. 7). Mean herbivore species richness across years was highest in ‘wildlife only’ areas (mean = 9.28; range = 6.11 to 11.11), followed by ‘ranching and wildlife’ (7.66; 1.73 to 12.00), ‘ranching’ (7.09; 3.73 to 11.17) and ‘pastoralist areas’ (5.92; 1.53 to 10.05). While richness was consistent for wildlife only properties, it was more variable in other land uses, particularly in pastoralist and ranching and wildlife areas (Fig. 7, Fig. 8). Estimated herbivore species richness slightly decreased with increasing distance from water (Fig. 8a), increased toward intermediate NDVI values and decreased again at higher NDVI values (Fig. 8b), and decreased sharply with increasing livestock relative abundance (Fig. 8c). See Fig. 9,10 and 11 for specific year trends.
Fig. 7. Violin plots summarize temporal trends of herbivore species richness (mean, standard deviation and density curve of data) for four different land uses and the total trend.
#Estimated species richness
Nsite <- model1$sims.list$Nsite
site.richness.matrix01 <- cbind(apply(Nsite[,,1],2,mean))
site.richness.matrix04 <- cbind(apply(Nsite[,,2],2,mean))
site.richness.matrix06 <- cbind(apply(Nsite[,,3],2,mean))
site.richness.matrix08 <- cbind(apply(Nsite[,,4],2,mean))
site.richness.matrix10 <- cbind(apply(Nsite[,,5],2,mean))
site.richness.matrix12 <- cbind(apply(Nsite[,,6],2,mean))
site.richness.matrix15 <- cbind(apply(Nsite[,,7],2,mean))
site.richness.matrix16 <- cbind(apply(Nsite[,,8],2,mean))
#Violin plot
violplotEstSR<-rbind(data.frame(SR=site.richness.matrix01,LandUse=data$LU100, Year=c("01")), data.frame(SR=site.richness.matrix04,LandUse=data$LU100, Year=c("04")), data.frame(SR=site.richness.matrix06,LandUse=data$LU100, Year=c("06")),
data.frame(SR=site.richness.matrix08,LandUse=data$LU100, Year=c("08")), data.frame(SR=site.richness.matrix10,LandUse=data$LU100, Year=c("10")), data.frame(SR=site.richness.matrix12,LandUse=data$LU100, Year=c("12")), data.frame(SR=site.richness.matrix15,LandUse=data$LU100, Year=c("15")),
data.frame(SR=site.richness.matrix16,LandUse=data$LU100, Year=c("16")))
violplotEstSR<-subset(violplotEstSR, !LandUse=="MixedFarmi") #Remove mixed farming cells
violplotEstSR<-subset(violplotEstSR, !LandUse=="NA") #Remove NA cells
violplotEstSR$LandUse<-factor(violplotEstSR$LandUse)
violplotEstSR$Year<-factor(violplotEstSR$Year)
levels(violplotEstSR$LandUse) <-c("Pastoralist", "Ranching","Ranching and wildlife","Wildlife only")
violplotEstSRTot<-violplotEstSR
violplotEstSRTot$LandUse<-c("Total")
violplot<-rbind(violplotEstSR,violplotEstSRTot)
viol <- ggplot(violplot, aes(x=Year, y=SR)) + geom_violin() + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), geom="pointrange", color = "black") + facet_wrap (~LandUse, ncol=5)
viol + theme_gray() + labs(x = "Year", y = "Estimated cell\nherbivore sp. richness") + theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15), strip.text = element_text(size=15))
Fig. 8. Maps illustrate spatiotemporal distribution of cell-specific estimates of herbivore species richness across Laikipia, Kenya, for 2001-2016.
SR<-readOGR(dsn = ".", layer = "EstiSpRichness")
intervals <- classIntervals(0:13, n = 6)
year<-c("X2001","X2004","X2006","X2008","X2010","X2012","X2015","X2016")
labelyear<-c(2001,2004,2006,2008,2010,2012,2015,2016)
filename<-c("01.jpg","04.jpg","06.jpg","08.jpg","10.jpg","12.jpg","15.jpg","16.jpg")
for (i in 1:8){
mymap<-tm_shape(prop) + tm_polygons("LU", palette = "Greys", border.col = "white", title="Land use") +
tm_legend(legend.show=T) +
tm_layout(legend.outside = T, scale = 1.3, title = labelyear[i], title.size = 2, frame = F, legend.position = c("left","center")) +
tm_scale_bar(position = c("left","bottom"), size = 0.5) + tm_compass (position = c("left","top")) +
tm_shape(SR) +
tm_bubbles(size = year[i], col= year[i], border.col="white", palette = RdYlGn (6), scale = 1.1, size.max=15, title.size = "Sp richness", title.col = "Sp richness", style = "fixed", breaks = c(0,0.99,3.9,5.9,7.9,9.9,13), legend.col.show = FALSE, legend.size.show = FALSE) +
tm_add_legend("symbol", col = RdYlGn (6), size = c(0.4,0.5,0.6,0.7,0.8,0.9), labels = c("0", "1-4","4-6","6-8","8-10","10-13"), border.col="white")
print(mymap)
save_tmap(tm=mymap, filename = filename[i])
}
library(magick)
list <- list.files(path="./LaikipiaMaps",pattern = "*.jpg", full.names = T)
image1<-image_read(list)
gif<-image_animate(image1, fps = 1)
print(gif)
image_write(gif, "SpRichnessLaikipia.gif")
Fig. 8. Mean estimated cell-specific herbivore species richness in relation to (a) distance to water, (b) NDVI and (c) livestock relative abundance (number of individuals) across rangelands of Laikipia County, Kenya. Dots indicate the land use of each site. Blue lines show the smoothing spline trend. See Appendix G for detail on yearly responses.
#Plot relationship between species richness and WT per management type
#Dataset for plotting
forplotWT<-rbind(data.frame(SR=site.richness.matrix01,WT=data$DistWater,LandUse=data$LU100, Year=2001),
data.frame(SR=site.richness.matrix04,WT=data$DistWater,LandUse=data$LU100, Year=2004),
data.frame(SR=site.richness.matrix06,WT=data$DistWater,LandUse=data$LU100, Year=2006),
data.frame(SR=site.richness.matrix08,WT=data$DistWater,LandUse=data$LU100, Year=2008),
data.frame(SR=site.richness.matrix10,WT=data$DistWater,LandUse=data$LU100, Year=2010),
data.frame(SR=site.richness.matrix12,WT=data$DistWater,LandUse=data$LU100, Year=2012),
data.frame(SR=site.richness.matrix15,WT=data$DistWater,LandUse=data$LU100, Year=2015),
data.frame(SR=site.richness.matrix16,WT=data$DistWater,LandUse=data$LU100, Year=2016))
forplotWT<-subset(forplotWT, !LandUse=="MixedFarmi")
forplotWT<-subset(forplotWT, !LandUse=="NA")
forplotWT$LandUse<-factor(forplotWT$LandUse)
forplotWT$Year<-factor(forplotWT$Year)
levels(forplotWT$LandUse) <-c("Pastoralist", "Ranching","Ranching and wildlife","Wildlife only")
#Total
p1 <- ggplot(forplotWT, aes(x=WT, y=SR)) + geom_point(aes(group=LandUse, color=LandUse)) +
stat_smooth(aes(x=WT, y=SR, group = 1), method = "gam", formula = y ~ s(x), se=F) +
scale_color_manual(name = "Land use", values = c("#FF7F00", "#8B0000", "#66CD00","#006400")) +
labs(x = "Distance to water (m)", y = "Estimated cell herbivore species richness") +
theme(legend.justification=c(1,0), legend.position="none")+
theme(axis.text = element_text(size = 14), axis.title = element_text(size=15)) +
annotate("text", x=12500,y=12, label="a", size=12)
#Plot relationship between species richness and NDVI per management type
#Dataset for plotting
forplotNDVI<-rbind(data.frame(SR=site.richness.matrix01,NDVI=data$NDVI01me_1,LandUse=data$LU100, Year=2001),
data.frame(SR=site.richness.matrix04,NDVI=data$NDVI04me_1,LandUse=data$LU100, Year=2004),
data.frame(SR=site.richness.matrix06,NDVI=data$NDVI06me_1,LandUse=data$LU100, Year=2006),
data.frame(SR=site.richness.matrix08,NDVI=data$NDVI08me_1,LandUse=data$LU100, Year=2008),
data.frame(SR=site.richness.matrix10,NDVI=data$NDVI10me_1,LandUse=data$LU100, Year=2010),
data.frame(SR=site.richness.matrix12,NDVI=data$NDVI12me_1,LandUse=data$LU100, Year=2012),
data.frame(SR=site.richness.matrix15,NDVI=data$NDVI15me_1,LandUse=data$LU100, Year=2015),
data.frame(SR=site.richness.matrix16,NDVI=data$NDVI16me_1,LandUse=data$LU100, Year=2016))
forplotNDVI<-subset(forplotNDVI, !LandUse=="MixedFarmi")
forplotNDVI<-subset(forplotNDVI, !LandUse=="NA")
forplotNDVI$LandUse<-factor(forplotNDVI$LandUse)
forplotNDVI$Year<-factor(forplotNDVI$Year)
levels(forplotNDVI$LandUse) <-c("Pastoralist", "Ranching","Ranching and wildlife","Wildlife only")
#Total
p2 <- ggplot(forplotNDVI, aes(x=NDVI, y=SR)) + geom_point(aes(group=LandUse, color=LandUse)) +
stat_smooth(aes(x=NDVI, y=SR, group = 1), method = "gam", formula = y ~ s(x), se=F) +
scale_color_manual(name = "Land use", values = c("#FF7F00", "#8B0000", "#66CD00","#006400")) +
labs(x = "NDVI", y = "") +
#theme_set(theme_bw()) +
theme(legend.justification=c(1,0), legend.position="none") +
theme(axis.text = element_text(size = 14), axis.title = element_text(size=15))+
annotate("text", x=0.77,y=12, label="b", size=12)
#Plot relationship between species richness and livestock per management type
#Dataset for plotting
forplot<-rbind(data.frame(SR=site.richness.matrix01,DA= AbunDom$`2001`,LandUse=data$LU100, Year=2001),
data.frame(SR=site.richness.matrix04,DA= AbunDom$`2004`,LandUse=data$LU100, Year=2004),
data.frame(SR=site.richness.matrix06,DA= AbunDom$`2006`,LandUse=data$LU100, Year=2006),
data.frame(SR=site.richness.matrix08,DA= AbunDom$`2008`,LandUse=data$LU100, Year=2008),
data.frame(SR=site.richness.matrix10,DA= AbunDom$`2010`,LandUse=data$LU100, Year=2010),
data.frame(SR=site.richness.matrix12,DA= AbunDom$`2012`,LandUse=data$LU100, Year=2012),
data.frame(SR=site.richness.matrix15,DA= AbunDom$`2015`,LandUse=data$LU100, Year=2015),
data.frame(SR=site.richness.matrix16,DA= AbunDom$`2016`,LandUse=data$LU100, Year=2016))
forplot<-subset(forplot, !LandUse=="MixedFarmi")
forplot<-subset(forplot, !LandUse=="NA")
forplot$LandUse<-factor(forplot$LandUse)
forplot$Year<-factor(forplot$Year)
levels(forplot$LandUse) <-c("Pastoralist", "Ranching","Ranching and wildlife","Wildlife only")
#Total
p3 <- ggplot(forplot, aes(x=DA, y=SR)) + geom_point(aes(group=LandUse, col=LandUse)) +
stat_smooth(aes(x=DA, y=SR, group = 1), method = "gam", formula = y ~ s(x), se=F) +
scale_color_manual(name = "Land use", values = c("#FF7F00", "#8B0000", "#66CD00","#006400")) +
labs(x = "Livestock relative abundance", y = "") +
#theme_set(theme_bw()) +
theme(legend.justification=c(1,0), legend.position=c(0.90,0.65),legend.title=element_text(size=14),legend.text=element_text(size=12)) +
theme(axis.text = element_text(size = 14), axis.title = element_text(size=15))+
annotate("text", x=1500,y=12, label="c", size=12)
#Plot all figures and save it
grid.arrange(p1, p2, p3, ncol = 3)
Figure 9. Mean estimated site-specific species richness in relation to distance to water (m) in the rangelands of Laikipia County, Kenya, for the period 2001-2016. Blue lines show the overall trend (Generalized additive model with a spline function mean and standard error for visualization), while other lines show linear trends colored by land use type.
p12 <- ggplot(forplotWT, aes(x=WT, y=SR)) + facet_wrap(~ Year, ncol=3)+ geom_point(aes(group=LandUse, color=LandUse))
p12 + #geom_smooth(method = "gam", formula = y ~ x, se = FALSE) +
stat_smooth(aes(x=WT, y=SR, group = 1), method = "gam", formula = y ~ s(x), se=F) +
scale_color_manual(name = "Land use", values = c("#FF7F00", "#8B0000", "#66CD00","#006400")) +
labs(x = "Distance to water (m)", y = "Estimated cell herbivore species richness") +
theme_set(theme_bw()) +
theme(legend.justification=c(1,0), legend.position=c(1,0.05),legend.title=element_text(size=14),legend.text=element_text(size=12)) +
theme(axis.text = element_text(size = 11), axis.title = element_text(size = 15), strip.text = element_text(size=15))
Figure 10. Mean estimated site-specific species richness in relation to NDVI in the rangelands of Laikipia County, Kenya, for the period 2001-2016. Blue lines show the overall trend (Generalized additive model with a spline function mean and standard error for visualization), while other lines show linear trends colored by land use type.
p22 <- ggplot(forplotNDVI, aes(x=NDVI, y=SR)) + facet_wrap(~ Year, ncol=3)+ geom_point(aes(group=LandUse, color=LandUse))
p22 + #geom_smooth(method = "gam", formula = y ~ x, se = FALSE) +
stat_smooth(aes(x=NDVI, y=SR, group = 1), method = "gam", formula = y ~ s(x), se=F) +
scale_color_manual(name = "Land use", values = c("#FF7F00", "#8B0000", "#66CD00","#006400")) +
labs(x = "NDVI", y = "Estimated cell herbivore species richness") +
theme_set(theme_bw()) +
theme(legend.justification=c(1,0), legend.position=c(1,0.05),legend.title=element_text(size=14),legend.text=element_text(size=12)) +
theme(axis.text = element_text(size = 11), axis.title = element_text(size = 15), strip.text = element_text(size=15))
Figure 11. Mean estimated site-specific species richness in relation to livestock abundance (number of individuals) in the rangelands of Laikipia County, Kenya, for the period 2001-2016. Blue lines show the overall trend (Generalized additive model with a spline function mean and standard error for visualization), while other lines show linear trends colored by land use type.
p23 <- ggplot(forplot, aes(x=DA, y=SR)) + facet_wrap(~ Year, ncol=3)+ geom_point(aes(group=LandUse, color=LandUse))
p23 + #geom_smooth(method = "gam", formula = y ~ x, se = FALSE) +
stat_smooth(aes(x=DA, y=SR, group = 1), method = "gam", formula = y ~ s(x), se=F) +
scale_color_manual(name = "Land use", values = c("#FF7F00", "#8B0000", "#66CD00","#006400")) +
labs(x = "Livestock relative abundance", y = "Estimated cell herbivore species richness") +
theme_set(theme_bw()) +
theme(legend.justification=c(1,0), legend.position=c(1,0.05),legend.title=element_text(size=14),legend.text=element_text(size=12)) +
theme(axis.text = element_text(size = 11), axis.title = element_text(size = 15), strip.text = element_text(size=15))
Occupancy probability summaries
Median site occupancy probability varied widely among species and years, ranging from 0.01 to 0.94. For most species, the averaged proportion of sites occupied (PSO) remained constant across years (Fig. 12). Median PSO for elephants increased from 2001 to 2016. Thomson’s gazelle showed a U-shaped pattern, with mean PSO decreasing from 2001 to 2008 and increasing from 2008 to 2016. Plains zebra occupied the greatest number of sites across years (Median PSO: 0.91), whereas species such as reticulated giraffe, Grevy’s zebra, impala, buffalo, oryx, waterbuck and gerenuk consistently occupied on average < 50% of studied sites (Fig. 12).
Fig. 12. Mean estimated proportion of area occupied (black dots; ± 95% credible intervals (CRI)) for 15 species of herbivores across the Laikipia plateau, Kenya, for the period of 2001-2016. Black triangles indicate the naïve occupancy estimated from the aerial survey dataset prior to the incorporation of detectability.
###############
## Occupancy ##
###############
#Naive occupancy
naiveocc<-matrix(NA, 15,8)
for (t in 1:8){
for (i in 1:dim(DetHist)[3]){
x<-DetHist[,,i,t]
sum1<-apply(x, 1, sum,na.rm=TRUE)
sum2<-replace(sum1,which(sum1>0),1)
sum3<-sum(sum2)
naiveocc[i,t]<-sum3/J
}}
rownames(naiveocc) <- SpList
colnames(naiveocc)<- c(2001,2004,2006,2008,2010,2012,2015,2016)
#Estimated occupancy
Occ <- model1$sims.list$PAO
occ.matrix<-matrix(NA,15,8)
for(t in 1:8){
occ.matrix[,t] <- cbind(apply(Occ[,,t],2,mean))
}
rownames(occ.matrix) <- SpList
colnames(occ.matrix)<- c(2001,2004,2006,2008,2010,2012,2015,2016)
occ.credible_lower<-matrix(NA,15,8)
for(t in 1:8){
occ.credible_lower[,t] <- cbind(apply(Occ[,,t], MARGIN = 2, quantile, prob = 0.025))
}
rownames(occ.credible_lower) <- SpList
colnames(occ.credible_lower)<- c(2001,2004,2006,2008,2010,2012,2015,2016)
occ.credible_upper<-matrix(NA,15,8)
for(t in 1:8){
occ.credible_upper[,t] <- cbind(apply(Occ[,,t], MARGIN = 2, quantile, prob = 0.975))
}
rownames(occ.credible_upper) <- SpList
colnames(occ.credible_upper)<- c(2001,2004,2006,2008,2010,2012,2015,2016)
#Plot
par(mfrow=c(3,1),mar=c(2.3,6,3,2))
plot(0:43,0:43, ylab = "Proportion of sites\noccupied (95% CRI)", ylim=c(0,1),
xlab= "", xlim = c(0.5,42.5), xaxt="n", col="white", bty='l', cex.lab=1.4, cex.axis=1.3)
segments(0:7, occ.credible_lower[1,],0:7,occ.credible_upper[1,], col = "gray55")
points(0:7,occ.matrix[1,], pch = 16)
points(0:7,naiveocc[1,], pch = 17, col = "gray28")
axis(1, 0:7,cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(9:16, occ.credible_lower[2,],9:16,occ.credible_upper[2,], col = "gray55")
points(9:16,occ.matrix[2,], pch = 16)
points(9:16,naiveocc[2,], pch = 17, col = "gray28")
axis(1, 9:16, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(18:25, occ.credible_lower[3,],18:25,occ.credible_upper[3,], col = "gray55")
points(18:25,occ.matrix[3,], pch = 16)
points(18:25,naiveocc[3,], pch = 17, col = "gray28")
axis(1, 18:25, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(27:34, occ.credible_lower[4,],27:34,occ.credible_upper[4,], col = "gray55")
points(27:34,occ.matrix[4,], pch = 16)
points(27:34,naiveocc[4,], pch = 17, col = "gray28")
axis(1, 27:34, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(36:43, occ.credible_lower[5,],36:43,occ.credible_upper[5,], col = "gray55")
points(36:43,occ.matrix[5,], pch = 16)
points(36:43,naiveocc[5,], pch = 17, col = "gray28")
axis(1, 36:43, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
mtext(c("Elephant","Reticulated giraffe","Plains zebra","Grevy's zebra",
"Thomson's gazelle"), side = 3, line = 0, at=c(3.5,12.5,21.5,30.5,39.5), cex=1)
abline(v=8, lty = 2, col = "grey")
abline(v=17, lty = 2, col = "grey")
abline(v=26, lty = 2, col = "grey")
abline(v=35, lty = 2, col = "grey")
plot(0:43,0:43, ylab = "Proportion of sites\noccupied (95% CRI)", ylim=c(0,1),
xlab= "", xlim = c(0.5,42.5), xaxt="n", col="white", bty='l', cex.lab=1.4, cex.axis=1.3)
segments(0:7, occ.credible_lower[6,],0:7,occ.credible_upper[6,], col = "gray55")
points(0:7,occ.matrix[6,], pch = 16)
points(0:7,naiveocc[6,], pch = 17, col = "gray28")
axis(1, 0:7,cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(9:16, occ.credible_lower[7,],9:16,occ.credible_upper[7,], col = "gray55")
points(9:16,occ.matrix[7,], pch = 16)
points(9:16,naiveocc[7,], pch = 17, col = "gray28")
axis(1, 9:16, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(18:25, occ.credible_lower[8,],18:25,occ.credible_upper[8,], col = "gray55")
points(18:25,occ.matrix[8,], pch = 16)
points(18:25,naiveocc[8,], pch = 17, col = "gray28")
axis(1, 18:25, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(27:34, occ.credible_lower[9,],27:34,occ.credible_upper[9,], col = "gray55")
points(27:34,occ.matrix[9,], pch = 16)
points(27:34,naiveocc[9,], pch = 17, col = "gray28")
axis(1, 27:34, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(36:43, occ.credible_lower[10,],36:43,occ.credible_upper[10,], col = "gray55")
points(36:43,occ.matrix[10,], pch = 16)
points(36:43,naiveocc[10,], pch = 17, col = "gray28")
axis(1, 36:43, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
mtext(c("Grant's gazelle","Hartebeest","Impala","African buffalo","Eland"),
side = 3, line = 0, at=c(3.5,12.5,21.5,30.5,39.5), cex=1)
abline(v=8, lty = 2, col = "grey")
abline(v=17, lty = 2, col = "grey")
abline(v=26, lty = 2, col = "grey")
abline(v=35, lty = 2, col = "grey")
plot(0:43,0:43, ylab = "Proportion of sites\noccupied (95% CRI)", ylim=c(0,1),
xlab= "", xlim = c(0.5,42.5), xaxt="n", col="white", bty='l', cex.lab=1.4, cex.axis=1.3)
segments(0:7, occ.credible_lower[11,],0:7,occ.credible_upper[11,], col = "gray55")
points(0:7,occ.matrix[11,], pch = 16)
points(0:7,naiveocc[11,], pch = 17, col = "gray28")
axis(1, 0:7,cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(9:16, occ.credible_lower[12,],9:16,occ.credible_upper[12,], col = "gray55")
points(9:16,occ.matrix[12,], pch = 16)
points(9:16,naiveocc[12,], pch = 17, col = "gray28")
axis(1, 9:16,cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(18:25, occ.credible_lower[13,],18:25,occ.credible_upper[13,], col = "gray55")
points(18:25,occ.matrix[13,], pch = 16)
points(18:25,naiveocc[13,], pch = 17, col = "gray28")
axis(1, 18:25, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(27:34, occ.credible_lower[14,],27:34,occ.credible_upper[14,], col = "gray55")
points(27:34,occ.matrix[14,], pch = 16)
points(27:34,naiveocc[14,], pch = 17, col = "gray28")
axis(1, 27:34, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
segments(36:43, occ.credible_lower[15,],36:43,occ.credible_upper[15,], col = "gray55")
points(36:43,occ.matrix[15,], pch = 16)
points(36:43,naiveocc[15,], pch = 17, col = "gray28")
axis(1, 36:43, cex.axis=0.9, labels=c("01","04","06","08","10","12","15","16"))
mtext(c("Ostrich","Common warthog","Beisa oryx","Defassa waterbuck","Gerenuk"),
side = 3, line = 0, at=c(3.5,12.5,21.5,30.5,39.5), cex=1)
abline(v=8, lty = 2, col = "grey")
abline(v=17, lty = 2, col = "grey")
abline(v=26, lty = 2, col = "grey")
abline(v=35, lty = 2, col = "grey")
We found evidence that the average occupancy probability of the herbivore assemblage was negatively related to distance to water (0.85 probability; Fig. 13). At the species level, we found that five species (giraffe, impala, buffalo, warthog and waterbuck) were attracted to available water (>0.95 probability), whereas occupancy probability for plains zebra and both gazelle species (Thomson’s and Grant’s) increased with increasing distance from water (> 0.95 probability), decreasing slightly again at the largest distances (Fig. 13; Fig. 14). Average occupancy probability for the herbivore assemblage was highest at intermediate levels of NDVI (i.e, lower occupancy at low and high NDVI values). However, there was a large amount of variability in occupancy at low and high NDVI values (Fig. 13; Fig. 14). Nine species (elephant, plains zebra, Thomson’s gazelle, hartebeest, impala, buffalo, eland, warthog and waterbuck) had a strong positive response to NDVI, with occupancy probability increasing with increasing NDVI values (0.99 probability). Gerenuk and Grevy’s zebra occupancy decreased with increasing NDVI (0.99 probability), suggesting that both species concentrate in more arid areas. We found strong evidence (0.99 probability) of a negative effect of livestock relative abundance on the herbivore assemblage average occupancy probability. Nearly every species, with the exception of Thomson’s gazelle, had a negative relationship of occupancy with livestock relative abundance (> 0.97 probability). Occupancy probability for Thomson’s gazelle increased with increasing livestock relative abundance (0.92 probability; Fig. 13; Fig. 14).
Fig 13. Posterior summaries of parameters effects at community and species levels on occupancy and detection probabilities for an herbivore community in Laikipia, Kenya. Distance to water, NDVI and livestock abundance are parameters affecting occupancy probability; whereas, woody vegetation biomass affects detection probability. Gray vertical lines show posterior mean (solid) and 95% credible intervals (CRI; dashed) of the community level hyperparameter. Black dots indicate posterior mean (±95 % CRI; Black bars indicate parameters are different from zero) at species level.
par(mfrow=c(2,3),mar=c(4,10,2,1))
DWeffect<-model1$summary[297:311,c(1:3,7)]
plot(DWeffect[,1],15:1, xlim = c(-1.5,1), xlab = "",
ylab= "", yaxt="n", pch = 16, main = expression("Distance to water"), cex.main=1.4, cex.lab=2, cex.axis=1.1)
abline(v=0, lwd = 1, col = "black")
abline(v=model1$summary[47,c(3,7)], lwd = 1, lty = 2, col = "darkgray")
abline(v=model1$summary[47,1], lwd = 1, col = "darkgray")
segments(DWeffect[,3],15:1, DWeffect[,4],15:1, col="grey", lwd=1)
segments(DWeffect[c(2,3,8,9,14),3],c(14,13,8,7,2), DWeffect[c(2,3,8,9,14),4],c(14,13,8,7,2), col="black", lwd=1)
points(DWeffect[,1],15:1, pch = 16)
axis(2,at=15:1, labels = c("Elephant","Reticulated giraffe","Plains zebra","Grevy's zebra",
"Thomson's gazelle","Grant's gazelle","Hartebeest","Impala",
"African buffalo","Eland","Ostrich","Common warthog","Beisa oryx",
"Defassa waterbuck","Gerenuk"), cex.axis=1.2, las=2)
DWeffect2<-model1$summary[312:326,c(1:3,7)]
plot(DWeffect2[,1],15:1, xlim = c(-1,1), xlab = "",
ylab= "", yaxt="n", pch = 16, main = expression("Distance to water"^"2"), cex.main=1.4, cex.lab=2, cex.axis=1.1)
abline(v=0, lwd = 1, col = "black")
abline(v=model1$summary[48,c(3,7)], lwd = 1, lty = 2, col = "darkgray")
abline(v=model1$summary[48,1], lwd = 1, col = "darkgray")
segments(DWeffect2[,3],15:1, DWeffect2[,4],15:1, col="grey", lwd=1)
points(DWeffect2[,1],15:1, pch = 16)
axis(2,at=15:1, labels = c("Elephant","Reticulated giraffe","Plains zebra","Grevy's zebra",
"Thomson's gazelle","Grant's gazelle","Hartebeest","Impala",
"African buffalo","Eland","Ostrich","Common warthog","Beisa oryx",
"Defassa waterbuck","Gerenuk"), cex.axis=1.2, las=2)
NDVIeffect<-model1$summary[327:341,c(1:3,7)]
plot(NDVIeffect[,1],15:1, xlim = c(-3.5,3.5), xlab = "",
ylab= "", yaxt="n", pch = 16, main = expression("NDVI"), cex.main=1.4, cex.lab=1.7, cex.axis=1.1)
abline(v=0, lwd = 1, col = "black")
abline(v=model1$summary[49,c(3,7)], lwd = 1, lty = 2, col = "darkgray")
abline(v=model1$summary[49,1], lwd = 1, col = "darkgray")
segments(NDVIeffect[,3],15:1, NDVIeffect[,4],15:1, col="grey", lwd=1)
segments(NDVIeffect[c(1,3,4,5,7,8,9,10,12,14,15),3],c(15,13,12,11,9,8,7,6,4,2,1), NDVIeffect[c(1,3,4,5,7,8,9,10,12,14,15),4],c(15,13,12,11,9,8,7,6,4,2,1), col="black", lwd=1)
points(NDVIeffect[,1],15:1, pch = 16)
axis(2,at=15:1, labels = c("Elephant","Reticulated giraffe","Plains zebra","Grevy's zebra",
"Thomson's gazelle","Grant's gazelle","Hartebeest","Impala",
"African buffalo","Eland","Ostrich","Common warthog","Beisa oryx",
"Defassa waterbuck","Gerenuk"), cex.axis=1.2, las=2)
NDVIeffect2<-model1$summary[342:356,c(1:3,7)]
plot(NDVIeffect2[,1],15:1, xlim = c(-1,1), xlab = "Parameter estimate",
ylab= "", yaxt="n", pch = 16, main = expression("NDVI"^"2"), cex.main=1.4, cex.lab=1.7, cex.axis=1.1)
abline(v=0, lwd = 1, col = "black")
abline(v=model1$summary[50,c(3,7)], lwd = 1, lty = 2, col = "darkgray")
abline(v=model1$summary[50,1], lwd = 1, col = "darkgray")
segments(NDVIeffect2[,3],15:1, NDVIeffect2[,4],15:1, col="grey", lwd=1)
segments(NDVIeffect2[c(5,9),3],c(11,7), NDVIeffect2[c(5,9),4],c(11,7), col="black", lwd=1)
points(NDVIeffect2[,1],15:1, pch = 16)
axis(2,at=15:1, labels = c("Elephant","Reticulated giraffe","Plains zebra","Grevy's zebra",
"Thomson's gazelle","Grant's gazelle","Hartebeest","Impala",
"African buffalo","Eland","Ostrich","Common warthog","Beisa oryx",
"Defassa waterbuck","Gerenuk"), cex.axis=1.2, las=2)
livestockeffect<-model1$summary[357:371,c(1:3,7)]
plot(livestockeffect[,1],15:1, xlim = c(-1.5,1), xlab = "Parameter estimate",
ylab= "", yaxt="n", pch = 16, main = expression("Livestock abundance"), cex.main=1.4, cex.lab=1.7, cex.axis=1.1)
abline(v=0, lwd = 1, col = "black")
abline(v=model1$summary[51,c(3,7)], lwd = 1, lty = 2, col = "darkgray")
abline(v=model1$summary[51,1], lwd = 1, col = "darkgray")
segments(livestockeffect[,3],15:1, livestockeffect[,4],15:1, col="grey", lwd=1)
segments(livestockeffect[c(1,2,3,4,5,6,8,9,10,11,12,13,14,15),3],c(15,14,13,12,11,10,8,7,6,5,4,3,2,1), livestockeffect[c(1,2,3,4,5,6,8,9,10,11,12,13,14,15),4],c(15,14,13,12,11,10,8,7,6,5,4,3,2,1), col="black", lwd=1)
points(livestockeffect[,1],15:1,pch = 16)
axis(2,at=15:1, labels = c("Elephant","Reticulated giraffe","Plains zebra","Grevy's zebra",
"Thomson's gazelle","Grant's gazelle","Hartebeest","Impala",
"African buffalo","Eland","Ostrich","Common warthog","Beisa oryx",
"Defassa waterbuck","Gerenuk"), cex.axis=1.2, las=2)
vegcovereffect<-model1$summary[372:386,c(1:3,7)]
plot(vegcovereffect[,1],15:1, xlim = c(-2.5,2), xlab = "Parameter estimate",
ylab= "", yaxt="n", pch = 16, main = expression("Woody vegetation biomass"), cex.main=1.4, cex.lab=1.7, cex.axis=1.1)
abline(v=0, lwd = 1, col = "black")
abline(v=model1$summary[52,c(3,7)], lwd = 1, lty = 2, col = "darkgray")
abline(v=model1$summary[52,1], lwd = 1, col = "darkgray")
segments(vegcovereffect[,3],15:1, vegcovereffect[,4],15:1, col="grey", lwd=1)
segments(vegcovereffect[c(2,3,4,5,6,7,8,10,11,13,14),3],c(14,13,12,11,10,9,8,6,5,3,2), vegcovereffect[c(2,3,4,5,6,7,8,10,11,13,14),4],c(14,13,12,11,10,9,8,6,5,3,2), col="black", lwd=1)
points(vegcovereffect[,1],15:1,pch = 16)
axis(2,at=15:1, labels = c("Elephant","Reticulated giraffe","Plains zebra","Grevy's zebra",
"Thomson's gazelle","Grant's gazelle","Hartebeest","Impala",
"African buffalo","Eland","Ostrich","Common warthog","Beisa oryx",
"Defassa waterbuck","Gerenuk"), cex.axis=1.2, las=2)
Figure 14. Mean occupancy probability for 15 species of herbivores as a function of distance to water (m), NDVI and livestock abundance and detection probabilities as a function of above ground woody vegetation biomass in Laikipia County, Kenya, for the period 2001-2016. Black dashed lines show posterior means and 95% credible intervals for community level hyperparameters.
colforplot<-colors()[c(34,56,657,72,30,153,552,91,142,259,258,109,84,456,208)]
DWater<-seq(1,12000, by= 50)
DWSt <- (DWater - mean.DW) / sd.DW
DWSt2 <- DWSt*DWSt
NDVIseq <- seq(0,1, by=0.01)
NDVISt <- (NDVIseq - mean.NDVI) / sd.NDVI
NDVISt2<-NDVISt*NDVISt
Abund<-seq(0,1500, by=10)
AbundSt <- (Abund - mean.DA) / sd.DA
mean.mu.u = model1$sims.list$mean.u
mean.u = model1$sims.list$mu.u
mu.a1 = model1$sims.list$mua1
mu.a2 = model1$sims.list$mua2
mu.a3 = model1$sims.list$mua3
mu.a4 = model1$sims.list$mua4
mu.a5 = model1$sims.list$mua5
a1 = model1$sims.list$a1
a2 = model1$sims.list$a2
a3 = model1$sims.list$a3
a4 = model1$sims.list$a4
a5 = model1$sims.list$a5
P=seq(1,80, by=1)
Pst <- (P - mean.COV) / sd.COV
mean.mu.v<-model1$sims.list$mean.v
mean.v <- model1$sims.list$mu.v
mean.b1<-model1$sims.list$mub1
b1 <- model1$sims.list$b1
#Dw
psi.mean1 <- plogis(mean(mean.mu.u) + mean(mu.a1)* DWSt + mean(mu.a2)*DWSt2)
psi.mean_dist1 <- matrix(NA, 15000, ncol = length(DWSt))
for (i in 1:15000){
psi.mean_dist1[i,] <- plogis(mean(mean.mu.u[i]) + mean(mu.a1[i])* DWSt + mean(mu.a2[i])*DWSt2)
}
credible_lowerS1 <- apply(psi.mean_dist1, MARGIN = 2, quantile, prob = 0.025)
credible_upperS1 <- apply(psi.mean_dist1, MARGIN = 2, quantile, prob = 0.975)
#NDVI
psi.mean2 <- plogis(mean(mean.mu.u) + mean(mu.a3)* NDVISt + mean(mu.a4)*NDVISt2)
psi.mean_dist2 <- matrix(NA, 15000, ncol = length(NDVISt))
for (i in 1:15000){
psi.mean_dist2[i,] <- plogis(mean(mean.mu.u[i]) + mean(mu.a3[i])* NDVISt + mean(mu.a4[i])*NDVISt2)
}
credible_lowerS2 <- apply(psi.mean_dist2, MARGIN = 2, quantile, prob = 0.025)
credible_upperS2 <- apply(psi.mean_dist2, MARGIN = 2, quantile, prob = 0.975)
#LA
psi.mean3 <- plogis(mean(mean.mu.u) + mean(mu.a5)* AbundSt)
psi.mean_dist3 <- matrix(NA, 15000, ncol = length(AbundSt))
for (i in 1:15000){
psi.mean_dist3[i,] <- plogis(mean(mean.mu.u[i]) + mean(mu.a5[i])* AbundSt)
}
credible_lowerS3 <- apply(psi.mean_dist3, MARGIN = 2, quantile, prob = 0.025)
credible_upperS3 <- apply(psi.mean_dist3, MARGIN = 2, quantile, prob = 0.975)
#detectability
pp.mean <- plogis(mean(mean.mu.v) + mean(mean.b1)*Pst)
pp.mean_dist <- matrix(NA, 15000, ncol = 80)
for (i in 1:15000){
pp.mean_dist[i,] <- plogis(mean(mean.mu.v[i]) + mean(mean.b1[i])*Pst)
}
credible_lowerS4 <- apply(pp.mean_dist, MARGIN = 2, quantile, prob = 0.025)
credible_upperS4 <- apply(pp.mean_dist, MARGIN = 2, quantile, prob = 0.975)
###############################
#PLOT
par(mfrow=c(2,2), mar=c(4,4.2,1,1))
plot(DWater,DWater, type="l", ylim=c(0,1), xlim=c(1,12000), ylab="Occupancy probability",cex.axis=1.2, cex.main=1.5,
xlab="Distance to water (m)", main="", cex.lab = 1.2,
col="white", bty="n")
for (i in 1:n) {
COVpred <- plogis(mean(mean.u[,i]) + mean(a1[,i])*DWSt + mean(a2[,i]*DWSt2))
lines(DWater,COVpred, type="l", ylim=c(0,1), xlim=c(0,1),col = colforplot[i], lwd=2)
}
polygon(c(DWater, rev(DWater)), c(credible_lowerS1, rev(credible_upperS1)), col = adjustcolor("gray", alpha.f = 0.5), border = NA)
lines(DWater, psi.mean1, lwd=2, lty=2)
plot(NDVIseq, NDVIseq, type="l", ylim=c(0,1), xlim=c(0,1), ylab="Occupancy probability",cex.axis=1.2, cex.main=1.5,
xlab="NDVI", main="", cex.lab = 1.2,
col="white", bty="n")
for (i in 1:n) {
COVpred <- plogis(mean(mean.u[,i]) + mean(a3[,i])*NDVISt + + mean(a4[,i])*NDVISt2)
lines(NDVIseq,COVpred, type="l", ylim=c(0,1), xlim=c(0,1),col = colforplot[i], lwd=2)
}
polygon(c(NDVIseq, rev(NDVIseq)), c(credible_lowerS2, rev(credible_upperS2)), col = adjustcolor("gray", alpha.f = 0.5), border = NA)
lines(NDVIseq, psi.mean2, lwd=2, lty=2)
plot(Abund,Abund, type="l", ylim=c(0,1), xlim=c(0,1500), ylab="Occupancy probability",cex.axis=1.2, cex.main=1.5,
xlab="Livestock abundance", main="", cex.lab = 1.2,
col="white", bty="n")
for (i in 1:n) {
COVpred <- plogis(mean(mean.u[,i]) + mean(a5[,i])*AbundSt)
lines(Abund,COVpred, type="l", ylim=c(0,1), xlim=c(0,1500),col = colforplot[i], lwd=2)
}
polygon(c(Abund, rev(Abund)), c(credible_lowerS3, rev(credible_upperS3)), col = adjustcolor("gray", alpha.f = 0.5), border = NA)
lines(Abund, psi.mean3, lwd=2, lty=2)
plot(P,P, type="l", ylim=c(0,0.6), xlim=c(0,80), xlab="Woody vegetation biomass", cex.axis=1.2, cex.main=1.5,
ylab="Detection probability", cex.lab = 1.2,
col="white",bty="n")
for (i in 1:n) {
COVpred <- mean(mean.v[,i]) + mean(b1[,i])*Pst
lines(P,plogis(COVpred), type="l", ylim=c(0,0.7), xlim=c(0,80), col = colforplot[i], lwd=2)
}
polygon(c(P, rev(P)), c(credible_lowerS4, rev(credible_upperS4)), col = adjustcolor("gray", alpha.f = 0.5), border = NA)
lines(P, pp.mean, lwd=2, lty=2)
legend(20,0.6, legend = c("Elephant","Reticulated giraffe","Plains zebra","Grevy's zebra",
"Thomson's gazelle","Grant's gazelle","Hartebeest","Impala",
"African buffalo","Eland","Ostrich","Common warthog","Beisa oryx",
"Defassa waterbuck","Gerenuk"),
lwd = 5, col = colforplot , cex=0.7 , bty = "n", ncol = 2)
Detection probability summaries
We found evidence that species-specific detectability increased with increasing mean group size (0.85 probability, Fig. 15). Median posterior estimates of detection probability for all species and all years were generally low (< 0.6). We found strong evidence (0.99 probability) that detectability at the herbivore assemblage level was negatively affected by woody vegetation biomass (WVB). Evidence of a negative relationship at the species detectability level with WVB was strong (>0.95 probability) for 12 species (Fig. 13; Fig. 14), except for elephant, buffalo and gerenuk. Each of these three species was unresponsive to variation in WVB (i.e., parameter estimates near zero; Fig. 13).
Fig. 15. Estimated mean detection probability and 95% credible intervals (CRI) as a function of group size in Laikipia County, Kenya.
grpsiz<-seq(1,15,by=1)
nsamp<-model1$mcmc.info$n.samples
predgz<-array(NA,dim=c(length(grpsiz),nsamp,2))
for(i in 1:nsamp){
predgz[,i,1]<-plogis(model1$sims.list$delta0.p[i]+model1$sims.list$delta1.p[i]*grpsiz)
predgz[,i,2]<-exp(model1$sims.list$phi0.p[i]+model1$sims.list$phi1.p[i]*grpsiz)
}
# Plot posterior mean
plot(grpsiz, grpsiz, ylab = "Detection probability (95% CRI)", xlab = "Group size", type = "l", lty = 1, xlim = c(0, 16), ylim = c(0, 0.5), frame = F, col = "blue", cex.axis=1.2, cex.main=1.5, cex.lab=1.2)
CLpredgz <- apply(predgz[,,1], MARGIN = 1, quantile, prob = 0.025)
CUpredgz <- apply(predgz[,,1], MARGIN = 1, quantile, prob = 0.975)
polygon(c(grpsiz, rev(grpsiz)), c(CLpredgz, rev(CUpredgz)), col = "gray", border = NA)
lines(grpsiz, apply(predgz[,,1], 1, mean), lwd = 3, col = "blue")
Discussion
Spatiotemporal community dynamics
The spatiotemporal patterns we observed reinforce findings of earlier analyses showing that wildlife friendly properties support a richer wildlife community than other land uses, with pastoral areas supporting, on average, the lowest herbivore species richness of all land uses (Georgiadis et al., 2007; Kinnaird and O’Brien, 2012). This is mainly driven by the negative responses of wild herbivores to livestock abundance. However, averaged trends across the different land uses can mask underlying variability. Our model indicates that spatial variability in herbivore species richness was high, particularly in pastoral areas, due to the high variability in livestock abundance and forage productivity. This high variability suggests that certain pastoral areas can sustain herbivore species richness levels that are comparable to areas dedicated exclusively to wildlife protection (although, wild species may still exist at low densities). Similar patterns have been reported for other regions of southern Kenya (Russell et al., 2018; Tyrrell et al., 2017). These patterns are supported by the notion that planned grazing management in pastoral rangelands can promote vegetation productivity and favor wild species richness compared to unmanaged areas (Odadi et al., 2017).
Species-specific occupancy responses
Echoing observed increases in elephant abundance in Kenya for the last two decades (Ogutu et al., 2016), the total area occupied by elephants expanded over time. This echoes successful conservation measures for elephants across the region (Litoroh et al., 2012). The variation in Thomson’s gazelle occupancy resembles high fluctuations in their numbers observed in southern Kenya and may be related to drought dynamics (Ogutu et al., 2014a). The high occupancy of plains zebra across time coincides with the high abundance reported for Laikipia (Georgiadis et al., 2007; Ogutu et al., 2016), supporting observed trends. We found a high proportion of sites occupied for hartebeest. This result differs from Kinnaird and O’Brien (2012), in which this species was among the ones with the lowest occupancy probability. Our estimates, however, present high uncertainties. More and immediate focused research on the abundance of this sub-species is needed given its limited population size (c. 1000 individuals) and primarily restriction to Laikipia County (Butynski & de Yong, in prep.). A relatively low proportion of the total area was occupied by reticulated giraffe, Grevy’s zebra, oryx, impala, buffalo, waterbuck, and gerenuk. Reticulated giraffe, Grevy’s zebra, and oryx are all endangered species, whereas gerenuk and the local sub-species of waterbuck are near threatened (IUCN, 2019). Populations of these species have been declining in Kenya over the past few decades (Ogutu et al., 2016). Conservation measures are highly important for the protection of each of these endangered species and sub-species given that Laikipia represents one of their last remaining refuges (e.g., 63% of the Grevy’s zebra and 13% of reticulated giraffe global populations occur in Laikipia (Rubenstein et al., 2018)).
All species occupancy probabilities responded negatively to livestock relative abundance except Thomson’s gazelle, which are better adapted to areas with short grasses generated by high livestock grazing pressure (Bhola et al., 2012; Georgiadis et al., 2007; Ogutu et al., 2014a). The increase in occupancy probability of several species at intermediate levels of NDVI is consistent with the preference of African ungulates for areas of intermediate-biomass to maximize their rate of intake of digestible energy (Ogutu et al., 2010), leading to movements of animals toward patches with the highest quality forage available during the dry season (Fynn and Bonyongo, 2010; Tyrrell et al., 2017). Two species, Grevy’s zebra and gerenuk, responded negatively to NDVI, which could be related to their preference for arid areas with low vegetation productivity (Rubenstein et al., 2016), or where competition with livestock is low. Similar to grasslands in southern Kenya, we found that impala, buffalo and waterbuck occurred close to water sources, whereas plains zebras were located farther from water (Ogutu et al., 2014b; Sitters et al., 2009; Tyrrell et al., 2017). We did not find a strong quadratic response to distance to water, unlike what has been observed for the Mara region (Ogutu et al., 2014b, 2010). This could be related to the relatively coarse spatial scale of our study compared with the finer spatial resolution used for the Mara. Besides these general findings at the species level, it is important to note that other factors not included in our models, such as interspecific competition, facilitation among herbivores and predation risk, can affect broad-scale distribution patterns (Bhola et al., 2012; Ogutu and Owen-Smith, 2005).
Implications for management and conservation
Our results highlight that high relative abundance of livestock during the dry season has strong negative effects on wild herbivores occurrence and species richness. For wildlife to persist, the landscape needs to maintain a functional heterogeneity (Fynn et al., 2016), by ensuring continued spatiotemporal heterogeneity in plant communities, which is dictated primarily by grazing pressure, to ensure the diversity and population stability of large herbivores. This will include ensuring a minimum of well-managed ‘wildlife-friendly’ areas with low livestock stocking levels where wild herbivores have access to high-quality foraging opportunities to enable reproduction and population growth (Georgiadis et al., 2007). During the rainy season, more favorable conditions may exist for coexistence of livestock and wild herbivores (Bhola et al., 2012), with livestock potentially facilitating wildlife (Odadi et al., 2011). However, limited forage and water availability during the dry season drives population dynamics in rangelands (Illius and O’Connor, 1999). Therefore, well-managed ‘wildlife-friendly’ areas where wild herbivores have access to reserve forage and that act as refugia during the dry seasons for maintaining population stability, are also important (Fynn and Bonyongo, 2010). Specifically, under scenarios of severe droughts when high livestock stocking rates can amplify negative impacts on wildlife (Ogutu et al., 2014a). A certain amount of land dedicated exclusively to ranching with high livestock stocking densities can be sustained in the matrix of land uses. However, there may exist a threshold to the amount of land dedicated to ranching that may be detrimental at the landscape-level for wildlife. Future research should determine the level and type of livestock and the amount of land dedicated to intense livestock production that does not compromise the wildlife community. Pastoralist areas are highly important for the prosperity of local communities and wildlife across African rangelands (Fynn and Bonyongo, 2010; Ogutu et al., 2014a; Reid et al., 2008; Russell et al., 2018). Pastoralists have evolved to move across the heterogeneous landscape, tracking the spatiotemporal changes in vegetation productivity (Tyrrell et al., 2017). The increasing human population size and high livestock abundances, together with a process of fencing and sedentarization, will inevitably lead to land degradation (Ogutu et al., 2014a; Western et al., 2009a) and overall wildlife decline (du Toit, et al., 2017; Ogutu et al., 2016, 2014a). Livestock numbers need to be controlled, with practices that favors vegetation growth (Odadi et al., 2017). It is also important to provide pastoral communities (and their livestock) the flexibility to move across the landscape to access vegetation reservoirs in wildlife friendly properties during dry periods (Fynn et al., 2016; Ogutu et al., 2014a; Russell et al., 2018). Such seasonal livestock grazing regimes benefit livestock survival, but also maintain a heterogeneous landscape, with greater variation in vegetation structures and increasing nutrient hotspots (glades) that benefit wildlife populations (Fynn et al., 2016). For this flexibility and mobility to persist, deep changes are needed in governance and policy, with well-regulated regional migrations of livestock and with better representation of pastoralists in decision making and broader benefits to those communities that the decisions most directly affect (Ogutu et al., 2014a). The challenge remains on obtaining regional social stability that can prevent deadly conflicts between social groups during times of severe drought (Kessing et al., 2018).
Model insights and limitations
Hierarchical multi-species models provide tools to incorporate rare species that frequently have too few sightings for individual analysis (Zipkin et al., 2009). Even though additional data is needed to improve precision of estimates, long-term trends and inferences on broader herbivore communities and species responses to different variables for species such as Grevy’s zebra, oryx, warthog and gerenuk, species that have been historically difficult to analyze using the aerial survey data, can provide valuable information for conservation.
In accordance with previous research on aerial surveys, our results suggest that species were imperfectly detected (Jachmann, 2002; Schlossberg et al., 2018). Species forming larger group sizes are more likely to be detected than those occurring in smaller group sizes. In addition, vegetation concealment can obstruct visibility and reduce detectability of animals by observers flying at approximately 200 km/h (Jachmann, 2002; Schlossberg et al., 2018). Accounting for the detection process allowed us to estimate species occupancy, minimizing underestimation of site-level richness (Royle and Dorazio, 2008) (see Fig. 16 for a comparison of raw and estimated herbivore species richness).
Our modeling approach on this dataset, however, does still have some limitations. Substituting space for time is a highly debated topic in occupancy models (Guillera-Arroita, 2011; Kery and Royle, 2016; Kendal and White, 2009; Whittington et al., 2015). There is a risk of confusing non-detection with species being absent from segments, biasing occupancy estimates upwards (Guillera-Arroita, 2011; Kendal and White, 2009). Bias is not expected when the animals studied are highly mobile, when there is enough time between replicate surveys and when the number of sites is high (Guillera-Arroita, 2011; Kendal and White, 2009; Whittington et al., 2015). However, in our study, there is the risk of overestimating occupancy for species that may not be mobile enough, and for which the model predicted low detectability by confounding non-detection with absence. Given the low abundance reported in previous studies for hartebeest and ostrich, it is likely that the proportions of sites occupied for these two species was overestimated (Butynski & de Yong, in prep; Ogutu et al., 2016). The occupancy results for these two species should be interpreted with caution. More research is needed to validate our findings, however, incorporating temporal replication into future aerial surveys designs could make similar analyses more robust. Such, designs could be done to maintain the same amount of effort, while providing comparability with historic data (e.g., reducing by half the amount of transects and repeating them twice).
Figure 16. Maps showing raw species richness calculated from the aerial surveys, and estimated posterior mean and standard deviation species richness using a hierarchical community- occupancy model for Laikipia County, Kenya, for the period 2001-2016.
Conclusion
We provide the first approximation of spatiotemporal trends of the herbivore assemblage across the rangelands of Laikipia County, Kenya. Laikipia exemplifies how private and communal land can play a critical role in wildlife conservation across rangelands when livestock abundances are regulated (Keesing et al., 2018; Kinnaird and O’Brien, 2012; Sundaresan and Riginos, 2010). Areas dedicated to wildlife protection are crucial for conservation across these regions and will play a central role in protecting endangered species. These properties, however, are not the only land-use that support wildlife communities. If managed properly, rangelands where pastoralism is practiced can also support rich herbivore communities and may be essential to maintain landscape connectivity. Nevertheless, human population growth and the associated livestock pose large and mounting challenges for conservation across the region as high livestock abundance has severe adverse effects on the wild herbivore community.
Aerial surveys are a common methodology implemented across Africa to study and monitor wildlife (Chase et al., 2016; Jachmann, 2002; Ogutu et al., 2016; Schlossberg et al., 2018). Using multispecies occupancy models that account for detection probability and incorporate rare species into analyses offer new and valuable information for decision makers. Stronger sampling designs incorporating this technique should be considered in wildlife monitoring programs. With most protected areas being too small to sustain their wildlife populations year-round (Western et al., 2009b), repeated and accurate surveys are critical to monitoring changes to wildlife communities and to understand the drivers of biodiversity change across increasingly anthropogenically disturbed landscapes outside formal protected areas.
Acknowledgments
We thank the Directorate of Resource Surveys and Remote Sensing of Kenya (DRSRS) for providing access to the dataset. We thank Grant Connette for helpful discussions on data analysis. We also thank three anonymous reviewers, whose comments helped to improve an earlier version of this paper. Support for this work comes from the Breakefield Critical Conservation Fund, Saudi-Aramco, the German Research Foundation (DFG, Grant No. OG 83/1-2), the European Union’s Horizon 2020 research and innovation programme under Grant Agreement No. 641918 through the AfricanBioServices Project, the Mpala Wildlife Foundation, National Geographic Society and Princeton University. Readers interested in visualizing occupancy probability maps for each species across the 8 sampling years can contact the lead author RDC.
References
Augustine, D.J., McNaughton, S.J., Frank, D.A., 2003. Feedbacks between soil nutrients and large herbivores in a managed savanna ecosystem. Ecol. Appl. 13: 1325-1337.
Bhola, N., Ogutu, J.O., Piepho, H.P., Said, M.Y., Reid, R.S., Hobbs, N.T., Olff, H., 2012. Comparative changes in density and demography of large herbivores in the Masai Mara Reserve and its surrounding human-dominated pastoral ranches in Kenya. Biodivers. Conserv. 21, 1509–1530. doi:10.1007/s10531-012-0261-y
Bouvet, A., Mermoz, S., Le Toan, T., Villard, L., Mathieu, R., Naidoo, L., Asner, G.P., 2018. An above-ground biomass map of African savannahs and woodlands at 25 m resolution derived from ALOS PALSAR. Remote Sens. Environ. 206, 156–173. doi:10.1016/j.rse.2017.12.030
Broms, K.M., Hooten, M.B., Fitzpatrick, R.M., 2016. Model selection and assessment for multi‐species occupancy models. Ecology. 97, 1759–1770
Butynski, T.M., de Jong, Y.A., 2014. Primate conservation in the rangeland agroecosystem of Laikipia county, central Kenya. Primate Conserv. 28, 117–128. doi:10.1896/052.028.0104
Caro, T., Paul, S., 2007. Policy piece: When protection falters. Afr. J. Ecol. 45, 279–281. doi:10.1111/j.1365-2028.2006.00706.x.Lejju
Ceballos, G., Ehrlich, P.R., Soberón, J., Salazar, I., Fay, J.P., 2005. Ecology: Global mammal conservation: What must we manage? Science 309, 603–607. https://doi.org/10.1126/science.1114015
Chase, M.J., Schlossberg, S., Griffin, C.R., Bouché, P.J.C., Djene, S.W., Elkan, P.W., Ferrerira, S., Grossman, F., Kohi, E.M., Landen, K., Omondi, P., Peltier, A., Selier, S.A.J., Sutcliffe, R., 2016. Continent-wide survey reveals massive decline in African savannah elephants. PeerJ 4, e2354. doi:10.7717/peerj.2354
County Government of Laikipia. 2018. Second County integrated development plan 2018-2022. Republic of Kenya. 220 pp. Available at: https://laikipia.go.ke/assets/file/ce3addf6-cidp-2018-2022.pdf
Dorazio, R.M., Royle, J.A., Söderstrom, B., Glimskär, A., 2006. Estimating species richness and accumulation by modeling species occurrence and detectability. Ecology 87, 842–854. doi:10.1002/ardp.200500215
Drescher, M., Brenner, J.C., 2018. The practice and promise of private land conservation. Ecol. Soc. 23, 3. doi:10.5751/ES-10020-230203
du Toit, J.T., Cumming, D.H.M., 1999. Functional significance of ungulate diversity in African savannas and the ecological implications of the spread of pastoralism. Biodivers. Conserv. 8, 1643–1661. https://doi.org/10.1023/A
du Toit, J.T., Cross, P.C., Valeix, M., 2017. Managing the livestock–wildlife interface on rangelands. In: Briske D. (Eds.) Rangeland systems. Springer, Cham, pp. 395-425.
Frank, D.A., McNaughton, S.J., Tracy, B.F., 1998. The ecology of the Earth’s grazing ecosystems. Bioscience 48, 513–521. doi:10.2307/1313313
Fynn, R.W.S., Bonyongo, M.C., 2010. Functional conservation areas and the future of Africa’s wildlife. Afr. J. Ecol. 49, 175–188.
Fynn, R.W.S., Augustine, D.J., Peel, M.J.S., de Garine-Wichatitsky, M., 2016. REVIEW: Strategic management of livestock to improve biodiversity conservation in African savannahs: A conceptual basis for wildlife-livestock coexistence. J. Appl. Ecol. 53, 388–397. doi:10.1111/1365-2664.12591
Gelman, A., Rubin, D., 1992. Inference from iterative simulation using multiple sequences. Stat. Sci. 7, 457–511.
Georgiadis, N.J., Olwero, J.G.N., Ojwang’, G., Romañach, S.S., 2007. Savanna herbivore dynamics in a livestock-dominated landscape: I. Dependence on land use, rainfall, density, and time. Biol. Conserv. 137, 461–472. doi:10.1016/j.biocon.2007.03.005
Gimenez, O., Morgan, B.J.T., Brooks, S.P., 2009. Weak identifiability in models for mark-recapture-recovery data, in: Thomson, D.L., Cooch, E.G., Conroy, M.J. (Eds.), Modeling Demographic Processes in Marked Populations. Springer, New York, USA., pp. 1055–1067.
Goijman, A.P., Conroy, M.J., Bernardos, J.N., Zaccagnini, M.E., 2015. Multi-season regional analysis of multi-species occupancy: Implications for bird conservation in agricultural lands in east-central Argentina. PLoS One 10. doi:10.1371/journal.pone.0130874
Guillera-Arroita, G., 2011. Impact of sampling with replacement in occupancy studies with spatial replication. Methods Ecol. Evol. 2, 401–406. doi:10.1111/j.2041-210X.2011.00089.x
Huntly, N., 1991. Herbivores and the dynamics of communities and ecosystems. Annu. Rev. Ecol. Syst. 22, 477–503. Illius, A.W., O’Connor, T.G., 1999. On the relevance of nonequilibrium concepts to arid and semiarid grazing systems. Ecol. Appl. 9, 799–813.
IUCN, 2019. The IUCN Red List of Threatened Species. Version 2019-1. Jachmann, H., 2002. Comparison of aerial counts with ground counts for large African herbivores. J. Appl. Ecol. 39, 841–852. doi:10.1046/j.1365-2664.2002.00752.x
Keesing, F., Ostfeld, R.S., Okanga, S., Huckett, S., Bayles, B.R., Chaplin-Kramer, R., Fredericks, L.P., Hedlund, T., Kowal, V., Tallis, H., Warui, C.M., Wood, S.A., Allan, B.F., 2018. Integrating livestock and wildlife in an African savanna. Nat. Sustain. 1, 566–573. doi:10.1038/s41893-018-0149-2
Kendall, W.L., White, G.C., 2009. A cautionary note on substituting spatial subunits for repeated temporal sampling in studies of site occupancy. J. Appl. Ecol. 46, 1182–1188 doi: 10.1111/j.1365-2664.2009.01732.x
Kery, M., Royle, J.A., 2016. Applied hierarchical modeling in ecology. Analysis and distribution, abundance and species richness in R and BUGS. Academic Press, London, UK.
Kinnaird, M.F., O’Brien, T.G., 2012. Effects of private-land use, livestock management, and human tolerance on diversity, distribution, and abundance of large African mammals. Conserv. Biol. 26, 1026–1039. doi:10.1111/j.1523-1739.2012.01942.x
Knapp, A.K., Blair, J.M., Briggs, J.M., Collins, S.L., Hartnett, D.C., Johnson, L.C., Towne, E.G., 1999. The keystone role of bison in North American tallgrass prairie. BioScience 49, 39–50.
Letai, J., Lind, J., 2013. Squeezed from all sides: Changing resource tenure and pastoralist innovation on the Laikipia Plateau, Kenya. In: Catley, A., Lind, J., Scoones, I. (Eds.), Pastoralism and development in Africa: Dynamic change at the margins. Routledge and Earthscan, London, pp. 164-176.
Litoroh, M., Omondi, P., Kock, R., Amin, R., 2012. Conservation and management strategy for the elephants in Kenya. 2012-2021. The Kenya Wildlife Service (KWS), Nairobi, Kenya.
Mortensen, B., Danielson, B., Harpole, W.S., Alberti, J., Arnillas, C.A., Biederman, L., Borer, E.T., Cadotte, M.W., Dwyer, J.M., Hagenah, N., Hautier, Y., Peri, P.L., Seabloom, E.W., 2018. Herbivores safeguard plant diversity by reducing variability in dominance. J. Ecol. 106, 101–112. doi:10.1111/1365-2745.12821
Nelson, F., 2008. Are large mammal declines in Africa inevitable? Afr. J. Ecol. 46, 3–4. doi:10.1111/j.1365-2028.2007.00906.x
Newmark, W.D., 2008. Isolation of African protected areas. Front. Ecol. Environ. 6, 321–328. doi:10.1890/070003 Odadi, W.O., Fargione, J., Rubenstein, D.I., 2017. Vegetation, wildlife, and livestock responses to planned grazing management in an African pastoral landscape. L. Degrad. Dev. 28, 2030–2038. doi:10.1002/ldr.2725
Odadi, W.O., Karachi, M.K., Abdulrazak, S.A., Young, T.P., 2011. African wild ungulates compete with or facilitate cattle depending on season. Science 333, 1753–1755. doi:10.1126/science.1208468
Ogutu, J.O., Owen-Smith, N., 2005. Oscillations in large mammal populations: are they related to predation or rainfall? Afr. J. Ecol. 43, 332–339.
Ogutu, J.O., Piepho, H.-P.P., Reid, R.S., Rainy, M.E., Kruska, R.L., Worden, J.S., Nyabenge, M., Hobbs, N.T., Piepho, H.-P.P., Hobbs, N.T., Reid, R.S., Rainy, M.E., Kruska, R.L., Worden, J.S., Nyabenge, M., Hobbs, N.T., 2010. Large herbivore responses to water and settlements in savannas. Ecol. Monogr. 80, 241–266. doi:10.1890/09-0439.1
Ogutu, J.O., Piepho, H.P., Said, M.Y., Kifugo, S.C., 2014a. Herbivore dynamics and range contraction in Kajiado County Kenya: Climate and land use changes, population pressures, governance, policy and human-wildlife conflicts. Open Ecol. J. 7, 9–31. doi:10.2174/1874213001407010009
Ogutu, J.O., Piepho, H.P., Said, M.Y., Ojwang, G.O., Njino, L.W., Kifugo, S.C., Wargute, P.W., 2016. Extreme wildlife declines and concurrent increase in livestock numbers in Kenya: What are the causes? PLoS One 11, 1–46. doi:10.1371/journal.pone.0163249
Ogutu, J.O., Reid, R.S., Piepho, H.P., Hobbs, N.T., Rainy, M.E., Kruska, R.L., Worden, J.S., Nyabenge, M., 2014b. Large herbivore responses to surface water and land use in an East African savanna: Implications for conservation and human-wildlife conflicts. Biodivers. Conserv. 23, 573–596. doi:10.1007/s10531-013-0617-y
Owen-Smith, N., 2004. Functional heterogenity in resourses within landscapes and herbivore population dynamics. Land. Ecol. 19, 761–771.
Owen-Smith, N., 2014. Spatial ecology of large herbivore populations. Ecography 37, 416–430. Pekel, J.F., Cottam, A., Gorelick, N., Belward, A.S., 2016. High-resolution mapping of global surface water and its long-term changes. Nature 540, 418–422. doi:10.1038/nature20584
Pettorelli, N., Ryan, S., Mueller, T., Bunnefeld, N., Jedrzejewska, B., Lima, M., Kausrud, K., 2011. The Normalized Difference Vegetation Index (NDVI): unforeseen successes in animal ecology. Clim. Res. 46:15-27. https://doi.org/10.3354/cr00936
Plummer, M., 2016. JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. (Version 4.2).
Post, E., 2013. Erosion of community diversity and stability by herbivore removal under warming. Proc. R. Soc. B Biol. Sci. 280. doi:10.1098/rspb.2012.2722
Prins, H.H.T., 2000. Competition between wildlife and livestock in Africa, in: Prins, H.H.T., Grootenhuis, J.G., Dolan, T.T. (Eds.), Wildlife Conservation by Sustainable Use. Kluwer Academic Publishers, Boston, MA, USA, pp. 51–80. doi:10.1007/978-94-011-4012-6
QGIS Development Team, 2018. Geographic Information System.
R Development Core Team, 2016. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
Redfern, J. V, Grant, Ri., Biggs, H., Getz, W.M., 2003. Surface-water contrainsts on herbivore foraging in the Kruger National Park, South Africa. Ecology 84, 2092–2107.
Reid, R.S., Galvin, K.A., Kruska, R.S., 2008. Global significance of extensive grazing lands and pastoral societies: An introduction, in: Galvin, K.A., Reid, R.S., Behnke Jr, R.H. (Eds.), Fragmentation in Semi-Arid and Arid Landscapes. Consequences for Human and Natural Systems. Springer, Dordrecht, The Netherlands., pp. 1–24.
Royle, J.A., Dorazio, R.M., 2008. Hierarchical modeling and inference in ecology: the analysis of data from populations, metapopulations and communities. Academic Press, San Diego, California, USA.
Rubenstein, D., Low Mackey, B., Davidson, Z., Kebede, F., King, S.R.B., 2016. Equus grevyi. The IUCN Red List of Threatened Species 2016: e.T7950A89624491. http://dx.doi.org/10.2305/IUCN.UK.2016-3.RLTS.T7950A89624491.en
Rubenstein, D., Parham, J., Stewart, C., Berger-wolf, T., Holmberg, J., Crall, J., Mackey, B.L., Funnel, S., Cockerill, K., Zeke, D., Mate, L., Nzomo, C., Warungu, R., Martins, D., Ontita, V., Omulupi, J., Weston, J., Anyona, G., Chege, G., Kimiti, D., Kaia, T., Gersick, A.S., Rubenstein, N., 2018. The State of Kenya’ s Grevy’ s Zebras and Reticulated Giraffes: Results of the Great Grevy’ s Rally 2018.
Russell, S., Tyrrell, P., Western, D., 2018. Seasonal interactions of pastoralists and wildlife in relation to pasture in an African savanna ecosystem. J. Arid Environ. 154, 70–81. doi:10.1016/j.jaridenv.2018.03.007
Schlossberg, S., Chase, M.J., Griffin, C.R., 2018. Using species traits to predict detectability of animals on aerial surveys: Ecol. Appl. 28, 106–118. doi:10.1002/eap.1632
Schmocker, J., Liniger, H.P., Ngeru, J.N., Brugnara, Y., Auchmann, R., 2015. Trends in mean and extreme precipitation in the Mount Kenya region from observations and reanalyses. International Journal of Climatology 36, 1500–1514. doi:10.1002/joc.4438
Sitters, J., Heitkönig, I.M.A., Holmgren, M., Ojwang’, G.S.O., 2009. Herded cattle and wild grazers partition water but share forage resources during dry years in East African savannas. Biol. Conserv. 142, 738–750. doi:10.1016/j.biocon.2008.12.001
Sundaresan, S.R., Riginos, C., 2010. Lessons learned from biodiversity conservation in the private lands of Laikipia, Kenya. Gt. Plains Res. 20, 17–27.
Tack, J.D., Jakes, A.F., Jones, P.F., Smith, J.T., Newton, R.E., Martin, B.H., Hebblewhite, M., Naugle, D.E., 2019. Beyond protected areas: Private lands and public policy anchor intact pathways for multi-species wildlife migration. Biol. Conserv. 234, 18–27. doi:10.1016/j.biocon.2019.03.017
Tallis, H., Kowal, V.A., Schieltz, J., Chaplin-Kramer, R., Wood, S.A., Musengezi, J., Keesing, F., Okanga, S., Huckett, S., Warui, C.M., Allan, B.F., Ostfeld, R.S., 2017. Can integrating wildlife and livestock enhance ecosystem services in central Kenya? Front. Ecol. Environ. 15, 328–335. doi:10.1002/fee.1501
Tyrrell, P., Russell, S., Western, D., 2017. Seasonal movements of wildlife and livestock in a heterogenous pastoral landscape: Implications for coexistence and community based conservation. Glob. Ecol. Conserv. 12, 59–72. doi:10.1016/j.gecco.2017.08.006
Watson, J.E.M., Dudley, N., Segan, D.B., Hockings, M., 2014. The performance and potential of protected areas. Nature 515, 67–73. doi:10.1038/nature13947
Western, D., Groom, R., Worden, J., 2009a. The impact of subdivision and sedentarization of pastoral lands on wildlife in an African savanna ecosystem. Biol. Conserv. 142, 2538–2546. doi:10.1016/j.biocon.2009.05.025
Western, D., Russell, S., Cuthil, I., 2009b. The status of wildlife in protected areas compared to non-protected areas of Kenya. PLoS One 4, e6140. doi:10.1371/journal.pone.0006140
Whittington, J., Heuer, K., Hunt, B., Hebblewhite, M., Lukacs, P.M., 2015. Estimating occupancy using spatially and temporally replicated snow surveys. Anim. Conserv. 18, 92–101. doi:10.1111/acv.12140
Young, T.P., Porensky, L.M., Riginos, C., Veblen, K.E., Odadi, W.O., Kimuyu, D.M., Charles, G.K., Young, H.S., 2018. Relationships Between Cattle and Biodiversity in Multiuse Landscape Revealed by Kenya Long-Term Exclosure Experiment. Rangel. Ecol. Manag. 71, 281–291. doi:10.1016/j.rama.2018.01.005
Young, T.P., Stubblefield, C.H., Isbell, L.A., 1997. Ants on swollen-thorn acacias: Species coexistence in a simple system. Oecologia 109, 98–107.
Zipkin, E.F., Dewan, A., Andrew Royle, J., 2009. Impacts of forest fragmentation on species richness: A hierarchical approach to community modelling. J. Appl. Ecol. 46, 815–822. doi:10.1111/j.1365-2664.2009.01664.x