Greater Gliders in Wombat State Forest

Using Mark-Recapture Distance Sampling to Estimate Abundance and Distribution of Greater Gliders in Wombat State Forest and Lerderderg State Park

Greater Glider

Show code

        function docReady(fn) {
            // see if DOM is already available
            if (document.readyState === "complete" || document.readyState === "interactive") {
                // call on next available tick
                setTimeout(fn, 1);
            } else {
                document.addEventListener("DOMContentLoaded", fn);
           }
        }    
        docReady(function() {
            history.replaceState({}, null, "https://justincally.github.io/blog/posts/2022-07-19-greater-gliders-in-wombat-sf/")
        });

Setup

During setup we load packages used in the analysis as well as the core data (greater glider observations and habitat variables). Basic data wrangling and data cleaning is also conducted.

Show code

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  cache = FALSE,
  message = FALSE, 
  warning = FALSE, 
  tidy.opts = list(width.cutoff = 60), tidy = TRUE
)
library(mapview)
mapview::mapviewOptions(fgb = FALSE)
library(cmdstanr)
library(data.table)
library(dplyr)
library(ggplot2)
library(bayesplot)
bayesplot_theme_set(theme_minimal())
library(cowplot)
library(VicmapR)
library(readr)
library(sf)
library(terra)
library(galah)
library(ggplot2)
library(ggspatial)
library(viridis)
library(readxl)
library(loo)
library(maptiles)
library(terrainr)
library(kableExtra)
library(pointblank)
library(tidyr)
library(stringr)
library(tidyterra)
library(ggridges)
library(maptiles)
library(tidyterra)
library(ggpmisc)
check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)

source(here::here("R", 'weighting_functions.R'))
source(here::here("_posts", '2022-07-19-greater-gliders-in-wombat-sf/GG_DataQuality/data-quality-functions.R'))

delwp_cols <- c(Teal = "#00B2A9",
                Navy = "#201547",
                Environment = "#CEDC00",
                `Climate Change` = "#FDDA24",
                Energy = "#0072CE",
                Water = "#71C5E8",
                Planning = "#642667",
                Infrastructure = "#AF272F",
                FFR = "#E57200",
                Corporate = "#201547")

delwp_cols_seq <- lapply(delwp_cols, function(x) {
  tinter::tinter(x, direction = "shades", steps = 6, adjust = 0.4)
})

read_excel_allsheets <- function(filename, tibble = FALSE) {
    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
}
galah_config(email = "justin.cally@delwp.vic.gov.au")
check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
register_knitr_engine(override = FALSE)

greater_glider_habitat_excel <- read_excel_allsheets("data/GG_Spotlight_H&D_Data.xlsx")

greater_glider_observations <- greater_glider_habitat_excel[[1]] %>%
  filter(Animal_no == 0 | !is.na(Distance_to_animal)) %>% # duplicate, no distance data
  left_join(greater_glider_habitat_excel[[2]], by = c("OldSiteID", "NewSiteID", "VisitNo")) %>%
  mutate(SiteNo = as.integer(as.factor(NewSiteID)), 
         `GG_height (estimate)` = as.numeric(`GG_height (estimate)`)) %>%
  rename(SiteID = NewSiteID) %>%
  arrange(SiteNo, VisitNo)

#### Greter Glider Storm Data ####
greater_glider_storm_data <- read_excel_allsheets("data/GG_ Spotlight_StormImpacts_Data.xlsx")

greater_glider_post_storm_obs <- greater_glider_storm_data[["Post-Storm GG Data"]] %>%
  filter(Animal_no == 0 | !is.na(Distance_to_animal)) %>% # duplicate, no distance data
  left_join(greater_glider_storm_data[["Site Data"]], by = c("OldSiteID", "NewSiteID", "VisitNo")) %>%
  mutate(SiteNo = as.integer(as.factor(NewSiteID))) %>%
  rename(SiteID = NewSiteID) %>%
  arrange(SiteNo, VisitNo)

# Add a study period column so we can make better arrays
greater_glider_all_sites  <- bind_rows(greater_glider_habitat_excel[[2]] %>%
                                               mutate(StudyPeriod = 1L), 
                                             greater_glider_storm_data[["Site Data"]] %>%
                                               mutate(StudyPeriod = 2L)) %>%
  arrange(NewSiteID, VisitNo, StudyPeriod) %>%
  mutate(SiteID = NewSiteID, 
         SiteNo = as.integer(as.factor(SiteID)))

greater_glider_all_observations <- bind_rows(greater_glider_observations %>%
                                               mutate(StudyPeriod = 1L), 
                                             greater_glider_post_storm_obs %>%
                                               mutate(StudyPeriod = 2L, 
                                                      ObserverX = as.numeric(ObserverX),
                                                      ObserverY = as.numeric(ObserverY),
                                                      GGy = as.numeric(GGy), 
                                                      GGx = as.numeric(GGx))) %>%
  dplyr::select(-SiteNo) %>%
  left_join(greater_glider_all_sites %>% 
              dplyr::select(SiteNo, SiteID, VisitNo, StudyPeriod))

greater_glider_all_site_simp <- greater_glider_all_sites %>%
  group_by(NewSiteID) %>%
  slice(1) %>%
  dplyr::select(SiteID = NewSiteID,
                `Transect start x`, `Transect start y`,
                `Transect end x`, `Transect end y`) %>%
  mutate(across(starts_with(c("Transect start", "Transect end")), ~as.numeric(gsub("°", "", .x)))) %>%
  mutate(`Transect end y` = -1*`Transect end y`,
         `Transect start y` = -1*`Transect start y`) %>%
  na.omit()

greater_glider_all_site_starts <- greater_glider_all_site_simp %>%
  sf::st_as_sf(coords = c("Transect start x", "Transect start y"))

greater_glider_all_site_ends <- greater_glider_all_site_simp %>%
  sf::st_as_sf(coords = c("Transect end x", "Transect end y"))

greater_glider_all_site_transects <- bind_rows(greater_glider_all_site_starts, greater_glider_all_site_ends) %>%
  group_by(SiteID) %>%
  summarise() %>%
  st_cast("LINESTRING") %>%
  sf::`st_crs<-`(4283) %>%
  ungroup() %>%
  left_join(greater_glider_all_sites %>%
              dplyr::select(SiteID, SiteNo) %>%
              unique()) %>%
  arrange(SiteNo)


greater_glider_observations_count <- greater_glider_all_site_transects %>% 
  st_drop_geometry() %>% 
  left_join(greater_glider_all_observations) %>%
  group_by(SiteID) %>% 
  mutate(vis = paste0(StudyPeriod, VisitNo)) %>%
  summarise(Count = sum(Count), 
            CountperVisit = Count/length(unique(vis))) %>%
  left_join(greater_glider_all_site_transects, by = "SiteID") %>%
  sf::st_as_sf() %>%
  sf::st_centroid() %>%
  mutate(Presence = case_when(Count > 0 ~ "Present", 
                              TRUE ~ "Absent"))

greater_glider_observations_count[is.na(greater_glider_observations_count)] <- 0

wombat_habitat_data <- read_excel_allsheets("data/GG_HabitatAssessment_H&D_data .xlsx")

# Missing sites for habitat data: 
# G12 and G54 
# G12 = G37
# G54 = G29
#### Temporary Fix ####

basal_sweep_summary <- wombat_habitat_data[["BasalSweep_All"]] %>%
  group_by(SiteID = NewSiteID, Species) %>%
  summarise(BasalArea = sum(`Basal Sweep`, na.rm = TRUE)) %>%
  dcast(SiteID ~ Species, value.var = "BasalArea")

colnames(basal_sweep_summary)[-1] <- paste("Basal Sweep", colnames(basal_sweep_summary)[-1])
colnames(basal_sweep_summary) <- gsub(" ", replacement = "_", colnames(basal_sweep_summary))

wombat_habitat_data_site_summary <- wombat_habitat_data$HabitatTree_All_ %>%
  rename(SiteID = NewSiteID) %>%
  mutate(DEADWOOD = case_when(DEADWOOD == ">30" ~ 30, 
                              TRUE ~ as.numeric(DEADWOOD)),
        EucRadiataDBH = case_when(SPECIES == "E. radiata" ~ DBH, 
                                  TRUE ~ 0)) %>%
  add_count(SiteID, SPECIES) %>%
  group_by(SiteID) %>%
  summarise(Trees = n(), 
            LivingTrees = sum(`DEAD/ALIVE` == "A", na.rm = T), 
            DeadTrees = sum(`DEAD/ALIVE` == "D", na.rm = T), 
            PropDead = DeadTrees/Trees, 
            MeanHeight = mean(HEIGHT, na.rm = T), 
            MedianDBH = median(DBH, na.rm = T), 
            SumDBH =  sum(DBH, na.rm = T),
            MedianDeadwood = median(DEADWOOD, na.rm = T), 
            SumDeadwood = sum(DEADWOOD, na.rm = T), 
            SumHollowSmall = sum(`HOLLOWS <100mm`, na.rm = T), 
            SumHollowMed = sum(`HOLLOWS 1-200mm`, na.rm = T), 
            SumHollowLarge = sum(`HOLLOWS >200mm`, na.rm = T), 
            SumHollow = sum(SumHollowSmall, SumHollowMed, SumHollowLarge),
            SumHollowMedLarge = sum(SumHollowMed, SumHollowLarge),
            EucRadiataProp = sum(SPECIES == "E. radiata")/Trees,
            EucObliquaProp = sum(SPECIES == "E. obliqua")/Trees,
            EucDalrympleanaProp = sum(SPECIES == "E. dalrympleana")/Trees,
            EucRadiataDBH = sum(EucRadiataDBH, na.rm = T),
            EucRadiataObliquaProp = EucRadiataProp + EucObliquaProp,
            EucRadiataDalrympleanaProp = EucRadiataProp + EucDalrympleanaProp,
            DominantSpecies = SPECIES[n == max(n)][1], 
            DominantSpecies2 = SPECIES[n == max(n[n!=max(n)])][1]) %>%
  ungroup() %>%
  mutate(DominantSpecies2 = case_when(is.na(DominantSpecies2) ~ DominantSpecies, 
                                      TRUE ~ DominantSpecies2), 
         DominantSpeciesGroup  = case_when(DominantSpecies %in% c("E. radiata", "E. obliqua") ~ DominantSpecies, 
                                      TRUE ~ "Other")) %>% 
  left_join(basal_sweep_summary, by = "SiteID")

wombat_habitat_data_site_summary[is.na(wombat_habitat_data_site_summary)] <- 0

WombatSF_SA <- readRDS("data/WombatSF_SA.rds") 
  
wombat_cov_stack <- rast("data/WombatSF_covariates.tif")

add_cov_stack <- rast(list.files("data/ClippedRasters", full.names = TRUE)) %>%
  project(wombat_cov_stack)

treecover <- rast("data/MOD44B.006_Percent_Tree_Cover_doy2020065_aid0001.tif")
treecover[treecover == 200] <- NA
treecover_proj <- project(treecover, wombat_cov_stack) %>%
  mask(terra::vect(WombatSF_SA)) %>%
  `names<-`("TreeDensity")

wombat_cov_stack[["TreeDensity"]] <- treecover_proj

storm_overlay_raster <- rast("data/storm_overlay_raster.tif")

fire_raster <- rast("data/fire_raster.tif")

VCF_diff <- rast("data/MOD44B/VCF_diff.tif")

wombat_cov_stack <- c(wombat_cov_stack, add_cov_stack, storm_overlay_raster, fire_raster, VCF_diff)

Wombat_EFG <- sf::st_read("data/Wombat_EFG_LMU", quiet = TRUE) %>%
  st_transform(3111)

area <- WombatSF_SA %>% 
  st_convex_hull() %>%
  st_transform(4283) %>%
  st_as_sf()

greater_glider_historical_records <- galah_call() %>%
  galah_identify(c("Petauroides volans")) %>%
  galah::galah_geolocate(area) %>%
  galah_filter(year >= 2012) %>%
  atlas_occurrences() %>%
  filter(!is.na(decimalLongitude) & !is.na(decimalLatitude)) %>%
  sf::st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4283) %>%
  mutate(Presence = "Historical Presence > 2012")

# Storm severity calculation 
storm_severity_site <- readRDS("data/storm_severity_site.rds")

The data quality for habitat variables is generally very good with only minor discrepancies for one record (dead or alive category not filled out). Several trees were included that do not necessarily match not historic records of species in the area (row 7), however these were checked and validated.

Site locations and historic records

Show code

wommap <- maptiles::get_tiles(
    WombatSF_SA %>% st_transform(4283),
    provider = "CartoDB.Positron",
    12,
    crop = TRUE,
    verbose = FALSE,
    forceDownload = FALSE
)



