Telemetry data

Map of sampling

library(tidyverse)
library(ggmap) 
library(maps)
library(mapdata)
register_google(key = "")
#library(ggsn) #scalebar

data <- read_csv("data/Carrizo_telemetry_with_shrub_density.csv")
data <- data %>% 
  drop_na(mesohabitat)

#get map rasters
cali <- get_map(location = c(lon = -119.62, lat = 35.11), zoom = 14, source="stamen", maptype= "terrain", crop= FALSE)
#cali <- get_map(location = c(lon = -119.62, lat = 35.11), zoom = 14, source="stamen", maptype= "toner", crop= FALSE)
p <-ggmap(cali)
#cali <- get_map(location = c(lon = -119.6, lat = 35.11), zoom = 14, source="google", maptype= "terrain", crop= FALSE)
#cali <- get_map(location = c(lon = -119.62, lat = 35.11), zoom = 14, source="google", maptype= "satellite", crop= FALSE)
cali <- get_map(location = c(lon = -119.62, lat = 35.11), zoom = 14, source="google", maptype= "hybrid", crop= FALSE)
q <-ggmap(cali)

#add points
p + geom_point(data = data, aes(x=long, y=lat, color = mesohabitat), alpha = .25, size = 0.5) + 
  scale_color_brewer(type = 'seq', palette = "Set1") +
  labs(x = "latitude", y = "longitude", color = "")

q + geom_point(data = data, aes(x=long, y=lat, color = mesohabitat), alpha = .25, size = 0.5) + 
  scale_color_brewer(type = 'seq', palette = "Set1") +
  labs(x = "latitude", y = "longitude", color = "")

meso_data <- data %>%
  group_by(time.class, ground, mesohabitat) %>%
  summarise(n = n(), density = mean(shrub.density.20))
meso_data

Telemetry stats

library(tidyverse)
data <- read_csv("data/Carrizo_telemetry_with_shrub_density.csv")
data <- data %>% 
  drop_na(mesohabitat)

meso_data <- data %>%
  group_by(time.class, ground, mesohabitat) %>%
  summarise(n = n(), density = mean(shrub.density.20))
meso_data
## # A tibble: 16 x 5
## # Groups:   time.class, ground [8]
##    time.class ground mesohabitat     n density
##    <chr>      <chr>  <chr>       <int>   <dbl>
##  1 am         above  open           69    3.25
##  2 am         above  shrub          14    6.14
##  3 am         below  open          214    3.14
##  4 am         below  shrub          33    7.27
##  5 AM         above  open          449    4.96
##  6 AM         above  shrub         189    7.42
##  7 AM         below  open          526    5.38
##  8 AM         below  shrub          40    8.38
##  9 pm         above  open           72    2.82
## 10 pm         above  shrub          36    7.94
## 11 pm         below  open          151    2.81
## 12 pm         below  shrub          57    6.81
## 13 PM         above  open          350    5.07
## 14 PM         above  shrub         322    6.53
## 15 PM         below  open          929    4.58
## 16 PM         below  shrub         101    9.44
m1 <- glm(n~mesohabitat*ground*time.class, data = meso_data, family = poisson)
summary(m1)
## 
## Call:
## glm(formula = n ~ mesohabitat * ground * time.class, family = poisson, 
##     data = meso_data)
## 
## Deviance Residuals: 
##  [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 
## Coefficients:
##                                            Estimate Std. Error z value
## (Intercept)                                4.234107   0.120386  35.171
## mesohabitatshrub                          -1.595049   0.293123  -5.442
## groundbelow                                1.131870   0.138440   8.176
## time.classAM                               1.872916   0.129306  14.484
## time.classpm                               0.042560   0.168469   0.253
## time.classPM                               1.623827   0.131719  12.328
## mesohabitatshrub:groundbelow              -0.274419   0.347703  -0.789
## mesohabitatshrub:time.classAM              0.729773   0.305679   2.387
## mesohabitatshrub:time.classpm              0.901902   0.357195   2.525
## mesohabitatshrub:time.classPM              1.511668   0.303124   4.987
## groundbelow:time.classAM                  -0.973591   0.152624  -6.379
## groundbelow:time.classpm                  -0.391256   0.199191  -1.964
## groundbelow:time.classPM                  -0.155694   0.151984  -1.024
## mesohabitatshrub:groundbelow:time.classAM -1.436727   0.394102  -3.646
## mesohabitatshrub:groundbelow:time.classpm -0.006662   0.432123  -0.015
## mesohabitatshrub:groundbelow:time.classPM -1.861187   0.371265  -5.013
##                                           Pr(>|z|)    
## (Intercept)                                < 2e-16 ***
## mesohabitatshrub                          5.28e-08 ***
## groundbelow                               2.94e-16 ***
## time.classAM                               < 2e-16 ***
## time.classpm                              0.800557    
## time.classPM                               < 2e-16 ***
## mesohabitatshrub:groundbelow              0.429975    
## mesohabitatshrub:time.classAM             0.016969 *  
## mesohabitatshrub:time.classpm             0.011571 *  
## mesohabitatshrub:time.classPM             6.13e-07 ***
## groundbelow:time.classAM                  1.78e-10 ***
## groundbelow:time.classpm                  0.049504 *  
## groundbelow:time.classPM                  0.305643    
## mesohabitatshrub:groundbelow:time.classAM 0.000267 ***
## mesohabitatshrub:groundbelow:time.classpm 0.987699    
## mesohabitatshrub:groundbelow:time.classPM 5.36e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3.4561e+03  on 15  degrees of freedom
## Residual deviance: 3.8258e-13  on  0  degrees of freedom
## AIC: 138.49
## 
## Number of Fisher Scoring iterations: 3
anova(m1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##                               Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
## NULL                                             15     3456.1          
## mesohabitat                    1  1154.41        14     2301.6 < 2.2e-16
## ground                         1    85.51        13     2216.1 < 2.2e-16
## time.class                     3  1641.35        10      574.8 < 2.2e-16
## mesohabitat:ground             1   341.82         9      233.0 < 2.2e-16
## mesohabitat:time.class         3    36.37         6      196.6 6.240e-08
## ground:time.class              3   136.30         3       60.3 < 2.2e-16
## mesohabitat:ground:time.class  3    60.29         0        0.0 5.096e-13
##                                  
## NULL                             
## mesohabitat                   ***
## ground                        ***
## time.class                    ***
## mesohabitat:ground            ***
## mesohabitat:time.class        ***
## ground:time.class             ***
## mesohabitat:ground:time.class ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
require(emmeans)
em.1 <- emmeans(m1, pairwise~mesohabitat*ground, adjust="tukey")
em.1
## $emmeans
##  mesohabitat ground   emmean         SE  df asymp.LCL asymp.UCL
##  open        above  5.118932 0.04573426 Inf  5.029295  5.208570
##  shrub       above  4.309719 0.08200722 Inf  4.148988  4.470450
##  open        below  5.870666 0.02986740 Inf  5.812127  5.929205
##  shrub       below  3.960890 0.07191483 Inf  3.819939  4.101840
## 
## Results are averaged over the levels of: time.class 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                    estimate         SE  df z.ratio p.value
##  open,above - shrub,above   0.8092135 0.09389785 Inf   8.618  <.0001
##  open,above - open,below   -0.7517343 0.05462310 Inf -13.762  <.0001
##  open,above - shrub,below   1.1580425 0.08522538 Inf  13.588  <.0001
##  shrub,above - open,below  -1.5609477 0.08727683 Inf -17.885  <.0001
##  shrub,above - shrub,below  0.3488290 0.10907303 Inf   3.198  0.0076
##  open,below - shrub,below   1.9097768 0.07787043 Inf  24.525  <.0001
## 
## Results are averaged over the levels of: time.class 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
plot(em.1, comparison = TRUE)

m2 <- glm(n~mesohabitat*ground*density, data = meso_data, family = poisson)
summary(m2)
## 
## Call:
## glm(formula = n ~ mesohabitat * ground * density, family = poisson, 
##     data = meso_data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -15.247   -3.840   -1.928    2.981   15.327  
## 
## Coefficients:
##                                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)                           1.70395    0.21688   7.856 3.95e-15
## mesohabitatshrub                      4.61052    0.47063   9.796  < 2e-16
## groundbelow                           2.34998    0.24326   9.660  < 2e-16
## density                               0.85246    0.04551  18.732  < 2e-16
## mesohabitatshrub:groundbelow         -7.07670    0.72109  -9.814  < 2e-16
## mesohabitatshrub:density             -1.04944    0.07541 -13.917  < 2e-16
## groundbelow:density                  -0.36550    0.05142  -7.108 1.18e-12
## mesohabitatshrub:groundbelow:density  0.86599    0.10178   8.509  < 2e-16
##                                         
## (Intercept)                          ***
## mesohabitatshrub                     ***
## groundbelow                          ***
## density                              ***
## mesohabitatshrub:groundbelow         ***
## mesohabitatshrub:density             ***
## groundbelow:density                  ***
## mesohabitatshrub:groundbelow:density ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3456.05  on 15  degrees of freedom
## Residual deviance:  902.47  on  8  degrees of freedom
## AIC: 1025
## 
## Number of Fisher Scoring iterations: 6
em.2 <- emmeans(m2, pairwise~mesohabitat*density, adjust="tukey")
em.2
## $emmeans
##  mesohabitat  density   emmean         SE  df asymp.LCL asymp.UCL
##  open        5.746557 6.727436 0.03437817 Inf  6.660056  6.794816
##  shrub       5.746557 4.257148 0.09698186 Inf  4.067067  4.447229
## 
## Results are averaged over the levels of: ground 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                       estimate        SE  df
##  open,5.74655707677096 - shrub,5.74655707677096 2.470288 0.1028948 Inf
##  z.ratio p.value
##   24.008  <.0001
## 
## Results are averaged over the levels of: ground 
## Results are given on the log (not the response) scale.
plot(em.2, comparison = TRUE)

#more appropriate estimate
ggplot(meso_data, aes(density, n, color = ground)) +
  stat_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1) +
  xlim(0, 10) +
  ylim(0,1000) +
  scale_color_brewer(palette = "Set1") +
  labs(color = "")

