Show code
pre, code {white-space:pre !important; overflow-x:auto}About this document
This following document has been written in RMarkdown and contains embedded analyses using R and STAN programming languages. Using this RMarkdown enables us to keep all the analyses used for the report in a single repository and file. This RMarkdown has been rendered to HTML and is available at https://justincally.github.io/statewide-deer-analysis/. Analyses in this document are carried out sequentially, meaning that some objects, functions or variables needed in later code chunks were made earlier. This repository and RMarkdown also uses a database connection to a secure DEECA database and thus cannot be reproduced without database access. However the main data used in our abundance model has been packaged up as data/multispecies_data.rds and is available in the github repository for users wanting to re-run our model.
An interactive map with estimated average densities for deer in Victoria is also available here: interactive map
Show code
knitr::include_url("interactive-map.html")Methods used to estimate deer density
Camera detection
The density of deer at a given site can be estimated using camera-trap distance sampling (CTDS). CTDS is a modified form of distance-sampling that allows us to infer the probability a given individual will be detected within the survey area (area in front of the camera). This detection probability is a function of the distance of the individual from the camera, whereby individuals entering the camera field of view further from the camera have a lower detection probability up to a given truncation distance from the camera where detection probability is near-zero.
An underlying assumption about CTDS is that the probability a deer will be available for detection at any given point location within the camera field of view is equal. Under this assumption, for a point transect, we take into account the total area for each distance-bin area, which increases at further distance bins. However, in this study we implement a novel method that considers group size of the detected species in the availability calculations. For larger groups, CTDS should account for the availability of the closest individual rather than the availability of all individuals. When group size increases it is more likely the the closest individual from the group is closer to the camera. If we assume that the triggering of the camera is dependent upon the closest individual than we must adjust our estimated availability to account for variable group sizes. If we do not adjust for group size and only use the distance to the closest individual for our DS models then we will likely underestimate the detection rate as distances will be smaller than reality. Alternatively, if we record distances to multiple individuals in the same photo and take an average or model them independently we would likely overestimate detection probability because individuals at further distances are recorded because a closer individual has triggered the camera trap. For more information on how group size is likely to effect CTDS abundance estimates a simulation study has been written here: https://justincally.github.io/blog/posts/2022-11-10-ctds-for-groups/
Supplementing the camera trap data with transect-level data
Our surveys involved the deployment of cameras at a site for usually 6-8 weeks as well as 150m transect searches that were walked up and back by a single individual for two transects, and one transect which was walked once by two observers. This gives us an approximate survey area at each site of 150m x 2 x 3 = 900m. Deer sign was recorded on these transects with pellets, footprints, rubbings, and wallows recorded as either present or absent.
This data provides supplemental presence-absence data that can be integrated into our model to help inform likely densities at sites where cameras did not record observations but one or more deer signs were detected along the transects.
The foundation for this integration is a Royle-Nichols model (Royle and Nichols 2003) that models the frequency of detection as a mixture of detection probability and abundance. Given that we also collect absolute abundance estimates via the camera trap data, relative abundance measures from this output can be transformed to absolute abundance and thus provide an avenue to integrate camera and transect-based data into a single model estimating abundance.
Setup
Load in relevant packages for analysis, additionally, connect to the database. Camera trap and site data is stored on the database.
Show code
library(cmdstanr)
library(dplyr)
library(sf)
library(bayesplot)
library(loo)
library(gt)
library(terra)
library(caret)
library(tidyterra)
library(tidyr)
library(VicmapR)
library(kableExtra)
library(brms)
library(Distance)
library(ggrepel)
library(pROC)
options(brms.backend = "cmdstanr")
check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
register_knitr_engine(override = FALSE)
options(mc.cores=8)
#### Database connection ####
con <- weda::weda_connect(password = keyring::key_get(service = "ari-dev-weda-psql-01",
username = "psql_user"), username = "psql_user")Custom Functions
Additional functions used in the data preparation, modelling and analysis are available in the /functions directory. These include functions relating to plotting, organizing model draws and calculating activity and availability of deer for use within the model.
Show code
# Source internal functions
sapply(list.files("functions", full.names = T), source, verbose = F)Prepare Data
In this section we wrangle and format data for the STAN model. The structure of the data entering the STAN model is specified in the data block of the STAN model.
Scope of models
Outline which species should be modeled, and which projects (state-wide and hog deer) to source data from.
Show code
# Species to run model for.
deer_species_all <- c("Cervus unicolor", "Dama dama", "Cervus elaphus", "Axis porcinus")
species_names <- c("Sambar Deer", "Fallow Deer", "Red Deer", "Hog Deer")
# Projects to select. used in querying database
project_short_name <- c("hog_deer_2023", "StatewideDeer")
# Buffer for data extraction.
spatial_buffer <- 1000
# Covariate Rasters
raster_files <- "/Volumes/DeerVic\ Photos/Processed_Rasters"
prediction_raster <- "data/prediction_raster/statewide_raster.tif"
# For the integrated model we place limits on the maximum density of deer to integrate over. In cases where there are no detections on the camera this is limited to 15 per km2. In cases where deer were detected on the camera this is expanded to 30. We believe this is sufficiently high. Very high values of these will be less efficient.
n_max_no_det <- 15
n_max_det <- 30Camera locations
Download the camera locations from the database, this table outlines the locations and the deployment history of the cameras.
Show code
cams_curated <- tbl(con, dbplyr::in_schema("camtrap", "curated_camtrap_operation")) %>%
dplyr::filter(ProjectShortName %in% !!project_short_name) %>% # Only retrieve cam records for hog deer 2023
dplyr::collect() %>%
dplyr::arrange(SiteID) %>%
sf::st_as_sf(., coords = c("Longitude", "Latitude"), crs = 4283)
n_site <- nrow(cams_curated)Formulas for detection and abundance
Here we outline formulas to be used in the models that we compare. The formulas account for the various fixed-effect parameters, random-effects are hard-coded within the STAN model.
Show code
#### Model formulas ####
#### Transect Formula: Survey Only ####
# Formula used in the RN model
transect_formula <- ~Survey
#### Abundance Formula Options ####
ab_formula_1 <- ~ scale(sqrt(BareSoil)) + scale(sqrt(Nitrogen)) + scale(log(PastureDistance)) + scale(BIO15) + scale(sqrt(ForestEdge)) +
scale(sqrt(MRVBF)) + scale(sqrt(SLOPE)) + scale(sqrt(NonPhotosyntheticVeg))
ab_formula_2 <- ~ scale(sqrt(BareSoil)) + scale(sqrt(Nitrogen)) + scale(log(PastureDistance)) + scale(BIO15) + scale(sqrt(ForestEdge)) + scale(sqrt(MRVBF))
ab_formula_3 <- ~ scale(sqrt(BareSoil)) + scale(sqrt(Nitrogen)) + scale(log(PastureDistance)) + scale(BIO15) + scale(sqrt(ForestEdge))
ab_formula_4 <- ~ 1
#### Detection Formula: Distance-sampling ####
det_formula <- ~ scale(HerbaceousUnderstoryCover) # average detection across all sites Create model data
Using the prepare_model_data() function we generate the data for the STAN model. This function will:
1. Download data from the database
2. Format that data to match the distance-sampling bins
3. Organised the counts into sites and group sizes
4. Generate model matrix of the various sub-models (distance-sampling, abundance and transect) 5. Generate data for the prediction process 6. Generate data for the random effect (bioregion)
7. Generate data for regional predictions (indexing based on DEECA regions)
Show code
if(!file.exists("data/multispecies_data.rds")) {
crown_land <- readRDS("data/crown_land.rds")
multispecies_data <- prepare_model_data_multispecies(species = deer_species_all,
projects = project_short_name,
buffer = spatial_buffer,
detection_formula = det_formula,
abundance_formula = ab_formula_1,
transect_formula = transect_formula,
con = con,
raster_dir = raster_files,
prediction_raster = prediction_raster,
n_max_no_det = n_max_no_det,
n_max_det = n_max_det,
crown_land = crown_land,
evaltransects = TRUE,
filter_behaviour = TRUE)
multispecies_data$keyfun <- 1 # 0 = HN, 1 = HZ
multispecies_data$raw_data[is.na(multispecies_data$raw_data)] <- 0
multispecies_data$transects[is.na(multispecies_data$transects)] <- 0
evc_groups <- readRDS("data/evc_groupnames.rds")
multispecies_data <- c(multispecies_data, evc_groups)
saveRDS(multispecies_data, "data/multispecies_data.rds")
} else {
multispecies_data <- readRDS("data/multispecies_data.rds")
}Note on assigning species to deer sign detections
The recorded deer sign data (pellets, footprints, rubbings and wallows) were not identifiable to deer species. Fresh pellets can be identified to species using genetic sequencing. It is also possible to identify deer species from pellets using visual structure, but this is challenging. Given the variation in surveyors and the uncertainty in visual assignment, we instead assigned species post-hoc via a set of logical steps to minimise false detections of species outside their known range. The known range was determined by creating a kernel density estimate (KDE) of species records over the past 30 years from incidental sightings recorded in the Atlas of Living Australia (ALA). The KDE had a bandwidth determined via the normal reference distribution and was thresholded so that 99.5% of observations fell within the KDE polygon. For a given deer sign detection at a site, we assigned that observation to a given species based on the following stepwise logic:
- species seen on camera = present
- species not seen on camera and another species seen on camera = absent
- species not seen on camera and no others seen, but within KDE distribution = present
- species not seen on camera, and not within KDE distribution = absent
This method may result in false positives or false negatives. At this stage we have not attempted to model those parameters. However, given our use of KDE restrictions we believe that false positives outside the known/existing range will be minimised and will hopefully prevent wider estimated distributions than the actual distribution.
Model Execution
MCMC settings
Below we list the MCMC setting for our model. We run models on eight parallel chains for 400 iterations each (200 warm-up and 200 sampling). These setting provide us with 1,600 posterior draws (8 x 200).
Show code
# STAN settings
ni <- 200 # sampling iterations
nw <- 200 # warm-up iterations
nc <- 8 # number of chainsRead in Models
These are the models used in the analysis. The first is a model using bioregion as a random effect for the abundance process.
Show code
functions {
/* Half-normal function
* Args:
* sigma: sigma
* midpoints: midpoints
* Returns:
* detection probability
*/
vector halfnorm(real sigma, vector midpoints) {
int bins = rows(midpoints);
vector[bins] p_raw; // detection probability
p_raw = exp(-(midpoints)^2/(2*sigma^2));
return p_raw;
}
/* Hazard function
* Args:
* sigma: sigma
* theta: theta
* midpoints: midpoints
* Returns:
* detection probability
*/
vector hazard(real sigma, real theta, vector midpoints) {
int bins = rows(midpoints);
vector[bins] p_raw; // detection probability
p_raw = 1 - exp(-(midpoints/sigma)^(-theta));
return p_raw;
}
vector prob_dist(real sigma, real theta, int keyfun, vector midpoints){
int bins = rows(midpoints);
vector[bins] out; // detection probability
if(keyfun == 0){
out = halfnorm(sigma, midpoints);
} else if(keyfun == 1){
out = hazard(sigma, theta, midpoints);
}
return out;
}
}
data {
int<lower=0> N; // number of observations
int<lower=0> S; // number of species
real delta; // bin width
int<lower=1> n_site; // site
int<lower=1> n_distance_bins; // distance bins
int<lower=1> n_gs; // number of group sizes
vector[n_gs] gs; // group sizes
vector[n_distance_bins] midpts; // midpt of each interval
real<lower=1> max_distance; // truncation distance
int<lower=1> max_int_dist; // max distance as an integer
real<lower=0> theta_frac; // fraction of camera view
array[n_site] int effort; // effort
array[n_site, n_gs, S] int n_obs; //number of observations
array[n_site, n_distance_bins, n_gs] int y; // observations matrix
// summary of whether species is known to be present at each site
int<lower = 0, upper = 1> any_seen[S, n_site];
// number of surveys at each site
int<lower = 0> n_survey[n_site];
// availability prior information
real<lower=0> bshape; // availability shape
real<lower=0> bscale; // availability scale
// camera trap detection parameters
int<lower=0> det_ncb; // number of covariates for detection model
matrix[n_site, det_ncb] det_model_matrix; // detection model matrix
array[n_gs, n_distance_bins] real pa ; // averaged availability for multiple individuals
// Abundance/occupancy model matrix
int<lower = 1> m_psi;
matrix[n_site, m_psi] X_psi;
// negbinom scale
// real reciprocal_phi_scale;
//transect level information
int<lower=1> trans; // total number of transects across all sites for all methods
array[S, trans] int<lower = 0, upper = 1> y2; // transect based binomial detections
int<lower = 0, upper = trans> start_idx[n_site];
int<lower = 0, upper = trans> end_idx[n_site];
int<lower=1> trans_det_ncb; // number of covariates for transect detection model
matrix[trans, trans_det_ncb] trans_pred_matrix; // transect detection model matrix
int<lower=1> n_max[n_site, S]; // max for poisson RN
// Prediction data
int<lower=1> npc; // number of prediction grid cells
matrix[npc, m_psi] X_pred_psi; // pred matrix
vector[npc] prop_pred; //offset
// bioregion RE
int<lower=1> np_bioreg;
int<lower=1> site_bioreg[n_site];
int<lower=1> pred_bioreg[npc];
// region data
int<lower=1> np_reg;
int<lower=1> site_reg[n_site];
int<lower=1> pred_reg[npc];
// key function, 0 = halfnorm
int keyfun;
}
transformed data {
// vector[n_distance_bins] pi; // availability for point
vector[n_site] survey_area;
array[S, n_site] real cam_seen;
for(i in 1:n_site) {
survey_area[i] = theta_frac * effort[i] * pi() * (max_distance/1000)^2;
for(s in 1:S) {
cam_seen[s,i] = sum(n_obs[i,,s]);
}
}
}
parameters {
// abundance parameters
//array[n_site] simplex[n_gs] eps_ngs; // random variation in group size
array[S] vector[m_psi] beta_psi;
vector[det_ncb] beta_det;
real log_theta;
// transect detection parameters
vector[trans_det_ncb] beta_trans_det;
// temporal availability parameters
real<lower=0, upper=1> activ;
// bioregion RE
array[S] real<lower=0> bioregion_sd;
array[S] vector[np_bioreg] bioregion_raw;
// eps group size params
array[S] vector[n_gs] zeta;
array[S] matrix[n_site, n_gs] eps_raw;
array[S] real<lower=0> grp_sd;
// od
array[S] real od_mu;
real odRNmu;
}
transformed parameters {
// random effects
array[S] vector[np_bioreg] eps_bioregion; // bioregion random effect
// distance parameters
array[n_site] real log_sigma;
array[n_site] real sigma;
array[n_site] vector[n_distance_bins] p_raw; // detection probability
array[n_site, n_distance_bins, n_gs] real p_raw_scale; // detection probability for point independence model and accounting for availability
array[n_site, n_distance_bins, n_gs] real<upper=0> log_p_raw; // log detection probability accounting for availability
array[n_site, n_gs] real log_p; // log probability with summed exponential for the multinomial logit model
array[n_site, n_gs] real<lower=0,upper=1> p; // log probability with summed exponential for the multinomial logit model
// abundance parameters
array[S, n_site, n_gs] real<lower=0> lambda;
// activity parameters
real log_activ = log(activ);
vector[trans] logit_trans_p = trans_pred_matrix * beta_trans_det; // observation process model
// lp_site for RN model
array[S, n_site] real lp_site;
vector<lower=0,upper=1>[trans] r = inv_logit(logit_trans_p);
array[S] vector[n_site] log_lambda_psi;
// eps group size
array[S, n_site] simplex[n_gs] eps_ngs; // random variation in group size
array[S, n_site] vector[n_gs] epsi;
array[S] matrix[n_site, n_gs] eps_site;
real<lower=0> theta = exp(log_theta);
array[S] real od; // bioregion random effect
real odRN = exp(odRNmu); // bioregion random effect
array[S, n_site] real<lower=0, upper=1> pbar;
for(s in 1:S) {
for(b in 1:np_bioreg) {
eps_bioregion[s,b] = bioregion_sd[s] * bioregion_raw[s,b];
}
od[s] = exp(odRNmu + od_mu[s]);
}
for(n in 1:n_site) {
log_sigma[n] = det_model_matrix[n,] * beta_det;
sigma[n] = exp(log_sigma[n]);
p_raw[n] = prob_dist(sigma[n], theta, keyfun, midpts);
for (i in 1:n_distance_bins) {
// assuming a half-normal detection fn from line
// p_raw[n,i] = exp(-(midpts[i])^2/(2*sigma[n]^2)); // half-normal function (pg 419 of AHM)
// hazard function
// p_raw[n,i] = 1 - exp(-(midpts[i]/theta)^(-sigma[n])); // hazard function (pg 507 of AHM)
for(j in 1:n_gs) {
p_raw_scale[n,i,j] = p_raw[n,i]*pa[j,i]; // pr(animal occurs and is detected in bin i)
log_p_raw[n,i,j] = log(p_raw_scale[n,i,j]);
}
}
for(s in 1:S) {
vector[n_max[n,s]+1] Nlik;
vector[n_max[n,s]+1] gN;
vector[n_max[n,s]+1] pcond;
log_lambda_psi[s,n] = X_psi[n,] * beta_psi[s] + eps_bioregion[s,site_bioreg[n]];
for(j in 1:n_gs) {
eps_site[s, n,j] = grp_sd[s] * eps_raw[s,n,j];
epsi[s,n,j] = exp(zeta[s,j] + eps_site[s,n,j]);
}
eps_ngs[s,n,] = epsi[s,n,]/sum(epsi[s,n,]);
// p-bar
for(k in 1:(n_max[n,s]+1))
Nlik[k] = exp(neg_binomial_2_log_lpmf(k-1 | log_lambda_psi[s,n], odRN));
gN = Nlik/sum(Nlik);
for(k in 1:(n_max[n,s]+1))
pcond[k] = (1 - (1 - r[start_idx[n]])^(k-1)) * gN[k];
pbar[s,n] = sum(pcond);
}
for(j in 1:n_gs) {
log_p[n,j] = log_sum_exp(log_p_raw[n,,j]);
p[n,j] = exp(log_p[n,j]);
// model site abundance
for(s in 1:S) {
lambda[s,n,j] = exp(log_lambda_psi[s,n] + log_p[n,j] + log_activ + log(eps_ngs[s,n,j])) .* survey_area[n];
}
}
// Royle-Nichols implementation in STAN (looping over possible discrete values of N)
// https://discourse.mc-stan.org/t/royle-and-nichols/14150
// https://discourse.mc-stan.org/t/identifiability-across-levels-in-occupancy-model/5340/2
for(s in 1:S) {
if (n_survey[n] > 0) {
vector[n_max[n,s]] lp;
if(any_seen[s,n] == 0){ // not seen
lp[1] = log_sum_exp(neg_binomial_2_log_lpmf(0 | log_lambda_psi[s,n], odRN), neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]));
} else {
lp[1] = neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]);
}
for (k in 2:n_max[n,s]){
lp[k] = neg_binomial_2_log_lpmf(k | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | 1-(1-r[start_idx[n]:end_idx[n]])^k);
}
lp_site[s,n] = log_sum_exp(lp);
} else {
lp_site[s, n] = 0;
}
}
}
}
model {
for(s in 1:S) {
beta_psi[s] ~ normal(0, 3); // prior for intercept in poisson model
bioregion_sd[s] ~ normal(0,2);
bioregion_raw[s,] ~ normal(0,1);
to_vector(eps_raw[s,,]) ~ std_normal();
grp_sd[s] ~ normal(0, 1);
zeta[s,] ~ normal(0, 2);
// od
od_mu[s] ~ normal(0,1);
}
odRNmu ~ normal(0,1);
beta_trans_det ~ normal(0, 2); // beta for transect detection
beta_det ~ normal(0, 4); // prior for sigma
activ ~ beta(bshape, bscale); //informative prior
log_theta ~ normal(0,2); // prior for theta
for(n in 1:n_site) {
for(j in 1:n_gs) {
for(s in 1:S) {
target += neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s]);
}
y[n,,j] ~ multinomial_logit(to_vector(log_p_raw[n,,j]));
}
for(s in 1:S) {
target += lp_site[s, n];
}
}
}
generated quantities {
array[S, n_site,n_gs] real n_obs_pred;
array[S, n_site, n_gs] real n_obs_true;
array[S, n_site] real N_site;
array[S, n_site] real N_site_pred;
array[n_site, max_int_dist+1] real DetCurve;
array[n_site, n_gs] real log_lik1;
array[S, n_site, n_gs] real log_lik2;
array[n_site, n_gs] real log_lik2_site;
array[S, n_site] real log_lik2_species;
array[n_site] real log_lik;
array[n_site] real log_lik_det;
//array[S, n_site] real Site_lambda;
real av_gs[S];
array[S] simplex[n_gs] eps_gs_ave;
array[S, npc] real pred;
//array[S, npc] real pred_trunc;
array[S, np_reg] real Nhat_reg;
array[S, np_bioreg] real Nhat_bioreg;
// array[np_reg] real Nhat_reg_design;
real Nhat[S];
//real Nhat_trunc[S];
real Nhat_sum;
int trunc_counter[S];
trunc_counter[S] = 0;
for(s in 1:S) {
eps_gs_ave[s] = exp(zeta[s])/sum(exp(zeta[s]));
av_gs[s] = sum(gs .* eps_gs_ave[s]); //average group size
}
for(n in 1:n_site) {
for(j in 1:n_gs) {
log_lik1[n,j] = multinomial_logit_lpmf(y[n,,j] | to_vector(log_p_raw[n,,j])); //for loo
for(s in 1:S) {
log_lik2[s,n,j] = neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s]); //for loo
n_obs_true[s, n, j] = gs[j] * (neg_binomial_2_log_rng(log_lambda_psi[s,n] + log(eps_ngs[s,n,j]), od[s]));
n_obs_pred[s, n,j] = gs[j] * neg_binomial_2_rng(lambda[s,n,j], od[s]);
}
log_lik2_site[n, j] = log_sum_exp(log_lik2[,n,j]);
}
// get loglik on a site level
log_lik_det[n] = log_sum_exp(log_lik1[n,]);
log_lik[n] = log_sum_exp(log_sum_exp(log_lik_det[n],
log_sum_exp(log_lik2_site[n,])), log_sum_exp(lp_site[,n]));
for(s in 1:S) {
//Site_lambda[s,n] = exp(log_lambda_psi[s,n]);
log_lik2_species[s, n] = log_sum_exp(log_sum_exp(log_lik2[s,n,]), lp_site[s,n]);
N_site[s,n] = sum(n_obs_true[s,n,]);
N_site_pred[s,n] = sum(n_obs_pred[s,n,]);
}
// loop over distance bins
if(keyfun == 0) {
for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
DetCurve[n, j+1] = exp(-(j+0.5)^2/(2*sigma[n]^2)); // half normal
}
} else if(keyfun == 1) {
for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
DetCurve[n, j+1] = 1 - exp(-((j+0.5)/sigma[n])^(-theta)); //hazard rate
}
}
}
for(s in 1:S) {
for(i in 1:np_reg) Nhat_reg[s,i] = 0;
for(i in 1:np_bioreg) Nhat_bioreg[s,i] = 0;
for(i in 1:npc) {
pred[s,i] = neg_binomial_2_log_rng(X_pred_psi[i,] * beta_psi[s] + eps_bioregion[s, pred_bioreg[i]], od[s]) * prop_pred[i] * av_gs[s]; //offset
if(pred[s,i] > max(N_site[s,])) {
trunc_counter[s] += 1;
}
// for Hog Deer limit to Gippsland region
if(pred_reg[i] != 2 && s == 4) {
pred[s,i] = 0;
}
Nhat_reg[s,pred_reg[i]] += pred[s,i];
Nhat_bioreg[s,pred_bioreg[i]] += pred[s,i];
}
Nhat[s] = sum(pred[s,]);
//Nhat_trunc[s] = sum(pred_trunc[s,]);
}
Nhat_sum = sum(Nhat);
}The second model uses the combination of bioregion and EVC as random-effects on the abundance process.
Show code
functions {
/* Half-normal function
* Args:
* sigma: sigma
* midpoints: midpoints
* Returns:
* detection probability
*/
vector halfnorm(real sigma, vector midpoints) {
int bins = rows(midpoints);
vector[bins] p_raw; // detection probability
p_raw = exp(-(midpoints)^2/(2*sigma^2));
return p_raw;
}
/* Hazard function
* Args:
* sigma: sigma
* theta: theta
* midpoints: midpoints
* Returns:
* detection probability
*/
vector hazard(real sigma, real theta, vector midpoints) {
int bins = rows(midpoints);
vector[bins] p_raw; // detection probability
p_raw = 1 - exp(-(midpoints/sigma)^(-theta));
return p_raw;
}
vector prob_dist(real sigma, real theta, int keyfun, vector midpoints){
int bins = rows(midpoints);
vector[bins] out; // detection probability
if(keyfun == 0){
out = halfnorm(sigma, midpoints);
} else if(keyfun == 1){
out = hazard(sigma, theta, midpoints);
}
return out;
}
}
data {
int<lower=0> N; // number of observations
int<lower=0> S; // number of species
real delta; // bin width
int<lower=1> n_site; // site
int<lower=1> n_distance_bins; // distance bins
int<lower=1> n_gs; // number of group sizes
vector[n_gs] gs; // group sizes
vector[n_distance_bins] midpts; // midpt of each interval
real<lower=1> max_distance; // truncation distance
int<lower=1> max_int_dist; // max distance as an integer
real<lower=0> theta_frac; // fraction of camera view
array[n_site] int effort; // effort
array[n_site, n_gs, S] int n_obs; //number of observations
array[n_site, n_distance_bins, n_gs] int y; // observations matrix
// summary of whether species is known to be present at each site
int<lower = 0, upper = 1> any_seen[S, n_site];
// number of surveys at each site
int<lower = 0> n_survey[n_site];
// availability prior information
real<lower=0> bshape; // availability shape
real<lower=0> bscale; // availability scale
// camera trap detection parameters
int<lower=0> det_ncb; // number of covariates for detection model
matrix[n_site, det_ncb] det_model_matrix; // detection model matrix
array[n_gs, n_distance_bins] real pa ; // averaged availability for multiple individuals
// Abundance/occupancy model matrix
int<lower = 1> m_psi;
matrix[n_site, m_psi] X_psi;
// negbinom scale
// real reciprocal_phi_scale;
//transect level information
int<lower=1> trans; // total number of transects across all sites for all methods
array[S, trans] int<lower = 0, upper = 1> y2; // transect based binomial detections
int<lower = 0, upper = trans> start_idx[n_site];
int<lower = 0, upper = trans> end_idx[n_site];
int<lower=1> trans_det_ncb; // number of covariates for transect detection model
matrix[trans, trans_det_ncb] trans_pred_matrix; // transect detection model matrix
int<lower=1> n_max[n_site, S]; // max for poisson RN
// Prediction data
int<lower=1> npc; // number of prediction grid cells
matrix[npc, m_psi] X_pred_psi; // pred matrix
vector[npc] prop_pred; //offset
// bioregion RE
int<lower=1> np_bioreg;
int<lower=1> site_bioreg[n_site];
int<lower=1> pred_bioreg[npc];
// region data
int<lower=1> np_reg;
int<lower=1> site_reg[n_site];
int<lower=1> pred_reg[npc];
// evc data
int<lower=1> np_evc;
int<lower=1> site_evc[n_site];
int<lower=1> pred_evc[npc];
// key function, 0 = halfnorm
int keyfun;
}
transformed data {
// vector[n_distance_bins] pi; // availability for point
vector[n_site] survey_area;
array[S, n_site] real cam_seen;
for(i in 1:n_site) {
survey_area[i] = theta_frac * effort[i] * pi() * (max_distance/1000)^2;
for(s in 1:S) {
cam_seen[s,i] = sum(n_obs[i,,s]);
}
}
}
parameters {
// abundance parameters
//array[n_site] simplex[n_gs] eps_ngs; // random variation in group size
array[S] vector[m_psi] beta_psi;
vector[det_ncb] beta_det;
real log_theta;
// transect detection parameters
vector[trans_det_ncb] beta_trans_det;
// temporal availability parameters
real<lower=0, upper=1> activ;
// bioregion RE
array[S] real<lower=0> bioregion_sd;
array[S] vector[np_bioreg] bioregion_raw;
// evc RE
real<lower=0> evc_sd;
array[S] vector[np_evc] evc_raw;
// eps group size params
array[S] vector[n_gs] zeta;
array[S] matrix[n_site, n_gs] eps_raw;
array[S] real<lower=0> grp_sd;
// od
array[S] real od_mu;
real odRNmu;
}
transformed parameters {
// random effects
array[S] vector[np_bioreg] eps_bioregion; // bioregion random effect
array[S] vector[np_evc] eps_evc; // bioregion random effect
// distance parameters
array[n_site] real log_sigma;
array[n_site] real sigma;
array[n_site] vector[n_distance_bins] p_raw; // detection probability
array[n_site, n_distance_bins, n_gs] real p_raw_scale; // detection probability for point independence model and accounting for availability
array[n_site, n_distance_bins, n_gs] real<upper=0> log_p_raw; // log detection probability accounting for availability
array[n_site, n_gs] real log_p; // log probability with summed exponential for the multinomial logit model
array[n_site, n_gs] real<lower=0,upper=1> p; // log probability with summed exponential for the multinomial logit model
// abundance parameters
array[S, n_site, n_gs] real<lower=0> lambda;
// activity parameters
real log_activ = log(activ);
vector[trans] logit_trans_p = trans_pred_matrix * beta_trans_det; // observation process model
// lp_site for RN model
array[S, n_site] real lp_site;
vector<lower=0,upper=1>[trans] r = inv_logit(logit_trans_p);
array[S] vector[n_site] log_lambda_psi;
// eps group size
array[S, n_site] simplex[n_gs] eps_ngs; // random variation in group size
array[S, n_site] vector[n_gs] epsi;
array[S] matrix[n_site, n_gs] eps_site;
real<lower=0> theta = exp(log_theta);
array[S] real od; // bioregion random effect
real odRN = exp(odRNmu); // bioregion random effect
array[S, n_site] real<lower=0, upper=1> pbar;
for(s in 1:S) {
for(b in 1:np_bioreg) {
eps_bioregion[s,b] = bioregion_sd[s] * bioregion_raw[s,b];
}
od[s] = exp(odRNmu + od_mu[s]);
for(k in 1:np_evc) {
eps_evc[s,k] = evc_sd * evc_raw[s,k];
}
}
for(n in 1:n_site) {
log_sigma[n] = det_model_matrix[n,] * beta_det;
sigma[n] = exp(log_sigma[n]);
p_raw[n] = prob_dist(sigma[n], theta, keyfun, midpts);
for (i in 1:n_distance_bins) {
// assuming a half-normal detection fn from line
// p_raw[n,i] = exp(-(midpts[i])^2/(2*sigma[n]^2)); // half-normal function (pg 419 of AHM)
// hazard function
// p_raw[n,i] = 1 - exp(-(midpts[i]/theta)^(-sigma[n])); // hazard function (pg 507 of AHM)
for(j in 1:n_gs) {
p_raw_scale[n,i,j] = p_raw[n,i]*pa[j,i]; // pr(animal occurs and is detected in bin i)
log_p_raw[n,i,j] = log(p_raw_scale[n,i,j]);
}
}
for(s in 1:S) {
vector[n_max[n,s]+1] Nlik;
vector[n_max[n,s]+1] gN;
vector[n_max[n,s]+1] pcond;
log_lambda_psi[s,n] = X_psi[n,] * beta_psi[s] + eps_bioregion[s,site_bioreg[n]] + eps_evc[s,site_evc[n]];
for(j in 1:n_gs) {
eps_site[s, n,j] = grp_sd[s] * eps_raw[s,n,j];
epsi[s,n,j] = exp(zeta[s,j] + eps_site[s,n,j]);
}
eps_ngs[s,n,] = epsi[s,n,]/sum(epsi[s,n,]);
// p-bar
for(k in 1:(n_max[n,s]+1))
Nlik[k] = exp(neg_binomial_2_log_lpmf(k-1 | log_lambda_psi[s,n], odRN));
gN = Nlik/sum(Nlik);
for(k in 1:(n_max[n,s]+1))
pcond[k] = (1 - (1 - r[start_idx[n]])^(k-1)) * gN[k];
pbar[s,n] = sum(pcond);
}
for(j in 1:n_gs) {
log_p[n,j] = log_sum_exp(log_p_raw[n,,j]);
p[n,j] = exp(log_p[n,j]);
// model site abundance
for(s in 1:S) {
lambda[s,n,j] = exp(log_lambda_psi[s,n] + log_p[n,j] + log_activ + log(eps_ngs[s,n,j])) .* survey_area[n];
}
}
// Royle-Nichols implementation in STAN (looping over possible discrete values of N)
// https://discourse.mc-stan.org/t/royle-and-nichols/14150
// https://discourse.mc-stan.org/t/identifiability-across-levels-in-occupancy-model/5340/2
for(s in 1:S) {
if (n_survey[n] > 0) {
vector[n_max[n,s]] lp;
if(any_seen[s,n] == 0){ // not seen
lp[1] = log_sum_exp(neg_binomial_2_log_lpmf(0 | log_lambda_psi[s,n], odRN), neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]));
} else {
lp[1] = neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]);
}
for (k in 2:n_max[n,s]){
lp[k] = neg_binomial_2_log_lpmf(k | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | 1-(1-r[start_idx[n]:end_idx[n]])^k);
}
lp_site[s,n] = log_sum_exp(lp);
} else {
lp_site[s, n] = 0;
}
}
}
}
model {
for(s in 1:S) {
beta_psi[s] ~ normal(0, 3); // prior for intercept in poisson model
bioregion_sd[s] ~ normal(0,2);
bioregion_raw[s,] ~ normal(0,1);
to_vector(eps_raw[s,,]) ~ std_normal();
grp_sd[s] ~ normal(0, 1);
zeta[s,] ~ normal(0, 2);
// od
od_mu[s] ~ normal(0,1);
// evc re
evc_raw[s,] ~ normal(0,1);
}
evc_sd ~ normal(0,1);
odRNmu ~ normal(0,1);
beta_trans_det ~ normal(0, 2); // beta for transect detection
beta_det ~ normal(0, 4); // prior for sigma
activ ~ beta(bshape, bscale); //informative prior
log_theta ~ normal(0,2); // prior for theta
for(n in 1:n_site) {
for(j in 1:n_gs) {
for(s in 1:S) {
target += neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s]);
}
y[n,,j] ~ multinomial_logit(to_vector(log_p_raw[n,,j]));
}
for(s in 1:S) {
target += lp_site[s, n];
}
}
}
generated quantities {
array[S, n_site,n_gs] real n_obs_pred;
array[S, n_site, n_gs] real n_obs_true;
array[S, n_site] real N_site;
array[S, n_site] real N_site_pred;
array[n_site, max_int_dist+1] real DetCurve;
array[n_site, n_gs] real log_lik1;
array[S, n_site, n_gs] real log_lik2;
array[n_site, n_gs] real log_lik2_site;
array[S, n_site] real log_lik2_species;
array[n_site] real log_lik;
array[n_site] real log_lik_det;
//array[S, n_site] real Site_lambda;
real av_gs[S];
array[S] simplex[n_gs] eps_gs_ave;
array[S, npc] real pred;
//array[S, npc] real pred_trunc;
array[S, np_reg] real Nhat_reg;
// array[np_reg] real Nhat_reg_design;
real Nhat[S];
//real Nhat_trunc[S];
real Nhat_sum;
int trunc_counter[S];
trunc_counter[S] = 0;
for(s in 1:S) {
eps_gs_ave[s] = exp(zeta[s])/sum(exp(zeta[s]));
av_gs[s] = sum(gs .* eps_gs_ave[s]); //average group size
}
for(n in 1:n_site) {
for(j in 1:n_gs) {
log_lik1[n,j] = multinomial_logit_lpmf(y[n,,j] | to_vector(log_p_raw[n,,j])); //for loo
for(s in 1:S) {
log_lik2[s,n,j] = neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s]); //for loo
n_obs_true[s, n, j] = gs[j] * (neg_binomial_2_log_rng(log_lambda_psi[s,n] + log(eps_ngs[s,n,j]), od[s]));
n_obs_pred[s, n,j] = gs[j] * neg_binomial_2_rng(lambda[s,n,j], od[s]);
}
log_lik2_site[n, j] = log_sum_exp(log_lik2[,n,j]);
}
// get loglik on a site level
log_lik_det[n] = log_sum_exp(log_lik1[n,]);
log_lik[n] = log_sum_exp(log_sum_exp(log_lik_det[n],
log_sum_exp(log_lik2_site[n,])), log_sum_exp(lp_site[,n]));
for(s in 1:S) {
//Site_lambda[s,n] = exp(log_lambda_psi[s,n]);
log_lik2_species[s, n] = log_sum_exp(log_sum_exp(log_lik2[s,n,]), lp_site[s,n]);
N_site[s,n] = sum(n_obs_true[s,n,]);
N_site_pred[s,n] = sum(n_obs_pred[s,n,]);
}
// loop over distance bins
if(keyfun == 0) {
for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
DetCurve[n, j+1] = exp(-(j+0.5)^2/(2*sigma[n]^2)); // half normal
}
} else if(keyfun == 1) {
for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
DetCurve[n, j+1] = 1 - exp(-((j+0.5)/sigma[n])^(-theta)); //hazard rate
}
}
}
for(s in 1:S) {
for(i in 1:np_reg) Nhat_reg[s,i] = 0;
for(i in 1:npc) {
pred[s,i] = neg_binomial_2_log_rng(X_pred_psi[i,] * beta_psi[s] + eps_bioregion[s, pred_bioreg[i]] + eps_evc[s, pred_evc[i]], od[s]) * prop_pred[i] * av_gs[s]; //offset
if(pred[s,i] > max(N_site[s,])) {
//pred_trunc[s,i] = max(N_site[s,]);
trunc_counter[s] += 1;
}
// for Hog Deer limit to Gippsland region
if(pred_reg[i] != 2 && s == 4) {
pred[s,i] = 0;
}
Nhat_reg[s,pred_reg[i]] += pred[s,i];
}
Nhat[s] = sum(pred[s,]);
//Nhat_trunc[s] = sum(pred_trunc[s,]);
}
Nhat_sum = sum(Nhat);
}Fit models
Detection Models
Using the Distance R package (Miller et al. 2019) we compare several detection functions (hazard and half-normal) with the inclusion of a parameter (herbaceous understorey cover), which may impact detection rates. The top model (as chosen by AIC) from this selection will be included in our full STAN model. Note that we use this approach only on photos where group size = 1. However, the model then used with the STAN model uses data from all group sizes. This following section is just to help us select the best detection submodel to include within our integrated STAN model. Previous attempts to use Leave-One-Out Cross Validation (LOO-CV) (Vehtari et al. 2020) for a suite of Bayesian detection models were not informative, thus we use this more conventional approach.
Show code
filter_behaviour <- TRUE
snapshot_interval <- 2
# Camera information
theta <- 40 * pi / 180 # camera angle in radians
### Get camera information for the camera sites
### For cameras that have problems you need to use the date from when the problem started
cams_curated2 <- dplyr::tbl(con, dbplyr::in_schema("camtrap", "curated_camtrap_operation")) %>%
dplyr::filter(ProjectShortName %in% !!project_short_name) %>% # Only retrieve cam records for hog deer 2023
dplyr::collect() %>%
dplyr::mutate(DateTimeDeploy = as.POSIXct(DateTimeDeploy),
DateTimeRetrieve = as.POSIXct(DateTimeRetrieve),
Tk = as.numeric(DateTimeRetrieve - DateTimeDeploy, units = "secs"), #seconds
Tk_prob = dplyr::coalesce(as.numeric(as.POSIXct(Problem1_to,
format = "%Y-%m-%d %H:%M:%OS") -
as.POSIXct(Problem1_from, format = "%Y-%m-%d %H:%M:%OS"),
units = "secs"), 0),
Tk_adj = Tk-Tk_prob,
Tkt = Tk_adj / snapshot_interval, # snapshot moments: every second second
Effort = (Tk_adj*theta)/(snapshot_interval * pi * 2)) %>%
dplyr::arrange(SiteID)
### Get the deer records (hog deer)
camtrap_records_deer <- dplyr::tbl(con, dbplyr::in_schema(schema = "camtrap", table = "curated_camtrap_records")) %>%
dplyr::filter(scientific_name %in% deer_species_all & ProjectShortName %in% !!project_short_name) %>%
dplyr::select(Species = scientific_name, SiteID, Distance = metadata_Distance, size = metadata_Multiples, Date, Time, Behaviour = metadata_Behaviour) %>%
dplyr::collect() %>%
dplyr::mutate(Distance = dplyr::case_when(Distance == "NA" | is.na(Distance) ~ "999",
TRUE ~ Distance)) %>%
dplyr::mutate(Time = as.POSIXct(Time, format = "%H:%M:%OS"),
Time_n = as.numeric(Time, units = "secs"),
Behaviour = na_if(x = Behaviour, y = "NA")) %>%
dplyr::rowwise() %>%
dplyr::mutate(DistanceMod = list(stringr::str_split(Distance, pattern = "_&_")[[1]])) %>%
dplyr::mutate(Distance = DistanceMod[which.min(as.numeric(stringr::str_extract(DistanceMod, pattern = "[0-9]+")))]) %>%
dplyr::mutate(Distance = dplyr::case_when(Distance == "0-2.5" ~ "0 - 2.5",
Distance == "999" ~ NA_character_,
TRUE ~ Distance)) %>%
dplyr::ungroup() %>%
dplyr::filter(Time_n %% snapshot_interval == 0) %>% #& #snapshot moment interval of 2s
{if(filter_behaviour) dplyr::filter(., is.na(.data[["Behaviour"]])) else .} %>% # filter out behaviors such as camera or marker interaction
dplyr::group_by(SiteID, Time_n, Species) %>%
dplyr::slice(1) %>% # if two photos occur in one second take only one (snapshot moment = 2)
dplyr::ungroup()
### Tidy format for the records
dcount<- camtrap_records_deer %>%
dplyr::select(Species, SiteID, Distance, size, Date, Time, Behaviour) %>%
dplyr::full_join(cams_curated2 %>%
dplyr::select(SiteID, Tkt, Effort), by="SiteID") %>%
dplyr::mutate(object=1:nrow(.)) %>%
dplyr::mutate(size = if_else(is.na(size),0L, size)) %>%
dplyr::arrange(SiteID)
### Format data in a time format for availability
summarised_count <- dcount
summarised_count$Distance <- factor(summarised_count$Distance, levels = c("0 - 2.5", "2.5 - 5", "5 - 7.5", "7.5 - 10", "10+"))
summarised_count <- summarised_count %>%
mutate(distance = case_when(Distance == "0 - 2.5" ~ 1.25,
Distance == "2.5 - 5" ~ 3.75,
Distance == "5 - 7.5" ~ 6.25,
Distance == "7.5 - 10" ~ 8.75,
Distance == "10+" ~ 11.25)) %>%
filter(size == 1)
# site variables
site_vars <- cams_curated2 %>%
left_join(dplyr::tbl(con, dbplyr::in_schema("deervic", "curated_site_data")) %>%
dplyr::filter(SiteID %in% !!c(cams_curated2$SiteID, "47191")) %>%
dplyr::collect() %>%
dplyr::mutate(HerbaceousUnderstoryCover = scale(NNWHUCover + ENWHUCover),
WoodyUnderstoryCover = scale(NWUCover + EWUCover),
Trail = factor(Trail),
SiteID = dplyr::case_when(SiteID == "47191" & CameraID == "HO04101053" ~ "47191A",
SiteID == "47191" & CameraID != "HO04101053" ~ "47191B",
TRUE ~ SiteID))) %>% # native + exotic herbaceous cover
dplyr::arrange(SiteID) %>%
as.data.frame() %>%
select(-Time, -Tkt, -Date, -Effort)
ds_data <- summarised_count %>% left_join(site_vars) %>%
filter(!is.na(HerbaceousUnderstoryCover))
mybreaks <- c(0,2.5,5,7.5,10,12.5)
trunc.list <- list(left=0, right=12.5)
conversion <- convert_units("metre", "metre", "metre")
# Different models to test
hn0 <- ds(ds_data, transect = "point", key="hn", adjustment = NULL,
convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
hr0 <- ds(ds_data, transect = "point", key="hr", adjustment = NULL,
convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
hn1 <- ds(ds_data, transect = "point", key="hn", adjustment = NULL,
formula = ~HerbaceousUnderstoryCover,
convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
hr1 <- ds(ds_data, transect = "point", key="hr", adjustment = NULL,
formula = ~HerbaceousUnderstoryCover,
convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
kableExtra::kbl(summarize_ds_models(hn0, hr0, hn1, hr1, output = "latex") %>%
mutate(Model = stringr::str_extract(Model, pattern = '(?<=\\{)[^\\}]+')),
digits=3,
caption="Model selection table for deer detection", format = "html") %>%
kable_styling(bootstrap_options = "striped")| Model | Key function | Formula | \(\chi^2\) \(p\)-value | \(\hat{P_a}\) | se(\(\hat{P_a}\)) | \(\Delta\)AIC | |
|---|---|---|---|---|---|---|---|
| 4 | hr1 | Hazard-rate | ~HerbaceousUnderstoryCover | 0 | 0.230 | 0.008 | 0.000 |
| 2 | hr0 | Hazard-rate | ~1 | 0 | 0.228 | 0.008 | 3.702 |
| 1 | hn0 | Half-normal | ~1 | 0 | 0.336 | 0.004 | 185.904 |
| 3 | hn1 | Half-normal | ~HerbaceousUnderstoryCover | 0 | 0.336 | 0.004 | 186.859 |
Model Options
Below we list the models we are implementing, all have a hazard detection function, but we compare four different abundance formulas as well as the inclusion of the EVC random effect within the abundance process.
Show code
#Models to run
models_to_run <- data.frame(formula = c("ab_formula_1",
"ab_formula_2",
"ab_formula_3",
"ab_formula_4",
"ab_formula_3",
"ab_formula_4",
"ab_formula_1",
"ab_formula_2"),
keyfun = c(1,1,1,1,1,1,1,1),
evc = c(F,F,F,F,T,T,T,T))Full Models
We can fit models using the cmdstanr R package (Gabry and Češnovar 2022). Here we fit models using a negative-binomial distribution. The models we compare are:
- Model with several key abundance covariates and hazard-rate detection function and bioregion random effect
- Model with several key abundance covariates and hazard-rate detection function and bioregion + EVC random effect
To avoid re-running models this chunk is not evaluated by default and the models are saved as rds objects
Show code
# Run the models and assign them to an element in a list matching the rows for the 'models_to_run' data.frame
model_fits <- list()
for(i in 1:nrow(models_to_run)) {
# for each iteration run a different model
if(models_to_run$evc[i]) {
model_to_fit <- model_negbin_ms_evc
} else {
model_to_fit <- model_negbin_ms
}
form_to_use <- get(models_to_run$formula[i])
params_to_use <- c(TRUE, labels(terms(ab_formula_1)) %in% labels(terms(form_to_use)))
data_to_use <- multispecies_data
data_to_use$X_psi <- as.matrix(data_to_use$X_psi[,params_to_use])
data_to_use$X_pred_psi <- as.matrix(data_to_use$X_pred_psi[,params_to_use])
data_to_use$m_psi <- sum(params_to_use)
data_to_use$keyfun <- models_to_run$keyfun[i]
model_fits[[i]] <- model_to_fit$sample(data = data_to_use,
chains = nc,
parallel_chains = nc,
init = 0.1,
max_treedepth = 10,
refresh = 25,
show_messages = TRUE,
show_exceptions = FALSE,
save_warmup = FALSE,
iter_sampling = ni,
iter_warmup = nw)
model_fits[[i]]$save_object(paste0("outputs/models/fit_multispecies_",
i,".rds"))
}Once all models have been run and saved, you can avoid having to re-run the above chunk by reading in the saved model objects, draws will also be read into memory and improve analysis speed:
Model Evaluation
Our strategy for model evaluation is to determine the best performing model out of a limited set of options with different fixed-effects, and random-effects. Model evaluation should consider formal diagnostics such as LOO (see below), as well as checks on the sensibility of predictions across unsampled locations.
LOO (Leave-One-Out 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 models. Better performing models according to loo::loo() will have higher elpd values.
Show code
loos <- list()
for(i in 1:length(model_fits)) {
loos[[i]] <- model_fits[[i]]$loo(variables = "log_lik2_species", cores = 8)
}
names(loos) <- paste("model", 1:length(model_fits))
loo_compare_table <- loo::loo_compare(x = loos)
# Plot loo table
gt(loo_compare_table %>%
as.data.frame() %>%
tibble::rownames_to_column()) %>%
fmt_number(decimals = 2) %>%
tab_style(locations = cells_row_groups(), style = cell_text(weight = "bold"))| elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic | |
|---|---|---|---|---|---|---|---|---|
| model 7 | 0.00 | 0.00 | 2,120.29 | 7.66 | 1.24 | 0.06 | −4,240.59 | 15.32 |
| model 8 | −0.09 | 0.30 | 2,120.20 | 7.64 | 1.19 | 0.05 | −4,240.41 | 15.29 |
| model 5 | −0.20 | 0.32 | 2,120.09 | 7.64 | 1.17 | 0.05 | −4,240.19 | 15.27 |
| model 6 | −1.76 | 0.68 | 2,118.53 | 7.66 | 1.04 | 0.05 | −4,237.06 | 15.32 |
| model 1 | −1.94 | 0.62 | 2,118.35 | 7.65 | 1.04 | 0.05 | −4,236.70 | 15.29 |
| model 2 | −2.16 | 0.68 | 2,118.13 | 7.65 | 1.00 | 0.05 | −4,236.26 | 15.30 |
| model 3 | −2.40 | 0.68 | 2,117.89 | 7.64 | 0.98 | 0.05 | −4,235.79 | 15.28 |
| model 4 | −4.37 | 0.92 | 2,115.93 | 7.69 | 0.85 | 0.04 | −4,231.86 | 15.38 |
Given that models are similar in there predictive performance at sampled sites we also evaluate the performance of predictions across unsampled sites. Below we see that more complex models (usually having a slighly better elpd), have poorer performance to unsampled areas and lead to unrealistic predictions.
Show code
Nhat_CV <- list()
for(i in 1:length(model_fits)) {
Nhat_CV[[i]] <- model_fits[[i]]$summary("Nhat", CV = ~ sd(.x)/mean(.x)) %>%
mutate(Model = i, Species = species_names)
}
Nhat_CV_c <- bind_rows(Nhat_CV)
Nhat_CV_c %>%
select(Model, Species, CV) %>%
kbl("html", escape=FALSE, caption = "Coefficients of Variation for Total Abundance estimates for the models assessed in the analysis") %>%
kable_styling() %>%
kableExtra::column_spec(3, background = ifelse(Nhat_CV_c$CV > 0.5, "red", ifelse(Nhat_CV_c$CV > 0.2, "orange", "green")), color = "white") %>%
collapse_rows(1) | Model | Species | CV |
|---|---|---|
| 1 | Sambar Deer | 0.1886410 |
| Fallow Deer | 0.9422867 | |
| Red Deer | 3.4679094 | |
| Hog Deer | 0.6694140 | |
| 2 | Sambar Deer | 0.1595965 |
| Fallow Deer | 0.4671161 | |
| Red Deer | 1.1285972 | |
| Hog Deer | 0.5351523 | |
| 3 | Sambar Deer | 0.1589921 |
| Fallow Deer | 0.4140389 | |
| Red Deer | 0.8865843 | |
| Hog Deer | 0.5135488 | |
| 4 | Sambar Deer | 0.1516583 |
| Fallow Deer | 0.2145688 | |
| Red Deer | 0.6266303 | |
| Hog Deer | 0.3045017 | |
| 5 | Sambar Deer | 0.2132886 |
| Fallow Deer | 0.5446794 | |
| Red Deer | 3.3479264 | |
| Hog Deer | 0.7870094 | |
| 6 | Sambar Deer | 0.2062165 |
| Fallow Deer | 0.2894896 | |
| Red Deer | 1.5069138 | |
| Hog Deer | 0.4926303 | |
| 7 | Sambar Deer | 0.2586611 |
| Fallow Deer | 0.5460739 | |
| Red Deer | 4.9932521 | |
| Hog Deer | 0.9147719 | |
| 8 | Sambar Deer | 0.2183291 |
| Fallow Deer | 0.4281970 | |
| Red Deer | 4.1957471 | |
| Hog Deer | 1.0076050 |
Selecting the ‘top’ model
Based on the findings of loo and comparisons of other checks (posterior predictive checks, variation in posterior draws, plausibility of predictions), we see that all model perform similarly using loo. However, several models produce predictions with lower levels of variation in the posterior predictions. Based on this, we use model 3 hereafter for analysis. This model has fixed effects of scale(sqrt(BareSoil)), scale(sqrt(Nitrogen)), scale(log(PastureDistance)), scale(BIO15) and scale(sqrt(ForestEdge)). It also has a random effect of bioregion for each species. This model does not have as many fixed effects as models 1 and 2, and also does not have a second random effect of EVC alongside bioregion (model 6). Unfortunately, the additional random effect of EVC produces much more variable and over-dispersed results for Red Deer, so we opt to use model 3 instead.
We acknowledge that this model selection process is not solely derived from analytical ‘scores’ or rankings (e.g. LOO or AIC), however, given the multispecies approach, negative binomial distribution and large amount of area being predicted to model selection should consider the predictive performance of many species and try and avoid overfitting and excessive dispersion in predictions. We acknowledge that a single model formula may not be optimal for all species; however we believe the benefit gained via a multispecies approach outweigh the constraints of the multispecies approach.
Note that model exploration and tuning was much more extensive than presented here, several months of model revisions, where many alterations were tried led us to this small selection of candidate models. For context previous exploration investigated (but later abandoned the use of):
- Gaussian process spatial random effects
- Site-level random effects
- Poisson distributions
- Zero-inflation models
- Hurdle models
- Fixed-effect only models
- Single-species only models
- Models with half-normal detection functions
- Models with other fixed-effect parameters not considered here
Show code
# assign index for top model
top <- 3Posterior predictive checks
Posterior predictive checks allow us to compare the observed data to the model-generated data. For each species we undertake posterior predictive checks for summary statistics relating to the number of deer seen on the cameras at each site. Ideally a well-fit model is able to make predictions that match the observed data. Here counts are the number of snapshot moments with deer. The summary statistics we use for the posterior predictive checks are:
- Maximum counts of deer seen at a site on a camera
- Mean number of deer seen at a site on a camera
- Standard deviation of the counts of deer on the cameras
- Average (mean) counts at sites in a scatter plot
- Proportion of sites with zero counts (camera or transect)
- Observed vs expected proportion of counts
When inspecting these PPC’s we should also keep in mind that the abundance parameter (\(\lambda\)) is not solely derived from camera counts, but the RN transect integrated model as well. Thus, in some cases we would expect more variation in the predicted vs expected values. For instance, mean counts should often be predicted to be higher than the ‘expected’ values given that the model is taking into account the latent abundance estimated from the transects.
Sambar Deer
Based on the plot below, we see that the PPCs for Sambar perform sufficiently well.
Show code
q95 <- function(x) quantile(x, 0.95, na.rm = T)
# q25 <- function(x) quantile(x, 0.25, na.rm = T)
q75 <- function(x) quantile(x, 0.75, na.rm = T)
sd_90 <- function (x, na.rm = FALSE) {
quants <- quantile(x, c(0.05, 0.95), na.rm = T)
x <- x[x < quants[2] & x > quants[1]]
sqrt(var(if (is.vector(x) || is.factor(x)) x else as.double(x),
na.rm = na.rm))
}
bayes_theme <- c(delwp_cols_seq$Teal[6:8], delwp_cols_shades$Navy[1:3])
bayesplot::color_scheme_set(bayes_theme)
funs <- c(max, mean, sd, ppc_scatter_avg, prop_zero, ppc_dens_overlay)
titles <- c("Max", "Mean", "SD", "Average Site Counts", "Proportion Zeros", "Density of Counts")
ppc_sambar <- list()
# funs <- c(prop_zero, mean, q90, sd)
for(i in 1:length(funs)) {
ppc_sambar[[i]] <- posterior_checks_multispecies(model = model_fits[[top]],
model_data = multispecies_data, species_index = 1,
stat = funs[[i]], integrated = T, only_det = F,
title = titles[i])
}
cowplot::plot_grid(plotlist = ppc_sambar, labels = "AUTO", ncol = 2)
Figure 1: Posterior predictive checks for Sambar Deer
Fallow Deer
The PPC’s below show that fallow has more dispersion in camera counts than Sambar. While the observed counts fall within the intervals of our model, there is higher uncertainty. Some of these findings appear to be driven by a site with very high counts (over 2,000 snapshot moments of deer individuals), this is more than twice the second highest.
Show code
ppc_fallow <- list()
for(i in 1:length(funs)) {
ppc_fallow[[i]] <- posterior_checks_multispecies(model = model_fits[[top]],
model_data = multispecies_data, species_index = 2,
stat = funs[[i]], integrated = T, only_det = F,
title = titles[i])
}
cowplot::plot_grid(plotlist = ppc_fallow, labels = "AUTO", ncol = 2)
Figure 2: Posterior predictive checks for Fallow Deer
Red Deer
Red Deer have PPC’s that show a high congruence between predicted and expected counts.
Show code
Figure 3: Posterior predictive checks for Red Deer
Hog Deer
Hog Deer have PPC’s that show our model performs sufficiently well in explaining the patterns of observed counts.
Show code
Figure 4: Posterior predictive checks for Hog Deer
Model Convergence
We observed no divergences or any STAN sampling warnings/issues for our models. To visualise the convergence of the model we can observe the mixing of chains for key parameters below:
beta_det: Parameters for the distance sampling model (intercept and herbaceous understorey)
beta_psi: Parameters for the fixed effect abundance model
bioregion_sd: Parameter for the bioregion random effect standard deviation
odRN: Parameter for the over-dispersion in the Royle-Nichols negative binomial model
od_mu: Parameters for the mean over-dispersion for each species for the camera counts (combined with odRN).
Show code
bayesplot::color_scheme_set(unname(delwp_cols[1:6]))
convergence_params <- model_fits[[top]]$draws(c("beta_det",
"beta_psi",
"bioregion_sd",
"beta_trans_det",
"odRN",
"od_mu"))
mcmc_trace(convergence_params, facet_args = list(ncol = 4)) +
theme(panel.spacing = unit(0.5, "lines"))
Figure 5: Chain mixing of several key parameters relating to detection, abundance and dispersion
Model Predictions
Within the STAN model we generate predictions for sampled and un-sampled locations. This provides us with site-level abundance estimates as well as estimates across all (un-sampled) public forest.
Site-based Predictions
Site-based Detections Map
Below we display a detection/non-detection map for all sites. Detections are recorded from either camera or transect detections.
Show code
vic_regions <- vicmap_query("open-data-platform:delwp_region") %>%
collect() %>%
st_transform(3111) %>%
st_simplify(dTolerance = 500)
delwp_pal <- colorRampPalette(c("#B9C600",
delwp_cols[["Teal"]],
delwp_cols[["Navy"]]))
site_detection_plot <- function(data, regions, cams_curated, species_index) {
cam_data <- cams_curated %>% dplyr::select(SiteID)
cam_data$Detected <- factor(data$any_seen[species_index, ])
cam_data$cam_seen <- factor(as.integer(as.logical(rowSums(data$n_obs[,,species_index]))))
if(species_index == 4) {
regions <- regions %>% filter(delwp_region == "GIPPSLAND")
cam_data <- cam_data %>% st_filter(regions %>% st_transform(4283))
}
lvls <- length(unique(cam_data$Detected, na.rm = T))
plot <- ggplot2::ggplot(data = cam_data %>% st_transform(3111)) +
ggplot2::geom_sf(data = regions, alpha = 0.75, fill = "grey80") +
ggplot2::geom_sf(aes(fill = Detected), shape = 21, size = 2, alpha = 0.75) +
ggplot2::scale_fill_manual(values = delwp_pal(lvls),
na.value = "transparent", na.translate = F,
name = "Detected", guide = guide_legend(override.aes = list(size = 6))) +
# ggtitle(species_names[species_index])+
theme_bw() +
theme(legend.text = element_text(size = 18), legend.key.size = unit(1, "cm"),
title = element_text(size = 20))
return(plot)
}
detection_plot <- cowplot::plot_grid(site_detection_plot(multispecies_data, regions = vic_regions, cams_curated = cams_curated, species_index = 1),
site_detection_plot(multispecies_data, regions = vic_regions, cams_curated = cams_curated, species_index = 2),
site_detection_plot(multispecies_data, regions = vic_regions, cams_curated = cams_curated, species_index = 3),
site_detection_plot(multispecies_data, regions = vic_regions, cams_curated = cams_curated, species_index = 4), ncol = 2,
labels = c("A) Sambar deer", "B) Fallow deer", "C) Red deer", "D) Hog deer"), greedy = TRUE)
ggsave(plot = detection_plot, filename = "outputs/plots/all_detection_plot.png", width = 4800, height = 3200, units = "px", bg = "white")
detection_plot
Figure 6: Detections of deer (camera or transects) from across Victoria. Detections of Hog Deer are restricted to Gippsland.
Site-based Prediction Map
Our model generates predictions of density (per km2) at each site using the \(\lambda\) value at that site, average group size at that site, and the dispersion at that site (varies across species). For this quantity (N_site), we plot visualisations of mean estimates for the various deer species surveyed for in this study. We save the numerical values, as a spatial object (sf) at outputs/site_deer_predictions.rds.
Show code
site_preds <- function(model, cams_curated, species_index) {
rn_dens <- model$summary("N_site", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), CV = ~sd(.)/mean(.), posterior::default_summary_measures(), posterior::default_convergence_measures())
which_sp <- which(stringr::str_detect(string = rn_dens$variable,
pattern = paste0("N_site\\[",
species_index)))
density_at_sites_rn <- cbind(cams_curated, rn_dens[which_sp, ])
return(density_at_sites_rn)
}
site_density_plot <- function(densities, regions, species) {
densities$density <- cut(densities$mean,
breaks = c(0, 0.25, 0.5, 1, 3, 5, 10, max(densities$mean)),
labels = c("< 0.25",
"0.25 - 0.5",
"0.5 - 1",
"1 - 3",
"3 - 5",
"5 - 10",
"10 +"), include.lowest = T, right = T)
if(species == "Hog") {
regions <- regions %>% filter(delwp_region == "GIPPSLAND")
densities <- densities %>% st_filter(regions %>% st_transform(4283))
}
lvls <- length(unique(densities$density, na.rm = T))
plot <- ggplot2::ggplot(data = densities) +
ggplot2::geom_sf(data = regions, alpha = 0.75, fill = "grey80") +
ggplot2::geom_sf(aes(fill = density, alpha = mean), shape = 21, size = 3) +
ggplot2::scale_fill_manual(values = delwp_pal(lvls),
na.value = "transparent", na.translate = F,
name = "", guide = guide_legend(override.aes = list(size = 6))) +
scale_alpha_continuous(range = c(0.5,1), guide = "none") +
labs(title = paste0('Abundance of ',species ,' deer'), fill = bquote('Deer per'~km^2)) +
theme_bw() +
theme(legend.text = element_text(size = 18), legend.key.size = unit(1, "cm"),
title = element_text(size = 20))
return(plot)
}
sambar_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 1)
fallow_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 2)
red_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 3)
hog_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 4)
site_deer_predictions <- list(Sambar = sambar_preds,
Fallow = fallow_preds,
Red = red_preds,
Hog = hog_preds)
saveRDS(site_deer_predictions, "outputs/site_deer_predictions.rds")
cowplot::plot_grid(site_density_plot(sambar_preds, vic_regions, species = "Sambar"),
site_density_plot(fallow_preds, vic_regions, species = "Fallow"),
site_density_plot(red_preds, vic_regions, species = "Red"),
site_density_plot(hog_preds, vic_regions, species = "Hog"), ncol = 1)
Figure 7: Site-level density estimates for all sites sampled as part of statewide and hog deer surveys. Point-density estimates are from trimmed means of the posterior draws (trimming of the top and bottom 2.5% of draws)
Average density estimates varied across sites for all four species. Sambar deer were the most commonly detected species on the camera traps (104/317 sites), and density estimates at sites ranged from 0 to 28.31 deer/km2 (mean = 1.94, sd = 2.76), for fallow, fewer site detections (30), but more variability in site density gave a range between 0 and 23.25 deer/km2 (mean = 0.66, sd = 1.54). Red deer were only detected on 12 sites with densities ranging from 0 to 6.81 (mean = 0.12, sd = 0.53). Finally, Hog Deer, which were only detected on 22 cameras in South Gippsland, had site densities ranging from 0 to 16.9 deer/km2 (mean = 0.48, sd = 1.49).
Site-density Summaries
As a check on model expectations we compare the average modelled densities at sites with (i) no evidence of deer, (ii) evidence of deer present on only transects, (iii) evidence of deer present on only the cameras, and (iv) evidence of deer on both cameras and transect. We expect that average densities are generally higher at sites that detected some form of deer than sites that did not detect any sign of deer. Additionally, we would also expect average densities to be generally higher at sites that had detections on cameras than those with only detections from transects. The table below shows these expectations to be correct.
Show code
density_summary_table <- function(preds, model_data, species, species_index) {
cam_seen <- as.integer(as.logical(rowSums(model_data$n_obs[,,species_index])))
transect_seen <- model_data$transects %>%
filter(Species == !!deer_species_all[species_index]) %>%
mutate(NoCam = case_when(Survey == "Camera" ~ 0L,
TRUE ~ 1)) %>%
group_by(SiteID) %>%
summarise(TransectDet = max(NoCam*Presence)) %>%
pull(TransectDet)
preds_sum <- preds %>%
st_drop_geometry() %>%
mutate(`Species` = species,
`CameraDetection` = cam_seen,
`AnyDetection` = model_data$any_seen[species_index, ],
`OnlyTransect` = transect_seen,
`Detection` = case_when(CameraDetection == 0 & AnyDetection == 0 ~ "Not seen",
CameraDetection == 0 & AnyDetection == 1 ~ "Only detected on transects",
CameraDetection == 1 & OnlyTransect == 0 ~ "Only seen on cameras",
CameraDetection == 1 & OnlyTransect == 1 ~ "Seen on both camera and transect")) %>%
group_by(`Species`, `Detection`) %>%
summarise(`Number of Sites` = n(),
`Average Density` = round(mean(mean), 2),
`SD between sites` = round(sd(mean), 2)) %>%
ungroup()
return(preds_sum)
}
density_summary <- bind_rows(
density_summary_table(sambar_preds, multispecies_data,
species = "Sambar", species_index = 1),
density_summary_table(fallow_preds, multispecies_data,
species = "Fallow", species_index = 2),
density_summary_table(red_preds, multispecies_data,
species = "Red", species_index = 3),
density_summary_table(hog_preds, multispecies_data,
species = "Hog", species_index = 4))
readr::write_csv(x = density_summary, file = "outputs/tables/site_group_estimates.csv")
density_summary %>%
kbl(format = "html", digits = 2, caption = "Average (mean) density estimates at the various groups of sites (according to detections)") %>%
kable_styling("striped") %>%
collapse_rows(1)| Species | Detection | Number of Sites | Average Density | SD between sites |
|---|---|---|---|---|
| Sambar | Not seen | 182 | 0.81 | 1.29 |
| Only detected on transects | 31 | 2.46 | 1.71 | |
| Only seen on cameras | 50 | 3.34 | 4.71 | |
| Seen on both camera and transect | 54 | 4.16 | 2.51 | |
| Fallow | Not seen | 259 | 0.52 | 1.53 |
| Only detected on transects | 28 | 0.97 | 1.43 | |
| Only seen on cameras | 20 | 1.37 | 1.38 | |
| Seen on both camera and transect | 10 | 1.89 | 1.62 | |
| Red | Not seen | 304 | 0.05 | 0.19 |
| Only detected on transects | 1 | 1.23 | NA | |
| Only seen on cameras | 5 | 1.47 | 1.48 | |
| Seen on both camera and transect | 7 | 2.01 | 2.37 | |
| Hog | Not seen | 293 | 0.27 | 1.18 |
| Only detected on transects | 2 | 2.03 | 2.71 | |
| Only seen on cameras | 22 | 3.10 | 2.39 |
Regional and Statewide Abundance
Within our model we calculate abundance/density for deer in each of the 74,570 km2 of public land. Based on these spatial predictions we can estimate abundance at a regional level (6 DEECA regions) and across the whole state.
Statewide Maps
Using the model predictions (pred) for all suitable public land we generate a raster (1km2 resolution). We save the average spatial estimates under outputs/rasters (.tif files) and also provide binned plots (png files) in outputs/plots.
Show code
pred_raster_full <- terra::rast(prediction_raster)
pred_raster <- terra::app(pred_raster_full[[stringr::str_subset(
stringr::str_remove_all(labels(terms(ab_formula_1)),
"scale[(]|[)]|log[(]|sqrt[(]"),
pattern = "[*]", negate = T)]], mean)
mean_raster <- list()
sd_raster <- list()
mean_raster_discrete <- list()
# gp_preds_draws_all <- model_fits[[top]]$draws("pred", format = "matrix")
gp_preds_draws_all <- model_fits[[top]]$draws("pred", format = "matrix")
# gp_preds_summary_all <- model_fits[[top]]$summary("pred", prob_occ = ~ mean(. > 0))
for(i in 1:length(deer_species_all)) {
which_sp <- which(stringr::str_detect(string = colnames(gp_preds_draws_all),
pattern = paste0("pred\\[", i)))
gp_preds_draws_sp <- gp_preds_draws_all[,which_sp]
terra::values(pred_raster)[!is.na(terra::values(pred_raster))] <- apply(gp_preds_draws_sp,
2,
mean,
na.rm = T) #filter out highest draw
mean_raster[[i]] <- pred_raster
mean_raster_discrete[[i]] <- mean_raster[[i]]
max_pred <- max(values(mean_raster[[i]]), na.rm = T)
values(mean_raster_discrete[[i]]) <- cut(values(mean_raster_discrete[[i]]) ,
breaks = c(0,0.25,0.5,1,3,5,10,max_pred),
include.lowest = T, right = T,
labels = c("< 0.25",
"0.25 - 0.5",
"0.5 - 1",
"1 - 3",
"3 - 5",
"5 - 10",
"10 +"))
# occ_raster[[i]] <- pred_raster
# terra::values(occ_raster[[i]])[!is.na(terra::values(occ_raster[[i]]))] <- apply(gp_preds_draws_sp, 2, function(x) sum(x > 0)/length(x))
sd_raster[[i]] <- pred_raster
terra::values(sd_raster[[i]])[!is.na(terra::values(sd_raster[[i]]))] <- apply(gp_preds_draws_sp, 2, sd)
}
# combine mean rasters together
combined_raster <- rast(mean_raster)
names(combined_raster) <- c("Average Sambar Deer Density (per km2)",
"Average Fallow Deer Density (per km2)",
"Average Red Deer Density (per km2)",
"Average Hog Deer Density (per km2)")
combined_sd_raster <- rast(sd_raster)
names(combined_sd_raster) <- c("Sambar Deer SD",
"Fallow Deer SD",
"Red Deer SD",
"Hog Deer SD")
# combined_occupancy_raster <- rast(occ_raster)
# names(combined_occupancy_raster) <- c("Average Sambar Deer Occupancy Probability",
# "Average Fallow Deer Occupancy Probability",
# "Average Red Deer Occupancy Probability",
# "Average Hog Deer Occupancy Probability")
# writeRaster(combined_occupancy_raster, "outputs/rasters/combined_occupancy_raster.tif", overwrite = T)
writeRaster(combined_raster, "outputs/rasters/combined_deer_average_density.tif", overwrite = T)
writeRaster(combined_sd_raster, "outputs/rasters/combined_deer_sd.tif", overwrite = T)Show code
# reg <- VicmapR::vicmap_query("open-data-platform:delwp_region") %>%
# VicmapR::collect() %>%
# sf::st_transform(3111) %>%
# sf::st_simplify(dTolerance = 250)
state <- VicmapR::vicmap_query("open-data-platform:vmlite_victoria_polygon_su5") %>%
filter(state != "NSW" & state != "SA" & feature_type_code != "sea") %>%
VicmapR::collect() %>%
sf::st_transform(3111)
gippsland <- vic_regions %>% filter(delwp_region == "GIPPSLAND")
plot_abundance <- function(raster, state, species, crop = NULL) {
if(!is.null(crop)) {
raster <- terra::crop(raster, vect(crop), mask = T)
state <- crop
}
lvls <- length(unique(values(raster, na.rm = T)))
plot <- ggplot2::ggplot() +
ggplot2::geom_sf(data = state, alpha = 1, linewidth = 0.5, fill = "grey90") +
tidyterra::geom_spatraster(data = raster, na.rm = T) +
# tidyterra::scale_fill_terrain_d(na.translate = FALSE) +
ggplot2::scale_fill_manual(values = delwp_pal(lvls), na.value = "transparent", na.translate = F) +
ggplot2::labs(fill = bquote('Deer per'~km^2)) +
ggspatial::annotation_scale() +
delwp_theme()
return(plot)
}
sambar_abundance_plot <- plot_abundance(mean_raster_discrete[[1]],
state = state,
species = "Sambar Deer")
fallow_abundance_plot <- plot_abundance(mean_raster_discrete[[2]],
state = state,
species = "Fallow Deer")
red_abundance_plot <- plot_abundance(mean_raster_discrete[[3]],
state = state,
species = "Red Deer")
hog_abundance_plot <- plot_abundance(mean_raster_discrete[[4]],
state = state, crop = gippsland,
species = "Hog Deer")
ggsave(plot = sambar_abundance_plot, filename = "outputs/plots/sambar_abundance_plot.png",
width = 2400, height = 1600, units = "px")
ggsave(plot = fallow_abundance_plot, filename = "outputs/plots/fallow_abundance_plot.png",
width = 2400, height = 1600, units = "px")
ggsave(plot = red_abundance_plot, filename = "outputs/plots/red_abundance_plot.png",
width = 2400, height = 1600, units = "px")
ggsave(plot = hog_abundance_plot, filename = "outputs/plots/hog_abundance_plot.png",
width = 2400, height = 1600, units = "px")
ggsave(plot = cowplot::plot_grid(sambar_abundance_plot, fallow_abundance_plot, red_abundance_plot, hog_abundance_plot, ncol = 2, labels = c("A) Sambar deer", "B) Fallow deer", "C) Red deer", "D) Hog deer"), greedy = F), filename = "outputs/plots/all_abundance_plot.png", width = 4800, height = 3200, units = "px", bg = "white")
cowplot::plot_grid(sambar_abundance_plot, fallow_abundance_plot, red_abundance_plot, hog_abundance_plot, ncol = 1, labels = c("A) Sambar deer", "B) Fallow deer", "C) Red deer", "D) Hog deer"))
Figure 8: Abundance of Sambar, Fallow, Red and Hog Deer across Victoria, dark-grey area reflects area not included in predictions (i.e. not public land)
Regional Abundance Estimates
For each DEECA region we provide species-level estimates of abundance with 90 % confidence interval. We also calculate the model-based average density within each region based on the abundance and the total area of public land within each region.
Show code
Nhat_sum <- model_fits[[top]]$summary("Nhat", trimmed_mean = ~mean(., trim = 0.025), CV = ~sd(.)/mean(.), posterior::default_summary_measures(), posterior::default_convergence_measures())
Nhat_all <- model_fits[[top]]$summary("Nhat_sum", trimmed_mean = ~mean(., trim = 0.025), CV = ~sd(.)/mean(.), posterior::default_summary_measures(), posterior::default_convergence_measures())
format_nhat <- function(sp, data = Nhat_sum) {
fd <- data %>%
filter(variable == !!paste0("Nhat[", sp, "]"))
paste0("$\\hat N$ = ",
formatC(fd$trimmed_mean, digits = 0, big.mark = ",", format = "f"),
" [90% CI: ",
formatC(fd$q5, digits = 0, big.mark = ",", format = "f"),
" – ",
formatC(fd$q95, digits = 0, big.mark = ",", format = "f"),
"]")
}
reg <- VicmapR::vicmap_query("open-data-platform:delwp_region") %>%
VicmapR::collect() %>%
dplyr::group_by(delwp_region) %>%
dplyr::summarise(geometry = sf::st_combine(geometry)) %>%
dplyr::ungroup() %>%
dplyr::mutate(delwp_region_fact = as.integer(factor(delwp_region)))
abundance_table <- function(model, regions, pred_area, caption, species_index, return_data = F) {
reg <- model$summary("Nhat_reg", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), posterior::default_summary_measures(), posterior::default_convergence_measures())
which_sp <- which(stringr::str_detect(string = reg$variable,
pattern = paste0("Nhat_reg\\[", species_index)))
regional_abundance <- reg[which_sp,] %>%
dplyr::bind_rows(model$summary("Nhat", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), posterior::default_summary_measures(), posterior::default_convergence_measures())[species_index,]) %>%
dplyr::mutate(variable = c(regions$delwp_region, "TOTAL"),
`Area km2` = round(c(pred_area, sum(pred_area))),
`Density (mean)` = round(trimmed_mean/`Area km2`, 2),
`Density (5 %)` = round(q5/`Area km2`, 2),
`Density (95 %)` = round(q95/`Area km2`, 2),
`Density mean (90% CI)` = paste0(`Density (mean)`,
" (",
`Density (5 %)`,
", ",
`Density (95 %)`,
")"),
CV = round(sd/trimmed_mean, 2)) %>%
dplyr::transmute(Region = variable,
`Trimmed mean` = round(trimmed_mean),
# Median = round(median),
SD = round(sd),
CV,
# MAD = round(mad),
`5 %` = round(q5),
`95 %` = round(q95),
`Area km2`,
`Density mean (90% CI)`)
if(return_data) {
if(species_index == 4) {
regional_abundance <- regional_abundance %>%
filter(Region == "GIPPSLAND")
}
return(regional_abundance %>%
mutate(Species = species_names[species_index], across(2:8, ~ formatC(.x, digits = 0, format = "f", big.mark = ",") )) %>%
dplyr::select(Species, everything()))
}
kableExtra::kbl(regional_abundance, digits = 2, format = "html", caption = caption, format.args = list(big.mark = ",")) %>%
kableExtra::kable_styling("striped") %>%
kableExtra::column_spec(1, bold = TRUE) %>%
kableExtra::row_spec(6, hline_after = T) %>%
kableExtra::row_spec(7, background = "#A5A1B5", color = "black", bold = T, hline_after = T)
}
pred_a_df <- data.frame(pred_reg = multispecies_data$pred_reg,
area = multispecies_data$prop_pred) %>%
group_by(pred_reg) %>%
summarise(area = sum(area))
# Save Abundance Table
all_ab_tab <- bind_rows(abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, species_index = 1, return_data = T),
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, species_index = 2, return_data = T),
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, species_index = 3, return_data = T),
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, species_index = 4, return_data = T))
write.csv(all_ab_tab, "outputs/tables/regional_abundance_estimates.csv")
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Sambar Deer Density", species_index = 1)| Region | Trimmed mean | SD | CV | 5 % | 95 % | Area km2 | Density mean (90% CI) |
|---|---|---|---|---|---|---|---|
| BARWON SOUTH WEST | 139 | 277 | 1.99 | 35 | 430 | 4,745 | 0.03 (0.01, 0.09) |
| GIPPSLAND | 54,895 | 8,898 | 0.16 | 42,618 | 70,948 | 25,211 | 2.18 (1.69, 2.81) |
| GRAMPIANS | 1,619 | 2,032 | 1.25 | 716 | 3,266 | 9,896 | 0.16 (0.07, 0.33) |
| HUME | 59,732 | 12,360 | 0.21 | 43,334 | 82,188 | 16,890 | 3.54 (2.57, 4.87) |
| LODDON MALLEE | 1,094 | 1,761 | 1.61 | 460 | 2,589 | 15,597 | 0.07 (0.03, 0.17) |
| PORT PHILLIP | 5,200 | 1,334 | 0.26 | 3,367 | 7,694 | 2,231 | 2.33 (1.51, 3.45) |
| TOTAL | 123,061 | 19,657 | 0.16 | 96,200 | 157,638 | 74,570 | 1.65 (1.29, 2.11) |
Show code
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Fallow Deer Density", species_index = 2)| Region | Trimmed mean | SD | CV | 5 % | 95 % | Area km2 | Density mean (90% CI) |
|---|---|---|---|---|---|---|---|
| BARWON SOUTH WEST | 10,514 | 14,091 | 1.34 | 3,049 | 31,646 | 4,745 | 2.22 (0.64, 6.67) |
| GIPPSLAND | 12,305 | 5,001 | 0.41 | 6,914 | 21,794 | 25,211 | 0.49 (0.27, 0.86) |
| GRAMPIANS | 4,350 | 3,203 | 0.74 | 2,106 | 9,264 | 9,896 | 0.44 (0.21, 0.94) |
| HUME | 16,021 | 5,519 | 0.34 | 9,362 | 26,721 | 16,890 | 0.95 (0.55, 1.58) |
| LODDON MALLEE | 3,296 | 2,198 | 0.67 | 1,384 | 7,277 | 15,597 | 0.21 (0.09, 0.47) |
| PORT PHILLIP | 1,501 | 767 | 0.51 | 761 | 2,879 | 2,231 | 0.67 (0.34, 1.29) |
| TOTAL | 48,932 | 20,900 | 0.43 | 29,888 | 85,063 | 74,570 | 0.66 (0.4, 1.14) |
Show code
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Red Deer Density", species_index = 3)| Region | Trimmed mean | SD | CV | 5 % | 95 % | Area km2 | Density mean (90% CI) |
|---|---|---|---|---|---|---|---|
| BARWON SOUTH WEST | 5,645 | 8,018 | 1.42 | 1,559 | 19,399 | 4,745 | 1.19 (0.33, 4.09) |
| GIPPSLAND | 179 | 196 | 1.10 | 49 | 502 | 25,211 | 0.01 (0, 0.02) |
| GRAMPIANS | 5,134 | 5,709 | 1.11 | 1,705 | 14,695 | 9,896 | 0.52 (0.17, 1.48) |
| HUME | 1,126 | 739 | 0.66 | 420 | 2,638 | 16,890 | 0.07 (0.02, 0.16) |
| LODDON MALLEE | 103 | 1,524 | 14.83 | 1 | 794 | 15,597 | 0.01 (0, 0.05) |
| PORT PHILLIP | 31 | 83 | 2.67 | 5 | 103 | 2,231 | 0.01 (0, 0.05) |
| TOTAL | 12,672 | 12,275 | 0.97 | 4,719 | 35,465 | 74,570 | 0.17 (0.06, 0.48) |
Show code
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Hog Deer Density", species_index = 4) | Region | Trimmed mean | SD | CV | 5 % | 95 % | Area km2 | Density mean (90% CI) |
|---|---|---|---|---|---|---|---|
| BARWON SOUTH WEST | 0 | 0 | NaN | 0 | 0 | 4,745 | 0 (0, 0) |
| GIPPSLAND | 4,243 | 2,263 | 0.53 | 2,121 | 8,464 | 25,211 | 0.17 (0.08, 0.34) |
| GRAMPIANS | 0 | 0 | NaN | 0 | 0 | 9,896 | 0 (0, 0) |
| HUME | 0 | 0 | NaN | 0 | 0 | 16,890 | 0 (0, 0) |
| LODDON MALLEE | 0 | 0 | NaN | 0 | 0 | 15,597 | 0 (0, 0) |
| PORT PHILLIP | 0 | 0 | NaN | 0 | 0 | 2,231 | 0 (0, 0) |
| TOTAL | 4,243 | 2,263 | 0.53 | 2,121 | 8,464 | 74,570 | 0.06 (0.03, 0.11) |
Deer were found to be widely distributed across Victoria. Across approximately 74,570 km2 of public land we estimate total deer abundance of the four species investigated in this study (Sambar, Fallow, Red and Hog) at 191,153 (90 % CI: 146,732 - 255,490. Sambar deer were the most populous species across Victoria (\(\hat N\) = 123,061 [90% CI: 96,200 – 157,638]), followed by Fallow at \(\hat N\) = 48,932 [90% CI: 29,888 – 85,063], Red (\(\hat N\) = 12,672 [90% CI: 4,719 – 35,465]), and Hog (\(\hat N\) = 4,243 [90% CI: 2,121 – 8,464]).
Deer Abundance in High Conservation Areas
Below we show code for extracting average abundance and densities for areas listed under the National Parks Act 1975, and greater than 20km2.
Show code
parkres <- vicmap_query("parkres") %>%
collect(quiet = TRUE) %>%
st_transform(3111) %>%
st_filter(vic_regions)
parkres_grouped <- parkres %>%
group_by(name, name_short) %>%
summarise(area = sum(areasqm)/1e6,
geometry = st_combine(geometry)) %>%
filter(area > 20) %>%
ungroup()
ab_e <- terra::extract(combined_raster,
terra::vect(parkres_grouped),
fun = "sum", method = "simple", na.rm = T, touches = TRUE) %>%
`names<-`(c("ID", species_names))
# get abundance estimates
park_ab <- bind_cols(parkres_grouped, round(ab_e[,2:5],0)) %>%
na.omit() %>%
mutate(`All Deer` = `Sambar Deer` + `Fallow Deer` + `Red Deer` + `Hog Deer`)%>%
mutate(across(contains("Deer"), ~ round(.x/area,2), .names = "Density {.col}")) %>%
mutate(Sambar = paste0(formatC(`Sambar Deer`, format = "fg", big.mark = ","), " (", `Density Sambar Deer`, ")"),
Fallow = paste0(formatC(`Fallow Deer`, format = "fg", big.mark = ","), " (", `Density Fallow Deer`, ")"),
Red = paste0(formatC(`Red Deer`, format = "fg", big.mark = ","), " (", `Density Red Deer`, ")"),
Hog = paste0(formatC(`Hog Deer`, format = "fg", big.mark = ","), " (", `Density Hog Deer`, ")"),
All = paste0(formatC(`All Deer`, format = "fg", big.mark = ","), " (", `Density All Deer`, ")")) %>%
arrange(desc(`Density All Deer`))
park_averages <- park_ab %>%
mutate(name = paste0(name, " (", round(area, 0), "km2)")) %>%
dplyr::select(Reserve = name,
Sambar, Fallow, Red, Hog, All) %>%
st_drop_geometry()
# write.csv(park_averages, "outputs/tables/park_averages.csv", row.names = F)Statewide Distributions
Distribution from an abundance threshold
From our abundance estimates we can obtain estimates of the distribution of the deer species across Victoria.
Show code
inv.cloglog <- function(x) 1-exp(-x)
detections <- function(data, regions, cams_curated, species_index) {
cam_data <- cams_curated %>% dplyr::select(SiteID)
cam_data$Detected <- factor(data$any_seen[species_index, ])
cam_data$cam_seen <- factor(as.integer(as.logical(rowSums(data$n_obs[,,species_index]))))
# if(species_index == 4) {
# regions <- regions %>% filter(delwp_region == "GIPPSLAND")
# cam_data <- cam_data %>% st_filter(regions %>% st_transform(4283))
# }
return(cam_data)
}
joint_threshold <- list()
combined_threshold <- list()
for(i in 1:length(species_names)) {
dets <- detections(multispecies_data,
vic_regions,
cams_curated = cams_curated,
species_index = i) %>%
st_make_valid() %>%
st_transform(3111)
pointVals <- terra::extract(combined_raster[[i]],
terra::vect(dets %>% st_buffer(1000)), fun = "mean", na.rm = T)[,2]
rocData <- dets[which(!is.na(pointVals)),] %>%
mutate(lambda = na.omit(pointVals),
occ = inv.cloglog(lambda))
roc <- roc(rocData$Detected, rocData$lambda)
ci_coords <- ci.coords(roc, "best",
best.policy = "random",
best.method = "closest.topleft",
progress = "none",
conf.level = 0.9)
threshold_rast <- combined_raster[[i]]
threshold_rast_lci <- combined_raster[[i]]
threshold_rast_uci <- combined_raster[[i]]
threshold_rast[values(threshold_rast) < ci_coords$threshold[2]] <- 0
threshold_rast[values(threshold_rast) >= ci_coords$threshold[2]] <- 1
threshold_rast_lci[values(threshold_rast_lci) < ci_coords$threshold[3]] <- 0
threshold_rast_lci[values(threshold_rast_lci) >= ci_coords$threshold[3]] <- 1
threshold_rast_uci[values(threshold_rast_uci) < ci_coords$threshold[1]] <- 0
threshold_rast_uci[values(threshold_rast_uci) >= ci_coords$threshold[1]] <- 1
combined_threshold[[i]] <- c(threshold_rast, threshold_rast_lci, threshold_rast_uci) %>% `names<-`(c("Threshold Median", "Threshold LCI", "Threshold UCI"))
joint_threshold[[i]] <- app(combined_threshold[[i]], sum, na.rm = T)
d <- data.frame(id = 0:3, values = c(NA_character_, "UCI", "Median", "LCI"))
levels(joint_threshold[[i]]) <- d
}
joint_threshold <- rast(joint_threshold)
names(combined_threshold) <- species_names
names(joint_threshold) <- species_names
writeRaster(joint_threshold, "outputs/rasters/combined_threshold_raster.tif", overwrite = TRUE)
get_th_area <- function(stack) {
med <- sum(values(stack[[1]], na.rm = T) * multispecies_data$prop_pred)
LCI <- sum(values(stack[[2]], na.rm = T) * multispecies_data$prop_pred)
UCI <- sum(values(stack[[3]], na.rm = T) * multispecies_data$prop_pred)
paste0(formatC(med, format = "fg", big.mark = ","), "km^2^ [95% CI: ",
formatC(LCI, format = "fg", big.mark = ","), ", ",
formatC(UCI, format = "fg", big.mark = ","), "]")
}
plot_threshold <- function(raster, state, crop = NULL) {
if(!is.null(crop)) {
raster <- terra::crop(raster, vect(crop), mask = T)
state <- crop
}
# lvls <- length(unique(values(raster, na.rm = T)))
plot <- ggplot2::ggplot() +
ggplot2::geom_sf(data = state, alpha = 1, linewidth = 0.5, fill = "grey98") +
tidyterra::geom_spatraster(data = raster, aes(alpha = after_stat(value),
fill = after_stat(value)), na.rm = T) +
# tidyterra::scale_fill_terrain_d(na.translate = FALSE) +
ggplot2::scale_fill_manual(values = unname(rep(delwp_cols[2], 3)),
na.value = "transparent",
na.translate = F) +
ggplot2::scale_alpha_manual(values = c(0.25, 0.65, 0.85),
na.value = "transparent",
na.translate = F) +
ggplot2::labs(alpha = "Distribution \nrange", fill = "Distribution \nrange") +
ggspatial::annotation_scale() +
delwp_theme()
return(plot)
}
sambar_th <- plot_threshold(joint_threshold[[1]],
state = state)
fallow_th <- plot_threshold(joint_threshold[[2]],
state = state)
red_th <- plot_threshold(joint_threshold[[3]],
state = state)
hog_th <- plot_threshold(joint_threshold[[4]],
state = state, crop = gippsland)
threshold_combined_plot <- cowplot::plot_grid(sambar_th, fallow_th, red_th, hog_th, ncol = 2, labels = c("A) Sambar deer", "B) Fallow deer", "C) Red deer", "D) Hog deer"), greedy = F)
ggsave(plot = threshold_combined_plot, filename = "outputs/plots/all_threshold_plot.png", width = 4800, height = 3200, units = "px", bg = "white")
threshold_combined_plot
Figure 9: Threshold distributions of Sambar, Fallow, Red and Hog Deer across Victoria, white area reflects area not included in predictions (i.e. not public land), LCI and UCI confidence bounds refer to bootstrapped threshold confidence intervals (90%) for optimal sensitivity and specificity.
We obtained a threshold distribution based on the mean abundance values for each grid cell. Our range estimates suggest Sambar Deer have a range of 38,582km2 [95% CI: 33,605, 40,590]. For Fallow deer the range estimate has more variability 25,151km2 [95% CI: 20,666, 30,458], this is similarly the case for Red Deer (9,403km2 [95% CI: 2,465, 21,985]). Lastly, Hog Deer have a more isolated range estimate (2,176km2 [95% CI: 1,542, 2,219]).
Comparison of range size to existing distributions
Using range size maps based on presence-only records and VBA elicitation (Forsyth, Stamation, and Woodford 2015, 2016) we can generate estimated existing distributions for the deer species on public land:
Show code
existing_dist <- list()
existing_dist[["Sambar"]] <- st_read("data/existing_distributions/Sambar_Deer_Distributions_Final_20201203", quiet = TRUE)
existing_dist[["Fallow"]] <- st_read("data/existing_distributions/Fallow_Deer_Distributions_Final_20201203", quiet = TRUE)
existing_dist[["Red"]] <- st_read("data/existing_distributions/Red_Deer_Distributions_Final_20201203", quiet = TRUE)
existing_dist[["Hog"]] <- st_read("data/existing_distributions/Hog_Deer_Distribution_20160329", quiet = TRUE)
crown_land <- readRDS("data/crown_land.rds")
# Get overlapping CA
existing_dist_cl <- lapply(existing_dist, FUN = function(x) {
d <- st_transform(x, 3111)
st_intersection(d, crown_land) %>%
st_union()
})
existing_ranges <- sapply(existing_dist_cl, st_area) %>%
units::set_units("m2") %>%
units::set_units("km2") %>%
formatC(., digits = 0, format = "fg", big.mark = ",") %>%
as.character()
names(existing_ranges) <- species_names
existing_rangesSambar Deer Fallow Deer Red Deer Hog Deer
"43,779" "10,259" "2,630" "1,147" Model Covariates
Our model uses covariates to inform three separate sub-models:
- A distance-sampling detection process in front of the camera
- An imperfect detection/observation model based on the transect methods and the camera survey
- An abundance model based on spatially-derived covariates, as well as a bioregion random effect
The following section explores the effect of these covariates and sub-models on detection and abundance.
Detection Processes
At a site, there are three processes that may lead us to not record deer on a camera when in reality they occupy the site (whereby the site is the home range area surrounding the camera.) Firstly, deer may be available for detection and enter the sampling area in front of the camera. However, due to a function of their distance from the camera we fail to detect this individual. Secondly, they may occupy a site (known from transect signs) but never enter the camera field of view, in this case we estimate the various detection rates from the various methods. The detection rate of the camera in this method is regarded as the spatial availability of deer for camera trap sampling. Thirdly, we take into account temporal availability via activity time.
Distance-sampling
Below we show the average detection rate for a given group of deer in front of the camera (up to 12.5m).
Show code
det_rates <- model_fits[[top]]$summary("p") %>%
mutate(var = stringr::str_extract(variable, "(?<=\\[).+?(?=\\])")) %>%
separate(var, into = c("site", "Group Size"), sep = ",")
av_det_rates <- det_rates %>%
group_by(`Group Size`) %>%
summarise(mean = mean(mean)) %>%
ungroup() %>%
transmute(`Group Size`,
`Average Site Detection Probability` = mean)
av_det_rates %>%
bind_rows() %>%
arrange(`Group Size`) %>%
kbl("html", digits = 3, caption = "Average detection rates for different group sizes") %>%
kable_styling("striped") %>%
collapse_rows(1)| Group Size | Average Site Detection Probability |
|---|---|
| 1 | 0.297 |
| 2 | 0.379 |
| 3 | 0.439 |
| 4 | 0.485 |
| 5 | 0.523 |
We can generate detection function plots that show how detection rates fall-off over varying distances. Initial exploration showed that there may be some relationship between herbaceous understorey and detection probability for Sambar Deer. However, when pooled together we see that herbaceous understorey appears to have no effect on detection.
Show code
site_vars <- dplyr::tbl(con, dbplyr::in_schema("deervic", "curated_site_data")) %>%
dplyr::filter(SiteID %in% !!c(cams_curated$SiteID, "47191")) %>%
dplyr::collect() %>%
dplyr::mutate(HerbaceousUnderstoryCover = NNWHUCover + ENWHUCover,
SiteID = dplyr::case_when(SiteID == "47191" & CameraID == "HO04101053" ~ "47191A",
SiteID == "47191" & CameraID != "HO04101053" ~ "47191B",
TRUE ~ SiteID)) %>% # native + exotic herbaceous cover
dplyr::arrange(SiteID) %>%
as.data.frame()
n_distance_bins <- 5
delta <- 2.5
midpts <- c(1.25, 3.75, 6.25, 8.75, 11.25)
max_distance <- 12.5
gs <- 1:5
pa <- matrix(data = NA, nrow = max(gs), ncol = length(midpts))
for(j in 1:max(gs)) {
pa[j, ] <- sapply(midpts, function(x) {pr_closest(x-(0.5*delta),
x+(0.5*delta),
max_distance = max_distance,
n = j)})
}
det_curve <- model_fits[[top]]$draws("DetCurve", format = "draws_matrix") %>%
as.data.frame() %>%
# head(250) %>%
pivot_longer(cols = everything())
det_curve_wr <- det_curve %>%
mutate(var = stringr::str_extract(name, "(?<=\\[).+?(?=\\])")) %>%
separate(var, into = c("s", "Distance"), sep = ",") %>%
mutate(Distance = as.numeric(Distance)-0.5)
det_vars_pred <- site_vars %>%
mutate(s = as.character(1:nrow(.)),
herbaceouslvl = cut(HerbaceousUnderstoryCover,
breaks = c(0, 25, 50, 75, 105),
labels = c("0 - 25%",
"25 - 50%",
"50 - 75%",
"75 - 100%"),
include.lowest = T, right = FALSE))
det_curve_sum <- det_curve_wr %>%
mutate(Distance = as.numeric(Distance)) %>%
left_join(det_vars_pred) %>%
group_by(herbaceouslvl, Distance) %>%
summarise(median = quantile(value, 0.5),
q5 = quantile(value, 0.05),
q95 = quantile(value, 0.95))
scaleHN <- det_curve_sum %>% filter(Distance == 1.5) %>% pull(median) %>% mean() # Approx scaling
y_comb_list <- list()
for(i in 1:max(gs)) {
y_comb_list[[i]] <- colSums(multispecies_data$y[,,i]) %>%
as.data.frame() %>%
`colnames<-`("Count") %>%
mutate(Distance = midpts,
Bin = 1:5,
gs = i,
Prop = Count/(max(Count)),
CountS = Count/pa[Bin,i])
}
y_combined <- bind_rows(y_comb_list) %>%
group_by(Distance) %>%
mutate(SumCount = sum(Count)) %>%
ungroup() %>%
mutate(propC = Count/SumCount,
PropSM = (CountS/max(CountS)) * propC) %>%
group_by(Distance) %>%
summarise(PropS = sum(PropSM) * scaleHN)
detection_plot_HN <- ggplot(aes(x = Distance), data = det_curve_sum) +
geom_col(aes(y = PropS), fill = "grey70", colour = "grey20", width = 2.5, data = y_combined, alpha = 0.7) +
# geom_error bar(aes(ymin = PropSq5, ymax = PropSq95), data = y_combined) +
# geom_line(aes(y = HNS)) +
geom_ribbon(aes(ymin = q5, ymax = q95, fill = herbaceouslvl), alpha = 0.5) +
geom_line(aes(y = median, colour = herbaceouslvl)) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), limits = c(0,1)) +
scale_fill_manual(values = unname(delwp_cols[c(1,2,3,8)])) +
scale_colour_manual(values = unname(delwp_cols[c(1,2,3,8)])) +
ylab("Detection probability (p)") +
labs(fill = "Herbaceous\nunderstorey\ncover", colour = "Herbaceous\nunderstorey\ncover") +
theme_classic()
ggsave(plot = detection_plot_HN, filename = "outputs/plots/hazard_detection_plot.png", width = 2400, height = 1600, units = "px")
detection_plot_HN
Figure 10: Distance-sampling detection process for Deer (group size = 1-5). Herbaceous understorey cover in the area surrounding appears to have some impact on detection probability
Transect-surveys
For Sambar, Fallow and Red Deer we were able to estimate detection rates of the various methods used to detect deer (camera, footprints, pellets, rubbings and wallows).
Show code
all_draws <- model_fits[[top]]$draws("beta_trans_det", format = "df")
# det_marginal_effects <- list()
# det_plot <- list()
obs_vars <- unlist(stringr::str_remove_all(string = labels(terms(transect_formula)), pattern = "log|scale|\\)|\\("))
obs_cols <- unlist(stringr::str_remove_all(string = colnames(multispecies_data$trans_pred_matrix), pattern = "log|scale|\\)|\\("))
obs_logs <- unlist(stringr::str_detect(string = colnames(multispecies_data$trans_pred_matrix), pattern = "log\\("))
obs_labs <- c("Survey Method")
fac <- c(TRUE, TRUE, TRUE, TRUE, TRUE)
fac2 <- c(TRUE)
det_obs <- multispecies_data$transects %>%
mutate(Survey = factor(Survey))
det_marginal_effects <- list()
det_plot <- list()
params_w_levs <- levels(det_obs$Survey)
for(i in 1:(length(obs_cols))) {
det_marginal_effects[[i]] <- marginal_effects_cmd(all_draws,
param = "beta_trans_det",
param_number = i, log = obs_logs[i],
model_data = det_obs,
model_column = obs_vars[c(1,attr(multispecies_data$trans_pred_matrix, "assign")[-1])[i]],
transition = FALSE) %>%
mutate(group = param,
param = params_w_levs[i],
factor = fac[i],
variable = as.numeric(variable))
if(fac[i]) {
det_marginal_effects[[i]] <- det_marginal_effects[[i]]
}
}
marginal_prob <- function(x, pwr = 3) {
xm <- 1-x
return(1-(xm^pwr))
}
det_marginal_effects_bind <- bind_rows(det_marginal_effects) %>%
rowwise() %>%
mutate(value = case_when(param != "Camera" ~ marginal_prob(value),
TRUE ~ value))
det_marginal_effects_split <- split(det_marginal_effects_bind, det_marginal_effects_bind$group)
for(i in 1:length(det_marginal_effects_split)) {
if(fac2[i]) {
plot_data <- det_marginal_effects_split[[i]] %>%
mutate(variable = param,
param = group)
} else {
plot_data <- det_marginal_effects_split[[i]]
}
det_plot[[i]] <- plot_data
}
data_for_plot <- bind_rows(det_plot) %>%
mutate(species = "All (combined)") %>%
mutate(Density = "1 deer/km^2^")
data_d3 <- data_for_plot %>%
mutate(value = 1 - (1 - value)^3,
Density = "3 deer/km^2^")
# d3 averages
data_d3_av <- data_d3 %>%
group_by(Method = variable) %>%
summarise(mean = round(mean(value), 2),
q5 = round(quantile(value, 0.05), 2),
q95 = round(quantile(value, 0.95), 2)) %>%
ungroup() %>%
mutate(label = paste0("*p* = ", mean, " [90% CI: ", q5, ", ", q95, "]")) %>%
arrange(desc(mean))
data_d5 <- data_for_plot %>%
mutate(value = 1 - (1 - value)^5,
Density = "5 deer/km^2^")
data_combined <- bind_rows(data_for_plot, data_d3, data_d5)
transect_det_plot <- marginal_effects_plot_cmd_all(data_combined,
factor = TRUE,
ylab = "Detection probability") +
xlab("Survey method") +
scale_y_continuous(limits = c(0,1)) +
theme(legend.position = "top") +
scale_fill_manual(values = unname(delwp_cols[c(1,2,3)])) +
# scale_x_discrete(expand = c(0, 0)) +
theme_classic() +
theme(legend.text = ggtext::element_markdown())
ggsave(plot = transect_det_plot, filename = "outputs/plots/transect_detection_plot.png", width = 2400, height = 1400, units = "px")
transect_det_plot
Figure 11: Unconditional detection probabilities for the various methods of survey, across various deer densities. Camera trap detection probability is based on average deployment length (53 days), while transects are based on 3 transect searches
From these plots we can compare the detection rates between methods. For instance, when deer are present at a density of 3 per km2, than the highest probability of detection is from a Camera (p = 0.73 [90% CI: 0.65, 0.81]), followed by Pellet (p = 0.71 [90% CI: 0.63, 0.78]), Footprint (p = 0.64 [90% CI: 0.56, 0.72]), Rubbing (p = 0.36 [90% CI: 0.29, 0.44]), and lastly Wallow (p = 0.03 [90% CI: 0.01, 0.05]).
Activity
The proportion of time that deer were active and available for detection was 0.52 (90% CI: 0.48, 0.56). This proportion can be visualised in the following plot:
Show code
activity_modeled <-model_fits[[top]]$summary("activ")$mean
#### Availability ####
#### Determine availability prior ####
trigger.events <- camtrap_records_deer %>%
# filter(is.na(Behaviour)) %>% ### Remove biased behaviour
dplyr::mutate(Sample.Label = as.character(SiteID))
trigger.events$rtime <- activity::gettime(x = trigger.events$Time,
tryFormats = "%Y-%m-%d %H:%M:%OS")
act_result <- activity::fitact(trigger.events$rtime,
sample="model",
reps=50,
bw = 30,
show = FALSE)
avail <- list(creation=data.frame(rate = act_result@act[1],
SE = act_result@act[2]))
act_data <- data.frame(x = 12*act_result@data/pi)
act_model <- as.data.frame(12*act_result@pdf/pi)
ggplot() +
geom_histogram(data = act_data, aes(x = x, after_stat(ncount)),
breaks = seq(0, 24), width = 2, colour = "grey",
fill = delwp_cols[["Navy"]], alpha = 0.75) +
geom_ribbon(data = act_model, aes(ymin = lcl, ymax = ucl, x = x),
fill = "grey40", alpha = 0.5) +
geom_line(data = act_model, aes(y = y, x = x), colour = delwp_cols[["Navy"]], linewidth = 1) +
# coord_polar(start = 0) +
scale_fill_brewer() +
ylab("Proportional Activity") +
xlab("Time of Day") +
scale_x_continuous("", limits = c(0, 24),
breaks = 2*seq(0, 12),
labels = 2*seq(0, 12)) +
delwp_theme()
Figure 12: Activity times for all deer shows variable activity time throughout the day
Seasonal deployments
Our surveys took place in Spring, Summer and Autumn. Below we show the camera deployment dates for the 317 sites.
Show code
seasons <- data.frame(Date=seq.Date(from = as.Date("2021-09-01"),
to = as.Date("2023-05-01"), by = "month"),
Season = c('Spring', 'Spring', 'Spring',
'Summer', 'Summer', 'Summer',
'Autumn', 'Autumn', 'Autumn',
'Winter', 'Winter', 'Winter',
'Spring', 'Spring', 'Spring',
'Summer', 'Summer', 'Summer',
'Autumn', 'Autumn', 'Autumn')) %>%
mutate(DateEnd = seq.Date(from = as.Date("2021-10-01"),
to = as.Date("2023-06-01"), by = "month"),
Season = factor(Season, levels = c("Summer", "Autumn", "Winter", "Spring"))) %>%
cross_join(data.frame(SiteID = cams_curated$SiteID)) %>%
group_by(Date, DateEnd, Season) %>%
summarise(SiteIDmin = min(SiteID),
SiteIDmax = max(SiteID)) %>%
ungroup()
cam_op_plot <- cams_curated %>%
# mutate(DateDeploy = as.Date(format(DateDeploy, "%m-%d"),"%m-%d"),
# DateRetrieve = as.Date(format(DateRetrieve, "%m-%d"), "%m-%d")) %>%
ggplot() +
geom_rect(data = seasons, aes(xmin=Date, xmax=DateEnd,
ymin=SiteIDmin, ymax=SiteIDmax,
fill=Season
), alpha = 0.5) +
geom_segment(aes(x = DateDeploy, xend = DateRetrieve, y = SiteID, yend = SiteID), linewidth = 0.8, colour = "grey30") +
scale_x_date(date_breaks = "6 months", date_labels="%b-%y", name = "Deployment date") +
scale_fill_manual(values = c("#AF272F", "#E57200", "#0072CE", "#CEDC00")) +
delwp_theme() +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank())
ggsave(plot = cam_op_plot, filename = "outputs/plots/cam_op_plot.png", height = 2400, width = 2000, units = "px")
cam_op_plot
Figure 13: Deployment times of the cameras across the survey period
Abundance Processes
Bioregion
We used a spatially-derived random-effect of the bioregion of the sampled site. This random-effect allows us the ability to make predictions that include the variance associated with this random effect. This random-effect also minimises predictions in areas without deer detection’s (e.g. the mallee).
Show code
bioregion_sd <- model_fits[[top]]$summary("bioregion_sd")
bayesplot::color_scheme_set(bayes_theme)
bioregion_contribution <- function(model, data, species, species_index) {
eps_bioregion <- model$draws("eps_bioregion", format = "matrix")
which_sp <- which(stringr::str_detect(string = colnames(eps_bioregion),
pattern = paste0("eps_bioregion\\[", i)))
eps_bioregion <- eps_bioregion[,which_sp]
bioregion_data <- data[["bioreg_sf"]]
colnames(eps_bioregion) <- bioregion_data$bioregion
plot <- mcmc_areas(eps_bioregion, prob_outer = 0.9) +
ggtitle(species)
return(plot)
}
bio_plot_data <- list()
for(i in 1:length(species_names)) {
bio_plot_data[[i]] <- bioregion_contribution(model = model_fits[[top]],
data = multispecies_data,
species = species_names[i], species_index = i)
}
cowplot::plot_grid(plotlist = bio_plot_data, ncol = 1)
Figure 14: Influence of the bioregion random effect on abundance (log-scale).
Bioregion usually had a substantial impact on abundance estimates across the state. The variance (standard deviation) associated with the bioregion random effect was large for Sambar (\(\sigma\) = 3.96, 90% CI: 2.8 - 5.38), but less for Fallow (\(\sigma\) = 1.67, 90% CI: 0.85 - 2.73). Red Deer also exhibited large variance in the bioregion random effect (\(\sigma\) = 3.7, 90% CI: 2.48 - 5.23), as well Hog (\(\sigma\) = 4.4, 90% CI: 3.11 - 5.93). The larger variances for Hog, Red and Sambar reflect their more refined range boundaries, wheras Fallow show less clear range boundaries and thus show less variance across the state. Below we plot the relative contribution to abundance from each bioregion on a map of Victoria:
Show code
# Bioregion map
bioregion_eff <- model_fits[[top]]$summary("eps_bioregion", posterior::default_summary_measures())
bioreg <- VicmapR::vicmap_query("open-data-platform:vbioreg100") %>%
VicmapR::collect() %>%
group_by(bioregion, bioregion_code) %>%
summarise(geometry = st_union(geometry)) %>%
ungroup() %>%
st_transform(3111)
bioregion_nhat <- model_fits[[top]]$summary("Nhat_bioreg", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), posterior::default_summary_measures())
bioregion_density <- list()
for(i in 1:length(species_names)) {
which_sp <- which(stringr::str_detect(string = bioregion_nhat$variable,
pattern = paste0("Nhat_bioreg\\[", i)))
bioregion_nhat_sp <- bioregion_nhat[which_sp,]
bioregion_nhat_sp$bioregion <- multispecies_data$bioreg_sf$bioregion
bioregion_nhat_sp$Area <- as.numeric(table(multispecies_data$pred_bioreg))
bioregion_nhat_sp <- bioregion_nhat_sp %>%
mutate(`Average Density` = trimmed_mean/Area,
`Density (5%)` = q5/Area,
`Density (95%)` = q95/Area,
`CV` = sd/trimmed_mean,
label = paste0(bioregion, " [D = "
,round(`Average Density`,2),
" (",round(`Density (5%)`,2)," - ",round(`Density (95%)`, 2),")]"))
bioregion_density[[i]] <- left_join(bioreg, bioregion_nhat_sp, by = "bioregion")
}
bioregion_effects <- list()
for(i in 1:length(species_names)) {
which_sp <- which(stringr::str_detect(string = bioregion_eff$variable,
pattern = paste0("eps_bioregion\\[", i)))
bioregion_eff_sp <- bioregion_eff[which_sp,]
intercept <- model_fits[[top]]$summary(paste0("beta_psi[",i,",1]"), mean = ~ mean(.))$mean
bioregion_eff_sp$bioregion <- multispecies_data$bioreg_sf$bioregion
bioregion_eff_sp$Area <- as.numeric(table(multispecies_data$pred_bioreg))
bioregion_eff_sp <- bioregion_eff_sp %>%
mutate(`Average Effect` = mean)
bioregion_effects[[i]] <- left_join(bioreg, bioregion_eff_sp, by = "bioregion")
}
bioregion_ggplot <- function(effect_df, title) {
ggplot(data = effect_df) +
geom_sf(data = state, linewidth = 0.5, fill = "grey90") +
geom_sf(aes(fill = `Average Effect`)) +
geom_label_repel(data = effect_df, aes(label = bioregion, fill = `Average Effect`, geometry = geometry),
size = 2.5, stat = "sf_coordinates", colour = "white") +
scale_fill_gradient2(low = "#B9C600",
mid = delwp_cols[["Teal"]],
high = delwp_cols[["Navy"]], guide = guide_coloursteps(title = "Average effect\n (log-scale)")) +
labs(x = "", y = "") +
ggspatial::annotation_scale() +
delwp_theme()
}
# cl_br <- st_intersection(crown_land, bioreg %>% filter(bioregion == "Bridgewater") %>% st_transform(3111))
cpb <- cowplot::plot_grid(bioregion_ggplot(bioregion_effects[[1]], "Sambar"),
bioregion_ggplot(bioregion_effects[[2]], "Fallow"),
bioregion_ggplot(bioregion_effects[[3]], "Red"),
bioregion_ggplot(bioregion_effects[[4]], "Hog"), labels = c("A) Sambar deer", "B) Fallow deer", "C) Red deer", "D) Hog deer"), ncol = 2)
ggsave(plot = cpb, filename = "outputs/plots/bioregion_effects_plot.png", width = 4800, height = 3200, units = "px", bg = "white")
cpb
Figure 15: Effect of Victorian Bioregions on abundance of deer
Fixed-parameter Effects
Five fixed-effects were included in our multispecies model, with each fixed-effect having a sperately estimated coefficient for each species. Below we shouw the code used to extract and summarise the beta_psi coefficients and plot their marginal contribution to abundance.
Show code
# beta_summaries
beta_summaries <- model_fits[[top]]$summary("beta_psi")
beta_table_summary <- beta_summaries %>%
mutate(species = factor(rep(species_names, times = 6),levels = species_names),
covariate = rep(c("(Intercept)",
"Bare soil (%)",
"Nitrogen (%)",
"Distance to pastural land (m)",
"Precipitation seasonality",
'Forest edge per km^2^ (m)'), each = 4)) %>%
arrange(species) %>%
select(species, covariate, mean, sd, q5, q95)
format_cov <- function(sp, b, data = beta_summaries, variable = "beta_psi") {
fd <- data %>%
filter(variable == !!paste0(variable, "[", sp, ",",b+1,"]"))
paste0("$\\beta$ = ",
round(fd$mean, 2),
" [90% CI: ",
round(fd$q5, 2),
" – ",
round(fd$q95, 2),
"]")
}
ab_joined_list <- list()
for(j in 1:length(species_names)) {
beta_draws <- model_fits[[top]]$draws("beta_psi", format = "df")
# which_sp <- which(stringr::str_detect(string = colnames(beta_draws),
# pattern = paste0("beta_psi\\[", i)))
#
# beta_draws <- beta_draws[,which_sp]
ab_marginal_effects <- list()
ab_plot <- list()
ab_to_use <- ab_formula_3
phi_vars <- unlist(stringr::str_remove_all(string = labels(terms(ab_to_use)), pattern = "log|scale|sqrt|\\)|\\("))
phi_logs <- unlist(stringr::str_detect(string = labels(terms(ab_to_use)), pattern = "log\\("))
phi_sqrts <- c(0.5,0.5,1,1,0.5)
phi_labs <- c("Bare soil (%)",
"Nitrogen (%)",
"Distance to pastural land (m)",
"Precipitation seasonality",
'Forest edge per km^2^ (m)')
fac <- c(FALSE, FALSE, FALSE, FALSE, FALSE)
for(i in 1:length(phi_vars)) {
ab_marginal_effects[[i]] <- marginal_effects_cmd(beta_draws,
param = "beta_psi", species_index = j,
param_number = i+1, log = phi_logs[i],
model_data = multispecies_data[["raw_data"]] %>% mutate(Nitrogen = Nitrogen*100),
abundance = TRUE,
pwr = phi_sqrts[i],
model_column = phi_vars[i],
transition = FALSE) %>%
mutate(species = species_names[j])
}
ab_joined_list[[j]] <- bind_rows(ab_marginal_effects)
}
all_me_data <- bind_rows(ab_joined_list) %>%
mutate(param = case_when(param == "BareSoil" ~ "Bare soil (%)",
param == "Nitrogen" ~ "Nitrogen (%)",
param == "PastureDistance" ~ "Distance to pastural land (m)",
param == "BIO15" ~ "Precipitation seasonality",
param == "ForestEdge" ~ 'Forest edge per km^2^ (m)'))
cov_plot <- marginal_effects_plot_cmd_all(all_me_data %>% mutate(species = factor(stringr::str_to_sentence(species), levels = stringr::str_to_sentence(species_names)),
param = factor(param, levels = phi_labs)),
col = "DarkGreen",
factor = FALSE,
ylab = "Contribution to abundance (log-scale)") +
ggplot2::facet_grid(species~param, scales = "free") +
delwp_theme() +
theme(strip.text = ggtext::element_markdown()) +
xlab("Covariate value")
ggsave(plot = cov_plot, filename = "outputs/plots/abundance_covariate_effects.png", width = 3000, height = 2000, units = "px")
cov_plot
Figure 16: Marginal response curves for the various fixed-effect parameters used in the model.
Of the five fixed-effect covariates used in the top model we found that all had tangible effects for one or more species. The direction and magnitude of these effects often varied between species. Bare soil had a slight positive effect on Sambar abundance (\(\beta\) = 0.37 [90% CI: 0.02 – 0.72]), negligible effect on Fallow (\(\beta\) = -0.1 [90% CI: -0.63 – 0.4]) and Red Deer (\(\beta\) = -0.73 [90% CI: -2.01 – 0.54]), and a slight negative effect on Hog Deer abundance (\(\beta\) = -0.72 [90% CI: -1.57 – 0.12]). Nitrogen generally had weaker, more variable effects. The effect of nitrogen was weakly positive for Sambar (\(\beta\) = 0.31 [90% CI: -0.08 – 0.68]), and Hog Deer (\(\beta\) = 0.53 [90% CI: -0.54 – 1.6]); but weakly negative for Fallow (\(\beta\) = 0.22 [90% CI: -0.36 – 0.72]), and Red Deer (\(\beta\) = -0.45 [90% CI: -1.48 – 0.63]). The effect of distance to pasture was consistently negative across all four species, which can also be interpreted as a prefence of deer to be ‘closer’ to pastural areas. However the effect of this covariate varied in magnitude for Sambar (\(\beta\) = -0.53 [90% CI: -0.83 – -0.22]), Fallow (\(\beta\) = -0.89 [90% CI: -1.35 – -0.43]), Red (\(\beta\) = -1.61 [90% CI: -2.71 – -0.48]), and Hog Deer (\(\beta\) = -0.3 [90% CI: -0.8 – 0.22]). Sambar Deer weakly favored areas with less seasonal variation in precipitation (\(\beta\) = -0.23 [90% CI: -0.57 – 0.12]), with Fallow (\(\beta\) = 0.39 [90% CI: 0.02 – 0.8]), Red (\(\beta\) = 1.05 [90% CI: 0.05 – 2.02]), and Hog Deer (\(\beta\) = 0.77 [90% CI: 0.21 – 1.31]) all showing positive relationships between precipitation seasonality and abundance. Lastly, we saw that the amount of forest edge in the landscape positively impacted abundance for Sambar (\(\beta\) = 0.17 [90% CI: -0.06 – 0.4]), and Fallow Deer (\(\beta\) = 0.61 [90% CI: 0.29 – 0.93]). This relationship was weakly negative for Red Deer (\(\beta\) = -0.7 [90% CI: -1.39 – 0.01]), and weakly positive for Hog Deer (\(\beta\) = 0.32 [90% CI: -0.07 – 0.69]). We can also show these ciefficients in a table:
Show code
which_sig <- which(beta_table_summary$q5 * beta_table_summary$q95 > 0)
write.csv(beta_table_summary, "outputs/tables/beta_coefs.csv")
beta_table_summary %>%
kbl("html", digits = 3, caption = "beta-coeffiient estimates for the drivers of species abundance, coefficient estimates not-overalpping zero are shown in bold") %>%
kable_styling("striped") %>%
collapse_rows(1) %>%
row_spec(row = which_sig, bold = TRUE)| species | covariate | mean | sd | q5 | q95 |
|---|---|---|---|---|---|
| Sambar Deer | (Intercept) | -3.615 | 0.956 | -5.182 | -2.166 |
| Bare soil (%) | 0.368 | 0.215 | 0.023 | 0.723 | |
| Nitrogen (%) | 0.313 | 0.230 | -0.079 | 0.683 | |
| Distance to pastural land (m) | -0.535 | 0.189 | -0.832 | -0.221 | |
| Precipitation seasonality | -0.230 | 0.208 | -0.573 | 0.122 | |
| Forest edge per km2 (m) | 0.166 | 0.137 | -0.057 | 0.396 | |
| Fallow Deer | (Intercept) | -1.680 | 0.513 | -2.563 | -0.930 |
| Bare soil (%) | -0.103 | 0.319 | -0.630 | 0.403 | |
| Nitrogen (%) | 0.220 | 0.335 | -0.359 | 0.723 | |
| Distance to pastural land (m) | -0.892 | 0.278 | -1.345 | -0.432 | |
| Precipitation seasonality | 0.392 | 0.239 | 0.023 | 0.799 | |
| Forest edge per km2 (m) | 0.607 | 0.192 | 0.288 | 0.931 | |
| Red Deer | (Intercept) | -6.781 | 1.170 | -8.861 | -4.997 |
| Bare soil (%) | -0.725 | 0.793 | -2.013 | 0.543 | |
| Nitrogen (%) | -0.447 | 0.635 | -1.481 | 0.628 | |
| Distance to pastural land (m) | -1.613 | 0.683 | -2.705 | -0.482 | |
| Precipitation seasonality | 1.049 | 0.599 | 0.052 | 2.021 | |
| Forest edge per km2 (m) | -0.699 | 0.426 | -1.386 | 0.005 | |
| Hog Deer | (Intercept) | -8.297 | 1.307 | -10.503 | -6.305 |
| Bare soil (%) | -0.717 | 0.512 | -1.567 | 0.120 | |
| Nitrogen (%) | 0.526 | 0.651 | -0.543 | 1.601 | |
| Distance to pastural land (m) | -0.298 | 0.312 | -0.799 | 0.222 | |
| Precipitation seasonality | 0.772 | 0.334 | 0.209 | 1.312 | |
| Forest edge per km2 (m) | 0.318 | 0.231 | -0.068 | 0.690 |
Dispersion
Below we show the parameter estimates for the dispersion parameters estimated within our negative binomial models for camera counts and Royle-Nichols model.
Show code
dispersion_data <- model_fits[[top]]$summary(c("od", "odRN")) %>%
mutate(species = c(species_names, "All"),
model = c(rep("Camera count", 4), "Royle-Nichols"),
parameter = c("$\\phi_1$", "$\\phi_2$", "$\\phi_3$", "$\\phi_4$", "$\\bar \\phi$")) %>%
select(model, species, parameter, mean, sd, q5, q95)
format_od <- function(data) {
paste0(dispersion_data$parameter, " = ", round(dispersion_data$mean, 3), " [90 % CI: ", round(dispersion_data$q5,3), " -- ", round(dispersion_data$q95,3), "]")
}
formatred_od <- format_od(dispersion_data)
dispersion_data %>%
kbl("html", digits = 3, caption = "Estimates of the dispersion parameter within negative binomial submodels") %>%
kable_styling("striped") %>%
collapse_rows(1)| model | species | parameter | mean | sd | q5 | q95 |
|---|---|---|---|---|---|---|
| Camera count | Sambar Deer | \(\phi_1\) | 0.114 | 0.014 | 0.093 | 0.137 |
| Fallow Deer | \(\phi_2\) | 0.021 | 0.004 | 0.015 | 0.030 | |
| Red Deer | \(\phi_3\) | 0.047 | 0.016 | 0.025 | 0.076 | |
| Hog Deer | \(\phi_4\) | 0.064 | 0.018 | 0.038 | 0.095 | |
| Royle-Nichols | All | \(\bar \phi\) | 0.790 | 0.140 | 0.591 | 1.041 |
Estimates of abundance were also dependent upon overdispersion (excess variation) in the camera counts and transect detections. Sambar (\(\phi_1\) = 0.114 [90 % CI: 0.093 – 0.137]), Fallow (\(\phi_2\) = 0.021 [90 % CI: 0.015 – 0.03]), Red (\(\phi_3\) = 0.047 [90 % CI: 0.025 – 0.076]), and Hog Deer (\(\phi_4\) = 0.064 [90 % CI: 0.038 – 0.095]) all had large amounts of overdispersion (where lower values (<1) equate to more variance). The Royle-Nichols model had less variance (than camera counts) associated with transect-level detections/non-detections (\(\bar \phi\) = 0.79 [90 % CI: 0.591 – 1.041]).
Relating abundance to vegetation
From our point density estimates of deer, we can model relationships between deer densities at sites and a range of structural vegetation measures collected during the surveys. For this part of the analysis we implement a multivariate regression.
Firstly the response variables we model are:
- Bare Ground Cover (Beta Distribution)
- Native Woody Understorey Cover (Beta Distribution)
- Native Non-woody Understorey Cover (Beta Distribution)
- Presence of seedlings (Negative Binomial Distribution)
- Presence of saplings (Negative Binomial Distribution)
- Presence of exotic flora (Bernoulli Distribution)
While our main predictor of interest is the density of deer at a site, we also control for predictors that are expected to have an impact on vegetation growth and composition. The predictors controlled for include features relating to climate, ecology, topography, fire history and forest structure. They are the following predictors:
- BIO01: Annual Mean Temperature
- BIO04: Temperature Seasonality
- BIO06: Min Temperature of Coldest Month
- BIO12: Annual Precipitation
- BIO15: Precipitation Seasonality (Coefficient of Variation)
- Soil Nitrogen
- Time since last fire (factor grouping)
- TWIND (Topographic wetness index)
- Forest Edge
- Canopy Cover
- Top Height
- Wallaby Presence
- EVC Group (random effect with 19 groups)
Data Preparation
We can use the sum of estimated deer density for the four species at each site as our key predictor variable. Given that we already extracted a range of bioclimatic and environmental covariates for the original model, we can combine those with the data available here.
Note that this model is run on a subset of sites (sites surveyed as part of the statewide deer survey), as the hog deer sites did not have complete vegetation metrics recorded.
Show code
density <- model_fits[[top]]$draws("N_site", format = "matrix")
site_estimates <- list()
site_density <- list()
for(i in 1:n_site) {
site_estimates[[i]] <- character()
for(j in 1:length(deer_species_all)) {
site_estimates[[i]][j] <- paste0("N_site[", j, ",", i, "]")
}
site_density[[i]] <- sum(colMeans(density)[site_estimates[[i]]])
}
site_density_ul <- unlist(site_density)
site_density_df <- cams_curated %>%
mutate(AllDeerDensity = site_density_ul,
Bioregion = factor(multispecies_data$site_bioreg),
EVC = factor(multispecies_data$site_evc)) %>%
filter(ProjectShortName == "StatewideDeer") %>%
left_join(multispecies_data[["raw_data"]] %>% select(-EVC, - BIOREGION)) %>%
mutate(TSLF_bin = as.factor(round(TSLF_bin)))
wallaby_pa <- tbl(con, dbplyr::in_schema("camtrap", "processed_site_substation_presence_absence")) %>%
filter(scientific_name == "Wallabia bicolor" & ProjectShortName %in% project_short_name) %>%
collect() %>%
transmute(SiteID,
Wallaby = as.logical(Presence))
site_vars <- dplyr::tbl(con, dbplyr::in_schema("deervic", "curated_site_data")) %>%
dplyr::filter(SiteID %in% !!c(cams_curated$SiteID, "47191")) %>%
dplyr::collect() %>%
dplyr::left_join(wallaby_pa) %>%
dplyr::mutate(HerbaceousUnderstoryCover = NNWHUCover + ENWHUCover,
Exotics = EWUCover + ENWHUCover,
ExoticPresence = as.integer(as.logical(Exotics)),
SeedlingsOrSaplings = as.integer(as.logical(Seedlings+Saplings)),
SeedlingsPresence = as.integer(as.logical(Seedlings)),
SaplingsPresence = as.integer(as.logical(Saplings)),
SiteID = dplyr::case_when(SiteID == "47191" & CameraID == "HO04101053" ~ "47191A",
SiteID == "47191" & CameraID != "HO04101053" ~ "47191B",
TRUE ~ SiteID)) %>% # native + exotic herbaceous cover
dplyr::arrange(SiteID) %>%
as.data.frame()
site_density_df_joined <- site_density_df %>%
left_join(site_vars) %>%
mutate(BGroundCover = ifelse(BGroundCover == 0, 0.125, BGroundCover)/100,
NWUCover = ifelse(NWUCover == 0, 0.125, NWUCover)/100,
NNWHUCover = ifelse(NNWHUCover == 0, 0.125, NNWHUCover)/100)Model Execution
We use the brms package to run our model in a Bayesian framework (Bürkner 2017, 2018, 2021). Our model is run in parallel on four cores over four chains for 1,000 sampling iterations per chain (1,000 warm up iterations).
Show code
bform1 <- bf(mvbind(BGroundCover,
NWUCover,
NNWHUCover,
Seedlings,
Saplings,
ExoticPresence) ~ AllDeerDensity
+ scale(BIO01)
+ scale(BIO04)
+ scale(BIO06)
+ scale(BIO12)
+ scale(BIO15)
+ scale(Nitrogen)
+ TSLF_bin
+ scale(TWIND)
+ scale(sqrt(ForestEdge))
+ scale(sqrt(CanopyCov))
+ scale(sqrt(TopHeight))
+ Wallaby
+ (1|EVC))
# Specifiy model families and links
vegfit <- brm(bform1,
data = site_density_df_joined,
family = list(Beta(link = "logit"),
Beta(link = "logit"),
Beta(link = "logit"),
negbinomial(),
negbinomial(),
bernoulli(link = "logit")),
inits = "0",
chains = 4,
cores = 4,
iter = 2000,
warmup = 1000,
backend = "cmdstanr")
saveRDS(vegfit, "outputs/models/vegfit.rds")Model Results
Model Summary
Below we show a complete summary of the model, including all coefficients estimated for all the parameters and response variables. Standard deviations for the EVC random effect and the estimates of the \(\phi\) parameter in the beta-distribution models, and the shape parameter in the negative binomial models are also shown. Note that all Rhats > 1.05 and all ESS > 200, suggesting good chain mixing and convergence.
Show code
vegfit <- readRDS("outputs/models/vegfit.rds")
veg_coefs <- fixef(vegfit, probs = c(0.05, 0.95), pars = c("BGroundCover_AllDeerDensity",
"NWUCover_AllDeerDensity",
"NNWHUCover_AllDeerDensity",
"Seedlings_AllDeerDensity",
"Saplings_AllDeerDensity",
"ExoticPresence_AllDeerDensity")) %>%
as.data.frame() %>%
mutate(Effect = paste0(c("Bare Ground Cover",
"Native Woody Understorey Cover",
"Native Non-woody Understorey Cover",
"Presence of seedlings",
"Presence of saplings",
"Presence of exotic flora"), " ~ Deer Density"))
all_veg_coef <- fixef(vegfit, probs = c(0.05, 0.95))
print_conditional <- function(x, i, d = 3, OR = FALSE) {
row <- x[i,]
name <- "\\beta"
if(OR) {
row <- row %>%
mutate(across(where(is.numeric), exp))
name <- "OR"
}
paste0("$",name,"$ = ", round(row[["Estimate"]], d), " [90% CI: ", round(row[["Q5"]], d), ", ", round(row[["Q95"]], d), "]")
}
veg_coefs_save <- veg_coefs %>%
transmute(effect = Effect, mean = Estimate, sd = Est.Error, `5 %` = Q5, `95 %` = Q95) %>%
mutate(across(where(is.numeric), round, 3))
write.csv(veg_coefs_save, "outputs/tables/veg_coefs.csv", row.names = F)
options(width=600)
summary(vegfit) Family: MV(beta, beta, beta, negbinomial, negbinomial, bernoulli)
Links: mu = logit; phi = identity
mu = logit; phi = identity
mu = logit; phi = identity
mu = log; shape = identity
mu = log; shape = identity
mu = logit
Formula: BGroundCover ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + Wallaby + (1 | EVC)
NWUCover ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + Wallaby + (1 | EVC)
NNWHUCover ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + Wallaby + (1 | EVC)
Seedlings ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + Wallaby + (1 | EVC)
Saplings ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + Wallaby + (1 | EVC)
ExoticPresence ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + Wallaby + (1 | EVC)
Data: site_density_df_joined (Number of observations: 250)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~EVC (Number of levels: 19)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(BGroundCover_Intercept) 0.26 0.12 0.06 0.51 1.00 1351 1309
sd(NWUCover_Intercept) 0.13 0.10 0.01 0.36 1.00 1422 1735
sd(NNWHUCover_Intercept) 0.23 0.13 0.02 0.51 1.00 1039 1630
sd(Seedlings_Intercept) 1.16 0.39 0.52 2.06 1.00 1138 1804
sd(Saplings_Intercept) 0.83 0.33 0.23 1.52 1.00 974 829
sd(ExoticPresence_Intercept) 0.63 0.41 0.04 1.58 1.01 757 1330
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
BGroundCover_Intercept -1.93 0.23 -2.39 -1.49 1.00 2558 2063
NWUCover_Intercept -0.52 0.20 -0.90 -0.14 1.00 3158 3469
NNWHUCover_Intercept -0.54 0.22 -0.96 -0.10 1.00 3306 3105
Seedlings_Intercept 3.13 0.61 1.93 4.30 1.00 1950 2638
Saplings_Intercept 3.45 0.52 2.45 4.47 1.00 1870 2448
ExoticPresence_Intercept -1.93 0.61 -3.15 -0.78 1.00 3061 3025
BGroundCover_AllDeerDensity -0.04 0.03 -0.09 0.02 1.00 4508 3244
BGroundCover_scaleBIO01 0.42 0.40 -0.35 1.20 1.00 1794 2119
BGroundCover_scaleBIO04 0.01 0.26 -0.48 0.50 1.00 1809 2187
BGroundCover_scaleBIO06 -0.29 0.46 -1.18 0.59 1.00 1807 2079
BGroundCover_scaleBIO12 -0.29 0.15 -0.57 0.00 1.00 4731 3442
BGroundCover_scaleBIO15 -0.20 0.08 -0.36 -0.04 1.00 3442 3334
BGroundCover_scaleNitrogen -0.07 0.16 -0.40 0.24 1.00 4304 3195
BGroundCover_TSLF_bin3 0.07 0.22 -0.35 0.50 1.00 2371 2631
BGroundCover_TSLF_bin4 0.05 0.25 -0.43 0.55 1.00 2395 2597
BGroundCover_TSLF_bin5 -0.24 0.25 -0.72 0.25 1.00 2512 2389
BGroundCover_scaleTWIND -0.02 0.11 -0.23 0.21 1.00 5371 3386
BGroundCover_scalesqrtForestEdge -0.03 0.07 -0.18 0.11 1.00 5161 2869
BGroundCover_scalesqrtCanopyCov 0.00 0.07 -0.14 0.14 1.00 4824 3191
BGroundCover_scalesqrtTopHeight -0.24 0.08 -0.39 -0.09 1.00 5686 3329
BGroundCover_WallabyTRUE -0.28 0.13 -0.53 -0.03 1.00 5501 2942
NWUCover_AllDeerDensity -0.05 0.03 -0.10 0.00 1.00 4897 3169
NWUCover_scaleBIO01 0.11 0.43 -0.76 0.94 1.00 1815 2598
NWUCover_scaleBIO04 -0.06 0.27 -0.57 0.48 1.00 1867 2634
NWUCover_scaleBIO06 -0.11 0.50 -1.07 0.90 1.00 1845 2493
NWUCover_scaleBIO12 0.20 0.15 -0.10 0.49 1.00 3308 3480
NWUCover_scaleBIO15 0.20 0.08 0.03 0.36 1.00 3546 3404
NWUCover_scaleNitrogen 0.01 0.17 -0.32 0.34 1.00 3779 3252
NWUCover_TSLF_bin3 -0.35 0.20 -0.73 0.06 1.00 2942 3117
NWUCover_TSLF_bin4 -0.57 0.23 -1.02 -0.13 1.00 2645 3202
NWUCover_TSLF_bin5 -0.89 0.25 -1.37 -0.42 1.00 2692 3335
NWUCover_scaleTWIND -0.00 0.13 -0.25 0.25 1.00 4692 3080
NWUCover_scalesqrtForestEdge 0.05 0.08 -0.10 0.20 1.00 5489 3338
NWUCover_scalesqrtCanopyCov -0.05 0.07 -0.19 0.09 1.00 4882 3279
NWUCover_scalesqrtTopHeight -0.15 0.08 -0.30 0.01 1.00 5121 3044
NWUCover_WallabyTRUE -0.02 0.14 -0.28 0.25 1.00 5602 3357
NNWHUCover_AllDeerDensity 0.06 0.03 0.01 0.12 1.00 5173 3576
NNWHUCover_scaleBIO01 0.16 0.50 -0.79 1.11 1.00 1577 2415
NNWHUCover_scaleBIO04 -0.29 0.31 -0.88 0.30 1.00 1654 2287
NNWHUCover_scaleBIO06 -0.48 0.57 -1.55 0.61 1.00 1603 2152
NNWHUCover_scaleBIO12 0.40 0.16 0.09 0.71 1.00 4099 3204
NNWHUCover_scaleBIO15 0.10 0.09 -0.08 0.29 1.00 3835 3330
NNWHUCover_scaleNitrogen -0.28 0.18 -0.63 0.07 1.00 3966 3262
NNWHUCover_TSLF_bin3 -0.34 0.22 -0.78 0.10 1.00 2738 3115
NNWHUCover_TSLF_bin4 -0.17 0.24 -0.66 0.31 1.00 2763 3086
NNWHUCover_TSLF_bin5 -0.19 0.25 -0.69 0.30 1.00 2745 3143
NNWHUCover_scaleTWIND -0.08 0.13 -0.35 0.19 1.00 3984 3175
NNWHUCover_scalesqrtForestEdge -0.00 0.08 -0.17 0.16 1.00 5110 3185
NNWHUCover_scalesqrtCanopyCov 0.03 0.08 -0.12 0.18 1.00 4914 2904
NNWHUCover_scalesqrtTopHeight 0.02 0.09 -0.16 0.18 1.00 4177 3033
NNWHUCover_WallabyTRUE 0.15 0.14 -0.14 0.43 1.00 5899 2878
Seedlings_AllDeerDensity 0.05 0.07 -0.09 0.21 1.00 4741 3409
Seedlings_scaleBIO01 0.14 0.91 -1.58 1.95 1.00 1849 2668
Seedlings_scaleBIO04 -0.10 0.59 -1.25 1.07 1.00 1886 2770
Seedlings_scaleBIO06 -0.57 1.00 -2.56 1.37 1.00 1817 2614
Seedlings_scaleBIO12 0.74 0.43 -0.10 1.58 1.00 1931 2620
Seedlings_scaleBIO15 -0.23 0.28 -0.76 0.32 1.00 2364 2890
Seedlings_scaleNitrogen -0.91 0.41 -1.69 -0.12 1.00 3201 2476
Seedlings_TSLF_bin3 -1.15 0.53 -2.17 -0.11 1.00 2205 2761
Seedlings_TSLF_bin4 -1.21 0.54 -2.26 -0.14 1.00 2233 2994
Seedlings_TSLF_bin5 -1.64 0.61 -2.80 -0.48 1.00 2022 2834
Seedlings_scaleTWIND -0.03 0.25 -0.53 0.47 1.00 3402 3122
Seedlings_scalesqrtForestEdge -0.39 0.19 -0.75 -0.02 1.00 4342 3329
Seedlings_scalesqrtCanopyCov -0.10 0.18 -0.46 0.27 1.00 4931 3322
Seedlings_scalesqrtTopHeight -0.94 0.21 -1.36 -0.54 1.00 4254 3074
Seedlings_WallabyTRUE 0.20 0.37 -0.55 0.90 1.00 4623 3433
Saplings_AllDeerDensity -0.01 0.06 -0.12 0.11 1.00 4054 3228
Saplings_scaleBIO01 0.16 0.86 -1.54 1.82 1.00 1383 2339
Saplings_scaleBIO04 0.15 0.53 -0.87 1.21 1.00 1432 2168
Saplings_scaleBIO06 0.27 0.95 -1.55 2.17 1.00 1351 2305
Saplings_scaleBIO12 0.23 0.32 -0.40 0.86 1.00 2140 3204
Saplings_scaleBIO15 -0.42 0.19 -0.78 -0.05 1.00 2955 2892
Saplings_scaleNitrogen 0.19 0.34 -0.47 0.85 1.00 3760 3011
Saplings_TSLF_bin3 -0.65 0.47 -1.57 0.30 1.00 2153 2776
Saplings_TSLF_bin4 -1.86 0.50 -2.85 -0.88 1.00 2264 2653
Saplings_TSLF_bin5 -1.47 0.53 -2.50 -0.38 1.00 1988 2685
Saplings_scaleTWIND 0.12 0.25 -0.38 0.61 1.00 2973 2628
Saplings_scalesqrtForestEdge -0.36 0.16 -0.68 -0.04 1.00 4175 3295
Saplings_scalesqrtCanopyCov 0.04 0.14 -0.23 0.31 1.00 3750 3096
Saplings_scalesqrtTopHeight -0.56 0.17 -0.90 -0.23 1.00 2626 2847
Saplings_WallabyTRUE 0.22 0.30 -0.35 0.80 1.00 4187 3579
ExoticPresence_AllDeerDensity 0.15 0.07 0.02 0.30 1.00 3798 2622
ExoticPresence_scaleBIO01 0.12 1.13 -2.12 2.42 1.00 1627 2337
ExoticPresence_scaleBIO04 0.07 0.71 -1.35 1.46 1.00 1673 2661
ExoticPresence_scaleBIO06 -0.37 1.31 -3.00 2.24 1.00 1573 2156
ExoticPresence_scaleBIO12 0.33 0.40 -0.46 1.12 1.00 3842 2956
ExoticPresence_scaleBIO15 0.32 0.21 -0.09 0.72 1.00 3219 3083
ExoticPresence_scaleNitrogen -1.21 0.47 -2.18 -0.32 1.00 3210 2717
ExoticPresence_TSLF_bin3 0.51 0.61 -0.71 1.67 1.00 1974 3007
ExoticPresence_TSLF_bin4 0.57 0.67 -0.75 1.85 1.00 1848 2894
ExoticPresence_TSLF_bin5 0.82 0.68 -0.54 2.13 1.00 1916 2616
ExoticPresence_scaleTWIND -0.39 0.32 -1.02 0.22 1.00 4482 3319
ExoticPresence_scalesqrtForestEdge 0.11 0.19 -0.25 0.50 1.00 4407 3114
ExoticPresence_scalesqrtCanopyCov 0.23 0.19 -0.15 0.61 1.00 5851 3125
ExoticPresence_scalesqrtTopHeight 0.12 0.21 -0.27 0.52 1.00 4431 2990
ExoticPresence_WallabyTRUE 0.46 0.35 -0.23 1.14 1.00 4962 3125
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi_BGroundCover 6.95 0.72 5.63 8.47 1.00 4130 2870
phi_NWUCover 4.03 0.35 3.41 4.72 1.00 5426 3127
phi_NNWHUCover 3.05 0.26 2.57 3.57 1.00 6043 2941
shape_Seedlings 0.29 0.03 0.23 0.36 1.00 3399 2982
shape_Saplings 0.37 0.04 0.30 0.46 1.00 4142 2742
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Deer Density Conditional Effects
Given our main interest in the model is deer density and how it impacts the suite of response variables, we generate marginal/conditional effects plots for the relationship between deer and the vegetation metric. Based on these plots we can see a “significant” (confidence intervals not-overlapping zero) effect of deer density on: bare ground (-ve), native woody understorey (-ve), native non-woody/herbaceous understorey (+ve), and presence of exotic species (+ve):
Show code
plot_list <- conditional_effects(vegfit, effects = "AllDeerDensity", prob = 0.9)
plot_list_50 <- conditional_effects(vegfit, effects = "AllDeerDensity", prob = 0.5)
plot_deer_relationship <- function(plot_data, plot_data_50, yaxis) {
ggplot(plot_data) +
geom_ribbon(aes(x = AllDeerDensity, ymin = `lower__`, ymax = `upper__`),
fill = delwp_cols[["Navy"]], alpha = 0.3) +
geom_ribbon(aes(x = AllDeerDensity, ymin = `lower__`, ymax = `upper__`),
fill = delwp_cols[["Navy"]], alpha = 0.3, data = plot_data_50) +
geom_smooth(aes(x = AllDeerDensity, y = `estimate__`),
colour = delwp_cols[["Navy"]], method = "loess", formula = y~x, linewidth = 2) +
ylab(paste(yaxis)) +
xlab("Deer density (per km^2^)") +
delwp_theme() +
theme(axis.title.x = ggtext::element_markdown())
}
vegplot <- cowplot::plot_grid(plot_deer_relationship(plot_list[[1]], plot_list_50[[1]], "Bare ground cover"),
plot_deer_relationship(plot_list[[2]], plot_list_50[[2]], "Native woody understorey cover"),
plot_deer_relationship(plot_list[[3]], plot_list_50[[3]], "Native herbaceous understorey cover"),
plot_deer_relationship(plot_list[[4]], plot_list_50[[4]], "Seedling counts") + scale_y_sqrt(),
plot_deer_relationship(plot_list[[5]], plot_list_50[[5]], "Sapling counts"),
plot_deer_relationship(plot_list[[6]], plot_list_50[[6]], "Presence of exotic flora"),
ncol = 3, labels = "AUTO")
ggsave(plot = vegplot, filename = "outputs/plots/vegetation_effects_plot.png", width = 3400, height = 2200, units = "px")
vegplot
Figure 17: Impact of deer density on a range of vegetation measure, including (A) bare ground cover, (B) native woody understorey, (C) native herbaceous understorey cover, (D) seedlings, (E) saplings, and (F) th presence of exotic species. The line shows the mean estimate, with the shading representing 50% and 90% confidence intervals. Inner and outer shading represents 50% and 90% confidence levels respectively.
We found deer density to have several plausible relationships with several structural vegetation componentsw, including negative relationships with bare ground cover (\(\beta\) = -0.035 [90% CI: -0.083, 0.009]), native woody understorey cover (\(\beta\) = -0.046 [90% CI: -0.089, -0.003]), and a positive relationship with native non-woody/herbaceous understorey cover (\(\beta\) = 0.064 [90% CI: 0.017, 0.111]). We also found strong evidence suggesting the presence of weeds is more likely at sites with higher densities of deer (\(\beta\) = 0.149 [90% CI: 0.035, 0.276], \(OR\) = 1.161 [90% CI: 1.036, 1.318]). We did not find tangible relationships between deer abundance and seedling or sapling counts.
Analysis Summary
Our studies detected deer at 148/317 cameras across Victoria and some form of deer sign (camera or transect) at 186/317 sites. Our average estimates of deer density at the 317 sites ranged from 0 to 28.3 deer per km2. Across approximately 74,570 km2 of public land in Victoria we estimate total deer abundance of the four species investigated in this study (Sambar, Fallow, Red and Hog) at 191,153 (90 % CI: 146,732 - 255,490. We also find plausible links between deer abundance and increased presence of exotic weed species.
Comparisons to other studies
Below we show a table of studies that have previously estimated deer density in Australia.
Show code
abundance_studies <- readr::read_csv("data/Australian_Deer_Studies_Density_Estimates.csv") %>%
mutate(`Density LB` = as.numeric(stringr::str_trim(`Density Min`)),
`Density UB` = as.numeric(stringr::str_trim(`Density Max`)))
abundance_table_present <- abundance_studies %>%
transmute(State, Species, Title, Author, `Study Sites`, `Year Start`,
`Year End`, `Density Mean`, `Density LB`,
`Density UB`) %>%
arrange(State, Species, Title, `Study Sites`)
write.csv(abundance_table_present, "outputs/tables/other_studies.csv", row.names = FALSE)
abundance_table_present %>%
kbl("html", digits = 3, caption = "Other studies estimating deer density in Australia") %>%
kable_styling("striped", font_size = 8) %>%
collapse_rows(1:3) %>%
scroll_box(height = "800px", width = "800px")| State | Species | Title | Author | Study Sites | Year Start | Year End | Density Mean | Density LB | Density UB |
|---|---|---|---|---|---|---|---|---|---|
| NSW | Fallow | Estimating deer density and abundance using spatial mark–resight models with camera trap data | Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Montane Woodland/Fores in Kosciuszko NP | 2017 | 2019 | 0.29 | 0.12 | 0.53 |
| Sambar | Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Montane Forest in Kosciuszko NP | 2017 | 2019 | 0.73 | 0.33 | 1.35 | ||
| Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Montane Woodland/Fores in Kosciuszko NP | 2017 | 2019 | 2.49 | 1.67 | 3.50 | |||
| Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Wetland Reserve | 2017 | 2019 | 0.48 | 0.24 | 0.84 | |||
| QLD | Red | I just want to count them! Considerations when choosing a deer population monitoring method | Amos, M., G. Baxter, N. Finch, A. Lisle, and P. Murray | Cressbrook Dam catchment | 2010 | 2012 | 45.75 | 39.80 | 51.70 |
| Amos, M., G. Baxter, N. Finch, A. Lisle, and P. Murray | Cressbrook Dam catchment | 2010 | 2012 | 26.50 | 23.70 | 29.30 | |||
| Amos, M., G. Baxter, N. Finch, A. Lisle, and P. Murray | Cressbrook Dam catchment | 2011 | 2011 | 26.26 | 12.92 | 53.30 | |||
| Rusa | Estimating deer density and abundance using spatial mark–resight models with camera trap data | Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Central Coastal QLD Woodland | 2017 | 2019 | 10.34 | 7.84 | 13.32 | |
| Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Central Coastal QLD Woodland | 2017 | 2019 | 0.68 | 0.21 | 1.77 | |||
| Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Samsonvale Resivoir | 2017 | 2019 | 3.11 | 1.76 | 5.07 | |||
| TAS | Fallow | Report of state-wide census of wild fallow deer in Tasmania project: Part A: Baseline aerial survey of fallow deer population, central and north-eastern Tasmania. | Lethbridge, M. R., Stead, M. G., & Wells, C. | Northern Midlands | 2019 | 2019 | 2.70 | 1.67 | 3.71 |
| VIC | Estimating deer density and abundance using spatial mark–resight models with camera trap data | Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Cardinia Resivoir Woodland | 2017 | 2019 | 2.09 | 1.46 | 2.92 | |
| Hog | The influence of evolutionary history and body size on partitioning of habitat resources by mammalian herbivores in south-eastern Australia | Davis, N. E., Gordon, I. R., & Coulson, G | Yanicke | 2003 | 2004 | 3.01 | 1.73 | 4.29 | |
| Red | Estimating deer density and abundance using spatial mark–resight models with camera trap data | Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Sugarloaf Resivoir Woodland | 2017 | 2019 | 24.57 | 19.79 | 30.64 | |
| Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Yan Yean Woodland | 2017 | 2019 | 19.76 | 17.58 | 22.17 | |||
| Sambar | BBRR Feral Pig Survey Project Delicknora 2023 | Durand M., and Tame, A. | Delicknora | NA | NA | 1.70 | 1.70 | 1.70 | |
| BBRR Feral Pig Survey Project Wulgulmerang 2022 | Durand M., and Tame, A. | Wulgulmerang | 2022 | 2022 | 3.29 | 3.29 | 3.29 | ||
| BBRR Feral Pig Survey Project Wulgulmerang 2023 | Durand M., and Tame, A. | Wulgulmerang | NA | NA | 6.40 | 6.40 | 6.40 | ||
| Estimating deer density and abundance using spatial mark–resight models with camera trap data | Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Cardinia Resivoir Woodland | 2017 | 2019 | 11.94 | 8.44 | 16.48 | ||
| Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T | Yan Yean Woodland | 2017 | 2019 | 3.93 | 2.53 | 6.29 | |||
| The application of catch–effort models to estimate the efficacy of aerial shooting operations on sambar deer | Ramsey David S. L., McMaster Damien, Thomas Elaine | Alpine NP - Bogong | 2020 | 2022 | 0.89 | 0.70 | 1.10 | ||
| Ramsey David S. L., McMaster Damien, Thomas Elaine | Alpine NP - Eastern Alps | 2020 | 2022 | 1.40 | 1.08 | 1.76 | |||
| Ramsey David S. L., McMaster Damien, Thomas Elaine | Burrowa-Pine | 2020 | 2022 | 0.40 | 0.27 | 0.56 | |||
| Ramsey David S. L., McMaster Damien, Thomas Elaine | Coopracambra | 2020 | 2022 | 1.12 | 0.04 | 0.20 | |||
| Ramsey David S. L., McMaster Damien, Thomas Elaine | Croajingolong NP | 2020 | 2022 | 0.61 | 0.44 | 0.81 | |||
| Ramsey David S. L., McMaster Damien, Thomas Elaine | Errinundra NP | 2020 | 2022 | 0.80 | 0.58 | 1.05 | |||
| Ramsey David S. L., McMaster Damien, Thomas Elaine | Mount Mitta Mitta RP | 2020 | 2022 | 0.30 | 0.14 | 0.49 | |||
| Ramsey David S. L., McMaster Damien, Thomas Elaine | Mt Buffalo NP | 2020 | 2022 | 0.78 | 0.64 | 0.94 | |||
| Ramsey David S. L., McMaster Damien, Thomas Elaine | Snowy River NP | 2020 | 2022 | 1.08 | 0.86 | 1.28 | |||
| Ramsey David S. L., McMaster Damien, Thomas Elaine | Wabba WP | 2020 | 2022 | 1.30 | 1.02 | 1.60 | |||
| The role of wild deer in the transmission of diseases of livestock | Pacioni, Ramsey et al | Gembrook | 2018 | 2020 | 5.90 | 5.70 | 6.10 | ||
| Pacioni, Ramsey et al | Kinglake | 2018 | 2020 | 6.80 | 3.70 | 10.90 | |||
| Pacioni, Ramsey et al | Willow Grove | 2018 | 2020 | 6.20 | 5.90 | 6.50 |
Presence-absence records
The camera traps used in this study also detected numerous other species. We can obtain a table of the presence-absences of these species using the following code:
Show code
# Get view
pa <- tbl(con, dbplyr::in_schema("camtrap", "processed_site_substation_presence_absence")) %>%
dplyr::filter(ProjectShortName %in% !!project_short_name) %>% # Only retrieve cam records for hog deer 2023
dplyr::collect() %>%
dplyr::arrange(SiteID)
# Get pa table
pa_table <- pa %>%
group_by(Species = common_name) %>%
summarise(`Sites Detected` = sum(Presence)) %>%
arrange(desc(`Sites Detected`))
write.csv(pa_table, "outputs/tables/other_species.csv", row.names = F)
pa_table %>%
kbl("html", digits = 3, caption = "Other species and the number of sites (out of 317) that they were recorded at.") %>%
kable_styling("striped", font_size = 8, full_width = F) %>%
scroll_box()| Species | Sites Detected |
|---|---|
| Black-tailed Wallaby | 214 |
| Eastern Grey Kangaroo | 129 |
| Sambar Deer | 104 |
| Bare-nosed Wombat | 91 |
| Red Fox | 90 |
| Common Brush-tailed Possum | 80 |
| Emu | 46 |
| Dingo | 38 |
| Superb Lyrebird | 38 |
| Pied Currawong | 36 |
| Laughing Kookaburra | 35 |
| Short-beaked Echidna | 32 |
| Domestic Cat (feral) | 30 |
| Fallow Deer | 30 |
| Red-necked Wallaby | 25 |
| Ravens and Crows | 24 |
| Hog Deer | 22 |
| Grey Currawong | 19 |
| Australian Magpie | 18 |
| Southern Long-nosed Bandicoot | 17 |
| Western Grey Kangaroo | 17 |
| Long-footed Potoroo | 16 |
| Rabbits and Hares | 16 |
| European Rabbit | 15 |
| Red Deer | 12 |
| Red Wattlebird | 11 |
| Eastern Ring-tailed Possum | 9 |
| Koala | 9 |
| White-winged Chough | 9 |
| Wonga Pigeon | 9 |
| Grey Shrike-thrush | 8 |
| Lace Monitor | 8 |
| Crimson Rosella | 6 |
| Goat (feral) | 6 |
| Australian Raven | 5 |
| Bassian Thrush | 3 |
| Bush Rat | 3 |
| Common Blackbird | 3 |
| Horse (feral) | 3 |
| Leadbeater’s Possum | 3 |
| Magpie-lark | 3 |
| Mountain Brush-tailed Possum | 3 |
| Pacific Black Duck | 3 |
| Satin Bowerbird | 3 |
| Australian Wood Duck | 2 |
| Black Rat | 2 |
| Eastern Whipbird | 2 |
| Flame Robin | 2 |
| Long-nosed Potoroo | 2 |
| Pig (feral) | 2 |
| White-browed Scrubwren | 2 |
| White-eared Honeyeater | 2 |
| Yellow-tailed Black-Cockatoo | 2 |
| Antechinus | 1 |
| Australian Owlet-nightjar | 1 |
| Black Swan | 1 |
| Black-faced Cuckoo-shrike | 1 |
| Brush Bronzewing | 1 |
| Brushtail Possums | 1 |
| Cattle (feral) | 1 |
| Chestnut Teal | 1 |
| Dingo & Dog (feral) | 1 |
| Eastern Rosella | 1 |
| Eastern Yellow Robin | 1 |
| Galah | 1 |
| Grey Butcherbird | 1 |
| Grey Fantail | 1 |
| Morepork | 1 |
| Red Kangaroo | 1 |
| Rufous Fantail | 1 |
| Stubble Quail | 1 |
| Sulphur-crested Cockatoo | 1 |
| White-faced Heron | 1 |