Granby 2015 Population Estimate and Drivers of Density

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

Granby 2015 Population Estimate and Drivers of Density

Clayton T. Lamb, PhD Candidate, University of Alberta

30 October, 2017

Load Data

Summary of Detections and Spatial Recaptures

Create Mask to Integrate Density Over

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?

Fit Models to Data

Estimate Abundance Within Kettle-Granby GBPU

Plot Density

Explore Density-Covariate Relationships