Resource selection function at the population level.

Analyses conducted at 95% MCP for habitat availability. 100% calculations are available in Scripts/RSF.R

#import relocations
reloc <- st_read("GIS/Lizard_Points/telemetryexportPR.shp")
## Reading layer `telemetryexportPR' from data source `F:\School\Carrizo_tel\GIS\Lizard_Points\telemetryexportPR.shp' using driver `ESRI Shapefile'
## Simple feature collection with 3545 features and 46 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 259157.3 ymin: 3888238 xmax: 262883 ymax: 3890634
## epsg (SRID):    26911
## proj4string:    +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
#path <- file.path("https://ndownloader.figshare.com/files/15042206") #this kind of works
#reloc <- st_read(path)
reloc$id <- "id"

#I want to make new datasets that include the survey dates
grid20 <- st_read("GIS/Site4/Fishnet20m.shp")
## Reading layer `Fishnet20m' from data source `F:\School\Carrizo_tel\GIS\Site4\Fishnet20m.shp' using driver `ESRI Shapefile'
## Simple feature collection with 81008 features and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 256931.5 ymin: 3886651 xmax: 263571.5 ymax: 3891531
## epsg (SRID):    26911
## proj4string:    +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
#plot(grid20)
j <- st_join(grid20, reloc)
j <- dplyr::select(j, uniID, relocation, lizard, year, date, time, 14:17,behavior)
data20 <- read.csv("data/mcp20m.csv")

j2 <- left_join(j, data20, by = "uniID")

#need to filter out observations outside the mcp
j2 <- filter(j2, !is.na(x_utm))
ck <- j2 %>% group_by(lizard, year) %>% count()
## Warning: Factor `lizard` contains implicit NA, consider using
## `forcats::fct_explicit_na`

## Warning: Factor `lizard` contains implicit NA, consider using
## `forcats::fct_explicit_na`

## Warning: Factor `lizard` contains implicit NA, consider using
## `forcats::fct_explicit_na`

## Warning: Factor `lizard` contains implicit NA, consider using
## `forcats::fct_explicit_na`

## Warning: Factor `lizard` contains implicit NA, consider using
## `forcats::fct_explicit_na`
j2 <- mutate(j2, aspect.cat = ifelse(Aspect >= 0 & Aspect < 22.5, "north", ifelse(Aspect >= 22.5 & Aspect < 67.5, "northeast", ifelse(Aspect >= 67.5 & Aspect < 112.5, "east", ifelse(Aspect >= 112.5 & Aspect < 157.5, "southeast", ifelse(Aspect >= 157.5 & Aspect < 202.5, "south", ifelse(Aspect >= 202.5 & Aspect < 247.5, "southwest", ifelse(Aspect >= 247.5 & Aspect < 292.5, "west", ifelse(Aspect >= 292.5 & Aspect < 337.5, "northwest", ifelse(Aspect >= 337.5 & Aspect < 360, "north", "flat"))))))))))

#separate datasets, used/available

#For each relocation, choose x available sites, add lizard name, rbind to available set. because used sample was done iwth replacement, need to include used and unused sites.

#functions

#this is the RSF function that takes the dataframe and selects random paired points. nu can be multiplied by any integer to change observed to random ratio. Y is ratio of available to used. 

rsf <- function(x, y){
  us <- filter(x, SumAllYr > 0) #filter cells with lizards
  av <- filter(j2, is.na(lizard))#filter cells with no lizards from full dataset
 # z <- us[!duplicated(us$uniID),]#remove duplicate lizards cells
  all <- rbind(us, av)
  #print(sum(all$SumAllYr))
  counts <- us %>% group_by(lizard, year) %>% count()
  counts <- as.data.frame(counts)
  empty <- all[FALSE,]
  for (i in 1:nrow(counts)) {
    nu <- counts[i, "n"]
    liz <- counts[i, "lizard"]
    yr <- counts[i, "year"]
    size <- 1 * nu
    rows <- all[sample(1:nrow(all), size, replace = TRUE),]
    rows$lizard <- liz
    rows$year <- yr
    empty <- rbind(empty, rows)
}
  empty$pres <- 0
  us$pres <- 1
  s_mcp95 <- rbind(empty, us)
  print(sum(s_mcp95$pres))
  s_mcp95
}

These boxplots compare the lizard use points (pres = 1) to randomly selected points across the MCP

#run RSF on telemetry points falling within 95% MCP

mcp95 <- rsf(j2, 1)
## [1] 3354
sum(mcp95$pres)
## [1] 3354
mcp95 %>% group_by(mesohabita) %>% summarise(sum = sum(pres))
## Warning: Factor `mesohabita` contains implicit NA, consider using
## `forcats::fct_explicit_na`

## Warning: Factor `mesohabita` contains implicit NA, consider using
## `forcats::fct_explicit_na`
## Simple feature collection with 3 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 259491.5 ymin: 3888431 xmax: 262271.5 ymax: 3890511
## epsg (SRID):    26911
## proj4string:    +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
## # A tibble: 3 x 3
##   mesohabita   sum                                                 geometry
##   <fct>      <dbl>                                       <MULTIPOLYGON [m]>
## 1 open        2606 (((260371.5 3888431, 260371.5 3888451, 260391.5 3888451~
## 2 shrub        747 (((260031.5 3888611, 260031.5 3888631, 260051.5 3888631~
## 3 <NA>           1 (((260331.5 3888451, 260331.5 3888471, 260351.5 3888471~
#ggplot(mcp95, aes(as.factor(pres), Shrub_Cov)) + geom_boxplot()+ theme_Publication() + xlab("Lizard Presence") + #ylab("Percent Shrub Cover")

#ggplot(mcp95, aes(as.factor(pres), Elev)) + geom_boxplot() + theme_Publication()+ xlab("Lizard Presence") + #ylab("Elevation")

#ggplot(mcp95, aes(as.factor(pres), Slope)) + geom_boxplot()+ theme_Publication() + xlab("Lizard Presence") + #ylab("Slope")

#ggplot(mcp95, aes(as.factor(pres), fill = aspect.cat)) + geom_bar(position = "dodge")+ theme_Publication()+ #xlab("Lizard Presence") + ylab("Aspect") + theme(legend.title = element_blank())

Modelling

These models use environmental covariates to predict the probility of lizard presence

m1 <- glmmTMB(pres ~ Shrub_Cov + (1|lizard), data = mcp95, family = binomial())
summary(m1)
##  Family: binomial  ( logit )
## Formula:          pres ~ Shrub_Cov + (1 | lizard)
## Data: mcp95
## 
##      AIC      BIC   logLik deviance df.resid 
##   9289.3   9309.8  -4641.7   9283.3     6705 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev.
##  lizard (Intercept) 5.904e-10 2.43e-05
## Number of obs: 6708, groups:  lizard, 78
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.09000    0.03329  -2.703  0.00687 ** 
## Shrub_Cov    3.94334    0.99332   3.970 7.19e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#results the same but variance of random effects is effectively 0

