multiscale-occupancy.R

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

psi <- .2
theta <- .3
p <- .5

nsite <- 100
ntime <- 4
nrep <- 3

z <- rbinom(nsite, 1, psi)
alpha <- matrix(rbinom(nsite * ntime, 1, theta), nrow = nsite)

y <- array(dim = c(nsite, ntime, nrep))
for (i in 1:nsite) {
  for (j in 1:ntime) {
    y[i, j, ] <- rbinom(nrep, 1, z[i] * alpha[i, j] * p)
  }
}

# potential combinations of alpha that can lead to all-zero capture history
alpha_potential <- expand.grid(rep(list(c(0, 1)), ntime))

stan_d <- list(nsite = nsite, 
               ntime = ntime, 
               nrep = nrep, 
               y = y, 
               n_possible = 2^ntime, 
               alpha_potential = alpha_potential)

m_init <- stan_model("multiscale-occupancy.stan")
m_fit <- sampling(m_init, data = stan_d)
pairs(m_fit, pars = c("psi", "theta", "p"))

multiscale-occupancy.stan

data {
  int<lower = 1> nsite;
  int<lower = 1> ntime;
  int<lower = 1> nrep;
  
  int<lower = 0, upper = 1> y[nsite, ntime, nrep];
  
  int<lower = 1> n_possible;
  matrix<lower = 0, upper = 1>[n_possible, ntime] alpha_potential;
}

transformed data {
  int<lower = 0, upper = 1> known_present[nsite];
  int<lower = 0, upper = 1> known_available[nsite, ntime];
  
  for (i in 1:nsite) {
    known_present[i] = 0;
    for (j in 1:ntime) {
      known_available[i, j] = 0;
      for (k in 1:nrep) {
        if (y[i, j, k] == 1) {
          known_present[i] = 1;
          known_available[i, j] = 1;
        }
      }
    }
  }
}

parameters {
  real<lower = 0, upper = 1> psi;
  real<lower = 0, upper = 1> theta;
  real<lower = 0, upper = 1> p;
}

transformed parameters {
  vector[nsite] log_lik;
  
  {
    vector[ntime] tmp_lp;
    matrix[n_possible, ntime] tmp_poss;
    vector[n_possible + 1] sum_poss;
    
    for (i in 1:nsite) {
      if (known_present[i]) {
        for (j in 1:ntime) {
          if (known_available[i, j]) {
            // present and available
            tmp_lp[j] = log(theta) + bernoulli_lpmf(y[i, j, ] | p);
          } else {
            // present, possibly unavailable
            tmp_lp[j] = log_sum_exp(
              log(theta) + bernoulli_lpmf(y[i, j, ] | p), 
              log1m(theta)
            );
          }
        }
        log_lik[i] = log(psi) + sum(tmp_lp);
      } else {
        // could be present or absent (was never detected)
        // and there are 2^ntime possible combinations
        // of alpha_{i, t} that are relevant if z_i = 1
        for (jj in 1:n_possible) {
          for (j in 1:ntime) {
            if (alpha_potential[jj, j] == 0) {
              // not available
              tmp_poss[jj, j] = log1m(theta);
            } else {
              // available but not detected
              tmp_poss[jj, j] = log(theta) + bernoulli_lpmf(y[i, j, ] | p);
            }
          }
          sum_poss[jj] = log(psi) + sum(tmp_poss[jj, ]);
        }
        sum_poss[n_possible + 1] = log1m(psi);
        log_lik[i] = log_sum_exp(sum_poss);
      }
    }
  }
}

model {
  target += sum(log_lik);
}