Data Processing
In this analysis data is present for 2016/27, 2019 (some sites pre-fire), 2020 (some sites post-fire). The following sections contain the code that was used to process and wrangle the data. Ultimately the data has been written to the Wildlife Ecology’s PostgreSQL database (ari-dev-weda-psql-01) under the southernark schema.
Load and Process Camera Data
Show code
#### 2016 Data: Sess 1 and 2 #### (Already modified by Paul)
# Load mammal camtrap data
CWRMammals <- readRDS('data/CWRMammals.RDS') %>%
mutate(Species = toupper(Species))
# Load in Feral Data
Ferals <- readRDS('data/Ferals.RDS') %>%
mutate(Species = toupper(Species))
# Get site covariates
siteCovs <- CWRMammals %>%
dplyr::select(Station, Group, FoxDen, Geo, Block) %>%
unique() %>%
arrange(Station)
#### Other Data ####
## Add the following element specific variables
# Species
## Add the following time based columns
# Season
# Visit
# Iteration
# Round
## Can link to the following site variables
# Group
# FoxDen
# Geo
# Block
#### 2016 Data (pre-fire) ####
data_2016 <- lapply(list.files(paste0(data_dir,"/2016-Sess-1"), full.names = TRUE), read.csv)
names(data_2016) <- substring(list.files(paste0(data_dir,"/2016-Sess-1")), 1, 3) %>% toupper()
for(x in names(data_2016)) {
colnames(data_2016[[x]]) <- gsub("o", "D", colnames(data_2016[[x]]))
data_2016[[x]]$Species <- x
# data_2016[[x]]$Season <- 1L
# data_2016[[x]]$Visit <- 1L # Only Iteration A before fires
# data_2016[[x]]$Iteration <- "A"
# data_2016[[x]]$Round <- "4.1"
data_2016[[x]] <- dplyr::rename(data_2016[[x]], Station = X) %>%
mutate(across(starts_with("D"), as.integer))
data_2016[[x]] <- tidyr::separate(data_2016[[x]], col = "Station", into = c("Station", "Sess", "Season", "Visit"))
data_2016[[x]]$Iteration <- LETTERS[as.integer(data_2016[[x]]$Visit)]
data_2016[[x]]$Round <- paste(data_2016[[x]]$Season, data_2016[[x]]$Visit, sep = ".")
data_2016[[x]]$Season <- as.integer(data_2016[[x]]$Season)
data_2016[[x]]$Visit <- as.integer(data_2016[[x]]$Visit)
colnames_data <- data_2016[[x]] %>%
dplyr::select(starts_with("D")) %>%
colnames()
data_2016[[x]][,colnames_data] <- as.data.frame(t(apply(data_2016[[x]][,colnames_data],
1,
function(x) {
c(x[!is.na(x)], x[is.na(x)])
})))
data_2016[[x]] <- data_2016[[x]][,colSums(is.na(data_2016[[x]]))<nrow(data_2016[[x]])] %>%
dplyr::select(-Sess)
}
#### 2016 Data (pre-fire) ####
data_2017 <- lapply(list.files(paste0(data_dir,"/2016-Sess-2"), full.names = TRUE), read.csv)
names(data_2017) <- substring(list.files(paste0(data_dir,"/2016-Sess-2")), 1, 3) %>% toupper()
for(x in names(data_2017)) {
colnames(data_2017[[x]]) <- gsub("o", "D", colnames(data_2017[[x]]))
data_2017[[x]]$Species <- x
# data_2017[[x]]$Season <- 1L
# data_2017[[x]]$Visit <- 1L # Only Iteration A before fires
# data_2017[[x]]$Iteration <- "A"
# data_2017[[x]]$Round <- "4.1"
data_2017[[x]] <- dplyr::rename(data_2017[[x]], Station = X) %>%
mutate(across(starts_with("D"), as.integer))
data_2017[[x]] <- tidyr::separate(data_2017[[x]], col = "Station", into = c("Station", "Sess", "Season", "Visit"))
data_2017[[x]]$Iteration <- LETTERS[as.integer(data_2017[[x]]$Visit)]
data_2017[[x]]$Round <- paste(data_2017[[x]]$Season, data_2017[[x]]$Visit, sep = ".")
data_2017[[x]]$Season <- as.integer(data_2017[[x]]$Season)
data_2017[[x]]$Visit <- as.integer(data_2017[[x]]$Visit)
colnames_data <- data_2017[[x]] %>%
dplyr::select(starts_with("D")) %>%
colnames()
data_2017[[x]][,colnames_data] <- as.data.frame(t(apply(data_2017[[x]][,colnames_data],
1,
function(x) {
c(x[!is.na(x)], x[is.na(x)])
})))
data_2017[[x]] <- data_2017[[x]][,colSums(is.na(data_2017[[x]]))<nrow(data_2017[[x]])] %>%
dplyr::select(-Sess)
}
#### 2019 Data (pre-fire) ####
data_2019 <- lapply(list.files(paste0(data_dir,"/2019_Green_A_r1_detections"), full.names = TRUE), read.csv)
names(data_2019) <- substring(list.files(paste0(data_dir,"/2019_Green_A_r1_detections")), 1, 3) %>% toupper()
for(x in names(data_2019)) {
colnames(data_2019[[x]]) <- gsub("o", "D", colnames(data_2019[[x]]))
data_2019[[x]]$Species <- x
data_2019[[x]]$Season <- 4L
data_2019[[x]]$Visit <- 1L # Only Iteration A before fires
data_2019[[x]]$Iteration <- "A"
data_2019[[x]]$Round <- "4.1"
data_2019[[x]] <- rename(data_2019[[x]], Station = X) %>%
mutate(across(starts_with("D"), as.integer))
colnames_data <- data_2019[[x]] %>%
dplyr::select(starts_with("D")) %>%
colnames()
data_2019[[x]][,colnames_data] <- as.data.frame(t(apply(data_2019[[x]][,colnames_data],
1,
function(x) {
c(x[!is.na(x)], x[is.na(x)])
})))
data_2019[[x]] <- data_2019[[x]][,colSums(is.na(data_2019[[x]]))<nrow(data_2019[[x]])]
}
#### 2020 Data (post-fire) ####
data_2020 <- list()
for(i in 1:3) {
data_2020[[i]] <- lapply(list.files(paste0(data_dir,"/2020_", i), full.names = TRUE), read.csv)
names(data_2020[[i]]) <- substring(list.files(paste0(data_dir,"/2020_", i)), 1, 3) %>% toupper()
for(x in names(data_2020[[i]])) {
colnames(data_2020[[i]][[x]]) <- gsub("o", "D", colnames(data_2020[[i]][[x]]))
data_2020[[i]][[x]]$Species <- x
data_2020[[i]][[x]]$Season <- 5L
data_2020[[i]][[x]]$Visit <- i
data_2020[[i]][[x]]$Iteration <- LETTERS[i]
data_2020[[i]][[x]]$Round <- paste0("5.", i)
data_2020[[i]][[x]] <- rename(data_2020[[i]][[x]], Station = 1) %>%
mutate(across(starts_with("D"), as.integer))
colnames_data <- data_2020[[i]][[x]] %>%
dplyr::select(starts_with("D")) %>%
colnames()
data_2020[[i]][[x]][,colnames_data] <- as.data.frame(t(apply(data_2020[[i]][[x]][,colnames_data],
1,
function(x) {
c(x[!is.na(x)], x[is.na(x)])
})))
data_2020[[i]][[x]] <- data_2020[[i]][[x]][,colSums(is.na(data_2020[[i]][[x]]))<nrow(data_2020[[i]][[x]])]
}
}
## 2020 data has empty species records depending on the time. Fill these with zeros
all_species <- sapply(data_2020, names) %>% unlist() %>% unique()
for(i in 1:3) {
sp_r <- names(data_2020[[i]])
sp_np <- setdiff(all_species, sp_r)
for(x in sp_np) {
data_2020[[i]][[x]] <- data_2020[[i]][[1]] %>%
mutate(Species = x) %>%
mutate(across(starts_with("D"), ~ replace(., . == 1L, 0L)))
}
}
# In 2016 and 2017 a predator bait was used for iteration D.
bound_2016_17_noD <- bind_rows(bind_rows(data_2016), bind_rows(data_2017)) %>%
filter(Iteration != "D")
bound_2016_17_D <- bind_rows(bind_rows(data_2016), bind_rows(data_2017)) %>%
filter(Iteration == "D") %>%
mutate(Visit = 3L) %>%
dplyr::select(-Round, - Iteration) %>%
rename_at(paste0("D", 1:52), ~paste0("D", 53:104))
bound_2016_2017 <- left_join(bound_2016_17_noD, bound_2016_17_D)
data_16_17_19_20 <- bind_rows(bound_2016_2017,
bind_rows(data_2019),
bind_rows(lapply(data_2020, bind_rows)))
# Write data to DB
colnames_to_write <- c("Species", "Station", "Season", "Round", "Iteration", "Visit")
presence_data_2016_2020 <- data_16_17_19_20 %>%
dplyr::select(!!colnames_to_write, everything()) %>%
arrange(Species, Season, Visit, Station) %>%
mutate(rowID = 1:n())
presence_data_2016_2020 <- presence_data_2016_2020 %>%
mutate(Season = case_when(Season == 1 ~ 2,
TRUE ~ as.numeric(Season)))
#### Read and write tables ####
DBI::dbWriteTable(con, Id(schema = "southernark",
table = "presence_data_2016_2020"),
presence_data_2016_2020,
overwrite = TRUE)
Process Detection Method Data
Two variables are calculated to account for imperfect detection: Days since deployment and type of bait. Days will be from 1 to 52 and bait type will be either nonPred or Pred, indicating the type of lure used with the camera trap.
Show code
days_data_2016_2020 <- presence_data_2016_2020 %>%
dplyr::select(starts_with("D"))
for(i in 1:52) {
i2 <- i + 52
non_NAs_first <- which(!is.na(days_data_2016_2020[[i]]))
non_NAs_second <- which(!is.na(days_data_2016_2020[[i2]]))
days_data_2016_2020[non_NAs_first,i] <- i
days_data_2016_2020[non_NAs_second,i2] <- i
}
bait_data_2016_2020 <- presence_data_2016_2020 %>%
dplyr::select(starts_with("D")) %>%
mutate_all(as.character)
for(i in 1:52) {
i2 <- i + 52
non_NAs_first <- which(!is.na(bait_data_2016_2020[[i]]))
non_NAs_second <- which(!is.na(bait_data_2016_2020[[i2]]))
bait_data_2016_2020[non_NAs_first,i] <- "nonPred"
bait_data_2016_2020[non_NAs_second,i2] <- "Pred"
}
#### Read and write tables ####
DBI::dbWriteTable(con, Id(schema = "southernark",
table = "days_data_2016_2020"),
days_data_2016_2020,
overwrite = TRUE)
DBI::dbWriteTable(con, Id(schema = "southernark",
table = "bait_data_2016_2020"),
bait_data_2016_2020,
overwrite = TRUE)
Show code
read_excel_allsheets <- function(filename, tibble = FALSE) {
# I prefer straight data.frames
# but if you like tidyverse tibbles (the default with read_excel)
# then just pass tibble = TRUE
sheets <- readxl::excel_sheets(filename)
x <- lapply(sheets, function(X) readxl::read_excel(filename, sheet = X))
if(!tibble) x <- lapply(x, as.data.frame)
names(x) <- sheets
x
}
# For each camera deployment get the day and the
dates <- read_excel_allsheets(paste0(data_dir, "/Camera_Dates.xlsx"))
dates[[1]] <- dates[[1]] %>%
rename(Station = SiteID, Setup_date = DateDeploy, Retrieval_date = DateRetrie) %>%
tidyr::separate(col = "round", into = c("Season", "Visit")) %>%
mutate(across(c(Season, Visit), as.integer))
dates[[3]] <- dates[[3]] %>%
tidyr::separate(col = "round", into = c("Season", "Visit")) %>%
mutate(across(c(Season, Visit), as.integer),
Season = 5L)
dates[[2]] <- dates[[2]] %>%
mutate(Season = 4L,
Visit = 1L)
dates_joined <- bind_rows(dates)
# convert to sf
dates_sf <- st_as_sf(dates_joined, coords = c("Easting", "Northing")) %>%
sf::`st_crs<-`(28355) %>%
sf::st_transform(4283) %>%
mutate(Station_Season = paste(Station, Season, sep = "_"))
dates_long <- dates_sf %>%
rowwise() %>%
do(data.frame(Station=.$Station,
Season = .$Season,
Visit = .$Visit,
day=seq(.$Setup_date,.$Retrieval_date,by="day"),
day_deploy = (seq(.$Setup_date,.$Retrieval_date,by="day")-.$Setup_date)/86400)) %>%
ungroup() %>%
mutate(day_deploy = as.integer(case_when(Visit == 4 ~ day_deploy + 52L,
TRUE ~ day_deploy)))
dates_long_sf <- dates_long %>%
dplyr::ungroup() %>%
left_join(southernark_site_data %>% dplyr::select(Station) %>% distinct()) %>%
mutate(Year = lubridate::year(day),
Day = lubridate::yday(day)) %>%
st_as_sf() %>%
filter(!st_is_empty(geometry))
# load in the grided daily data
daily_data <- list()
nc_files <- list.files(paste0(data_dir, "/daily_temps"))
years_nc <- gsub(x = nc_files, pattern = "tmax.|.nc", replacement = "")
daily_data <- lapply(paste(data_dir, "daily_temps", nc_files, sep = "/"), raster::brick)
names(daily_data) <- years_nc
temp_max <- numeric()
dates_long_no_sf <-sf::st_drop_geometry(dates_long_sf)
for(i in 1:nrow(dates_long_sf)) {
raster_year <- as.character(dates_long_no_sf[i, "Year"])
raster_day <- as.integer(dates_long_no_sf[i, "Day"])
temp_max[i] <- raster::extract(x = daily_data[[raster_year]][[raster_day]], dates_long_sf[i,], fun = mean, na.rm =TRUE, method = "bilinear")
}
dates_long_sf$tmax <- temp_max
max_narm <- function(x) {
max(x, na.rm = TRUE)
}
dates_wide <- dcast(as.data.table(dates_long_sf), formula = Station+Season ~ day_deploy, value.var = "tmax", fun.aggregate = max_narm)
invisible(lapply(names(dates_wide),function(.name) set(dates_wide, which(is.infinite(dates_wide[[.name]])), j = .name,value =NA)))
temp_wide <- dates_wide
colnames(temp_wide)[3:ncol(temp_wide)] <- paste0("D", 1+as.integer(colnames(temp_wide)[3:ncol(temp_wide)]))
temp_wide$D104 <- NA
temp_wide_ordered <- presence_data_2016_2020 %>%
dplyr::select(Station, Season, rowID) %>%
left_join(temp_wide %>%
mutate(Season = case_when(Season == 1 ~ 2L,
TRUE ~ Season)),
by = c("Station", "Season")) %>%
arrange(rowID) %>%
dplyr::select(starts_with("D"))
#### Read and write tables ####
DBI::dbWriteTable(con, Id(schema = "southernark",
table = "temp_days_data_2016_2020"),
temp_wide_ordered,
overwrite = TRUE)
Load and Process Spatial Data
We include site covariates on a seasonal basis for:
- Landscape fire severity (Department of Environment 2021)
- Bioclim variables (Karger et al. 2017):
- BIO1: Temperature (Mean)
- BIO12: Precipitation (Annual)
- BIO15: Precipitation (Seasonality)
- BIO1: Temperature (Mean)
Site-specific covariates extracted various environmental variables from a range of open-access data sources. For the landscape variables of fire severity we use a distance and severity weighted calculation from Lindenmayer et al. (2021). This method allows for the effects of disturbance at the camera trap location and surrounding area (250m radius) to be included as a continuous variable. Because the fire severity data is composed of 5 severity classes related to metrics like the percentage of crown scorch (Department of Environment 2021) the fire severity at a site will be weighted accordingly. The spatio-severity weighting of fire can be visualised through the a simulated contour plot, where the higher weighted parameter space is coloured with darker reds (Figure 1). The calculation of the fire severity score weighting for each raster pixel within a 250m radius of the site follows the following gaussian kernel estimate:
\[fire\ severity = \sum\limits_{i=1}^n e^{-ϕ_{1} d_{n}^2} \times (1 - e^{-ϕ_{2} s_{n}^2})\] Where:
- \(phi_1\) is the distance scale parameter (73)
- \(d_{n}\) is the distance of the raster cell n from the site (in km)
- \(phi_2\) is the fire severity scale parameter (0.0865)
- \(s_{n}\) is the fire severity of the raster cell n
Bioclim variables were obtained from CHELSA (Climatologies at high resolution for the earth’s land surface areas) (Karger et al. 2017). A temperature-based (BIO01) variable and two precipitation-based (BIO12 and BIO15) variables were used. These variables represent annual mean temperature (BIO01), annual precipitation (BIO12) and precipitation seasonality (BIO15). Bioclimatic variables were extracted at each station point location.
Show code
# Logging example
# distance <- seq(from = 0, to = 0.25, length.out = 100)
# time <- seq(from = 0, to = 30)
#
# wl <- matrix(dimnames = list(distance, time), nrow = length(distance), ncol = length(time))
#
# for(d in 1:length(distance)) {
# for(t in 1:length(time)) {
# wl[d,t] <- wk(dist = distance[d], lag = time[t], dp=2, tp = 2, ssp = 73, tsp = 0.0050)
# }
# }
#############
# Fire example
# distance <- seq(from = 0, to = 0.25, length.out = 100)
sev <- seq(from = 0, to = 5, length.out = 30)
distance <- seq(from = 0, to = 0.25, length.out = 100)
w <- matrix(dimnames = list(distance, sev), nrow = length(distance), ncol = length(sev))
for(d in 1:length(distance)) {
for(t in 1:length(sev)) {
w[d,t] <- wk(dist = distance[d], lag = sev[t], dp=2, tp = 2, ssp = 73, tsp = 0.0865, inverse_lag = TRUE)
}
}
filled.contour(x = distance, y = sev, z = w,
plot.title = title(main = "Spatio-severity weighting for fire impact",
xlab = "Distance from site (km)",
ylab = "Fire severity"),
key.title = title(main = "Weight"))
Figure 1: Weightings of fire severity at a particular station/site in 2020 were calculated based on the fire severity measured at the site and in the surrounding area (250m radius)
Show code
# filled.contour(x = distance, y = time, z = wl,
# plot.title = title(main = "Spatio-temporal weighting for logging impact",
# xlab = "Distance from site (km)",
# ylab = "Years since harvest"),
# key.title = title(main = "Weight"))
Show code
camera_locations <- read.delim('data/Cam_Locations.txt') %>%
dplyr::select(Station = FULLID, X, Y) %>%
sf::st_as_sf(., coords = c("X", "Y"), crs = 28355) %>% #utm 55h
sf::st_transform(4283)
# Pull down entire sf layer for last log in east gippsland
LASTLOG25_EAST_GIPPSLAND <- st_read(dsn = con,
layer = DBI::Id(schema = "external", table = "LASTLOG25_EAST_GIPPSLAND"))
ex <- extent(sf::st_buffer(camera_locations %>%
sf::st_transform(3111), dist = 1000))
#### Elevation #####
elev <- raster::raster(here::here("_posts",
"2021-09-21-southernarkfire",
"data",
"vmelev_dem10m.tif")) #Dem10
# tri <- tri(elev, exact = TRUE)
#
# raster::writeRaster(tri, filename = here::here("_posts",
# "2021-09-21-southernarkfire",
# "data",
# "vmelev_tri.tif"))
tri <- raster::raster(here::here("_posts",
"2021-09-21-southernarkfire",
"data",
"vmelev_tri.tif"))
FIRE_SEVERITY_2020 <- rpostgis::pgGetRast(con,
name = c('external','FIRE_SEVERITY_2020'),
boundary = c(ex[4], ex[3], ex[2], ex[1]))
FIRE_SEVERITY_RESCALED <- FIRE_SEVERITY_2020
FIRE_SEVERITY_RESCALED <- FIRE_SEVERITY_RESCALED - 1
FIRE_SEVERITY_RESCALED[FIRE_SEVERITY_RESCALED <= 0] <- NA
# FIRE_SEVERITY_RESCALED <- 4 - FIRE_SEVERITY_RESCALED
gippy_rast <- raster(gippy %>% sf::st_transform(3111), res = 20)
# sample 1km buffers
gippy_points <- camera_locations %>%
sf::st_transform(3111) %>%
sf::st_buffer(dist = 1000) %>%
sf::st_transform(4283)
## Road Connectivity
# Extract roads for the points
spatial_data <- camera_locations %>%
left_join(data.frame(Station = rep(camera_locations$Station,4),
Season = rep(c(1,2,4,5), each = nrow(camera_locations))), by = "Station")
for(i in 1:nrow(gippy_points)) {
gippy_points[i, "road_density_1000"] <- get_road_density(gippy_points[i,])
}
for(i in which(is.na(gippy_points$road_density_1000))) {
gippy_points[i, "road_density_1000"] <- get_road_density(gippy_points[i,])
}
spatial_data <- left_join(spatial_data, gippy_points %>% st_drop_geometry(), by = "Station")
## Generate Disturbance Scores
# For gippy points get areas that have been logged in the last 30 years (1km buffer)
# For each point get the weighting of the amount of disturbance at the site
spatial_data_list <- split(spatial_data, spatial_data$Season)
spatial_data_list <- mclapply(spatial_data_list, function(x) {
current_year <- unique(x$Season) + 2015L
# Last log raster
last_log_raster <- LASTLOG25_EAST_GIPPSLAND %>%
st_transform(3111) %>%
filter(year(ENDDATE) > 1990 & year(ENDDATE) < current_year) %>%
mutate(years_since = current_year - year(ENDDATE)) %>%
fasterize::fasterize(raster = gippy_rast, background = NA, field = "years_since", fun = "max")
x$logging_density <- landscape_disturbance(sites = x,
disturbance_raster = last_log_raster,
buffer = 250, # In m
dp=2, tp = 2, ssp = 73, tsp = 0.0050)
return(x)
}, mc.cores = 4)
spatial_data_list[["5"]]$fire_severity <- landscape_disturbance(sites = spatial_data_list[["5"]],
disturbance_raster = FIRE_SEVERITY_RESCALED,
buffer = 250, # In m 5 ha area
dp=2,
tp = 2,
ssp = 73, # 250m = 0.1 weighting
tsp = 0.0865, #10 fold between highest severity and lowest severity burn (5 to 1)
inverse_lag = TRUE)# change for fire severity: inverse because higher values are worse:
for (i in 1:3) {
spatial_data_list[[i]]$fire_severity <- as.numeric(0)
}
spatial_data <- bind_rows(spatial_data_list)
#Join this data to the other site covariates
site_data <- left_join(spatial_data %>% sf::st_transform(4283), siteCovs, by = "Station")
site_data$scaled_road_density <- scales::rescale(site_data$road_density_1000)
site_data$scaled_logging <- scales::rescale(site_data$logging_density)
site_data$scaled_fire_severity <- scales::rescale(site_data$fire_severity)
# Read in rasters for bioclim
bioclims <- list.files("data/bioclim", full.names = FALSE)
bio_name <- gsub(pattern = ".tif", replacement = "", x = bioclims)
bioclim_site_data <- lapply(bioclims, function(x) {
sf_df <- site_data %>%
dplyr::select(Station)
bio_raster <- raster::raster(paste0("data/bioclim/", x))
raster::extract(bio_raster, y = sf_df)
}) %>% as.data.frame() %>% `colnames<-`(bio_name)
site_data <- bind_cols(site_data, bioclim_site_data)
southernark_site_data <- presence_data_2016_2020 %>%
dplyr::select(Station, Season, Species, rowID, Visit) %>%
mutate(Group = substr(Station, 1,4)) %>%
left_join(site_data %>%
dplyr::select(Station, Season, Geo, Block, road_density_1000, scaled_road_density, logging_density, scaled_logging, fire_severity, scaled_fire_severity, BIO01, BIO04, BIO12, BIO15), by = c("Station", "Season")) %>%
st_as_sf()
# Add in fox density
fox_den_files <- list.files('data/fox_density', full.names = TRUE)
fox_den <- raster::stack(fox_den_files)
names(fox_den) <- paste0("FoxDensity_", 2015:2019)
crs(fox_den) <- 3111
# extract based on sites
fox_den_sites <- bind_cols(st_drop_geometry(southernark_site_data[which(sf::st_is_empty(southernark_site_data) == FALSE), c("Station","Season") ]), as.data.frame(extract(fox_den, y = southernark_site_data[which(sf::st_is_empty(southernark_site_data) == FALSE),])))
fox_den_year_sites <- melt(as.data.table(fox_den_sites), id.vars = c("Station", "Season"), value.name = "FoxDensity") %>%
mutate(year_fd = as.integer(substr(variable, start = 12, 15))) %>%
filter(Season + 2015L == year_fd | (Season == 5 & year_fd == 2019)) %>%
dplyr::select(-year_fd, -variable) %>%
unique()
southernark_site_data <- left_join(southernark_site_data, fox_den_year_sites, by = c("Station", "Season")) %>%
arrange(rowID)
# Add elevation
elev_df <- southernark_site_data %>%
filter(!sf::st_is_empty(geometry)) %>%
dplyr::select(Station) %>%
distinct()
elev_df$elevation <- raster::extract(elev, elev_df)
elev_df$tri <- raster::extract(tri, elev_df)
elev_df <- elev_df %>% st_drop_geometry()
southernark_site_data <- left_join(southernark_site_data, elev_df, by = c("Station")) %>%
arrange(rowID)
# Correct Seasons to 16/17 (season 2)
southernark_site_data <- southernark_site_data %>%
mutate(Season = case_when(Season == 1 ~ 2,
TRUE ~ Season))
lat_lon <- as.data.frame(st_coordinates(southernark_site_data)) %>%
rename(lat = Y, lon = X)
southernark_site_data <- bind_cols(southernark_site_data, lat_lon)
# Write sf spatial data to db
st_write(obj = southernark_site_data,
dsn = con,
layer = Id(schema = "southernark", table = "southernark_site_data"),
delete_layer = TRUE)
Data Summary
Total Detections for Each Species, Each Season
A total of 10 species were monitored on the camera traps in all four years (2016, 2017, 2019 and 2020). Three of them being introduced predators (cats, dogs and foxes) while seven were native species (common brushtail possum, common ringtail possum, goanna, long-footed potoroo, long-nosed bandicoot, long-nosed potoroo and southern brown bandicoot). In 2019 27 sites recorded an unidentifiable potoroo spp (Table 1).
Show code
summarised_detections_all <- presence_data_2016_2020 %>%
as.data.table() %>%
melt(id.vars = c("Species", "Season"), measure.vars = paste0("D", 1:104), value.name = "detection") %>%
mutate(Fire = case_when(Season < 5 ~ "Pre-fire",
TRUE ~ "Post-fire")) %>%
group_by(Species, Season) %>%
summarise(across(starts_with("D"))) %>%
summarise(`0` = sum(detection == 0, na.rm = TRUE),
`1` = sum(detection == 1, na.rm = TRUE),
`NA` = sum(is.na(detection))) %>%
ungroup() %>%
mutate(Detection_Rate = `1`/`0`)
summarised_detections_site_season <- presence_data_2016_2020 %>%
as.data.table() %>%
melt(id.vars = c("Species", "Season", "Station"), measure.vars = paste0("D", 1:104), value.name = "detection") %>%
mutate(Fire = case_when(Season < 5 ~ "Pre-fire",
TRUE ~ "Post-fire")) %>%
group_by(Species, Season, Station) %>%
# summarise(across(starts_with("D"))) %>%
summarise(detections = sum(detection == 1, na.rm = TRUE),
non_detections = sum(detection == 0, na.rm = TRUE),
detection_rate = sum(detection == 1, na.rm = TRUE)/(sum(detection == 1, na.rm = TRUE) + sum(detection == 0, na.rm = TRUE)),
Presence = max(detection, na.rm = TRUE)) %>%
ungroup()
summarised_season_aggregates <- summarised_detections_site_season %>%
group_by(Species, Season) %>%
summarise(Stations = n(),
`Stations Present` = sum(Presence, na.rm = TRUE),
`Presence Percentage` = paste0(round(100*`Stations Present`/Stations, 2), "%"))
summarised_season_aggregates %>%
kbl(format = "html", caption = "Total for the target species over the study period") %>%
kable_styling(full_width = FALSE) %>%
column_spec(1:3, bold = TRUE) %>%
collapse_rows(columns = 1:2)
| Species | Season | Stations | Stations Present | Presence Percentage |
|---|---|---|---|---|
| CAT | 2 | 685 | 200 | 29.2% |
| 4 | 238 | 36 | 15.13% | |
| 5 | 250 | 47 | 18.8% | |
| CBP | 2 | 685 | 145 | 21.17% |
| 4 | 238 | 76 | 31.93% | |
| 5 | 250 | 68 | 27.2% | |
| CRP | 2 | 685 | 21 | 3.07% |
| 4 | 238 | 6 | 2.52% | |
| 5 | 250 | 4 | 1.6% | |
| DOG | 2 | 685 | 27 | 3.94% |
| 4 | 238 | 14 | 5.88% | |
| 5 | 250 | 4 | 1.6% | |
| FOX | 2 | 685 | 21 | 3.07% |
| 4 | 238 | 10 | 4.2% | |
| 5 | 250 | 13 | 5.2% | |
| GOA | 2 | 685 | 168 | 24.53% |
| 4 | 238 | 60 | 25.21% | |
| 5 | 250 | 2 | 0.8% | |
| LFP | 2 | 685 | 208 | 30.36% |
| 4 | 238 | 86 | 36.13% | |
| 5 | 250 | 136 | 54.4% | |
| LNB | 2 | 685 | 154 | 22.48% |
| 4 | 238 | 21 | 8.82% | |
| 5 | 250 | 21 | 8.4% | |
| LNP | 2 | 685 | 24 | 3.5% |
| 4 | 238 | 3 | 1.26% | |
| 5 | 250 | 2 | 0.8% | |
| POT | 4 | 238 | 27 | 11.34% |
| SBB | 2 | 685 | 6 | 0.88% |
| 4 | 238 | 2 | 0.84% | |
| 5 | 250 | 1 | 0.4% |
Effect of Fire Severity in 2020
Fire severity varied throughout the landscape and across sites. Comparing occupied vs non-occupied sites in season 5 (2020), we see that there is some reason to believe that fire severity negatively impacted occupancy; with unoccupied sites generally having greater landscape fire severity scores (Figure 2). However this difference is not clear-cut and appears to vary between species. Hence a mixed effect occupancy model can be used.
Show code
fire_presence_2020 <- left_join(summarised_detections_site_season, southernark_site_data)
fire_pres_plot <- fire_presence_2020 %>%
filter(Season == 5) %>%
mutate(`Site Presence` = as.character(Presence)) %>%
ggplot(aes(x = scaled_fire_severity, y = Species)) +
# geom_col(aes(fill = `Site Presence`)) +
ggridges::geom_density_ridges(aes(fill = `Site Presence`), alpha = 0.7, panel_scaling = FALSE) +
scale_fill_brewer(palette = "Set2") +
xlab("Scaled fire severity") +
theme_minimal()
fire_pres_plot
Figure 2: Comparison of fire severity scores in 2020 for stations that recorded either presences or absences.
Visualising Site Data
The sites shown here are those that have pre-fire and post-fire surveys. The raster overlay depicts the severity and extent of the 2019/2020 bushfires. The majority of sites in the analysis were impacted by the fires, however the severity varied.
Show code
season_5_map_data <- southernark_site_data %>%
filter(Season == 5) %>%
dplyr::select(Station, scaled_fire_severity) %>%
distinct()
# Fire data
load('data/fire_lite.rda')
gen_pal <- colorNumeric("inferno", domain = c(0,6), na.color = "transparent", reverse = TRUE)
base_map <- leaflet() %>%
addProviderTiles(providers$CartoDB.Positron, group = "Map") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satellite") %>%
# addWMSTiles(baseUrl = "http://data.gov.au/geoserver/fy-bushfire-boundaries-201920/ows?SERVICE=WMS",
# layers = "ckan_f91c2ffd_a0a0_4042_a8ac_71757a6ef727",
# group = "wms_crown_land",
# WMSTileOptions(format = "image/png", transparent = TRUE),
# attribution = "© DELWP") %>%
leaflet.extras::addFullscreenControl() %>%
leafem::addMouseCoordinates()
# Create a continuous palette function
firepal <- colorNumeric(
palette = "inferno",
domain = c(0,1), reverse = TRUE)
map <- base_map %>%
addCircles(data = season_5_map_data,
radius = 500,
label = ~ paste(Station, round(scaled_fire_severity, 2), sep = ": "),
fillColor = ~firepal(scaled_fire_severity), fill = TRUE,
stroke = FALSE, fillOpacity = 0.8,
group = "Stations") %>%
addRasterImage(fire_lite, colors = gen_pal, opacity = 0.75, group = "Fire Severity", project = FALSE) %>%
addLegend(pal = gen_pal, values = c(2,6),
title = "Fire Severity") %>%
addLayersControl(overlayGroups = c("Stations", "Fire Severity"), baseGroups = c("Map", "Satellite"))
map
Potoroo Problem
In 2019 27 stations/sites detected a potoroo species (long-nosed potoroo or long-footed potoroo. Because those records could not be identified to a species level the presence of LNP and LFP during that season may be biased in that records are underestimated. In order to correct this the ‘potoroo’ records are reclassed as either LNP or LFP based on i) previous site detections of either species or if no records were recorded in 206/2017 then (ii) mean detetection rates of LFP vs LNP.
Show code
# Potoroos had issues of being identified in 2019 with detections being classed as "potoroos spp"
# If we exclude these records we risk estimates for occupancy in 2019 (pre-fire) being underestimated... Potentially makeing the effect of fire appear less negative or even positive.
# Get site specific potoroo detection rates in 2016/2017
lnp_lfp_detection_rates <- summarised_detections_site_season %>%
filter(Species %in% c("LNP", "LFP") & Season < 3) %>%
group_by(Species, Station) %>%
summarise(detection_rate = sum(detections)/(sum(detections) + sum(non_detections))) %>%
as.data.table() %>%
dcast(Station ~ Species, value.var = "detection_rate")
mean_lnp_rate <- mean(lnp_lfp_detection_rates$LNP)
mean_lfp_rate <- mean(lnp_lfp_detection_rates$LFP)
mean_lfp_prob <- mean_lfp_rate/(mean_lfp_rate+mean_lnp_rate)
lnp_lfp_detection_rates <- lnp_lfp_detection_rates %>%
mutate(lfp_prob = case_when(LFP == 0 & LNP == 0 ~ mean_lfp_rate,
TRUE ~ LFP/(LFP+LNP)))
# For the Potoroo detections assign them as LFP or LNPs on the same day
potoroo_data <- presence_data_2016_2020 %>%
filter(Species == "POT") %>%
dplyr::select(Station, starts_with("D")) %>%
left_join(lnp_lfp_detection_rates %>% dplyr::select(Station, lfp_prob)) %>%
mutate(lfp_prob = tidyr::replace_na(lfp_prob, mean_lfp_prob))
potoroo_data_obs <- potoroo_data %>%
as.data.table() %>%
melt(id.vars = c("Station", "lfp_prob")) %>%
filter(value == 1) %>%
dplyr::rowwise() %>%
mutate(Species = case_when(value == 1 ~ ifelse(rbinom(1,1,lfp_prob) == 1, "LFP", "LNP"),
TRUE ~ NA_character_)) %>%
ungroup()
# fill presence data
presence_data_2016_2020_mod <- presence_data_2016_2020
for(i in 1:nrow(potoroo_data_obs)) {
presence_data_2016_2020_mod[presence_data_2016_2020_mod$Season == 4 &
presence_data_2016_2020_mod$Visit == 1 &
presence_data_2016_2020_mod$Station == potoroo_data_obs[i, ]$Station &
presence_data_2016_2020_mod$Species == potoroo_data_obs[i,]$Species,
potoroo_data_obs[i,]$variable] <- 1
}
Modelling Approach
The main goal of this study is to compare the effect of the 2019/2020 bushfires on the occupancy of several native species. The occupancy rate is expected to decline after the fires. Especially in areas severely burnt. The model will also take into account the effect of other variables (e.g. predation, climate, disturbances (roads and logging)). With regards to predation we use a precursor model to generate a series of variables for cat, dog and fox occupancy based on predictions.
Predator Presence Model: GLMM
Our key modelling approach uses a Bayesian approach (via the ubms R package). However, courser maximum likelihood methods can also be used to assess the drivers of species occupancy. The GLMM allows for random effects but does not include a detection submodel, alternatively the unmarked GLM contains detection observation submodels but no random effects. The limitations of these modeling tools from the lme4 and unmarked modelling R packages places the use of ubms or JAGS/BUGS as ideal packages/languages for implementing our analysis. However, for efficient (fast) model selection we used GLMM’s with occupancy data across days condensed into presence (1) or absence (0). Parameter estimates in this approach should be similar albeit more marginal with lower predicted occupancy rates as imperfect detection is not accounted for.
Show code
Predators <- c("CAT", "FOX", "DOG")
# Restrict Data to only the sites that were consistently visited
sites_visited_post <- presence_data_2016_2020_mod %>%
filter(Season == 5) %>%
pull(Station) %>%
unique()
sites_visited_pre_post <- presence_data_2016_2020_mod %>%
filter(Season < 5 & Station %in% sites_visited_post) %>%
pull(Station) %>%
unique()
#Drop rows with no geometry
no_geo <- southernark_site_data %>%
filter(st_is_empty(geometry)) %>%
pull(rowID)
# Restrict the model to data
predator_model_rows <- presence_data_2016_2020_mod %>%
filter(Species %in% Predators & !(rowID %in% no_geo)) %>% #& Station %in% sites_visited_pre_post
pull(rowID)
predator_site_data <- presence_data_2016_2020_mod %>%
filter(rowID %in% predator_model_rows) %>%
dplyr::select(Species, Station, Season, Visit) %>%
distinct() %>%
left_join(sf::st_drop_geometry(southernark_site_data))
data_for_predator_glmm <- predator_site_data %>%
left_join(summarised_detections_site_season, by = c("Station", "Species", "Season"))
umf_predators <- unmarkedFrameOccu(y = presence_data_2016_2020_mod[predator_model_rows,] %>%
dplyr::select(starts_with("D")),
siteCovs = predator_site_data,
obsCovs = list(days = days_data_2016_2020[predator_model_rows,],
baitMethod = bait_data_2016_2020[predator_model_rows,],
temp = temp_days_data_2016_2020[predator_model_rows,]))
if(!file.exists("models/mod_dredge_predator.rds")) {
global_model_predator <- lme4::glmer(Presence~Species
+ Species*scaled_fire_severity # species reactions to fire
+ Species*scale(BIO01) # Climate
# + Species*scale(BIO04) # Climate
+ Species*scale(BIO12) # Climate
+ Species*scale(BIO15) # Climate
+ (1|Season)
+ (1|Group), # Site specific effects\
family = binomial,
data = data_for_predator_glmm,
control=lme4::glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)), na.action = "na.pass")
subsetExpr <- makeRule(scale(BIO01), scale(BIO12), scale(BIO15))
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", 8), type = clusterType))
clusterExport(clust, "data_for_predator_glmm")
clusterExport(clust, c("glmer", "glmerControl"), envir = as.environment("package:lme4"))
mod_dredge_predator <- MuMIn::pdredge(global_model_predator,
fixed = ~ Species + Species*scaled_fire_severity + (1|Season) + (1|Group),
subset = subsetExpr)
saveRDS(mod_dredge_predator, "models/mod_dredge_predator.rds", compress = "xz")
saveRDS(global_model_predator, "models/global_model_predator.rds", compress = "xz")
} else {
mod_dredge_predator <- readRDS("models/mod_dredge_predator.rds")
global_model_predator <- readRDS("models/global_model_predator.rds")
}
mod_dredge_predator %>%
as.data.frame() %>%
dplyr::select(-`(Intercept)`, -`scale(BIO01)`, -`scale(BIO12)`, -`scale(BIO15)`, -`scaled_fire_severity`) %>%
kbl(format = "html", caption = "Model selection for the predator models using GLMM", digits = 3, escape = FALSE) %>%
kable_styling(full_width = FALSE) %>%
row_spec(1, bold = TRUE) %>%
scroll_box(width = "800px", height = "400px")
| Species | scale(BIO01):Species | scale(BIO12):Species | scale(BIO15):Species | scaled_fire_severity:Species | df | logLik | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|
| 63 |
|
|
|
|
|
17 | -997.853 | 2029.881 | 0.000 | 9.776036e-01 |
| 27 |
|
|
|
NA |
|
14 | -1004.987 | 2038.094 | 8.213 | 1.609772e-02 |
| 54 |
|
NA |
|
|
|
14 | -1006.748 | 2041.616 | 11.735 | 2.767199e-03 |
| 9 |
|
|
NA | NA |
|
11 | -1009.944 | 2041.963 | 12.082 | 2.325954e-03 |
| 45 |
|
|
NA |
|
|
14 | -1008.207 | 2044.534 | 14.653 | 6.432954e-04 |
| 18 |
|
NA |
|
NA |
|
11 | -1011.711 | 2045.498 | 15.617 | 3.972549e-04 |
| 0 |
|
NA | NA | NA |
|
8 | -1015.840 | 2047.721 | 17.840 | 1.307142e-04 |
| 36 |
|
NA | NA |
|
|
11 | -1014.161 | 2050.397 | 20.516 | 3.429638e-05 |
Occupancy Predator Model: Stacked Approach
Model
The aim of the predator occupancy model is to generate accurate predictions of the occupancy of the three species of predators at each location in each season. The chosen model for based on the model selection above (Table 2) is:
Detection Submodel: ~ Baiting Method x days Occupancy Submodel: ~ Species + Species x Fire + Species x BIO01 + Species x BIO12 + Species x BIO15 + (1|Group) + (1|Season)
Show code
if(!file.exists("models/predator_global_ubms.rds")) {
predator_global_ubms <- stan_occu(~Species+baitMethod*days # affecting detection
~Species
+ Species*scaled_fire_severity # species reactions to fire
+ Species*scale(BIO01) # Climate
# + Species*scale(BIO04) # Climate
+ Species*scale(BIO12) # Climate
+ Species*scale(BIO15) # Climate
+ (1|Group) # Site specific effects
+ (1|Season), # Seasonal effects
data = umf_predators, chains=4, iter=4000)
# predator_global_unmarked <- occu(~baitMethod+days # affecting detection
# ~Species
# + scaled_fire_severity # species reactions to fire
# + scaled_road_density # predator access on road
# + scaled_logging # habitat modifications
# + BIO01 # Climate
# + BIO04 # Climate
# + BIO12 # Climate
# + BIO15 # Climate
# # + (1|Station) # Site specific effects
# + as.factor(Season), # Seasonal effects
# data = umf_predators)
saveRDS(predator_global_ubms, "models/predator_global_ubms.rds", compress = "xz")
} else {
predator_global_ubms <- readRDS("models/predator_global_ubms.rds")
}
Model Summary
The summary for the top model (global) is broken up into the Occupancy sub-model and the Detection sub-model. The estimates for the model parameters, standard deviation (SD), 95% confidence intervals, effective sample size and Rhat is presented below (Table (tab:modelOutputPred)). n_eff > 200 and Rhat < 1.05 are general indicators the model has reached convergence. The random effect of ‘Group’ has low effective sample size due to the number of levels and the lack of repeat measurements. The effective sample size was even lower if station was used instead of group.
Show code
predator_summary <- summary(predator_global_ubms, submodel = "state")
predator_summary_sig <- which((predator_summary$`2.5%` * predator_summary$`97.5%`) > 0)
predator_summary%>%
dplyr::select(-`25%`, -`50%`, -`75%`) %>%
kbl("html", caption = "Model estimates showing the effects of various covariates on the probability of occupancy for predator species. Estimates with a 95% CI not overlapping zero are shown in bold. The reference level (Intercept) is 'CAT'.", digits = 3) %>%
kable_styling(full_width = FALSE) %>%
row_spec(predator_summary_sig, bold = TRUE)
| mean | se_mean | sd | 2.5% | 97.5% | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| (Intercept) | -0.173 | 0.011 | 0.315 | -0.768 | 0.559 | 886.269 | 1.002 |
| SpeciesDOG | -1.867 | 0.010 | 0.341 | -2.489 | -1.177 | 1107.115 | 1.003 |
| SpeciesFOX | -2.905 | 0.010 | 0.272 | -3.446 | -2.375 | 760.574 | 1.004 |
| scaled_fire_severity | -1.136 | 0.018 | 0.632 | -2.546 | -0.026 | 1270.626 | 1.001 |
| scale(BIO01) | -0.309 | 0.002 | 0.119 | -0.550 | -0.082 | 2816.664 | 1.001 |
| scale(BIO12) | 0.024 | 0.004 | 0.131 | -0.232 | 0.272 | 1290.406 | 1.002 |
| scale(BIO15) | -0.063 | 0.003 | 0.124 | -0.308 | 0.182 | 1614.376 | 1.000 |
| SpeciesDOG:scaled_fire_severity | -0.605 | 0.017 | 0.956 | -2.637 | 1.140 | 3057.998 | 1.001 |
| SpeciesFOX:scaled_fire_severity | 0.983 | 0.017 | 0.762 | -0.509 | 2.432 | 2097.096 | 1.001 |
| SpeciesDOG:scale(BIO01) | 0.153 | 0.005 | 0.220 | -0.274 | 0.587 | 2055.440 | 1.003 |
| SpeciesFOX:scale(BIO01) | -0.245 | 0.004 | 0.202 | -0.638 | 0.150 | 3192.690 | 1.001 |
| SpeciesDOG:scale(BIO12) | -0.400 | 0.007 | 0.258 | -0.883 | 0.118 | 1246.663 | 1.001 |
| SpeciesFOX:scale(BIO12) | -0.852 | 0.006 | 0.246 | -1.337 | -0.380 | 1948.567 | 1.000 |
| SpeciesDOG:scale(BIO15) | 0.156 | 0.007 | 0.249 | -0.334 | 0.633 | 1262.307 | 1.001 |
| SpeciesFOX:scale(BIO15) | 0.792 | 0.005 | 0.236 | 0.334 | 1.255 | 1904.120 | 1.001 |
| sigma [1|Group] | 0.379 | 0.029 | 0.185 | 0.075 | 0.738 | 40.085 | 1.083 |
| sigma [1|Season] | 0.373 | 0.012 | 0.385 | 0.018 | 1.389 | 957.279 | 1.001 |
Show code
predator_summary_det <- summary(predator_global_ubms, submodel = "det")
predator_summary_det_sig <- which((predator_summary_det$`2.5%` * predator_summary_det$`97.5%`) > 0)
predator_summary_det%>%
dplyr::select(-`25%`, -`50%`, -`75%`) %>%
kbl("html", caption = "Model estimates showing the effect of covariates on the probability of detection for predator species. Estimates with a 95% CI not overlapping zero are shown in bold. The reference level (Intercept) is 'nonPred'.", digits = 3) %>%
kable_styling(full_width = FALSE) %>%
row_spec(predator_summary_det_sig, bold = TRUE)
| mean | se_mean | sd | 2.5% | 97.5% | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| (Intercept) | -4.008 | 0.003 | 0.119 | -4.247 | -3.780 | 1730.650 | 1.001 |
| SpeciesDOG | -0.705 | 0.006 | 0.272 | -1.276 | -0.211 | 2051.467 | 1.001 |
| SpeciesFOX | 0.233 | 0.004 | 0.184 | -0.139 | 0.586 | 2728.877 | 1.000 |
| baitMethodPred | 0.929 | 0.005 | 0.188 | 0.557 | 1.299 | 1538.164 | 1.004 |
| days | 0.001 | 0.000 | 0.004 | -0.007 | 0.010 | 2340.507 | 1.001 |
| baitMethodPred:days | -0.025 | 0.000 | 0.009 | -0.043 | -0.008 | 1769.780 | 1.004 |
Based on the \(\sigma\) values for the random effects (shown above), there is a non-zero estimated effect of group (a set of 3 neighboring sites) and season in explaining the overall variation. The estimates for the random effects can also be viewed and are as follows:
Show code
ran <- ubms::ranef(predator_global_ubms, submodel="state", summary=TRUE)
ran[["Group"]][[1]] %>%
kbl("html", caption = "Model estimates for random effect levels of 'Group' on the probability of occupancy for predator species.", digits = 3) %>%
kable_styling(full_width = FALSE) %>%
kableExtra::scroll_box(width = "800px", height = "600px")
| Estimate | SD | 2.5% | 97.5% | |
|---|---|---|---|---|
| G076 | -0.294 | 0.503 | -1.354 | 0.716 |
| G077 | -0.304 | 0.525 | -1.408 | 0.740 |
| G084 | -0.125 | 0.506 | -1.119 | 0.948 |
| G112 | 0.089 | 0.524 | -0.819 | 1.288 |
| G113 | -0.164 | 0.505 | -1.168 | 0.866 |
| G114 | -0.040 | 0.515 | -0.989 | 1.100 |
| G115 | -0.303 | 0.520 | -1.405 | 0.714 |
| G116 | -0.287 | 0.522 | -1.387 | 0.731 |
| G121 | -0.188 | 0.500 | -1.186 | 0.874 |
| G122 | -0.008 | 0.516 | -0.942 | 1.143 |
| G123 | -0.313 | 0.524 | -1.404 | 0.720 |
| G124 | -0.314 | 0.515 | -1.439 | 0.681 |
| G125 | -0.063 | 0.500 | -1.026 | 1.024 |
| G130 | -0.147 | 0.510 | -1.169 | 0.936 |
| G132 | -0.141 | 0.495 | -1.105 | 0.920 |
| G148 | -0.341 | 0.522 | -1.416 | 0.696 |
| G149 | -0.071 | 0.501 | -1.020 | 1.019 |
| G150 | 0.076 | 0.505 | -0.798 | 1.178 |
| G156 | -0.271 | 0.488 | -1.290 | 0.722 |
| G157 | -0.173 | 0.495 | -1.197 | 0.833 |
| G158 | 0.061 | 0.523 | -0.847 | 1.230 |
| G159 | -0.319 | 0.528 | -1.438 | 0.723 |
| G160 | -0.141 | 0.482 | -1.099 | 0.867 |
| G161 | -0.150 | 0.482 | -1.065 | 0.877 |
| G162 | -0.153 | 0.500 | -1.127 | 0.889 |
| G163 | -0.077 | 0.497 | -1.019 | 1.004 |
| G165 | 0.040 | 0.525 | -0.872 | 1.243 |
| G166 | -0.185 | 0.518 | -1.213 | 0.867 |
| G167 | -0.178 | 0.509 | -1.182 | 0.868 |
| G168 | -0.323 | 0.508 | -1.410 | 0.665 |
| G170 | -0.047 | 0.514 | -1.008 | 1.093 |
| G171 | -0.313 | 0.520 | -1.394 | 0.710 |
| G184 | -0.075 | 0.504 | -1.029 | 1.025 |
| G185 | -0.329 | 0.505 | -1.379 | 0.642 |
| G186 | -0.181 | 0.500 | -1.156 | 0.867 |
| G187 | 0.041 | 0.514 | -0.886 | 1.220 |
| G192 | -0.302 | 0.526 | -1.429 | 0.714 |
| G193 | -0.159 | 0.478 | -1.084 | 0.841 |
| G194 | -0.374 | 0.531 | -1.539 | 0.617 |
| G195 | -0.323 | 0.527 | -1.430 | 0.707 |
| G196 | -0.319 | 0.515 | -1.381 | 0.699 |
| G197 | -0.251 | 0.489 | -1.257 | 0.724 |
| G198 | -0.171 | 0.506 | -1.154 | 0.893 |
| G200 | -0.192 | 0.503 | -1.182 | 0.831 |
| G201 | -0.040 | 0.514 | -0.978 | 1.137 |
| G202 | -0.274 | 0.535 | -1.383 | 0.796 |
| G203 | -0.159 | 0.507 | -1.154 | 0.883 |
| G204 | -0.073 | 0.505 | -1.019 | 1.046 |
| G206 | -0.180 | 0.522 | -1.190 | 0.902 |
| G220 | -0.364 | 0.513 | -1.468 | 0.633 |
| G221 | -0.203 | 0.491 | -1.184 | 0.838 |
| G222 | 0.136 | 0.529 | -0.763 | 1.320 |
| G223 | 0.171 | 0.552 | -0.766 | 1.403 |
| G224 | -0.267 | 0.489 | -1.266 | 0.722 |
| G225 | -0.171 | 0.509 | -1.180 | 0.864 |
| G226 | -0.320 | 0.520 | -1.463 | 0.698 |
| G228 | -0.019 | 0.509 | -0.928 | 1.252 |
| G229 | -0.126 | 0.478 | -1.054 | 0.876 |
| G230 | -0.084 | 0.497 | -1.027 | 0.992 |
| G231 | -0.383 | 0.509 | -1.506 | 0.600 |
| G232 | -0.359 | 0.514 | -1.451 | 0.637 |
| G234 | -0.258 | 0.487 | -1.264 | 0.718 |
| G235 | -0.257 | 0.497 | -1.263 | 0.716 |
| G236 | -0.013 | 0.496 | -0.910 | 1.141 |
| G237 | -0.206 | 0.514 | -1.354 | 0.829 |
| G238 | -0.206 | 0.492 | -1.191 | 0.780 |
| G239 | -0.189 | 0.501 | -1.152 | 0.860 |
| G240 | 0.102 | 0.542 | -0.824 | 1.314 |
| G241 | -0.305 | 0.522 | -1.432 | 0.723 |
| G242 | -0.045 | 0.526 | -1.013 | 1.233 |
| G243 | -0.161 | 0.504 | -1.158 | 0.884 |
| G244 | -0.268 | 0.519 | -1.337 | 0.771 |
| G245 | -0.181 | 0.506 | -1.178 | 0.871 |
| G246 | -0.057 | 0.512 | -1.023 | 1.088 |
| G247 | -0.171 | 0.507 | -1.204 | 0.889 |
| G256 | -0.311 | 0.514 | -1.350 | 0.710 |
| G257 | -0.227 | 0.493 | -1.209 | 0.787 |
| G258 | 0.044 | 0.516 | -0.871 | 1.167 |
| G259 | -0.192 | 0.497 | -1.170 | 0.844 |
| G260 | -0.312 | 0.521 | -1.419 | 0.702 |
| G261 | -0.217 | 0.490 | -1.191 | 0.797 |
| G262 | -0.355 | 0.511 | -1.474 | 0.615 |
| G263 | -0.142 | 0.481 | -1.048 | 0.876 |
| G264 | -0.188 | 0.499 | -1.191 | 0.881 |
| G265 | -0.239 | 0.491 | -1.262 | 0.754 |
| G266 | 0.018 | 0.501 | -0.885 | 1.169 |
| G267 | -0.139 | 0.488 | -1.099 | 0.903 |
| G268 | -0.098 | 0.501 | -1.046 | 0.999 |
| G269 | -0.331 | 0.519 | -1.432 | 0.699 |
| G271 | -0.027 | 0.522 | -0.983 | 1.142 |
| G272 | 0.033 | 0.519 | -0.903 | 1.150 |
| G273 | -0.169 | 0.496 | -1.144 | 0.893 |
| G274 | -0.030 | 0.501 | -0.970 | 1.050 |
| G275 | -0.191 | 0.497 | -1.200 | 0.831 |
| G276 | -0.166 | 0.503 | -1.177 | 0.861 |
| G277 | -0.190 | 0.494 | -1.190 | 0.834 |
| G278 | -0.070 | 0.496 | -1.022 | 1.002 |
| G279 | -0.085 | 0.504 | -1.012 | 0.996 |
| G281 | -0.171 | 0.499 | -1.142 | 0.867 |
| G292 | -0.075 | 0.499 | -1.026 | 0.990 |
| G293 | -0.155 | 0.516 | -1.167 | 0.950 |
| G294 | -0.177 | 0.512 | -1.209 | 0.855 |
| G295 | -0.039 | 0.514 | -0.993 | 1.095 |
| G296 | -0.240 | 0.485 | -1.212 | 0.767 |
| G297 | -0.081 | 0.498 | -0.998 | 1.004 |
| G298 | 0.072 | 0.516 | -0.818 | 1.213 |
| G299 | -0.163 | 0.484 | -1.125 | 0.841 |
| G300 | -0.261 | 0.492 | -1.286 | 0.762 |
| G302 | -0.380 | 0.518 | -1.504 | 0.614 |
| G303 | -0.261 | 0.498 | -1.260 | 0.733 |
| G304 | -0.223 | 0.502 | -1.245 | 0.805 |
| G305 | -0.047 | 0.482 | -0.950 | 1.008 |
| G306 | -0.333 | 0.514 | -1.397 | 0.673 |
| G307 | -0.286 | 0.522 | -1.357 | 0.787 |
| G308 | -0.153 | 0.515 | -1.138 | 0.889 |
| G309 | -0.320 | 0.524 | -1.442 | 0.704 |
| G310 | -0.181 | 0.491 | -1.179 | 0.867 |
| G311 | -0.217 | 0.511 | -1.281 | 0.820 |
| G312 | -0.166 | 0.497 | -1.159 | 0.861 |
| G313 | -0.198 | 0.501 | -1.221 | 0.827 |
| G314 | -0.312 | 0.513 | -1.396 | 0.698 |
| G315 | -0.069 | 0.492 | -0.981 | 1.026 |
| G316 | -0.067 | 0.496 | -0.988 | 1.027 |
| G317 | -0.299 | 0.521 | -1.403 | 0.690 |
| G319 | -0.156 | 0.531 | -1.193 | 1.130 |
| G331 | 0.261 | 0.557 | -0.642 | 1.466 |
| G332 | -0.333 | 0.518 | -1.419 | 0.670 |
| G333 | 0.049 | 0.502 | -0.828 | 1.160 |
| G335 | -0.055 | 0.480 | -0.956 | 0.967 |
| G336 | -0.161 | 0.507 | -1.166 | 0.893 |
| G338 | -0.365 | 0.522 | -1.472 | 0.648 |
| G342 | -0.395 | 0.527 | -1.540 | 0.599 |
| G343 | -0.188 | 0.496 | -1.214 | 0.834 |
| G344 | -0.192 | 0.512 | -1.202 | 0.835 |
| G345 | -0.335 | 0.526 | -1.460 | 0.685 |
| G346 | -0.033 | 0.501 | -0.973 | 1.033 |
| G347 | -0.206 | 0.504 | -1.175 | 0.877 |
| G349 | -0.205 | 0.494 | -1.197 | 0.784 |
| G350 | -0.201 | 0.493 | -1.203 | 0.805 |
| G351 | -0.213 | 0.480 | -1.181 | 0.782 |
| G352 | -0.183 | 0.495 | -1.152 | 0.843 |
| G353 | -0.182 | 0.487 | -1.133 | 0.828 |
| G355 | -0.200 | 0.488 | -1.183 | 0.795 |
| G368 | 0.142 | 0.544 | -0.757 | 1.381 |
| G369 | -0.068 | 0.479 | -0.961 | 0.976 |
| G370 | -0.142 | 0.495 | -1.095 | 0.883 |
| G371 | -0.030 | 0.488 | -0.926 | 1.035 |
| G373 | 0.034 | 0.524 | -0.884 | 1.238 |
| G375 | -0.156 | 0.489 | -1.125 | 0.850 |
| G377 | -0.387 | 0.520 | -1.479 | 0.635 |
| G378 | -0.216 | 0.500 | -1.230 | 0.790 |
| G379 | 0.023 | 0.511 | -0.869 | 1.164 |
| G380 | -0.220 | 0.520 | -1.242 | 0.836 |
| G381 | -0.098 | 0.488 | -0.994 | 0.937 |
| G382 | -0.226 | 0.492 | -1.212 | 0.796 |
| G383 | -0.314 | 0.518 | -1.420 | 0.700 |
| G385 | -0.228 | 0.500 | -1.231 | 0.824 |
| G386 | -0.341 | 0.521 | -1.474 | 0.679 |
| G387 | -0.189 | 0.497 | -1.164 | 0.853 |
| G388 | -0.161 | 0.504 | -1.146 | 0.882 |
| G389 | -0.294 | 0.509 | -1.379 | 0.728 |
| G390 | -0.208 | 0.508 | -1.226 | 0.847 |
| G391 | -0.191 | 0.489 | -1.191 | 0.821 |
| G405 | -0.312 | 0.487 | -1.326 | 0.683 |
| G406 | 0.116 | 0.502 | -0.761 | 1.224 |
| G407 | -0.197 | 0.494 | -1.147 | 0.825 |
| G408 | -0.114 | 0.503 | -1.084 | 0.926 |
| G409 | -0.215 | 0.476 | -1.151 | 0.792 |
| G410 | -0.384 | 0.526 | -1.526 | 0.646 |
| G411 | -0.058 | 0.505 | -1.001 | 1.014 |
| G413 | -0.214 | 0.491 | -1.207 | 0.793 |
| G414 | -0.294 | 0.503 | -1.304 | 0.733 |
| G416 | -0.163 | 0.493 | -1.140 | 0.883 |
| G417 | -0.344 | 0.516 | -1.416 | 0.676 |
| G421 | -0.103 | 0.485 | -1.026 | 0.943 |
| G422 | -0.339 | 0.516 | -1.470 | 0.667 |
| G423 | -0.094 | 0.493 | -1.043 | 0.945 |
| G424 | -0.313 | 0.540 | -1.448 | 0.734 |
| G425 | -0.038 | 0.526 | -0.999 | 1.136 |
| G426 | -0.032 | 0.504 | -0.952 | 1.089 |
| G427 | -0.178 | 0.499 | -1.163 | 0.854 |
| G441 | -0.157 | 0.481 | -1.084 | 0.823 |
| G443 | 0.083 | 0.499 | -0.793 | 1.202 |
| G444 | -0.135 | 0.499 | -1.098 | 0.923 |
| G445 | -0.461 | 0.546 | -1.810 | 0.564 |
| G446 | -0.430 | 0.529 | -1.563 | 0.580 |
| G447 | -0.224 | 0.486 | -1.190 | 0.795 |
| G449 | 0.072 | 0.526 | -0.842 | 1.238 |
| G450 | -0.079 | 0.504 | -1.049 | 0.991 |
| G451 | -0.278 | 0.515 | -1.364 | 0.753 |
| G452 | -0.163 | 0.482 | -1.118 | 0.849 |
| G456 | -0.342 | 0.525 | -1.497 | 0.692 |
| G457 | -0.107 | 0.496 | -1.067 | 0.962 |
| G460 | -0.183 | 0.506 | -1.216 | 0.873 |
| G461 | -0.314 | 0.516 | -1.426 | 0.691 |
| G477 | -0.321 | 0.507 | -1.376 | 0.658 |
| G478 | -0.098 | 0.497 | -1.031 | 0.993 |
| G479 | -0.143 | 0.494 | -1.069 | 0.943 |
| G481 | -0.357 | 0.520 | -1.467 | 0.636 |
| G484 | -0.247 | 0.506 | -1.399 | 0.774 |
| G486 | -0.113 | 0.521 | -1.132 | 1.029 |
| G487 | -0.240 | 0.491 | -1.242 | 0.771 |
| G489 | -0.032 | 0.494 | -0.945 | 1.069 |
| G490 | -0.350 | 0.517 | -1.435 | 0.670 |
| G491 | 0.110 | 0.516 | -0.783 | 1.268 |
| G515 | -0.323 | 0.498 | -1.364 | 0.676 |
| G517 | -0.311 | 0.505 | -1.361 | 0.687 |
| G518 | -0.049 | 0.482 | -0.948 | 1.007 |
| G519 | -0.236 | 0.540 | -1.498 | 0.850 |
| G520 | -0.194 | 0.479 | -1.161 | 0.792 |
| G521 | 0.011 | 0.509 | -0.901 | 1.158 |
| G522 | -0.095 | 0.482 | -1.011 | 0.921 |
| G523 | -0.320 | 0.496 | -1.375 | 0.669 |
| G524 | -0.253 | 0.505 | -1.286 | 0.788 |
| G551 | -0.035 | 0.490 | -0.931 | 1.028 |
| G552 | -0.309 | 0.523 | -1.414 | 0.716 |
| G553 | -0.285 | 0.480 | -1.268 | 0.683 |
| G554 | 0.059 | 0.516 | -0.832 | 1.227 |
| G555 | -0.260 | 0.494 | -1.264 | 0.755 |
| G556 | 0.071 | 0.502 | -0.810 | 1.180 |
| G557 | -0.108 | 0.479 | -1.015 | 0.926 |
| G558 | -0.165 | 0.482 | -1.109 | 0.842 |
| G559 | -0.171 | 0.489 | -1.133 | 0.856 |
| G560 | -0.277 | 0.507 | -1.348 | 0.730 |
| G591 | 0.002 | 0.516 | -0.938 | 1.135 |
| G592 | -0.215 | 0.489 | -1.233 | 0.784 |
| G593 | -0.242 | 0.510 | -1.303 | 0.804 |
| G594 | -0.405 | 0.535 | -1.540 | 0.626 |
| G595 | -0.147 | 0.494 | -1.122 | 0.894 |
| G596 | -0.046 | 0.499 | -0.958 | 1.085 |
| G597 | -0.243 | 0.505 | -1.252 | 0.793 |
| G627 | -0.305 | 0.493 | -1.318 | 0.691 |
| G628 | 0.122 | 0.529 | -0.793 | 1.362 |
| G631 | -0.149 | 0.513 | -1.292 | 0.919 |
| G662 | -0.220 | 0.514 | -1.284 | 0.781 |
| G663 | -0.187 | 0.492 | -1.187 | 0.820 |
| G664 | -0.182 | 0.491 | -1.133 | 0.838 |
| G699 | -0.223 | 0.536 | -1.302 | 0.885 |
| G734 | -0.225 | 0.527 | -1.319 | 0.862 |
Show code
ran[["Season"]][[1]] %>%
kbl("html", caption = "Model estimates for random effect levels of 'Seasom' on the probability of occupancy for predator species.", digits = 3) %>%
kable_styling(full_width = FALSE)
| Estimate | SD | 2.5% | 97.5% | |
|---|---|---|---|---|
| 2 | -0.218 | 0.137 | -0.474 | 0.064 |
| 4 | -0.314 | 0.211 | -0.744 | 0.095 |
| 5 | 0.006 | 0.349 | -0.540 | 0.820 |
MCMC Convergence and Diagnostics
Given that we implemented a Bayesian approach here the following section presents several model diagnostics used to ensure model convergence and suitability. The traceplots show how that the chains have mixed well and the number of iterations used in the model is sufficient in reaching convergence.
Show code
ubms::traceplot(predator_global_ubms, pars=c("beta_state", "beta_det"), ncol = 3)
Figure 3: Traceplots of the four chains run for 3000 iterations in the global model. Plots show good evidence that chain-mixing and convergence has been achieved
Model Fit
Residuals
The residuals for the occupancy model can be calculated and plotted. The ubms package implements a method from Wright, Irvine, and Higgs (2019) to calculate residuals, where residuals are calculated independently for the observation and detection submodels. The diagnostic information from visualising these plots can be summarised as follows: if 95 % of the binned residual points fall within the grey shading then it is an indication the model fits the data well (Kellner 2021).
Show code
plot_residuals(predator_global_ubms, submodel="state")
Figure 4: Residual plot for the predator model, showing that binned residuals largely fall within the margins of error a model with good predictive power would have
Posterior Predictive Tests
We can predict the posterior distribution of the model and compare those predictions to the actual data. Here we show that the posterior prediction distribution matching the real data at a generalized level.
Show code
sim_y <- posterior_predict(predator_global_ubms, "y", draws=250)
prop0 <- apply(sim_y, 1, function(x) mean(x==0, na.rm=TRUE))
prop1 <- apply(sim_y, 1, function(x) mean(x==1, na.rm=TRUE))
actual_prop0 <- mean(getY(predator_global_ubms) == 0, na.rm=TRUE)
actual_prop1 <- mean(getY(predator_global_ubms) == 1, na.rm=TRUE)
#Compare 0
par(mfrow = c(1,2))
hist(prop0, col='gray', breaks = 15)
abline(v=actual_prop0, col='red', lwd=2)
#Compare 1
hist(prop1, col='gray', breaks = 15)
abline(v=actual_prop1, col='red', lwd=2)
Figure 5: Comparison of posterior predictions from the predator model to true mean values of absences and presences (0 and 1)
Predator Predictions
Using the predictions from the predator model we can generate site and season specific predictions of occupancy for cats, dogs and foxes (Figure (fig:predatorPredictions)). The distribution of these from the Bayesian ubms model is species dependent but relatively normal and suitable to use as a predictor variable for a native species model.
Show code
# At every site and season get mean predictions of predator occupancy
# predator_predictions_glmm <- bind_cols(predator_site_data,
# data.frame(Predicted = predict(predator_glmm,
# type = "response"))) %>%
# as.data.table() %>%
# dcast(Station + Season + FoxDensity ~ Species, value.var = "Predicted", fun.aggregate = mean)
predator_predictions_ubms <- bind_cols(predator_site_data,
predict(predator_global_ubms,
submodel = "state")) %>%
as.data.table() %>%
dcast(Station + Season + FoxDensity ~ Species, value.var = "Predicted", fun.aggregate = mean) %>%
mutate(FOX_scaled = scales::rescale(FOX),
FoxDen_scaled = scales::rescale(FoxDensity)) %>%
mutate(FoxScore = FOX_scaled + FoxDen_scaled)
pred_compare <- bind_cols(predator_site_data %>%
dplyr::select(Species, Station, Season, Visit),
predict(predator_global_ubms,
submodel = "state"))
pred_compare_plot <- pred_compare %>%
ggplot2::ggplot(aes(x = Predicted, y = Species)) +
ggridges::geom_density_ridges(aes(fill = Species), alpha = 0.75)+
scale_fill_manual(values = unname(delwp_cols[1:3]))+
# xlab("Predicted predator occupancy (ubms/bayesian model)") +
ylab("Predicted predator occupancy") +
# geom_point(shape = 21) +
# facet_wrap(~ Species, scales = "free") +
theme_minimal()
pred_compare_plot
Figure 6: Distribution of predator predicted occupancies for the ‘ubms’ model
Show code
# predator_predictions_ubms %>%
# ggplot2::ggplot(aes(x = FoxDensity, y = FOX)) +
# geom_point(shape = 21) +
# # facet_wrap(~ Species, scales = "free") +
# theme_minimal()
# # Comparison of predictions for
Occupancy Native Mammal Model: Stacked Approach
The below code shows how the data was curated into the unmarked frame for use within the native species model.
Show code
# combined probability
cp <- function(p)
{
ev <- do.call(expand.grid,replicate(length(p),0:1,simplify=FALSE))
pe <- apply(ev,1,function(x) prod(p*(x==1)+(1-p)*(x==0)))
tapply(pe,rowSums(ev),sum)
}
Natives <- c("LNP", "LNB", "LFP", "CBP", "GOA")
# Natives <- c("LNP","CBP")
# Restrict the model to data
natives_model_rows <- presence_data_2016_2020_mod %>%
filter(Species %in% Natives & !(rowID %in% no_geo)) %>% #Station %in% sites_visited_pre_post &
pull(rowID)
natives_site_data <- presence_data_2016_2020_mod %>%
filter(rowID %in% natives_model_rows) %>%
dplyr::select(Species, Station, Season, Visit) %>%
distinct() %>%
left_join(sf::st_drop_geometry(southernark_site_data)) %>%
left_join(predator_predictions_ubms) %>%
rowwise() %>%
mutate(TotalPred = sum(cp(c(FOX, CAT, DOG))[2:4])+FoxDensity) %>%
ungroup()
umf_natives <- unmarkedFrameOccu(y = presence_data_2016_2020_mod %>%
filter(rowID %in% natives_model_rows)%>%
dplyr::select(starts_with("D")),
siteCovs = natives_site_data,
obsCovs = list(days = days_data_2016_2020[natives_model_rows,],
baitMethod = bait_data_2016_2020[natives_model_rows,],
temp = temp_days_data_2016_2020[natives_model_rows,]))
Model Selection
The global model used for model selection is as follows, with all parameters having a fixed interaction with ‘Species’:
- fire severity (fixed)
- Fox Density: Because we had fox density estimates from a previously constructed predator model we oped to use the more informative fox density score for each site as opposed to the occupancy. This is justified as while the fox density and fox occupancy scores were not correlated statistically it would make it difficult for land managers to interpret both occupancy and density as seperate parameters. In essence, these variables are both proxy measures for the effect of fox predation on native species but fox density is more informative than fox occupancy. (fixed)
- cat occupancy
- dog occupancy
- BIO01
- BIO12
- BIO15
Elevation, BIO04 and latitude were removed as potential predictor variables due to their high correlation with mean temperature (BIO01). No two parameters were included when correlation coefficients exceeded 0.75 (Figure 7) except for cat occupancy and BIO01 (r=0.8) due to their key interest in this study.
Show code
# For our model we find that the variables of 'CAT', 'DOG' and 'FOX' are quite correlated and a better fit may drop some (e.g. cat and dog) or combine the occupancy rates to obtain a total predator score. Combining them into a single metric (TotalPred) is the conditional probability that a predator (fox, dog or cat) is present based on predicted occupancy plus the modelled fox abundance.
panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y, use="pairwise.complete.obs"))
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(natives_site_data[,c("scaled_fire_severity", "CAT", "DOG", "FOX", "FoxDensity", "BIO01", "BIO04", "BIO12", "BIO15", "elevation", "lat")], lower.panel=panel.smooth, upper.panel=panel.cor)
Figure 7: Correlation coefficients between possible predictor variables within the model
Show code
data_for_native_glmm <- natives_site_data %>%
left_join(summarised_detections_site_season, by = c("Station", "Species", "Season")) %>%
dplyr::mutate(GroupSpecies = paste0(Species, Group),
SeasonSpecies = paste0(Species,Season),
StationSpecies = paste0(Species,Station))
if(!file.exists("models/mod_dredge_native.rds")) {
global_model_natives <- lme4::glmer(Presence~Species
+ Species*scaled_fire_severity # species reactions to fire
+ Species*FoxDensity
+ Species*CAT
+ Species*DOG
+ Species*scale(BIO01) # Climate
# + Species*scale(BIO04) # Climate
+ Species*scale(BIO12) # Climate
+ Species*scale(BIO15) # Climate
+ (1|Season)
+ (1|Group), # Site specific effects\
family = binomial,
data = data_for_native_glmm,
control=lme4::glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)), na.action = "na.pass")
subsetExpr <- makeRule(CAT,DOG,scale(BIO01), scale(BIO12), scale(BIO15))
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", 8), type = clusterType))
clusterExport(clust, "data_for_native_glmm")
clusterExport(clust, c("glmer", "glmerControl"), envir = as.environment("package:lme4"))
mod_dredge_native <- MuMIn::pdredge(global_model_natives, cluster = clust,
fixed = ~ Species:scaled_fire_severity + Species:FoxDensity + (1|Season) + (1|Group),
subset = subsetExpr)
top_model_natives <- MuMIn::get.models(mod_dredge_native, subset = 1)[[1]]
saveRDS(mod_dredge_native, "models/mod_dredge_native.rds", compress = "xz")
saveRDS(global_model_natives, "models/global_model_natives.rds", compress = "xz")
saveRDS(top_model_natives, "models/top_model_natives.rds", compress = "xz")
} else {
mod_dredge_native <- readRDS("models/mod_dredge_native.rds")
global_model_natives <- readRDS("models/global_model_natives.rds")
top_model_natives <- readRDS("models/top_model_natives.rds")
}
mod_dredge_native %>%
as.data.frame() %>%
dplyr::select(-`(Intercept)`, -`scale(BIO01)`, -`scale(BIO12)`, -`scale(BIO15)`, -`scaled_fire_severity`, -CAT, -DOG, -FoxDensity) %>%
kbl(format = "html", caption = "Model selection for the predator models using GLMM", digits = 3, escape = FALSE) %>%
kable_styling(full_width = FALSE) %>%
row_spec(2, bold = TRUE) %>%
scroll_box(width = "800px", height = "800px")
| Species | CAT:Species | DOG:Species | FoxDensity:Species | scale(BIO01):Species | scale(BIO12):Species | scale(BIO15):Species | scaled_fire_severity:Species | df | logLik | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 4063 |
|
|
|
|
|
|
NA |
|
37 | -2291.328 | 4657.140 | 0.000 | 7.207261e-01 |
| 3549 |
|
|
NA |
|
|
|
NA |
|
32 | -2298.352 | 4661.068 | 3.928 | 1.011058e-01 |
| 8191 |
|
|
|
|
|
|
|
|
42 | -2288.760 | 4662.142 | 5.002 | 5.910985e-02 |
| 3292 |
|
NA | NA |
|
|
|
NA |
|
27 | -2304.070 | 4662.399 | 5.259 | 5.197185e-02 |
| 1999 |
|
|
|
|
|
NA | NA |
|
32 | -2300.053 | 4664.469 | 7.329 | 1.846002e-02 |
| 3806 |
|
NA |
|
|
|
|
NA |
|
32 | -2300.197 | 4664.756 | 7.617 | 1.599054e-02 |
| 7420 |
|
NA | NA |
|
|
|
|
|
32 | -2300.496 | 4665.355 | 8.215 | 1.185581e-02 |
| 7677 |
|
|
NA |
|
|
|
|
|
37 | -2295.511 | 4665.506 | 8.366 | 1.099200e-02 |
| 6127 |
|
|
|
|
|
NA |
|
|
37 | -2296.163 | 4666.810 | 9.670 | 5.726733e-03 |
| 7934 |
|
NA |
|
|
|
|
|
|
37 | -2297.198 | 4668.880 | 11.741 | 2.033951e-03 |
| 1742 |
|
NA |
|
|
|
NA | NA |
|
27 | -2307.519 | 4669.297 | 12.157 | 1.651515e-03 |
| 5870 |
|
NA |
|
|
|
NA |
|
|
32 | -2303.947 | 4672.258 | 15.118 | 3.758378e-04 |
| 1485 |
|
|
NA |
|
|
NA | NA |
|
27 | -2321.360 | 4696.980 | 39.840 | 1.609195e-09 |
| 5356 |
|
NA | NA |
|
|
NA |
|
|
27 | -2321.586 | 4697.432 | 40.292 | 1.283935e-09 |
| 5613 |
|
|
NA |
|
|
NA |
|
|
32 | -2316.855 | 4698.073 | 40.933 | 9.317288e-10 |
| 1228 |
|
NA | NA |
|
|
NA | NA |
|
22 | -2327.445 | 4699.064 | 41.924 | 5.675608e-10 |
| 7159 |
|
|
|
|
NA |
|
|
|
37 | -2314.730 | 4703.945 | 46.805 | 4.946162e-11 |
| 3031 |
|
|
|
|
NA |
|
NA |
|
32 | -2327.202 | 4718.767 | 61.627 | 2.990215e-14 |
| 967 |
|
|
|
|
NA | NA | NA |
|
27 | -2337.124 | 4728.508 | 71.368 | 2.292953e-16 |
| 2517 |
|
|
NA |
|
NA |
|
NA |
|
27 | -2340.149 | 4734.557 | 77.417 | 1.113923e-17 |
| 5095 |
|
|
|
|
NA | NA |
|
|
32 | -2335.571 | 4735.506 | 78.366 | 6.930822e-18 |
| 6645 |
|
|
NA |
|
NA |
|
|
|
32 | -2338.338 | 4741.040 | 83.900 | 4.356241e-19 |
| 453 |
|
|
NA |
|
NA | NA | NA |
|
22 | -2362.241 | 4768.656 | 111.516 | 4.389630e-25 |
| 4581 |
|
|
NA |
|
NA | NA |
|
|
27 | -2358.224 | 4770.707 | 113.567 | 1.574050e-25 |
| 6902 |
|
NA |
|
|
NA |
|
|
|
32 | -2394.030 | 4852.424 | 195.284 | 2.834213e-43 |
| 2774 |
|
NA |
|
|
NA |
|
NA |
|
27 | -2442.223 | 4938.705 | 281.565 | 5.207356e-62 |
| 6388 |
|
NA | NA |
|
NA |
|
|
|
27 | -2453.462 | 4961.183 | 304.044 | 6.847864e-67 |
| 2260 |
|
NA | NA |
|
NA |
|
NA |
|
22 | -2490.758 | 5025.689 | 368.549 | 6.735199e-81 |
| 4838 |
|
NA |
|
|
NA | NA |
|
|
27 | -2501.486 | 5057.232 | 400.092 | 9.523846e-88 |
| 710 |
|
NA |
|
|
NA | NA | NA |
|
22 | -2509.880 | 5063.934 | 406.794 | 3.339328e-89 |
| 4324 |
|
NA | NA |
|
NA | NA |
|
|
22 | -2513.657 | 5071.488 | 414.348 | 7.644102e-91 |
| 196 |
|
NA | NA |
|
NA | NA | NA |
|
17 | -2527.069 | 5088.242 | 431.102 | 1.758267e-94 |
From the above model selection (Table 7) we see that generally additional terms are not penalised (more than they improve the fit) and that the best model fit (according to AICc). We chose the top model to use for the occupancy model. Based on a GLMM this top model was a relatively good ROC indicating a decent level of predictive power (Figure ??)
Show code
# LFP_rows <- which(data_for_native_glmm$Species == "LNB")
rocobj <- plot.roc(data_for_native_glmm$Presence, predict(top_model_natives),
main="Confidence intervals", percent=TRUE,
ci=TRUE, # compute AUC (of AUC by default)
print.auc=TRUE) # print the AUC (will contain the CI)
ciobj <- ci.se(rocobj, # CI of sensitivity
specificities=seq(0, 100, 5)) # over a select set of specificities
plot(ciobj, type="shape", col="#1c61b6AA") # plot as a blue shape
plot(ci(rocobj, of="thresholds", thresholds="best")) # add one threshold
Figure 8: Response to the Operating Curve (ROC) for the top GLMM, indicating high degrees of specificity and sensitivity, suggesting a good model fit of the top model.
Chosen Model (Bayesian)
Detection Submodel: ~ Baiting Method x Days + Species*Temperature Occupancy Submodel: ~ (Species x Landscape Fire Severity) + (Species x Fox Density) + (Species x Cat Occupancy) + (Species x Dog Occupancy) + (Species x BIO01) + (Species x BIO12) + (1|Group) + (1|Season)
Show code
# options(mc.cores=8) # increase for this model
if(!file.exists("models/natives_global_ubms.rds")) {
natives_global_ubms <- stan_occu(~Species*temp+baitMethod*days
~Species
+ Species*scaled_fire_severity # species reactions to fire
+ Species*FoxDensity
+ Species*CAT
+ Species*DOG
+ Species*scale(BIO01) # Climate
+ Species*scale(BIO12) # Climate
+ (1|Season)
+ (1|Group), # Seasonal effects
data = umf_natives, chains=4, iter=4000)
saveRDS(natives_global_ubms, "models/natives_global_ubms.rds", compress = "xz")
} else {
natives_global_ubms <- readRDS("models/natives_global_ubms.rds")
}
# options(mc.cores=4) # reset for this model
Model Summary
The summary for the top model (global) is broken up into the Occupancy sub-model and the Detection sub-model. The estimates for the model parameters, standard deviation (SD), 95% confidence intervals, effective sample size and Rhat is presented below. n_eff > 200 and Rhat < 1.05 are general indicators the model has reached convergence.
Show code
native_summary <- summary(natives_global_ubms, submodel = "state")
native_summary_sig <- which((native_summary$`2.5%` * native_summary$`97.5%`) > 0)
native_summary%>%
dplyr::select(-`25%`, -`50%`, -`75%`) %>%
kbl("html", caption = "Model estimates showing the effects of various covariates on the probability of occupancy for native species. Estimates with a 95% CI not overlapping zero are shown in bold. The reference level (Intercept) is 'BTP'.", digits = 3) %>%
kable_styling(full_width = FALSE) %>%
row_spec(native_summary_sig, bold = TRUE)
| mean | se_mean | sd | 2.5% | 97.5% | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| (Intercept) | -1.564 | 0.008 | 0.603 | -2.753 | -0.419 | 5863.509 | 1.000 |
| SpeciesGOA | 0.304 | 0.007 | 0.639 | -0.913 | 1.635 | 8062.450 | 1.000 |
| SpeciesLFP | 1.350 | 0.008 | 0.649 | 0.019 | 2.559 | 6385.894 | 1.001 |
| SpeciesLNB | -0.314 | 0.008 | 0.632 | -1.627 | 0.906 | 6878.865 | 1.000 |
| SpeciesLNP | -2.104 | 0.009 | 0.774 | -3.602 | -0.540 | 7948.920 | 1.001 |
| scaled_fire_severity | -0.154 | 0.006 | 0.500 | -1.132 | 0.801 | 7246.667 | 1.000 |
| FoxDensity | 0.506 | 0.003 | 0.297 | -0.087 | 1.087 | 7642.236 | 1.000 |
| CAT | 0.689 | 0.015 | 1.179 | -1.580 | 3.047 | 6441.594 | 1.000 |
| DOG | -0.021 | 0.014 | 1.467 | -2.951 | 2.936 | 11511.676 | 1.000 |
| scale(BIO01) | 0.516 | 0.002 | 0.134 | 0.253 | 0.789 | 6250.900 | 1.000 |
| scale(BIO12) | -0.410 | 0.001 | 0.101 | -0.607 | -0.213 | 6984.902 | 1.000 |
| SpeciesGOA:scaled_fire_severity | -5.333 | 0.011 | 1.175 | -7.869 | -3.242 | 11698.691 | 1.000 |
| SpeciesLFP:scaled_fire_severity | -0.966 | 0.005 | 0.467 | -1.881 | -0.044 | 7671.095 | 1.000 |
| SpeciesLNB:scaled_fire_severity | -1.891 | 0.007 | 0.636 | -3.169 | -0.671 | 8975.476 | 1.000 |
| SpeciesLNP:scaled_fire_severity | -2.765 | 0.014 | 1.341 | -5.736 | -0.486 | 9654.119 | 1.000 |
| SpeciesGOA:FoxDensity | 0.157 | 0.004 | 0.432 | -0.663 | 1.023 | 9835.799 | 1.000 |
| SpeciesLFP:FoxDensity | -8.654 | 0.008 | 0.936 | -10.521 | -6.894 | 12148.437 | 1.000 |
| SpeciesLNB:FoxDensity | -0.865 | 0.005 | 0.462 | -1.811 | 0.021 | 9822.855 | 1.001 |
| SpeciesLNP:FoxDensity | -0.575 | 0.007 | 0.753 | -2.166 | 0.771 | 11304.681 | 1.000 |
| SpeciesGOA:CAT | -0.591 | 0.015 | 1.389 | -3.427 | 2.088 | 8419.263 | 1.000 |
| SpeciesLFP:CAT | 1.957 | 0.018 | 1.473 | -0.717 | 5.126 | 7024.700 | 1.000 |
| SpeciesLNB:CAT | 0.913 | 0.016 | 1.376 | -1.695 | 3.729 | 7773.338 | 1.000 |
| SpeciesLNP:CAT | -0.576 | 0.017 | 1.611 | -3.970 | 2.483 | 8763.270 | 1.001 |
| SpeciesGOA:DOG | 0.182 | 0.015 | 1.733 | -3.259 | 3.748 | 13890.089 | 1.000 |
| SpeciesLFP:DOG | -1.351 | 0.018 | 1.808 | -5.348 | 1.803 | 10523.167 | 1.000 |
| SpeciesLNB:DOG | 0.861 | 0.015 | 1.715 | -2.296 | 4.611 | 12588.329 | 1.000 |
| SpeciesLNP:DOG | -0.183 | 0.015 | 1.752 | -3.782 | 3.274 | 13793.131 | 1.000 |
| SpeciesGOA:scale(BIO01) | 0.739 | 0.003 | 0.218 | 0.307 | 1.164 | 7017.232 | 1.000 |
| SpeciesLFP:scale(BIO01) | -1.379 | 0.002 | 0.172 | -1.715 | -1.039 | 6515.736 | 1.000 |
| SpeciesLNB:scale(BIO01) | -0.077 | 0.002 | 0.178 | -0.421 | 0.276 | 7485.517 | 1.000 |
| SpeciesLNP:scale(BIO01) | 0.566 | 0.004 | 0.352 | -0.086 | 1.288 | 9269.766 | 1.000 |
| SpeciesGOA:scale(BIO12) | 0.542 | 0.001 | 0.144 | 0.263 | 0.825 | 9240.857 | 1.000 |
| SpeciesLFP:scale(BIO12) | 0.391 | 0.001 | 0.125 | 0.139 | 0.629 | 7510.863 | 1.000 |
| SpeciesLNB:scale(BIO12) | 0.607 | 0.001 | 0.137 | 0.343 | 0.879 | 8706.726 | 1.000 |
| SpeciesLNP:scale(BIO12) | 0.887 | 0.002 | 0.226 | 0.452 | 1.338 | 9819.703 | 1.000 |
| sigma [1|Group] | 0.758 | 0.010 | 0.592 | 0.051 | 2.269 | 3302.146 | 1.002 |
| sigma [1|Season] | 0.541 | 0.002 | 0.066 | 0.415 | 0.676 | 1232.133 | 1.001 |
Show code
native_summary_det <- summary(natives_global_ubms, submodel = "det")
native_summary_sig_det <- which((native_summary_det$`2.5%` * native_summary_det$`97.5%`) > 0)
native_summary_det %>%
dplyr::select(-`25%`, -`50%`, -`75%`) %>%
kbl("html", caption = "Model estimates showing the effects of covariates on the probability of detection for native species. Estimates with a 95% CI not overlapping zero are shown in bold. The reference level (Intercept) is 'nonPred'.", digits = 3) %>%
kable_styling(full_width = FALSE) %>%
row_spec(native_summary_sig_det, bold = TRUE)
| mean | se_mean | sd | 2.5% | 97.5% | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| (Intercept) | -1.983 | 0.002 | 0.115 | -2.209 | -1.750 | 4191.553 | 1.000 |
| SpeciesGOA | -1.365 | 0.003 | 0.249 | -1.860 | -0.874 | 7012.408 | 1.000 |
| SpeciesLFP | 0.828 | 0.002 | 0.130 | 0.575 | 1.085 | 4253.309 | 1.001 |
| SpeciesLNB | -1.573 | 0.003 | 0.225 | -2.008 | -1.130 | 6835.344 | 1.000 |
| SpeciesLNP | -0.600 | 0.004 | 0.381 | -1.377 | 0.134 | 8460.647 | 1.000 |
| temp | -0.003 | 0.000 | 0.006 | -0.014 | 0.008 | 4153.719 | 1.001 |
| baitMethodPred | -0.559 | 0.001 | 0.090 | -0.737 | -0.386 | 11110.226 | 1.000 |
| days | -0.015 | 0.000 | 0.001 | -0.018 | -0.012 | 9598.819 | 1.000 |
| SpeciesGOA:temp | 0.017 | 0.000 | 0.011 | -0.005 | 0.039 | 6478.581 | 1.001 |
| SpeciesLFP:temp | -0.008 | 0.000 | 0.007 | -0.020 | 0.005 | 3968.598 | 1.001 |
| SpeciesLNB:temp | 0.036 | 0.000 | 0.011 | 0.016 | 0.057 | 6742.782 | 1.000 |
| SpeciesLNP:temp | 0.027 | 0.000 | 0.018 | -0.008 | 0.062 | 8328.035 | 1.000 |
| baitMethodPred:days | 0.005 | 0.000 | 0.004 | -0.003 | 0.013 | 11733.468 | 1.000 |
Based on the \(\sigma\) values for the random effects (shown above), there is an effect of group (a set of 3 neighboring sites) and season in explaining the overall variation. The estimates for the random effects can also be viewed and are as follows:
Show code
ran_n <- ubms::ranef(natives_global_ubms, submodel="state", summary=TRUE)
ran_n[["Group"]][[1]] %>%
kbl("html", caption = "Model estimates for random effect levels of 'Group' on the probability of occupancy for native species.", digits = 3) %>%
kable_styling(full_width = FALSE) %>%
kableExtra::scroll_box(width = "800px", height = "600px")
| Estimate | SD | 2.5% | 97.5% | |
|---|---|---|---|---|
| G076 | -1.312 | 0.735 | -2.724 | 0.157 |
| G077 | -1.483 | 0.732 | -2.949 | -0.067 |
| G084 | -2.355 | 0.987 | -4.494 | -0.683 |
| G112 | -1.567 | 0.750 | -3.050 | -0.134 |
| G113 | -1.585 | 0.729 | -3.029 | -0.161 |
| G114 | -1.935 | 0.743 | -3.399 | -0.488 |
| G115 | -1.573 | 0.720 | -3.008 | -0.196 |
| G116 | -1.524 | 0.721 | -2.960 | -0.142 |
| G121 | -2.205 | 0.735 | -3.683 | -0.820 |
| G122 | -1.722 | 0.726 | -3.148 | -0.302 |
| G123 | -1.848 | 0.717 | -3.255 | -0.465 |
| G124 | -1.738 | 0.708 | -3.134 | -0.349 |
| G125 | -1.473 | 0.721 | -2.909 | -0.096 |
| G130 | -1.025 | 0.728 | -2.464 | 0.365 |
| G132 | -1.396 | 0.712 | -2.803 | -0.025 |
| G148 | -1.592 | 0.719 | -2.996 | -0.196 |
| G149 | -1.879 | 0.743 | -3.344 | -0.455 |
| G150 | -1.500 | 0.743 | -2.982 | -0.057 |
| G156 | -1.415 | 0.688 | -2.748 | -0.073 |
| G157 | -1.343 | 0.734 | -2.815 | 0.093 |
| G158 | -1.847 | 0.739 | -3.312 | -0.397 |
| G159 | -1.257 | 0.708 | -2.649 | 0.119 |
| G160 | -1.229 | 0.699 | -2.609 | 0.110 |
| G161 | -1.788 | 0.697 | -3.177 | -0.442 |
| G162 | -1.993 | 0.730 | -3.428 | -0.583 |
| G163 | -1.652 | 0.728 | -3.116 | -0.268 |
| G165 | -1.679 | 0.733 | -3.135 | -0.278 |
| G166 | -1.950 | 0.736 | -3.417 | -0.549 |
| G167 | -1.802 | 0.719 | -3.236 | -0.402 |
| G168 | -1.965 | 0.718 | -3.394 | -0.601 |
| G170 | -2.282 | 0.749 | -3.788 | -0.847 |
| G171 | -1.920 | 0.723 | -3.347 | -0.520 |
| G184 | -1.398 | 0.759 | -2.886 | 0.095 |
| G185 | -1.039 | 0.720 | -2.460 | 0.343 |
| G186 | -1.239 | 0.728 | -2.695 | 0.144 |
| G187 | -1.343 | 0.754 | -2.852 | 0.120 |
| G192 | -1.467 | 0.715 | -2.896 | -0.079 |
| G193 | -1.394 | 0.697 | -2.766 | -0.046 |
| G194 | -1.570 | 0.670 | -2.926 | -0.288 |
| G195 | -1.495 | 0.716 | -2.916 | -0.116 |
| G196 | -1.671 | 0.715 | -3.088 | -0.300 |
| G197 | -1.202 | 0.691 | -2.587 | 0.129 |
| G198 | -1.524 | 0.720 | -2.972 | -0.123 |
| G200 | -1.616 | 0.717 | -3.026 | -0.217 |
| G201 | -1.753 | 0.738 | -3.217 | -0.307 |
| G202 | -1.975 | 0.718 | -3.416 | -0.581 |
| G203 | -2.079 | 0.739 | -3.553 | -0.675 |
| G204 | -1.934 | 0.734 | -3.401 | -0.503 |
| G206 | -2.129 | 0.737 | -3.637 | -0.723 |
| G220 | -1.653 | 0.742 | -3.119 | -0.200 |
| G221 | -1.141 | 0.730 | -2.546 | 0.306 |
| G222 | -1.464 | 0.762 | -2.959 | 0.011 |
| G223 | -1.761 | 0.753 | -3.238 | -0.289 |
| G224 | -1.709 | 0.694 | -3.111 | -0.384 |
| G225 | -1.933 | 0.730 | -3.395 | -0.502 |
| G226 | -1.454 | 0.710 | -2.860 | -0.070 |
| G228 | -1.015 | 0.708 | -2.451 | 0.337 |
| G229 | -0.871 | 0.701 | -2.241 | 0.482 |
| G230 | -1.239 | 0.707 | -2.627 | 0.143 |
| G231 | -1.421 | 0.691 | -2.806 | -0.103 |
| G232 | -1.408 | 0.681 | -2.757 | -0.084 |
| G234 | -1.686 | 0.701 | -3.080 | -0.327 |
| G235 | -1.401 | 0.683 | -2.768 | -0.093 |
| G236 | -1.442 | 0.715 | -2.865 | -0.053 |
| G237 | -1.802 | 0.726 | -3.247 | -0.401 |
| G238 | -1.995 | 0.728 | -3.488 | -0.625 |
| G239 | -2.116 | 0.738 | -3.614 | -0.675 |
| G240 | -1.500 | 0.735 | -2.966 | -0.055 |
| G241 | -1.246 | 0.705 | -2.638 | 0.102 |
| G242 | -1.766 | 0.731 | -3.241 | -0.354 |
| G243 | -1.770 | 0.728 | -3.189 | -0.348 |
| G244 | -2.065 | 0.757 | -3.606 | -0.585 |
| G245 | -1.768 | 0.725 | -3.199 | -0.346 |
| G246 | -2.299 | 0.753 | -3.778 | -0.850 |
| G247 | -2.151 | 0.746 | -3.655 | -0.729 |
| G256 | -1.667 | 0.752 | -3.146 | -0.201 |
| G257 | -1.332 | 0.748 | -2.760 | 0.107 |
| G258 | -1.405 | 0.753 | -2.878 | 0.051 |
| G259 | -1.619 | 0.719 | -3.020 | -0.239 |
| G260 | -1.163 | 0.710 | -2.542 | 0.232 |
| G261 | -0.967 | 0.715 | -2.376 | 0.389 |
| G262 | -0.950 | 0.706 | -2.337 | 0.430 |
| G263 | -1.050 | 0.712 | -2.467 | 0.337 |
| G264 | -1.487 | 0.718 | -2.918 | -0.095 |
| G265 | -1.116 | 0.685 | -2.471 | 0.231 |
| G266 | -1.284 | 0.728 | -2.716 | 0.129 |
| G267 | -1.599 | 0.715 | -3.042 | -0.232 |
| G268 | -1.130 | 0.711 | -2.527 | 0.240 |
| G269 | -1.825 | 0.719 | -3.248 | -0.449 |
| G271 | -1.506 | 0.730 | -2.943 | -0.093 |
| G272 | -1.690 | 0.731 | -3.131 | -0.247 |
| G273 | -1.764 | 0.731 | -3.217 | -0.348 |
| G274 | -1.629 | 0.741 | -3.102 | -0.184 |
| G275 | -2.079 | 0.737 | -3.549 | -0.675 |
| G276 | -1.778 | 0.722 | -3.198 | -0.381 |
| G277 | -2.144 | 0.736 | -3.610 | -0.750 |
| G278 | -1.630 | 0.735 | -3.097 | -0.214 |
| G279 | -1.948 | 0.734 | -3.415 | -0.532 |
| G281 | -1.672 | 0.750 | -3.167 | -0.211 |
| G292 | -1.821 | 0.766 | -3.311 | -0.336 |
| G293 | -1.340 | 0.755 | -2.842 | 0.140 |
| G294 | -1.488 | 0.753 | -2.963 | -0.028 |
| G295 | -1.636 | 0.739 | -3.112 | -0.214 |
| G296 | -1.202 | 0.694 | -2.580 | 0.114 |
| G297 | -1.792 | 0.738 | -3.242 | -0.368 |
| G298 | -0.769 | 0.736 | -2.219 | 0.681 |
| G299 | -1.341 | 0.709 | -2.734 | 0.038 |
| G300 | -1.128 | 0.689 | -2.504 | 0.179 |
| G302 | -1.316 | 0.724 | -2.751 | 0.110 |
| G303 | -1.687 | 0.707 | -3.101 | -0.304 |
| G304 | -1.442 | 0.705 | -2.834 | -0.073 |
| G305 | -1.904 | 0.730 | -3.382 | -0.515 |
| G306 | -1.557 | 0.701 | -2.956 | -0.203 |
| G307 | -1.645 | 0.712 | -3.047 | -0.276 |
| G308 | -2.319 | 0.736 | -3.806 | -0.899 |
| G309 | -1.912 | 0.734 | -3.374 | -0.535 |
| G310 | -1.246 | 0.707 | -2.621 | 0.144 |
| G311 | -1.363 | 0.718 | -2.793 | 0.016 |
| G312 | -1.400 | 0.721 | -2.853 | -0.002 |
| G313 | -1.855 | 0.731 | -3.289 | -0.454 |
| G314 | -2.045 | 0.722 | -3.434 | -0.644 |
| G315 | -1.764 | 0.739 | -3.261 | -0.298 |
| G316 | -2.210 | 0.760 | -3.728 | -0.741 |
| G317 | -2.098 | 0.742 | -3.596 | -0.679 |
| G319 | -1.764 | 0.748 | -3.267 | -0.326 |
| G331 | -1.278 | 0.753 | -2.775 | 0.201 |
| G332 | -1.377 | 0.709 | -2.765 | -0.027 |
| G333 | -1.162 | 0.732 | -2.621 | 0.259 |
| G335 | -0.941 | 0.730 | -2.346 | 0.474 |
| G336 | -0.823 | 0.725 | -2.262 | 0.581 |
| G338 | -1.340 | 0.721 | -2.763 | 0.047 |
| G342 | -1.614 | 0.695 | -2.992 | -0.282 |
| G343 | -1.230 | 0.719 | -2.634 | 0.170 |
| G344 | -1.769 | 0.733 | -3.206 | -0.377 |
| G345 | -1.737 | 0.705 | -3.122 | -0.386 |
| G346 | -1.442 | 0.735 | -2.881 | -0.007 |
| G347 | -1.877 | 0.748 | -3.320 | -0.461 |
| G349 | -2.163 | 0.746 | -3.650 | -0.729 |
| G350 | -1.689 | 0.732 | -3.120 | -0.267 |
| G351 | -1.370 | 0.729 | -2.818 | 0.035 |
| G352 | -1.867 | 0.748 | -3.373 | -0.431 |
| G353 | -1.580 | 0.734 | -3.027 | -0.176 |
| G355 | -1.521 | 0.729 | -2.951 | -0.093 |
| G368 | -1.112 | 0.744 | -2.588 | 0.324 |
| G369 | -0.497 | 0.727 | -1.913 | 0.921 |
| G370 | -0.728 | 0.719 | -2.170 | 0.649 |
| G371 | -1.183 | 0.717 | -2.600 | 0.206 |
| G373 | -1.045 | 0.730 | -2.508 | 0.346 |
| G375 | -1.357 | 0.728 | -2.811 | 0.063 |
| G377 | -1.297 | 0.689 | -2.668 | 0.070 |
| G378 | -1.963 | 0.729 | -3.414 | -0.558 |
| G379 | -1.629 | 0.735 | -3.072 | -0.196 |
| G380 | -1.356 | 0.725 | -2.785 | 0.039 |
| G381 | -1.387 | 0.724 | -2.823 | -0.013 |
| G382 | -1.499 | 0.722 | -2.922 | -0.098 |
| G383 | -1.758 | 0.726 | -3.218 | -0.358 |
| G385 | -2.002 | 0.760 | -3.517 | -0.540 |
| G386 | -2.135 | 0.754 | -3.665 | -0.687 |
| G387 | -2.062 | 0.740 | -3.563 | -0.628 |
| G388 | -1.779 | 0.732 | -3.227 | -0.370 |
| G389 | -1.855 | 0.746 | -3.346 | -0.409 |
| G390 | -0.987 | 0.731 | -2.424 | 0.398 |
| G391 | -1.401 | 0.741 | -2.855 | 0.034 |
| G405 | -0.775 | 0.681 | -2.110 | 0.554 |
| G406 | -0.457 | 0.738 | -1.891 | 0.969 |
| G407 | -1.224 | 0.715 | -2.627 | 0.174 |
| G408 | -1.658 | 0.722 | -3.089 | -0.272 |
| G409 | -1.423 | 0.694 | -2.807 | -0.109 |
| G410 | -1.467 | 0.691 | -2.851 | -0.138 |
| G411 | -1.503 | 0.736 | -2.991 | -0.083 |
| G413 | -1.076 | 0.691 | -2.443 | 0.260 |
| G414 | -1.551 | 0.710 | -2.957 | -0.177 |
| G416 | -1.520 | 0.713 | -2.938 | -0.127 |
| G417 | -1.399 | 0.736 | -2.835 | 0.011 |
| G421 | -1.449 | 0.745 | -2.936 | -0.034 |
| G422 | -1.544 | 0.730 | -2.984 | -0.126 |
| G423 | -1.541 | 0.743 | -3.034 | -0.109 |
| G424 | -1.629 | 0.739 | -3.090 | -0.204 |
| G425 | -1.803 | 0.756 | -3.313 | -0.350 |
| G426 | -1.389 | 0.765 | -2.909 | 0.131 |
| G427 | -1.326 | 0.740 | -2.788 | 0.123 |
| G441 | -1.002 | 0.728 | -2.441 | 0.428 |
| G443 | -1.047 | 0.739 | -2.539 | 0.391 |
| G444 | -0.759 | 0.724 | -2.195 | 0.652 |
| G445 | -1.238 | 0.672 | -2.546 | 0.087 |
| G446 | -1.361 | 0.685 | -2.732 | -0.055 |
| G447 | -1.710 | 0.734 | -3.168 | -0.285 |
| G449 | -1.236 | 0.737 | -2.675 | 0.207 |
| G450 | -1.493 | 0.750 | -3.004 | -0.060 |
| G451 | -1.287 | 0.744 | -2.734 | 0.143 |
| G452 | -1.781 | 0.740 | -3.217 | -0.367 |
| G456 | -1.890 | 0.747 | -3.392 | -0.430 |
| G457 | -1.821 | 0.760 | -3.370 | -0.360 |
| G460 | -1.500 | 0.740 | -2.975 | -0.069 |
| G461 | -1.593 | 0.728 | -3.053 | -0.167 |
| G477 | -1.089 | 0.705 | -2.444 | 0.294 |
| G478 | -1.360 | 0.741 | -2.827 | 0.056 |
| G479 | -0.933 | 0.726 | -2.384 | 0.485 |
| G481 | -1.968 | 0.738 | -3.478 | -0.548 |
| G484 | -1.150 | 0.737 | -2.600 | 0.283 |
| G486 | -1.864 | 0.766 | -3.380 | -0.404 |
| G487 | -1.800 | 0.720 | -3.197 | -0.374 |
| G489 | -1.691 | 0.738 | -3.152 | -0.245 |
| G490 | -1.061 | 0.721 | -2.493 | 0.330 |
| G491 | -1.753 | 0.759 | -3.293 | -0.264 |
| G515 | -1.201 | 0.719 | -2.596 | 0.205 |
| G517 | -1.793 | 0.721 | -3.208 | -0.419 |
| G518 | -1.304 | 0.736 | -2.728 | 0.132 |
| G519 | -1.586 | 0.784 | -3.165 | -0.041 |
| G520 | -1.671 | 0.738 | -3.157 | -0.275 |
| G521 | -1.781 | 0.748 | -3.258 | -0.339 |
| G522 | -1.649 | 0.745 | -3.111 | -0.230 |
| G523 | -2.056 | 0.755 | -3.579 | -0.576 |
| G524 | -1.571 | 0.746 | -3.084 | -0.117 |
| G551 | -1.388 | 0.737 | -2.834 | 0.085 |
| G552 | -1.477 | 0.758 | -2.964 | 0.019 |
| G553 | -1.676 | 0.714 | -3.093 | -0.272 |
| G554 | -1.513 | 0.742 | -2.966 | -0.088 |
| G555 | -1.588 | 0.721 | -3.036 | -0.209 |
| G556 | -1.689 | 0.743 | -3.179 | -0.246 |
| G557 | -2.032 | 0.770 | -3.559 | -0.525 |
| G558 | -2.479 | 0.776 | -4.051 | -1.009 |
| G559 | -2.138 | 0.775 | -3.676 | -0.637 |
| G560 | -2.104 | 0.749 | -3.597 | -0.668 |
| G591 | -1.891 | 0.757 | -3.394 | -0.430 |
| G592 | -1.597 | 0.748 | -3.081 | -0.131 |
| G593 | -1.836 | 0.745 | -3.295 | -0.397 |
| G594 | -1.698 | 0.747 | -3.164 | -0.242 |
| G595 | -1.531 | 0.749 | -3.003 | -0.062 |
| G596 | -1.973 | 0.776 | -3.491 | -0.458 |
| G597 | -1.903 | 0.767 | -3.416 | -0.432 |
| G627 | -2.184 | 0.750 | -3.686 | -0.750 |
| G628 | -1.491 | 0.767 | -3.024 | -0.014 |
| G631 | -1.275 | 0.759 | -2.772 | 0.181 |
| G662 | -1.501 | 0.782 | -3.044 | 0.008 |
| G663 | -1.897 | 0.752 | -3.409 | -0.446 |
| G664 | -1.457 | 0.749 | -2.964 | 0.006 |
| G699 | -1.659 | 0.806 | -3.248 | -0.088 |
| G734 | -1.631 | 0.805 | -3.217 | -0.027 |
Show code
ran_n[["Season"]][[1]] %>%
kbl("html", caption = "Model estimates for random effect levels of 'Season' on the probability of occupancy for native species.", digits = 3) %>%
kable_styling(full_width = FALSE)
| Estimate | SD | 2.5% | 97.5% | |
|---|---|---|---|---|
| 2 | -1.803 | 0.531 | -2.874 | -0.805 |
| 4 | -1.714 | 0.513 | -2.736 | -0.737 |
| 5 | -1.352 | 0.595 | -2.547 | -0.184 |
MCMC Convergence and Diagnostics
Given that we implemented a Bayesian approach here the following section presents several model diagnostics used to ensure model convergence and suitability. The traceplots show how that the chains have mixed well and the number of iterations used in the model is sufficient in reaching convergence.
Show code
ubms::traceplot(natives_global_ubms, pars=c("beta_state", "beta_det"), ncol = 3)
Figure 9: Traceplots of the four chains run for 3000 iterations in the natives model. Plots show good evidence that chain-mixing and convergence has been achieved
Model Fit
Residuals
The residuals for the occupancy model can be calculated and plotted. The ubms package implements a method from Wright, Irvine, and Higgs (2019) to calculate residuals, where residuals are calculated independently for the observation and detection submodels. The diagnostic information from visualising these plots can be summarised as follows: if 95 % of the binned residual points fall within the grey shading then it is an indication the model fits the data well (Kellner 2021).
Show code
plot_residuals(natives_global_ubms, submodel="state")
Figure 10: Residual plot for the natives model, showing that binned residuals largely fall within the margins of error a model with good predictive power would have
Posterior Predictive Tests
We can predict the posterior distribution of the model and compare those predictions to the actual data. Here we show that the posterior prediction distribution matching the real data at a generalized level (Figure 11). However, given the slight overestimate in the median simulate posterior value to the true value we investigate the predictive power at a species level for both the power to predict occupancy at a given site (Figure 12A) and the number of times that the species is detected at a station when occupied (Figure 12B).
Show code
n_draws <- 250
nrows <- nrow(umf_natives@siteCovs)
species <- unique(umf_natives@siteCovs$Species)
sim_y <- posterior_predict(natives_global_ubms, "y", draws=n_draws, re.form = NULL)
prop0 <- apply(sim_y, 1, function(x) mean(x==0, na.rm=TRUE))
prop1 <- apply(sim_y, 1, function(x) mean(x==1, na.rm=TRUE))
actual_prop0 <- mean(getY(natives_global_ubms) == 0, na.rm=TRUE)
actual_prop1 <- mean(getY(natives_global_ubms) == 1, na.rm=TRUE)
#Compare 0
par(mfrow = c(1,2))
hist(prop0, col='gray', breaks = 15)
abline(v=actual_prop0, col='red', lwd=2)
#Compare 1
hist(prop1, col='gray', breaks = 15)
abline(v=actual_prop1, col='red', lwd=2)
Figure 11: Comparison of posterior predictions from the native species model to true mean values of absences and presences (0 and 1).
Show code
sim_y_props <- list()
# Species level
sim_y_props <-list()
for(x in species) {
sim_y_species_rows <- which(umf_natives@siteCovs$Species == x)
vec <- integer()
for(n in sim_y_species_rows) {
vec <- c(vec, ((n*104)-103):(n*104))
}
prop0 <- numeric()
prop1 <- numeric()
for(i in 1:n_draws) {
mat <- apply(matrix(sim_y[i,vec], ncol = 104, byrow = TRUE), 1, function(x) max(x, na.rm = TRUE))
mat[which(!is.finite(mat))] <- NA
prop1[i] <- mean(mat, na.rm = TRUE)
prop0[i] <- 1 - prop1[i]
}
p1 <- mean(apply(getY(natives_global_ubms)[sim_y_species_rows,], 1, function(x) {max(x,na.rm = T)}))
sim_y_props[[x]] <- data.frame(Species = x,
prop0 = prop0,
prop1 = prop1,
actual_prop0 = 1-p1,
actual_prop1 = p1)
}
bound_sim2 <- bind_rows(sim_y_props)
# sim_y_species <- sim_y_props %>%
# as.data.table() %>%
# melt.data.table(id.vars = c("Species", "Season"), measure.vars = cols_to_sel) %>%
# mutate(prop = substr(variable, 1, 5))
# group_by(Species, Season, prop) %>%
# summarise(Proportion = mean(value))
station_props <- bound_sim2 %>%
ggplot(aes(x = prop1)) +
geom_density(fill = "grey90") +
geom_vline(aes(xintercept = actual_prop1), colour = "red") +
xlim(0, 0.4) +
xlab("Proportion of stations\n where species were recorded")+
theme_bw() +
facet_wrap(~ Species, ncol = 1, scales = "free")
sim_y_props <-list()
for(x in species) {
sim_y_species_rows <- which(umf_natives@siteCovs$Species == x)
vec <- integer()
for(n in sim_y_species_rows) {
vec <- c(vec, ((n*104)-103):(n*104))
}
prop0 <- numeric()
prop1 <- numeric()
for(i in 1:n_draws) {
prop1[i] <- mean(apply(matrix(sim_y[i,vec], ncol = 104, byrow = TRUE), 1, function(x) max(x, na.rm = TRUE)))
prop0[i] <- 1 - prop1[i]
}
p1 <- mean(apply(getY(natives_global_ubms)[sim_y_species_rows,], 1, function(x) {max(x,na.rm = T)}))
sim_y_props[[x]] <- data.frame(Species = x,
prop0 = apply(sim_y[,vec], 1, function(x) mean(x==0, na.rm=TRUE)),
prop1 = apply(sim_y[,vec], 1, function(x) mean(x==1, na.rm=TRUE)),
actual_prop0 = mean(getY(natives_global_ubms)[sim_y_species_rows,] == 0, na.rm=TRUE),
actual_prop1 = mean(getY(natives_global_ubms)[sim_y_species_rows,] == 1, na.rm=TRUE))
}
bound_sim <- bind_rows(sim_y_props)
# sim_y_species <- sim_y_props %>%
# as.data.table() %>%
# melt.data.table(id.vars = c("Species", "Season"), measure.vars = cols_to_sel) %>%
# mutate(prop = substr(variable, 1, 5))
# group_by(Species, Season, prop) %>%
# summarise(Proportion = mean(value))
all_props <- bound_sim %>%
ggplot(aes(x = prop1)) +
geom_density(fill = "grey90") +
geom_vline(aes(xintercept = actual_prop1), colour = "red") +
xlim(0, 0.1) +
theme_bw() +
xlab("Proportion of days during deployment\n where species were recorded")+
facet_wrap(~ Species, ncol = 1, scales = "free")
cowplot::plot_grid(station_props, all_props, nrow = 1, labels = "AUTO")
Figure 12: Comparison of posterior predictions from the native species model to true mean values of presences, grouped for each species and station (A) as well as total counts of presences across the duration of camera trap deployment (max 104 days). The results indicate that predictive power is high for being able to predict whether or not a species is present at a site/station (A) and how many times that species is likely to be detected when present (B).
Methods Summary
In order to investigate the impact of fire and predation on native species occupancy we used a Bayesian approach to model species occupancy (Kellner 2021) while accounting for imperfect detection. Firstly, we generated predictions of occupancy for three invasive predators (foxes, cats and dogs) that were recorded in 2016/2017, 2019 and 2020. These site and season specific model predictions were then use as covariates within a native species occupancy model (albeit fox occupancy was replaced by fox density CITATION). We used model selection on model selection on frequentist mixed-effects models of predator and native species occupancies (Bates et al. 2015), with the top model (lowest AICc) from model selection (Tables 7, 2; Bartoń (2020)) then used in the Bayesian occupancy model, which accounted for imperfect detection due to the inclusion of a detection submodel/sub-process.
The top predator model included landscape fire severity (Department of Environment (2021); but see Figure 1), mean annual temperature, annual precipitation and precipitation seasonality (BIO01, BIO12 and BIO15) (Karger et al. 2017). Because stations/sites were deployed in nearby clusters of three we used this cluster of three (Group) as a random effect as inclusion of the station/site resulted in poor model convergence. Additionally, the season (2016/17, 2019 or 2020) was included as second random effect. Dynamic occupancy models were not used due to missing data (2018) and inconsistent measurements of occupancy at a given station (partly due to the 2019/2020 bushfires). All model covariates were fixed with a Species interaction effect so that the impacts of the covariates can be investigated on the species specific level.
The top native species model included landscape fire severity [Department of Environment (2021); but see Supplementary Information], cat occupancy, dog occupancy and fox density (CITATION), mean annual temperature and annual precipitation (BIO01 and BIO12) (Karger et al. 2017). Similarly to the predator model ‘group’ and ‘season’ were included as random effects and all covariates were fixed with a ‘species’ interaction term for the five native species modelled here: Long-footed potoroo, Long-nosed potoroo, Long-nosed bandicoot, Common brushtail possum and Goanna sp..
In both predator and native species models our detection submodel/subprocess included an interaction between the bait type (predator lure or native mammal lure) and days since deployment (1 to 52). Additionally, species was included as a fixed effect in the detection submodel as there was evidence to suggest either detectibility varied between species or species existed at different densities meaning the frequency of observations would differ even when the true latent state is the same (i.e. the species is present).
We investigated the model fit of our occupancy models through several diagnostic and posterior predictive tests: all suggesting convergence and mixing of chains, good model fits and accurate posterior predictive power (see Supplementary Material for extensive checks and relevent R code).
Finally, we compared whether predicted occupancy rates in fire affected areas may be lower for some species or at some stations due to the likelihood of fire refugia or landscape topography that may allow for greater persistence post-fire. We did this by determining the standardised mean change (cohen’s d) in occupancy at sites pre and post fire (Cohen 2013) against Terrain Ruggedness Index (TRI), which is a measure of the variability (ruggedness) of the elevation surrounding a given point (Riley, DeGloria, and Elliot 1999).
Results
Detection probabilities for predators and native mammals
Effect of bait type and days
Show code
predator_det_data <- list()
for(i in 1:3) {
predator_det_data[[i]] <- condition_effects_plot(model = predator_global_ubms,
submodel = "det",
covariate = "days",
facet = "baitMethod",
level = i,
return_data = TRUE) %>%
mutate(model = "Predator model")
}
predator_det_data <- bind_rows(predator_det_data)
native_det_data <- list()
for(i in 1:5) {
native_det_data[[i]] <- condition_effects_plot(model = natives_global_ubms,
submodel = "det",
covariate = "days",
facet = "baitMethod",
level = i,
return_data = TRUE) %>%
mutate(model = "Native species model")
}
native_det_data <- bind_rows(native_det_data)
det_data <- bind_rows(predator_det_data, native_det_data) %>%
mutate(`Bait lure` = case_when(baitMethod == "nonPred" ~ "Non predator",
baitMethod == "Pred" ~ "Predator")) %>%
arrange(desc(model)) %>%
unique() %>%
mutate(Species = case_when(Species == "LFP" ~ "Long-footed Potoroo",
Species == "LNP" ~ "Long-nosed Potoroo",
Species == "LNB" ~ "Long-nosed Bandicoot",
Species == "CBP" ~ "Common Brushtail Possum",
Species == "GOA" ~ "Goanna",
Species == "DOG" ~ "Dingo/Wild dog",
Species == "CAT" ~ "Feral cat",
Species == "FOX" ~ "Fox"))
det_plot <- ggplot(data=det_data, aes(x=days, y=Predicted, colour = `Bait lure`, fill = `Bait lure`)) +
geom_ribbon(aes(ymin=`2.5%`, ymax=`97.5%`), alpha=0.5) +
geom_line() +
labs(x = "Days", y = "Detection Probability (0 – 1)") +
# ylim(0, 0.15) +
scale_fill_manual(values = c(Predator = delwp_cols[["Navy"]], `Non predator` = delwp_cols[["Teal"]])) +
scale_colour_manual(values = c(Predator = delwp_cols[["Navy"]], `Non predator` = delwp_cols[["Teal"]])) +
facet_wrap( ~ factor(model, levels = c("Predator model", "Native species model"))+Species, scales = "free_y", ncol = 3) +
theme_bw() +
theme(panel.spacing = unit(25, "pt"))
day1_pred_pred <- predict(predator_global_ubms, submodel = "det", newdata = data.frame(days = 1, baitMethod = "Pred", Species = "CAT"))
day1_pred_nonpred <- predict(predator_global_ubms, submodel = "det", newdata = data.frame(days = 1, baitMethod = "nonPred", Species = "CAT"))
day1_nonpred_pred <- predict(natives_global_ubms, submodel = "det", newdata = data.frame(days = 1, baitMethod = "Pred", Species = "LFP", temp = median(umf_natives@obsCovs$temp)))
day1_nonpred_nonpred <- predict(natives_global_ubms, submodel = "det", newdata = data.frame(days = 1, baitMethod = "nonPred", Species = "LFP", temp = median(umf_natives@obsCovs$temp)))
shift_legend2(det_plot)
Figure 13: Estimated detection probabilities were dependent on the type of lure used (predator vs non-predator) as well as the number of days since the lure was set-up. Detection probability was highest initially after set-up in cases where the correct lure was used for the target species (e.g. predator lure for predators or non-predator lure for native species/mammals). The intercept of the detection probability varied between species.
Effect of bait and temperature
Show code
native_det_temp_data <- list()
for(i in 1:5) {
native_det_temp_data[[i]] <- condition_effects_plot(model = natives_global_ubms,
submodel = "det",
covariate = "temp",
facet = "baitMethod",
level = i,
return_data = TRUE) %>%
mutate(model = "Native species model")
}
native_det_data <- bind_rows(native_det_temp_data)
det_data_temp <- bind_rows(native_det_temp_data) %>%
filter(baitMethod == "nonPred") %>%
mutate(`Bait lure` = case_when(baitMethod == "nonPred" ~ "Non predator",
baitMethod == "Pred" ~ "Predator")) %>%
arrange(desc(model)) %>%
unique() %>%
mutate(Species = case_when(Species == "LFP" ~ "Long-footed Potoroo",
Species == "LNP" ~ "Long-nosed Potoroo",
Species == "LNB" ~ "Long-nosed Bandicoot",
Species == "CBP" ~ "Common Brushtail Possum",
Species == "GOA" ~ "Goanna",
Species == "DOG" ~ "Dingo/Wild dog",
Species == "CAT" ~ "Feral cat",
Species == "FOX" ~ "Fox"))
det_plot_temp <- ggplot(data=det_data_temp, aes(x=temp, y=Predicted, colour = `Bait lure`, fill = `Bait lure`)) +
geom_ribbon(aes(ymin=`2.5%`, ymax=`97.5%`), alpha=0.5) +
geom_line() +
labs(x = "Daily Maximum Temperature", y = "Detection Probability (0 – 1)") +
# ylim(0, 0.15) +
scale_fill_manual(values = c(`Non predator` = delwp_cols[["Teal"]])) +
scale_colour_manual(values = c(`Non predator` = delwp_cols[["Teal"]])) +
facet_wrap( ~ factor(model, levels = c("Predator model", "Native species model"))+Species, scales = "free_y", ncol = 3) +
theme_bw() +
theme(panel.spacing = unit(25, "pt"))
shift_legend2(det_plot_temp)
Figure 14: Estimated detection probabilities in the native species model were dependent on the nightly temperature. Detection probabilities varied according to nightly temperature. However the direction and magnitude of this effect differed between species. Goannas and Long-nosed potoroos showed some evidence that higher detection on warmer days/nights detection probabilities were higher. Long-nosed bandicoots showed stronger evidence that detection probabilities were higher on warmer days/nights. However, for Long-footed Potoroos this relationship was negative and for Common brushtail possums there was no discernable effect.
Show code
cumulative_nd <- data.frame(temp = mean(natives_global_ubms@data@obsCovs$temp, na.rm = TRUE),
days = seq(1, to = max(natives_global_ubms@data@obsCovs$days, na.rm = TRUE)),
baitMethod = "nonPred") %>%
left_join(data.frame(Species = unique(natives_global_ubms@data@siteCovs$Species)), by = character())
cumulative_det <- bind_cols(cumulative_nd, predict(natives_global_ubms,
newdata = cumulative_nd, submodel = "det")) %>%
group_by(Species) %>%
mutate(nd = 1-Predicted,
cum = cumprod(nd),
cum_det = 1-cum,
nd_2.5 = 1-(cumprod(1-`2.5%`)),
nd_97.5 = 1-(cumprod(1-`97.5%`)),
threshold_75 = cum_det > 0.75,
threshold_90 = cum_det > 0.9,
delta = cum_det - lag(cum_det)) %>%
mutate(Species = case_when(Species == "LFP" ~ "Long-footed Potoroo",
Species == "LNP" ~ "Long-nosed Potoroo",
Species == "LNB" ~ "Long-nosed Bandicoot",
Species == "CBP" ~ "Common Brushtail Possum",
Species == "GOA" ~ "Goanna",
Species == "DOG" ~ "Dingo/Wild dog",
Species == "CAT" ~ "Feral cat",
Species == "FOX" ~ "Fox"))
cumulative_det_threshold_90 <- cumulative_det %>%
filter(threshold_90) %>%
group_by(Species) %>%
slice(1)
cumulative_det_threshold_75 <- cumulative_det %>%
filter(threshold_75) %>%
group_by(Species) %>%
slice(1)
cumulative_det_delta1 <- cumulative_det %>%
filter(delta < 0.01) %>%
group_by(Species) %>%
slice(1)
ggplot(data = cumulative_det, aes(x = days, y = cum_det)) +
geom_ribbon(aes(ymin = nd_2.5, ymax = nd_97.5), fill = delwp_cols["Planning"], alpha = 0.5) +
geom_line() +
geom_vline(data = cumulative_det_threshold_75, aes(xintercept = days), linetype = 2) +
geom_vline(data = cumulative_det_threshold_90, aes(xintercept = days), linetype = 4, colour = "#a50f15") +
geom_label(data = cumulative_det_threshold_75, aes(x = days, y = 0.1, label = ">75% p"), size = 4.5) +
geom_label(data = cumulative_det_threshold_90, aes(x = days, y = 0.5, label = ">90% p"), size = 4.5, colour = "#a50f15") +
ylab("Detection probability (p)") +
xlab("Days since deployment") +
facet_wrap(~ Species, scales = "free_x") +
theme_bw() +
theme(panel.spacing = unit(25, "pt"))
Figure 15: Cumulative detection probabilities over the course of deployment varies between species. Dpending on target species, cameras may wish to be deployed for varying time periods. While all species have a 75% probability of detection within 50 days only Long-footed Potoroos, Common Brushtail Possums and Long-nosed Bandicoots achieve 90% cumulative detection probabilities within the 50 days. These cumulative detection probabilities are based on fixed mean maximum temperatures of 18.8°C
Show code
# ggplot(data = cumulative_det, aes(x = days, y = delta)) +
# geom_line() +
# ylab("Detection probability (p) change") +
# xlab("Days since deployment") +
# facet_wrap(~ Species, scales = "free_x") +
# theme_bw() +
# theme(panel.spacing = unit(25, "pt"))
Detection probabilities varied for native mammal species and for predators. Both the type of bait/lure used and the number of days since the lure was set up predicted the detection probability.
For the predator model (the model where the response variable was the detection of cats, dogs or foxes) we found (i) an effect of lure type and (ii) an interaction with days. Firstly, predator lures had a higher detection probability than non-predator lures when baits were initially deployed (Day 1 cat detection \(p_{predator\ lure}\) = 0.043 [95% CI = 0.032, 0.057], Day 1 cat detection \(p_{native\ lure}\) = 0.018 [95% CI = 0.014, 0.022]). However, the detection probability for the predator lure decreased over time with the baits having a similar detection probability after 25-30 days (Figure 13).
For the native species model (the model where the response variable was the detection of bandicoots, potoroos, possums and goannas) we found (i) an effect of lure type, (ii) a decrease in detection over time and (iii) no discernible interaction between lure type and days on detection probability. Specifically, non predator lures were more effective at detective native species (Day 1 long-footed potoroo detection \(p_{native\ lure}\) = NA [95% CI = NA, NA], Day 1 long-footed potoroo detection \(p_{predator\ lure}\) = NA [95% CI = NA, NA]), with the probability of detection reducing over time for both lures (Figure 13). In both predator and native species models the species detected had a fixed effect (Table 4, 9) in changing the detection probability; however no interactions with species were modelled thus ‘species’ only affects the intercept and not the slope of decreased detectibility over time, or differences between lure types.
Given the variability in detection probabilities we generated cumulative probability of detection values across the total observation period (maximum of 52 days). The results of this exercise shown in Figure 15.
Additionally, we investigated the role that maximum daily temperature has on detectibility. It is predicted that detectibility for reptiles (goannas) will increase with higher maximum daily temperatures. To assess this relationship we obtained estimated daily maximum temperatures from a 0.5 degree x 0.5 degree raster grid (Oceanic and Labratory) 2021). See Table 12) for the conditional effects of temperature and days for different species and bait methods and Figure 14 for the graphical representation of detection probabilities across temperature values.
Show code
preds <- c("temp")
table_conditional_det <- conditional_effects_ubms(natives_global_ubms,
submodel = "det",
categoricalParam = "Species",
continuousParam = "temp", all = FALSE)
table_condition_bait <- conditional_effects_ubms(natives_global_ubms,
submodel = "det",
categoricalParam = "baitMethod",
continuousParam = "days", all = FALSE)
table_conditional_det %>%
dplyr::select(variable, Category = Species, mean, sd, `2.5%`, `97.5%`) %>%
bind_rows(table_condition_bait %>%
dplyr::select(variable, Category = baitMethod, mean, sd, `2.5%`, `97.5%`)) %>%
kbl("html", caption = "Conditional effects from model estimates showing the effects of various covariates on the probability of detection (logit-scale) for native species. Estimates with a 95% CI not overlapping zero are shown in bold.", digits = 3) %>%
kable_styling(full_width = FALSE) %>%
row_spec(which(table_conditional_det$Sig), bold = TRUE) %>%
column_spec(1, bold = TRUE) %>%
# row_spec(which(table_conditional$Species == "ALL"), background = "#f0f0f0") %>%
collapse_rows(1:2)
| variable | Category | mean | sd | 2.5% | 97.5% |
|---|---|---|---|---|---|
| temp | CBP | -0.003 | 0.006 | -0.014 | 0.008 |
| GOA | 0.014 | 0.010 | -0.005 | 0.033 | |
| LFP | -0.011 | 0.003 | -0.017 | -0.004 | |
| LNB | 0.034 | 0.009 | 0.016 | 0.051 | |
| LNP | 0.024 | 0.017 | -0.009 | 0.058 | |
| days | nonPred | -0.015 | 0.001 | -0.018 | -0.012 |
| Pred | -0.010 | 0.004 | -0.018 | -0.002 |
Predictions of effects on native mammal occupancy
Given that we modelled species as a categorical variable and fitted it with various interaction terms (for fire, predators, and bioclimatic variables), it can be difficult to interpret the marginal effects of fire on each species from the model estimates. Instead we can obtain the conditional effects by combining the estimate of the base effect (e.g. fire severity) with the estimate of the interaction term (e.g. species X fire severity). This is done by extracting the posterior distributions from the fitted model. Table 13 shows the conditional effects for the native species model while Table 14 shows the conditional effects for the precoursor predator model.
Show code
preds <- c("scaled_fire_severity", "FoxDensity", "CAT", "DOG", "scale(BIO01)", "scale(BIO12)")
table_conditional <- lapply(preds, function(x) {
conditional_effects_ubms(natives_global_ubms,
submodel = "state",
categoricalParam = "Species",
continuousParam = x, all = FALSE)
}) %>% bind_rows()
table_conditional %>%
dplyr::select(variable, Species, mean, sd, `2.5%`, `97.5%`) %>%
kbl("html", caption = "Conditional effects from model estimates showing the effects of various covariates on the probability of occupancy (logit-scale) for native species. Estimates with a 95% CI not overlapping zero are shown in bold.", digits = 3) %>%
kable_styling(full_width = FALSE) %>%
row_spec(which(table_conditional$Sig), bold = TRUE) %>%
column_spec(1, bold = TRUE) %>%
# row_spec(which(table_conditional$Species == "ALL"), background = "#f0f0f0") %>%
collapse_rows(1:2)
| variable | Species | mean | sd | 2.5% | 97.5% |
|---|---|---|---|---|---|
| scaled_fire_severity | CBP | -0.154 | 0.500 | -1.132 | 0.801 |
| GOA | -5.487 | 1.241 | -8.135 | -3.255 | |
| LFP | -1.120 | 0.554 | -2.202 | -0.065 | |
| LNB | -2.045 | 0.705 | -3.473 | -0.702 | |
| LNP | -2.919 | 1.401 | -6.012 | -0.499 | |
| FoxDensity | CBP | 0.506 | 0.297 | -0.087 | 1.087 |
| GOA | 0.663 | 0.376 | -0.062 | 1.423 | |
| LFP | -8.148 | 0.900 | -9.925 | -6.448 | |
| LNB | -0.359 | 0.410 | -1.194 | 0.414 | |
| LNP | -0.069 | 0.737 | -1.644 | 1.227 | |
| CAT | CBP | 0.689 | 1.179 | -1.580 | 3.047 |
| GOA | 0.098 | 1.737 | -3.335 | 3.492 | |
| LFP | 2.646 | 1.724 | -0.589 | 6.126 | |
| LNB | 1.602 | 1.712 | -1.680 | 5.074 | |
| LNP | 0.112 | 1.942 | -3.826 | 3.835 | |
| DOG | CBP | -0.021 | 1.467 | -2.951 | 2.936 |
| GOA | 0.161 | 2.222 | -4.246 | 4.682 | |
| LFP | -1.372 | 2.236 | -5.932 | 2.869 | |
| LNB | 0.840 | 2.172 | -3.228 | 5.332 | |
| LNP | -0.205 | 2.274 | -4.781 | 4.182 | |
| scale(BIO01) | CBP | 0.516 | 0.134 | 0.253 | 0.789 |
| GOA | 1.255 | 0.216 | 0.838 | 1.677 | |
| LFP | -0.862 | 0.160 | -1.177 | -0.552 | |
| LNB | 0.439 | 0.171 | 0.116 | 0.773 | |
| LNP | 1.082 | 0.351 | 0.445 | 1.797 | |
| scale(BIO12) | CBP | -0.410 | 0.101 | -0.607 | -0.213 |
| GOA | 0.132 | 0.131 | -0.119 | 0.393 | |
| LFP | -0.019 | 0.109 | -0.233 | 0.192 | |
| LNB | 0.197 | 0.123 | -0.042 | 0.441 | |
| LNP | 0.477 | 0.216 | 0.055 | 0.906 |
Show code
preds <- c("scaled_fire_severity", "scale(BIO01)", "scale(BIO12)", "scale(BIO15)")
pred_table_conditional <- lapply(preds, function(x) {
conditional_effects_ubms(predator_global_ubms,
submodel = "state",
categoricalParam = "Species",
continuousParam = x, all = FALSE)
}) %>% bind_rows()
pred_table_conditional %>%
dplyr::select(variable, Species, mean, sd, `2.5%`, `97.5%`) %>%
kbl("html", caption = "Conditional effects from model estimates showing the effects of various covariates on the probability of occupancy (logit-scale) for Predators. Estimates with a 95% CI not overlapping zero are shown in bold.", digits = 3) %>%
kable_styling(full_width = FALSE) %>%
collapse_rows(1:2) %>%
row_spec(which(pred_table_conditional$Sig), bold = TRUE) %>%
column_spec(1, bold = TRUE)
| variable | Species | mean | sd | 2.5% | 97.5% |
|---|---|---|---|---|---|
| scaled_fire_severity | CAT | -1.136 | 0.632 | -2.546 | -0.026 |
| DOG | -1.741 | 1.064 | -4.029 | 0.098 | |
| FOX | -0.153 | 0.826 | -1.836 | 1.326 | |
| scale(BIO01) | CAT | -0.309 | 0.119 | -0.550 | -0.082 |
| DOG | -0.156 | 0.193 | -0.517 | 0.233 | |
| FOX | -0.555 | 0.172 | -0.892 | -0.214 | |
| scale(BIO12) | CAT | 0.024 | 0.131 | -0.232 | 0.272 |
| DOG | -0.376 | 0.225 | -0.814 | 0.077 | |
| FOX | -0.829 | 0.215 | -1.252 | -0.417 | |
| scale(BIO15) | CAT | -0.063 | 0.124 | -0.308 | 0.182 |
| DOG | 0.093 | 0.220 | -0.337 | 0.525 | |
| FOX | 0.729 | 0.212 | 0.322 | 1.146 |
Show code
# row_spec(which(pred_table_conditional$Species == "ALL"), background = "#f0f0f0")
Effects of Fire
Show code
stations_used <- natives_site_data %>%
filter(Station %in% sites_visited_pre_post) %>%
pull(Station) %>%
unique()
season_4_sites <- natives_site_data %>%
filter(Season == 4) %>%
pull(Station) %>%
unique()
season_5_sites <- natives_site_data %>%
filter(Season == 5) %>%
pull(Station) %>%
unique()
new_data_pre_post_fire <- natives_site_data %>%
# filter(Season == 5) %>%
dplyr::select(Station,
Group,
Season,
Species,
scaled_fire_severity,
contains("BIO"),
CAT,
DOG,
FOX,
FoxDensity,
tri,
GroupSpecies,
SeasonSpecies) %>%
distinct() %>%
mutate(Time = case_when(Season < 5 ~ "Pre-fire",
Season == 5 ~ "Post-fire")) %>%
arrange(Time)
fire_preds_ag <- bind_cols(new_data_pre_post_fire,
predict(natives_global_ubms,
submodel = "state",
newdata = new_data_pre_post_fire))
new_data_4_5 <- new_data_pre_post_fire %>%
filter(Season > 3 & Station %in% season_4_sites & Station %in% season_5_sites)
post <- ubms:::sim_lp(natives_global_ubms,
submodel = "state",
transform = TRUE,
newdata = new_data_pre_post_fire,
samples = 1:100,
re.form = NULL) %>%
t()
fire_preds <- bind_cols(new_data_pre_post_fire %>%
dplyr::select(Station, Season, Species, scaled_fire_severity, tri),
post) %>%
filter(Station %in% sites_visited_pre_post) %>%
mutate(Time = case_when(Season < 5 ~ "Pre-fire",
Season == 5 ~ "Post-fire")) %>%
group_by(Station, Time, Species) %>%
filter(Season == max(as.numeric(Season), na.rm = TRUE)) %>%
ungroup() %>%
as.data.table() %>%
melt(id.vars = c("Station", "Time", "Season", "Species", "scaled_fire_severity", "tri"))
#### Predator predictions ####
new_data_pre_post_fire_pred <- predator_site_data %>%
# filter(Season == 5) %>%
dplyr::select(Station,
Group,
Season,
Species,
scaled_fire_severity,
contains("BIO"),
GroupSpecies,
SeasonSpecies) %>%
distinct() %>%
mutate(Time = case_when(Season < 5 ~ "Pre-fire",
Season == 5 ~ "Post-fire")) %>%
arrange(Time)
post_predators <- ubms:::sim_lp(predator_global_ubms,
submodel = "state",
transform = TRUE,
newdata = new_data_pre_post_fire_pred,
samples = 1:100,
re.form = NULL) %>%
t()
fire_preds_pred <- bind_cols(new_data_pre_post_fire_pred %>%
dplyr::select(Station, Season, Species, scaled_fire_severity),
post_predators) %>%
filter(Station %in% sites_visited_pre_post) %>%
mutate(Time = case_when(Season < 5 ~ "Pre-fire",
Season == 5 ~ "Post-fire")) %>%
group_by(Station, Time, Species) %>%
filter(Season == max(as.numeric(Season), na.rm = TRUE)) %>%
ungroup() %>%
as.data.table() %>%
melt(id.vars = c("Station", "Time", "Season", "Species", "scaled_fire_severity"))
fire_preds_pre_post_comparison_pred <- fire_preds_pred %>%
as.data.table() %>%
dcast(Station + Species ~ Time, value.var = "value", fun.aggregate = mean) %>%
as.data.frame() %>%
# rename(`2019` = `4`, `2020` = `5`) %>%
mutate(diff = `Post-fire`-`Pre-fire`,
perc_red = (1 - `Post-fire`/`Pre-fire`)) %>%
group_by(Species) %>%
summarise(`Mean Occupancy Change` = round(mean(diff), 3),
`NSites Decreased` = sum(diff < 0),
`Mean Percentage Reduction` = round(mean(perc_red)*100, 2),
`NSites Decreased 20%` = sum(perc_red > 0.2),
`total sites` = n())
#### End predator ####
station_plot <- fire_preds %>%
# mutate(Season = as.character(2015L + Season)) %>%
ggplot(aes(x = value, y = Station)) +
ggridges::geom_density_ridges(aes(fill = Time), alpha = 0.7) +
scale_fill_manual(values = c("Pre-fire" = delwp_cols[["Teal"]], "Post-fire" = delwp_cols[["FFR"]])) +
theme_minimal() +
theme(legend.position = "top") +
facet_wrap(~ Species, ncol = 5, scales = "free_x")
station_sum_plot <- fire_preds %>%
# mutate(Season = as.character(2015L + Season),
mutate(Station = "") %>%
ggplot(aes(x = value, y = Station)) +
scale_y_discrete(limits = "", expand = c(0,0)) +
theme_minimal() +
ylab("") +
ggridges::geom_density_ridges(aes(fill = Time), alpha = 0.7) +
scale_fill_manual(values = c("Pre-fire" = delwp_cols[["Teal"]], "Post-fire" = delwp_cols[["FFR"]]), guide = FALSE) +
theme(plot.margin = margin(l = 50)) +
facet_wrap(~ Species, ncol = 5, scales = "free")
fire_preds_pre_post_comparison <- fire_preds %>%
as.data.table() %>%
dcast(Station + Species ~ Time, value.var = "value", fun.aggregate = mean) %>%
as.data.frame() %>%
# rename(`2019` = `4`, `2020` = `5`) %>%
mutate(diff = `Post-fire`-`Pre-fire`,
perc_red = (1 - `Post-fire`/`Pre-fire`)) %>%
group_by(Species) %>%
summarise(`Mean Occupancy Change` = round(mean(diff), 3),
`NSites Decreased` = sum(diff < 0),
`Mean Percentage Reduction` = round(mean(perc_red)*100, 2),
`NSites Decreased 20%` = sum(perc_red > 0.2),
`total sites` = n())
fire_prediction_plot_multiseason <- fire_preds_ag %>%
mutate(Season = 2015L + Season) %>%
filter(Station %in% sites_visited_pre_post) %>%
# filter(!Species %in% c("CRP", "SBB")) %>%
ggplot(aes(x = Season, y = Predicted, fill = scaled_fire_severity)) +
scale_x_continuous(breaks = c(2017:2020)) +
scale_fill_distiller(palette = "YlOrRd", direction = 1,
guide = guide_legend(title = "Landscape fire severity (scaled)",
override.aes = list(size = 3))) +
geom_quasirandom(shape = 21, size = 1) +
facet_wrap(~ Species, ncol = 2, scales = "free_y") +
theme_bw()
Occupancy rates generally decreased after the 2019/2020 bushfires, and this effect was extenuated based on the fire severity within the surrounding landscape. However, the responses were species and site specific (Figure 16). If we compare stations that had pre-fire surveys (in 2016, 2017 or 2019) and a post-fire survey (2020) then we that Long-Nosed Potoroo, Long-nosed Bandicoot and Goanna predicted occupancy rates decreased the most (% reduction in probability of occupancy), with the majority of sites having mean decreases in occupancy. However, Long-footed Potoroos appeared to have felt a lesser reduction in occupancy after the fires and Common Brushtail Possums show no (or even a slight positive) change in occupancy when comparing pre and post fire data. The following data summarises this difference for the five native species. Note that predicted occupancy rates for pre-fire surveys are taken from the most recent pre-fire survey (e.g. if surveys were conducted in 2017 and 2019 only the 2019 data is used).
Show code
fire_preds_pre_post_comparison %>%
bind_rows(fire_preds_pre_post_comparison_pred) %>%
kbl(format = "html", caption = "Change in native species occupancy for sites that had both pre-fire and post-fire surveys.") %>%
kable_styling(full_width = FALSE)
| Species | Mean Occupancy Change | NSites Decreased | Mean Percentage Reduction | NSites Decreased 20% | total sites |
|---|---|---|---|---|---|
| CBP | 0.049 | 0 | -24.37 | 0 | 236 |
| GOA | -0.180 | 209 | 67.43 | 203 | 236 |
| LFP | -0.034 | 150 | 2.79 | 66 | 236 |
| LNB | -0.097 | 200 | 35.31 | 173 | 236 |
| LNP | -0.017 | 203 | 47.84 | 198 | 236 |
| CAT | -0.074 | 194 | 16.48 | 110 | 236 |
| DOG | -0.038 | 201 | 32.79 | 179 | 236 |
| FOX | 0.011 | 0 | -18.61 | 0 | 236 |
The site and specific predictions can be obtained via taking random draws of the posterior distribution for our model. Figure 16 shows the site-specific changes in predicted occupancy for the five native species.
Figure 16: Station and species specific predictions for the change in occupancy rates from pre-fire to post-fire. Aggregate distributions across all stations are shown at the bottom. Note that the aggregate values may not reflect the complexity of fire severity across the landscape/at various sites.
The effect of fire severity on occupancy can be visualised with the following plot. The effects of landscape fire severity show the largest effects for Long-Nosed Potoroo, Long-nosed Bandicoot and Goanna. More variable effects appear for Common Brushtail Possums and Long-footed Potoroos (Figure 17).
Show code
shift_legend2(fire_prediction_plot_multiseason)
Figure 17: Change in site occupancy for sites that recorded pre and post fire surveys. Plots are grouped by species and coloured according to the scaled landscape fre severity score. For surveys in 2016, 2017 and 2019 fire severity is zero as only the fire severity of the 2019/2020 fires were taken into account.
Show code
native_fire <- condition_effects_plot(model = natives_global_ubms,
submodel = "state",
covariate = "scaled_fire_severity",
facet = "Species",
ncol = 1, scales = "free_y", return_data = TRUE) %>%
mutate(Species = case_when(Species == "LFP" ~ "Long-footed Potoroo",
Species == "LNP" ~ "Long-nosed Potoroo",
Species == "LNB" ~ "Long-nosed Bandicoot",
Species == "CBP" ~ "Common Brushtail Possum",
Species == "GOA" ~ "Goanna",
Species == "DOG" ~ "Dingo/Wild dog",
Species == "CAT" ~ "Feral cat",
Species == "FOX" ~ "Fox")) %>%
ggplot(data=., aes_string(x="scaled_fire_severity", y="Predicted")) +
geom_ribbon(aes(ymin=`2.5%`, ymax=`97.5%`), alpha=0.3) +
geom_line() +
labs(x = "Scaled fire severity", y = "Occupancy (0 – 1)") +
facet_wrap(facets = "Species", ncol = 1, scales = "free_y") +
theme_bw() +
theme(panel.spacing = unit(25, "pt"))
pred_fire <- condition_effects_plot(model = predator_global_ubms,
submodel = "state",
covariate = "scaled_fire_severity",
facet = "Species",
ncol = 1, scales = "free_y") +
xlab("Scaled fire severity") #+
# theme(plot.margin = margin(5,5,348,5))
cowplot::plot_grid(native_fire, pred_fire, nrow = 1, labels = "AUTO")
Figure 18: Marginal effects of fire severity on native species (A) and predator species (B) occupancy rates. These are thentransformed values on the occupancy scale. Mean estimates are shown with the black line and 95% CI are shown with the grey bands.
Show code
log_plot <- plot_predicted_marginal(model = natives_global_ubms,
model_stack = umf_natives,
param = "scaled_logging",
groupings = c("Species", "Group", "Season"),
other_covs = c("scaled_fire_severity", "scaled_road_density", "BIO01", "BIO04", "BIO12", "BIO15", "CAT", "FoxScore", "DOG"), breaks = 10)
cat_plot <- plot_predicted_marginal(model = natives_global_ubms,
model_stack = umf_natives,
param = "CAT",
groupings = c("Species", "Group", "Season"),
other_covs = c("scaled_fire_severity", "scaled_road_density", "BIO01", "BIO04", "BIO12", "BIO15", "scaled_logging", "FoxScore", "DOG"), breaks = 10)
bio01_plot <- plot_predicted_marginal(model = natives_global_ubms,
model_stack = umf_natives,
param = "BIO01",
groupings = c("Species", "Group", "Season"),
other_covs = c("scaled_fire_severity", "scaled_road_density", "CAT", "BIO04", "BIO12", "BIO15", "scaled_logging", "FoxScore", "DOG"), breaks = 10)
bio04_plot <- plot_predicted_marginal(model = natives_global_ubms,
model_stack = umf_natives,
param = "BIO04",
groupings = c("Species", "Group", "Season"),
other_covs = c("scaled_fire_severity", "scaled_road_density", "CAT", "BIO01", "BIO12", "BIO15", "scaled_logging", "FoxScore", "DOG"), breaks = 10)
bio12_plot <- plot_predicted_marginal(model = natives_global_ubms,
model_stack = umf_natives,
param = "BIO12",
groupings = c("Species", "Group", "Season"),
other_covs = c("scaled_fire_severity", "scaled_road_density", "CAT", "BIO01", "BIO04", "BIO15", "scaled_logging", "FoxScore", "DOG"), breaks = 10)
bio15_plot <- plot_predicted_marginal(model = natives_global_ubms,
model_stack = umf_natives,
param = "BIO15",
groupings = c("Species", "Group", "Season"),
other_covs = c("scaled_fire_severity", "scaled_road_density", "CAT", "BIO01", "BIO04", "BIO12", "scaled_logging", "FoxScore", "DOG"), breaks = 10)
fox_plot <- plot_predicted_marginal(model = natives_global_ubms,
model_stack = umf_natives,
param = "FoxScore",
groupings = c("Species", "Group", "Season"),
other_covs = c("scaled_fire_severity", "scaled_road_density", "CAT", "BIO01", "BIO04", "BIO12", "scaled_logging", "BIO15", "DOG"), breaks = 10)
dog_plot <- plot_predicted_marginal(model = natives_global_ubms,
model_stack = umf_natives,
param = "DOG",
groupings = c("Species", "Group", "Season"),
other_covs = c("scaled_fire_severity", "scaled_road_density", "FoxScore", "BIO01", "BIO04", "BIO12", "scaled_logging", "BIO15", "CAT"), breaks = 10)
# Effect of Fire Severity on Occupancy
cowplot::plot_grid(log_plot, cat_plot,
fox_plot, dog_plot,
bio01_plot, bio04_plot,
bio12_plot, bio15_plot, nrow = 4, labels = "AUTO")
Effects of Predation
The relationship between predation and native species occupancy varied; but apart from two relationships (LFP - fox density and LFP - cat occupancy) there were no relationships where the 95 % CI intervals did not overlap zero (Table 13. Specifically, there was a sharp decline in LFP occupancy when fox density was great than zero; with predicted long-footed potoroo occupancy estimated at zero for fox densities above 0.5 (Figure 19). Conversely there was a weaker positive relationship between long-footed potoroo occupancy and cat occupancy; indicating cats were likely to occupy sites that were also occupied by long-footed potoroos (Figure 19). However this relationship is marred by both cats and long-footed potoroos tending to occupy areas with lower value of BIO01 (mean annual temperature) (Tables 13, 14). Suggesting possible climatic or environmental confounding factors for this association.
Show code
foxplot <- condition_effects_plot(model = natives_global_ubms,
submodel = "state",
covariate = "FoxDensity",
facet = "Species",
ncol = 1, scales = "free_y", return_data = TRUE) %>%
mutate(Species = case_when(Species == "LFP" ~ "Long-footed Potoroo",
Species == "LNP" ~ "Long-nosed Potoroo",
Species == "LNB" ~ "Long-nosed Bandicoot",
Species == "CBP" ~ "Common Brushtail Possum",
Species == "GOA" ~ "Goanna",
Species == "DOG" ~ "Dingo",
Species == "CAT" ~ "Feral cat",
Species == "FOX" ~ "Fox")) %>%
ggplot(data=., aes_string(x="FoxDensity", y="Predicted")) +
geom_ribbon(aes(ymin=`2.5%`, ymax=`97.5%`), alpha=0.3) +
geom_line() +
labs(x = "Fox density", y = "Occupancy (0 – 1)") +
facet_wrap(facets = "Species", ncol = 1, scales = "free_y") +
theme_bw() +
theme(panel.spacing = unit(25, "pt"))
catplot <- condition_effects_plot(model = natives_global_ubms,
submodel = "state",
covariate = "CAT",
facet = "Species",
ncol = 1, scales = "free_y", return_data = TRUE) %>%
mutate(Species = case_when(Species == "LFP" ~ "Long-footed Potoroo",
Species == "LNP" ~ "Long-nosed Potoroo",
Species == "LNB" ~ "Long-nosed Bandicoot",
Species == "CBP" ~ "Common Brushtail Possum",
Species == "GOA" ~ "Goanna",
Species == "DOG" ~ "Dingo",
Species == "CAT" ~ "Feral cat",
Species == "FOX" ~ "Fox")) %>%
ggplot(data=., aes_string(x="CAT", y="Predicted")) +
geom_ribbon(aes(ymin=`2.5%`, ymax=`97.5%`), alpha=0.3) +
geom_line() +
labs(x = "Feral cat occupancy", y = "Occupancy (0 – 1)") +
facet_wrap(facets = "Species", ncol = 1, scales = "free_y") +
theme_bw() +
theme(panel.spacing = unit(25, "pt"))
dogplot <- condition_effects_plot(model = natives_global_ubms,
submodel = "state",
covariate = "DOG",
facet = "Species",
ncol = 1, scales = "free_y", return_data = TRUE) %>%
mutate(Species = case_when(Species == "LFP" ~ "Long-footed Potoroo",
Species == "LNP" ~ "Long-nosed Potoroo",
Species == "LNB" ~ "Long-nosed Bandicoot",
Species == "CBP" ~ "Common Brushtail Possum",
Species == "GOA" ~ "Goanna",
Species == "DOG" ~ "Dingo",
Species == "CAT" ~ "Feral cat",
Species == "FOX" ~ "Fox")) %>%
ggplot(data=., aes_string(x="DOG", y="Predicted")) +
geom_ribbon(aes(ymin=`2.5%`, ymax=`97.5%`), alpha=0.3) +
geom_line() +
labs(x = "Dingo occupancy", y = "Occupancy (0 – 1)") +
facet_wrap(facets = "Species", ncol = 1, scales = "free_y") +
theme_bw() +
theme(panel.spacing = unit(25, "pt"))
cowplot::plot_grid(foxplot, catplot, dogplot, nrow = 1, labels = "AUTO")
Figure 19: Marginal effects of total predation on native species occupancy rates (A = Fox, B = Cat and C = Dog). These are the graphical and transformed values from the table above. Mean estimates are shown with the black line and 95% CI are shown with the grey bands.
Reasons for lesser impacts from fire for certain species
There are several reasons why certain species may have occupancy less decreased by fire. One is the potential for fire refugia to exist and be more utilised. Because fire severity was calculated on a landscape level, surrounding fire refugia to a site would presumably lead to a lower fire severity score. However, depending on how far fire refugia is from a station the effect on the fire severity score may not be discernable. Additionally, fire severity measurements (Department of Environment 2021) may not be accurate i detecting very small areas where fire was less severe.
We use another potential variable that may infer likelihood of fire refugia to see whether species less effected by fire severity (LFP) could have used them. Terrain Ruggedness Index (TRI) is a measure of the variability (ruggedness) of the elevation surrounding a given point (Riley, DeGloria, and Elliot 1999). If TRI increased the abundance of fire refugia we would expect to see locations where occupancy persisted post-fire having hire TRI than in locations where occupancy decreased. However, no such trend for LFP can be seen for this relationship (Figure 20). Goana sp. may show a slight positive relationship however, this species was severely fire affected with detection at only two sites in 2020.
Show code
d <- fire_preds %>% filter(Time == "Post-fire") %>% dplyr::select(Station, scaled_fire_severity)
dq <- quantile(x = d$scaled_fire_severity, 0.5)
fire_preds_tri <- fire_preds %>%
group_by(Station, Time, Species) %>%
summarise(mean = mean(value),
sd = sd(value),
TRI = unique(tri)) %>%
ungroup() %>%
as.data.table() %>%
dcast(Station + Species + TRI ~ Time, value.var = c("mean", "sd")) %>%
as.data.frame()
fire_preds_fs <- fire_preds %>%
filter(Time == "Post-fire") %>%
dplyr::select(Station, scaled_fire_severity) %>%
distinct()
fire_preds_tri <- fire_preds_tri %>%
left_join(fire_preds_fs) %>%
filter(scaled_fire_severity > dq)
for(i in 1:nrow(fire_preds_tri)) {
es <- compute.es::mes(m.1 = fire_preds_tri[i, "mean_Post-fire"] ,
m.2 = fire_preds_tri[i, "mean_Pre-fire"],
sd.1 = fire_preds_tri[i, "sd_Post-fire"],
sd.2 = fire_preds_tri[i, "sd_Pre-fire"],
n.1 = 100,
n.2 = 100, verbose = FALSE)
fire_preds_tri[i, "cohen's d"] <- es$d
fire_preds_tri[i, "cohen's d variance"] <- es$var.d
fire_preds_tri[i, "cohen's d LCI"] <- es$l.d
fire_preds_tri[i, "cohen's d UCI"] <- es$u.d
}
fire_preds_tri %>%
ggplot(aes(x = `cohen's d`, y = TRI)) +
geom_segment(aes(x = `cohen's d LCI`, xend = `cohen's d UCI`, yend = TRI)) +
geom_point(shape = 21, fill = delwp_cols[["Environment"]], size = 1.5) +
geom_smooth(colour = delwp_cols[["Navy"]], method = "gam") +
facet_wrap(~Species, ncol = 3, scales = "fixed") +
ylim(0, 20) +
xlab("Standardised change in mean \noccupancy before-after fires (cohen's d)") +
theme_bw()
Figure 20: Difference between mean occupancy rates at stations before and after fires. Values are standardised using cohen’s d, which accounts for variable standard deviations between pre and post fire occupancy predictions. Change in occupancy (cohen’s d) is plotted against TRI. No species showed trends that would suggest TRI would allow for increased or maintained occupancy for stations with high severity fires. Station’s (points) are filtered so that only those with a fire severity in the top 50 % worst burned sites are shown
Show code
mean_prepost <- fire_preds_ag %>%
group_by(Species, Time) %>%
summarise(median = quantile(Predicted, 0.5),
LCI95 = quantile(Predicted, 0.025),
UCI95 = quantile(Predicted, 0.975)) %>%
ungroup() %>%
as.data.table() %>%
dcast(Species ~ Time, value.var = c("median", "LCI95", "UCI95")) %>%
mutate(diff_av = round(`median_Post-fire` - `median_Pre-fire`,2),
diff_LCI = round(`LCI95_Post-fire` - `LCI95_Pre-fire`, 2),
diff_UCI = round(`UCI95_Post-fire` - `UCI95_Pre-fire`, 2),
pre_fire = paste0(round(`median_Pre-fire`, 2), " [", round(`LCI95_Pre-fire`, 2), ", ", round(`UCI95_Pre-fire`, 2), "]"),
post_fire = paste0(round(`median_Post-fire`, 2), " [", round(`LCI95_Post-fire`, 2), ", ", round(`UCI95_Post-fire`, 2), "]")) %>%
dplyr::select(Species, pre_fire, post_fire, diff_av)
d <- fire_preds_ag %>% filter(Time == "Post-fire") %>% dplyr::select(Station, scaled_fire_severity)
dq <- quantile(x = d$scaled_fire_severity, 0.5)
mean_sev <- fire_preds_ag %>%
filter(Time == "Post-fire") %>%
mutate(fire_class = case_when(scaled_fire_severity >= dq ~ "High",
scaled_fire_severity < dq ~ "Low")) %>%
group_by(Species, fire_class) %>%
summarise(median = quantile(Predicted, 0.5),
LCI95 = quantile(Predicted, 0.025),
UCI95 = quantile(Predicted, 0.975)) %>%
ungroup() %>%
as.data.table() %>%
dcast(Species ~ fire_class, value.var = c("median", "LCI95", "UCI95")) %>%
mutate(diff_av = round(`median_High` - `median_Low`,2),
High = paste0(round(`median_High`, 2), " [", round(`LCI95_High`, 2), ", ", round(`UCI95_High`, 2), "]"),
Low = paste0(round(`median_Low`, 2), " [", round(`LCI95_Low`, 2), ", ", round(`UCI95_Low`, 2), "]")) %>%
dplyr::select(Species, Low, High, diff_av)
fire_preds_pred_ag <- bind_cols(new_data_pre_post_fire_pred,
predict(predator_global_ubms,
submodel = "state",
newdata = new_data_pre_post_fire_pred))
mean_prepost_pred <- fire_preds_pred_ag %>%
group_by(Species, Time) %>%
summarise(median = quantile(Predicted, 0.5),
LCI95 = quantile(Predicted, 0.025),
UCI95 = quantile(Predicted, 0.975)) %>%
ungroup() %>%
as.data.table() %>%
dcast(Species ~ Time, value.var = c("median", "LCI95", "UCI95")) %>%
mutate(diff_av = round(`median_Post-fire` - `median_Pre-fire`,2),
diff_LCI = round(`LCI95_Post-fire` - `LCI95_Pre-fire`, 2),
diff_UCI = round(`UCI95_Post-fire` - `UCI95_Pre-fire`, 2),
pre_fire = paste0(round(`median_Pre-fire`, 2), " [", round(`LCI95_Pre-fire`, 2), ", ", round(`UCI95_Pre-fire`, 2), "]"),
post_fire = paste0(round(`median_Post-fire`, 2), " [", round(`LCI95_Post-fire`, 2), ", ", round(`UCI95_Post-fire`, 2), "]")) %>%
dplyr::select(Species, pre_fire, post_fire, diff_av)
mean_sev_pred <- fire_preds_pred_ag %>%
filter(Time == "Post-fire") %>%
mutate(fire_class = case_when(scaled_fire_severity >= dq ~ "High",
scaled_fire_severity < dq ~ "Low")) %>%
group_by(Species, fire_class) %>%
summarise(median = quantile(Predicted, 0.5),
LCI95 = quantile(Predicted, 0.025),
UCI95 = quantile(Predicted, 0.975)) %>%
ungroup() %>%
as.data.table() %>%
dcast(Species ~ fire_class, value.var = c("median", "LCI95", "UCI95")) %>%
mutate(diff_av = round(`median_High` - `median_Low`,2),
High = paste0(round(`median_High`, 2), " [", round(`LCI95_High`, 2), ", ", round(`UCI95_High`, 2), "]"),
Low = paste0(round(`median_Low`, 2), " [", round(`LCI95_Low`, 2), ", ", round(`UCI95_Low`, 2), "]")) %>%
dplyr::select(Species, Low, High, diff_av)
References
R Session information
Show code
sessionInfo() %>% pander()
R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8
attached base packages: parallel, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: pROC(v.1.18.0), spatialEco(v.1.3-7), jagsUI(v.1.5.2), leafem(v.0.1.6), leaflet.extras(v.1.0.0), leaflet.providers(v.1.9.0), leaflet(v.2.0.4.1), lme4(v.1.1-27.1), Matrix(v.1.3-3), raster(v.3.5-2), DBI(v.1.1.1), maptools(v.1.1-1), sp(v.1.4-5), sf(v.1.0-2), VicmapR(v.0.1.8), pander(v.0.6.4), kableExtra(v.1.3.4.9000), ggbeeswarm(v.0.6.0), ggplot2(v.3.3.5), data.table(v.1.14.0), ubms(v.1.0.2.9016), unmarked(v.1.1.1), lattice(v.0.20-44), dbplyr(v.2.1.1) and dplyr(v.1.0.7)
loaded via a namespace (and not attached): backports(v.1.2.1), Hmisc(v.4.5-0), systemfonts(v.1.0.2), plyr(v.1.8.6), splines(v.4.1.0), crosstalk(v.1.1.1), chk(v.0.7.0), rstantools(v.2.1.1), inline(v.0.3.19), digest(v.0.6.28), htmltools(v.0.5.1.1), viridis(v.0.6.1), fansi(v.0.5.0), checkmate(v.2.0.0), magrittr(v.2.0.1), cluster(v.2.1.2), RcppParallel(v.5.1.4), matrixStats(v.0.61.0), svglite(v.2.0.0), prettyunits(v.1.1.1), jpeg(v.0.1-8.1), colorspace(v.2.0-2), blob(v.1.2.1), rvest(v.1.0.0), xfun(v.0.24), rgdal(v.1.5-23), callr(v.3.7.0), crayon(v.1.4.1), jsonlite(v.1.7.2), keyring(v.1.2.0), survival(v.3.2-11), glue(v.1.4.2), compute.es(v.0.2-5), gtable(v.0.3.0), webshot(v.0.5.2), V8(v.3.4.2), pkgbuild(v.1.2.0), rstan(v.2.21.2), scales(v.1.1.1), Rcpp(v.1.0.7), RPostgres(v.1.3.3), htmlTable(v.2.2.1), viridisLite(v.0.4.0), units(v.0.7-2), foreign(v.0.8-81), bit(v.4.0.4), proxy(v.0.4-26), Formula(v.1.2-4), stats4(v.4.1.0), StanHeaders(v.2.21.0-7), htmlwidgets(v.1.5.3), httr(v.1.4.2), RColorBrewer(v.1.1-2), wk(v.0.5.0), ellipsis(v.0.3.2), pkgconfig(v.2.0.3), loo(v.2.4.1), farver(v.2.1.0), nnet(v.7.3-16), sass(v.0.4.0), utf8(v.1.2.2), here(v.1.0.1), tidyselect(v.1.1.1), labeling(v.0.4.2), rlang(v.0.4.11), munsell(v.0.5.0), tools(v.4.1.0), cli(v.3.0.1), generics(v.0.1.0), ggridges(v.0.5.3), evaluate(v.0.14), stringr(v.1.4.0), yaml(v.2.2.1), processx(v.3.5.2), knitr(v.1.33), bit64(v.4.0.5), purrr(v.0.3.4), s2(v.1.0.6), lemon(v.0.4.5), nlme(v.3.1-152), xml2(v.1.3.2), compiler(v.4.1.0), rstudioapi(v.0.13), beeswarm(v.0.4.0), curl(v.4.3.2), png(v.0.1-7), e1071(v.1.7-7), RPostgreSQL(v.0.7), tibble(v.3.1.5), bslib(v.0.2.5.1), stringi(v.1.7.3), highr(v.0.9), ps(v.1.6.0), RSpectra(v.0.16-0), classInt(v.0.4-3), tinter(v.0.1.0), nloptr(v.1.2.2.2), vctrs(v.0.3.8), pillar(v.1.6.3), lifecycle(v.1.0.1), jquerylib(v.0.1.4), cowplot(v.1.1.1), latticeExtra(v.0.6-29), R6(v.2.5.1), MuMIn(v.1.43.17), KernSmooth(v.2.23-20), gridExtra(v.2.3), vipor(v.0.4.5), rjags(v.4-10), codetools(v.0.2-18), distill(v.1.2), boot(v.1.3-28), MASS(v.7.3-54), assertthat(v.0.2.1), rprojroot(v.2.0.2), withr(v.2.4.2), mgcv(v.1.8-35), hms(v.1.1.0), terra(v.1.4-20), rpart(v.4.1-15), grid(v.4.1.0), tidyr(v.1.1.3), coda(v.0.19-4), class(v.7.3-19), minqa(v.1.2.4), rmarkdown(v.2.9), downlit(v.0.2.1), lubridate(v.1.7.10) and base64enc(v.0.1-3)