m1.1 <- glmmTMB(pres ~ Shrub_Cov + Elev + Slope + aspect.cat + (1|lizard), data = mcp95, family = binomial())
summary(m1.1)
##  Family: binomial  ( logit )
## Formula:          
## pres ~ Shrub_Cov + Elev + Slope + aspect.cat + (1 | lizard)
## Data: mcp95
## 
##      AIC      BIC   logLik deviance df.resid 
##   8995.1   9083.7  -4484.6   8969.1     6695 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  lizard (Intercept) 1.153e-09 3.395e-05
## Number of obs: 6708, groups:  lizard, 78
## 
## Conditional model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         12.217919   1.622542   7.530 5.07e-14 ***
## Shrub_Cov            3.588422   1.055116   3.401 0.000671 ***
## Elev                -0.016692   0.002353  -7.095 1.30e-12 ***
## Slope               -0.071160   0.017219  -4.133 3.58e-05 ***
## aspect.catflat       0.014296   0.228651   0.063 0.950147    
## aspect.catnorth     -0.975047   0.240405  -4.056 5.00e-05 ***
## aspect.catnortheast -0.140572   0.239884  -0.586 0.557875    
## aspect.catnorthwest -0.240329   0.172449  -1.394 0.163430    
## aspect.catsouth     -0.449271   0.170458  -2.636 0.008397 ** 
## aspect.catsoutheast -0.125606   0.182538  -0.688 0.491384    
## aspect.catsouthwest -0.240899   0.161982  -1.487 0.136963    
## aspect.catwest       0.119942   0.161151   0.744 0.456707    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(m1.1, type = 2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: pres
##             Chisq Df Pr(>Chisq)    
## Shrub_Cov  11.567  1  0.0006715 ***
## Elev       50.335  1  1.296e-12 ***
## Slope      17.079  1  3.585e-05 ***
## aspect.cat 81.250  8  2.738e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.2 <- glmmTMB(pres ~Shrub_Cov + Slope * Elev + aspect.cat  + (1|lizard), data = mcp95, family = binomial(link=logit))
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')
summary(m1.2)
##  Family: binomial  ( logit )
## Formula:          
## pres ~ Shrub_Cov + Slope * Elev + aspect.cat + (1 | lizard)
## Data: mcp95
## 
##      AIC      BIC   logLik deviance df.resid 
##       NA       NA       NA       NA     6694 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  lizard (Intercept) 1.005e-09 3.171e-05
## Number of obs: 6708, groups:  lizard, 78
## 
## Conditional model:
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)         -14.20672         NA      NA       NA
## Shrub_Cov             3.92338         NA      NA       NA
## Slope                 7.86481         NA      NA       NA
## Elev                  0.02019         NA      NA       NA
## aspect.catflat        0.26903         NA      NA       NA
## aspect.catnorth      -1.01557         NA      NA       NA
## aspect.catnortheast  -0.23474         NA      NA       NA
## aspect.catnorthwest  -0.43113         NA      NA       NA
## aspect.catsouth      -0.61635         NA      NA       NA
## aspect.catsoutheast  -0.24387         NA      NA       NA
## aspect.catsouthwest  -0.40834         NA      NA       NA
## aspect.catwest       -0.08412         NA      NA       NA
## Slope:Elev           -0.01097         NA      NA       NA
#warning message because random effects variance is estimate at zero
#the estimates are identical in these models so can just use glm?

