Hi! I am exploring how to incorporate imperfect detection into abundance joint species distribution modelling after reading Tobler et al. (2019), which then led me to Yamaura et al. (2012). This may be partly related to Spatial capture-recapture model or State-space models in greta .
Below are two code chunks, one to simulate abundance data and another to fit model (apologies if they look ugly!) To jump ahead, I ran into an error while trying to implement Yamaura et al. (2012)'s model:
Error: model contains a discrete random variable that doesn't have a fixed value, so cannot be sampled from
I think it is because I cannot use poisson() for unobserved variable (commented in code below)? This is also hinted in the helpfile. Has anyone come up with a alternative solution? It would be very helpful and thanks in advance!
Simulate data
# Simulate data to test fit_model.R script
# Following the data generative process in
# Yamaura et al. (2012) Biodiversity and Conservation
set.seed(123)
# number of species
n_spp <- 10
# number of sites
n_site <- 100
# number of sampling events
t <- 10
# probability of detection for each species
# assuming the same p_j across sites for now
alpha0_true <- rnorm(n_spp, 2, 1)
mu_true <- alpha0_true
sigma_true <- rexp(n_spp, 2)
p_true <- numeric(n_spp)
for (j in seq_len(n_spp)) {
logit_p_true <- rnorm(1, mu_true, sigma_true)
p_true[j] <- exp(logit_p_true) / (exp(logit_p_true) + 1)
}
# true abundance
# using Poisson for now (explore neg. binom. later)
beta0_true <- rnorm(n_spp, 2, 1)
lambda_true <- exp(rep(1, n_site) %*% t(beta0_true))
N_true <- matrix(rpois(n_site * n_spp, lambda_true), n_site, n_spp)
# observed
n_true <- array(NA, dim = c(n_site, n_spp, t))
for (i in seq_len(n_site)) {
for (j in seq_len(n_spp)) {
n_true[i, j, ] <- rbinom(t, N_true[i, j], p_true[j])
}
}
Fit model
library(greta)
# source("code/sim_data.R")
# data
n <- as_data(n_true)
# variables
X <- ones(n_site)
V <- ones(n_site)
# priors
sd_beta0 <- exponential(2)
beta0 <- normal(0, sd_beta0, dim = n_spp)
lambda <- exp(X %*% t(beta0))
N <- poisson(lambda) # I think this doesn't work as per helpfile
N <- do.call(abind, c(replicate(t, N, simplify = FALSE), along = 3))
sd_alpha0 <- exponential(2)
alpha0 <- normal(0, sd_alpha0, dim = n_spp)
mu <- V %*% t(alpha0)
sigma <- exponential(2, dim = n_spp)
sigma <- do.call(rbind, replicate(n_site, t(sigma), simplify = FALSE))
# p <- ilogit(multivariate_normal(mu, sigma))
p <- ilogit(normal(mu, sigma))
p <- do.call(abind, c(replicate(t, p, simplify = FALSE), along = 3))
# observational model
distribution(n) <- binomial(N, p)
mod <- model(N, p)