ContextPlot <- ggplot() +
  terrainr::geom_spatial_rgb(data = wommap, aes(x = x, y = y, r = red, g = green, b = blue)) +
  geom_sf(data = WombatSF_SA %>% st_transform(4283), alpha = 0.3, lwd = 0.1, stroke = 0.1, inherit.aes = FALSE) +
  # geom_sf(data = Wombat_EFG %>%
  #            filter(EFG_NAME %in% c("Forby Forest", "Moist Forest")) %>% 
  #           st_transform(4283), 
  #         aes(fill = EFG_NAME), lwd = 0.1, stroke = 0.1, inherit.aes = FALSE) +
  geom_sf(data = greater_glider_historical_records, aes(colour = Presence), 
          shape = 19, 
          size = 2, inherit.aes = FALSE) +
  geom_sf(data = greater_glider_observations_count, 
          aes(colour = Presence, size = CountperVisit), shape = 19, inherit.aes = FALSE) +
  scale_colour_manual(values = c(`Present` = "#99000d", 
                               `Absent` = "Black", 
                               `Historical Presence > 2012` = "orange")) +
  # scale_fill_manual(values = c(`Forby Forest` = delwp_cols[["Environment"]], 
  #                              `Moist Forest` = delwp_cols[["Teal"]])) +
  scale_size_continuous(name = "Average greater\ngliders per visit", breaks = c(2,4,6)) +
  # ylim(c(-37.57, -37.28)) +
  ylab("Latitude") +
  # xlim(c(144.07, 144.6)) +
  xlab("Longitude") +
  theme_classic()

ContextPlot

Figure 1: Distribution of greater glider sites surveyed and historical records. The average count of greater gliders detected per visit is indicated by the size of the point when present. Absent records are marked as a black dot, while historical records are shown as orange dots. The Grey area is the moist forest and forby forest that overlapped with the Greater Glider habitat distribution model considered in the sampling design. Thus the grey area is the scope of our analysis.

The total area available for the study is 374.55 km^2.

Formatting double observer data

The modelling of greater gliders using mark-recapture distance sampling requires several data objects to be constructed. Firstly, a data frame of individual greater gliders and whether they were detected by observer 1/2 as well as the perpendicular distance from the transect to the individual. We calculate the perpendicular distance through trigonometry functions using the bearing of the transect, bearing to the animal and the distance to the animal from the measurement point. Secondly, we generate a data frame for each visit with the number of detections/observations at that site. Thirdly, we group the records of greater gliders into distance bins in a matrix with n(site-visit) rows and n(distance-bins) columns. We set the parameters for the analysis (e.g., number of distance bins, maximum distance, breakpoints/midpoints of distance bins) below.

Below we also set the parameters for the distance-sampling model, such as number of distance bins and maximum distance.

Show code

n_distance_bins <- 8 # bins
max_distance <- 80 # you can bias the estimates if this is not set smartly. 
delta<- max_distance/n_distance_bins
# transect_length <- 0.5 # variable data in model
bin_breakpoints <- seq(0, max_distance, length.out = n_distance_bins + 1)
midpts <- seq(delta/2, max_distance, delta)

# survey info 
n_visit <- 2
n_studyperiod <- 2
n_site <- nrow(greater_glider_all_site_transects)

Show code

# Angle function
angle_diff <- function(theta1, theta2){
  theta <- abs(theta1 - theta2) %% 360 
  return(ifelse(theta > 180, 360 - theta, theta))
}


#################################################
# Greater gliders per site, per visit and per year 
GG_OP_all <- greater_glider_all_observations %>% 
  filter(Animal_no != 0) %>%
  mutate(GG_ID = paste(SiteNo, StudyPeriod, VisitNo, Animal_no, sep = "_"),
         RawAngleToAnimal = abs(angle_diff(as.numeric(Bearing_to_Animal), `Transect bearing`)), 
         AngleToAnimal = case_when(RawAngleToAnimal > 90 ~ RawAngleToAnimal-90,
                                   TRUE ~ RawAngleToAnimal), 
         HorizontalDistance = sin(AngleToAnimal * pi/180) * as.numeric(Distance_to_animal)) %>%
  group_by(SiteNo, StudyPeriod, VisitNo, Animal_no, GG_ID) %>%
  mutate(DReal = mean(as.numeric(Distance_to_animal), na.rm = T),
         DHReal = mean(HorizontalDistance)) %>%
  ungroup() %>%
  filter(DHReal <= max_distance) %>%
  dplyr::select(SiteID, SiteNo, StudyPeriod, VisitNo, Animal_no, GG_ID, DReal, M = Seen_by_Obs1, C = Seen_by_Obs2, DHReal) %>%
  group_by(GG_ID) %>%
  mutate(M = max(M), 
         C = max(C)) %>%
  distinct() %>%
  mutate(Union = 1) %>%
  arrange(SiteNo, VisitNo, StudyPeriod) %>%
  as.data.frame() 

GG_multi_obs <-  left_join(data.frame(SiteNo = unique(greater_glider_all_sites$SiteNo), 
                                      SiteID = unique(greater_glider_all_sites$SiteID)),
                           expand_grid(VisitNo = 1:n_visit, 
                                      StudyPeriod = 1:n_studyperiod), by = character()) %>%
  left_join(GG_OP_all, by = c("SiteID", "SiteNo", "StudyPeriod", "VisitNo")) %>%
  group_by(SiteID, SiteNo, StudyPeriod, VisitNo) %>%
  summarise(obs = sum(Union)) %>%
  ungroup() %>%
  left_join(greater_glider_all_sites %>%
              select(SiteNo, SiteID, VisitNo, StudyPeriod) %>%
              mutate(surveyed = 0)) %>%
  mutate(obs = coalesce(obs, surveyed)) %>%
  select(-surveyed) %>%
  arrange(SiteNo, VisitNo, StudyPeriod) %>%
  as.data.frame() %>%
  unique() 

############# TEMPORARY #############
# GG_multi_obs[is.na(GG_multi_obs)] <- 0

#####################################

all_site_n <- length(unique(GG_multi_obs$SiteNo))

y3d <- array(NA, dim = c(all_site_n, n_visit, n_studyperiod))

for(i in 1:n_visit) {
  y3d[,,i] <- GG_multi_obs %>% 
                      filter(StudyPeriod == i) %>%
                      as.data.table() %>%
                      dcast(SiteNo ~ VisitNo, value.var = "obs") %>%
    dplyr::select(-SiteNo) %>%
    as.matrix()
}

# derived parameter: number of site visits 
n_recs_all <- sum(!is.na(GG_multi_obs$obs))

# binned greater glider records 
GG_binned_all <- GG_multi_obs%>%
  dplyr::select(SiteNo, VisitNo, StudyPeriod) %>%
  unique() %>%
  left_join(GG_OP_all, Joining, by = c("SiteNo", "VisitNo", "StudyPeriod")) %>%
  mutate(Dbin = cut(DHReal, breaks=c(bin_breakpoints), include.lowest = FALSE, right = FALSE), 
         SiteVisit = paste(SiteNo, VisitNo, StudyPeriod, sep = "_")) %>%
  arrange(SiteNo)

y_all <- as.data.table(GG_binned_all) %>%
  dcast(SiteVisit + SiteNo + VisitNo + StudyPeriod ~ Dbin, value.var = "Union", fun.aggregate = sum) %>%
  arrange(SiteNo, VisitNo) %>%
  tibble::column_to_rownames("SiteVisit") %>%
  dplyr::select(-SiteNo, -VisitNo, -`NA`, -StudyPeriod) %>%
  `colnames<-`(midpts) %>%
  as.matrix()

y4d <- array(NA, c(all_site_n, length(midpts), n_visit, n_studyperiod))

for(yr in 1:n_studyperiod) {
  for(v in 1:n_visit) {
    ss <- paste(1:all_site_n, v, yr, sep ="_")
    y4d[,,v,yr] <- y_all[ss,]
  }
}
  
################################################

Exploration of abundance covariates

Habitat covariates

A range of data was collected at each site. This data includes key habitat variables on the stands/canopy trees at various points along each transect. Broadly, this data relates to the:

• senescence of the trees (living, dead amount of dead wood)
  
• tree age/maturity via diameter at breast height (DBH) and canopy height
  
• number of hollows (small, medium and large)
  
• tree species at each site
  
• basal sweeps al points along each transect

Below we group and summarise the habitat variables collected at each site. We visualise them by plotting the number of observations (y axis) against the value of the habitat variable (x axis).

Show code

# historic_disturbances_intersects <- readRDS("data/historic_disturbances_intersects.rds")
# Explore the relationship of GG obs and habitat covariates 
GG_obs_hab <- left_join(GG_multi_obs, wombat_habitat_data_site_summary, by = "SiteID") %>%
  arrange(SiteNo, VisitNo) #%>%
  # left_join(historic_disturbances_intersects %>% 
  #             dplyr::select(SiteID, LOGGED, FIRE))

GG_obs_hab_grp <- GG_obs_hab %>% 
  filter(StudyPeriod == 1) %>%
  dplyr::select(-StudyPeriod) %>%
  dplyr::select(-DominantSpecies, -DominantSpecies2, -DominantSpeciesGroup) %>%
  group_by(SiteID) %>%
  mutate(Siteobs = sum(obs)) %>%
  ungroup() %>%
  as.data.table() %>%
  melt(id.vars = c("SiteID", "SiteNo", "Siteobs")) %>%
  unique() %>%
  filter(!(variable %in% c("VisitNo", "obs", "Trees")))

GG_obs_hab_grp_cat <- GG_obs_hab %>% 
  filter(StudyPeriod == 1) %>%
  dplyr::select(-StudyPeriod) %>%
  dplyr::select(SiteID, SiteNo, DominantSpecies, DominantSpecies2, DominantSpeciesGroup, obs) %>%
  group_by(SiteID) %>%
  mutate(Siteobs = sum(obs)) %>%
  ungroup() %>%
  as.data.table() %>%
  melt(id.vars = c("SiteID", "SiteNo", "Siteobs")) %>%
  unique() %>%
  filter(!(variable %in% c("VisitNo", "obs", "Trees")))

contplot <- GG_obs_hab_grp %>%
  ggplot(aes(x = value, y = Siteobs)) +
  geom_point(shape = 21, fill = "DarkGreen", alpha = 0.7) +
  geom_smooth(method = "lm", level = 0.80, formula = y ~ x) +
  ylim(c(-2, 22)) +
  facet_wrap(~ variable, scales = "free_x", ncol = 3) +
  ylab("Greater gliders observed") +
  theme_bw()

catplot <- GG_obs_hab_grp_cat %>%
  ggplot(aes(x = value, y = Siteobs)) +
  geom_boxplot(fill = "DarkGreen", alpha = 0.7) +
  geom_point(shape = 21, fill = "DarkGreen", alpha = 0.7) +
  ylab("Greater gliders observed") +
  facet_wrap(~ variable, scales = "free_x", ncol = 3) +
  theme_bw() +
  theme(axis.text.x.bottom = element_text(angle = 270))

cowplot::plot_grid(contplot, catplot, ncol = 1, labels = "AUTO", rel_heights = c(0.7, 0.3))

Figure 2: Relationship between habitat variables collected at each site and the observed greater gliders at each site (summation over visits). (A) Shows variables that are on a continuous scale. A basic linear model line is fitted to help aid visualisation, this should not be statistically interpreted (grey bands represent 80% confidence intervals). (B) Shows the categorical variables (dominant and 2nd dominant species at each site).

Spatial covariates

In addition to the habitat covariates, there are a range of climate, topographic and modelled habitat covariates that are available on a wide spatial scale for the study area. As above, we investigate the relationship between these spatially extracted covariates and abundance.

Show code

GG_SpatialLocs <- left_join(GG_multi_obs, greater_glider_all_site_transects, by = "SiteNo") %>%
  sf::st_as_sf() %>%
  sf::st_transform(3111) 

# Add an 80m buffer to account for strip width of transect
# interpolate and take the mean
GG_SpatialData <- terra::extract(wombat_cov_stack, vect(GG_SpatialLocs %>% 
                                                          st_buffer(80, endCapStyle = "FLAT")), 
                                 na.rm = T, method = "bilinear", fun = "mean") %>% 
  mutate(storm_overlay = case_when(storm_overlay > 0 ~ TRUE, 
                                   TRUE ~ FALSE), 
         Fire_after_1982 = case_when(Fire_after_1982 == 0 ~ 0, 
                                      TRUE ~ 1))

GG_spatial_obs <- cbind(GG_multi_obs, GG_SpatialData)

GG_obs_spat_grp <- GG_spatial_obs %>% 
  filter(StudyPeriod == 1) %>%
  dplyr::select(-StudyPeriod) %>%
  dplyr::select(-TSLF_bin, -EVC) %>%
  group_by(SiteID) %>%
  mutate(Siteobs = sum(obs)) %>%
  ungroup() %>%
  as.data.table() %>%
  melt(id.vars = c("SiteID", "SiteNo", "Siteobs")) %>%
  unique() %>%
  filter(!(variable %in% c("VisitNo", "obs", "Trees", "ID")))

contplotspat <- GG_obs_spat_grp %>%
  ggplot(aes(x = value, y = Siteobs)) +
  geom_point(shape = 21, fill = "DarkGreen", alpha = 0.7) +
  geom_smooth(method = "lm", level = 0.80, formula = y ~ x) +
  ylim(c(-2, 22)) +
  facet_wrap(~ variable, scales = "free_x", ncol = 4) +
  ylab("Greater gliders observed") +
  theme_bw()

contplotspat

Figure 3: Relationships between greater glider counts at sites and a suite of spatially explicit covariates

Model building

Based on the above relationships between greater glider observations and habitat/spatial covariates we construct a suite of model formulas for the abundance component of the mark-recapture distance sampling model. The covariates the proved of interest were as follows…

Habitat covariates:

• PropDead: Proportion of measured trees that were dead
  
• MedianDBH: The median diameter at breast height of trees
  
• SumHollow: The total number of hollows on the measured trees
  
• MeanHeight: The average height for the PCQ (point centre quarter) measured trees
  
• EucRadiataDBH: The total DBH of trees measured (using PCQ) that were E. radiata
• SumDeadwood: Total amount of deadwood measure on the PCQ trees
  
