#################################
### SECR FOR GRANBY 2015
#################################
#
##Load Packages
library(secr)
library(ggplot2)
library(rgdal)
library(raster)
library(magrittr)
library(RColorBrewer)
###Load and Prepare SECR capture history
bearCH <- read.capthist(
"/Users/clayton.lamb/Dropbox/Work/gm_Granby grizzly/REPORT_Granby/U_Idado_Presentation/CaptHist/Cap_Granby_2015.txt", ##Individual Capture Data
"/Users/clayton.lamb/Dropbox/Work/gm_Granby grizzly/REPORT_Granby/U_Idado_Presentation/CaptHist/Trap_Granby_2015.txt", ##Trap data
detector ="proximity",
covname = 'Sex',
binary.usage = TRUE,
sep="\t")
## Summary of data
summary(bearCH)
## Object class capthist
## Detector type proximity
## Detector number 124
## Average spacing 4233.28 m
## x-range 356509 449370 m
## y-range 5428768 5557584 m
##
## Usage range by occasion
## 1 2 3 4
## min 0 0 1 1
## max 1 1 1 1
##
## Counts by occasion
## 1 2 3 4 Total
## n 36 21 36 42 135
## u 36 9 13 16 74
## f 32 25 15 2 74
## M(t+1) 36 45 58 74 74
## losses 0 0 0 0 0
## detections 44 25 51 57 177
## detectors visited 24 19 33 34 110
## detectors used 108 120 124 124 476
#################################################################
########## PREP MASK
#################################################################
GBmask <- make.mask (traps(bearCH), buffer = 20000, type = 'trapbuffer', spacing = 3000, poly.habitat = T)
##get rid of USA, sorry U of I!!!
GBmask<- GBmask[GBmask$y>5422268,]
##Plot
plot(GBmask, col="black")
plot(traps(bearCH), col="red", add=TRUE)
######################################################################################################
### SPATIAL COVARIATES
######################################################################################################
####################################
######## ADD road density
####################################
#Loading Raster
roadden_raw <- raster("/Users/clayton.lamb/Dropbox/Work/gm_Granby grizzly/REPORT_Granby/U_Idado_Presentation/Spatial/rd_dens_8km",
verbose=FALSE) %>%
projectRaster(crs="+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0",
method="bilinear")
#rescale (standardize, secr had convergence issues with unstandardized version)
roadden <- scale(roadden_raw, center=TRUE, scale=TRUE)
#turn to polygon, secr doens't like raster
roadden_poly <-as(roadden, "SpatialGridDataFrame")
##rename columns
names(roadden_poly) <- "rd_dens"
###Add road dens values to mask
GBmask <- addCovariates(GBmask , roadden_poly)
##PLOT
###Road Density (Standardized)
plot(GBmask, covariate="rd_dens")
####################################
######## ADD habitat
####################################
#Loading Raster
Hab_raw <- raster("/Users/clayton.lamb/Dropbox/Work/gm_Granby grizzly/REPORT_Granby/U_Idado_Presentation/Spatial/habitat_8km", verbose=FALSE) %>%
projectRaster(crs="+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0", method="bilinear")
###Reverse values so a high habitat value is good habitat
x <- values(Hab_raw)
values(Hab_raw) <- ((x-1)*-1)+6
#rescale
Habitat <- scale(Hab_raw, center=TRUE, scale=TRUE)
#turn to polygon, secr doens't like raster
Habitat_poly <- as(Habitat, "SpatialGridDataFrame")
##rename columns
names(Habitat_poly) <- "Habitat"
###Add BEI values to mask
GBmask <- addCovariates(GBmask, Habitat_poly)
##PLOT
##Habitat Quality
plot(GBmask, covariate="Habitat")
########################################################################
########################################################################
######## MODEL
########################################################################
########################################################################
##Model null
bear.fit.NULL <- secr.fit(bearCH, mask=GBmask, model=list(D~1, lambda0~1), detectfn = "HHN", trace=FALSE)
##Model 1
bear.fit1.spat <- secr.fit(bearCH, mask=GBmask, model=list(D~Habitat, lambda0~1), detectfn = "HHN", trace=FALSE)
##Model 2
bear.fit2.spat <- secr.fit(bearCH, mask=GBmask, model=list(D~rd_dens, lambda0~1), detectfn = "HHN", trace=FALSE)
##Model 3
bear.fit3.spat <- secr.fit(bearCH, mask=GBmask, model=list(D~Habitat + rd_dens, lambda0~1), detectfn = "HHN", trace=FALSE)
####################################
######## COMPARE MODELS
####################################
##Which model is more suported by AIC? Pick model with lowest AIC, retain as "Top_Model"
AIC<-AIC(bear.fit.NULL, bear.fit1.spat, bear.fit2.spat, bear.fit3.spat)
#Retain top model
Top_Model_spat<-get(row.names(AIC)[1])
##Return results
AIC
## model detectfn
## bear.fit3.spat D~Habitat + rd_dens lambda0~1 sigma~1 hazard halfnormal
## bear.fit1.spat D~Habitat lambda0~1 sigma~1 hazard halfnormal
## bear.fit2.spat D~rd_dens lambda0~1 sigma~1 hazard halfnormal
## bear.fit.NULL D~1 lambda0~1 sigma~1 hazard halfnormal
## npar logLik AIC AICc dAICc AICcwt
## bear.fit3.spat 5 -631.4346 1272.869 1273.752 0.000 0.9445
## bear.fit1.spat 4 -635.4198 1278.840 1279.419 5.667 0.0555
## bear.fit2.spat 4 -643.6017 1295.203 1295.783 22.031 0.0000
## bear.fit.NULL 3 -649.1635 1304.327 1304.670 30.918 0.0000
###PREP DATA
#
##create new mask, ~ size of grizzly bear population unit
popmask <- make.mask (traps(bearCH), buffer = 12000, type = 'trapbuffer',spacing=1000)
##get rid of USA, sorry U of I!!!
popmask<- popmask[popmask$y>5422268,]
##Add covariates to mask
popmask <- addCovariates(popmask, Habitat_poly)
popmask <- addCovariates(popmask, roadden_poly)
###PLOT
plot(popmask)
plot(traps(bearCH), add=TRUE)
####################################
## ABUNDANCE ESTIMATE
####################################
#
##Predict N, we retain "E.N" as this measure does not assume all individuals
##captured have a home range center inside study area
popest_all <- region.N(Top_Model_spat, popmask, spacing = 3000, se.N = TRUE)
popest_all[1,]
## estimate SE.estimate lcl ucl n
## E.N 95.90467 11.75789 75.4869 121.8451 74
****The KG GBPU population estimate is 95.9 (95% CIs: 75.5 - 121.8 ) grizzly bears***
########################################################################
######## SPATIALLY VIEW DENSITY WITHIN STUDY AREA
########################################################################
##Dens surface for top model
surfaceD <- predictDsurface(Top_Model_spat, mask=popmask, cl.D = TRUE)
##To Raster
surfaceD_rast <- raster(surfaceD, covariate="D.0")
## Warning in matrix(tmp, nrow = ny, ncol = nx, byrow = TRUE): data length
## [17199] is not a sub-multiple or multiple of the number of rows [153]
proj4string(surfaceD_rast) <- "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
##prep color
my.palette <- brewer.pal(n = 9, name = "Reds")
##PLOT
plot(surfaceD_rast *1E5 , axes=FALSE, col=my.palette, box=FALSE, horizontal=TRUE,
legend.args=list(text=as.expression(bquote("Grizzly Density ( / 1000" ~ km^2 ~")" )), side=1, font=2, line=-2.7, cex=1.2))
########################################################################
######## PLOT D ~ Covariates
########################################################################
######################
###### Habitat
######################
Habitat=min(covariates(GBmask)$Habitat):max(covariates(GBmask)$Habitat)
Habitat_seq <- seq(min(covariates(GBmask)$Habitat),max(covariates(GBmask)$Habitat),length.out = 10)
Habitat_replace <- seq(6, 1, length.out = 10)
##Prep input data
newdata3 = data.frame(session=as.factor("Granby_2015"),
t=factor("1", levels=c("1","2","3","4")),
bk=factor("0", levels=c("0", "1")),
PARK = factor("NOT_PARK", levels=c("PARK", "NOT_PARK")),
Habitat= Habitat_seq ,
rd_dens=0,
rd_dens_class= factor("1", levels=c("0", "1")),
rd_dens_class_RC= factor("1", levels=c("0", "1")),
Habitat_replace <- Habitat_replace,
MOTOR = factor("MOTOR", levels=c("MOTOR", "NOT_MOTOR")))
##predict
Habitat_resp_list <- predict(Top_Model_spat,newdata= newdata3 )
##extract data from list into df
Habitat_resp <- as.data.frame(t(data.frame(sapply(Habitat_resp_list,FUN= function(x) data.frame(D=x$estimate[1]*1E5,
SE=x$SE.estimate[1]*1E5)))))
##ensure variables are numeric
Habitat_resp$D <- as.numeric(Habitat_resp$D)
Habitat_resp$SE <- as.numeric(Habitat_resp$SE)
##Add predictor variable and rescale back to non-transformed scale
Habitat_resp$Predictor <- (((Habitat_seq * 2.6) + 5)/2) +0.7
##PLOT
plot1 <- ggplot(Habitat_resp, aes(x=Predictor, y=D))+
geom_ribbon(aes(x=Predictor, ymax= D + SE, ymin=D-SE ), fill = "grey70") +
geom_path()+
xlim(1,6)+
ylim(0,70)+
xlab("Habitat rating (1:low; 6:high)")+
ylab((bquote(paste("Grizzly density (/1000 km"^{2},")"))))+
theme_bw() +theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.title.x = element_text( vjust=-0.5, size=25),
axis.text.x = element_text( size=16),
axis.title.y = element_text( vjust=1.5, size=25),
axis.text.y = element_text( size=16),
strip.text.y = element_text(size=22, angle=270, face="bold"),
axis.line.x = element_line(colour = "black"),
axis.line.y = element_line(colour = "black"),
legend.position=c(0.2,0.35),
legend.text = element_text( size = 18),
legend.title = element_text(size = 25, face = "bold"),
legend.key = element_rect(colour = "black"),
legend.background = element_rect(fill="transparent"))
######################
###### Road Density
######################
rd_dens=min(covariates(GBmask)$rd_dens):max(covariates(GBmask)$rd_dens)
rd_dens_seq <- seq(min(covariates(GBmask)$rd_dens),max(covariates(GBmask)$rd_dens),length.out = 10)
rd_dens_replace <- seq(0, 6, length.out = 10)
##Prep input data
newdata3 = data.frame(session=as.factor("Granby_2015"),
rd_dens= rd_dens_seq,
rd_dens_replace=rd_dens_replace,
Habitat=1.7
)
##predict
rd_dens_resp_list <- predict(Top_Model_spat,newdata= newdata3 )
##extract data from list into df
rd_dens_resp <- as.data.frame(t(data.frame(sapply(rd_dens_resp_list,FUN= function(x) data.frame(D=x$estimate[1]*1E5,
SE=x$SE.estimate[1]*1E5)))))
##ensure variables are numeric
rd_dens_resp$D <- as.numeric(rd_dens_resp$D)
rd_dens_resp$SE <- as.numeric(rd_dens_resp$SE)
##Add predictor variable and rescale back to non-transformed scale
#rd_dens_resp$Predictor <- (((rd_dens_seq * 2.6) + 5)/2) +0.7
rd_dens_resp$Predictor <- rd_dens_replace
##PLOT
plot2 <- ggplot(rd_dens_resp, aes(x=Predictor, y=D))+
geom_ribbon(aes(x=Predictor, ymax= D + SE, ymin=D-SE ), fill = "grey70") +
geom_path()+
xlim(1,6)+
ylim(0,70)+
xlab((bquote(paste("Road density (km/km"^{2},")"))))+
ylab((bquote(paste("Grizzly density (/1000 km"^{2},")"))))+
theme_bw() +theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
axis.title.x = element_text( vjust=-0.5, size=25),
axis.text.x = element_text( size=16),
axis.title.y = element_text( vjust=1.5, size=25),
axis.text.y = element_text( size=16),
strip.text.y = element_text(size=22, angle=270, face="bold"),
axis.line.x = element_line(colour = "black"),
axis.line.y = element_line(colour = "black"),
legend.position=c(0.2,0.35),
legend.text = element_text( size = 18),
legend.title = element_text(size = 25, face = "bold"),
legend.key = element_rect(colour = "black"),
legend.background = element_rect(fill="transparent"))
plot1
plot2