packages
library(raster)
## Loading required package: sp
covariate for distribution of possums
tree_cover <- readRDS('data/nveg_ras.rds')
# extract center and scale for full region
tree_cover_scaled <- scale(getValues(tree_cover))
scale <- attr(tree_cover_scaled, 'scaled:scale')
center <- attr(tree_cover_scaled, 'scaled:center')
PRESENCE-ONLY DATA
po_raw <- read.csv('data/po.csv')
po_ras <- rasterize(po_raw, tree_cover,
fun = 'count', background = 0)
po <- getValues(po_ras)
table(po)
## po
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 4387 375 165 76 66 27 15 8 5 5 4 6 5 2 1 2
## 16 17 19 21 25 26 27 44 95
## 1 1 1 1 1 2 1 1 1
# vis data
plot(tree_cover)
points(SpatialPoints(po_raw), pch = 20)
# area for po data is sqkm
area_ras <- area(tree_cover)
area_po <- getValues(area_ras)
# covariates for po data
tree_cover_po <- as.vector(tree_cover_scaled)
access_ras <- readRDS('data/access_ras.rds')
access <- as.vector(scale(getValues(access_ras)))
COUNT DATA
counts <- read.csv('data/counts.csv')
counts_sp <- SpatialPoints(counts[, c('lon', 'lat')])
# area of count data is search radius 50m. convert to sqkm here
area_counts <- (pi * 50 ^ 2) / 1e6
tree_cover_counts <- as.vector(scale(
extract(tree_cover,
counts_sp),
center = center,
scale = scale))
DETECTION DATA
scat_det <- read.csv('data/scat_det.csv')
scat_det_sp <- SpatialPoints(scat_det[, c('lon', 'lat')])
# scat search radius 10m
area_scat_det <- (pi * 10 ^ 2) / 1e6
tree_cover_scat_det <- as.vector(scale(
extract(tree_cover,
scat_det_sp),
center = center,
scale = scale))
predictions across whole region
tree_cover_pred <- as.vector(tree_cover_scaled)
access_pred <- access
create jags data
win.data <- list(
# presence-only data
n_po = length(po),
po = po,
tree_cover_po = tree_cover_po,
area_po = area_po,
access = access,
# spotlight count data
n_counts = length(counts),
counts = counts$count,
tree_cover_counts = tree_cover_counts,
area_counts = area_counts,
# time to detection scat data
n_scat_det = length(scat_det),
tree_cover_scat_det = tree_cover_scat_det,
area_scat_det = area_scat_det,
# prediction
n_pred = length(tree_cover_pred),
tree_cover_pred = tree_cover_pred,
access_pred = access_pred
)
sink("jags_mod.txt")
cat("
# Define model and write into R working directory
model {
## Priors ##
# Priors for abundance of possums
alpha ~ dnorm(0, 0.1)
# Prior for slope of abundance predictor
beta ~ dnorm(0, 0.1)
# Prior for possum spatial sampling bias
alpha_bias ~ dnorm(0, 0.1)
# Prior for slope in bias predictors
beta_bias ~ dnorm(0, 0.1)
# Prior for scat detection
alpha_scat_det ~ dnorm(0, 0.1)
## Likelihood ##
# count data of possums
for (i in 1:n_counts){
counts[i] ~ dpois(lambda_1[i]*area_counts)
log(lambda_1[i]) <- alpha + beta*tree_cover_counts[i]
}
# presence-only possum data
for (i in 1:n_po){
po[i] ~ dpois(lambda_2[i]*area_po[i]*bias[i])
log(lambda_2[i]) <- alpha + beta*tree_cover_po[i]
# access is the covariate for bias
log(bias[i]) <- alpha_bias + beta_bias*access[i]
}
# scat detection data
for (i in 1:n_scat_det) {
scat_det[i] ~ dbern(prob_scat_det[i])
# prob_scat_det is prob osberving 1 or more scats given expected abund scats
prob_scat_det[i] <- 1 - exp(-abund_scat[i])
# alpha_scat_det scaling param to relate abund possums to abund scats
abund_scat[i] <- alpha_scat_det*abund_possum[i]
abund_possum[i] <- lambda_3[i]*area_scat_det
log(lambda_3[i]) <- alpha + beta*tree_cover_scat_det[i]
}
## Derived quantities ##
for (i in 1:n_pred){
log(lambda_pred[i]) <- alpha + beta*tree_cover_pred[i]
log(bias_pred[i]) <- alpha_bias + beta_bias*access_pred[i]
}
}
", fill = TRUE)
##
##
## # Define model and write into R working directory
##
## model {
##
## ## Priors ##
##
## # Priors for abundance of possums
## alpha ~ dnorm(0, 0.1)
##
## # Prior for slope of abundance predictor
## beta ~ dnorm(0, 0.1)
##
## # Prior for possum spatial sampling bias
## alpha_bias ~ dnorm(0, 0.1)
##
## # Prior for slope in bias predictors
## beta_bias ~ dnorm(0, 0.1)
##
## # Prior for scat detection
## alpha_scat_det ~ dnorm(0, 0.1)
##
## ## Likelihood ##
##
## # count data of possums
## for (i in 1:n_counts){
##
## counts[i] ~ dpois(lambda_1[i]*area_counts)
## log(lambda_1[i]) <- alpha + beta*tree_cover_counts[i]
##
## }
##
## # presence-only possum data
## for (i in 1:n_po){
##
## po[i] ~ dpois(lambda_2[i]*area_po[i]*bias[i])
## log(lambda_2[i]) <- alpha + beta*tree_cover_po[i]
##
## # access is the covariate for bias
## log(bias[i]) <- alpha_bias + beta_bias*access[i]
##
## }
##
## # scat detection data
## for (i in 1:n_scat_det) {
##
## scat_det[i] ~ dbern(prob_scat_det[i])
##
## # prob_scat_det is prob osberving 1 or more scats given expected abund scats
## prob_scat_det[i] <- 1 - exp(-abund_scat[i])
##
## # alpha_scat_det scaling param to relate abund possums to abund scats
## abund_scat[i] <- alpha_scat_det*abund_possum[i]
##
## abund_possum[i] <- lambda_3[i]*area_scat_det
## log(lambda_3[i]) <- alpha + beta*tree_cover_scat_det[i]
##
## }
##
## ## Derived quantities ##
##
## for (i in 1:n_pred){
##
## log(lambda_pred[i]) <- alpha + beta*tree_cover_pred[i]
## log(bias_pred[i]) <- alpha_bias + beta_bias*access_pred[i]
##
## }
##
## }
sink()
MCMC
# Run JAGS, check convergence and summarize posteriors:
# out <- jagsUI::jags(data = win.data, inits = NULL,
# parameters.to.save = c("beta", "alpha",
# "alpha_scat_det",
# "alpha_bias", "beta_bias",
# "lambda_pred"),
# model.file = "jags_mod.txt",
# n.chains = 3, n.thin = 1, n.iter = 5000, n.burnin = 2500,
# codaOnly = TRUE, parallel = FALSE)
# saveRDS(out, "jags_mod.rds")
out <- readRDS('jags_mod.rds')
prediction
pred_map <- tree_cover
values(pred_map) <- as.vector(out$mean$lambda_pred)
pred_map <- mask(pred_map, tree_cover)
plot(pred_map)