• Basal_Sweep_E._radiata: Total basal sweep area for E. radiata
  
• Basal_Sweep_E._dalrympleana: Total basal sweep area for E. dalrymleana
  
• Basal_Sweep_E._obliqua: Total basal sweep area for E. obliqua

Spatial covariates:

• BIO15: Precipitation seasonality (Karger et al. 2017)
  
• BIO04: Temperature seasonality (Karger et al. 2017)
  
• BIO06: Mean daily minimum air temperature of the coldest month (Karger et al. 2017)
  
• TreeDensity: Density of trees according to satellite image analysis (Sexton et al. 2013) 
• TWIND: Topographic wetness index (J. P. Guerschman and Hill 2018)
  
• AHMI: Annual heat moisture index (Wagner et al. 2020)
  
• cond: Number of days with condensation index (CI) < 0 (Wagner et al. 2020)
  
• BareSoil: Estimate bare soil proportion (J. Guerschman 2014)
• NDVI: Normalized Difference Vegetation Index (Didan 2015)
  
• Elevation: Average elevation of transect (Department of Environment 2019)
  
• Nitrogen: Average modelled nitrogen in top 15cm of soil (J. P. Guerschman and Hill 2018)
  
• Fire after 1982: Whether a fire (since 1982) overlapped the transect (Department of Environment 2022)

In order to use these covariates, we should the correlation of these variables:

Show code

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)
}

GGCombinedData <- cbind(GG_SpatialData, GG_obs_hab) %>%
  group_by(SiteID, StudyPeriod) %>%
  slice(1) %>%
  ungroup()

GGCombinedDataCor <- GGCombinedData %>%
  dplyr::select(PropDead, MeanHeight, SumDeadwood, MedianDBH, SumHollowMedLarge, EucRadiataDBH, `Basal_Sweep_E._radiata`, Basal_Sweep_E._dalrympleana, Basal_Sweep_E._obliqua, EucRadiataDalrympleanaProp, Basal_Sweep_E._rubida, BIO15, BIO04, BIO06, TreeDensity, TWIND, AHMI, cond,BareSoil, SLOPE, NDVI, Elevation, Nitrogen, Fire_after_1982)

pairs(GGCombinedDataCor, lower.panel=panel.smooth, upper.panel=panel.cor)

Show code

# p <- pairs(GGCombinedData %>%
#         select(-LivingTrees, -Trees, -DeadTrees, -DominantSpecies, -DominantSpecies2, -DominantSpeciesGroup,
#                -SiteID, -SiteNo, -StudyPeriod, -VisitNo, -obs), lower.panel=panel.smooth, upper.panel=panel.cor)

Based on the above correlation plot we should avoid using NDVI, rain, BIO04, Elevation and BIO06.

Using these covariates we generate a set of model formulas and model matrices for analysis:

Show code

abundance_formulas <- list()
abundance_formulas[["Habitat_Only"]] <- ~ scale(PropDead) + scale(MeanHeight) + scale(SumDeadwood) + scale(MedianDBH) + scale(SumHollowMedLarge) + scale(EucRadiataDBH) + scale(Basal_Sweep_E._radiata) + scale(Basal_Sweep_E._dalrympleana) + scale(Basal_Sweep_E._obliqua) + scale(Basal_Sweep_E._rubida) #+ scale(Basal_Sweep_E._viminalis) #+ scale(EucRadiataDalrympleanaProp)
abundance_formulas[["Spatial_Only"]] <- ~ scale(cond) + scale(BIO15) + scale(TreeDensity) + scale(TWIND) + scale(AHMI) + scale(BareSoil) + scale(NDVI) + scale(Nitrogen) + Fire_after_1982
abundance_formulas[["Null"]] <- ~ 1

#### all 2 
abundance_formulas[["All"]] <- ~ scale(cond) + scale(BIO15) + scale(TreeDensity) + scale(TWIND) + scale(AHMI) + scale(BareSoil) + scale(NDVI) + scale(Nitrogen) + Fire_after_1982 + scale(PropDead) + scale(MeanHeight) + scale(SumDeadwood) + scale(MedianDBH) + scale(SumHollowMedLarge) + scale(EucRadiataDBH) + scale(Basal_Sweep_E._radiata) + scale(Basal_Sweep_E._dalrympleana) + scale(Basal_Sweep_E._obliqua) + scale(Basal_Sweep_E._rubida) #+ scale(Basal_Sweep_E._viminalis) # + scale(EucRadiataDalrympleanaProp)

# abundance_formulas[["Hab_Spat_CombinedAll"]] <- ~ scale(cond) + scale(BIO15) + scale(TreeDensity) + scale(TWIND) + scale(AHMI) + scale(E_radiata_sl) + scale(PropDead) + scale(MedianDBH) + scale(SumHollow) + scale(EucRadiataProp)

transition_formula <- ~ sqrt(storm_severity)
transition_formula_null <- ~ 1
transition_formula_vcf <- ~ scale(DeltaTreeCover) + scale(MeanHeight) + scale(MedianDBH)
# transition_formula_vcf <- ~ scale(sqrt(storm_severity)) + scale(MeanHeight) + scale(MedianDBH)
transition_formula_vcf_spatial <- ~ scale(DeltaTreeCover) 

Show code

GGCombinedData$storm_overlay <- factor(GGCombinedData$storm_overlay)

GGCombinedData_mod <- GGCombinedData %>%
  select(-storm_overlay) %>%
  left_join(bind_rows(storm_severity_site %>%
                        mutate(StudyPeriod = 2), 
                      storm_severity_site %>%
                        mutate(StudyPeriod = 1, 
                               storm_severity = 0,
                               storm_overlay = as.factor(FALSE))))

GGCombinedData_mod$Fire_after_1982 <- factor(GGCombinedData_mod$Fire_after_1982)

aggPred <- aggregate(wombat_cov_stack, fact = 10, fun = "mean", na.rm = TRUE)

ss_raster_smooth <- rast("data/ss_raster_smooth.tif")

terra::`add<-`(aggPred, ss_raster_smooth)

# Area of trees: used in predicting density area 
GG_SpatialArea_Tree <- exactextractr::coverage_fraction(aggPred[[1]], WombatSF_SA, crop = TRUE)[[1]]
names(GG_SpatialArea_Tree) <- "Area"
terra::`add<-`(aggPred, GG_SpatialArea_Tree)

predDF <- aggPred %>%
  as.data.frame(., xy= TRUE, cells = T) %>%
  na.omit() %>%
  mutate(storm_overlay = case_when(storm_overlay > 0 ~ TRUE, 
                                   TRUE ~ FALSE), 
         Fire_after_1982 = case_when(Fire_after_1982 == 0 ~ 0, 
                                   TRUE ~ 1))

predDF$Fire_after_1982 <- factor(predDF$Fire_after_1982)

predDF$storm_overlay <- factor(predDF$storm_overlay)

ModelData <- sapply(abundance_formulas, function(x) {
  
  drows <- GGCombinedData_mod %>%
                        filter(StudyPeriod == 1)
  
  mm <- model.matrix(x, data = bind_rows(drows, predDF))
  
  # mm_list <- list()
  # for(i in 1:n_studyperiod) {
  #   mm_list[[i]] <- mm[GGCombinedData_mod$StudyPeriod == i,] %>%
  #     as.matrix()
  # }
  return(mm[1:nrow(drows),,drop = FALSE]) 
  }, simplify = FALSE)

drows <- GGCombinedData_mod %>%
                        filter(StudyPeriod == 1)

tm <- model.matrix(transition_formula, data = bind_rows(drows, predDF))
  
tm_list <- list()
  for(i in 1:(n_studyperiod-1)) {
    tm_list[[i]] <- tm[1:nrow(drows),,drop = FALSE] %>%
      as.matrix()
  }

tm_vcf <- model.matrix(transition_formula_vcf, data = bind_rows(drows, predDF))
  
tm_vcf_list <- list()
  for(i in 1:(n_studyperiod-1)) {
    tm_vcf_list[[i]] <- tm_vcf[1:nrow(drows),,drop = FALSE] %>%
      as.matrix()
  }

tm_null <- model.matrix(transition_formula_null, data = GGCombinedData_mod)
  
tm_null_list <- list()
  for(i in 1:(n_studyperiod-1)) {
    tm_null_list[[i]] <- tm_null[GGCombinedData_mod$StudyPeriod == i+1,] %>%
      as.matrix()
  }
  

# filler
predDF$`PropDead` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`SumHollowMedLarge` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`MedianDBH` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`EucRadiataDBH` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`Basal_Sweep_E._radiata` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`Basal_Sweep_E._dalrympleana` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`Basal_Sweep_E._obliqua` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`Basal_Sweep_E._rubida` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`Basal_Sweep_E._viminalis` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`EucRadiataDalrympleanaProp` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`MeanHeight` <- rnorm(nrow(predDF), mean = 10, sd = 1)
predDF$`SumDeadwood` <- rnorm(nrow(predDF), mean = 10, sd = 1)



# set data to zero for prediction so model works 
PredictionData <- lapply(abundance_formulas, function(x) {
  pred.dfmm <- predDF %>%
  model.matrix(x, data = .) %>%
  as.data.frame()
  
  for(c in which(colnames(pred.dfmm) %in% labels(terms(abundance_formulas[["Habitat_Only"]])))) {
    pred.dfmm[,c] <- 0
  }
  
  return(as.matrix(pred.dfmm))
})

TransitionData <- predDF %>%
  model.matrix(transition_formula, data = .) %>%
  as.matrix()
  
TransitionData <- rep(list(TransitionData), n_studyperiod-1)

TransitionData_Null <- predDF %>%
  model.matrix(transition_formula_null, data = .) %>%
  as.matrix()
  
TransitionData_Null <- rep(list(TransitionData_Null), n_studyperiod-1)

TransitionData_VCF <- predDF %>%
  model.matrix(transition_formula_vcf, data = .) %>%
  as.matrix()
  
TransitionData_VCF <- rep(list(TransitionData_VCF), n_studyperiod-1)

Show code

missing_svp <- GG_multi_obs %>% 
  filter(is.na(obs))

missing_svp_array <- array(0, dim= c(n_site,
                                      n_visit,
                                      n_studyperiod))

for(i in 1:nrow(missing_svp)) {
  missing_svp_array[missing_svp$SiteNo[i], missing_svp$VisitNo[i], missing_svp$StudyPeriod[i]] <- 1
}

# y3d[is.na(y3d)] <- 0 # NAs break stan
y3d_long <- reshape2::melt(y3d) %>%
  `colnames<-`(c("SiteNo", "VisitNo", "StudyPeriod", "obs")) %>%
  filter(!is.na(obs))

stan_data <- list()
for(i in 1:length(abundance_formulas)) {
  
stan_data[[i]] <- list(N1 = sum(GG_OP_all$C), 
                       N2 = sum(GG_OP_all$M), 
                       M = GG_OP_all[GG_OP_all$C == 1, "M"], 
                       C = GG_OP_all[GG_OP_all$M == 1, "C"], 
                       delta = delta, 
                       transect_length = array(as.numeric(st_length(greater_glider_all_site_transects)/1000),dim = c(n_site, n_visit, n_studyperiod)),
                       D1 = GG_OP_all[GG_OP_all$C == 1, "DHReal"], 
                       D2 = GG_OP_all[GG_OP_all$M == 1, "DHReal"], 
                       visits = n_visit, 
                       studyperiods = n_studyperiod,
                       transitions = n_studyperiod-1,
                       n_site = n_site,
                       n_distance_bins = n_distance_bins,
                       midpts = midpts,
                       y = y4d,
                       n_obs = y3d_long$obs,
                       surveys = nrow(y3d_long),
                       period_sizes = table(y3d_long$StudyPeriod) %>% as.integer(),
                       visit_sizes = table(y3d_long$VisitNo) %>% as.integer(),
                       site_sizes = table(y3d_long$SiteNo) %>% as.integer(),
                       max_int_dist = as.integer(ceiling(max_distance)), 
                       poisson_ncb = if(i == 0) 1 else ncol(ModelData[[i]]), 
                       poisson_model_matrix = ModelData[[i]],
                       phi_ncb = ncol(tm_list[[1]]),
                       phi_model_matrix= tm_list,
                       max_distance = max_distance, 
                       npc = nrow(PredictionData[[i]]),
                       pred_model_matrix = PredictionData[[i]], 
                       phi_pred_matrix = TransitionData,
                       pred_area = 0.25, 
                       native_area = predDF$Area,
                       missing_data = missing_svp_array,
                       key_fun = 0, # 0 is half-normal and 1 is hazard
                       hs_df = 1, 
                       hs_df_global = 0.25,
                       hs_scale_global = 1/sqrt(nrow(y3d_long)), # ratio of expected non-zero to zero divided by total observation as per brms convention
                       hs_scale_slab = 1,
                       hs_df_slab = 2) 
}

stan_data_vcf <- stan_data
for(i in 1:length(abundance_formulas)) {
  stan_data_vcf[[i]][["phi_ncb"]] <- ncol(tm_vcf_list[[1]])
  stan_data_vcf[[i]][["phi_model_matrix"]] <- tm_vcf_list
  stan_data_vcf[[i]][["phi_pred_matrix"]] = TransitionData_VCF
}

stan_data_vcf_spatial <- stan_data
for(i in 1:length(abundance_formulas)) {
  stan_data_vcf_spatial[[i]][["phi_ncb"]] <- 2
  stan_data_vcf_spatial[[i]][["phi_model_matrix"]] <- list(tm_vcf_list[[1]][,1:2])
  stan_data_vcf_spatial[[i]][["phi_pred_matrix"]] = list(TransitionData_VCF[[1]][,1:2])
}

