We first import the wolf kill sites and hunt-kills

getwd()
## [1] "C:/Users/6040901/OneDrive - Høgskolen i Innlandet/Dokumenter/New analysis II"
library(rgdal);library(amt);library(lme4);library(performance);library(MuMIn);library(dplyr); library(car);library(ggplot2);library(gridExtra);library(sjPlot);library(sjlabelled); library(snakecase);library(ggeffects);library(png);library(grid);library(forestmodel); library(raster);library(sf);library(sp);library(broom);library(ggsn);library(broom)

wolf.data<-read.csv("new_wolf3006.csv", sep = ",") 
hunt.data<-read.csv("hunt_sites.csv", sep = ",") 
study_area<-readOGR("study_area.shp") #This corresponds to the two wolf territories
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\6040901\OneDrive - Høgskolen i Innlandet\Dokumenter\New analysis II\study_area.shp", layer: "study_area"
## with 1 features
## It has 2 fields
wolf.data$actual_date_1<-as.POSIXct(wolf.data$actual_date_1, "%d/%m/%Y %H:%M", tz="UTC")

wolf.kills<-wolf.data%>% #select only the kill sites 
  filter(kill==1)

wolf.kills.fall<-wolf.kills%>% 
  filter(season=="fall")

wolf.kills.winter<-wolf.kills%>% 
  filter(season=="winter")

#Make a track using the amt package 
wolf_fall_tr<- mk_track(wolf.kills.fall, utmx, utmy, actual_date_1)

We then imported all the raster files that are our covariates

setwd("C:/Users/6040901/OneDrive - Høgskolen i Innlandet/Dokumenter/New Analysis II")

ruggedness.w<-raster("rug_wolf.tif") #terrain ruggedness
distmain.w<-raster("distm_wolf.tif") #distance to main roads
build.w<-raster("build_wolf.tif") #building density
distcl.w<-raster("distcl_wolf.tif") #distance to young forests and clearcuts 
distfor.w<-raster("distfor_wolf.tif") #distance to secondary roads 
moose.w<-raster("moose.tif") #moose density 
wolf.ud.w<-raster("wolf1819.tif") #wolf space use 
distb.w<-raster("distbog_wolf.tif") #distance to bogs 

crs(ruggedness.w)<-"+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"
crs(distmain.w)<-"+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"
crs(distcl.w)<-"+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"
crs(build.w)<-"+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"
crs(wolf.ud.w)<-"+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"
crs(moose.w)<-"+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"
crs(distfor.w)<-"+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"
crs(wolf.ud.w)<-"+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"
crs(distb.w)<-"+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"

#Give the same resolution to all rasters 
distcl.w<-resample(distcl.w,distmain.w, method="bilinear")
build.w<-resample (build.w, distmain.w, method="bilinear")
distfor.w<-resample(distfor.w, distmain.w, method="bilinear")
moose.w<-resample(moose.w, distmain.w, method="bilinear")
wolf.ud.w <- resample(wolf.ud.w,distmain.w, method="bilinear" )
rug.w<-resample(ruggedness.w, distmain.w, method="bilinear")
wolf.ud.w<-resample(wolf.ud.w, distmain.w, method="bilinear")
moose.w<-resample(moose.w, distmain.w, method="bilinear")
distb.w<-resample(distb.w, distmain.w, method="bilinear")

#And the same extent

moose.w<-setExtent(moose.w, distmain.w)
wolf.ud.w<-setExtent(wolf.ud.w, distmain.w)
distfor.w<- setExtent(distfor.w, distmain.w)
distcl.w <- setExtent(distcl.w, distmain.w)
build.w <- setExtent(build.w, distmain.w)
wolf.ud.w<-setExtent(wolf.ud.w, distmain.w)
rug.w<-setExtent(rug.w, distmain.w)
distb.w<-setExtent(distb.w, distmain.w)

Next we explore what number of random points we need to ensure parameter estimate stability.