m1 <- glm(pres ~Shrub_Cov + Slope + Elev + aspect.cat + NDVI + Solar, data = mcp95, family = binomial(link=logit))
summary(m1)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope + Elev + aspect.cat + 
##     NDVI + Solar, family = binomial(link = logit), data = mcp95)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7601  -1.1793   0.1573   1.1080   1.9071  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.532e+00  1.190e+01  -0.381 0.703338    
## Shrub_Cov            4.860e+00  1.081e+00   4.495 6.94e-06 ***
## Slope               -2.915e-02  2.030e-02  -1.436 0.150979    
## Elev                -1.809e-02  1.865e-03  -9.701  < 2e-16 ***
## aspect.catflat       8.263e-02  2.288e-01   0.361 0.718002    
## aspect.catnorth     -9.387e-01  2.430e-01  -3.862 0.000112 ***
## aspect.catnortheast -1.207e-01  2.424e-01  -0.498 0.618431    
## aspect.catnorthwest -2.328e-01  1.756e-01  -1.326 0.184948    
## aspect.catsouth     -4.697e-01  1.813e-01  -2.590 0.009596 ** 
## aspect.catsoutheast -1.109e-01  1.866e-01  -0.594 0.552179    
## aspect.catsouthwest -3.467e-01  1.707e-01  -2.031 0.042260 *  
## aspect.catwest       6.522e-02  1.624e-01   0.402 0.687955    
## NDVI                 7.623e+00  1.227e+00   6.211 5.25e-10 ***
## Solar                4.186e-05  3.316e-05   1.263 0.206754    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9299.3  on 6707  degrees of freedom
## Residual deviance: 8927.9  on 6694  degrees of freedom
## AIC: 8955.9
## 
## Number of Fisher Scoring iterations: 4
car::Anova(m1, type = 2)
## Analysis of Deviance Table (Type II tests)
## 
## Response: pres
##            LR Chisq Df Pr(>Chisq)    
## Shrub_Cov    20.255  1  6.778e-06 ***
## Slope         2.064  1     0.1508    
## Elev        102.107  1  < 2.2e-16 ***
## aspect.cat   74.193  8  7.154e-13 ***
## NDVI         40.436  1  2.032e-10 ***
## Solar         1.656  1     0.1981    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2 <- glm(pres ~ Shrub_Cov + Elev * Slope + aspect.cat + NDVI + Solar, data = mcp95, family = binomial(link=logit))
summary(m2)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Elev * Slope + aspect.cat + 
##     NDVI + Solar, family = binomial(link = logit), data = mcp95)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9699  -1.1718   0.3166   1.1009   3.2077  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          1.647e+01  1.520e+01   1.084  0.27850    
## Shrub_Cov            5.143e+00  1.117e+00   4.606 4.11e-06 ***
## Elev                 2.294e-02  3.920e-03   5.853 4.83e-09 ***
## Slope                7.720e+00  6.965e-01  11.084  < 2e-16 ***
## aspect.catflat       3.178e-01  2.300e-01   1.382  0.16698    
## aspect.catnorth     -1.112e+00  2.455e-01  -4.528 5.95e-06 ***
## aspect.catnortheast -3.269e-01  2.437e-01  -1.342  0.17973    
## aspect.catnorthwest -5.328e-01  1.786e-01  -2.982  0.00286 ** 
## aspect.catsouth     -4.136e-01  1.868e-01  -2.214  0.02681 *  
## aspect.catsoutheast -9.809e-02  1.891e-01  -0.519  0.60391    
## aspect.catsouthwest -3.263e-01  1.742e-01  -1.873  0.06114 .  
## aspect.catwest      -1.162e-01  1.636e-01  -0.710  0.47745    
## NDVI                 7.572e+00  1.267e+00   5.976 2.29e-09 ***
## Solar               -9.051e-05  4.301e-05  -2.104  0.03534 *  
## Elev:Slope          -1.076e-02  9.682e-04 -11.111  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9299.3  on 6707  degrees of freedom
## Residual deviance: 8765.1  on 6693  degrees of freedom
## AIC: 8795.1
## 
## Number of Fisher Scoring iterations: 5
car::Anova(m2, type =3)
## Analysis of Deviance Table (Type III tests)
## 
## Response: pres
##            LR Chisq Df Pr(>Chisq)    
## Shrub_Cov    21.406  1  3.716e-06 ***
## Elev         35.480  1  2.577e-09 ***
## Slope       161.106  1  < 2.2e-16 ***
## aspect.cat   69.667  8  5.726e-12 ***
## NDVI         37.324  1  1.001e-09 ***
## Solar         4.403  1    0.03587 *  
## Elev:Slope  162.792  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- glm(pres ~ Shrub_Cov * Slope + Elev + NDVI + Solar, data = mcp95, family = binomial(link=logit))
summary(m3)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov * Slope + Elev + NDVI + Solar, 
##     family = binomial(link = logit), data = mcp95)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8828  -1.1889   0.3014   1.1235   1.7646  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      2.317e+01  7.794e+00   2.973  0.00295 ** 
## Shrub_Cov        9.761e+00  1.941e+00   5.028 4.95e-07 ***
## Slope           -6.409e-03  1.990e-02  -0.322  0.74734    
## Elev            -1.579e-02  1.699e-03  -9.292  < 2e-16 ***
## NDVI             8.673e+00  1.206e+00   7.189 6.52e-13 ***
## Solar           -3.661e-05  2.176e-05  -1.682  0.09251 .  
## Shrub_Cov:Slope -1.560e+00  5.103e-01  -3.058  0.00223 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9299.3  on 6707  degrees of freedom
## Residual deviance: 8989.2  on 6701  degrees of freedom
## AIC: 9003.2
## 
## Number of Fisher Scoring iterations: 4
car::Anova(m3, type = 3)
## Analysis of Deviance Table (Type III tests)
## 
## Response: pres
##                 LR Chisq Df Pr(>Chisq)    
## Shrub_Cov         30.388  1  3.537e-08 ***
## Slope              0.104  1  0.7476179    
## Elev              88.958  1  < 2.2e-16 ***
## NDVI              54.741  1  1.375e-13 ***
## Solar              2.797  1  0.0944669 .  
## Shrub_Cov:Slope   12.836  1  0.0003401 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m4 <- glm(pres ~ Slope * Elev + NDVI + Solar*Shrub_Cov + aspect.cat, data = mcp95, family = binomial(link=logit))
summary(m4)
## 
## Call:
## glm(formula = pres ~ Slope * Elev + NDVI + Solar * Shrub_Cov + 
##     aspect.cat, family = binomial(link = logit), data = mcp95)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9072  -1.1762   0.3142   1.0997   3.2572  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          3.916e+01  1.626e+01   2.408 0.016046 *  
## Slope                8.019e+00  7.104e-01  11.288  < 2e-16 ***
## Elev                 2.381e-02  3.948e-03   6.032 1.62e-09 ***
## NDVI                 7.520e+00  1.281e+00   5.873 4.29e-09 ***
## Solar               -1.518e-04  4.577e-05  -3.317 0.000911 ***
## Shrub_Cov           -1.180e+03  3.003e+02  -3.930 8.50e-05 ***
## aspect.catflat       3.198e-01  2.301e-01   1.389 0.164693    
## aspect.catnorth     -1.069e+00  2.448e-01  -4.369 1.25e-05 ***
## aspect.catnortheast -2.509e-01  2.441e-01  -1.028 0.304128    
## aspect.catnorthwest -4.690e-01  1.793e-01  -2.616 0.008906 ** 
## aspect.catsouth     -4.075e-01  1.870e-01  -2.179 0.029329 *  
## aspect.catsoutheast -8.238e-02  1.894e-01  -0.435 0.663523    
## aspect.catsouthwest -3.181e-01  1.744e-01  -1.824 0.068176 .  
## aspect.catwest      -1.176e-01  1.637e-01  -0.718 0.472655    
## Slope:Elev          -1.117e-02  9.877e-04 -11.314  < 2e-16 ***
## Solar:Shrub_Cov      3.118e-03  7.900e-04   3.947 7.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9299.3  on 6707  degrees of freedom
## Residual deviance: 8749.6  on 6692  degrees of freedom
## AIC: 8781.6
## 
## Number of Fisher Scoring iterations: 5
car::Anova(m4, type = 3)
## Analysis of Deviance Table (Type III tests)
## 
## Response: pres
##                 LR Chisq Df Pr(>Chisq)    
## Slope            161.995  1  < 2.2e-16 ***
## Elev              37.690  1  8.291e-10 ***
## NDVI              36.058  1  1.915e-09 ***
## Solar             10.954  1   0.000934 ***
## Shrub_Cov         15.394  1  8.728e-05 ***
## aspect.cat        61.315  8  2.571e-10 ***
## Slope:Elev       163.848  1  < 2.2e-16 ***
## Solar:Shrub_Cov   15.527  1  8.134e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m5 <- glm(pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar*aspect.cat, data = mcp95, family = binomial(link=logit))
summary(m5)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar * 
##     aspect.cat, family = binomial(link = logit), data = mcp95)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0997  -1.1701   0.2736   1.1100   3.2572  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                6.909e+01  9.337e+01   0.740   0.4593    
## Shrub_Cov                  5.521e+00  1.128e+00   4.896 9.79e-07 ***
## Slope                      7.439e+00  7.202e-01  10.330  < 2e-16 ***
## Elev                       2.054e-02  4.133e-03   4.971 6.67e-07 ***
## NDVI                       8.301e+00  1.291e+00   6.429 1.28e-10 ***
## Solar                     -2.250e-04  2.468e-04  -0.911   0.3621    
## aspect.catflat            -3.490e+02  2.091e+02  -1.669   0.0951 .  
## aspect.catnorth           -3.063e+02  1.214e+02  -2.523   0.0116 *  
## aspect.catnortheast       -1.663e+02  1.297e+02  -1.283   0.1995    
## aspect.catnorthwest       -7.831e+01  9.671e+01  -0.810   0.4181    
## aspect.catsouth           -1.306e+02  1.010e+02  -1.294   0.1958    
## aspect.catsoutheast        1.129e+01  1.047e+02   0.108   0.9142    
## aspect.catsouthwest       -7.085e+01  9.647e+01  -0.735   0.4626    
## aspect.catwest            -1.080e+01  9.387e+01  -0.115   0.9084    
## Slope:Elev                -1.039e-02  1.000e-03 -10.383  < 2e-16 ***
## Solar:aspect.catflat       9.187e-04  5.500e-04   1.670   0.0949 .  
## Solar:aspect.catnorth      8.050e-04  3.199e-04   2.516   0.0119 *  
## Solar:aspect.catnortheast  4.376e-04  3.417e-04   1.281   0.2003    
## Solar:aspect.catnorthwest  2.050e-04  2.547e-04   0.805   0.4208    
## Solar:aspect.catsouth      3.419e-04  2.657e-04   1.287   0.1981    
## Solar:aspect.catsoutheast -2.945e-05  2.757e-04  -0.107   0.9149    
## Solar:aspect.catsouthwest  1.856e-04  2.540e-04   0.731   0.4648    
## Solar:aspect.catwest       2.824e-05  2.472e-04   0.114   0.9091    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9299.3  on 6707  degrees of freedom
## Residual deviance: 8736.0  on 6685  degrees of freedom
## AIC: 8782
## 
## Number of Fisher Scoring iterations: 5
car::Anova(m5, type = 3)
## Analysis of Deviance Table (Type III tests)
## 
## Response: pres
##                  LR Chisq Df Pr(>Chisq)    
## Shrub_Cov          24.211  1  8.633e-07 ***
## Slope             140.852  1  < 2.2e-16 ***
## Elev               25.452  1  4.536e-07 ***
## NDVI               43.309  1  4.673e-11 ***
## Solar               0.746  1  0.3878458    
## aspect.cat         29.288  8  0.0002823 ***
## Slope:Elev        143.388  1  < 2.2e-16 ***
## Solar:aspect.cat   29.118  8  0.0003024 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mnull <- glm(pres ~ 1, data = mcp95, family = binomial(link=logit))