STAN model

Below is the STAN code for the mark-recapture distance sampling model. The model estimates a set of mark-recapture, distance-sampling and abundance parameters as well as predicting density back to the site level and to a spatial scale. Models are coded in the STAN programming language (Gelman, Lee, and Guo 2015) and executed using the cmdstanr R package (Gabry and Češnovar 2022).

Update: On Oct 20, 2022 NASA released 2021-2022 data for MOD44B (vegetation continuous fields), enabling us to calculate the change in tree cover before and after the storm

Horseshoe priors model

This model uses a regularised horsehshoe prior for the abundance submodel (Piironen and Vehtari 2017). Horseshoe priors are applied to all abundance covariates apart from the intercept, following the convention and code implementation in brms (Bürkner 2017, 2018, 2021). The regularised horseshoe prior allows us to model a larger suite of covariates despite the small value of n (sample size). Considering we have 42 unique sites with different model covariates and 20 model covariates, the regularisation approach provides a more parsimonious approach than model selection. In practice the horseshoe prior will strongly concentrate (regularise) uninformative parameters around zero, with informative parameters moving out of this regularised space. This allows us to intuitively interpret the effect (or lack thereof of all assessed covariates).

Show code

# model_negbin <- cmdstan_model(here::here("_posts", '2022-07-19-greater-gliders-in-wombat-sf', "stan", "modhs_nb.stan"))
# model_negbin_re <- cmdstan_model(here::here("_posts", '2022-07-19-greater-gliders-in-wombat-sf', "stan", "modhs_nb_re.stan"))
# model_negbin_nohs <- cmdstan_model(here::here("_posts", '2022-07-19-greater-gliders-in-wombat-sf', "stan", "mod_nb.stan"))
modhs <- cmdstan_model(here::here("_posts", '2022-07-19-greater-gliders-in-wombat-sf', "stan", "modhs.stan"))

Using the above code we fit the suite of models defined by the model formulas. Model fitting is conducted using 1,000 sampling iterations (1,000 warmup iterations) across four parallel chains. A set of initial values for parameters is also set before model fitting.

Show code

# MCMC settings
ni <- 1000
nw <- 1000
nt <- 1
nb <- 1000
nc <- 8

# fitlisths <- list()
# if(!file.exists("outputs/fitlisths[[1]].rds")) {
# 
# 
# for(i in 1:length(stan_data)) {
# fitlisths[[i]] <- modhs$sample(data = stan_data[[i]], 
#                   parallel_chains = nc, 
#                   show_messages = TRUE, 
#                   save_warmup = FALSE,
#                   iter_sampling = ni, 
#                   iter_warmup = nw, 
#                   # init = inits,
#                   adapt_delta = 0.99)
# fitlisths[[i]]$save_object(paste0("outputs/fitlisths[[",i,"]].rds"), compress = "xz")
# }
# names(fitlisths) <- paste0(names(abundance_formulas), "_hs")
# } else {
#   for(i in 1:length(stan_data)) {
#     fitlisths[[i]] <- readRDS(paste0("outputs/fitlisths[[",i,"]].rds"))
#   }
#   names(fitlisths) <- paste0(names(abundance_formulas), "_hs")
# }

fitlisths_nt <- list()
stan_data_null <- stan_data
if(!file.exists("outputs/fitlisths_nt[[1]].rds")) {


for(i in 1:length(stan_data_null)) {

  stan_data_null[[i]]$phi_ncb <- ncol(tm_null_list[[1]])
  stan_data_null[[i]]$phi_model_matrix <- tm_null_list
  stan_data_null[[i]]$phi_pred_matrix <- TransitionData_Null
  
  
fitlisths_nt[[i]] <- modhs$sample(data = stan_data_null[[i]], 
                  parallel_chains = nc, 
                  show_messages = TRUE, 
                  save_warmup = FALSE,
                  iter_sampling = ni, 
                  iter_warmup = nw, 
                  # init = inits2,
                  adapt_delta = 0.99)
fitlisths_nt[[i]]$save_object(paste0("outputs/fitlisths_nt[[",i,"]].rds"), compress = "xz")
}
names(fitlisths_nt) <- paste0(names(abundance_formulas), "_transition_null_hs")
} else {
  for(i in 1:length(stan_data)) {
    fitlisths_nt[[i]] <- readRDS(paste0("outputs/fitlisths_nt[[",i,"]].rds"))
  }
  names(fitlisths_nt) <- paste0(names(abundance_formulas), "_transition_null_hs")
}

##### VCF Transition #### 
fitlisths_vcf <- vector("list", length(stan_data_vcf))

if(!file.exists("outputs/fitlisths_vcf[[1]].rds")) {

for(i in 1:length(stan_data_vcf)) {
  
fitlisths_vcf[[i]] <- modhs$sample(data = stan_data_vcf[[i]], 
                  parallel_chains = nc, 
                  show_messages = TRUE, 
                  save_warmup = FALSE,
                  iter_sampling = ni, 
                  iter_warmup = nw, 
                  # init = inits,
                  adapt_delta = 0.99)

fitlisths_vcf[[i]]$save_object(paste0("outputs/fitlisths_vcf[[",i,"]].rds"), compress = "xz")
}
names(fitlisths_vcf) <- paste0(names(abundance_formulas), "_transition_vcf_hs")
} else {
  for(i in 1:length(stan_data_vcf)) {
    fitlisths_vcf[[i]] <- readRDS(paste0("outputs/fitlisths_vcf[[",i,"]].rds"))
  }
  names(fitlisths_vcf) <- paste0(names(abundance_formulas), "_transition_vcf_hs")
}

##### VCF Transition: only spatial #### 
fitlisths_vcf_spatial <- vector("list", length(stan_data_vcf_spatial))

if(!file.exists("outputs/fitlisths_vcf_spatial[[1]].rds")) {

for(i in 1:length(stan_data_vcf_spatial)) {
  
fitlisths_vcf_spatial[[i]] <- modhs$sample(data = stan_data_vcf_spatial[[i]], 
                  parallel_chains = nc, 
                  show_messages = TRUE, 
                  save_warmup = FALSE,
                  iter_sampling = ni, 
                  iter_warmup = nw, 
                  # init = inits,
                  adapt_delta = 0.99)

fitlisths_vcf_spatial[[i]]$save_object(paste0("outputs/fitlisths_vcf_spatial[[",i,"]].rds"), compress = "xz")
}
names(fitlisths_vcf_spatial) <- paste0(names(abundance_formulas), "_transition_vcf_spatial_hs")
} else {
  for(i in 1:length(stan_data_vcf_spatial)) {
    fitlisths_vcf_spatial[[i]] <- readRDS(paste0("outputs/fitlisths_vcf_spatial[[",i,"]].rds"))
  }
  names(fitlisths_vcf_spatial) <- paste0(names(abundance_formulas), "_transition_vcf_spatial_hs")
}

##### Negbin VCF Transition #### 
# fitlisths_negbin <- vector("list", length(stan_data_vcf))
# 
# if(!file.exists("outputs/fitlisths_negbin[[1]].rds")) {
# 
# for(i in 1:length(stan_data_vcf)) {
#   
# fitlisths_negbin[[i]] <- model_negbin$sample(data = stan_data_vcf[[i]], 
#                   parallel_chains = nc, 
#                   show_messages = TRUE, 
#                   save_warmup = FALSE,
#                   iter_sampling = ni, 
#                   iter_warmup = nw, 
#                   # init = inits,
#                   adapt_delta = 0.99)
# 
# fitlisths_negbin[[i]]$save_object(paste0("outputs/fitlisths_negbin[[",i,"]].rds"), compress = "xz")
# }
# names(fitlisths_negbin) <- paste0(names(abundance_formulas), "_fitlisths_negbin_hs")
# } else {
#   for(i in 1:length(stan_data_vcf)) {
#     fitlisths_negbin[[i]] <- readRDS(paste0("outputs/fitlisths_negbin[[",i,"]].rds"))
#   }
#   names(fitlisths_negbin) <- paste0(names(abundance_formulas), "_fitlisths_negbin_hs")
# }

Cross validation

We use leave-one-out cross-validation (LOO-CV); a Bayesian model evaluation process (Vehtari, Gelman, and Gabry 2017; Vehtari et al. 2020) to compare the habitat-only, spatial-only, all parameters and null model. As long as the model with ‘all’ parameters performs similarly to the habitat-only and spatial-only models it is suitable to use it exclusively for all analyses and predictions because of our use of the horseshoe regularisation. Importantly, the cross-validation also gives us the ability to compare the predictive performance against a null model. In doing so we see that the model with abundance covariates performs better than a null model. Noting that the null model will still include a random site effect and random observation effect.

Show code

fitlistCombined <- c(fitlisths_nt, fitlisths_vcf, fitlisths_vcf_spatial)
loofits <- lapply(fitlistCombined, function(x) {
  x$loo()
})

modcomp <- loo_compare(loofits)
# loo::loo_model_weights(loofits, method = "stacking")
modcomp %>% 
  as.data.frame() %>%
  kbl(format = "html", digits = 2, row.names = TRUE, 
      caption = "Model comparison using leave-one-out (LOO) cross-validation") %>%
  kable_styling(bootstrap_options = c("striped"))

Table 1: Model comparison using leave-one-out (LOO) cross-validation

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
All_transition_null_hs	0.00	0.00	-343.18	34.63	39.17	8.04	686.36	69.25
All_transition_vcf_spatial_hs	-1.38	1.56	-344.55	35.07	40.89	8.29	689.11	70.13
All_transition_vcf_hs	-3.48	4.26	-346.66	35.17	46.44	9.29	693.32	70.33
Spatial_Only_transition_null_hs	-17.04	10.68	-360.22	35.31	25.39	4.86	720.44	70.62
Spatial_Only_transition_vcf_spatial_hs	-18.10	10.56	-361.28	35.51	26.88	5.01	722.55	71.03
Spatial_Only_transition_vcf_hs	-20.42	11.05	-363.60	35.97	31.79	5.65	727.20	71.95
Habitat_Only_transition_null_hs	-25.69	13.77	-368.87	37.40	27.68	7.21	737.74	74.80
Habitat_Only_transition_vcf_spatial_hs	-27.02	14.23	-370.19	37.99	31.22	7.93	740.39	75.98
Habitat_Only_transition_vcf_hs	-28.98	13.17	-372.16	37.21	36.18	8.77	744.32	74.41
Null_transition_null_hs	-48.96	14.20	-392.14	35.33	6.65	1.39	784.28	70.65
Null_transition_vcf_hs	-49.59	14.62	-392.77	36.00	12.26	2.24	785.54	72.01
Null_transition_vcf_spatial_hs	-50.48	14.57	-393.66	35.86	9.38	2.07	787.31	71.71

Based on the above model selection table, the model with all covariates performs sufficiently well for it to be subsequently used as the model for further analyses. The elpd differences between the models are relatively minor (apart from the null model) and thus suggest that a single model does not necessarily perform better. However, in order to get good prediction across the landscape, we will use the spatially explicit horseshoe model (Spatial_Only_hs) for predictions over the geographical area of interest.

We also can quickly check whether a hazard detection function performs better than a half-normal detection function model. Using loo, we compare a half-normal and hazard model, with all other parameters kept the same:

Show code

hz_data <- stan_data_vcf[[3]]
hz_data$key_fun <- 1

fithz <- modhs$sample(data = hz_data, 
                  parallel_chains = nc, 
                  show_messages = T, 
                  save_warmup = FALSE,
                  iter_sampling = ni, 
                  iter_warmup = nw, 
                  # init = inits,
                  adapt_delta = 0.99, refresh = 0)

#> Running MCMC with 4 chains, at most 8 in parallel...
#> 
#> Chain 3 finished in 8.9 seconds.
#> Chain 1 finished in 9.5 seconds.
#> Chain 2 finished in 9.9 seconds.
#> Chain 4 finished in 10.3 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 9.7 seconds.
#> Total execution time: 10.4 seconds.

Show code

# compare loo
modcompdet <- loo_compare(list(`Half-normal` = fitlistCombined$Null_transition_vcf_hs$loo(), 
                               `Hazard` = fithz$loo()))

#get average detection
av_det_keyfun <- bind_rows(fitlistCombined[["Null_transition_vcf_hs"]]$summary("log_p"),
            fithz$summary("log_p")) %>%
  dplyr::select(`Average detection probability (p hat)` = mean, sd, q5, q95) %>%
  mutate_all(exp)

modcompdet %>% 
  bind_cols(av_det_keyfun) %>%
  mutate(`Detection Function` = c("Half-normal", "Hazard")) %>%
  dplyr::select(`Detection Function`, everything()) %>%
  kbl(format = "html", digits = 2, row.names = F, full_width = T, 
      caption = "Model comparison using leave-one-out (LOO) cross-validation for half-normal vs hazard detection functions") %>%
  kable_styling(bootstrap_options = c("striped"))  

Table 2: Model comparison using leave-one-out (LOO) cross-validation for half-normal vs hazard detection functions

Detection Function	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic	Average detection probability (p hat)	sd	q5	q95
Half-normal	0.00	0.00	-392.77	36.0	12.26	2.24	785.54	72.01	0.50	1.06	0.45	0.56
Hazard	-1.45	1.51	-394.22	36.1	12.81	2.27	788.43	72.19	0.42	1.16	0.34	0.51