# Explore sensitivity of HSF coefficients to the number of available points 
n.frac <- c(1, 5, 10, 50) 
n.pts<- ceiling(nrow(wolf.kills.fall)*n.frac) 
n.rep <- 20

res1.wolf <- tibble( 
  n.pts = rep(n.pts, n.rep), 
  frac = rep(n.frac, n.rep), 
  res= map( n.pts, ~ wolf_fall_tr %>% 
              random_points(n = .x) %>% 
              extract_covariates(rug.w) %>% 
              extract_covariates(distmain.w) %>% 
              extract_covariates(build.w) %>% 
              extract_covariates(distcl.w) %>% 
              extract_covariates(distfor.w) %>%
              extract_covariates(distb.w) %>% 
              extract_covariates(moose.w) %>% 
              extract_covariates(wolf.ud.w) %>% 
              mutate(rug = scale(rug_wolf)[, 1], 
                     distm = scale(distm_wolf)[, 1],
                     build = scale(build_wolf) [,1], 
                     distcl = scale(distcl_wolf)[, 1], 
                     distf = scale(distfor_wolf) [,1], 
                     distb = scale(distbog_wolf) [,1], 
                     moose2 = scale(moose) [,1], 
                     wolf = scale(wolf1819) [,1], 
                     w = ifelse(case_, 1, 5000)) %>%
              glm(case_ ~ rug + distm + build + distcl + distf +moose2+distb+ wolf, 
                  weight = w, data = ., family = binomial()) %>% tidy()))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
res1.wolf %>% unnest(cols = res) %>% 
  mutate(kills.fall = dplyr::recode(term, "(Intercept)"="Intercept", 
                                    rug = "Ruggedness", 
                                    build = "Building density", 
                                    distcl = "Distance to clearcut", 
                                    distf = "Distance to forest roads", 
                                    distm = "Distance to main roads", 
                                    moose2="Moose density", 
                                    distb = "Distance to bogs", 
                                    wolf = "Wolf space use" )) %>% 
  ggplot(aes(factor(frac), y = estimate)) + geom_boxplot() + 
  facet_wrap(~ kills.fall, scale ="free") + geom_jitter(alpha = 0.2) + 
  labs(x = "Number of available points (multiplier of no. of used locations)", y = "Estimate") +
  theme_light()

Use the largest sample size here for the rest of the paper

wolf.fall.dataset <- wolf_fall_tr %>% 
  extract_covariates(rug.w) %>% 
  extract_covariates(distmain.w) %>% 
  extract_covariates(build.w) %>% 
  extract_covariates(distcl.w) %>% 
  extract_covariates(distfor.w) %>%
  extract_covariates(distb.w) %>% 
  extract_covariates(moose.w) %>% 
  extract_covariates(wolf.ud.w) %>% 
  mutate(rug = scale(rug_wolf)[, 1], 
         distm = scale(distm_wolf)[, 1], 
         build = scale(build_wolf) [,1], 
         distcl = scale(distcl_wolf)[, 1], 
         distf = scale(distfor_wolf) [,1], 
         distb = scale(distbog_wolf) [,1],
         moose2 = scale(moose) [,1], 
         wolf = scale(wolf1819) [,1])

wolf.fall.dataset$case_<-"TRUE"

wolf.fall.random<-random_points(study_area, n=max(n.pts)) 
wolf.fall.random<- wolf.fall.random %>%
  extract_covariates(rug.w) %>% 
  extract_covariates(distmain.w) %>% 
  extract_covariates(build.w) %>% 
  extract_covariates(distcl.w) %>% 
  extract_covariates(distfor.w) %>%
  extract_covariates(distb.w) %>% 
  extract_covariates(moose.w) %>% 
  extract_covariates(wolf.ud.w) %>% 
  mutate(rug = scale(rug_wolf)[, 1], 
         distm = scale(distm_wolf)[, 1], 
         build = scale(build_wolf) [,1], 
         distcl = scale(distcl_wolf)[, 1], 
         distf = scale(distfor_wolf) [,1], 
         distb = scale(distbog_wolf) [,1], 
         moose2 = scale(moose) [,1], 
         wolf = scale(wolf1819) [,1])