knitr::kable(AIC(m1, m2, m3, m4, m5, mnull))
df AIC
m1 14 8955.878
m2 15 8795.087
m3 7 9003.236
m4 16 8781.560
m5 23 8781.969
mnull 1 9301.263
confint(m5)
## Waiting for profiling to be done...
##                                   2.5 %        97.5 %
## (Intercept)               -1.352511e+02  2.348408e+02
## Shrub_Cov                  3.315726e+00  7.736871e+00
## Slope                      6.053674e+00  8.875471e+00
## Elev                       1.248956e-02  2.869113e-02
## NDVI                       5.790265e+00  1.085355e+01
## Solar                     -6.636102e-04  3.147853e-04
## aspect.catflat            -7.699248e+02  5.616967e+01
## aspect.catnorth           -5.408758e+02 -5.859124e+01
## aspect.catnortheast       -4.142031e+02  9.896793e+01
## aspect.catnorthwest       -2.519884e+02  1.312228e+02
## aspect.catsouth           -3.146591e+02  8.489209e+01
## aspect.catsoutheast       -1.809363e+02  2.338141e+02
## aspect.catsouthwest       -2.446275e+02  1.375026e+02
## aspect.catwest            -1.777705e+02  1.943006e+02
## Slope:Elev                -1.238376e-02 -8.463968e-03
## Solar:aspect.catflat      -1.470970e-04  2.025716e-03
## Solar:aspect.catnorth      1.522980e-04  1.423194e-03
## Solar:aspect.catnortheast -2.614504e-04  1.091008e-03
## Solar:aspect.catnorthwest -3.467807e-04  6.625337e-04
## Solar:aspect.catsouth     -2.252235e-04  8.262152e-04
## Solar:aspect.catsoutheast -6.151370e-04  4.763671e-04
## Solar:aspect.catsouthwest -3.629410e-04  6.431035e-04
## Solar:aspect.catwest      -5.118430e-04  4.679421e-04
knitr::kable(confint(m5))
## Waiting for profiling to be done...
2.5 % 97.5 %
(Intercept) -135.2510532 234.8407921
Shrub_Cov 3.3157262 7.7368710
Slope 6.0536737 8.8754711
Elev 0.0124896 0.0286911
NDVI 5.7902654 10.8535461
Solar -0.0006636 0.0003148
aspect.catflat -769.9247647 56.1696696
aspect.catnorth -540.8757903 -58.5912426
aspect.catnortheast -414.2031343 98.9679275
aspect.catnorthwest -251.9883593 131.2227662
aspect.catsouth -314.6591378 84.8920883
aspect.catsoutheast -180.9362610 233.8140801
aspect.catsouthwest -244.6275141 137.5026466
aspect.catwest -177.7704628 194.3006225
Slope:Elev -0.0123838 -0.0084640
Solar:aspect.catflat -0.0001471 0.0020257
Solar:aspect.catnorth 0.0001523 0.0014232
Solar:aspect.catnortheast -0.0002615 0.0010910
Solar:aspect.catnorthwest -0.0003468 0.0006625
Solar:aspect.catsouth -0.0002252 0.0008262
Solar:aspect.catsoutheast -0.0006151 0.0004764
Solar:aspect.catsouthwest -0.0003629 0.0006431
Solar:aspect.catwest -0.0005118 0.0004679
#Best model is m5

summary(m5)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar * 
##     aspect.cat, family = binomial(link = logit), data = mcp95)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0997  -1.1701   0.2736   1.1100   3.2572  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                6.909e+01  9.337e+01   0.740   0.4593    
## Shrub_Cov                  5.521e+00  1.128e+00   4.896 9.79e-07 ***
## Slope                      7.439e+00  7.202e-01  10.330  < 2e-16 ***
## Elev                       2.054e-02  4.133e-03   4.971 6.67e-07 ***
## NDVI                       8.301e+00  1.291e+00   6.429 1.28e-10 ***
## Solar                     -2.250e-04  2.468e-04  -0.911   0.3621    
## aspect.catflat            -3.490e+02  2.091e+02  -1.669   0.0951 .  
## aspect.catnorth           -3.063e+02  1.214e+02  -2.523   0.0116 *  
## aspect.catnortheast       -1.663e+02  1.297e+02  -1.283   0.1995    
## aspect.catnorthwest       -7.831e+01  9.671e+01  -0.810   0.4181    
## aspect.catsouth           -1.306e+02  1.010e+02  -1.294   0.1958    
## aspect.catsoutheast        1.129e+01  1.047e+02   0.108   0.9142    
## aspect.catsouthwest       -7.085e+01  9.647e+01  -0.735   0.4626    
## aspect.catwest            -1.080e+01  9.387e+01  -0.115   0.9084    
## Slope:Elev                -1.039e-02  1.000e-03 -10.383  < 2e-16 ***
## Solar:aspect.catflat       9.187e-04  5.500e-04   1.670   0.0949 .  
## Solar:aspect.catnorth      8.050e-04  3.199e-04   2.516   0.0119 *  
## Solar:aspect.catnortheast  4.376e-04  3.417e-04   1.281   0.2003    
## Solar:aspect.catnorthwest  2.050e-04  2.547e-04   0.805   0.4208    
## Solar:aspect.catsouth      3.419e-04  2.657e-04   1.287   0.1981    
## Solar:aspect.catsoutheast -2.945e-05  2.757e-04  -0.107   0.9149    
## Solar:aspect.catsouthwest  1.856e-04  2.540e-04   0.731   0.4648    
## Solar:aspect.catwest       2.824e-05  2.472e-04   0.114   0.9091    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9299.3  on 6707  degrees of freedom
## Residual deviance: 8736.0  on 6685  degrees of freedom
## AIC: 8782
## 
## Number of Fisher Scoring iterations: 5
knitr::kable(car::Anova(m5, type = 3))
LR Chisq Df Pr(>Chisq)
Shrub_Cov 24.211208 1 0.0000009
Slope 140.852318 1 0.0000000
Elev 25.451576 1 0.0000005
NDVI 43.309350 1 0.0000000
Solar 0.745685 1 0.3878458
aspect.cat 29.287682 8 0.0002823
Slope:Elev 143.387861 1 0.0000000
Solar:aspect.cat 29.117594 8 0.0003024
library(stargazer)
## Warning: package 'stargazer' was built under R version 3.4.4
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
#write model outputs to tables 

stargazer(confint(m5), type = "html", summary = FALSE, out = "tables/m5confint.doc")
## Waiting for profiling to be done...
## 
## <table style="text-align:center"><tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>2.5 %</td><td>97.5 %</td></tr>
## <tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">(Intercept)</td><td>-135.251</td><td>234.841</td></tr>
## <tr><td style="text-align:left">Shrub_Cov</td><td>3.316</td><td>7.737</td></tr>
## <tr><td style="text-align:left">Slope</td><td>6.054</td><td>8.875</td></tr>
## <tr><td style="text-align:left">Elev</td><td>0.012</td><td>0.029</td></tr>
## <tr><td style="text-align:left">NDVI</td><td>5.790</td><td>10.854</td></tr>
## <tr><td style="text-align:left">Solar</td><td>-0.001</td><td>0.0003</td></tr>
## <tr><td style="text-align:left">aspect.catflat</td><td>-769.925</td><td>56.170</td></tr>
## <tr><td style="text-align:left">aspect.catnorth</td><td>-540.876</td><td>-58.591</td></tr>
## <tr><td style="text-align:left">aspect.catnortheast</td><td>-414.203</td><td>98.968</td></tr>
## <tr><td style="text-align:left">aspect.catnorthwest</td><td>-251.988</td><td>131.223</td></tr>
## <tr><td style="text-align:left">aspect.catsouth</td><td>-314.659</td><td>84.892</td></tr>
## <tr><td style="text-align:left">aspect.catsoutheast</td><td>-180.936</td><td>233.814</td></tr>
## <tr><td style="text-align:left">aspect.catsouthwest</td><td>-244.628</td><td>137.503</td></tr>
## <tr><td style="text-align:left">aspect.catwest</td><td>-177.770</td><td>194.301</td></tr>
## <tr><td style="text-align:left">Slope:Elev</td><td>-0.012</td><td>-0.008</td></tr>
## <tr><td style="text-align:left">Solar:aspect.catflat</td><td>-0.0001</td><td>0.002</td></tr>
## <tr><td style="text-align:left">Solar:aspect.catnorth</td><td>0.0002</td><td>0.001</td></tr>
## <tr><td style="text-align:left">Solar:aspect.catnortheast</td><td>-0.0003</td><td>0.001</td></tr>
## <tr><td style="text-align:left">Solar:aspect.catnorthwest</td><td>-0.0003</td><td>0.001</td></tr>
## <tr><td style="text-align:left">Solar:aspect.catsouth</td><td>-0.0002</td><td>0.001</td></tr>
## <tr><td style="text-align:left">Solar:aspect.catsoutheast</td><td>-0.001</td><td>0.0005</td></tr>
## <tr><td style="text-align:left">Solar:aspect.catsouthwest</td><td>-0.0004</td><td>0.001</td></tr>
## <tr><td style="text-align:left">Solar:aspect.catwest</td><td>-0.001</td><td>0.0005</td></tr>
## <tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr></table>
stargazer(AIC(m1, m2, m3, m4, m5, mnull), type = "html", summary = FALSE, out = "tables/AIC.doc")
## 
## <table style="text-align:center"><tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>df</td><td>AIC</td></tr>
## <tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">m1</td><td>14</td><td>8,955.878</td></tr>
## <tr><td style="text-align:left">m2</td><td>15</td><td>8,795.087</td></tr>
## <tr><td style="text-align:left">m3</td><td>7</td><td>9,003.236</td></tr>
## <tr><td style="text-align:left">m4</td><td>16</td><td>8,781.560</td></tr>
## <tr><td style="text-align:left">m5</td><td>23</td><td>8,781.969</td></tr>
## <tr><td style="text-align:left">mnull</td><td>1</td><td>9,301.263</td></tr>
## <tr><td colspan="3" style="border-bottom: 1px solid black"></td></tr></table>