Based on the above comparison there is no strong evidence to prefer the hazard model over the half-normal model. The remaining parts of the analysis will use a half-normal detection function. There also appears to be some evidence that the estimated average detection rate for the half-normal function is more precise.

Top model

Show code

fit <- fitlistCombined[["All_transition_vcf_hs"]]

Diagnostic checks

Based on the below MCMC chains, convergence and chain mixing occurred.

Show code

color_scheme_set("viridis")
mcmc_trace(fit$draws(c("beta_poisson","log_p","sigma", "y0", "b_distance")), facet_args = list(ncol = 3)) + theme_bw()

Posterior predictive checks

We compare the observed values to the predicted model values for several key statistics (proportion zeros, mean, standard deviation and maximum) regarding the count (observed greater gliders) variable. Modelled estimates of n_obs_pred should equal observed variables.

Show code

color_scheme_set("teal")
empty_obs <- paste0(missing_svp$SiteNo, ",", missing_svp$VisitNo, ",", missing_svp$StudyPeriod)
cols <- paste0("n_obs_pred[", empty_obs, "]")
pred_rep <- fit$draws(variable = c("n_obs_pred"), format = "draws_matrix") %>%
  as_tibble() %>%
  dplyr::select(-one_of(cols)) %>%
  as.matrix()

obs_nanarm <- y3d_long$obs

prop_zero<- function(x) mean(x == 0, na.rm = T) 
trimmed_mean <- function(x) mean(x, trim = 0.25)
pp1<- ppc_stat(obs_nanarm, pred_rep, stat=prop_zero, binwidth=0.015) +xlim(0.4, 0.75) + 
  ggtitle("Proportion zeros") + theme_bw() + theme(legend.position = "top")

pp2<- ppc_stat(obs_nanarm, pred_rep, stat="mean", binwidth=0.03) + ggtitle("Mean") + theme_bw() + theme(legend.position = "top")

pp3<- ppc_stat(obs_nanarm, pred_rep, stat="sd", binwidth=0.05) + ggtitle("Standard deviation") + theme_bw() + theme(legend.position = "top")

pp4<- ppc_stat(obs_nanarm, pred_rep, stat="max", binwidth=1) + xlim(0, 35) + ggtitle("Maximum") + theme_bw() + theme(legend.position = "top")

post1 <- pp_check(obs_nanarm, pred_rep, fun="scatter_avg") + xlim(0, 12) + ylim(0, 12) + theme_bw() + theme(legend.position = "top")
post2 <- ppc_bars(obs_nanarm, pred_rep, freq=FALSE) + theme_bw() + theme(legend.position = "top")
post3 <- ppc_rootogram(obs_nanarm, pred_rep) + theme_bw() + theme(legend.position = "top")
post4 <- ppc_error_scatter_avg(obs_nanarm, pred_rep) + theme(legend.position = "top")


cowplot::plot_grid(pp1,pp2,pp3,pp4,post3,post2,post1,post4, ncol = 2, labels = "AUTO", align = "hv", axis = "b")

Figure 4: Posterior predictive checks for the greater glider density model. Observed statistics that are compared to the posterior distributions for the n_obs_pred derived quantity. The figures correspond to statistics of (A) proportion of zeros, (B) mean count, (C) standard deviation, (D) maximum count, (E) average count on a scatter plot, (F) proportion of each count value, (G) square-root of the true count against the expected square-root count and (H) the average predictive error across the different count values.

Site-visit specific estimates and total population size

We can estimate the total number of greater gliders at each site-visit by combining the number of observed greater gliders with the modelled unobserved number of greater gliders. Each site-visit will have a unique modelled range of greater gliders at the site due to the covariates governing the abundance process as well as an error term for each visit and a random effect for each site.

Show code

transect_area <- sum(array(as.numeric(st_length(greater_glider_all_site_transects)/1000),dim = c(n_site, n_visit, n_studyperiod))[,1,1]*((max_distance/1000)*2))

no_double_visits <- c(4,17,22,23,24,28,30,35,36,38)
cols_n_dv <- c(paste0("n_per_site[", no_double_visits, ",", c(1), "]"), 
               paste0("n_per_site[", no_double_visits, ",", c(2), "]"))

color_scheme_set("blue")
cols_n <- paste0("n[", empty_obs, "]")
modelled_n <- fit$draws('n', format = "draws_matrix", inc_warmup = FALSE) %>%
  as_tibble() %>%
  dplyr::select(-one_of(cols_n)) %>%
  tibble::rowid_to_column(var = "draw") %>%
  as.data.table() %>%
  melt(id.var = "draw") %>%
  group_by(draw) %>%
  summarise(PopSize = sum(value)) 

pop1 <- modelled_n %>%
  ggplot() +
  geom_histogram(aes(x = PopSize), colour = "DarkBlue", fill = "LightBlue", lwd = 0.2, bins = 40) +
  scale_y_discrete(expand = c(0, 0)) +
  ylab("Density") +
  xlab("Total estimated population in \nsurveyed area across visits and years") +
  theme_classic() 

pop1_s1 <- fit$draws('n_per_site', format = "draws_matrix", inc_warmup = FALSE) %>%
  as_tibble() %>%
  dplyr::select(contains("1]"), -one_of(cols_n_dv) ) %>%
  tibble::rowid_to_column(var = "draw") %>%
  as.data.table() %>%
  melt(id.var = "draw") %>%
  group_by(draw) %>%
  summarise(`Estimated population size in searched area` = sum(value, na.rm = T)) %>%
  mutate(Season = "Pre-storm")

pop1_s2 <- fit$draws('n_per_site', format = "draws_matrix", inc_warmup = FALSE) %>%
  as_tibble() %>%
  dplyr::select(contains("2]"), -one_of(cols_n_dv) ) %>%
  tibble::rowid_to_column(var = "draw") %>%
  as.data.table() %>%
  melt(id.var = "draw") %>%
  group_by(draw) %>%
  summarise(`Estimated population size in searched area` = sum(value, na.rm = T)) %>%
  mutate(Season = "Post-storm")

change_on_transects <- bind_cols(pop1_s1 %>% 
                                   ungroup() %>%
                                   select(S1 = `Estimated population size in searched area`), 
                                 pop1_s2 %>% 
                                   select(S2 = `Estimated population size in searched area`)) %>%
  mutate(Change = S2-S1, 
         Percentage = 100*Change/S1, 
         Type = "Change") %>%
  group_by(Type) %>%
  summarise(`Average Change` = mean(Change), 
            `Q5 Change` = quantile(Change, 0.05), 
            `Q95 Change` = quantile(Change, 0.95), 
            `Average Percentage` = mean(Percentage), 
            `Q5 Perc` = quantile(Percentage, 0.05), 
            `Q95 Perc` = quantile(Percentage, 0.95))

abundance_on_transects <- bind_rows(pop1_s1, pop1_s2) %>%
  mutate(Density = "Density")

abundance_on_transects_summary <- abundance_on_transects %>%
  group_by(Season) %>%
  summarise(`Median Estimated Population Size in Searched Area` = median(`Estimated population size in searched area`), 
            `5%` = quantile(`Estimated population size in searched area`, 0.05), 
            `95%` = quantile(`Estimated population size in searched area`, 0.95))

abundance_on_transects_plot <- abundance_on_transects %>%
  ggplot(aes(x = `Estimated population size in searched area`, y = Density, fill = Season)) +
  geom_density_ridges(alpha = 0.75) +
  scale_fill_manual(values = c("Pre-storm" = delwp_cols_seq[["Teal"]][1], "Post-storm" = delwp_cols_seq[["FFR"]][1]), name = "Year") +
  scale_y_discrete(expand = c(0, 0), name = "", labels = "") +
  scale_x_continuous(name = paste0("Estimated population size in searched area (", round(transect_area, 2),  " km^2^)")) +
  theme_minimal() +
  theme(axis.title.x = ggtext::element_markdown())



idxs <- apply(expand.grid( 1:n_site, 1:n_visit, 1:n_studyperiod), 1, paste, collapse = ",")
notidxs <- paste0(missing_svp$SiteNo, ",", missing_svp$VisitNo, ",", missing_svp$StudyPeriod)

idxsf <- setdiff(idxs, notidxs)

n_obs_df <- cbind(variable = paste0("n[",idxsf, "]"), y3d_long) %>%
  left_join(greater_glider_all_site_transects %>%
              st_drop_geometry() %>%
              select(SiteID, SiteNo)) %>%
  rename(true = obs) %>%
  mutate(label = paste0("Site:", SiteID, ", Visit: ", VisitNo, ", Year:", StudyPeriod))

pop2 <- bayesplot::mcmc_intervals(fit$draws("n", format = "draws_matrix", inc_warmup = FALSE) %>% 
                                    as_tibble() %>%
                                    dplyr::select(-one_of(cols_n)) %>%
                                    as.matrix()) +
    geom_point(data = n_obs_df, aes(x = true, y = variable), fill = "Red", shape = 21, size = 2, position =  position_nudge(y = 0)) +
  scale_y_discrete(expand = c(0.01, 0), labels = n_obs_df$label) 

n_site_total <- fit$draws('n_per_site', format = "draws_matrix", inc_warmup = FALSE) 


colnames(n_site_total) <- GG_multi_obs %>% 
  arrange(StudyPeriod, SiteNo) %>%
  mutate(label = paste0("Site:", SiteID, ", Year:", StudyPeriod)) %>%
  pull(label) %>%
  unique()

GG_obs_mean <- GG_multi_obs %>%
  mutate(label = paste0("Site:", SiteID, ", Year:", StudyPeriod)) %>%
  group_by(label) %>%
  summarise(true = mean(obs)) %>%
  na.omit()

pop3 <- bayesplot::mcmc_intervals(n_site_total %>% 
                                    as_tibble() %>%
                                    select_if(~ !all(is.na(.))) %>%
                                    as.matrix()) +
    geom_point(data = GG_obs_mean, aes(x = true, y = label), 
               fill = "Red", shape = 21, size = 2, position =  position_nudge(y = 0)) 

plot_grid(plot_grid(pop1, pop3, labels = c("A", "B"), ncol = 1, rel_heights = c(0.35, 0.65)), pop2, ncol = 2, labels = c("", "", "C"))

Figure 5: Estimates of greater glider abundance at sites and across visits (unobserved and observed). The total number of Greater Gliders across sites and visits (unobserved and observed) is shown in (A). (B) Shows the site average estimates of Greater Gliders and (C) shows the estimated number of Greater Gliders present at a given visit.

Change between pre-storm and post-storm

On a site-by-site basis it is evident that there was a general reduction in Greater Glider abundance in the post-storm suverys compared with the pre-storm surveys.

Show code

site_areas <- greater_glider_all_site_transects %>% mutate(area = as.numeric(st_length(geometry)*160/1e6)) %>% st_drop_geometry()

intervals_data <- bayesplot::mcmc_intervals(n_site_total %>% 
                                    as_tibble() %>%
                                    select_if(~ !all(is.na(.))) %>%
                                    as.matrix())$data %>%
  left_join(GG_obs_mean, by = c("parameter" = "label")) %>%
  separate(col = "parameter", into = c("Site", "Year"), sep = ", ") %>%
  mutate(Year = case_when(Year == "Year:1" ~ "Pre-storm", 
                          TRUE ~ "Post-storm"),
         Site = stringr::str_remove_all(Site, "Site:"), 
         Site2 = paste(Site, Year, sep = ": ")) %>%
  mutate(Site_no_g = as.integer(stringr::str_remove_all(Site, "G"))) %>% 
  arrange(desc(Site_no_g)) %>%
  left_join(site_areas %>% 
              dplyr::select(Site = SiteID, area)) %>%
  mutate(m = m/area, 
         ll = ll/area,
         l = l/area,
         h = h/area,
         hh = hh/area) # convert to density

intervals_data$Site <- factor(intervals_data$Site, levels = unique(intervals_data$Site))

raw_counts_diff <- GG_multi_obs %>%
  group_by(SiteID, StudyPeriod) %>% 
  summarise(Count = sum(obs)) %>%
  mutate(Presence = case_when(Count > 0 ~ "Present", 
                              TRUE ~ "Absent")) %>%
  dplyr::select(SiteID, StudyPeriod, Count) %>%
  pivot_wider(names_from = StudyPeriod, values_from = Count) %>%
  mutate(diff = `2` - `1`) %>%
  mutate(change_direction = case_when(diff < 0 ~ "Reduction",
                                      diff > 0 ~ "Increase", 
                                      diff == 0 ~ "No change")) %>%
  ungroup()

raw_counts_diff$Site <- factor(raw_counts_diff$SiteID, levels = unique(raw_counts_diff$SiteID))

intervals_wide <- intervals_data %>%
  select(Site, Year, m) %>%
  pivot_wider(names_from = Year, values_from = m) %>%
  dplyr::select(Site, `Pre-storm\nModelled\nDensity\n(GGs/km²)` = `Pre-storm`, `Post-storm\nModelled\nDensity\n(GGs/km²)` = `Post-storm`) %>% 
  left_join(raw_counts_diff %>% 
              dplyr::select(Site, `Pre-storm\nCounts` = `1`, `Post-storm\nCounts` = `2`)) %>%
  arrange(desc(Site)) %>%
  mutate_if(is.numeric, round)

# table(raw_counts_diff$change_direction)

