## Chapter 4 Analysis Script ## Josephine B. Smit #Part 1: Modeling elephant site use and detection at the coarse scale, for the entire Ruaha-Rungwa ecosystem #Part 2: Modeling elephant site use and detection at the fine scale, for protected areas in the Ruaha-Rungwa ecosystem ################################################ ### Part 1: Coarse scale, entire ecosystem ##### ################################################ library(unmarked) #version 1.2.5, dated 13th May 2022 library(MuMIn) #version 1.46.0, dated 24th February 2022 library(AICcmodavg) #version 2.3.1, dated 26th August 2020 library(dplyr) #version 1.0.9, dated 28th April 2022 library(magrittr) #version 2.0.3, dated 30th May 2022 library(ggplot2) #version 3.3.6, dated 3rd May 2022 #set working directory # Load detection history detection_history <- read.csv("detection_history.csv", # First variable ("Site") has row.names but not data row.names = "Site") # Load covariate data #detection covariates effort <- read.csv("effort.csv", # First variable ("Site") has row.names but not data row.names = "Site") substrate <- read.csv("Substrate.csv", # First variable ("Site") has row.names but not data row.names = "Site") PA <- read.csv("PA.csv", # First variable ("Site") has row.names but not data row.names = "Site") H <- read.csv("H.csv", # First variable ("Site") has row.names but not data row.names = "Site") C <- read.csv("C.csv", row.names = "Site") R <- read.csv("R.csv", row.names = "Site") W <- read.csv("W.csv", row.names = "Site") #site covariates (different combinations as some covariates are collinear) #CGFSA + H + W + R site_cov <- read.csv("site_cov.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov) #PA + H + W + R site_cov_pa <- read.csv("site_cov_PA.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_pa) #Pop + H + W + R site_cov_pop <- read.csv("site_cov_Pop.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_pop) #DistV + H + W + R site_cov_distv <- read.csv("site_cov_DistV.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_distv) #CGFSA + NDVI + W + R site_cov_ndvi <- read.csv("site_cov_NDVI.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_ndvi) #PA + C + W + R site_cov_pa_c <- read.csv("site_cov_PA_C.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_pa_c) #PA + NDVI + W + R site_cov_pa_ndvi <- read.csv("site_cov_PA_NDVI.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_pa_ndvi) #Pop + C + W + R site_cov_pop_c <- read.csv("site_cov_Pop_C.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_pop_c) #Pop + ndvi + W + R site_cov_pop_ndvi <- read.csv("site_cov_Pop_NDVI.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_pop_ndvi) #DistV + C + W + R site_cov_distv_c <- read.csv("site_cov_DistV_C.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_distv_c) #DistV + NDVI + W + R site_cov_distv_ndvi <- read.csv("site_cov_DistV_NDVI.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_distv_ndvi) #Settlement + H + W + R site_cov_settlement <- read.csv("site_cov_Settlement.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_settlement) #Settlement + C + W + R site_cov_settlement_c <- read.csv("site_cov_Settlement_C.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_settlement_c) #Settlement + NDVI + W + R site_cov_settlement_ndvi <- read.csv("site_cov_Settlement_NDVI.csv", # First variable ("Site") has row.names but not data row.names = "Site") head(site_cov_settlement_ndvi) ####################### #DETECTION MODEL####### ####################### ##Start modeling detection with fixed global model for occupancy #global model for occupancy: Crop + H + R + W #possible detection covariates are: #Effort (effort) #Substrate quality (substrate) #PA - whether or not sampling occasion is in a protected area #Dominant vegetation (H) #Tree cover (C) #Distance to riparian habitat (R) #Water availability (W) #H and C are collinear #PA # Build a new unmarkedFramOccu umf <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H, R = R, W = W, substrate = substrate), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov) #scale observation covariates obsCovs(umf) <- scale(obsCovs(umf)) #scale site covariates siteCovs(umf) <- scale(siteCovs(umf)) summary(umf) global.occu <- occu(~effort + PA + H + R + W + substrate ~CGFSA + W + R + H, umf) summary(global.occu) #run dredge next to find detection model while keeping site covariates fixed occu_dredge_det <- dredge(global.model = global.occu, rank = "AICc", subset=`psi(W)`&`psi(CGFSA)`&`psi(H)`&`psi(R)`) occu_dredge_det[1:10,] #rerun dredge with C instead of H umf <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, C = C, R = R, W = W, substrate = substrate), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov) #scale observation covariates obsCovs(umf) <- scale(obsCovs(umf)) #scale site covariates siteCovs(umf) <- scale(siteCovs(umf)) summary(umf) global.occu <- occu(~effort + PA + C + R + W + substrate ~CGFSA + W + R + H, umf) summary(global.occu) occu_dredge_det <- dredge(global.model = global.occu, rank = "AICc", subset=`psi(W)`&`psi(CGFSA)`&`psi(H)`&`psi(R)`) occu_dredge_det[1:10,] ############################################ #Now model site use #Use/fix best detection model: eff + PA + H ############################################ #site use covariates are: #Crop (CGFSA) #Human population density (Pop) #Proportion of a site that is protected area (PA) #Distance to settlement (DistV) #Building density (Settlement) #Water availability (W) #Distance to riparian habitat (R) #Dominant vegetation type (H) #Tree cover (C) #Vegetation productivity (NDVI) #Crop, PA, Pop,LnDistV, Settlement collinear #Crop & C collinear #H, C, NDVI collinear #Global site use model #Crop + H + R + W #Test with alternative collinear covariates #Crop (CGFSA) # Build a new unmarkedFramOccu umf <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov) # S4 class for occupancy model data summary(umf) #scale observation covariates obsCovs(umf) <- scale(obsCovs(umf)) #scale site covariates siteCovs(umf) <- scale(siteCovs(umf)) summary(umf) #global site use model with CGFSA occu.mglob.det <- occu(formula = ~effort + PA + H # detection formula first ~CGFSA + H + R + W, # occupancy formula second, data = umf) # Summarize summary(occu.mglob.det) # Do Mackenzie-Bailey goodness of fit test for single-season occupancy model mgdet_mb.gof.boot <- mb.gof.test(occu.mglob.det, nsim = 1000) # View Results mgdet_mb.gof.boot #c-hat= 1.1, P=0.3 #dredge CGFSA # Fit all sub-sets of the global site use model while keeping detection covariates fixed occu_dredge <- dredge(global.model = occu.mglob.det, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #Run with Pop instead of Crop umfpop <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_pop) #scale observation covariates obsCovs(umfpop) <- scale(obsCovs(umfpop)) #scale site covariates siteCovs(umfpop) <- scale(siteCovs(umfpop)) summary(umfpop) occu.m1.pop <- occu(formula = ~effort + PA + H # detection formula first ~Pop + H + R + W, # occupancy formula second, data = umfpop) occu_dredge <- dredge(global.model = occu.m1.pop, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #run with PA instead of Crop umfpa <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_pa) #scale observation covariates obsCovs(umfpa) <- scale(obsCovs(umfpa)) #scale site covariates siteCovs(umfpa) <- scale(siteCovs(umfpa)) summary(umfpa) occu.m1.pa <- occu(formula = ~effort + PA + H # detection formula first ~PA + H + R + W, # occupancy formula second, data = umfpa) occu_dredge <- dredge(global.model = occu.m1.pa, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #Run with DistV instead of Crop umfdistv <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_distv) #scale observation covariates obsCovs(umfdistv) <- scale(obsCovs(umfdistv)) #scale site covariates siteCovs(umfdistv) <- scale(siteCovs(umfdistv)) summary(umfdistv) occu.m1.distv <- occu(formula = ~effort + PA + H # detection formula first ~LnDistV + H + R + W, # occupancy formula second, data = umfdistv) occu_dredge <- dredge(global.model = occu.m1.distv, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #Crop + NDVI umfndvi <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_ndvi) #scale observation covariates obsCovs(umfndvi) <- scale(obsCovs(umfndvi)) #scale site covariates siteCovs(umfndvi) <- scale(siteCovs(umfndvi)) summary(umfndvi) occu.m1.ndvi <- occu(formula = ~effort + PA + H # detection formula first ~CGFSA + NDVI + R + W, # occupancy formula second, data = umfndvi) occu_dredge <- dredge(global.model = occu.m1.ndvi, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #Pop + C umfpopc <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov__pop_c) #scale observation covariates obsCovs(umfpopc) <- scale(obsCovs(umfpopc)) #scale site covariates siteCovs(umfpopc) <- scale(siteCovs(umfpopc)) summary(umfpopc) occu.m1.popc <- occu(formula = ~effort + PA + H # detection formula first ~Pop + C + R + W, # occupancy formula second, data = umfpopc) occu_dredge <- dredge(global.model = occu.m1.popc, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #Pop + NDVI umfpopndvi <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_pop_ndvi) #scale observation covariates obsCovs(umfpopndvi) <- scale(obsCovs(umfpopndvi)) #scale site covariates siteCovs(umfpopndvi) <- scale(siteCovs(umfpopndvi)) summary(umfpopndvi) occu.m1.popndvi <- occu(formula = ~effort + PA + H # detection formula first ~Pop + NDVI + R + W, # occupancy formula second, data = umfpopndvi) occu_dredge <- dredge(global.model = occu.m1.popndvi, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #PA + C umfpac <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_pa_c) #scale observation covariates obsCovs(umfpac) <- scale(obsCovs(umfpac)) #scale site covariates siteCovs(umfpac) <- scale(siteCovs(umfpac)) summary(umfpac) occu.m1.pac <- occu(formula = ~effort + PA + H # detection formula first ~PA + C + R + W, # occupancy formula second, data = umfpac) occu_dredge <- dredge(global.model = occu.m1.pac, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #PA + NDVI umfpandvi <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_pa_ndvi) #scale observation covariates obsCovs(umfpandvi) <- scale(obsCovs(umfpandvi)) #scale site covariates siteCovs(umfpandvi) <- scale(siteCovs(umfpandvi)) summary(umfpandvi) occu.m1.pandvi <- occu(formula = ~effort + PA + H # detection formula first ~PA + NDVI + R + W, # occupancy formula second, data = umfpandvi) occu_dredge <- dredge(global.model = occu.m1.pandvi, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #DistV + C umfdistvc <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_distv_c) #scale observation covariates obsCovs(umfdistvc) <- scale(obsCovs(umfdistvc)) #scale site covariates siteCovs(umfdistvc) <- scale(siteCovs(umfdistvc)) summary(umfdistvc) occu.m1.distvc <- occu(formula = ~effort + PA + H # detection formula first ~LnDistV + C + R + W, # occupancy formula second, data = umfdistvc) occu_dredge <- dredge(global.model = occu.m1.distvc, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #DistV + NDVI umfdistvndvi <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_distv_ndvi) #scale observation covariates obsCovs(umfdistvndvi) <- scale(obsCovs(umfdistvndvi)) #scale site covariates siteCovs(umfdistvndvi) <- scale(siteCovs(umfdistvndvi)) summary(umfdistvndvi) occu.m1.distvndvi <- occu(formula = ~effort + PA + H # detection formula first ~LnDistV + NDVI + R + W, # occupancy formula second, data = umfdistvndvi) occu_dredge <- dredge(global.model = occu.m1.distvndvi, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #Settlement umfsettlement <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_settlement) #scale observation covariates obsCovs(umfsettlement) <- scale(obsCovs(umfsettlement)) #scale site covariates siteCovs(umfsettlement) <- scale(siteCovs(umfsettlement)) summary(umfsettlement) occu.m1.settlement <- occu(formula = ~effort + PA + H # detection formula first ~Settlement + H + R + W, # occupancy formula second, data = umfsettlement) occu_dredge <- dredge(global.model = occu.m1.settlement, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #Settlement + C umfsettlementc <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_settlement_c) #scale observation covariates obsCovs(umfsettlementc) <- scale(obsCovs(umfsettlementc)) #scale site covariates siteCovs(umfsettlementc) <- scale(siteCovs(umfsettlementc)) summary(umfsettlementc) occu.m1.settlementc <- occu(formula = ~effort + PA + H # detection formula first ~Settlement + C + R + W, # occupancy formula second, data = umfsettlementc) occu_dredge <- dredge(global.model = occu.m1.settlementc, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] #Settlement + NDVI umfsettlementndvi <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_settlement_ndvi) #scale observation covariates obsCovs(umfsettlementndvi) <- scale(obsCovs(umfsettlementvndvi)) #scale site covariates siteCovs(umfsettlementvndvi) <- scale(siteCovs(umfsettlementndvi)) summary(umfsettlementndvi) occu.m1.settlementndvi <- occu(formula = ~effort + PA + H # detection formula first ~Settlement + NDVI + R + W, # occupancy formula second, data = umfsettlementndvi) occu_dredge <- dredge(global.model = occu.m1.settlementndvi, # the rank argument can use AICc, QAICc, or others (see help page) rank = "AICc", subset=`p(effort)`&`p(PA)`&`p(H)`) occu_dredge[1:10,] # 3. Create model list for those that are relatively well-supported umf <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov) #scale observation covariates obsCovs(umf) <- scale(obsCovs(umf)) #scale site covariates siteCovs(umf) <- scale(siteCovs(umf)) summary(umf) summary(umfndvi) summary(umfpop) # We will fit this model, then put the supported models in a list occu.m1 <- occu(formula = ~effort + PA + H # detection formula first ~CGFSA, # occupancy formula second, data = umfndvi) occu.m2 <- occu(formula = ~effort + PA + H # detection formula first ~Pop, # occupancy formula second, data = umfpop) occu.m3 <- occu(formula = ~effort + PA + H # detection formula first ~CGFSA + NDVI, # occupancy formula second, data = umfndvi) occu.m4 <- occu(formula = ~effort + PA + H # detection formula first ~CGFSA + H, # occupancy formula second, data = umf) occu.m5 <- occu(formula = ~effort + PA + H # detection formula first ~CGFSA + R, # occupancy formula second, data = umfndvi) summary(occu.m1) confint(occu.m1, type="det") confint(occu.m1, type="state") summary(occu.m2) confint(occu.m2, type="det") confint(occu.m2, type="state") summary(occu.m3) confint(occu.m3, type="det") confint(occu.m3, type="state") summary(occu.m4) confint(occu.m4, type="det") confint(occu.m4, type="state") summary(occu.m5) confint(occu.m5, type="det") confint(occu.m5, type="state") #global models occu.glob1 <- occu(formula = ~effort + PA + H # detection formula first ~CGFSA + R + NDVI + W, # occupancy formula second, data = umfndvi) occu.glob2 <- occu(formula = ~effort + PA + H # detection formula first ~Pop + R + H + W, # occupancy formula second, data = umfpop) occu.glob3 <- occu(formula = ~effort + PA + H # detection formula first ~CGFSA + R + H + W, # occupancy formula second, data = umf) # Make model list occu_model_list <- list(occ_m1 = occu.m1, occ_m2 = occu.m2, occ_m3 = occu.m3, occ_m4 = occu.m4, occ_m5 = occu.m5) ## ----modelaveragedpredictions---------- #also want to predict the probability of occupancy/use at each site #while accounting for model uncertainty using model-averaging #use the modavgPred function in the AICcmodavg package to predict values site_cov_occu <- read.csv("site_cov_occu.csv", row.names = "Site") umf <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, PA = PA, H = H), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_occu) #scale observation covariates obsCovs(umf) <- scale(obsCovs(umf)) #scale site covariates siteCovs(umf) <- scale(siteCovs(umf)) summary(umf) occu_modavg_psi_predict <- modavgPred(occu_model_list, # c.hat = 1, # to change variance inflation factor, default = 1) parm.type = "psi", # psi = occupancy, can also be "det" for detection probability newdata = umf@siteCovs)[c("mod.avg.pred", "lower.CL", "upper.CL")] ## Put predictions, CI, and all site covariates into one data frame occu_modavg_psi_predict_df <- data.frame(Predicted = occu_modavg_psi_predict$mod.avg.pred, lower = occu_modavg_psi_predict$lower.CL, upper = occu_modavg_psi_predict$upper.CL, site_cov_occu) # Look at first values head(occu_modavg_psi_predict_df) ## ----modelaveraging_covariates--------------------- #relationships with covariates #use model-averaging to look at how occupancy is predicted to change with covariates #predict the relationship between occupancy and covariates while accounting for model #uncertainty. # First, set-up a new dataframe to predict along a sequence of the covariate. # Predicting requires all covariates, so let's hold the other covariates constant at their mean value occu_CGFSA_newdata <- data.frame(CGFSA = seq(min(site_cov$CGFSA), max(site_cov$CGFSA), by = 0.5), H = mean(site_cov$H), # hold other variables constant R = mean(site_cov$R), W = mean(site_cov$W)) # hold other variables constant # Model-averaged prediction of occupancy and confidence interval occu_CGFSA_pred <- modavgPred(occu_model_list, # c.hat = # to change variance inflation factor, default = 1) parm.type = "psi", # psi = occupancy newdata = occu_CGFSA_newdata)[c("mod.avg.pred", "lower.CL", "upper.CL")] # Put prediction, confidence interval, and covariate values together in a data frame occu_CGFSA_pred_df <- data.frame(Predicted = occu_CGFSA_pred$mod.avg.pred, lower = occu_CGFSA_pred$lower.CL, upper = occu_CGFSA_pred$upper.CL, occu_CGFSA_newdata) # Plot the relationship occu_CGFSA_pred_plot <- ggplot(occu_CGFSA_pred_df, aes(x = CGFSA, y = Predicted)) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5, linetype = "dashed") + geom_path(size = 1) + labs(x = "Percentage cropland (standardized)", y = "Site use probability") + theme_classic() + coord_cartesian(ylim = c(0,1)) + theme(text = element_text(family = "HelveticaNeue", colour = "black"), axis.text = element_text(colour = "black")) occu_CGFSA_pred_plot #plot on original scale #Plot influence of cropland on site use probability in the original scale UN.CPorig<- seq(min(site_cov_ndvi$CGFSA),max(site_cov_ndvi$CGFSA), length.out = 185) ##use datafile with original scale UN.CPpred<- (UN.CPorig-mean(site_cov_ndvi$CGFSA))/sd(site_cov_ndvi$CGFSA) CPnewdata<- data.frame(CGFSA = UN.CPpred, NDVI = 0, R = 0, W = 0) ##use mean of standardized values for the other covariates in your model UN.CP <- predict(occu.m1, type="state", newdata=CPnewdata, appendData = TRUE) plot(UN.CPorig, UN.CP[,"Predicted"], xlab="Proportion of cropland", ylab="Site use probability", type="l", ylim = c(0,1), frame = F, lwd = 2) matlines(UN.CPorig, UN.CP[,3:4], lty = 1, col = "grey", lwd = 1) # plot occupancy as a function of H # First, set-up a new dataframe to predict along a sequence of the covariate. # Predicting requires all covariates, so let's hold the other covariates constant at their mean value occu_H_newdata <- data.frame(H = seq(min(site_cov$H), max(site_cov$H), by = 0.5), CGFSA = mean(site_cov$CGFSA), # hold other variables constant R = mean(site_cov$R), W = mean(site_cov$W)) # hold other variables constant # Model-averaged prediction of occupancy and confidence interval occu_H_pred <- modavgPred(occu_model_list, # c.hat = # to change variance inflation factor, default = 1) parm.type = "psi", # psi = occupancy newdata = occu_H_newdata)[c("mod.avg.pred", "lower.CL", "upper.CL")] # Put prediction, confidence interval, and covariate values together in a data frame occu_H_pred_df <- data.frame(Predicted = occu_H_pred$mod.avg.pred, lower = occu_H_pred$lower.CL, upper = occu_H_pred$upper.CL, occu_H_newdata) # Plot the relationship occu_H_pred_plot <- ggplot(occu_H_pred_df, aes(x = H, y = Predicted)) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5, linetype = "dashed") + geom_path(size = 1) + labs(x = "Percentage miombo (standardized)", y = "Site use probability") + theme_classic() + coord_cartesian(ylim = c(0,1)) + theme(text = element_text(family = "HelveticaNeue", colour = "black"), axis.text = element_text(colour = "black")) occu_H_pred_plot #plot on original scale #Plot influence of habitat on site use probability in the original scale UN.CPorig<- seq(min(site_cov$H),max(site_cov$H), length.out = 185) ##use datafile with original scale UN.CPpred<- (UN.CPorig-mean(site_cov$H))/sd(site_cov$H) CPnewdata<- data.frame(H = UN.CPpred, CGFSA = 0, R = 0, W = 0) ##use mean of standardized values for the other covariates in your model UN.CP <- predict(occu.m4, type="state", newdata=CPnewdata, appendData = TRUE) plot(UN.CPorig, UN.CP[,"Predicted"], xlab="Proportion miombo woodland", ylab="Site use probability", type="l", ylim = c(0,1), frame = F, lwd = 2) matlines(UN.CPorig, UN.CP[,3:4], lty = 1, col = "grey", lwd = 1) # plot occupancy as a function of R # First, set-up a new dataframe to predict along a sequence of the covariate. # Predicting requires all covariates, so let's hold the other covariates constant at their mean value occu_R_newdata <- data.frame(R = seq(min(site_cov$R), max(site_cov$R), by = 0.5), CGFSA = mean(site_cov$CGFSA), # hold other variables constant H = mean(site_cov$H), W = mean(site_cov$W)) # hold other variables constant # Model-averaged prediction of occupancy and confidence interval occu_R_pred <- modavgPred(occu_model_list, # c.hat = # to change variance inflation factor, default = 1) parm.type = "psi", # psi = occupancy newdata = occu_R_newdata)[c("mod.avg.pred", "lower.CL", "upper.CL")] # Put prediction, confidence interval, and covariate values together in a data frame occu_R_pred_df <- data.frame(Predicted = occu_R_pred$mod.avg.pred, lower = occu_R_pred$lower.CL, upper = occu_R_pred$upper.CL, occu_R_newdata) # Plot the relationship occu_R_pred_plot <- ggplot(occu_R_pred_df, aes(x = R, y = Predicted)) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5, linetype = "dashed") + geom_path(size = 1) + labs(x = "Distance to river (standardized)", y = "Site use probability") + theme_classic() + coord_cartesian(ylim = c(0,1)) + theme(text = element_text(family = "HelveticaNeue", colour = "black"), axis.text = element_text(colour = "black")) occu_R_pred_plot #plot on original scale #Plot influence of river on site use probability in the original scale UN.CPorig<- seq(min(site_cov$R),max(site_cov$R), length.out = 185) ##use datafile with original scale UN.CPpred<- (UN.CPorig-mean(site_cov$R))/sd(site_cov$R) CPnewdata<- data.frame(R = UN.CPpred, CGFSA = 0, NDVI = 0, W = 0) ##use mean of standardized values for the other covariates in your model UN.CP <- predict(occu.m5, type="state", newdata=CPnewdata, appendData = TRUE) plot(UN.CPorig, UN.CP[,"Predicted"], xlab="Distance to river (km)", ylab="Site use probability", type="l", ylim = c(0,1), frame = F, lwd = 2) matlines(UN.CPorig, UN.CP[,3:4], lty = 1, col = "grey", lwd = 1) ##Plot influence of NDVI on site use probability in the original scale UN.CPorig<- seq(min(site_cov_ndvi$NDVI),max(site_cov_ndvi$NDVI), length.out = 185) ##use datafile with original scale UN.CPpred<- (UN.CPorig-mean(site_cov_ndvi$NDVI))/sd(site_cov_ndvi$NDVI) CPnewdata<- data.frame(NDVI = UN.CPpred, CGFSA = 0, R = 0, W = 0) ##use mean of standardized values for the other covariates in your model UN.CP <- predict(occu.m3, type="state", newdata=CPnewdata, appendData = TRUE) plot(UN.CPorig, UN.CP[,"Predicted"], xlab="Average NDVI", ylab="Site use probability", type="l", ylim = c(0,1), frame = F, lwd = 2) matlines(UN.CPorig, UN.CP[,3:4], lty = 1, col = "grey", lwd = 1) ##Plot influence of Pop on site use probability in the original scale UN.CPorig<- seq(min(site_cov_pop$Pop),max(site_cov_pop$Pop), length.out = 185) ##use datafile with original scale UN.CPpred<- (UN.CPorig-mean(site_cov_pop$Pop))/sd(site_cov_pop$Pop) CPnewdata<- data.frame(Pop = UN.CPpred, H = 0, R = 0, W = 0) ##use mean of standardized values for the other covariates in your model UN.CP <- predict(occu.m2, type="state", newdata=CPnewdata, appendData = TRUE) plot(UN.CPorig, UN.CP[,"Predicted"], xlab="Human population density (per sq km)", ylab="Site use probability", type="l", ylim = c(0,1), frame = F, lwd = 2) matlines(UN.CPorig, UN.CP[,3:4], lty = 1, col = "grey", lwd = 1) #calculate proportion of area occupied using best ranked model. Need to divide by number of sites #on last row of code. re <- ranef(occu.m1) EBUP <- bup(re, stat="mean") CI <- confint(re, level=0.9) rbind(PAO = c(Estimate = sum(EBUP), colSums(CI)) / 185) #create fit stats object fitstats <- function(occu.m1, method = "nonparboot") { observed <- getY(occu.m1@data) expected <- fitted(occu.m1) resids <- residuals(occu.m1, method = "nonparboot") sse <- sum(resids^2, na.rm = TRUE) chisq <- sum((observed - expected)^2 / expected, na.rm = TRUE) freeTuke <- sum((sqrt(observed) - sqrt(expected))^2, na.rm = TRUE) out <- c(SSE = sse, Chisq = chisq, freemanTukey = freeTuke) return(out) } pb <- parboot(occu.m1, fitstats, nsim = 1000, report = TRUE, method = "nonparboot") pb #All tests are OK # Do Mackenzie-Bailey goodness of fit test for top model mgm1_mb.gof.boot <- mb.gof.test(occu.m1, # Demonstrate with small number of sims (10), # but then change to large number (e.g. 1000) nsim = 1000) # View Results mgm1_mb.gof.boot #chat=1.1, p=0.28 # Do Mackenzie-Bailey goodness of fit test for global model mglob1_mb.gof.boot <- mb.gof.test(occu.glob1, # Demonstrate with small number of sims (10), # but then change to large number (e.g. 1000) nsim = 1000) # View Results mglob1_mb.gof.boot #c-hat= 1.09, p=0.278 #predictions of p for each site p_top <- predict(occu.m1,type="det") #predictions of psi for each site psi_top <- predict(occu.m1,type="state") ################################################ ### Part 2: Coarse scale, protected areas ###### ################################################ # Load detection history detection_history <- read.csv("detection_history.csv", row.names = "Site") # Examine data head(detection_history) # Load covariate data - each observation level covariate is a different csv #sitecovs can be one csv effort <- read.csv("effort.csv", # First variable ("Site") has row.names but not data row.names = "Site") substrate <- read.csv("substrate.csv", row.names = "Site") H <- read.csv("H.csv", row.names = "Site") bound <- read.csv("Bound.csv", row.names = "Site") carcass <- read.csv("Carcass15.csv", row.names = "Site") hum <- read.csv("Hum.csv", row.names = "Site") post <- read.csv("Post.csv", row.names = "Site") R <- read.csv("R.csv", row.names = "Site") W <- read.csv("W.csv", row.names = "Site") C <- read.csv("C.csv", row.names = "Site") site_cov <- read.csv("site_cov.csv", row.names = "Site") #collinear covariates #H, NDVI & C #Hum, Bound, Carcass ######################### #Model Detection first ## ######################### #global model for occupancy/use is Carcass + H + W + R + Post #possible detection covariates are: #Effort (effort) #Substrate quality (substrate) #Dominant vegetation (H) #Tree cover (C) #Distance to riparian habitat (R) #Water availability (W) #Carcass occurrence probability (Carcass) #Probability of illegal human use (Hum) #Distance to PA boundary (Bound) #Distance to ranger post (Post) #fit global model umf <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, substrate = substrate, H = H, W = W, R = R, post = post, bound = bound), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov) # S4 class for occupancy model data summary(umf) #scale observation covariates obsCovs(umf) <- scale(obsCovs(umf)) #scale site covariates siteCovs(umf) <- scale(siteCovs(umf)) summary(umf) #global model for occupancy/use and detection occu.mglob <- occu(formula = ~effort + W + H + R + post + substrate + bound # detection formula first ~Carcass + H + W + R + Post, # occupancy formula second, data = umf) # Summarize summary(occu.mglob) # Do Mackenzie-Bailey goodness of fit test for single-season occupancy model mg_mb.gof.boot <- mb.gof.test(occu.mglob, # Demonstrate with small number of sims (10), # but then change to large number (e.g. 1000) nsim = 50) # View Results mg_mb.gof.boot #Fit all sub-sets of the global model #fix site covars when finding best detection model occu_dredge_det <- dredge(global.model = occu.mglob, rank = "AICc", subset=`psi(Carcass)`&`psi(H)`&`psi(W)`&`psi(R)`&`psi(Post)`) occu_dredge_det[1:10,] #redo dredge for detection model with Hum instead of Bound # Build unmarkedFramOccu umf1 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, substrate = substrate, H = H, W = W, R = R, post = post, hum = hum), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov) # S4 class for occupancy model data summary(umf1) #scale observation covariates obsCovs(umf1) <- scale(obsCovs(umf1)) #scale site covariates siteCovs(umf1) <- scale(siteCovs(umf1)) summary(umf1) #global model for occupancy and detection occu.mglob1 <- occu(formula = ~effort + W + H + R + post + substrate + hum # detection formula first ~Carcass + H + W + R + Post, # occupancy formula second, data = umf1) # Summarize summary(occu.mglob1) #fix site covars when finding best detection model occu_dredge_det_hum <- dredge(global.model = occu.mglob1, rank = "AICc", subset=`psi(Carcass)`&`psi(H)`&`psi(W)`&`psi(R)`&`psi(Post)`) occu_dredge_det_hum[1:10,] #redo dredge for detection model with Carcass instead of Bound # Build unmarkedFramOccu umf2 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, substrate = substrate, H = H, W = W, R = R, post = post, carcass = carcass), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov) # S4 class for occupancy model data summary(umf2) #scale observation covariates obsCovs(umf2) <- scale(obsCovs(umf2)) #scale site covariates siteCovs(umf2) <- scale(siteCovs(umf2)) summary(umf2) #global model for occupancy and detection occu.mglob2 <- occu(formula = ~effort + W + H + R + post + substrate + carcass # detection formula first ~Carcass + H + W + R + Post, # occupancy formula second, data = umf2) # Summarize summary(occu.mglob2) #fix site covars when finding best detection model occu_dredge_det_carcass <- dredge(global.model = occu.mglob2, rank = "AICc", subset=`psi(Carcass)`&`psi(H)`&`psi(W)`&`psi(R)`&`psi(Post)`) occu_dredge_det_carcass[1:10,] #best detection model so far has carcass, effort, H, post, R #redo dredge for this best detection model with C instead of H # Build unmarkedFramOccu umf3 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, substrate = substrate, C = C, W = W, R = R, post = post, carcass = carcass), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov) # S4 class for occupancy model data summary(umf3) #scale observation covariates obsCovs(umf3) <- scale(obsCovs(umf3)) #scale site covariates siteCovs(umf3) <- scale(siteCovs(umf3)) summary(umf3) #global model for occupancy and detection occu.mglob3 <- occu(formula = ~effort + W + C + R + post + substrate + carcass # detection formula first ~Carcass + H + W + R + Post, # occupancy formula second, data = umf3) # Summarize summary(occu.mglob3) #fix site covars when finding best detection model occu_dredge_det_C <- dredge(global.model = occu.mglob3, rank = "AICc", subset=`psi(Carcass)`&`psi(H)`&`psi(W)`&`psi(R)`&`psi(Post)`) occu_dredge_det_C[1:10,] ###################################################### ## Model occupancy/use with best detection model ##### ###################################################### #best detection model p(carcass + effort + H + post + R) #Site use covariates are: #Dominant vegetation (H) #Tree cover (C) #Distance to riparian habitat (R) #Water availability (W) #Carcass occurrence probability (Carcass) #Probability of illegal human use (Hum) #Distance to PA boundary (Bound) #Distance to ranger post (Post) #Vegetation productivity (NDVI) #global model with Crop # Build unmarkedFramOccu umf4 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, R = R, post = post, carcass = carcass), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov) # S4 class for occupancy model data summary(umf4) #scale observation covariates obsCovs(umf4) <- scale(obsCovs(umf4)) #scale site covariates siteCovs(umf4) <- scale(siteCovs(umf4)) summary(umf4) #fix det covars when finding best occupancy model #start with global occupancy model psi(Carcass + H + W + R + post) occu.mglob <- occu(formula = ~effort + H + R + post + carcass # detection formula first ~Carcass + H + W + R + Post, # occupancy formula second, data = umf4) # Summarize summary(occu.mglob) #dredge with fixed detection models occu_dredge_psi <- dredge(global.model = occu.mglob, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psi[1:10,] ## Hum instead of carcass site_cov_hum <- read.csv("site_cov_Hum.csv", row.names = "Site") # Build unmarkedFramOccu umf5 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_hum) #scale observation covariates obsCovs(umf5) <- scale(obsCovs(umf5)) #scale site covariates siteCovs(umf5) <- scale(siteCovs(umf5)) summary(umf5) #global model for occupancy and detection occu.mglob6 <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Hum + H + W + R + Post, # occupancy formula second, data = umf5) # Summarize summary(occu.mglob6) #dredge with fixed detection model occu_dredge_psi5 <- dredge(global.model = occu.mglob6, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psi5[1:10,] ## Bound instead of carcass site_cov_bound <- read.csv("site_cov_Bound.csv", row.names = "Site") # Build unmarkedFramOccu umf6 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_bound) #scale observation covariates obsCovs(umf6) <- scale(obsCovs(umf6)) #scale site covariates siteCovs(umf6) <- scale(siteCovs(umf6)) summary(umf6) #global model for occupancy and detection occu.mglob7 <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Bound + H + W + R + Post, # occupancy formula second, data = umf6) # Summarize summary(occu.mglob7) #dredge with fixed detection model occu_dredge_psi7 <- dredge(global.model = occu.mglob7, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psi7[1:10,] ## C instead of H site_cov_c <- read.csv("site_cov_C.csv", row.names = "Site") # Build unmarkedFramOccu umf7 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_c) #scale observation covariates obsCovs(umf7) <- scale(obsCovs(umf7)) #scale site covariates siteCovs(umf7) <- scale(siteCovs(umf7)) summary(umf7) #global model for occupancy and detection occu.mglob8 <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Carcass + C + W + R + Post, # occupancy formula second, data = umf7) # Summarize summary(occu.mglob8) ##dredge with fixed detection model occu_dredge_psi9 <- dredge(global.model = occu.mglob8, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psi9[1:10,] ## NDVI site_cov_ndvi <- read.csv("site_cov_NDVI.csv", row.names = "Site") # Build unmarkedFramOccu umf8 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_ndvi) #scale observation covariates obsCovs(umf8) <- scale(obsCovs(umf8)) #scale site covariates siteCovs(umf8) <- scale(siteCovs(umf8)) summary(umf8) #global model for occupancy and detection occu.mglob9 <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Carcass + NDVIave + W + R + Post, # occupancy formula second, data = umf8) # Summarize summary(occu.mglob9) ##dredge with fixed detection model occu_dredge_psi11 <- dredge(global.model = occu.mglob9, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psi11[1:10,] #NDVIstd site_cov_ndvistd <- read.csv("site_cov_NDVIstd.csv", row.names = "Site") # Build unmarkedFramOccu umfndvistd <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_ndvistd) #scale observation covariates obsCovs(umfndvistd) <- scale(obsCovs(umfndvistd)) #scale site covariates siteCovs(umfndvistd) <- scale(siteCovs(umfndvistd)) summary(umfndvistd) #global model for occupancy and detection occu.mglobndvistd <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Carcass + NDVIstd + W + R + Post, # occupancy formula second, data = umfndvistd) # Summarize summary(occu.mglobndvistd) ##dredge with fixed detection model occu_dredge_psindvistd <- dredge(global.model = occu.mglobndvistd, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psindvistd[1:10,] ## C and Hum site_cov_c_hum <- read.csv("site_cov_C_Hum.csv", row.names = "Site") # Build unmarkedFramOccu umf9 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_c_hum) #scale observation covariates obsCovs(umf9) <- scale(obsCovs(umf9)) #scale site covariates siteCovs(umf9) <- scale(siteCovs(umf9)) summary(umf9) #global model for occupancy and detection occu.mglob10 <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Hum + C + W + R + Post, # occupancy formula second, data = umf9) # Summarize summary(occu.mglob10) ##dredge with fixed detection model occu_dredge_psi13 <- dredge(global.model = occu.mglob10, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psi13[1:10,] ## C and Bound site_cov_c_bound <- read.csv("site_cov_C_Bound.csv", row.names = "Site") # Build unmarkedFramOccu umf10 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_c_bound) #scale observation covariates obsCovs(umf10) <- scale(obsCovs(umf10)) #scale site covariates siteCovs(umf10) <- scale(siteCovs(umf10)) summary(umf10) #global model for occupancy and detection occu.mglob11 <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Bound + C + W + R + Post, # occupancy formula second, data = umf10) # Summarize summary(occu.mglob11) ##dredge with fixed detection model occu_dredge_psi15 <- dredge(global.model = occu.mglob11, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psi15[1:10,] ## NDVI and Hum site_cov_ndvi_hum <- read.csv("site_cov_NDVI_Hum.csv", row.names = "Site") # Build unmarkedFramOccu umf11 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_ndvi_hum) #scale observation covariates obsCovs(umf11) <- scale(obsCovs(umf11)) #scale site covariates siteCovs(umf11) <- scale(siteCovs(umf11)) summary(umf11) #global model for occupancy and detection occu.mglob12 <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Hum + NDVIave + W + R + Post, # occupancy formula second, data = umf11) # Summarize summary(occu.mglob12) ##dredge with fixed detection model occu_dredge_psi17 <- dredge(global.model = occu.mglob12, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psi17[1:10,] ## NDVIstd and Hum site_cov_ndvistd_hum <- read.csv("site_cov_NDVIstd_Hum.csv", row.names = "Site") # Build unmarkedFramOccu umfndvistdhum <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_ndvistd_hum) #scale observation covariates obsCovs(umfndvistdhum) <- scale(obsCovs(umfndvistdhum)) #scale site covariates siteCovs(umfndvistdhum) <- scale(siteCovs(umfndvistdhum)) summary(umfndvistdhum) #global model for occupancy and detection occu.mglobndvistdhum <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Hum + NDVIstd + W + R + Post, # occupancy formula second, data = umfndvistdhum) # Summarize summary(occu.mglobndvistdhum) ##dredge with fixed detection model occu_dredge_psindvistdhum <- dredge(global.model = occu.mglobndvistdhum, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psindvistdhum[1:10,] ## NDVI and Bound site_cov_ndvi_bound <- read.csv("site_cov_NDVI_Bound.csv", row.names = "Site") # Build unmarkedFramOccu umf12 <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_ndvi_bound) #scale observation covariates obsCovs(umf12) <- scale(obsCovs(umf12)) #scale site covariates siteCovs(umf12) <- scale(siteCovs(umf12)) summary(umf12) #global model for occupancy and detection occu.mglob13 <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Bound + NDVIave + W + R + Post, # occupancy formula second, data = umf12) # Summarize summary(occu.mglob13) ##dredge with fixed detection model occu_dredge_psi19 <- dredge(global.model = occu.mglob13, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psi19[1:10,] #NDVIstd with Bound site_cov_ndvistd_bound <- read.csv("site_cov_NDVIstd_Bound.csv", row.names = "Site") # Build unmarkedFramOccu umfndvistdbound <- unmarkedFrameOccu( # y is a matrix with observed detection history # (0's and 1's, one row per site, one column per survey) y = as.matrix(detection_history), # obsCovs = observation covariates in a list, # each variable has site rows x survey columns obsCovs = list(effort = effort, H = H, carcass = carcass, R = R, post = post), # siteCovs = dataframe with site rows x column variables siteCovs = site_cov_ndvistd_bound) #scale observation covariates obsCovs(umfndvistdbound) <- scale(obsCovs(umfndvistdbound)) #scale site covariates siteCovs(umfndvistdbound) <- scale(siteCovs(umfndvistdbound)) summary(umfndvistdbound) #global model for occupancy and detection occu.mglobndvistdbound <- occu(formula = ~effort + H + carcass + post + R # detection formula first ~Bound + NDVIstd + W + R + Post, # occupancy formula second, data = umfndvistdbound) # Summarize summary(occu.mglobndvistdbound) ##dredge with fixed detection model occu_dredge_psindvistdbound <- dredge(global.model = occu.mglobndvistdbound, rank = "AICc", subset=`p(effort)`&`p(carcass)`&`p(H)`&`p(post)`&`p(R)`) occu_dredge_psindvistdbound[1:10,] #top model list occu.m1 <- occu(formula = ~effort + carcass + H + post + R ~Carcass + C + R, data = umf7) occu.m2 <- occu(formula = ~effort + carcass + H + post + R ~Carcass + C + W, data = umf7) occu.m3 <- occu(formula = ~effort + carcass + H + post + R ~Carcass + C + W + R, data = umf7) AICc(occu.m1) AICc(occu.m2) AICc(occu.m3) coef(occu.m1) coef(occu.m2) coef(occu.m3) occu_model_list <- list(occ_m1 = occu.m1, occ_m2 = occu.m2, occ_m3 = occu.m3) ## ----modelaveragedpredictions---------- occu_modavg_psi_predict <- modavgPred(occu_model_list, # c.hat = 1, # to change variance inflation factor, default = 1) parm.type = "psi", # psi = occupancy, can also be "det" for detection probability newdata = umf14@siteCovs)[c("mod.avg.pred", "lower.CL", "upper.CL")] ## Put predictions, CI, and all site covariates into one data frame occu_modavg_psi_predict_df <- data.frame(Predicted = occu_modavg_psi_predict$mod.avg.pred, lower = occu_modavg_psi_predict$lower.CL, upper = occu_modavg_psi_predict$upper.CL, site_cov_c) head(occu_modavg_psi_predict_df) ## ----modelaveraging_covariates--------------------- # First, set-up a new dataframe to predict along a sequence of the covariate. # Predicting requires all covariates, so let's hold the other covariates constant at their mean value occu_carcass_newdata <- data.frame(Carcass = seq(min(site_cov_c$Carcass), max(site_cov_c$Carcass), by = 0.001), C = mean(site_cov_c$C), # hold other variables constant R = mean(site_cov_c$R), W = mean(site_cov_c$W), Post = mean(site_cov_c$Post)) # hold other variables constant # Model-averaged prediction of occupancy and confidence interval occu_carcass_pred <- modavgPred(occu_model_list, # c.hat = # to change variance inflation factor, default = 1) parm.type = "psi", # psi = occupancy newdata = occu_carcass_newdata)[c("mod.avg.pred", "lower.CL", "upper.CL")] # Put prediction, confidence interval, and covariate values together in a data frame occu_carcass_pred_df <- data.frame(Predicted = occu_carcass_pred$mod.avg.pred, lower = occu_carcass_pred$lower.CL, upper = occu_carcass_pred$upper.CL, occu_carcass_newdata) # Plot the relationship occu_carcass_pred_plot <- ggplot(occu_carcass_pred_df, aes(x = Carcass, y = Predicted)) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5, linetype = "dashed") + geom_path(size = 1) + labs(x = "Probability of carcass occurrence (standardized)", y = "Estimated site use probability") + theme_classic() + coord_cartesian(ylim = c(0,1)) + theme(text = element_text(family = "HelveticaNeue", colour = "black"), axis.text = element_text(colour = "black")) occu_carcass_pred_plot #Plot influence of carcass on site use probability in the original scale UN.CPorig<- seq(min(site_cov_c$Carcass),max(site_cov_c$Carcass), length.out = 139) ##use datafile with original scale UN.CPpred<- (UN.CPorig-mean(site_cov_c$Carcass))/sd(site_cov_c$Carcass) CPnewdata<- data.frame(Carcass = UN.CPpred, R = 0, C = 0, W = 0, Post = 0) ##use mean of standardized values for the other covariates in your model UN.CP <- predict(occu.m1, type="state", newdata=CPnewdata, appendData = TRUE) plot(UN.CPorig, UN.CP[,"Predicted"], xlab="Elephant carcass occurrence probability", ylab="Estimated probability of site use", type="l", ylim = c(0,1), frame = F, lwd = 2) matlines(UN.CPorig, UN.CP[,3:4], lty = 1, col = "grey", lwd = 1) #C # Predicting requires all covariates, so let's hold the other covariates constant at their mean value occu_c_newdata <- data.frame(C = seq(min(site_cov_c$C), max(site_cov_c$C), by = 0.5), Carcass = mean(site_cov_c$Carcass), # hold other variables constant R = mean(site_cov_c$R), W = mean(site_cov_c$W), Post = mean(site_cov_c$Post)) # hold other variables constant # Model-averaged prediction of occupancy and confidence interval occu_c_pred <- modavgPred(occu_model_list, # c.hat = # to change variance inflation factor, default = 1) parm.type = "psi", # psi = occupancy newdata = occu_c_newdata)[c("mod.avg.pred", "lower.CL", "upper.CL")] # Put prediction, confidence interval, and covariate values together in a data frame occu_c_pred_df <- data.frame(Predicted = occu_c_pred$mod.avg.pred, lower = occu_c_pred$lower.CL, upper = occu_c_pred$upper.CL, occu_c_newdata) # Plot the relationship occu_c_pred_plot <- ggplot(occu_c_pred_df, aes(x = C, y = Predicted)) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5, linetype = "dashed") + geom_path(size = 1) + labs(x = "Percentage cover (standardized)", y = "Estimated site use probability") + theme_classic() + coord_cartesian(ylim = c(0,1)) + theme(text = element_text(family = "HelveticaNeue", colour = "black"), axis.text = element_text(colour = "black")) occu_c_pred_plot #Plot influence of population density on site use probability in the original scale UN.CPorig<- seq(min(site_cov_c$C),max(site_cov_c$C), length.out = 139) ##use datafile with original scale UN.CPpred<- (UN.CPorig-mean(site_cov_c$C))/sd(site_cov_c$C) CPnewdata<- data.frame(C = UN.CPpred, R = 0, Carcass = 0, W = 0, Post = 0) ##use mean of standardized values for the other covariates in your model UN.CP <- predict(occu.m1, type="state", newdata=CPnewdata, appendData = TRUE) plot(UN.CPorig, UN.CP[,"Predicted"], xlab="Percentage tree cover", ylab="Estimated probability of site use", type="l", ylim = c(0,1), frame = F, lwd = 2) matlines(UN.CPorig, UN.CP[,3:4], lty = 1, col = "grey", lwd = 1) #R # Predicting requires all covariates, so let's hold the other covariates constant at their mean value occu_r_newdata <- data.frame(R = seq(min(site_cov_c$R), max(site_cov_c$R), by = 0.5), Carcass = mean(site_cov_c$Carcass), # hold other variables constant C = mean(site_cov_c$C), W = mean(site_cov_c$W), Post = mean(site_cov_c$Post)) # hold other variables constant # Model-averaged prediction of occupancy and confidence interval occu_r_pred <- modavgPred(occu_model_list, # c.hat = # to change variance inflation factor, default = 1) parm.type = "psi", # psi = occupancy newdata = occu_r_newdata)[c("mod.avg.pred", "lower.CL", "upper.CL")] # Put prediction, confidence interval, and covariate values together in a data frame occu_r_pred_df <- data.frame(Predicted = occu_r_pred$mod.avg.pred, lower = occu_r_pred$lower.CL, upper = occu_r_pred$upper.CL, occu_r_newdata) # Plot the relationship occu_r_pred_plot <- ggplot(occu_r_pred_df, aes(x = R, y = Predicted)) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5, linetype = "dashed") + geom_path(size = 1) + labs(x = "Distance to river (standardized)", y = "Estimated site use probability") + theme_classic() + coord_cartesian(ylim = c(0,1)) + theme(text = element_text(family = "HelveticaNeue", colour = "black"), axis.text = element_text(colour = "black")) occu_r_pred_plot #Plot influence of river on site use probability in the original scale UN.CPorig<- seq(min(site_cov_c$R),max(site_cov_c$R), length.out = 139) ##use datafile with original scale UN.CPpred<- (UN.CPorig-mean(site_cov_c$R))/sd(site_cov_c$R) CPnewdata<- data.frame(R = UN.CPpred, C = 0, Carcass = 0, W = 0, Post = 0) ##use mean of standardized values for the other covariates in your model UN.CP <- predict(occu.m1, type="state", newdata=CPnewdata, appendData = TRUE) plot(UN.CPorig, UN.CP[,"Predicted"], xlab="Distance to river (km)", ylab="Estimated probability of site use", type="l", ylim = c(0,1), frame = F, lwd = 2) matlines(UN.CPorig, UN.CP[,3:4], lty = 1, col = "grey", lwd = 1) #W # Predicting requires all covariates, so let's hold the other covariates constant at their mean value occu_w_newdata <- data.frame(W = seq(min(site_cov_c$W), max(site_cov_c$W), by = 0.5), Carcass = mean(site_cov_c$Carcass), # hold other variables constant C = mean(site_cov_c$C), R = mean(site_cov_c$R), Post = mean(site_cov_c$Post)) # hold other variables constant # Model-averaged prediction of occupancy and confidence interval occu_w_pred <- modavgPred(occu_model_list, # c.hat = # to change variance inflation factor, default = 1) parm.type = "psi", # psi = occupancy newdata = occu_w_newdata)[c("mod.avg.pred", "lower.CL", "upper.CL")] # Put prediction, confidence interval, and covariate values together in a data frame occu_w_pred_df <- data.frame(Predicted = occu_w_pred$mod.avg.pred, lower = occu_w_pred$lower.CL, upper = occu_w_pred$upper.CL, occu_w_newdata) # Plot the relationship occu_w_pred_plot <- ggplot(occu_w_pred_df, aes(x = W, y = Predicted)) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5, linetype = "dashed") + geom_path(size = 1) + labs(x = "Water availability (standardized)", y = "Estimated site use probability") + theme_classic() + coord_cartesian(ylim = c(0,1)) + theme(text = element_text(family = "HelveticaNeue", colour = "black"), axis.text = element_text(colour = "black")) occu_w_pred_plot #Plot influence of water on site use probability in the original scale UN.CPorig<- seq(min(site_cov_c$W),max(site_cov_c$W), length.out = 139) ##use datafile with original scale UN.CPpred<- (UN.CPorig-mean(site_cov_c$W))/sd(site_cov_c$W) CPnewdata<- data.frame(W = UN.CPpred, C = 0, Carcass = 0, R = 0, Post = 0) ##use mean of standardized values for the other covariates in your model UN.CP <- predict(occu.m2, type="state", newdata=CPnewdata, appendData = TRUE) plot(UN.CPorig, UN.CP[,"Predicted"], xlab="Water availability", ylab="Estimated probability of site use", type="l", ylim = c(0,1), frame = F, lwd = 2) matlines(UN.CPorig, UN.CP[,3:4], lty = 1, col = "grey", lwd = 1) #request confidence intervals using asympotic normal confint(occu.m1, type='det') confint(occu.m1, type='state') #generate PAO and 95% CI s<-nrow(detection_history) re <- ranef(occu.m1) EBUP <- bup(re, stat="mode") CI <- confint(re, level=0.95) print("95 % EB interval on number sites occupied") rbind(s_occup = c(Estimate = sum(EBUP), colSums(CI)) ) print("95 % EB interval on proportion of sites occupied") rbind(PAO = c(Estimate = sum(EBUP), colSums(CI))/s ) #extract estimated coefficients coef(occu.m1) coef(occu.m1, type = "state") coef(occu.m1, type = "det") #To check which types are available for a model, use the names method. names(occu.m1) #to extract the covariance matrix of the estimates vcov(occu.m1, type = "det") vcov(occu.m1, type = "state") coef(occu.m2, type = "state") coef(occu.m2, type = "det") confint(occu.m2, type='det') confint(occu.m2, type='state') coef(occu.m3, type = "state") coef(occu.m3, type = "det") confint(occu.m3, type='det') confint(occu.m3, type='state') #predictions of p for each site p_top <- predict(occu.m1,type="det") #predictions of lambda for each site psi_top <- predict(occu.m1,type="state") # Do Mackenzie-Bailey goodness of fit test for single-season occupancy model occum1_mb.gof.boot <- mb.gof.test(occu.m1, # Demonstrate with small number of sims (10), # but then change to large number (e.g. 1000) nsim = 1000) # View Results occum1_mb.gof.boot #estimate c-hat: 1.06, p=value is 0.35 # Do Mackenzie-Bailey goodness of fit test for single-season occupancy model occumglob_mb.gof.boot <- mb.gof.test(occu.mglob, # Demonstrate with small number of sims (10), # but then change to large number (e.g. 1000) nsim = 1000) # View Results occumglob_mb.gof.boot #estimate c-hat: 1.06, p-value:0.346 occu.mglob