# write_xlsx(wolf.fall.random, "wolf_random_fall.xlsx") # We exported them into arcGIS to check that everything looked okay 
# write_xlsx(wolf.fall.dataset, "wolf_actual_fall.xlsx")

#Re-import the file 
wolf.hsf.fall<-read.csv("wolf_rsf_fall.csv", sep = ",") 
wolf.hsf.fall$w <- ifelse(wolf.hsf.fall$case_, 1, 5000)

wolf.hsf.fall<-wolf.hsf.fall[complete.cases(wolf.hsf.fall$distf),]
wolf.hsf.fall<-wolf.hsf.fall[complete.cases(wolf.hsf.fall$build),]
wolf.hsf.fall<-wolf.hsf.fall[complete.cases(wolf.hsf.fall$distm),]
wolf.hsf.fall<-wolf.hsf.fall[complete.cases(wolf.hsf.fall$rug),]
wolf.hsf.fall<-wolf.hsf.fall[complete.cases(wolf.hsf.fall$moose),]
wolf.hsf.fall<-wolf.hsf.fall[complete.cases(wolf.hsf.fall$wolf),]

hsf.wolf.fall <- glm(case_ ~ rug + build + distm+ distf+distb + wolf + distcl +moose, data = wolf.hsf.fall, weight = w, family = binomial(link = "logit"), na.action="na.fail")

dred.wolf_fall<-dredge(hsf.wolf.fall, rank = "AIC") 
summary(get.models(dred.wolf_fall, 1)[[1]])
## 
## Call:
## glm(formula = case_ ~ build + distb + distcl + wolf + 1, family = binomial(link = "logit"), 
##     data = wolf.hsf.fall, weights = w, na.action = "na.fail")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7712  -0.2142  -0.1769  -0.1434   5.2428  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.255e+01  2.052e-01 -61.153  < 2e-16 ***
## build       -1.775e-01  8.051e-02  -2.204  0.02752 *  
## distb        1.386e-03  5.004e-04   2.769  0.00562 ** 
## distcl      -3.564e-04  1.446e-04  -2.466  0.01367 *  
## wolf         6.174e+02  9.722e+01   6.351 2.15e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2685.2  on 5084  degrees of freedom
## Residual deviance: 2626.1  on 5080  degrees of freedom
## AIC: 2636.1
## 
## Number of Fisher Scoring iterations: 9

Next we plot the relative risk of wolf predation within our study area following the method used by Kauffmann et al. (2007)

#First we import the average raster for each of our covariate variables 
rug.wolf.av.hunting19<-raster("av_rug2_wolf.tif")
build.wolf.av.hunting19<-raster("av_build_wolf.tif")
distf.wolf.av.hunting19<-raster("av_distf_wolf.tif")
distcl.wolf.av.hunting19<-raster("av_distcl_wolf.tif")
wolf.ud.av.hunting19<-raster("av_wolf_hunt19.tif")
moose.av.hunting19<-raster("av_moose19_wolf.tif")
distm.av.hunting19<-raster("av_distm_wolf.tif")
distb.av.hunting19<-raster("av_distbog_wolf.tif")

distcl.wolf.av.hunting19<-resample(distcl.wolf.av.hunting19, distmain.w, method="bilinear")
build.wolf.av.hunting19<-resample(build.wolf.av.hunting19, distmain.w, method="bilinear")
moose.wolf.av.hunting19<-resample(moose.av.hunting19, distmain.w, method="bilinear")
distf.wolf.av.hunting19<-resample(distf.wolf.av.hunting19, distmain.w, method="bilinear")
wolf.ud.av.hunting19<-resample(wolf.ud.av.hunting19, distmain.w, method="bilinear")
distm.wolf.av.hunting19<-resample(distm.av.hunting19,distmain.w, method="bilinear" )
rug.wolf.av.hunting19<-resample(rug.wolf.av.hunting19,distmain.w, method="bilinear" )
distb.wolf.av.hunting19<-resample(distb.av.hunting19, distmain.w, method="bilinear")