intervals_diff <- intervals_data %>% 
  group_by(Site) %>%
  arrange(Year) %>%
  summarise(mean_change = m[1] - m[2], 
            max_m = max(m, na.rm = T), 
            perc_change = 100*(m[2] - m[1])/m[2]) %>%
  ungroup() %>%
  mutate(direction = case_when(mean_change > 0 ~ "Increase", 
                               mean_change < 0 ~ "Decrease", 
                               mean_change == 0 & max_m == 0 ~ "No change, abundance = 0",
                               mean_change == 0 & max_m > 0 ~ "No change, abundance > 0",
                               is.na(mean_change) ~ "Incomplete data"))

intervals_diff$perc_change[is.infinite(intervals_diff$perc_change)] <- 100

intervals_diff_summary <- intervals_diff %>%
  group_by(direction) %>%
  summarise(`Average change in density (GG/km2)` = mean(mean_change), 
            `Average percentage change (%)` = mean(perc_change, na.rm = T),
            `Sites` = n(), 
            `Percentage of Total` = 100*n()/nrow(intervals_diff))

base_change_effect <- fit$summary("beta_phi[1]") %>%
  dplyr::select(mean, q5, q95) %>%
  as.numeric() %>% 
  exp() %>%
  round(2)

storm_change_effect <- fit$summary("beta_phi[2]") %>%
  dplyr::select(mean, q5, q95) %>%
  as.numeric() %>% 
  exp() %>%
  round(2)

height_change_effect <- fit$summary("beta_phi[3]") %>%
  dplyr::select(mean, q5, q95) %>%
  as.numeric() %>% 
  exp() %>%
  round(2)

dbh_change_effect <- fit$summary("beta_phi[4]") %>%
  dplyr::select(mean, q5, q95) %>%
  as.numeric() %>% 
  exp() %>%
  round(2)

dfnpc <- tibble(x = 500, y = 42, tb = list(intervals_wide))

intervalsplot <- intervals_data %>%
  ggplot() +
  geom_segment(aes(x = ll, xend = hh, y = Site, yend = Site, colour = Year), alpha = 0.8, size = 0.8) +
  geom_segment(aes(x = l, xend = h, y = Site, yend = Site, colour = Year), size = 3, alpha = 0.8) +
  geom_point(aes(x = m, y = Site, fill = Year, shape = Year, size = Year), colour = "gray70", alpha = 1) +
  geom_table(data = dfnpc,
             aes(x = x, y = y, label = tb), nudge_x = 220, nudge_y = 3.64,
             table.theme = ttheme_gtdefault(base_size = 12, padding = unit(c(0.015, 0.0148), "npc"))) +
  scale_fill_manual(values = c("Pre-storm" = delwp_cols_seq[["Teal"]][1], "Post-storm" = delwp_cols_seq[["FFR"]][1])) +
  scale_colour_manual(values = c("Pre-storm" = delwp_cols_seq[["Teal"]][3], "Post-storm" = delwp_cols_seq[["FFR"]][3])) +
  scale_shape_manual(values = c("Pre-storm" = 21, "Post-storm" = 23)) +
  scale_size_manual(values = c("Pre-storm" = 5.2, "Post-storm" = 4)) +
  scale_x_continuous(name = "Estimated density of greater gliders at each site (gliders per km^2^)", limits = c(0, 720)) +
  ggtitle("Change in site abundance", subtitle = "Modelled estimates") +
  theme_classic() +
  theme(panel.grid.major.y = element_line(), 
        axis.title.x = ggtext::element_markdown(), legend.position = "bottom", 
        plot.title = element_text(size =20, face = "bold", vjust = -14, hjust = 0.05), 
        plot.subtitle = element_text(size =18, face = "italic", vjust = -15, hjust = 0.045))

ggsave(plot = intervalsplot, filename = "outputs/intervalsplot.png", height = 18, width = 13, device = "png", units = "in")

ggsave(plot = intervalsplot, filename = "outputs/intervalsplot.pdf", height = 18, width = 13, device = "pdf", units = "in")

knitr::include_graphics("outputs/intervalsplot.png")

Figure 6: Change in greater glider estimates between pre- and post-storm surveys. Median estimates are shown with a circle point with 50% and 90 % confidence intervals shown as the bar and lines respectively. The table shows the median estimate from the modelled results as well as the raw counts for 2 visits to each site.

Show code

intervals_diff_summary %>%
  kbl("html", digits = 0, caption = "Summary of modelled site-specific changes to abundance between pre- and post-storm surveys") %>%
  kable_styling(bootstrap_options = c("striped"))

Table 3: Summary of modelled site-specific changes to abundance between pre- and post-storm surveys

direction	Average change in density (GG/km2)	Average percentage change (%)	Sites	Percentage of Total
Decrease	-34	59	20	48
Incomplete data	NA	NaN	10	24
Increase	11	-27	5	12
No change, abundance = 0	0	NaN	4	10
No change, abundance > 0	0	0	3	7

We estimated the intercept of the transition parameter at 0.55 [90% CI: 0.41 - 0.73]. And the change in tree cover transition effect as 0.85 [90% CI: 0.67 - 1.07]; which is a negligable relationship between increased tree cover at a given site and abundance. The height of trees at the sites also had a negligable effect on the change in Greater Glider abundance (0.97 [90% CI: 0.72 - 1.3]). However the median DBH of trees at the sites did have a relationship, siggesting a protective effect of mature habitat: 1.47 [90% CI: 1.06 - 2.04]

Detection functions

The mark-recapture distance sampling model relies on two mark recapture models (for observer 1 given observer 2 detection and observer 2 given observer 1 detection). The predicted detection recapture rate is plotted against the model covariate of perpendicular distance. The combined probability that at least observer 1 or 2 detected an individual can then be estimated. The intercept of this value then scales the half-normal detection function according to the point-independence method.

Show code

det_prob_pred <- fit$summary(c("MRUnion", "MR1", "MR2"))
det_prob_pred$Distance <- rep(seq(from = 0, to = max_distance, by = 1), 3)
det_prob_pred$Observer <- rep(c("Obs1 or Obs2", "Obs1", "Obs2"), each = max_distance+1)

g0_table <- det_prob_pred %>% filter(Distance == 0) %>% as.data.frame()

b_dist_effect <- fit$summary(c("b_distance")) %>% 
  dplyr::select(mean, q5, q95) %>%
  as.numeric() %>% 
  round(3)

sigma_effect <- fit$summary(c("sigma")) %>% 
  dplyr::select(mean, q5, q95) %>%
  as.numeric() %>% 
  round(1)

p_hat_effect <- fit$summary(c("log_p")) %>% 
  dplyr::select(mean, q5, q95) %>%
  as.numeric() %>% 
  exp() %>%
  round(3)

det_prob_MR <- det_prob_pred %>%
  ggplot(aes(x = Distance, y = mean, fill = Observer, colour = Observer)) +
  geom_ribbon(aes(ymin = q5, ymax = q95), alpha = 0.5) +
  geom_line() +
  viridis::scale_fill_viridis(discrete = T, option = "G", end = 0.8, guide = guide_legend(title.position = "top")) +
  viridis::scale_colour_viridis(discrete = T, option = "G", end = 0.8, guide = guide_legend(title.position = "top")) +
  ylim(c(0,1)) +
  scale_x_continuous(expand = c(0, 0)) +
  ggtitle("Mark Recapture Model") +
  ylab("Recapture probability [g()]") +
  xlab("Distance (m)") +
  theme_classic() +
  theme(legend.position = "top", legend.title.align = 0.5)

sigma <- fit$summary("sigma")$mean
# beta <- fit$summary("beta")$mean
# hnormsim <- rep(NA, 99)
# for (b in 1:99) {
#   hnormsim[b] <- exp(-(b+1)*b/(2*sigma^2))
# }
# hzsim <- rep(NA, 99)
# for (b in 1:99) {
#   hzsim[b] <- 1 - exp(-((b + -.5)/sigma)^(-beta))
# }

y0_summary <- fit$summary("y0")
det_curve <- fit$summary("DetCurve")

y_combined <- colSums(y_all) %>%
  as.data.frame() %>%
  `colnames<-`("Count") %>%
  tibble::rownames_to_column("Distance") %>%
  mutate(Distance = as.integer(Distance), 
         Prop = Count/(max(Count)),
         PropS = Prop*y0_summary$mean,
         PropSq5 = Prop*y0_summary$q5,
         PropSq95 = Prop*y0_summary$q95)


hnormsim_scaled <- data.frame(MRDS = det_curve$mean,
                              MRDS_LCI = det_curve$q5,
                              MRDS_UCI = det_curve$q95, 
                              Distance = 0:max_distance)


detection_plot_HN <- ggplot(aes(x = Distance), data = hnormsim_scaled) +
  geom_col(aes(y = PropS), fill = "grey70", colour = "grey20", width = delta, data = y_combined, alpha = 0.7) +
  # geom_errorbar(aes(ymin = PropSq5, ymax = PropSq95),  data = y_combined) +
  # geom_line(aes(y = HNS)) +
  geom_ribbon(aes(ymin = MRDS_LCI, ymax = MRDS_UCI), alpha = 0.5, fill = "Red") +
  geom_line(aes(y = MRDS), colour = "Red") +
  ggtitle("MRDS Model", subtitle = "Red = g(0) x p(dist)") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0,1)) +
  ylab("Detection probability (p)") +
  xlab("Distance (m)") +
  theme_classic()

det_plot_c <- plot_grid(det_prob_MR, detection_plot_HN, ncol = 2, labels = "AUTO")

ggsave(filename = "outputs/det_figure.png", plot = det_plot_c, width = 7, height = 4)

det_plot_c

Figure 7: Comparison of logistic regression for the mark-recapture model (A) and the combined mark-recapture distance-sampling model, using the point independence method (B). For (B) The histogram shows the proportional detection rates from the simulated data. The red line shows the model predictions for the detection function; accounting for g(0), or the detection probability being less than 1 at distance 0

Show code

g0_table %>%
  dplyr::select(Observer, Distance, mean, sd, q5, q95) %>% 
  kbl("html", digits = 2, caption = "Detection probabilities for observers 1, 2 and their union at g(0)/along the transect") %>%
  kable_styling(bootstrap_options = c("striped"))

Table 4: Detection probabilities for observers 1, 2 and their union at g(0)/along the transect

Observer	Distance	mean	sd	q5	q95
Obs1 or Obs2	0	0.96	0.01	0.93	0.98
Obs1	0	0.83	0.04	0.77	0.89
Obs2	0	0.76	0.05	0.68	0.83

Interestingly, from the mark-recapture model distance had no discernable effect on detection rates, with the b_distance coefficient estimated at -0.005 [90% CI: -0.017 - 0.007]. For the half-normal distance sampling model, sigma was estimated at 33.8 [90% CI: 30.3 - 38]. Ultimately, this led to an average detection probability for the search area of 0.498 [95% CI: 0.451 - 0.553].

Although not used in the top model the hazard model would produce a detection curve of:

Show code

det_curve_hz <- fithz$summary("DetCurve")

hazard_scaled <- data.frame(MRDS = det_curve_hz$mean,
                              MRDS_LCI = det_curve_hz$q5,
                              MRDS_UCI = det_curve_hz$q95, 
                              Distance = 0:max_distance)


detection_plot_HZ <- ggplot(aes(x = Distance), data = hazard_scaled) +
  geom_col(aes(y = PropS), fill = "grey70", colour = "grey20", width = delta, data = y_combined, alpha = 0.7) +
  # geom_errorbar(aes(ymin = PropSq5, ymax = PropSq95),  data = y_combined) +
  # geom_line(aes(y = HNS)) +
  geom_ribbon(aes(ymin = MRDS_LCI, ymax = MRDS_UCI), alpha = 0.5, fill = "Blue") +
  geom_line(aes(y = MRDS), colour = "Blue") +
  ggtitle("MRDS Model (Hazard function)", subtitle = "Blue = g(0) x p(dist)") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0,1)) +
  ylab("Detection probability (p)") +
  xlab("Distance (m)") +
  theme_classic()

detection_plot_HZ

Figure 8: Detection curve for a hazard model (not used in the analysis, for comparison only)

Key model parameters

Key model parameters from the top MRDS model are plotted below. The model parameters from the mark-recapture model are the detection probabilities of observer 1, 2 and their combined detection probability at a perpendicular distance of 0 (g(0)). The sigma value (shape of half-normal fitted function) is presented in the distance-sampling submodel section. The average detection probability at a site from 0m to the max distance (80m) on either side is shown as the p_hat value. Finally, the model estimates for the covariates governing the abundance process are presented alongside the intercept (mean log lambda).

Show code

vars <- unlist(stringr::str_remove_all(string = labels(terms(abundance_formulas[["All"]])), pattern =  "scale|\\)|\\("))
transition_vars <- unlist(stringr::str_extract_all(labels(terms(transition_formula_vcf)),  "(?<=\\().+?(?=\\))"))
# transition_vars <- "DeltaTreeCover"
posterior_vars <- list()
posterior_vars[["Mark recapture submodel"]] <- fit$draws(c("y01", "y02", "y0"))
dimnames(posterior_vars[["Mark recapture submodel"]])$variable <- c("p(Observer 1)", "p(Observer 2)", "p(Observer 1 U Observer 2)")
posterior_vars[["Distance sampling submodel"]] <- fit$draws("sigma")
dimnames(posterior_vars[["Distance sampling submodel"]])$variable <- c("sigma (half-normal scale term)") 
posterior_vars[["MRDS combined model"]] <- exp(fit$draws("log_p"))
dimnames(posterior_vars[["MRDS combined model"]])$variable <- c("p_hat (average detection probability)")
posterior_vars[["Abundance submodel"]] <- fit$draws("beta_poisson")
dimnames(posterior_vars[["Abundance submodel"]])$variable <- c("Intercept", vars)
posterior_vars[["Transition submodel"]] <- fit$draws("beta_phi")
dimnames(posterior_vars[["Transition submodel"]])$variable <- c("Intercept", transition_vars)