Model Cross-Validation

Cross-validation to assess model fit

library(boot)
cv.glm(mcp95, m1, K = 10)$delta
## [1] 0.2380534 0.2379921
cv.glm(mcp95, m2, K = 10)$delta
## [1] 0.2332428 0.2331842
cv.glm(mcp95, m3, K = 10)$delta
## [1] 0.2395499 0.2395288
cv.glm(mcp95, m4, K = 10)$delta
## [1] 0.2328606 0.2328015
cv.glm(mcp95, m5, K = 10)$delta
## [1] 0.2324763 0.2323958

Mesohabitat

Create models for shrub and open mesohabitats separately

df_list <- split(j2, j2$mesohabita)
meso <- map2(df_list, 2, rsf)
open <- meso$open
shrub <- meso$shrub

#open mesohabitat
mes1 <- glm(pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar*aspect.cat, data = open, family = binomial(link=logit))
summary(mes1)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar * 
##     aspect.cat, family = binomial(link = logit), data = open)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1747  -1.1631   0.2586   1.1003   3.3110  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                6.888e+01  1.317e+02   0.523  0.60109    
## Shrub_Cov                  3.731e+00  1.342e+00   2.781  0.00542 ** 
## Slope                      8.246e+00  8.474e-01   9.731  < 2e-16 ***
## Elev                       2.739e-02  4.924e-03   5.562 2.67e-08 ***
## NDVI                       7.271e+00  1.460e+00   4.982 6.29e-07 ***
## Solar                     -2.363e-04  3.483e-04  -0.678  0.49750    
## aspect.catflat            -5.502e+02  3.003e+02  -1.832  0.06694 .  
## aspect.catnorth           -4.626e+02  1.782e+02  -2.597  0.00941 ** 
## aspect.catnortheast       -3.040e+02  1.910e+02  -1.591  0.11154    
## aspect.catnorthwest       -6.128e+01  1.340e+02  -0.457  0.64754    
## aspect.catsouth           -9.597e+01  1.370e+02  -0.701  0.48351    
## aspect.catsoutheast       -3.951e+01  1.421e+02  -0.278  0.78092    
## aspect.catsouthwest       -1.517e+01  1.334e+02  -0.114  0.90946    
## aspect.catwest             1.500e+01  1.319e+02   0.114  0.90951    
## Slope:Elev                -1.155e-02  1.179e-03  -9.793  < 2e-16 ***
## Solar:aspect.catflat       1.448e-03  7.900e-04   1.833  0.06676 .  
## Solar:aspect.catnorth      1.217e-03  4.695e-04   2.591  0.00956 ** 
## Solar:aspect.catnortheast  8.010e-04  5.037e-04   1.590  0.11178    
## Solar:aspect.catnorthwest  1.594e-04  3.530e-04   0.452  0.65150    
## Solar:aspect.catsouth      2.512e-04  3.606e-04   0.697  0.48609    
## Solar:aspect.catsoutheast  1.039e-04  3.740e-04   0.278  0.78120    
## Solar:aspect.catsouthwest  3.955e-05  3.513e-04   0.113  0.91036    
## Solar:aspect.catwest      -3.969e-05  3.475e-04  -0.114  0.90907    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 7225.4  on 5211  degrees of freedom
## Residual deviance: 6749.1  on 5189  degrees of freedom
## AIC: 6795.1
## 
## Number of Fisher Scoring iterations: 5
knitr::kable(car::Anova(mes1, type = 3))
LR Chisq Df Pr(>Chisq)
Shrub_Cov 7.7623006 1 0.0053348
Slope 124.4450936 1 0.0000000
Elev 32.1761861 1 0.0000000
NDVI 25.7646475 1 0.0000004
Solar 0.4595259 1 0.4978456
aspect.cat 40.2119470 8 0.0000029
Slope:Elev 126.9397996 1 0.0000000
Solar:aspect.cat 39.9578227 8 0.0000033
#shrub mesohabitat
mes2 <- glm(pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar*aspect.cat, data = shrub, family = binomial(link=logit))
summary(mes2)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar * 
##     aspect.cat, family = binomial(link = logit), data = shrub)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6007  -1.0737   0.1716   1.0738   1.8180  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -1.385e+02  2.881e+02  -0.481   0.6308    
## Shrub_Cov                  2.014e+01  2.434e+00   8.278  < 2e-16 ***
## Slope                      6.502e+00  1.418e+00   4.586 4.52e-06 ***
## Elev                       1.279e-02  8.186e-03   1.563   0.1181    
## NDVI                       6.273e+00  2.846e+00   2.204   0.0275 *  
## Solar                      3.356e-04  7.598e-04   0.442   0.6587    
## aspect.catflat            -4.177e+02  5.335e+02  -0.783   0.4336    
## aspect.catnorth           -2.086e+02  3.309e+02  -0.630   0.5284    
## aspect.catnortheast       -7.053e+01  3.547e+02  -0.199   0.8424    
## aspect.catnorthwest        6.438e+00  2.930e+02   0.022   0.9825    
## aspect.catsouth            3.714e+01  2.973e+02   0.125   0.9006    
## aspect.catsoutheast        3.322e+02  3.309e+02   1.004   0.3154    
## aspect.catsouthwest        3.519e+01  2.919e+02   0.121   0.9041    
## aspect.catwest             1.253e+02  2.895e+02   0.433   0.6651    
## Slope:Elev                -9.060e-03  1.959e-03  -4.624 3.76e-06 ***
## Solar:aspect.catflat       1.099e-03  1.403e-03   0.783   0.4336    
## Solar:aspect.catnorth      5.516e-04  8.715e-04   0.633   0.5268    
## Solar:aspect.catnortheast  1.887e-04  9.340e-04   0.202   0.8399    
## Solar:aspect.catnorthwest -1.679e-05  7.711e-04  -0.022   0.9826    
## Solar:aspect.catsouth     -1.000e-04  7.822e-04  -0.128   0.8982    
## Solar:aspect.catsoutheast -8.727e-04  8.703e-04  -1.003   0.3160    
## Solar:aspect.catsouthwest -9.387e-05  7.683e-04  -0.122   0.9028    
## Solar:aspect.catwest      -3.295e-04  7.620e-04  -0.432   0.6655    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2071.1  on 1493  degrees of freedom
## Residual deviance: 1846.0  on 1471  degrees of freedom
## AIC: 1892
## 
## Number of Fisher Scoring iterations: 5
knitr::kable(car::Anova(mes2, type = 3))
LR Chisq Df Pr(>Chisq)
Shrub_Cov 76.4924014 1 0.0000000
Slope 24.9024213 1 0.0000006
Elev 2.4707692 1 0.1159812
NDVI 5.0164561 1 0.0251075
Solar 0.2171558 1 0.6412156
aspect.cat 11.9657488 8 0.1527386
Slope:Elev 25.5347362 1 0.0000004
Solar:aspect.cat 11.9516447 8 0.1533745

Ground

split by above and below ground

df_list <- split(j2, j2$ground)
ground <- map2(df_list, 2, rsf)
below <- ground$below
above <- ground$above

