Load Data

#################################
### 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 Detections and Spatial Recaptures

## 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

Create Mask to Integrate Density Over

#################################################################
########## 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)

Add Spatial Covariates to Mask

Your goal is to answer the Q, what is more important, habitat or roads, bottom-up or top-down pressue? or Both?

######################################################################################################
### 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")

Fit Models to Data

########################################################################
########################################################################
######## 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

Estimate Abundance Within Kettle-Granby GBPU

###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***

Plot Density

########################################################################
######## 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))

Explore Density-Covariate Relationships

########################################################################
######## 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