posterior_plots <- list()

for(i in 1:length(posterior_vars)) {

  color_scheme_set(delwp_cols_seq[[i]][1:6])
posterior_plots[[i]] <- mcmc_areas(
  posterior_vars[[i]], 
  # pars = c("cyl", "drat", "am", "sigma"),
  prob = 0.9, # 80% intervals
  prob_outer = 0.95, # 99%
  point_est = "mean"
) +
  labs(title = names(posterior_vars)[i])
}

posterior_plots[[1]] <- posterior_plots[[1]] + xlim(c(0,1))
posterior_plots[[2]] <- posterior_plots[[2]] + xlim(c(20,45))
posterior_plots[[3]] <- posterior_plots[[3]] + xlim(c(0,1))

height_scale <- unlist(lapply(posterior_vars, function(x) length(dimnames(x)$variable))) + 1

cowplot::plot_grid(plotlist = posterior_plots, ncol = 1, rel_heights = height_scale/sum(height_scale), align = "v")

Figure 9: Estimates of key parameters from the mark-recapture distance sampling model. Mean estimates are shown with a vertical line with 90% confidence intervals shown with darker shading. Tails are trimmed at 99% of the core density.

Show code

covnames <- c('Intercept', 'Condensing atmospheric conditions', 'BIO15: Precipitation seasonality', 'Tree cover (%)', 'TWIND', 'AHMI', 'Bare soil', 'NDVI', 'Soil nitrogen', 'Fire after 1982', 'Proportion of dead trees', 'Mean height', 'Sum Deadwood', 'Median DBH', 'Total hollows > 100 mm', 'DBH of _E. radiata_', 'Basal Sweep of _E. radiata_', 'Basal Sweep of _E. dalrympleana/cypellocarpa_', 'Basal Sweep of _E. obliqua_', 'Basal Sweep of _E. rubida_', 'Intercept', 'Change in Tree Cover %', 'Mean tree height pre-storm', 'Median DBH pre-storm')
beta_summary <- fit$summary(c("beta_poisson", "beta_phi")) %>% 
  mutate(Model = c(rep("Abundance", 20), rep("Transition", 4)), 
         Covariate = c("Intercept", vars, "Intercept", transition_vars)) %>%
  select(Model, Covariate, mean, median, sd, q5, q95) %>% 
  mutate(Covariate = covnames)

which_sig <- which((beta_summary$q5 * beta_summary$q95) > 0)

write.csv(beta_summary, "outputs/beta_summary.csv")

kbl(beta_summary, caption = "Model covariate estimates", digits = 2) %>% 
  kable_styling(bootstrap_options = "basic") %>%
  row_spec(which_sig, bold = T) %>%
  collapse_rows(1:2)

Table 5: Model covariate estimates

Model	Covariate	mean	median	sd	q5	q95
Abundance	Intercept	3.39	3.39	0.20	3.05	3.72
	Condensing atmospheric conditions	-0.24	-0.21	0.23	-0.66	0.06
	BIO15: Precipitation seasonality	0.18	0.14	0.20	-0.07	0.55
	Tree cover (%)	0.76	0.77	0.30	0.26	1.25
	TWIND	0.23	0.23	0.13	0.02	0.44
	AHMI	0.04	0.00	0.32	-0.45	0.62
	Bare soil	0.01	0.00	0.16	-0.24	0.28
	NDVI	-0.41	-0.39	0.31	-0.96	0.03
	Soil nitrogen	0.12	0.08	0.16	-0.08	0.43
	Fire after 1982	-1.84	-1.82	0.47	-2.64	-1.10
	Proportion of dead trees	0.42	0.42	0.10	0.25	0.59
	Mean height	0.09	0.07	0.13	-0.08	0.32
	Sum Deadwood	0.03	0.02	0.11	-0.14	0.22
	Median DBH	-0.07	-0.05	0.16	-0.36	0.17
	Total hollows > 100 mm	-0.01	0.00	0.09	-0.17	0.14
	DBH of E. radiata	0.38	0.39	0.18	0.07	0.68
	Basal Sweep of E. radiata	0.22	0.21	0.13	0.01	0.43
	Basal Sweep of E. dalrympleana/cypellocarpa	0.38	0.38	0.16	0.11	0.65
	Basal Sweep of E. obliqua	0.23	0.21	0.21	-0.04	0.60
	Basal Sweep of E. rubida	-0.45	-0.45	0.14	-0.68	-0.21
Transition	Intercept	-0.60	-0.59	0.17	-0.89	-0.32
	Change in Tree Cover %	-0.16	-0.16	0.14	-0.40	0.07
	Mean tree height pre-storm	-0.03	-0.03	0.18	-0.33	0.27
	Median DBH pre-storm	0.39	0.38	0.20	0.06	0.71

Abundance covariate effects

The effect of the abundance model covariates on greater glider density are best visualised by plotting the response (density) against the range of values of the covariate (e.g. mean height of trees).

Show code

marginal_effects_cmd <- function(draws, param, param_number, model_data, model_column, pwr = 1) {
  intercept = paste0(param, "[1]")
  poi = paste0(param, "[", param_number, "]")
  model_draws <- draws[, c(intercept, poi)]
  
  
  
  if(class(model_data[[model_column]]) == "factor") {
    sr <- c(0,1)
  } else {
    sr <- seq(from = -1, to = 1, by = 0.1)
  }
  
  effect <- matrix(ncol = length(sr), 
                   nrow = nrow(model_draws))
  c <- 1
  for(i in sr) {
  effect[,c] <- exp(model_draws[[1]] * 1 + model_draws[[2]] * i) # scaled so no need to add additional mean vals
  c = c+1
  }
  
  colnames(effect) <- sr
  effect <- as.data.table(effect)
  effect <- melt(effect, id.vars = NULL, measure.vars = colnames(effect))
    if(!class(model_data[[model_column]]) == "factor") {
  effect[, variable := scales::rescale(as.numeric(variable)^pwr, to = c(min(model_data[[model_column]]), 
                                                                    max(model_data[[model_column]])))]
    }
  effect$param <- model_column
  return(effect)
}

marginal_effects_plot_cmd <- function(data, col, factor = FALSE) {
  data_sum <- data %>%
    group_by(variable) %>%
    summarise(median = median(value), 
              q5 = quantile(value, 0.05), 
              q25 = quantile(value, 0.25), 
              q75 = quantile(value, 0.75), 
              q95 = quantile(value, 0.95))
  
  if(factor) {
    plot <- data %>%
    ggplot(aes(x = variable, y = value)) +
    geom_violin(draw_quantiles = c(0.05, 0.5, 0.95), fill = col, alpha = 0.5) +
    ylab("Density (gliders per km^2^)") +
    # ylim(c(0,47)) +
    theme_classic() + 
    theme(axis.title.y = ggtext::element_markdown())
  } else{
  
  plot <- data_sum %>%
    ggplot(aes(x = variable, y = median)) +
    geom_ribbon(aes(ymin = q5, ymax = q95), alpha = 0.25, fill = col) +
    geom_ribbon(aes(ymin = q25, ymax = q75), alpha = 0.25, fill = col) +
    geom_line(size = 1.5, colour = col) +
    ylab("Density (gliders per km^2^)") +
    # ylim(c(0,47)) +
    theme_classic() + 
    theme(axis.title.y = ggtext::element_markdown())
  }
  return(plot)
}

all_draws <- tidybayes::tidy_draws(fit)

model_marginal_effects <- list()
model_marginal_effects_transition <- list()

for(i in 2:(length(vars)+1)) {
model_marginal_effects[[vars[i-1]]] <- marginal_effects_cmd(all_draws, 
                                     param = "beta_poisson", 
                                     param_number = i, 
                                     model_data = GGCombinedData_mod, 
                                     model_column = vars[i-1]) 

}

for(i in 2:(length(transition_vars)+1)) {
model_marginal_effects_transition[[transition_vars[i-1]]] <- marginal_effects_cmd(all_draws, 
                                     param = "beta_phi", 
                                     param_number = i, 
                                     model_data = GGCombinedData_mod, 
                                     model_column = transition_vars[i-1], pwr = 1) # for sqrt transform

}

fire_diff <- model_marginal_effects[["Fire_after_1982"]] %>% 
  group_by(variable, param) %>% summarise(Median = median(value), 
                                          Q5 = quantile(value, 0.05), 
                                          Q95 = quantile(value, 0.95))

marg_plots <- list()
marg_plots[["cond"]] <- marginal_effects_plot_cmd(model_marginal_effects[["cond"]], 
                          col = delwp_cols[[1]]) +
  xlab("Condensing atmospheric conditions")

marg_plots[["bio15"]] <- marginal_effects_plot_cmd(model_marginal_effects[["BIO15"]], 
                          col = delwp_cols[[1]]) +
  xlab("BIO15: Precipitation seasonality")

marg_plots[["td"]] <- marginal_effects_plot_cmd(model_marginal_effects[["TreeDensity"]], 
                          col = delwp_cols[[1]]) +
  xlab("Tree cover (%)")
  
marg_plots[["twind"]] <- marginal_effects_plot_cmd(model_marginal_effects[["TWIND"]], 
                          col = delwp_cols[[1]]) +
  xlab("TWIND")

marg_plots[["ahmi"]] <- marginal_effects_plot_cmd(model_marginal_effects[["AHMI"]], 
                          col = delwp_cols[[1]]) +
  xlab("AHMI")

# marg_plots[["eradpred"]] <- marginal_effects_plot_cmd(model_marginal_effects[["E_radiata_sl"]], 
#                           col = delwp_cols[[1]]) +
#   xlab("Modelled _E. radiata_ habitat suitability") + 
#   theme(axis.title.x = ggtext::element_markdown())

marg_plots[["npv"]] <- marginal_effects_plot_cmd(model_marginal_effects[["BareSoil"]], 
                          col = delwp_cols[[1]]) +
  xlab("Bare soil")

marg_plots[["ndvi"]] <- marginal_effects_plot_cmd(model_marginal_effects[["NDVI"]], 
                          col = delwp_cols[[1]]) +
  xlab("NDVI")

marg_plots[["nitrogen"]] <- marginal_effects_plot_cmd(model_marginal_effects[["Nitrogen"]], 
                          col = delwp_cols[[1]]) +
  xlab("Soil nitrogen")

marg_plots[["fire"]] <- marginal_effects_plot_cmd(model_marginal_effects[["Fire_after_1982"]], 
                          col = delwp_cols[[1]], factor = TRUE) +
  xlab("Fire after 1982")

marg_plots[["dead"]] <- marginal_effects_plot_cmd(model_marginal_effects[["PropDead"]], 
                          col = delwp_cols[[2]]) +
  xlab("Proportion of dead trees")

marg_plots[["height"]] <- marginal_effects_plot_cmd(model_marginal_effects[["MeanHeight"]], 
                          col = delwp_cols[[2]]) +
  xlab("Mean height")

marg_plots[["sumdeadwood"]] <- marginal_effects_plot_cmd(model_marginal_effects[["SumDeadwood"]], 
                          col = delwp_cols[[2]]) +
  xlab("Sum deadwood")

marg_plots[["dbh"]] <- marginal_effects_plot_cmd(model_marginal_effects[["MedianDBH"]], 
                          col = delwp_cols[[2]]) +
  xlab("Median DBH")

marg_plots[["hol"]]<- marginal_effects_plot_cmd(model_marginal_effects[["SumHollowMedLarge"]], 
                          col = delwp_cols[[2]]) +
  xlab("Total hollows")

marg_plots[["erp"]] <- marginal_effects_plot_cmd(model_marginal_effects[["EucRadiataDBH"]],
                          col = delwp_cols[[2]]) +
  xlab("DBH of _E. radiata_") +
  theme(axis.title.x = ggtext::element_markdown())

marg_plots[["er"]] <- marginal_effects_plot_cmd(model_marginal_effects[["Basal_Sweep_E._radiata"]],
                          col = delwp_cols[[2]]) +
  xlab("Basal sweep of _E. radiata_") +
  theme(axis.title.x = ggtext::element_markdown())

marg_plots[["ed"]] <- marginal_effects_plot_cmd(model_marginal_effects[["Basal_Sweep_E._dalrympleana"]],
                          col = delwp_cols[[2]]) +
  xlab("Basal sweep of _E. dalrympleana/cypellocarpa_") +
  theme(axis.title.x = ggtext::element_markdown())

marg_plots[["eo"]] <- marginal_effects_plot_cmd(model_marginal_effects[["Basal_Sweep_E._obliqua"]],
                          col = delwp_cols[[2]]) +
  xlab("Basal sweep of _E. obliqua_") +
  theme(axis.title.x = ggtext::element_markdown())

marg_plots[["eru"]] <- marginal_effects_plot_cmd(model_marginal_effects[["Basal_Sweep_E._rubida"]],
                          col = delwp_cols[[2]]) +
  xlab("Basal sweep of _E. rubida_") +
  theme(axis.title.x = ggtext::element_markdown())