#below ground
gr1 <- glm(pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar*aspect.cat, data = below, family = binomial(link=logit))
summary(gr1)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar * 
##     aspect.cat, family = binomial(link = logit), data = below)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2263  -1.1560   0.2505   1.1006   3.2445  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.949e+02  1.717e+02   2.300 0.021463 *  
## Shrub_Cov                  7.720e+00  1.614e+00   4.783 1.73e-06 ***
## Slope                      7.144e+00  9.036e-01   7.906 2.65e-15 ***
## Elev                       2.223e-02  5.467e-03   4.065 4.80e-05 ***
## NDVI                       5.583e+00  1.641e+00   3.402 0.000669 ***
## Solar                     -1.085e-03  4.537e-04  -2.391 0.016802 *  
## aspect.catflat            -1.789e+02  2.524e+02  -0.709 0.478486    
## aspect.catnorth           -7.209e+02  2.249e+02  -3.205 0.001350 ** 
## aspect.catnortheast       -6.706e+02  2.377e+02  -2.821 0.004789 ** 
## aspect.catnorthwest       -3.565e+02  1.738e+02  -2.051 0.040255 *  
## aspect.catsouth           -4.537e+02  1.770e+02  -2.563 0.010369 *  
## aspect.catsoutheast       -3.981e+02  1.828e+02  -2.178 0.029422 *  
## aspect.catsouthwest       -4.165e+02  1.731e+02  -2.405 0.016155 *  
## aspect.catwest            -3.132e+02  1.720e+02  -1.821 0.068534 .  
## Slope:Elev                -1.002e-02  1.255e-03  -7.983 1.43e-15 ***
## Solar:aspect.catflat       4.725e-04  6.641e-04   0.711 0.476786    
## Solar:aspect.catnorth      1.895e-03  5.927e-04   3.198 0.001384 ** 
## Solar:aspect.catnortheast  1.767e-03  6.269e-04   2.820 0.004808 ** 
## Solar:aspect.catnorthwest  9.373e-04  4.579e-04   2.047 0.040639 *  
## Solar:aspect.catsouth      1.193e-03  4.661e-04   2.560 0.010469 *  
## Solar:aspect.catsoutheast  1.049e-03  4.813e-04   2.180 0.029250 *  
## Solar:aspect.catsouthwest  1.096e-03  4.560e-04   2.404 0.016229 *  
## Solar:aspect.catwest       8.252e-04  4.530e-04   1.822 0.068509 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5428.7  on 3915  degrees of freedom
## Residual deviance: 5075.3  on 3893  degrees of freedom
## AIC: 5121.3
## 
## Number of Fisher Scoring iterations: 5
knitr::kable(car::Anova(gr1, type = 3))
LR Chisq Df Pr(>Chisq)
Shrub_Cov 23.349781 1 0.0000014
Slope 79.726881 1 0.0000000
Elev 17.039882 1 0.0000366
NDVI 11.835989 1 0.0005810
Solar 6.774784 1 0.0092455
aspect.cat 28.120007 8 0.0004520
Slope:Elev 81.888865 1 0.0000000
Solar:aspect.cat 27.997670 8 0.0004747
#above ground
gr2 <- glm(pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar*aspect.cat, data = above, family = binomial(link=logit))
summary(gr2)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar * 
##     aspect.cat, family = binomial(link = logit), data = above)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1620  -1.1456   0.1957   1.0820   2.2384  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -3.900e+02  1.971e+02  -1.979 0.047869 *  
## Shrub_Cov                  9.610e+00  1.744e+00   5.511 3.56e-08 ***
## Slope                      8.990e+00  1.150e+00   7.818 5.36e-15 ***
## Elev                       2.172e-02  6.431e-03   3.378 0.000731 ***
## NDVI                       1.189e+01  2.114e+00   5.622 1.88e-08 ***
## Solar                      9.791e-04  5.202e-04   1.882 0.059837 .  
## aspect.catflat             1.268e+02  3.103e+02   0.408 0.682944    
## aspect.catnorth           -1.387e+02  2.460e+02  -0.564 0.572924    
## aspect.catnortheast        8.325e+00  2.534e+02   0.033 0.973793    
## aspect.catnorthwest        3.325e+02  2.010e+02   1.655 0.097971 .  
## aspect.catsouth            2.575e+02  2.034e+02   1.266 0.205546    
## aspect.catsoutheast        5.332e+02  2.163e+02   2.465 0.013702 *  
## aspect.catsouthwest        3.856e+02  1.989e+02   1.939 0.052516 .  
## aspect.catwest             4.668e+02  1.976e+02   2.362 0.018166 *  
## Slope:Elev                -1.249e-02  1.599e-03  -7.812 5.64e-15 ***
## Solar:aspect.catflat      -3.357e-04  8.164e-04  -0.411 0.680921    
## Solar:aspect.catnorth      3.668e-04  6.483e-04   0.566 0.571553    
## Solar:aspect.catnortheast -2.032e-05  6.677e-04  -0.030 0.975724    
## Solar:aspect.catnorthwest -8.761e-04  5.291e-04  -1.656 0.097752 .  
## Solar:aspect.catsouth     -6.806e-04  5.354e-04  -1.271 0.203632    
## Solar:aspect.catsoutheast -1.404e-03  5.692e-04  -2.466 0.013663 *  
## Solar:aspect.catsouthwest -1.016e-03  5.235e-04  -1.942 0.052184 .  
## Solar:aspect.catwest      -1.230e-03  5.202e-04  -2.364 0.018101 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3870.5  on 2791  degrees of freedom
## Residual deviance: 3519.4  on 2769  degrees of freedom
## AIC: 3565.4
## 
## Number of Fisher Scoring iterations: 5
knitr::kable(car::Anova(gr2, type = 3))
LR Chisq Df Pr(>Chisq)
Shrub_Cov 31.405592 1 0.0000000
Slope 81.306158 1 0.0000000
Elev 11.744565 1 0.0006102
NDVI 34.160149 1 0.0000000
Solar 3.838218 1 0.0500968
aspect.cat 48.337414 8 0.0000001
Slope:Elev 81.698813 1 0.0000000
Solar:aspect.cat 48.296152 8 0.0000001

Year

df_list <- split(j2, j2$year)
year <- map2(df_list, 2, rsf)
## [1] 1104
## [1] 576
## [1] 1674
year1 <- year$`2016`
year2 <- year$`2017`
year3 <- year$`2018`