distf.wolf.av.hunting19 <- setExtent(distf.wolf.av.hunting19, distmain.w)
distcl.wolf.av.hunting19 <- setExtent(distcl.wolf.av.hunting19, distmain.w)
moose.wolf.av.hunting19<-setExtent(moose.wolf.av.hunting19, distmain.w)
build.wolf.av.hunting19<-setExtent(build.wolf.av.hunting19, distmain.w)
wolf.ud.av.hunting19<-setExtent(wolf.ud.av.hunting19, distmain.w)
distm.wolf.av.hunting19<-setExtent(distm.wolf.av.hunting19,distmain.w)
rug.wolf.av.hunting19<-setExtent(rug.wolf.av.hunting19,distmain.w)
distb.wolf.av.hunting19<-setExtent(distb.wolf.av.hunting19, distmain.w)

options(scipen=1000)
summary(get.models(dred.wolf_fall, 1)[[1]])
## 
## Call:
## glm(formula = case_ ~ build + distb + distcl + wolf + 1, family = binomial(link = "logit"), 
##     data = wolf.hsf.fall, weights = w, na.action = "na.fail")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7712  -0.2142  -0.1769  -0.1434   5.2428  
## 
## Coefficients:
##                Estimate  Std. Error z value             Pr(>|z|)    
## (Intercept) -12.5485956   0.2052008 -61.153 < 0.0000000000000002 ***
## build        -0.1774522   0.0805092  -2.204              0.02752 *  
## distb         0.0013859   0.0005004   2.769              0.00562 ** 
## distcl       -0.0003564   0.0001446  -2.466              0.01367 *  
## wolf        617.3961108  97.2192592   6.351       0.000000000215 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2685.2  on 5084  degrees of freedom
## Residual deviance: 2626.1  on 5080  degrees of freedom
## AIC: 2636.1
## 
## Number of Fisher Scoring iterations: 9
wolf.fall.plotting<-glm(case_ ~ build + distb + wolf + distcl +moose, data = wolf.hsf.fall, weight = w, family = binomial(link = "logit"), na.action="na.fail")
summary(wolf.fall.plotting)
## 
## Call:
## glm(formula = case_ ~ build + distb + wolf + distcl + moose, 
##     family = binomial(link = "logit"), data = wolf.hsf.fall, 
##     weights = w, na.action = "na.fail")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7940  -0.2144  -0.1772  -0.1425   5.2686  
## 
## Coefficients:
##                Estimate  Std. Error z value             Pr(>|z|)    
## (Intercept) -12.3099883   0.3521206 -34.960 < 0.0000000000000002 ***
## build        -0.1838005   0.0834654  -2.202              0.02766 *  
## distb         0.0014190   0.0005000   2.838              0.00454 ** 
## wolf        612.3342034  96.7377747   6.330       0.000000000245 ***
## distcl       -0.0003833   0.0001479  -2.591              0.00956 ** 
## moose        -0.4212533   0.5155146  -0.817              0.41384    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2685.2  on 5084  degrees of freedom
## Residual deviance: 2625.4  on 5079  degrees of freedom
## AIC: 2637.4
## 
## Number of Fisher Scoring iterations: 9
wolf.risk.hunting.season19<- exp((612.3342034*(wolf.ud.w-wolf.ud.av.hunting19))+
                    (-0.4212533*(moose.w-moose.wolf.av.hunting19))+
                    (-0.0003833*(distcl.w - distcl.wolf.av.hunting19))+
                    (-0.1838005*(build.w-build.wolf.av.hunting19))+
                    (0.0014190*(distb.w-distb.wolf.av.hunting19)))

plot(wolf.risk.hunting.season19)

We carried out the same analyses with the wolf kills after the hunting season had ended and with the hunter-killed moose locations.