# marg_plots[["ev"]] <- marginal_effects_plot_cmd(model_marginal_effects[["Basal_Sweep_E._viminalis"]],
#                           col = delwp_cols[[2]]) +
#   xlab("Basal sweep of _E. viminalis_") +
#   theme(axis.title.x = ggtext::element_markdown())

# marg_plots[["erdp"]] <- marginal_effects_plot_cmd(model_marginal_effects[["EucRadiataDalrympleanaProp"]],
#                           col = delwp_cols[[2]]) +
#   xlab("Proportion of _E. dalrympleana_ and _E. radiata_") +
#   theme(axis.title.x = ggtext::element_markdown())

marg_plots[["ss"]] <- marginal_effects_plot_cmd(model_marginal_effects_transition[["DeltaTreeCover"]],
                          col = delwp_cols[[9]]) +
  labs(y = "Relative change in density", x = "Change in tree cover %")

marg_plots[["dthps"]] <- marginal_effects_plot_cmd(model_marginal_effects_transition[["MeanHeight"]],
                          col = delwp_cols[[9]]) +
  labs(y = "Relative change in density", x = "Mean tree height pre-storm")

marg_plots[["mdbhps"]] <- marginal_effects_plot_cmd(model_marginal_effects_transition[["MedianDBH"]],
                          col = delwp_cols[[9]]) +
  labs(y = "Relative change in density", x = "Median DBH pre-storm")



ggsave(plot = cowplot::plot_grid(plotlist = marg_plots, labels = "AUTO", ncol = 5, align = "h", greedy = TRUE, hjust = -0.3), filename = "outputs/effects_plot.pdf", height = 13, width = 22, device = "pdf")
ggsave(plot = cowplot::plot_grid(plotlist = marg_plots, labels = "AUTO", ncol = 5, align = "h", greedy = TRUE, hjust = -0.3), filename = "outputs/effects_plot.svg", height = 13, width = 22, device = "svg")
ggsave(plot = cowplot::plot_grid(plotlist = marg_plots, labels = "AUTO", ncol = 5, align = "h", greedy = TRUE, hjust = -0.3), filename = "outputs/effects_plot.png", height = 13, width = 22, device = "png")

cowplot::plot_grid(plotlist = marg_plots, labels = "AUTO", ncol = 5, align = "h")

Figure 10: Effect of spatial covariates (A - J) and habitat covariates (K - S) on greater glider density km²., as well as the effect of transition covariates (T) on the proportional change in Greater Glider density

Model predictions

Using the spatial covariates in the model we can predict the density of greater gliders across the study area of interest. Currently, the best performing model combines habitat and spatial covariates, however to predict across a landscape we use a regularised horseshoe model of six spatially explicit covariates. Thus predictions can only be based on the spatial covariates. That is, covariates that exist for all grid cells across the study area. This method allows us to calculate total greater glider abundance for the study area.

Show code

cells_n <- nrow(predDF)

yr1_preds <- paste0("n_pred[", 1:cells_n, ",1]")
yr2_preds <- paste0("n_pred[", 1:cells_n, ",2]")
preds_yr1 <- fitlistCombined$Spatial_Only_transition_vcf_spatial_hs$summary(yr1_preds, trimmed_mean = ~ mean(., trim = 0.05))
preds_yr2 <- fitlistCombined$Spatial_Only_transition_vcf_spatial_hs$summary(yr2_preds, trimmed_mean = ~ mean(., trim = 0.05))

preds_locations <- cbind(preds_yr1, predDF) %>%
  mutate(StudyPeriod = "Pre-storm") %>%
  bind_rows(cbind(preds_yr2, predDF) %>%
  mutate(StudyPeriod = "Post-storm")) %>%
  mutate(density = trimmed_mean*4, 
         density = case_when(density > 100 ~ density, 
                             TRUE~ density)) 

preds_locations$StudyPeriod <- factor(preds_locations$StudyPeriod, levels = c("Pre-storm", "Post-storm"))

wombat_tiles <- maptiles::get_tiles(
  WombatSF_SA %>% st_transform(3111) %>% st_buffer(1000) %>% st_transform(4283),
  provider = "OpenStreetMap",
  11,
  crop = TRUE,
  verbose = FALSE,
  forceDownload = FALSE
) 


predplot <- WombatSF_SA %>%
  st_transform(3111) %>%
  ggplot() +
  # geom_spatraster_rgb(data = wombat_tiles, inheret.aes = F) +
  geom_tile(aes(x=x, y=y, fill=density), data=preds_locations) + 
  geom_sf(fill=NA, lwd = 0.1) +
  annotation_scale(location = "br", width_hint = 0.25) +
  annotation_north_arrow(location = "tr", which_north = "true", height=unit(0.5, "in"),
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  scale_fill_viridis(discrete=FALSE, direction=-1, option="D", trans = "sqrt") +
  labs(x="Longitude",y="Latitude",fill="Density (km2)") +
  facet_wrap(~StudyPeriod, ncol = 1) +
  theme_bw() +
  # geom_sf(data = greater_glider_observations_count, 
  #         aes(colour = Presence, size = CountperVisit), shape = 19) +
  # scale_colour_manual(values = c(`Present` = "#99000d", 
  #                              `Absent` = "Black")) +
  # scale_size_continuous(name = "Average counts\nper visit") +
  theme(axis.title.x = element_text(face="bold", size=12),
        axis.title.y = element_text(face="bold", size=12),
        plot.subtitle = element_text(face="bold", size=12),
        legend.title = element_text(face="bold"))

predplot

Figure 11: Spatially explicit density of Greater Glider’s across the study area in Wombat State Forest and Lerderderg State Park.

Show code

# Difference 
preds_locations_diff <- preds_locations %>%
  arrange(StudyPeriod) %>%
  group_by(cell, x, y) %>%
  summarise(density_delta = density[2] - density[1], 
            density_perc = 100*density_delta/density[1])

predplot_diff <- WombatSF_SA %>%
  st_transform(3111) %>%
  ggplot() +
  geom_spatraster_rgb(data = rast(wombat_tiles)) +
  geom_tile(aes(x=x, y=y, fill=density_delta), data=preds_locations_diff) + 
  geom_sf(fill=NA, lwd = 0.1) +
  annotation_scale(location = "br", width_hint = 0.25) +
  annotation_north_arrow(location = "tr", which_north = "true", height=unit(0.5, "in"),
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  scale_fill_viridis(discrete=FALSE, direction=-1, option="A", limits = c(-110,5), guide = guide_colorbar(title.position = "top",barwidth = 12)) +
  labs(x="Longitude",y="Latitude",fill="Density change \\(gliders per km^2^)") +
  theme_bw() +
  # geom_sf(data = greater_glider_observations_count, 
  #         aes(colour = Presence, size = CountperVisit), shape = 19) +
  # scale_colour_manual(values = c(`Present` = "#99000d", 
  #                              `Absent` = "Black")) +
  # scale_size_continuous(name = "Average counts\nper visit") +
  theme(axis.title.x = element_text(face="bold", size=12),
        axis.title.y = element_text(face="bold", size=12),
        plot.subtitle = element_text(face="bold", size=12),
        legend.title = ggtext::element_markdown(face="bold"), 
        legend.position = "bottom")

predplot_diff

Figure 12: Spatially explicit changes in density of Greater Glider’s across the study area in Wombat State Forest and Lerderderg State Park

Show code

npredhat <- fitlistCombined$Spatial_Only_transition_vcf_hs$draws(yr1_preds, format = "draws_matrix") %>% 
  rowSums()

npredhat2 <- fitlistCombined$Spatial_Only_transition_vcf_hs$draws(yr2_preds, format = "draws_matrix") %>% 
  rowSums()

npreddiff <- npredhat2 - npredhat

data_for_change <- data.frame(Draw = 1:length(npredhat), 
                              Original = round(npredhat),
                              Change = round(npreddiff), 
                              name = factor("Change")) %>%
  # filter(Change > quantile(npreddiff, 0.005) & Change < quantile(npreddiff, 0.995)) %>%
  mutate(Direction = case_when(Change >= 0 ~ "Increase", 
                   TRUE ~ "Decrease"), 
         Percentage = 100*Change/Original) 

# quantile(data_for_change$Percentage, c(0.5,0.05,0.95))

percentile_zero <- ecdf(npreddiff)(0)

nph_df <- tibble(Time = c("Pre-storm", "Post-storm", "Change (post - pre)"), 
                 Median = round(c(median(npredhat), median(npredhat2), median(npreddiff))),
                     `5 % CI` = round(c(quantile(npredhat, 0.05), 
                                  quantile(npredhat2, 0.05), 
                                  quantile(npreddiff, 0.05))),
                     `95 % CI` = round(c(quantile(npredhat, 0.95), 
                                   quantile(npredhat2, 0.95), 
                                   quantile(npreddiff, 0.95))), 
                 SD = round(c(sd(npredhat), sd(npredhat2), sd(npreddiff))), 
                 CV = round(c(sd(npredhat)/mean(npredhat), 
                        sd(npredhat2)/mean(npredhat), 
                        sd(npreddiff)/mean(npredhat)), 2))

nph_df %>%
  kbl(format = "html", 
      caption = "Abundance estimates (with 5 % and 95 % confidence intervals) for Greater Gliders in the survey area of Wombat State Forest/Lerderderg State Park", row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("basic", "bordered"))

Table 6: Abundance estimates (with 5 % and 95 % confidence intervals) for Greater Gliders in the survey area of Wombat State Forest/Lerderderg State Park

Time	Median	5 % CI	95 % CI	SD	CV
Pre-storm	13549	10495	17776	2227	0.16
Post-storm	8560	6489	11565	1572	0.11
Change (post - pre)	-4966	-8158	-2139	1861	0.14

Show code

abundance_change <- data.frame(Draw = rep(1:length(npredhat), 2), 
                               Abundance = c(round(npredhat), round(npredhat2)),
                               Year = c(rep("Pre-storm", length(npredhat)), 
                                        rep("Post-storm", length(npredhat))), 
                               Density = "Density")

abundance_change_plot <- abundance_change %>%
  ggplot(aes(x = Abundance, y = Density, fill = Year)) +
  geom_density_ridges(alpha = 0.75) +
  scale_fill_manual(values = c("Pre-storm" = delwp_cols_seq[["Teal"]][1], "Post-storm" = delwp_cols_seq[["FFR"]][1]), name = "Year") +
  scale_y_discrete(expand = c(0, 0), name = "", labels = "") +
  scale_x_continuous(limits = c(quantile(abundance_change$Abundance, 0.000), quantile(abundance_change$Abundance, 0.999)), name = "Estimated population size in study area (374.9 km^2^)", labels = scales::comma) +
  theme_minimal() + 
  theme(axis.title.x = ggtext::element_markdown())

combined_abundance_plot <- cowplot::plot_grid(abundance_on_transects_plot, 
                                              abundance_change_plot, 
                                              align = "v", nrow = 2, labels = "AUTO")

ggsave(plot = combined_abundance_plot, filename = "outputs/abundance.png", height = 8, width = 8, device = "png", units = "in", bg = "white")
ggsave(plot = combined_abundance_plot, filename = "outputs/abundance.pdf", height = 8, width = 8, device = "pdf", units = "in")

We can also view the probable distribution of the change in Greater Glider abundance between pre- and post-storm by sampling the posterior distribution. Based on this data it is 99.8 % likely that abundance decreased.

Show code

data_for_change %>%
  ggplot(aes(x = Change, y = name, fill = factor(stat(quantile)))) +
  stat_density_ridges(
    geom = "density_ridges_gradient", calc_ecdf = TRUE, scale = 100, bandwidth = 200,
    quantiles = percentile_zero, quantile_lines = TRUE
  ) +
  scale_fill_manual(name = "Direction", values = c(`1` = "DarkRed", `2` = "DarkGreen"), labels = c("Decrease", "Increase")) +
  geom_vline(xintercept = 0, linetype = "solid", colour = "black") +
  scale_y_discrete(expand = c(0, 0), name = "", labels = "") +
  scale_x_continuous(limits = c(quantile(npreddiff, 0.00), quantile(npreddiff, 1)), name = "Change in Greater Glider abundance from pre- to post-storm", labels = scales::comma) +
  theme_minimal()

Figure 13: Posterior distribution of change in abundnace of Greater Gliders from pre- to post-storm. Green areas are posterior samples that suggest population increase and red are posterior samples that suggest population decrease.

Show code

gg_count_summary <- GGCombinedData_mod %>%
  select(SiteID, StudyPeriod, obs) %>%
  pivot_wider(id_cols = 1, names_from = 2, values_from = 3, names_prefix = "GG obs Y")
# Site habitat summary data 
gg_hab_summary <- left_join(gg_count_summary, GGCombinedData_mod %>% 
  select(SiteID, PropDead, MeanHeight, SumDeadwood, MedianDBH, SumHollowMedLarge, EucRadiataDBH, `Basal_Sweep_E._radiata`, `Basal_Sweep_E._dalrympleana`, `Basal_Sweep_E._obliqua`, `Basal_Sweep_E._rubida`) %>%
  distinct()) %>%
  mutate(No = str_remove(SiteID, "G") %>% as.numeric(), 
         PropDead = PropDead*100) %>%
  arrange(No) %>% 
  select(-No, -SiteID) %>%
  pivot_longer(everything()) %>%
  group_by(name) %>%
  summarise(across(everything(), list(mean = mean, sd = sd, min = min, max = max), na.rm = TRUE)) %>% mutate_if(is.numeric, ~round(.x, 1))

write_csv(gg_hab_summary, "outputs/habitat_summary.csv")