#Year 1
y1 <- glm(pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar*aspect.cat, data = year1, family = binomial(link=logit))
summary(y1)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar * 
##     aspect.cat, family = binomial(link = logit), data = year1)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.41018  -1.13489   0.09786   1.08899   1.97322  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                4.803e+01  2.404e+02   0.200  0.84161    
## Shrub_Cov                  1.225e+01  1.983e+00   6.177 6.53e-10 ***
## Slope                      8.257e+00  1.250e+00   6.605 3.98e-11 ***
## Elev                       4.684e-02  7.095e-03   6.603 4.04e-11 ***
## NDVI                       5.774e+00  1.773e+00   3.257  0.00112 ** 
## Solar                     -2.179e-04  6.340e-04  -0.344  0.73110    
## aspect.catflat            -1.709e+03  1.019e+03  -1.678  0.09335 .  
## aspect.catnorth           -4.627e+02  2.903e+02  -1.594  0.11102    
## aspect.catnortheast       -2.696e+02  2.825e+02  -0.954  0.33989    
## aspect.catnorthwest       -2.659e+02  2.457e+02  -1.082  0.27907    
## aspect.catsouth           -3.133e+02  2.486e+02  -1.260  0.20766    
## aspect.catsoutheast       -3.100e+02  2.624e+02  -1.181  0.23742    
## aspect.catsouthwest       -7.464e+01  2.421e+02  -0.308  0.75786    
## aspect.catwest            -3.494e+01  2.414e+02  -0.145  0.88494    
## Slope:Elev                -1.143e-02  1.720e-03  -6.648 2.96e-11 ***
## Solar:aspect.catflat       4.497e-03  2.679e-03   1.678  0.09331 .  
## Solar:aspect.catnorth      1.217e-03  7.647e-04   1.591  0.11152    
## Solar:aspect.catnortheast  7.114e-04  7.442e-04   0.956  0.33911    
## Solar:aspect.catnorthwest  6.995e-04  6.467e-04   1.082  0.27943    
## Solar:aspect.catsouth      8.191e-04  6.542e-04   1.252  0.21053    
## Solar:aspect.catsoutheast  8.103e-04  6.901e-04   1.174  0.24034    
## Solar:aspect.catsouthwest  1.950e-04  6.372e-04   0.306  0.75951    
## Solar:aspect.catwest       9.121e-05  6.355e-04   0.144  0.88588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3060.9  on 2207  degrees of freedom
## Residual deviance: 2797.7  on 2185  degrees of freedom
## AIC: 2843.7
## 
## Number of Fisher Scoring iterations: 5
knitr::kable(car::Anova(y1, type = 3))
LR Chisq Df Pr(>Chisq)
Shrub_Cov 40.1184629 1 0.0000000
Slope 54.3672680 1 0.0000000
Elev 48.0808420 1 0.0000000
NDVI 10.8478704 1 0.0009891
Solar 0.1181141 1 0.7310890
aspect.cat 44.9359358 8 0.0000004
Slope:Elev 55.5118090 1 0.0000000
Solar:aspect.cat 44.8027609 8 0.0000004
#Year 2
y2 <- glm(pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar*aspect.cat, data = year2, family = binomial(link=logit))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(y2)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar * 
##     aspect.cat, family = binomial(link = logit), data = year2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.96393  -1.14264   0.00065   1.05996   2.18541  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -1.374e+02  2.585e+02  -0.532 0.594877    
## Shrub_Cov                  8.251e+00  3.054e+00   2.702 0.006898 ** 
## Slope                      9.247e+00  1.829e+00   5.056 4.29e-07 ***
## Elev                       3.925e-02  1.064e-02   3.690 0.000224 ***
## NDVI                       2.227e+01  3.514e+00   6.337 2.34e-10 ***
## Solar                      2.714e-04  6.841e-04   0.397 0.691638    
## aspect.catflat            -8.987e+01  5.253e+02  -0.171 0.864147    
## aspect.catnorth           -1.474e+02  3.609e+02  -0.408 0.682999    
## aspect.catnortheast       -1.342e+05  4.173e+06  -0.032 0.974347    
## aspect.catnorthwest       -1.530e+01  2.689e+02  -0.057 0.954620    
## aspect.catsouth           -9.247e+01  2.725e+02  -0.339 0.734349    
## aspect.catsoutheast       -4.452e+01  3.037e+02  -0.147 0.883455    
## aspect.catsouthwest       -1.143e+01  2.617e+02  -0.044 0.965172    
## aspect.catwest             1.703e+02  2.602e+02   0.654 0.512799    
## Slope:Elev                -1.295e-02  2.532e-03  -5.116 3.13e-07 ***
## Solar:aspect.catflat       2.398e-04  1.381e-03   0.174 0.862125    
## Solar:aspect.catnorth      3.892e-04  9.504e-04   0.410 0.682147    
## Solar:aspect.catnortheast  3.530e-01  1.098e+01   0.032 0.974347    
## Solar:aspect.catnorthwest  4.443e-05  7.072e-04   0.063 0.949899    
## Solar:aspect.catsouth      2.444e-04  7.162e-04   0.341 0.732901    
## Solar:aspect.catsoutheast  1.212e-04  7.981e-04   0.152 0.879297    
## Solar:aspect.catsouthwest  3.247e-05  6.880e-04   0.047 0.962358    
## Solar:aspect.catwest      -4.450e-04  6.842e-04  -0.650 0.515449    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1597.0  on 1151  degrees of freedom
## Residual deviance: 1417.5  on 1129  degrees of freedom
## AIC: 1463.5
## 
## Number of Fisher Scoring iterations: 19
knitr::kable(car::Anova(y2, type = 3))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
LR Chisq Df Pr(>Chisq)
Shrub_Cov 7.4662601 1 0.0062866
Slope 31.0780428 1 0.0000000
Elev 14.2729149 1 0.0001581
NDVI 46.9698942 1 0.0000000
Solar 0.1626145 1 0.6867607
aspect.cat 16.9889504 8 0.0302243
Slope:Elev 32.0912954 1 0.0000000
Solar:aspect.cat 16.9679917 8 0.0304441
#Year 3
y3 <- glm(pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar*aspect.cat, data = year3, family = binomial(link=logit))
summary(y3)
## 
## Call:
## glm(formula = pres ~ Shrub_Cov + Slope * Elev + NDVI + Solar * 
##     aspect.cat, family = binomial(link = logit), data = year3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1833  -1.1177   0.2151   1.0104   3.5098  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.992e+02  1.663e+02   1.197  0.23112    
## Shrub_Cov                  4.875e+00  1.857e+00   2.625  0.00866 ** 
## Slope                      6.699e+00  1.193e+00   5.617 1.94e-08 ***
## Elev                       1.598e-03  7.003e-03   0.228  0.81955    
## NDVI                       1.577e+01  2.530e+00   6.231 4.63e-10 ***
## Solar                     -5.365e-04  4.400e-04  -1.219  0.22273    
## aspect.catflat            -6.754e+02  3.331e+02  -2.028  0.04261 *  
## aspect.catnorth           -3.750e+02  1.938e+02  -1.935  0.05304 .  
## aspect.catnortheast       -4.621e+02  2.450e+02  -1.886  0.05928 .  
## aspect.catnorthwest       -1.167e+02  1.691e+02  -0.690  0.49020    
## aspect.catsouth           -2.557e+00  1.790e+02  -0.014  0.98860    
## aspect.catsoutheast       -1.867e+02  1.938e+02  -0.963  0.33532    
## aspect.catsouthwest       -4.148e+01  1.703e+02  -0.244  0.80757    
## aspect.catwest            -1.371e+02  1.670e+02  -0.821  0.41189    
## Slope:Elev                -9.437e-03  1.675e-03  -5.633 1.77e-08 ***
## Solar:aspect.catflat       1.778e-03  8.762e-04   2.030  0.04239 *  
## Solar:aspect.catnorth      9.858e-04  5.108e-04   1.930  0.05363 .  
## Solar:aspect.catnortheast  1.216e-03  6.458e-04   1.883  0.05964 .  
## Solar:aspect.catnorthwest  3.055e-04  4.455e-04   0.686  0.49289    
## Solar:aspect.catsouth      6.918e-06  4.713e-04   0.015  0.98829    
## Solar:aspect.catsoutheast  4.922e-04  5.100e-04   0.965  0.33453    
## Solar:aspect.catsouthwest  1.093e-04  4.485e-04   0.244  0.80747    
## Solar:aspect.catwest       3.613e-04  4.400e-04   0.821  0.41162    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4641.3  on 3347  degrees of freedom
## Residual deviance: 4006.8  on 3325  degrees of freedom
## AIC: 4052.8
## 
## Number of Fisher Scoring iterations: 6
knitr::kable(car::Anova(y3, type = 3))
LR Chisq Df Pr(>Chisq)
Shrub_Cov 6.9758424 1 0.0082617
Slope 39.1401559 1 0.0000000
Elev 0.0520799 1 0.8194829
NDVI 42.2642812 1 0.0000000
Solar 1.5542334 1 0.2125116
aspect.cat 21.2574172 8 0.0064942
Slope:Elev 39.4966687 1 0.0000000
Solar:aspect.cat 21.2267260 8 0.0065689

Figures

library(sjPlot)

s1 <- glm(pres ~ Shrub_Cov, family = binomial(link=logit), data = mcp95)
newdata <- data.frame(Shrub_Cov = seq(0,1, by = 0.002))
pred.v <- predict(s1, newdata, type = "response", se.fit = TRUE)
pred.v <- as.data.frame(pred.v)
pred.v <- cbind(pred.v, newdata)
#shrub cover only as a predictor i.e. without the influence of the other predictors
#a <- ggplot(pred.v, aes(Shrub_Cov, fit)) + 
  #geom_line() + theme_Publication() + 
  #geom_ribbon(aes(Shrub_Cov, ymin = pred.v$fit - se.fit, ymax = pred.v$fit + se.fit), alpha = 0.3) + 
  #labs(x = "shrub cover", y = "presence")

#a

Predicted values from the full model (m2) i.e. influence of all covariates included. Marginal effects

#full model plot with covariates
a <- plot_model(m5, type = "pred", terms = c("Shrub_Cov")) + 
  theme_Publication() + 
  labs(x = "shrub cover", y = "presence") + 
  ggtitle("") +
  scale_y_continuous(limits = c(0, 1), breaks = c(0,0.25, 0.5, 0.75, 1), label = c("0.0", "0.25", "0.50", "0.75", "1.0")) +
  scale_x_continuous(breaks = c(0, 0.1, 0.2, 0.3), label = c("0", "10", "20", "30"))
## Warning: package 'ggthemes' was built under R version 3.4.4
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
a

b <- plot_model(m5, type = "pred", terms = c("Slope", "Elev")) + 
  theme_Publication()+ 
  labs(x = "slope", y = "presence", color = "elevation") + 
  ggtitle("") + 
  scale_y_continuous(breaks = c(0,0.25, 0.5, 0.75, 1), label = c("0.0", "0.25", "0.50", "0.75", "1.0")) + scale_color_discrete(label = c("low", "medium", "high"))
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
b