Methods
Data Wrangling
The following code shows the process we used to wrangle the data into a format suitable for the modelling approach used here. The final data object (umf_stack), is the data format used by unmarked and ubms.
Show code
# Data wrangling #
library(readxl)
library(dplyr)
library(ubms)
library(data.table)
#### Data Read in ####
mueller <- read_excel("data/Mueller_allSpp.xlsx")
other_sites <- read_excel("data/Drummer_Timboon_Genoa_allSpp.xlsx")
all_sites <- bind_rows(mueller, other_sites)
cov.names <- readxl::excel_sheets("data/Covariates.xlsx")
covariates <- lapply(cov.names, function(x) {
read_excel(path = "data/Covariates.xlsx", sheet = x)
})
names(covariates) <- cov.names
TSLF <- as.data.table(covariates[[1]]) %>%
melt(id.vars = c("Location", "Site")) %>%
rename(Year = variable, TSLF = value)
#### Site Covariate Data ####
# Get site covariates in a stack
sa_site_covs <- purrr::reduce(covariates[c(2,3,4,6)], left_join, by = c("Site", "Location")) %>%
mutate(id = Site) %>%
select(id, Location, Site, FireFreq, Elevation, Slope_Cat = CLASS, EFG_NAME)
rainfall_baiting <- purrr::reduce(covariates[c(7,8)], left_join, by = c("Year", "Location")) %>%
mutate(Year = Year - 2007)
# number of years 2008-2017 (scaled 1-10)
sa_site_covs_stack <- sapply(sa_site_covs, rep.int, times=5) %>% as.data.frame()
sa_site_covs_stack$Year <- rep(c(1L:4L,10L), each = nrow(sa_site_covs))
sa_site_covs_stack <- sa_site_covs_stack %>%
left_join(rainfall_baiting, by = c("Location", "Year"))
# Rep again for species
sa_site_covs_stack_species <- sapply(sa_site_covs_stack, rep.int, times=4) %>% as.data.frame()
sa_site_covs_stack_species$Species <- rep(c("BTP", "LNB", "LNP", "SBB"), each = nrow(sa_site_covs_stack))
sa_site_covs_stack_species_dm <- sapply(sa_site_covs_stack_species, rep.int, times=2) %>% as.data.frame()
sa_site_covs_stack_species_dm$Detection_method <- rep(c("cam", "ht"), each = nrow(sa_site_covs_stack_species))
sa_site_covs_stack_species_dm <- left_join(sa_site_covs_stack_species_dm, TSLF, by = c("Location", "Year", "Site"))
#### Presence-Absence Detection ####
setDT(all_sites)
y_melt <- melt(all_sites, id.vars = c("Location", "Site", "Species")) %>%
tidyr::separate(col = variable, into = c("UID", "obs")) %>%
tidyr::separate(col = "UID", into = c("Detection_method", "Year"), sep = -2) %>%
mutate(obs = as.integer(obs),
obs_group = case_when(Detection_method == "cam" & obs < 8 ~ 5L,
Detection_method == "cam" & obs >= 8 & obs < 15 ~ 6L,
Detection_method == "cam" & obs >= 15 & obs < 22 ~ 7L,
Detection_method == "cam" & obs >= 22 ~ 8L, #1-4 observations can also be used
TRUE ~ as.integer(obs)),
Year = as.integer(Year)-7) %>%
# dplyr::rowwise() %>%
# mutate(detection = case_when(is.na(value) ~ NA_integer_,
# value == 0 ~ 0L,
# value > 0 ~ 1L)) %>%
group_by(Location, Site, Species, Year, Detection_method, obs_group) %>%
summarise(detection = max(as.numeric(value), na.rm = TRUE)) %>%
ungroup() %>%
distinct() %>%
arrange(obs_group) %>%
mutate(obs_group = paste0("obs_", obs_group),
detection = as.integer(detection)) %>%
mutate(detection = case_when(detection > 0 ~ 1L,
TRUE ~ detection))
setDT(y_melt)
y_stack <- data.table::dcast(y_melt, Location + Site + Species + Year ~ obs_group, value.var = "detection")
#### Probability of Detections ####
# Detection
method_stack <- y_stack %>%
mutate(dplyr::across(c("obs_1", "obs_2", "obs_3", "obs_4"), ~ "ht")) %>%
mutate(dplyr::across(c("obs_5", "obs_6", "obs_7", "obs_8"), ~ "cam"))
#### Data for Model ####
# Data same lengths
join_cols <- c("Location", "Site", "Year", "Species")
#### Inner join ####
sa_site_covs_stack_species_dm <- sa_site_covs_stack_species_dm %>% select(-Detection_method) %>% distinct()
joined_data_stack_covs <- left_join(sa_site_covs_stack_species_dm %>%
mutate(Year = as.numeric(Year)), y_stack, by = join_cols) %>%
dplyr::select(id, FireFreq, Elevation, Slope_Cat, EFG_NAME, Site, Location, Year, Species, TSLF, Baiting, `6mthlag`) %>%
mutate(Elevation = as.numeric(Elevation), Year = as.numeric(Year), TSLF = as.numeric(TSLF), `6mthlag` = as.numeric(`6mthlag`),
Baiting = as.numeric(Baiting), FireFreq = as.numeric(FireFreq))
joined_data_stack_y <- left_join(sa_site_covs_stack_species_dm %>%
mutate(Year = as.numeric(Year)), y_stack, by = join_cols) %>%
dplyr::select(contains("obs_"))
joined_data_stack_method <- left_join(sa_site_covs_stack_species_dm %>%
mutate(Year = as.numeric(Year)), method_stack, by = join_cols) %>%
dplyr::select(contains("obs_"))
#### Final Formatting ####
umf_stack <- unmarkedFrameOccu(y=joined_data_stack_y, siteCovs=joined_data_stack_covs,
obsCovs=list(method=joined_data_stack_method))
saveRDS(umf_stack, "data/umf_stack.rds")
saveRDS(y_melt, "data/y_melt.rds")
Statistical approach
The approach we took to account for the imperfect detections of species that arise through hair tube and camera trap sampling was to use an occupancy modelling framework. This approach considers the probability that a species will be present at some sites but will go undetected, and this probability is explicitly incorporated into the presence estimates.
Using the ubms R package (Kellner 2021), we fit a “stacked” model, which is a single-season model with pseudo-replication (i sites will be repeated in the data over j years). The effects of pseudo-replication are controlled by having site (n = 150) and location (n = 4) as random variables. This “stacked” approach was chosen over a multi-season dynamic occupancy model (Darryl I. MacKenzie et al. 2003) for several reasons. Firstly, multi-season models can require large data sets, given no data was available between 2012 and 2016, multi-season models may not be suitable for use after 2011. Secondly, the primary aim of this study is to assess the effect of fox baiting on the occupancy of the target native mammals; hence the transition probabilities that are generated in multi-season models are not of direct interest for this study.
The data used in this analysis comprised of three components. (i) Covariate data of various environmental and ecological variables for given sites and for given time periods (years). (ii) Detection probability data comprised of a categorical variable (hair tube and camera traps) that was used in the estimation of detection probability as part of the modelling approach. (iii) Presence-absence data for each of the species, across sites and over observation periods. The presence-absence data was simplified to fit a structure that is consistent across the methods. That is, counts were grouped by week and transformed to presence (x > 0) or absence (x = 0). This ensured that each method had four observation periods per deployment (weeks 1,2,3 and 4).
The models developed in this analysis are standard occupancy models based on a zero-inflated binomial distribution (Darryl I. MacKenzie et al. 2017). The models consist of two processes: (i) an occupancy state process that models the observed occupancy at a site (\(z_i\)) based on the true probability of occupancy for that site \(\psi_i\):
\[z_i \sim Bernoulli(\psi_i)\] (ii) an observation process modelling the observation of a species at a site and a detection probability (\(y_{ij}\)), which is dependent on probability of detection at a given site and condition \(p_{ij}\):
\[y_{ij} | z_j \sim Bernoulli(z_i \times p_{ij})\] The site covariates of elevation, species, baiting, year, fire frequency and rainfall were modelled as fixed effects in the occupancy state process. Site and location were included as random effects within the occupancy state process. The method of detection (camera trap or hair tube) covariates were included in modelling the observation process. Models were developed using the ubms R package (Kellner 2021); which fits Bayesian hierarchical models with Stan (Carpenter et al. 2017).
We used model selection to compare the goodness of fit of three candidate models. One global model contained the covariates metioned above, while two simplified models removed (i) fire frequency and rainfall and (ii) fire frequency, rainfall and the interaction effect of baiting and species The global and simplified models can be compared based on their predictive accuracy (Vehtari, Gelman, and Gabry 2017). Given that we are implementing a Bayesian approach to the occupancy models, the computational costs of extensive model selection and dredging is high. Here we use model selection to find a suitable model that is well-founded in ecological theory from a select few well-reasoned models. The model selection implemented here found that the global model was the best-fitting model (additional details in Supplementary Material). The model formula for the top model (in ubms syntax) is:
~method ~ scale(Elevation)*Species
+ scale(Baiting)*Species
+ scale(Year)*Species
+ scale(FireFreq)*Species
+ scale(Rainfall)
+ (1|Location)
+ (1|Site)This syntax corresponds to a the following models for the occupancy and detection sub-process, where \(\alpha\) is the random effect of site \(j\) and \(\gamma\) is the random effect of location \(k\).
\[\alpha_j \sim normal(\mu_{site}, \tau_{site})\] \[\gamma_k \sim normal(\mu_{location}, \tau_{location})\] \[logit(\psi_{ijk}) \sim \beta_1Species_i + \beta_2Elevation_i + \beta_3(Elevation_i \times Species_i) \\+ \beta_4Baiting_i + \beta_5(Baiting_i \times Species_i) \\+ \beta_6Year_i + \beta_7(Year_i \times Species_i) \\+ \beta_8Firefreq_i + \beta_9(Firefreq_i \times Species_i) \\+ \beta_{10}Rainfall_i \\+ \alpha_{j} + \gamma_{k}\] \[logit(p_i) \sim \beta_{11}Method_i\]
Data Summary
The presence-absence data, summarised by Year, detection method and species is as follows. Each row should total 960, which is the product of 4 locations, 60 sites per locations and 4 observations per site (4 X 60 X 4).
Show code
PresenceAbsenceData <- y_melt[, c('detection', 'Detection_method', 'Year', 'Species')] %>%
mutate(Year = Year +2007L,
Species = as.character(Species)) %>%
group_by(Detection_method, Species, Year) %>%
summarise(`0` = sum(detection == 0, na.rm = TRUE),
`1` = sum(detection == 1, na.rm = TRUE),
`NA` = sum(is.na(detection)))
PresenceAbsenceData %>%
kbl(format = "html", caption = "Total presence and absence counts for the four target species over the study period") %>%
kable_styling(full_width = FALSE) %>%
column_spec(1:3, bold = TRUE) %>%
collapse_rows(columns = 1:2)
| Detection_method | Species | Year | 0 | 1 | NA |
|---|---|---|---|---|---|
| cam | BTP | 2010 | 180 | 12 | 768 |
| 2011 | 131 | 13 | 816 | ||
| 2017 | 854 | 66 | 40 | ||
| LNB | 2010 | 172 | 20 | 768 | |
| 2011 | 135 | 9 | 816 | ||
| 2017 | 906 | 14 | 40 | ||
| LNP | 2010 | 156 | 36 | 768 | |
| 2011 | 116 | 28 | 816 | ||
| 2017 | 903 | 17 | 40 | ||
| SBB | 2010 | 186 | 6 | 768 | |
| 2011 | 136 | 8 | 816 | ||
| 2017 | 917 | 3 | 40 | ||
| ht | BTP | 2008 | 711 | 249 | 0 |
| 2009 | 703 | 257 | 0 | ||
| 2010 | 676 | 284 | 0 | ||
| 2011 | 622 | 338 | 0 | ||
| LNB | 2008 | 955 | 5 | 0 | |
| 2009 | 957 | 3 | 0 | ||
| 2010 | 958 | 2 | 0 | ||
| 2011 | 951 | 9 | 0 | ||
| LNP | 2008 | 721 | 239 | 0 | |
| 2009 | 736 | 224 | 0 | ||
| 2010 | 778 | 182 | 0 | ||
| 2011 | 744 | 216 | 0 | ||
| SBB | 2008 | 873 | 87 | 0 | |
| 2009 | 881 | 79 | 0 | ||
| 2010 | 919 | 41 | 0 | ||
| 2011 | 910 | 50 | 0 |
Fit models
Fit a global model with numerous parameters. From this broad global model, several more simplified models can be constructed and compared using leave-one-out cross-validation (LOO) or based on the estimated widely applicable information criterion (WAIC). The covariates included in the model and the reasons are as follows:
Detection submodel
- Detection method: Hair-tubes and camera traps are different methods of detecting species presence and are expected to detect species at different probabilities.
Occupancy submodel
Baiting:Species: This is the key variable in the study, baiting was implemented at the four locations in different years and if fox baiting had an impact on occupancy of target species this should be reflected in the estimate. If fox baiting is successful we would expect a non-negative or positive estimate for baiting. The baiting variable is a continuous variable of how many years baiting at the site has been undertaken for as we expect a cumulative effect of baiting over years (e.g. the longer the baiting is conducted the larger the effect). The effect of baiting on occupancy is likely speces-specific; hence an interaction term is fitted
Elevation:Species: Elevation is a likely factor driving species abundance and occupancy for the four target species. The effect of elevation on species occupancy is expected to be dependent on the target species. Hence an interaction effect is modelled.
Year:Species: year was scaled so that 1 reflected 2008 and 10 was 2017. As such, the effect of year on occupancy will indicate the trend of each species over time.
FireFreq:Species: Frequency of recent fires may also impact occupancy in species-specific ways. The ecological fire group could also be used but FireFreq allows for a continuous scale to be used in the analysis.
6mthlag: This is the rainfall associated with each site for a six-month period and may impact species occupancy. This is a broad environmental variable that may signify a range of ecological responses to climate and environmental conditions.
(1|Location): A random-effect variable for location is implemented to control for and estimate the variability in occupancy patterns accross locations, as there may well be lrge amounts of variation between locations that cannot be explained by the model parameters implemented here.
(1|Site): Similar to the location random effect, the site random effect variable should explain variability between sites. Unfortunately
ubmsdoes not allow for nested random effects; ideally sites should be nested within location.
Global Model
The following model contains all covariates listed above
Show code
# cache to avoid re-running
if(!file.exists("models/fit_global.rds")) {
fit_global <- stan_occu(~method ~ scale(Elevation)*Species
+ scale(as.numeric(Baiting))*Species
+ scale(FireFreq)*Species
+ scale(`6mthlag`)*Species
+ (1|Location)
+ (1|Site),
data=umf_stack, chains=4, iter=2000)
saveRDS(fit_global, "models/fit_global.rds", compress = "xz")
} else {
fit_global <- readRDS("models/fit_global.rds")
}
Simplified Models
It is anticipated not all covariates have a strong relationship with occupancy of the target species. Based on the estimates of the global model several covariates were removed based on weak and/or variable effect sizes. Fire frquency and rainfall were omitted in the following model. These site/location specific covarites may have relationships with species occupancy, however could not be detected in the above model and are likely more complex than what the available data allows. Variation between sites and locations will still be estimated through the sites and location random effects; which may reflect a suite of ecological and environmental drivers of species occupancy (including fire, rainfall, temperature vegetation class).
Show code
# cache to avoid re-running
if(!file.exists("models/fit_simple1.rds")) {
fit_simple1 <- stan_occu(~method ~ Baiting
+ scale(Elevation)*Species
+ (1|Location)
+ (1|Site),
data=umf_stack, chains=4, iter=2000)
saveRDS(fit_simple1, "models/fit_simple1.rds", compress = "xz")
} else {
fit_simple1 <- readRDS("models/fit_simple1.rds")
}
# cache to avoid re-running
if(!file.exists("models/fit_simple2.rds")) {
fit_simple2 <- stan_occu(~method ~ scale(Elevation)*Species
+ scale(as.numeric(Baiting))*Species
+ (1|Site),
data=umf_stack, chains=4, iter=2000)
# fit_null <- stan_occu(~1~1, data=umf_stack, chains=6, iter=4000)
saveRDS(fit_simple2, "models/fit_simple2.rds", compress = "xz")
} else {
fit_simple2 <- readRDS("models/fit_simple2.rds")
}
if(!file.exists("models/fit_simple3.rds")) {
fit_simple3 <- stan_occu(~method ~ scale(as.numeric(Baiting))*Species
+ (1|Location)
+ (1|Site),
data=umf_stack, chains=4, iter=2000)
saveRDS(fit_simple3, "models/fit_simple3.rds", compress = "xz")
} else {
fit_simple3 <- readRDS("models/fit_simple3.rds")
}
if(!file.exists("models/fit_null.rds")) {
fit_null <- stan_occu(~1 ~ 1,
data=umf_stack, chains=4, iter=2000)
saveRDS(fit_null, "models/fit_null.rds", compress = "xz")
} else {
fit_null <- readRDS("models/fit_null.rds")
}
# cache to avoid re-running
if(!file.exists("models/fit_simple5.rds")) {
fit_simple5 <- stan_occu(~method ~ scale(Elevation)*Species
+ scale(as.numeric(Baiting))*Species*Location
+ scale(FireFreq)*Species
+ scale(`6mthlag`)*Species
+ (1|Site),
data=umf_stack, chains=4, iter=2000)
# fit_null <- stan_occu(~1~1, data=umf_stack, chains=6, iter=4000)
saveRDS(fit_simple5, "models/fit_simple5.rds", compress = "xz")
} else {
fit_simple5 <- readRDS("models/fit_simple5.rds")
}
Model Selection
The global and simplified models can be compared based on their predictive accuracy (elpd). Given that we are implementing a Bayesian approach to the occupancy models the computational costs of extensive model selection and dredging is high. Here we use model selection to find a suitable model that is well-founded in ecological theory from a select few well-reasoned models. The model selection table shows that the global model is the best performing model, with the next best performing model having a lower elpd. Given this, fit_simple5 will be the exclusive model used in further analyses.
Show code
mods <- fitList(fit_global, fit_simple1, fit_simple2, fit_simple3, fit_simple5, fit_null)
formulas <- lapply(mods@models, function(x) {
data.frame(formula = x@call[["formula"]] %>% capture.output() %>% trimws() %>% paste(., collapse = ""))
}) %>% bind_rows()
formulas$model <- names(mods)
modelSelection <- modSel(mods)
fit_top <- fit_simple5
round(modelSelection, 2) %>% tibble::rownames_to_column(var = "model") %>%
left_join(formulas) %>%
kableExtra::kbl("html") %>%
kable_styling()
| model | elpd | nparam | elpd_diff | se_diff | weight | formula |
|---|---|---|---|---|---|---|
| fit_simple5 | -4510.12 | 161.69 | 0.00 | 0.00 | 0.97 |
~method ~ scale(Elevation) * Species + scale(as.numeric(Baiting)) Species Location + scale(FireFreq) * Species + scale(6mthlag) *Species + (1 | Site)
|
| fit_global | -4638.91 | 137.52 | -128.79 | 16.16 | 0.00 |
~method ~ scale(Elevation) * Species + scale(as.numeric(Baiting)) Species + scale(FireFreq) Species + scale(6mthlag) *Species + (1 | Location) + (1 | Site)
|
| fit_simple2 | -4638.95 | 137.02 | -128.82 | 16.59 | 0.00 | ~method ~ scale(Elevation) * Species + scale(as.numeric(Baiting)) *Species + (1 | Site) |
| fit_simple1 | -4666.08 | 129.08 | -155.96 | 18.18 | 0.00 | ~method ~ Baiting + scale(Elevation) * Species + (1 | Location) +(1 | Site) |
| fit_simple3 | -4747.64 | 128.36 | -237.51 | 21.31 | 0.01 | ~method ~ scale(as.numeric(Baiting)) * Species + (1 | Location) +(1 | Site) |
| fit_null | -5470.25 | 3.11 | -960.13 | 40.70 | 0.02 | ~1 ~ 1 |
Model Diagnostics
The summary for the top model (global) is broken up into the Occupancy sub-model and the Detection sub-model. The estimates for the model paramters, standard deviation (SD), 95% confidence intervals, effective sample size and Rhat is presented below. n_eff > 200 and Rhat < 1.05 are general indicators the model has reached convergence.
Show code
fit_top
Call:
stan_occu(formula = ~method ~ scale(Elevation) * Species + scale(as.numeric(Baiting)) *
Species * Location + scale(FireFreq) * Species + scale(`6mthlag`) *
Species + (1 | Site), data = umf_stack, chains = 4, iter = 2000)
Occupancy (logit-scale):
Estimate SD
(Intercept) 0.41264 0.2334
scale(Elevation) 0.30731 0.1740
SpeciesLNB -5.23548 0.5984
SpeciesLNP -1.93652 0.2842
SpeciesSBB -2.22498 0.2954
scale(as.numeric(Baiting)) 0.67040 0.3071
LocationGenoa 1.32509 0.4924
LocationMueller -0.05302 0.4209
LocationTamboon -0.64394 0.4333
scale(FireFreq) 0.16534 0.1476
scale(`6mthlag`) 0.00269 0.0930
scale(Elevation):SpeciesLNB -0.29664 0.3002
scale(Elevation):SpeciesLNP -1.33357 0.2191
scale(Elevation):SpeciesSBB -0.44475 0.2190
SpeciesLNB:scale(as.numeric(Baiting)) -0.22856 0.4989
SpeciesLNP:scale(as.numeric(Baiting)) -1.15741 0.4143
SpeciesSBB:scale(as.numeric(Baiting)) -2.19429 0.5411
scale(as.numeric(Baiting)):LocationGenoa -0.07301 0.5555
scale(as.numeric(Baiting)):LocationMueller 1.46509 0.6311
scale(as.numeric(Baiting)):LocationTamboon -2.20607 0.4374
SpeciesLNB:LocationGenoa -0.04155 0.7727
SpeciesLNP:LocationGenoa -1.21664 0.5916
SpeciesSBB:LocationGenoa -2.53411 0.6381
SpeciesLNB:LocationMueller 3.30193 0.7379
SpeciesLNP:LocationMueller 2.30824 0.5033
SpeciesSBB:LocationMueller 0.09484 0.4938
SpeciesLNB:LocationTamboon -0.09797 0.9175
SpeciesLNP:LocationTamboon 0.35779 0.5300
SpeciesSBB:LocationTamboon -1.35950 0.7457
SpeciesLNB:scale(FireFreq) 0.23380 0.2431
SpeciesLNP:scale(FireFreq) 0.48616 0.1776
SpeciesSBB:scale(FireFreq) 0.31984 0.1833
SpeciesLNB:scale(`6mthlag`) 0.32491 0.2466
SpeciesLNP:scale(`6mthlag`) 0.17544 0.1416
SpeciesSBB:scale(`6mthlag`) -0.03986 0.1526
SpeciesLNB:scale(as.numeric(Baiting)):LocationGenoa 0.30629 0.7450
SpeciesLNP:scale(as.numeric(Baiting)):LocationGenoa 0.15503 0.7001
SpeciesSBB:scale(as.numeric(Baiting)):LocationGenoa 1.92481 0.7947
SpeciesLNB:scale(as.numeric(Baiting)):LocationMueller -1.13600 0.7853
SpeciesLNP:scale(as.numeric(Baiting)):LocationMueller -1.45725 0.7647
SpeciesSBB:scale(as.numeric(Baiting)):LocationMueller -0.80716 0.7897
SpeciesLNB:scale(as.numeric(Baiting)):LocationTamboon 2.73927 0.8365
SpeciesLNP:scale(as.numeric(Baiting)):LocationTamboon 0.63454 0.6313
SpeciesSBB:scale(as.numeric(Baiting)):LocationTamboon 1.14372 0.9749
sigma [1|Site] 1.05487 0.0968
2.5% 97.5%
(Intercept) -0.0437 0.8757
scale(Elevation) -0.0343 0.6468
SpeciesLNB -6.4851 -4.1577
SpeciesLNP -2.5108 -1.3978
SpeciesSBB -2.8333 -1.6748
scale(as.numeric(Baiting)) 0.1262 1.2936
LocationGenoa 0.4110 2.3573
LocationMueller -0.8583 0.7666
LocationTamboon -1.5013 0.2006
scale(FireFreq) -0.1229 0.4519
scale(`6mthlag`) -0.1799 0.1857
scale(Elevation):SpeciesLNB -0.8970 0.2904
scale(Elevation):SpeciesLNP -1.7701 -0.9113
scale(Elevation):SpeciesSBB -0.8710 -0.0125
SpeciesLNB:scale(as.numeric(Baiting)) -1.2429 0.6912
SpeciesLNP:scale(as.numeric(Baiting)) -2.0169 -0.3961
SpeciesSBB:scale(as.numeric(Baiting)) -3.2886 -1.1928
scale(as.numeric(Baiting)):LocationGenoa -1.0738 1.1221
scale(as.numeric(Baiting)):LocationMueller 0.2794 2.7538
scale(as.numeric(Baiting)):LocationTamboon -3.0840 -1.3936
SpeciesLNB:LocationGenoa -1.5307 1.5217
SpeciesLNP:LocationGenoa -2.4314 -0.0694
SpeciesSBB:LocationGenoa -3.8583 -1.3224
SpeciesLNB:LocationMueller 1.9217 4.8381
SpeciesLNP:LocationMueller 1.3722 3.3217
SpeciesSBB:LocationMueller -0.8442 1.0858
SpeciesLNB:LocationTamboon -1.9286 1.7092
SpeciesLNP:LocationTamboon -0.6878 1.4060
SpeciesSBB:LocationTamboon -2.9202 0.0116
SpeciesLNB:scale(FireFreq) -0.2250 0.7148
SpeciesLNP:scale(FireFreq) 0.1378 0.8350
SpeciesSBB:scale(FireFreq) -0.0345 0.6740
SpeciesLNB:scale(`6mthlag`) -0.1440 0.8441
SpeciesLNP:scale(`6mthlag`) -0.0974 0.4478
SpeciesSBB:scale(`6mthlag`) -0.3284 0.2631
SpeciesLNB:scale(as.numeric(Baiting)):LocationGenoa -1.2397 1.7502
SpeciesLNP:scale(as.numeric(Baiting)):LocationGenoa -1.3085 1.4529
SpeciesSBB:scale(as.numeric(Baiting)):LocationGenoa 0.3461 3.4796
SpeciesLNB:scale(as.numeric(Baiting)):LocationMueller -2.6866 0.4298
SpeciesLNP:scale(as.numeric(Baiting)):LocationMueller -2.9691 0.0746
SpeciesSBB:scale(as.numeric(Baiting)):LocationMueller -2.3968 0.7335
SpeciesLNB:scale(as.numeric(Baiting)):LocationTamboon 1.1365 4.4261
SpeciesLNP:scale(as.numeric(Baiting)):LocationTamboon -0.5951 1.8468
SpeciesSBB:scale(as.numeric(Baiting)):LocationTamboon -0.7427 3.0149
sigma [1|Site] 0.8768 1.2563
n_eff Rhat
(Intercept) 1114 1.01
scale(Elevation) 1507 1.00
SpeciesLNB 1420 1.00
SpeciesLNP 1702 1.00
SpeciesSBB 1691 1.00
scale(as.numeric(Baiting)) 1255 1.00
LocationGenoa 1027 1.00
LocationMueller 1072 1.00
LocationTamboon 1245 1.00
scale(FireFreq) 1426 1.00
scale(`6mthlag`) 2102 1.00
scale(Elevation):SpeciesLNB 2696 1.00
scale(Elevation):SpeciesLNP 2126 1.00
scale(Elevation):SpeciesSBB 2143 1.00
SpeciesLNB:scale(as.numeric(Baiting)) 1693 1.00
SpeciesLNP:scale(as.numeric(Baiting)) 1490 1.00
SpeciesSBB:scale(as.numeric(Baiting)) 1723 1.00
scale(as.numeric(Baiting)):LocationGenoa 1175 1.00
scale(as.numeric(Baiting)):LocationMueller 1200 1.00
scale(as.numeric(Baiting)):LocationTamboon 1528 1.00
SpeciesLNB:LocationGenoa 1790 1.00
SpeciesLNP:LocationGenoa 1376 1.00
SpeciesSBB:LocationGenoa 1461 1.00
SpeciesLNB:LocationMueller 1335 1.00
SpeciesLNP:LocationMueller 1471 1.00
SpeciesSBB:LocationMueller 1337 1.00
SpeciesLNB:LocationTamboon 1829 1.00
SpeciesLNP:LocationTamboon 1924 1.00
SpeciesSBB:LocationTamboon 1825 1.00
SpeciesLNB:scale(FireFreq) 2996 1.00
SpeciesLNP:scale(FireFreq) 2712 1.00
SpeciesSBB:scale(FireFreq) 2614 1.00
SpeciesLNB:scale(`6mthlag`) 3732 1.00
SpeciesLNP:scale(`6mthlag`) 2587 1.00
SpeciesSBB:scale(`6mthlag`) 2634 1.00
SpeciesLNB:scale(as.numeric(Baiting)):LocationGenoa 1488 1.00
SpeciesLNP:scale(as.numeric(Baiting)):LocationGenoa 1330 1.00
SpeciesSBB:scale(as.numeric(Baiting)):LocationGenoa 1729 1.00
SpeciesLNB:scale(as.numeric(Baiting)):LocationMueller 1574 1.00
SpeciesLNP:scale(as.numeric(Baiting)):LocationMueller 1301 1.00
SpeciesSBB:scale(as.numeric(Baiting)):LocationMueller 1500 1.00
SpeciesLNB:scale(as.numeric(Baiting)):LocationTamboon 2177 1.00
SpeciesLNP:scale(as.numeric(Baiting)):LocationTamboon 1702 1.00
SpeciesSBB:scale(as.numeric(Baiting)):LocationTamboon 2004 1.00
sigma [1|Site] 776 1.01
Detection (logit-scale):
Estimate SD 2.5% 97.5% n_eff Rhat
(Intercept) -1.79 0.0913 -1.97 -1.61 2141 1
methodht 1.76 0.0979 1.57 1.95 2157 1
LOOIC: 9020.248
Runtime: 10.747 minBased on the \(\sigma\) values for the random effects (shown above), there is an effect of site and location in explaining the overall variation. The estimates for the random effects can also be viewed and are as follows:
Show code
ran <- ranef(fit_top, submodel="state", summary=TRUE)
ran[[1]][[1]] %>%
kbl("html", digits = 3, caption = "Site random effects") %>%
kable_styling(full_width = FALSE)
| Estimate | SD | 2.5% | 97.5% | |
|---|---|---|---|---|
| D01 | 0.411 | 1.095 | -1.732 | 2.546 |
| D02 | 0.403 | 1.088 | -1.777 | 2.561 |
| D03 | 0.412 | 1.109 | -1.836 | 2.609 |
| D04 | 0.376 | 1.068 | -1.679 | 2.462 |
| D05 | 0.423 | 1.104 | -1.709 | 2.554 |
| D06 | 0.394 | 1.093 | -1.766 | 2.488 |
| D07 | 0.446 | 1.087 | -1.725 | 2.641 |
| D08 | 0.392 | 1.079 | -1.801 | 2.481 |
| D09 | 0.409 | 1.097 | -1.746 | 2.550 |
| D10 | -0.341 | 0.758 | -1.830 | 1.180 |
| D11 | 1.004 | 0.710 | -0.390 | 2.397 |
| D12 | 1.033 | 0.699 | -0.365 | 2.375 |
| D13 | -0.691 | 0.739 | -2.199 | 0.745 |
| D14 | -0.253 | 0.721 | -1.698 | 1.130 |
| D15 | -0.546 | 0.736 | -2.021 | 0.862 |
| D16 | 0.772 | 0.712 | -0.635 | 2.160 |
| D17 | -0.124 | 0.715 | -1.531 | 1.288 |
| D18 | 0.383 | 0.695 | -0.984 | 1.809 |
| D19 | -0.058 | 0.701 | -1.471 | 1.287 |
| D20 | -0.521 | 0.732 | -1.990 | 0.903 |
| D21 | -0.015 | 0.732 | -1.491 | 1.448 |
| D22 | 0.124 | 0.745 | -1.330 | 1.563 |
| D23 | 1.302 | 0.668 | -0.016 | 2.587 |
| D24 | -0.671 | 0.715 | -2.127 | 0.682 |
| D25 | -0.642 | 0.662 | -1.966 | 0.643 |
| D26 | -0.056 | 0.650 | -1.358 | 1.171 |
| D27 | 1.602 | 0.663 | 0.322 | 2.957 |
| D28 | 0.906 | 0.617 | -0.319 | 2.071 |
| D29 | -0.165 | 0.688 | -1.593 | 1.177 |
| D30 | 0.669 | 0.644 | -0.638 | 1.914 |
| D31 | 0.079 | 0.621 | -1.166 | 1.283 |
| D32 | -1.203 | 0.779 | -2.816 | 0.231 |
| D33 | -0.644 | 0.675 | -2.032 | 0.642 |
| D34 | 0.861 | 0.659 | -0.396 | 2.180 |
| D35 | 0.535 | 0.629 | -0.722 | 1.739 |
| D36 | 0.333 | 0.667 | -0.959 | 1.621 |
| D37 | 0.890 | 0.668 | -0.476 | 2.217 |
| D38 | 0.343 | 0.658 | -0.973 | 1.616 |
| D39 | 0.818 | 0.594 | -0.325 | 1.940 |
| D40 | 1.175 | 0.610 | -0.015 | 2.346 |
| D41 | -0.126 | 0.631 | -1.401 | 1.113 |
| D42 | 0.681 | 0.617 | -0.567 | 1.831 |
| D43 | 0.072 | 0.648 | -1.239 | 1.303 |
| D44 | 1.042 | 0.596 | -0.106 | 2.245 |
| D45 | 1.427 | 0.578 | 0.276 | 2.586 |
| D46 | 1.095 | 0.586 | -0.083 | 2.217 |
| D47 | 1.182 | 0.587 | 0.012 | 2.310 |
| D48 | 1.604 | 0.643 | 0.351 | 2.854 |
| D49 | 0.547 | 0.622 | -0.666 | 1.780 |
| D50 | 1.050 | 0.632 | -0.131 | 2.289 |
| D51 | 1.165 | 0.653 | -0.052 | 2.501 |
| D52 | 1.635 | 0.662 | 0.367 | 2.947 |
| D53 | -0.205 | 0.688 | -1.578 | 1.128 |
| D54 | 0.804 | 0.617 | -0.400 | 2.022 |
| D55 | -0.759 | 0.727 | -2.242 | 0.652 |
| D56 | 0.382 | 0.736 | -1.053 | 1.803 |
| D57 | 0.823 | 0.704 | -0.559 | 2.156 |
| D58 | 0.332 | 0.747 | -1.154 | 1.768 |
| D59 | -0.076 | 0.774 | -1.611 | 1.424 |
| D60 | 1.481 | 0.714 | 0.084 | 2.881 |
| G01 | 0.408 | 1.096 | -1.742 | 2.550 |
| G02 | 0.387 | 1.071 | -1.796 | 2.404 |
| G03 | 0.407 | 1.070 | -1.710 | 2.505 |
| G04 | 0.426 | 1.110 | -1.796 | 2.617 |
| G05 | 0.416 | 1.101 | -1.712 | 2.586 |
| G06 | 0.416 | 1.102 | -1.815 | 2.590 |
| G07 | 0.425 | 1.078 | -1.737 | 2.585 |
| G08 | 0.415 | 1.087 | -1.666 | 2.531 |
| G09 | 0.424 | 1.100 | -1.706 | 2.620 |
| G10 | 0.630 | 0.731 | -0.806 | 2.066 |
| G14 | -0.366 | 0.778 | -1.891 | 1.115 |
| G15 | -0.074 | 0.769 | -1.538 | 1.469 |
| G16 | 1.053 | 0.727 | -0.376 | 2.462 |
| G17 | 0.561 | 0.746 | -0.852 | 2.056 |
| G18 | 2.072 | 0.724 | 0.674 | 3.473 |
| G19 | 0.419 | 0.784 | -1.127 | 1.917 |
| G20 | 0.496 | 0.796 | -1.062 | 2.024 |
| G21 | 0.437 | 0.792 | -1.162 | 1.948 |
| G22 | -1.589 | 0.865 | -3.290 | 0.085 |
| G23 | 0.255 | 0.825 | -1.333 | 1.918 |
| G24 | 0.940 | 0.748 | -0.540 | 2.394 |
| G25 | 0.580 | 0.794 | -0.947 | 2.126 |
| G26 | 0.102 | 0.790 | -1.419 | 1.629 |
| G27 | 0.436 | 0.783 | -1.090 | 1.941 |
| G28 | 0.433 | 0.769 | -1.138 | 1.932 |
| G29 | 0.467 | 0.746 | -1.005 | 1.916 |
| G30 | 0.932 | 0.746 | -0.545 | 2.393 |
| G31 | 0.312 | 0.756 | -1.224 | 1.776 |
| G32 | 0.767 | 0.706 | -0.640 | 2.104 |
| G33 | 0.123 | 0.732 | -1.338 | 1.531 |
| G34 | 0.857 | 0.692 | -0.518 | 2.196 |
| G35 | 0.496 | 0.780 | -1.036 | 1.995 |
| G36 | 0.701 | 0.865 | -0.956 | 2.399 |
| G37 | 1.012 | 0.741 | -0.438 | 2.443 |
| G38 | 0.025 | 0.807 | -1.582 | 1.618 |
| G39 | 0.101 | 0.754 | -1.374 | 1.553 |
| G40 | 0.444 | 0.779 | -1.071 | 1.984 |
| G41 | 0.147 | 0.726 | -1.282 | 1.559 |
| G42 | 0.179 | 0.770 | -1.317 | 1.685 |
| G43 | 0.696 | 0.791 | -0.837 | 2.298 |
| G44 | -0.012 | 0.749 | -1.432 | 1.513 |
| G45 | 0.445 | 0.766 | -1.047 | 1.979 |
| G46 | 1.012 | 0.711 | -0.366 | 2.443 |
| G47 | 0.222 | 0.717 | -1.194 | 1.618 |
| G48 | 0.827 | 0.704 | -0.559 | 2.176 |
| G49 | 0.451 | 0.699 | -0.929 | 1.871 |
| G50 | 0.872 | 0.677 | -0.426 | 2.188 |
| G51 | 0.130 | 0.825 | -1.434 | 1.749 |
| G52 | 0.234 | 0.742 | -1.240 | 1.693 |
| G53 | 0.236 | 0.801 | -1.257 | 1.779 |
| G54 | -0.859 | 0.814 | -2.449 | 0.770 |
| G55 | 0.138 | 0.782 | -1.337 | 1.702 |
| G56 | 0.157 | 0.798 | -1.397 | 1.734 |
| G57 | 1.655 | 0.740 | 0.214 | 3.131 |
| G58 | 0.434 | 0.772 | -1.039 | 1.962 |
| G59 | 0.406 | 0.757 | -1.065 | 1.858 |
| G60 | 0.430 | 0.728 | -0.984 | 1.840 |
| M01 | 0.406 | 1.071 | -1.755 | 2.508 |
| M02 | 0.426 | 1.073 | -1.676 | 2.505 |
| M03 | 0.418 | 1.077 | -1.644 | 2.573 |
| M04 | 0.402 | 1.079 | -1.728 | 2.459 |
| M05 | 0.411 | 1.086 | -1.712 | 2.579 |
| M06 | 0.407 | 1.103 | -1.709 | 2.552 |
| M07 | 0.424 | 1.108 | -1.760 | 2.543 |
| M08 | 0.410 | 1.078 | -1.666 | 2.487 |
| M09 | 0.402 | 1.079 | -1.710 | 2.524 |
| M10 | -0.455 | 0.746 | -1.967 | 0.992 |
| M11 | 0.509 | 0.656 | -0.770 | 1.784 |
| M12 | 1.847 | 0.616 | 0.687 | 3.090 |
| M13 | 0.656 | 0.598 | -0.497 | 1.827 |
| M14 | 0.882 | 0.594 | -0.299 | 2.027 |
| M15 | 1.867 | 0.607 | 0.720 | 3.051 |
| M16 | 1.066 | 0.586 | -0.128 | 2.201 |
| M17 | 1.475 | 0.602 | 0.319 | 2.669 |
| M18 | 1.053 | 0.607 | -0.107 | 2.232 |
| M19 | 0.819 | 0.597 | -0.337 | 2.001 |
| M20 | 0.211 | 0.616 | -0.983 | 1.430 |
| M21 | 1.196 | 0.603 | 0.018 | 2.369 |
| M22 | 0.803 | 0.613 | -0.368 | 2.001 |
| M23 | 0.646 | 0.600 | -0.513 | 1.851 |
| M24 | 1.084 | 0.607 | -0.092 | 2.277 |
| M25 | 0.170 | 0.629 | -1.059 | 1.413 |
| M26 | 0.398 | 0.646 | -0.890 | 1.666 |
| M27 | 0.014 | 0.675 | -1.320 | 1.302 |
| M28 | 0.572 | 0.610 | -0.641 | 1.753 |
| M29 | 0.301 | 0.661 | -1.023 | 1.578 |
| M30 | 0.217 | 0.658 | -1.081 | 1.469 |
| M31 | -0.622 | 0.750 | -2.144 | 0.822 |
| M32 | 0.913 | 0.605 | -0.278 | 2.119 |
| M33 | 0.703 | 0.639 | -0.570 | 1.984 |
| M34 | -0.525 | 0.636 | -1.787 | 0.678 |
| M35 | -0.568 | 0.666 | -1.893 | 0.708 |
| M36 | -0.257 | 0.673 | -1.607 | 1.029 |
| M37 | -0.312 | 0.648 | -1.566 | 0.953 |
| M38 | 0.420 | 0.631 | -0.819 | 1.621 |
| M39 | 0.235 | 0.629 | -1.037 | 1.443 |
| M40 | 1.293 | 0.620 | 0.090 | 2.495 |
| M41 | -0.729 | 0.704 | -2.146 | 0.573 |
| M42 | -0.690 | 0.662 | -2.054 | 0.583 |
| M43 | 0.143 | 0.652 | -1.161 | 1.398 |
| M44 | -0.347 | 0.654 | -1.628 | 0.918 |
| M45 | 0.852 | 0.620 | -0.390 | 2.059 |
| M46 | 0.296 | 0.639 | -0.983 | 1.509 |
| M47 | 0.206 | 0.616 | -1.051 | 1.404 |
| M48 | 0.532 | 0.613 | -0.744 | 1.735 |
| M49 | 0.289 | 0.639 | -0.954 | 1.531 |
| M50 | 0.411 | 0.628 | -0.822 | 1.638 |
| M51 | 0.543 | 0.626 | -0.697 | 1.739 |
| M52 | -0.405 | 0.693 | -1.832 | 0.889 |
| M53 | 0.282 | 0.658 | -0.975 | 1.571 |
| M54 | 0.407 | 0.682 | -0.944 | 1.766 |
| M55 | -1.311 | 0.783 | -2.815 | 0.205 |
| M56 | 0.639 | 0.773 | -0.864 | 2.184 |
| M57 | 0.467 | 0.732 | -0.983 | 1.886 |
| M58 | 0.553 | 0.603 | -0.646 | 1.721 |
| M59 | 1.047 | 0.572 | -0.061 | 2.154 |
| M60 | 1.522 | 0.579 | 0.393 | 2.659 |
| T01 | 0.409 | 1.063 | -1.657 | 2.544 |
| T02 | 0.422 | 1.057 | -1.713 | 2.549 |
| T03 | 0.411 | 1.094 | -1.750 | 2.493 |
| T04 | 0.409 | 1.091 | -1.753 | 2.590 |
| T05 | 0.417 | 1.079 | -1.667 | 2.522 |
| T06 | 0.409 | 1.081 | -1.730 | 2.518 |
| T07 | 0.407 | 1.050 | -1.675 | 2.430 |
| T08 | 0.434 | 1.065 | -1.648 | 2.600 |
| T09 | 0.412 | 1.089 | -1.674 | 2.529 |
| T10 | -0.977 | 0.703 | -2.401 | 0.381 |
| T11 | 0.262 | 0.710 | -1.067 | 1.710 |
| T12 | -1.772 | 0.777 | -3.349 | -0.340 |
| T13 | -1.701 | 0.740 | -3.235 | -0.326 |
| T14 | -1.253 | 0.776 | -2.832 | 0.248 |
| T15 | 2.850 | 0.736 | 1.414 | 4.347 |
| T16 | -1.334 | 0.798 | -2.943 | 0.178 |
| T17 | 0.395 | 0.734 | -1.091 | 1.811 |
| T18 | -1.563 | 0.765 | -3.151 | -0.145 |
| T19 | -0.744 | 0.686 | -2.111 | 0.597 |
| T20 | -0.862 | 0.670 | -2.178 | 0.400 |
| T21 | 1.078 | 0.656 | -0.163 | 2.401 |
| T22 | 1.171 | 0.668 | -0.083 | 2.550 |
| T23 | 1.108 | 0.626 | -0.083 | 2.372 |
| T24 | 1.036 | 0.650 | -0.219 | 2.316 |
| T25 | 1.677 | 0.674 | 0.423 | 3.039 |
| T26 | 0.440 | 0.640 | -0.798 | 1.683 |
| T27 | -0.864 | 0.686 | -2.230 | 0.436 |
| T28 | 1.653 | 0.664 | 0.419 | 2.987 |
| T29 | 0.205 | 0.645 | -1.051 | 1.472 |
| T30 | 1.423 | 0.688 | 0.122 | 2.791 |
| T31 | 0.744 | 0.640 | -0.469 | 1.996 |
| T32 | 1.395 | 0.655 | 0.166 | 2.708 |
| T33 | 1.008 | 0.628 | -0.191 | 2.261 |
| T34 | 0.748 | 0.626 | -0.451 | 2.001 |
| T35 | 0.236 | 0.652 | -1.012 | 1.530 |
| T36 | 1.151 | 0.647 | -0.070 | 2.477 |
| T37 | 0.693 | 0.677 | -0.609 | 2.007 |
| T38 | 2.666 | 0.701 | 1.392 | 4.127 |
| T39 | 0.584 | 0.670 | -0.705 | 1.936 |
| T40 | 2.026 | 0.681 | 0.728 | 3.444 |
| T41 | 0.658 | 0.718 | -0.710 | 2.116 |
| T42 | 1.380 | 0.671 | 0.076 | 2.686 |
| T43 | -0.042 | 0.651 | -1.305 | 1.264 |
| T44 | -0.344 | 0.664 | -1.606 | 0.944 |
| T45 | -1.595 | 0.774 | -3.175 | -0.132 |
| T46 | -0.969 | 0.662 | -2.287 | 0.324 |
| T47 | -1.574 | 0.699 | -2.975 | -0.282 |
| T48 | 0.241 | 0.640 | -0.977 | 1.505 |
| T49 | -0.998 | 0.667 | -2.311 | 0.323 |
| T50 | -1.529 | 0.755 | -3.070 | -0.126 |
| T51 | 1.598 | 0.668 | 0.363 | 2.964 |
| T52 | 0.995 | 0.645 | -0.250 | 2.256 |
| T53 | 0.512 | 0.697 | -0.828 | 1.928 |
| T54 | 1.030 | 0.672 | -0.231 | 2.368 |
| T55 | 0.942 | 0.654 | -0.316 | 2.231 |
| T56 | 1.624 | 0.666 | 0.378 | 2.976 |
| T57 | 1.008 | 0.630 | -0.204 | 2.262 |
| T58 | 1.701 | 0.668 | 0.431 | 3.029 |
| T59 | 1.088 | 0.648 | -0.174 | 2.363 |
| T60 | 1.166 | 0.654 | -0.030 | 2.487 |
Given that we implemented a Bayesian approach here the following section presents several model diagnostics used to ensure model convergence and suitability.
MCMC Convergence and Diagnostics
The traceplots how that the chains have mixed well and the number of iterations used in the model is sufficient in reaching convergence.
Show code
traceplot(fit_top, pars=c("beta_state", "beta_det"), ncol = 3)
Figure 1: Traceplots of the four chains run for 2000 iterations in the global model. Plots show good evidence that chain-mixing and convergence has been achieved
Model Fit
Residuals
The residuals for the occupancy model can be calculated and plotted. The ubms package implements a method from Wright, Irvine, and Higgs (2019) to calculate residuals, where residuals are calculated independently for the observation and detection submodels. The diagnostic information from visualising these plots can be summarised as follows: if 95 % of the binned residual points fall within the grey shading then it is an indication the model fits the data well (Kellner 2021).
Show code
plot_residuals(fit_top, submodel="state")
Figure 2: Residual plot for the global model, showing that binned residuals largely fall within the margins of error a model with good predictive power would have
Posterior Predictive Tests
We can predict the posterior distribution of the model and compare those predictions to the actual data. Here we show that the posterior prediction distribution matching the real data at a generalized level.
Show code
sim_y <- posterior_predict(fit_top, "y", draws=500)
prop0 <- apply(sim_y, 1, function(x) mean(x==0, na.rm=TRUE))
prop1 <- apply(sim_y, 1, function(x) mean(x==1, na.rm=TRUE))
actual_prop0 <- mean(getY(fit_top) == 0, na.rm=TRUE)
actual_prop1 <- mean(getY(fit_top) == 1, na.rm=TRUE)
#Compare 0
par(mfrow = c(1,2))
hist(prop0, col='gray')
abline(v=actual_prop0, col='red', lwd=2)
#Compare 1
hist(prop1, col='gray')
abline(v=actual_prop1, col='red', lwd=2)
Figure 3: Comparison of posterior predictions from the global model to true mean values of absences and presences (0 and 1)
Model Results
Marginal Covariate Effects for Occupancy
Showing the generalised marginal effects of the covariates on the probability of occupancy is difficult for the model used. This is primarily because of interaction terms used in the model. The interaction terms ensure that estimates and predictions are generated for each species. Using the ubms::plot_marginal() function does not take these interactions into effect, hence the marginal effects shown are only for the reference level of species (Brushtail Possum). In order to more accurately visualise these effects we used model predictions; to see these ouputs, navigate to the ‘Predictions’ section below.
Show code
plot_marginal(fit_top, "state")
Marginal Covariate Effects for Detection
Comparing the methods of detection is a bit easier because of the simple submodel process for determining the probability of detection (solely dependent on method). Based on the results from the detection submodel we see that the hair tubes appear to have a higher probability of detection than camera traps.
Show code
detection_predictions <- ubms::posterior_linpred(fit_top,
submodel="det",
transform = TRUE,
draws = 1000,
newdata=data.frame(method = c("cam", "ht"))) %>%
as.data.frame() %>%
`colnames<-`(c("Camera", "Hair Tube")) %>% melt()
ggplot(data = detection_predictions,
aes(x = variable, y = value)) +
tidybayes::geom_eye() +
ylim(0,1) +
xlab("Detection method") +
ylab("Probability of detection (0-1)") +
theme_minimal()
Figure 4: Posterior predictions (grey shading) for the probability of detection using camera traps and hair tubes. Hair tubes show a higher detection probability than camera traps. Mean values from 1,000 draws of the posterior are shown with the black dot, lines represent 95 % CI’s.
Show code
# plot_marginal(fit_top, "det") + theme_minimal()
Predictions
Sometimes marginal estimates can be difficult to interpret in isolation because of factor levels and various interactions. Using model predictions can allow for the model to be projected to a ‘newdata’ set that matches the covariates of interest. Here we use model predictions to visualise the impact of baiting on the probability of occupancy for the four target species across the four areas of interest
Key Comparisons to make based on predictions:
- All species occupancy over time (years after baiting), for this prediction the following parameters are used:
- mean elevation
- mean fire frequency
- mean rainfall
- All locations and sites
- mean elevation
- If baiting did not occur vs baiting occurred: Hypothetical showing the predicted effect on occupancy if baiting did not occur. For this prediction the following parameters are used:
- mean elevation
- mean fire frequency
- mean rainfall
- All locations and average of sites within a location
- mean elevation
Additional Marginal Effects on a Species Basis. In these cases predictions are made across all locations for each species with other covariate values set to the mean:
Fire Frequency
Elevation
- Rainfall (6mthlag)
Show code
newdata <- umf_stack@siteCovs %>%
# filter(Location != "Tamboon") %>%
group_by(Site, Location, Species) %>%
summarise(Elevation = mean(Elevation),
FireFreq = mean(FireFreq),
`6mthlag` = mean(`6mthlag`),
Year = mean(Year))
baiting_newdata <- sapply(newdata, rep.int, times=11, simplify = FALSE)
baiting_newdata$Baiting <- rep(c(0L:10L), each = nrow(newdata))
baiting_newdata <- bind_cols(baiting_newdata)
baiting_predictions <- predict(fit_simple5, submodel="state", newdata=baiting_newdata, level = c(0.5, 0.95))
baiting_predictions <- bind_cols(baiting_newdata, baiting_predictions) %>%
group_by(Location, Baiting, Species) %>%
mutate(av_pred = mean(Predicted),
av_LCI_95 = mean(`2.5%`),
av_UCI_95 = mean(`97.5%`),
av_LCI_50 = mean(`25%`),
av_UCI_50 = mean(`75%`)) %>%
ungroup()
Show code
bait_prob <- baiting_predictions %>%
# filter(Location == "Mueller") %>%
dplyr::select(Species, av_pred, av_LCI_95, av_UCI_95, Baiting, Location) %>%
distinct() %>%
group_by(Species, Location) %>%
mutate(av_change = av_pred - lag(av_pred),
av_changeLCI = av_LCI_95 - lag(av_LCI_95),
av_change_UCI = av_UCI_95 - lag(av_UCI_95))
# group_by(Species) %>%
# summarise(av_change = mean(av_change, na.rm = TRUE),
# av_change_LCI = mean(av_LCI, na.rm = TRUE),
# av_change_UCI = mean(av_UCI, na.rm = TRUE))
# bait_prob %>%
# filter(Baiting %in% c(0,10)) %>%
# select(1:6) %>%
# as.data.table() %>%
# dcast(Species + Location ~ Baiting, value.var = c("av_pred", "av_LCI_95", "av_UCI_95")) %>%
# mutate(Start = paste0(round(`av_pred_0`, 2), " (", round(`av_LCI_95_0`, 2), ", ", round(`av_UCI_95_0`, 2), ")"),
# End = paste0(round(`av_pred_10`, 2), " (", round(`av_LCI_95_10`, 2), ", ", round(`av_UCI_95_10`, 2), ")"),
# diff = round(av_pred_10 - av_pred_0,2)) %>%
# select(Location, Species, Start, End, diff) %>%
# arrange(Location) %>%
# kbl(digits = 3)
ggplot(data = baiting_predictions %>% mutate(Baiting = as.integer(Baiting)),
aes(x = Baiting, y = Predicted, fill = Location)) +
# geom_ribbon(aes(ymin = av_LCI_50, ymax = av_UCI_50), alpha = 0.3) +
geom_ribbon(aes(ymin = av_LCI_95, ymax = av_UCI_95), alpha = 0.5) +
# geom_quasirandom(shape = 21, alpha = 0.3) +
geom_line(aes(y = av_pred, colour = Location), lwd = 1) +
# annotate("rect", xmin = 5, xmax = 9, ymin = 0, ymax = 1, fill = "grey60", alpha = 0.5) +
scale_fill_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
scale_color_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
scale_x_continuous(limits = c(-0.5, 10.5), breaks = 0:10) +
ylim(0,1) +
xlab("Years after baiting began") +
ylab("Model predictions of occupancy probability (0–1)") +
facet_wrap(~Species+Location, ncol = 4, nrow = 4) +
theme_bw()
Figure 5: Effect of Baiting over Time on Site Occupancy: Model predictions for the site occupancy of four species: Brushtail Possums (BTP), Long-nosed Bandicoots (LNB), Long-nosed Potoroos (LNP), and Southern Brown Bandicoots (SBB). Model predictions are plotted for each of the four study locations (60 sites at each location). Shading represents 95 % confidence intervals
Show code
newdata_year <- umf_stack@siteCovs %>%
# filter(Location != "Tamboon") %>%
group_by(Site, Location, Species, Year) %>%
summarise(Elevation = mean(Elevation),
FireFreq = mean(FireFreq),
`6mthlag` = mean(`6mthlag`))
no_baiting <- newdata_year %>%
mutate(Baiting = 0L)
no_baiting_predictions <- bind_cols(no_baiting, predict(fit_top, submodel="state", newdata=no_baiting)) %>%
group_by(Location, Year, Species) %>%
mutate(av_pred = mean(Predicted),
av_LCI = mean(`2.5%`),
av_UCI = mean(`97.5%`)) %>%
ungroup() %>%
mutate(Baited = "Never Baited")
yes_baiting <- umf_stack@siteCovs %>%
# filter(Location != "Tamboon") %>%
group_by(Site, Location, Species, Year, Baiting) %>%
summarise(Elevation = mean(Elevation),
FireFreq = mean(FireFreq),
`6mthlag` = mean(`6mthlag`))
yes_baiting_predictions <- bind_cols(yes_baiting, predict(fit_top, submodel="state", newdata=yes_baiting)) %>%
group_by(Location, Year, Species) %>%
mutate(av_pred = mean(Predicted),
av_LCI = mean(`2.5%`),
av_UCI = mean(`97.5%`)) %>%
ungroup() %>% mutate(Baited = "As Conducted")
baiting_comparison <- bind_rows(no_baiting_predictions, yes_baiting_predictions)
# baiting implemented as planned
when_baiting_started <- umf_stack@siteCovs %>%
# filter(Location != "Tamboon") %>%
group_by(Location) %>%
filter(Baiting > 0) %>%
summarise(Baiting_Start = min(Year)+2007) %>%
ungroup()
real_data <- umf_stack@siteCovs %>% bind_cols(as.data.frame(umf_stack@y)) %>%
group_by(Year, Location, Species) %>%
summarise(Occupancy_Prob = mean(c(obs_1, obs_2, obs_3, obs_4), na.rm = TRUE)) %>%
ungroup()
Show code
ggplot(data = baiting_comparison %>% mutate(Year = Year + 2007L),
aes(x = Year, y = Predicted, fill = Baited, colour = Baited, linetype = Baited)) +
geom_vline(aes(xintercept = Baiting_Start), data = when_baiting_started, colour = "DarkGreen") +
# geom_ribbon(aes(ymin = 0, ymax = Baiting, x = Year), inherit.aes = FALSE, fill = "grey90", alpha = 0.5) +
geom_ribbon(aes(ymin = av_LCI, ymax = av_UCI), alpha = 0.3) +
# geom_quasirandom(shape = 21, alpha = 0.3) +
geom_line(aes(y = av_pred), lwd = 1) +
# geom_point(data = real_data %>% mutate(Year = Year + 2007L),
# aes(x = Year, y = Occupancy_Prob), shape = 23, fill = "grey", inherit.aes = FALSE, size = 2) +
# annotate("rect", xmin = 5, xmax = 9, ymin = 0, ymax = 1, fill = "grey60", alpha = 0.5) +
scale_fill_brewer(palette = "Set1") +
scale_color_brewer(palette = "Set1") +
scale_x_continuous(limits = c(2008, 2017), n.breaks = 10, labels = c("08", "09", "10", "11", "12", "13", "14", "15", "16", "17")) +
ylim(0,1) +
xlab("Year") +
ylab("Model predictions of occupancy probability (0–1)") +
facet_wrap(~Species+Location, ncol = 4, nrow = 4) +
theme_minimal()
Figure 6: Comparison of the Conducted Baiting Regime to a Hypothetical Scenario where no Baiting Occured: Model predictions for the site occupancy of four species: Brushtail Possums (BTP), Long-nosed Bandicoots (LNB), Long-nosed Potoroos (LNP), and Southern Brown Bandicoots (SBB). The two predicted scenarios show the predicted occupancy rate for the ‘as conducted’ baiting regime (red solid line), where the green vertical line indicates the baiting start date. In comparison the alternative no baiting scenario (blue dashed line) are based on model predictions where baiting was fixed at zero over the years. 95% CI are shown with shading for each respective scenario.
Show code
plot_predicted_marginal <- function(model, model_stack, param, groupings, other_covs, breaks = 10) {
# range
if(class(model_stack@siteCovs[[param]]) %in% c("numeric", "integer")) {
min_param <- min(model_stack@siteCovs[[param]])
max_param <- max(model_stack@siteCovs[[param]])
}
seq_el <- seq.int(from = min_param, to = max_param, length.out = breaks)
newdata_r <- umf_stack@siteCovs %>%
group_by_at(groupings) %>%
summarise_at(.vars = other_covs, .funs = mean)
newdata <- sapply(newdata_r, rep.int, times=length(seq_el), simplify = FALSE)
newdata[[param]] <- rep(seq_el, each = nrow(newdata_r)) %>% as.numeric()
final_predictions <- bind_cols(newdata, predict(model,
submodel="state",
newdata=newdata)) %>%
group_by_at(c("Species", param)) %>%
mutate(av_pred = mean(Predicted),
av_LCI = mean(`2.5%`),
av_UCI = mean(`97.5%`)) %>%
ungroup()
ggplot(data = final_predictions,
aes(x = !!rlang::sym(param), y = Predicted)) +
geom_ribbon(aes(ymin = av_LCI, ymax = av_UCI), alpha = 0.3) +
geom_quasirandom(shape = 21, alpha = 0.3) +
geom_line(aes(y = av_pred), lwd = 1) +
# annotate("rect", xmin = 5, xmax = 9, ymin = 0, ymax = 1, fill = "grey60", alpha = 0.5) +
scale_fill_brewer(palette = "Set2") +
scale_color_brewer(palette = "Set2") +
scale_x_continuous() +
ylim(0,1) +
xlab(param) +
ylab("Model Predictions of Occupancy Probability (0-1)") +
facet_wrap(~Species, ncol = 1, nrow = 4) +
theme_minimal()
}
Show code
elev_data <- plot_predicted_marginal(model = fit_top,
model_stack = umf_stack,
param = "Elevation",
groupings = c("Species", "Location", "Site"),
other_covs = c("FireFreq", "Year", "Baiting", "6mthlag"), breaks = 10, return_data = TRUE) %>%
group_by(Species, Elevation) %>%
mutate(av_pred = mean(Predicted),
av_LCI = mean(`2.5%`),
av_UCI = mean(`97.5%`)) %>%
ungroup() %>%
group_by(Species, Elevation, Location) %>%
mutate(av_pred_loc = mean(Predicted),
av_LCI_loc = mean(`2.5%`),
av_UCI_loc = mean(`97.5%`)) %>%
ungroup()
elev_plot <- elev_data %>%
ggplot(aes(x = Elevation, y = Predicted)) +
geom_ribbon(aes(ymin = av_LCI, ymax = av_UCI, x = Elevation), inherit.aes = FALSE, alpha = 0.3) +
# ggbeeswarm::geom_quasirandom(aes(fill = Location, x = Elevation, y = Predicted), inherit.aes = FALSE, shape = 21, alpha = 0.3) +
geom_line(aes(y = av_pred, x = Elevation),inherit.aes = FALSE, lwd = 2) +
geom_line(aes(y = av_pred_loc, x = Elevation, colour = Location),inherit.aes = FALSE, lwd = 1, linetype = "dashed") +
# annotate("rect", xmin = 5, xmax = 9, ymin = 0, ymax = 1, fill = "grey60", alpha = 0.5) +
# scale_fill_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
scale_color_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
scale_x_continuous() +
ylim(0,1) +
xlab("Elevation") +
ylab("Model Predictions of Occupancy Probability (0-1)") +
facet_wrap(facets = c("Species"), ncol = 1) +
theme_bw() +
theme(legend.position = "none")
fire_data <- plot_predicted_marginal(model = fit_top,
model_stack = umf_stack,
param = "FireFreq",
groupings = c("Species", "Location", "Site"),
other_covs = c("Elevation", "Year", "Baiting", "6mthlag"), breaks = 10, return_data = TRUE) %>%
group_by(Species, FireFreq) %>%
mutate(av_pred = mean(Predicted),
av_LCI = mean(`2.5%`),
av_UCI = mean(`97.5%`)) %>%
ungroup() %>%
group_by(Species, FireFreq, Location) %>%
mutate(av_pred_loc = mean(Predicted),
av_LCI_loc = mean(`2.5%`),
av_UCI_loc = mean(`97.5%`)) %>%
ungroup()
fire_plot <- fire_data %>%
ggplot(aes(x = FireFreq, y = Predicted)) +
geom_ribbon(aes(ymin = av_LCI, ymax = av_UCI, x = FireFreq), inherit.aes = FALSE, alpha = 0.3) +
# ggbeeswarm::geom_quasirandom(aes(fill = Location, x = Elevation, y = Predicted), inherit.aes = FALSE, shape = 21, alpha = 0.3) +
geom_line(aes(y = av_pred, x = FireFreq),inherit.aes = FALSE, lwd = 2) +
geom_line(aes(y = av_pred_loc, x = FireFreq, colour = Location),inherit.aes = FALSE, lwd = 1, linetype = "dashed") +
# annotate("rect", xmin = 5, xmax = 9, ymin = 0, ymax = 1, fill = "grey60", alpha = 0.5) +
# scale_fill_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
scale_color_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
scale_x_continuous() +
ylim(0,1) +
xlab("FireFreq") +
ylab("Model Predictions of Occupancy Probability (0-1)") +
facet_wrap(facets = c("Species"), ncol = 1) +
theme_bw() +
theme(legend.position = "none")
# year_data <- plot_predicted_marginal(model = fit_top,
# model_stack = umf_stack,
# param = "Year",
# groupings = c("Species", "Location", "Site"),
# other_covs = c("Elevation", "FireFreq", "Baiting", "6mthlag"), breaks = 10, return_data = TRUE) %>%
# group_by(Species, Year) %>%
# mutate(av_pred = mean(Predicted),
# av_LCI = mean(`2.5%`),
# av_UCI = mean(`97.5%`)) %>%
# ungroup() %>%
# group_by(Species, Year, Location) %>%
# mutate(av_pred_loc = mean(Predicted),
# av_LCI_loc = mean(`2.5%`),
# av_UCI_loc = mean(`97.5%`)) %>%
# ungroup()
#
# year_plot <- year_data %>%
# ggplot(aes(x = Year, y = Predicted)) +
# geom_ribbon(aes(ymin = av_LCI, ymax = av_UCI, x = Year), inherit.aes = FALSE, alpha = 0.3) +
# # ggbeeswarm::geom_quasirandom(aes(fill = Location, x = Elevation, y = Predicted), inherit.aes = FALSE, shape = 21, alpha = 0.3) +
# geom_line(aes(y = av_pred, x = Year),inherit.aes = FALSE, lwd = 2) +
# geom_line(aes(y = av_pred_loc, x = Year, colour = Location),inherit.aes = FALSE, lwd = 1, linetype = "dashed") +
# # annotate("rect", xmin = 5, xmax = 9, ymin = 0, ymax = 1, fill = "grey60", alpha = 0.5) +
# # scale_fill_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
# scale_color_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
# scale_x_continuous() +
# ylim(0,1) +
# xlab("Year") +
# ylab("Model Predictions of Occupancy Probability (0-1)") +
# facet_wrap(facets = c("Species"), ncol = 1) +
# theme_bw() +
# theme(legend.position = "none")
Rainfall_data <- plot_predicted_marginal(model = fit_top,
model_stack = umf_stack,
param = "6mthlag",
groupings = c("Species", "Location", "Site"),
other_covs = c("Elevation", "FireFreq", "Baiting", "Year"), breaks = 10, return_data = TRUE) %>%
group_by(Species, `6mthlag`) %>%
mutate(av_pred = mean(Predicted),
av_LCI = mean(`2.5%`),
av_UCI = mean(`97.5%`)) %>%
ungroup() %>%
group_by(Species, `6mthlag`, Location) %>%
mutate(av_pred_loc = mean(Predicted),
av_LCI_loc = mean(`2.5%`),
av_UCI_loc = mean(`97.5%`)) %>%
ungroup()
Rainfall_plot <- Rainfall_data %>%
ggplot(aes(x = `6mthlag`, y = Predicted)) +
geom_ribbon(aes(ymin = av_LCI, ymax = av_UCI, x = `6mthlag`), inherit.aes = FALSE, alpha = 0.3) +
# ggbeeswarm::geom_quasirandom(aes(fill = Location, x = Elevation, y = Predicted), inherit.aes = FALSE, shape = 21, alpha = 0.3) +
geom_line(aes(y = av_pred, x = `6mthlag`),inherit.aes = FALSE, lwd = 2) +
geom_line(aes(y = av_pred_loc, x = `6mthlag`, colour = Location),inherit.aes = FALSE, lwd = 1, linetype = "dashed") +
# annotate("rect", xmin = 5, xmax = 9, ymin = 0, ymax = 1, fill = "grey60", alpha = 0.5) +
# scale_fill_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
scale_color_manual(values = delwp_cols[c(1:3, 5)] %>% `names<-`(unique(baiting_predictions$Location))) +
scale_x_continuous() +
ylim(0,1) +
xlab("Rainfall") +
ylab("Model Predictions of Occupancy Probability (0-1)") +
facet_wrap(facets = c("Species"), ncol = 1) +
theme_bw()
cowplot::plot_grid(elev_plot, fire_plot, Rainfall_plot, nrow = 1, labels = "AUTO", rel_widths = c(0.8,0.8,1.1))
Figure 7: Effect of Covariates on Site Occupancy: Model predictions for the site occupancy of four species: Brushtail Possums (BTP), Long-nosed Bandicoots (LNB), Long-nosed Potoroos (LNP), and Southern Brown Bandicoots (SBB). Model predictions are plotted for each species with individual points showing the prediction for a site at that location (60 sites at each location). The four panes (A-D) show the effects of Elevation (A), Fire frequency (B), Year (C) and 6-month lag Rainfall (D). Across these four model covariates, the probability of occupancy varies across species.
Summary of Results
Show code
model_summary <- summary(fit_top, "state") %>% round(3)
baiting_effects <- conditional_effects_ubms(fit_top,
submodel = "state",
categoricalParam = "Species",
continuousParam = "scale(as.numeric(Baiting))", all = TRUE)
# detection
detection_coefs <- predict(fit_top, "det", newdata = data.frame(method = c("cam", "ht"))) %>%
`rownames<-`(c("cam", "ht")) %>% round(3)
We investigated the drivers of site occupancy for four native mammal species: Brushtail Possums (BTP), Long-Nosed Bandicoots (LNB), Long-nosed Potoroos (LNP) and Southern Brown Bandicoots (SBB). We found that fox baiting had variable effect on the site occupancy of the target native species, with the mean slope of the ‘years of baiting’ being slightly negative (\(\beta_{ALL}\) = -0.057 [95% CI = -0.577, 0.471]). The effect of fox baiting had variable magnitudes depending on the target species detected (\(\beta_{BTP}\) = 0.67 [95% CI = 0.126, 1.294], \(\beta_{LNB}\) = 0.442 [95% CI = -0.484, 1.25], \(\beta_{LNP}\) = -0.487 [95% CI = -1.16, 0.113], \(\beta_{SBB}\) = -1.524 [95% CI = -2.499, -0.658]). Figure 5 shows the predicted occupancy rates over different lengths of time since baiting began; these predictions are partitioned across the four species and four sites. These plots show no discernable effect of baiting over time as well as the high variability in occupancy across the sites and locations.
Model covariates used in this study had variable associations with site occupancy depending on the species. Elevation had a positive association with occupancy for Brushtail Possums and a negative association for Long-nosed Potoroos (Figure 7A). Fire frequency had weak positive trends for Brushtail Possums, Long-nosed Potoroos and Southern Brown Bandicoots (Figure 7B). Year had a positive association with occupancy for Long-nosed Bandicoots but a negative association for Southern Brown Bandicoots (Figure 7C). Lastly, rainfall (6-month lag) did not show associations with occupancy for any species (Figure 7D). Model estimates for all covariates are shown in Table 3; estimates with a 95% CI not overlapping zero are shown in bold.
Show code
model_summary_sig <- which((model_summary$`2.5%` * model_summary$`97.5%`) > 0)
model_summary %>%
select(-`25%`, -`50%`, -`75%`) %>%
kbl("html", caption = "Model estimates showing the effects of various covariates on the probability of occupancy. Estimates with a 95% CI not overlapping zero are shown in bold.") %>%
kable_styling(full_width = FALSE) %>%
row_spec(model_summary_sig, bold = TRUE)
| mean | se_mean | sd | 2.5% | 97.5% | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| (Intercept) | 0.413 | 0.007 | 0.233 | -0.044 | 0.876 | 1113.909 | 1.005 |
| scale(Elevation) | 0.307 | 0.004 | 0.174 | -0.034 | 0.647 | 1507.286 | 1.001 |
| SpeciesLNB | -5.235 | 0.016 | 0.598 | -6.485 | -4.158 | 1420.348 | 1.003 |
| SpeciesLNP | -1.937 | 0.007 | 0.284 | -2.511 | -1.398 | 1701.692 | 1.003 |
| SpeciesSBB | -2.225 | 0.007 | 0.295 | -2.833 | -1.675 | 1690.528 | 1.002 |
| scale(as.numeric(Baiting)) | 0.670 | 0.009 | 0.307 | 0.126 | 1.294 | 1255.406 | 1.002 |
| LocationGenoa | 1.325 | 0.015 | 0.492 | 0.411 | 2.357 | 1027.400 | 1.005 |
| LocationMueller | -0.053 | 0.013 | 0.421 | -0.858 | 0.767 | 1071.870 | 1.004 |
| LocationTamboon | -0.644 | 0.012 | 0.433 | -1.501 | 0.201 | 1244.605 | 1.003 |
| scale(FireFreq) | 0.165 | 0.004 | 0.148 | -0.123 | 0.452 | 1425.681 | 1.003 |
scale(6mthlag)
|
0.003 | 0.002 | 0.093 | -0.180 | 0.186 | 2102.119 | 1.000 |
| scale(Elevation):SpeciesLNB | -0.297 | 0.006 | 0.300 | -0.897 | 0.290 | 2696.384 | 1.000 |
| scale(Elevation):SpeciesLNP | -1.334 | 0.005 | 0.219 | -1.770 | -0.911 | 2125.767 | 1.000 |
| scale(Elevation):SpeciesSBB | -0.445 | 0.005 | 0.219 | -0.871 | -0.012 | 2143.398 | 1.000 |
| SpeciesLNB:scale(as.numeric(Baiting)) | -0.229 | 0.012 | 0.499 | -1.243 | 0.691 | 1692.905 | 1.001 |
| SpeciesLNP:scale(as.numeric(Baiting)) | -1.157 | 0.011 | 0.414 | -2.017 | -0.396 | 1489.892 | 1.002 |
| SpeciesSBB:scale(as.numeric(Baiting)) | -2.194 | 0.013 | 0.541 | -3.289 | -1.193 | 1722.536 | 1.000 |
| scale(as.numeric(Baiting)):LocationGenoa | -0.073 | 0.016 | 0.556 | -1.074 | 1.122 | 1174.589 | 1.001 |
| scale(as.numeric(Baiting)):LocationMueller | 1.465 | 0.018 | 0.631 | 0.279 | 2.754 | 1200.436 | 1.002 |
| scale(as.numeric(Baiting)):LocationTamboon | -2.206 | 0.011 | 0.437 | -3.084 | -1.394 | 1527.951 | 1.001 |
| SpeciesLNB:LocationGenoa | -0.042 | 0.018 | 0.773 | -1.531 | 1.522 | 1789.815 | 1.001 |
| SpeciesLNP:LocationGenoa | -1.217 | 0.016 | 0.592 | -2.431 | -0.069 | 1375.589 | 1.002 |
| SpeciesSBB:LocationGenoa | -2.534 | 0.017 | 0.638 | -3.858 | -1.322 | 1460.551 | 1.001 |
| SpeciesLNB:LocationMueller | 3.302 | 0.020 | 0.738 | 1.922 | 4.838 | 1334.972 | 1.003 |
| SpeciesLNP:LocationMueller | 2.308 | 0.013 | 0.503 | 1.372 | 3.322 | 1471.276 | 1.001 |
| SpeciesSBB:LocationMueller | 0.095 | 0.014 | 0.494 | -0.844 | 1.086 | 1337.034 | 1.003 |
| SpeciesLNB:LocationTamboon | -0.098 | 0.021 | 0.918 | -1.929 | 1.709 | 1829.438 | 1.003 |
| SpeciesLNP:LocationTamboon | 0.358 | 0.012 | 0.530 | -0.688 | 1.406 | 1924.188 | 1.003 |
| SpeciesSBB:LocationTamboon | -1.359 | 0.017 | 0.746 | -2.920 | 0.012 | 1825.075 | 1.005 |
| SpeciesLNB:scale(FireFreq) | 0.234 | 0.004 | 0.243 | -0.225 | 0.715 | 2996.276 | 1.000 |
| SpeciesLNP:scale(FireFreq) | 0.486 | 0.003 | 0.178 | 0.138 | 0.835 | 2711.774 | 1.001 |
| SpeciesSBB:scale(FireFreq) | 0.320 | 0.004 | 0.183 | -0.035 | 0.674 | 2613.937 | 1.001 |
SpeciesLNB:scale(6mthlag)
|
0.325 | 0.004 | 0.247 | -0.144 | 0.844 | 3731.526 | 1.000 |
SpeciesLNP:scale(6mthlag)
|
0.175 | 0.003 | 0.142 | -0.097 | 0.448 | 2587.056 | 1.000 |
SpeciesSBB:scale(6mthlag)
|
-0.040 | 0.003 | 0.153 | -0.328 | 0.263 | 2633.687 | 1.000 |
| SpeciesLNB:scale(as.numeric(Baiting)):LocationGenoa | 0.306 | 0.019 | 0.745 | -1.240 | 1.750 | 1487.602 | 1.001 |
| SpeciesLNP:scale(as.numeric(Baiting)):LocationGenoa | 0.155 | 0.019 | 0.700 | -1.308 | 1.453 | 1330.325 | 1.001 |
| SpeciesSBB:scale(as.numeric(Baiting)):LocationGenoa | 1.925 | 0.019 | 0.795 | 0.346 | 3.480 | 1729.234 | 1.000 |
| SpeciesLNB:scale(as.numeric(Baiting)):LocationMueller | -1.136 | 0.020 | 0.785 | -2.687 | 0.430 | 1574.030 | 1.000 |
| SpeciesLNP:scale(as.numeric(Baiting)):LocationMueller | -1.457 | 0.021 | 0.765 | -2.969 | 0.075 | 1300.810 | 1.003 |
| SpeciesSBB:scale(as.numeric(Baiting)):LocationMueller | -0.807 | 0.020 | 0.790 | -2.397 | 0.734 | 1499.980 | 1.000 |
| SpeciesLNB:scale(as.numeric(Baiting)):LocationTamboon | 2.739 | 0.018 | 0.837 | 1.137 | 4.426 | 2177.108 | 1.000 |
| SpeciesLNP:scale(as.numeric(Baiting)):LocationTamboon | 0.635 | 0.015 | 0.631 | -0.595 | 1.847 | 1702.388 | 1.002 |
| SpeciesSBB:scale(as.numeric(Baiting)):LocationTamboon | 1.144 | 0.022 | 0.975 | -0.743 | 3.015 | 2003.509 | 1.003 |
| sigma [1|Site] | 1.055 | 0.003 | 0.097 | 0.877 | 1.256 | 775.845 | 1.006 |
In order to evaluate the benefits of the fox baiting program, the effects of fox baiting on native mammal survival and occupancy need to be made clear and contrasted against a null model (no intervention). Using model predictions we show that the effect that fox baiting has had on small native mammal occupancy cannot be discerned from the available data used in the monitoring of this program (Figure 6).
From the analysis of the observation submodel/process we also found that hair tubes (\(p_{ht}\) = 0.491 [95% CI = 0.474 - 0.509]) have a higher estimated detection probability than camera traps (\(p_{cam}\) = 0.143 [95% CI = 0.122 - 0.167]). Detection probabilities are shown in Figure 4.
References
R Session information
Show code
sessionInfo() %>% pander()
R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8
attached base packages: stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: readxl(v.1.3.1), pander(v.0.6.4), kableExtra(v.1.3.4.9000), ggbeeswarm(v.0.6.0), ggplot2(v.3.3.5), data.table(v.1.14.0), ubms(v.1.0.2.9016), unmarked(v.1.1.1), lattice(v.0.20-44) and dplyr(v.1.0.7)
loaded via a namespace (and not attached): minqa(v.1.2.4), colorspace(v.2.0-2), ellipsis(v.0.3.2), rgdal(v.1.5-23), rprojroot(v.2.0.2), htmlTable(v.2.2.1), base64enc(v.0.1-3), rstudioapi(v.0.13), farver(v.2.1.0), rstan(v.2.21.2), svUnit(v.1.0.6), RSpectra(v.0.16-0), fansi(v.0.5.0), lubridate(v.1.7.10), xml2(v.1.3.2), codetools(v.0.2-18), splines(v.4.1.0), downlit(v.0.2.1), knitr(v.1.33), Formula(v.1.2-4), jsonlite(v.1.7.2), nloptr(v.1.2.2.2), cluster(v.2.1.2), png(v.0.1-7), ggdist(v.3.0.0), compiler(v.4.1.0), httr(v.1.4.2), backports(v.1.2.1), assertthat(v.0.2.1), Matrix(v.1.3-3), cli(v.3.0.1), htmltools(v.0.5.1.1), prettyunits(v.1.1.1), tools(v.4.1.0), coda(v.0.19-4), gtable(v.0.3.0), glue(v.1.4.2), posterior(v.1.0.1), reshape2(v.1.4.4), V8(v.3.4.2), Rcpp(v.1.0.7), cellranger(v.1.1.0), jquerylib(v.0.1.4), raster(v.3.5-2), vctrs(v.0.3.8), svglite(v.2.0.0), nlme(v.3.1-152), tensorA(v.0.36.2), xfun(v.0.24), stringr(v.1.4.0), tinter(v.0.1.0), ps(v.1.6.0), lme4(v.1.1-27.1), rvest(v.1.0.0), lifecycle(v.1.0.1), tidybayes(v.3.0.0), terra(v.1.4-20), MASS(v.7.3-54), scales(v.1.1.1), parallel(v.4.1.0), inline(v.0.3.19), RColorBrewer(v.1.1-2), yaml(v.2.2.1), curl(v.4.3.2), gridExtra(v.2.3), loo(v.2.4.1), StanHeaders(v.2.21.0-7), sass(v.0.4.0), rpart(v.4.1-15), distill(v.1.2), latticeExtra(v.0.6-29), stringi(v.1.7.3), highr(v.0.9), checkmate(v.2.0.0), boot(v.1.3-28), pkgbuild(v.1.2.0), rlang(v.0.4.11), pkgconfig(v.2.0.3), systemfonts(v.1.0.2), matrixStats(v.0.61.0), distributional(v.0.2.2), evaluate(v.0.14), purrr(v.0.3.4), htmlwidgets(v.1.5.3), rstantools(v.2.1.1), chk(v.0.7.0), labeling(v.0.4.2), cowplot(v.1.1.1), processx(v.3.5.2), tidyselect(v.1.1.1), here(v.1.0.1), plyr(v.1.8.6), magrittr(v.2.0.1), R6(v.2.5.1), Hmisc(v.4.5-0), generics(v.0.1.0), DBI(v.1.1.1), foreign(v.0.8-81), pillar(v.1.6.3), withr(v.2.4.2), nnet(v.7.3-16), survival(v.3.2-11), abind(v.1.4-5), sp(v.1.4-5), tibble(v.3.1.5), crayon(v.1.4.1), arrayhelpers(v.1.1-0), utf8(v.1.2.2), rmarkdown(v.2.9), jpeg(v.0.1-8.1), grid(v.4.1.0), callr(v.3.7.0), digest(v.0.6.28), webshot(v.0.5.2), tidyr(v.1.1.3), RcppParallel(v.5.1.4), stats4(v.4.1.0), munsell(v.0.5.0), beeswarm(v.0.4.0), viridisLite(v.0.4.0), vipor(v.0.4.5) and bslib(v.0.2.5.1)