Spatial capture-recapture model

#1

I thought I’d post this code for fitting a spatial capture-recapture model in case it proves useful to someone else.

Simulate some data:

## Simulate activity centers for N=100 individuals
N <- 100
xlim <- ylim <- c(0,1)
s <- cbind(runif(N, xlim[1], xlim[2]),
           runif(N, ylim[1], ylim[2]))

## Trap locations
x0 <- seq(0.25, 0.75, length=7)
x <- cbind(rep(x0, each=7), rep(x0, times=7))
J <- nrow(x) ## nTraps

if(1==2) {
    plot(s, xlim=c(0,1), ylim=c(0,1))
    points(x, pch=3)
}

## Parameters of encounter model
lam0 <- 0.1
sigma <- 0.1

K <- 5 ## nOccasions
lambda <- dist <- matrix(NA, N, J)

## Encounter histories for all N individuals
yall <- array(NA, c(N, J)) ## Shout-out from Georgia
for(j in 1:J) {
    dist[,j] <- sqrt((s[,1]-x[j,1])^2 +
                     (s[,2]-x[j,2])^2)
    lambda[,j] <- lam0*exp(-dist[,j]^2 / (2*sigma^2))*K
    yall[,j] <- rpois(N, lambda[,j])
}

## Discard data on individuals not detected
detected.all <- rowSums(yall)>0
y <- yall[detected.all,]  ## Observed encounter histories
(n <- nrow(y))            ## nCaptured

Specify the greta model. Because greta can’t (yet?) sample from discrete distributions, the indicator variable z, which is used in data augmentation, must be marginalized from the likelihood. I’ve done this using the zeros trick.

library(greta)

x <- as_data(x)
y <- as_data(y)
M <- 200 ## Data augmentation size
z1.dat <- as_data(rep(1L, n))          ## We know z=1 for captured individuals
zeros.trick <- as_data(rep(0L, M-n))   ## Used to marginalize z out of likelihood

lam0 <- normal(0, 2, dim=1, truncation=c(0, Inf))
sigma <- normal(0, 2, dim=1, truncation=c(0, Inf))
psi <- beta(1, 1, dim=1)

## Activity centers
s <- cbind(uniform(xlim[1], xlim[2], dim=M),
           uniform(ylim[1], ylim[2], dim=M))

## Encounter rate
lambda <- lam0*exp(-rdist(s, x)^2 / (2*sigma*sigma))*K

## Marginal (neg log) likelihood for augmented individuals
nll <- -log(exp(-rowSums(lambda[(n+1):M,]))*psi + 1-psi)

## Distributions for the data
distribution(y) <- poisson(lambda[1:n,])
distribution(z1.dat) <- bernoulli(psi)
distribution(zeros.trick) <- poisson(nll)

system.time({
    gm <- model(lam0, sigma, psi)
}) ## Should be <1s

plot(gm)


system.time({
    mc <- mcmc(gm, n_samples=5000, warmup=1000, chains=1, n_cores=16)
}) ## Could take a few minutes


plot(mc)

coda::effectiveSize(mc)

summary(window(mc, start=1001))

Since z was marginalized, we have to recover it after the fact in order to get the posterior for N.

## First, compute the conditional probability: Pr(z=1|y=0)
cp.z1 <- exp(-rowSums(lambda[(n+1):M,]))*psi / exp(-nll)
cp.z1.post <- as.matrix(calculate(cp.z1, mc))

## Now sample z for the augmented individuals (who implicitly have y=0)
z.post <- matrix(NA_integer_, nrow(cp.z1.post), ncol(cp.z1.post))
z.post[] <- rbinom(n=prod(dim(z.post)), 1, cp.z1.post)

## Here's N
N.post <- rowSums(z.post)+n
hist(N.post)

This resolves the issue that I posted earlier. Apparently, the model compiles and runs much much faster if you avoid for loops.

If anyone sees opportunities for speeding this up further, please let me know.

Greta Sampling seems very slow
Imperfect detection abundance model
#2

Ha! I should have read this before replying over there. Exactly the same approach though.

Only recommendation for speedup is to run multiple chains. greta can run those simultaneously, and therefore run them in parallel too. So that’s a quicker way to get lots of samples. In fact running lots of chains in greta often leads to faster tuning of sampler parameters (greta shares information between samplers during the warmup phase only), and faster sampling.

There’s an identifiability issue with the activity centre locations with multiple chains, but that shouldn’t affect convergence diagnostics on the other parameters.

#3

Thanks Nick. It’s great how easy it is to run multiple chains in parallel. I might be able to speed things up even more by integrating out the activity centers and by avoiding data augmentation, but that would make it harder to extend the model and embed more interesting point processes.

Can you explain the identifiability issue related to the activity centers? I just monitored them using a model with a single chain and a model with multiple chains, and everything looked fine to me.

Also, are you still planning on adding functionality to sample discrete latent variables? There are many types of SCR models (especially those involving unmarked individuals) in which there are complex dependencies among the latent variables, making it impossible to compute the marginal likelihood. JAGS can handle these models, but it’s quiet slow and I’d rather use greta.

#4

Sorry, I was thinking about unmarked models when I mentioned the identifiability issue. In those models you can swap any two activity centre locations and get the same result, but here you have marked individuals so no issue.

Yes, sampling discrete variables is on the to do list, and hopefully something I’ll have a prototype for in the next 3 or 4 weeks (see dev sprint plan here - just finishing off the simulation interface now). I’m planning to introduce some generic samplers for discrete variables (e.g. a discrete slice sampler), and handle models with both discrete and continuous variables by alternating between the discrete sampling and HMC for continuous parameters. That won’t retain all the nice geometric properties (and therefore efficient sampling) of HMC, but it should still be more efficient than Gibbs-sampling everything.

I’ve also got a prototype user interface for marginalising discrete variables that I’m planning to release soonish. Would be keen to get your feedback on that too!

Nick

#5

That’s great news. Let me know if you need any help testing the discrete variable functionality.