The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
In this vignette, we will use the nimbleSCR (Bischof et al. 2020) package and NIMBLE (de Valpine et al. 2017; NIMBLE Development Team 2020) to simulate and fit an efficient Bayesian SCR model.
For this example, we create a \(12*12\) habitat grid represented by grid cell centers. We center a \(8*8\) trapping grid leaving a buffer area of 2 units around it.
# CREATE HABITAT GRID
coordsHabitatGridCenter <- cbind(rep(seq(11.5, 0.5, by=-1), 12),
sort(rep(seq(0.5, 11.5, by=1), 12)))
colnames(coordsHabitatGridCenter) <- c("x","y")
# CREATE TRAP GRID
trapCoords <- cbind(rep(seq(2.5, 9.5,by=1),8),
sort(rep(seq(2.5, 9.5,by=1),8)))
colnames(trapCoords) <- c("x","y")
# PLOT CHECK
plot(coordsHabitatGridCenter[,"y"] ~ coordsHabitatGridCenter[,"x"], pch = 1, cex = 1.5) #pch=16)
points(trapCoords[,"y"] ~ trapCoords[,"x"], col="red", pch=16 )
par(xpd=TRUE)
legend(x = 7, y = 13,
legend=c("Habitat grid centers", "Traps"),
pt.cex = c(1.5,1),
horiz = T,
pch=c(1,16),
col=c("black", "red"),
bty = 'n')
Here we rescale the trap coordinates to the habitat to allow a fast habitat cell ID lookup and the local evaluation (see Milleret et al. (2019) and Turek et al. (2021) for further details).
ScaledtrapCoords <- scaleCoordsToHabitatGrid(coordsData = trapCoords,
coordsHabitatGridCenter = coordsHabitatGridCenter)$coordsDataScaled
habitatMask <- matrix(1, nrow = 12, ncol= 12, byrow = TRUE)
We also set up the objects necessary to perform the local evaluation using the ‘getLocalObjects’ function. Special care should be taken when choosing ‘dmax’ relative to \(\sigma\) (Milleret et al. 2019). Here we are using a value \(>2*\sigma\) (see below for the \(\sigma\) chosen.
We use the model code described in Turek et al. (2021). Note that we do not provide ‘detNums’ and ‘detIndices’ arguments in the ‘dbinomLocal_normal’ function as we wish to use this model to simulate data (see ‘?dbinomLocal_normal’ for further details).
modelCode <- nimbleCode({
## priors
psi ~ dunif(0, 1)
sigma ~ dunif(0, 50)
p0 ~ dunif(0, 1)
## loop over individuals
for(i in 1:M) {
## AC coordinates
sxy[i,1] ~ dunif(0, x.max)
sxy[i,2] ~ dunif(0, y.max)
## habitat constraint
ones[i] ~ dHabitatMask( s = sxy[i,1:2],
xmin = lowerCoords[1],
xmax = upperCoords[1],
ymin = lowerCoords[2],
ymax = upperCoords[2],
habitat = habitat.mx[1:y.max,1:x.max])
## latent dead/alive indicators
z[i] ~ dbern(psi)
## likelihood
y[i, 1:lengthYCombined] ~ dbinomLocal_normal(size = trials[1:n.traps],
p0 = p0,
s = sxy[i,1:2],
sigma = sigma,
trapCoords = trapCoords[1:n.traps,1:2],
localTrapsIndices = trapIndex[1:n.cells,1:maxNBDets],
localTrapsNum = nTraps[1:n.cells],
habitatGrid = habitatIDDet[1:y.maxDet,1:x.maxDet],
indicator = z[i],
lengthYCombined = lengthYCombined)
}
## derived quantity: total population size
N <- sum(z[1:M])
})
The model formulation uses data augmentation to derive N estimates (Royle and Dorazio 2012). We therefore need to choose the total number of individuals M (detected + augmented). Here we used 200. The expected total number of individuals present in the population is ‘M * psi’.
When simulating detections using this formulation of the SCR model in NIMBLE, all the information about detections (where and how many) is stored in ‘y’ in that order (See ?dbinomLocal_normal for more details.):
We now need to provide the maximum number of spatial recaptures that can be simulated per individual. We recommend using ‘trapLocal$numlocalindicesmax’ that defines the maximum number of traps available for detections when local evaluation is used. This will enable the simulation of as many spatial detections as allowed by the restrictions imposed by the local evaluation (defined by the ‘dmax’ argument from ‘getLocalObjects’). This means that the length of the ‘y’ observation vector for each individual is equal to the length of \(c(detNums, x, detIndices)\) and is therefore equal to \(lengthYCombined = 1+ trapLocal\$numLocalIndicesMax * 2\).
nimConstants <- list(M = M,
n.traps = dim(ScaledtrapCoords)[1],
y.max = dim(habitatMask)[1],
x.max = dim(habitatMask)[2],
y.maxDet = dim(trapLocal$habitatGrid)[1],
x.maxDet = dim(trapLocal$habitatGrid)[2],
n.cells = dim(trapLocal$localIndices)[1],
maxNBDets = trapLocal$numLocalIndicesMax,
trapIndex = trapLocal$localIndices,
nTraps = trapLocal$numLocalIndices,
habitatIDDet = trapLocal$habitatGrid,
lengthYCombined = lengthYCombined)
nimData <- list(trapCoords = ScaledtrapCoords,
habitat.mx = habitatMask,
ones = rep(1, nimConstants$M),
lowerCoords = c(min(coordsHabitatGridCenter[,1]) - 0.5, min(coordsHabitatGridCenter[,2]) - 0.5),
upperCoords = c(max(coordsHabitatGridCenter[,1]) + 0.5, max(coordsHabitatGridCenter[,2]) + 0.5),
trials = rep(1, dim(ScaledtrapCoords)[1])
)
# We set the parameter values as inits
nimInits <- list(p0 = p0,
psi = psi,
sigma = sigma)
We first need to obtain the list of nodes that will be simulated. We used the ‘getDependencies’ function from NIMBLE. Using the ‘simulate’ function from NIMBLE, we will then simulate the activity center (AC) locations (‘sxy’), the state of the individual (‘z’) and SCR observation data (‘y’) given the values we provided for ‘p0’, ‘sigma’ and ‘psi’.
# FIRST WE GET THE NODES TO SIMULATE
nodesToSim <- model$getDependencies(c("sxy", "z"), self=T)
# THEN WE SIMULATE THOSE NODES
set.seed(100)
model$simulate(nodesToSim, includeData = FALSE)
After running ‘simulate’, the simulated data are stored in the ‘model’ object. For example, we can access the simulated ‘z’ and check how many individuals were considered present:
We have simulated 108 individuals present in the population of which 86 were detected.
Here we also make some plots to check where are located the simulated detections in relation with the simulated location of the activity center for a given individual.
## [1] 6
#LET'S PLOT THEM
plot(coordsHabitatGridCenter[,"y"] ~ coordsHabitatGridCenter[,"x"], pch = 1, cex = 1.5) #pch=16)
points(trapCoords[,"y"] ~ trapCoords[,"x"], col="red", pch=16 )
# PLOT ITS ACTIVITY CETENR
points(model$sxy[i, 2] ~ model$sxy[i, 1],col="orange", pch=16)
whichdets <- model$y[i, (trapLocal$numLocalIndicesMax + 2):(trapLocal$numLocalIndicesMax+model$y[i, 1] + 1)]
points(ScaledtrapCoords[whichdets, "y"] ~
ScaledtrapCoords[whichdets, "x"], col="blue", pch=16)
par(xpd=TRUE)
legend(x = -1, y = 13,
legend=c("Habitat grid centers", "Traps", "Simulated AC", "Detections"),
pt.cex = c(1.5,1),
horiz = T,
pch=c(1, 16, 16, 16),
col=c("black", "red", "orange", "blue"),
bty = 'n')
Here, we build the NIMBLE model again using the simulated ‘y’ as data. For simplicity, we used the simulated ‘z’ as initial values. Then we can fit the SCR model with the simulated ‘y’ data set.
nimData1 <- nimData
nimData1$y <- model$y
nimInits1 <- nimInits
nimInits1$z <- model$z
nimInits1$sxy <- model$sxy
# CREATE AND COMPILE THE NIMBLE MODEL
model <- nimbleModel( code = modelCode,
constants = nimConstants,
data = nimData1,
inits = nimInits1,
check = F,
calculate = F)
model$calculate()
## [1] -1963.177
MCMCconf <- configureMCMC(model = model,
monitors = c("N", "sigma", "p0","psi"),
control = list(reflective = TRUE),
thin = 1)
## ===== Monitors =====
## thin = 1: N, p0, psi, sigma
## ===== Samplers =====
## binary sampler (200)
## - z[] (200 elements)
## RW sampler (403)
## - psi
## - sigma
## - p0
## - sxy[] (400 elements)
## Add block sampling of sxy coordinates see Turek et al 2021 for further details
MCMCconf$removeSamplers("sxy")
ACnodes <- paste0("sxy[", 1:nimConstants$M, ", 1:2]")
for(node in ACnodes) {
MCMCconf$addSampler(target = node,
type = "RW_block",
control = list(adaptScaleOnly = TRUE),
silent = TRUE)
}
MCMC <- buildMCMC(MCMCconf)
cMCMC <- compileNimble(MCMC, project = model)
niter <- 5000
nburnin <- 1000
nchains <- 3
# RUN THE MCMC
MCMCRuntime <- system.time(myNimbleOutput <- runMCMC( mcmc = cMCMC,
niter = niter,
nburnin = nburnin,
nchains = nchains,
samplesAsCodaMCMC = TRUE))
Now we imagine a scenario where we would like to build a SCR model with trap covariates on baseline detection probability (\(p_0\)). Compared to the model in (1.3), the baseline detection probability (\(p_0\)) is linearly modelled as a function of two covariates on the logistic scale such as:
\(logit(p0Traps{_j}) = logit(p_0) + BetaTraps_1 * trapCovs1_j + BetaTraps_2 * trapCovs2_j\)
First, we simulate values for two independent trap covariates.
set.seed(1)
trapCovs <- cbind( runif(nimConstants$n.traps,-1, 1),
runif(nimConstants$n.traps,-1, 1))
We now add the trap covariates to the data object.
We choose the simulated betaTraps values.
Because we are modelling the effect of trap covariates on baseline detection probability (\(p_0\)), this means that we have a \(p0Traps\) that varies for each trap j. When this is the case, we need to provide the ‘p0Traps’ argument for the ‘dbinomLocal_normal’ function, because the ‘p0’ argument only accepts a scalar.
modelCodeTrap <- nimbleCode({
## priors
psi ~ dunif(0, 1)
sigma ~ dunif(0, 50)
p0 ~ dunif(0, 1)
#trap covariates
betaTraps[1] ~ dunif(-5,5)
betaTraps[2] ~ dunif(-5,5)
p0Traps[1:n.traps] <- ilogit(logit(p0) + betaTraps[1] * trapCovs[1:n.traps, 1] +
betaTraps[2] * trapCovs[1:n.traps, 2])
## loop over individuals
for(i in 1:M) {
## AC coordinates
sxy[i,1] ~ dunif(0, x.max)
sxy[i,2] ~ dunif(0, y.max)
## habitat constraint
ones[i] ~ dHabitatMask( s = sxy[i,1:2],
xmin = lowerCoords[1],
xmax = upperCoords[1],
ymin = lowerCoords[2],
ymax = upperCoords[2],
habitat = habitat.mx[1:y.max,1:x.max])
## latent dead/alive indicators
z[i] ~ dbern(psi)
## likelihood
y[i, 1:lengthYCombined] ~ dbinomLocal_normal( size = trials[1:n.traps],
p0Traps = p0Traps[1:n.traps],
s = sxy[i,1:2],
sigma = sigma,
trapCoords = trapCoords[1:n.traps,1:2],
localTrapsIndices = trapIndex[1:n.cells,1:maxNBDets],
localTrapsNum = nTraps[1:n.cells],
habitatGrid = habitatIDDet[1:y.maxDet,1:x.maxDet],
indicator = z[i],
lengthYCombined = lengthYCombined)
}
## derived quantity: total population size
N <- sum(z[1:M])
})
# FIRST WE GET THE NODES TO SIMULATE
nodesToSim <- model$getDependencies(c("sxy", "z"), self=T)
# THEN WE SIMULATE THOSE NODES
set.seed(100)
model$simulate(nodesToSim, includeData = FALSE)
We have simulated 108 individuals present in the population of which 93 were detected.
Here we make some plots to check where are located the detections in relation with the location of activity center.
## [1] 3
#LET'S PLOT THEM
plot(coordsHabitatGridCenter[,"y"] ~ coordsHabitatGridCenter[,"x"], pch = 1, cex = 1.5) #pch=16)
points(trapCoords[,"y"] ~ trapCoords[,"x"], col="red", pch=16 )
# PLOT ITS ACTIVITY CETENR
points(model$sxy[i, 2] ~ model$sxy[i, 1],col="orange", pch=16)
whichdets <- model$y[i, (trapLocal$numLocalIndicesMax + 2):(trapLocal$numLocalIndicesMax+model$y[i, 1] + 1)]
points(ScaledtrapCoords[whichdets, "y"] ~
ScaledtrapCoords[whichdets, "x"], col="blue", pch=16)
par(xpd=TRUE)
legend(x = -1, y = 13,
legend=c("Habitat grid centers", "Traps", "Activity center", "Detections"),
pt.cex = c(1.5,1),
horiz = T,
pch=c(1, 16, 16, 16),
col=c("black", "red", "orange", "blue"),
bty = 'n')
nimData1 <- nimData
nimData1$y <- model$y
nimInits1 <- nimInits
nimInits1$z <- model$z
nimInits1$sxy <- model$sxy
# CREATE AND COMPILE THE NIMBLE MODEL
model <- nimbleModel( code = modelCodeTrap,
constants = nimConstants,
data = nimData1,
inits = nimInits1,
check = F,
calculate = F)
model$calculate()
## [1] -1976.644
MCMCconf <- configureMCMC(model = model,
monitors = c("N", "sigma", "p0","psi", "betaTraps"),
control = list(reflective = TRUE),
thin = 1)
## ===== Monitors =====
## thin = 1: N, betaTraps, p0, psi, sigma
## ===== Samplers =====
## binary sampler (200)
## - z[] (200 elements)
## RW sampler (405)
## - psi
## - sigma
## - p0
## - betaTraps[] (2 elements)
## - sxy[] (400 elements)
## Add block sampling of sxy coordinates see Turek et al 2021 for further details
MCMCconf$removeSamplers("sxy")
ACnodes <- paste0("sxy[", 1:nimConstants$M, ", 1:2]")
for(node in ACnodes) {
MCMCconf$addSampler(target = node,
type = "RW_block",
control = list(adaptScaleOnly = TRUE),
silent = TRUE)
}
MCMC <- buildMCMC(MCMCconf)
cMCMC <- compileNimble(MCMC, project = model, resetFunctions = TRUE)
# RUN THE MCMC
MCMCRuntime <- system.time(myNimbleOutput <- runMCMC( mcmc = cMCMC,
niter = niter,
nburnin = nburnin,
nchains = nchains,
samplesAsCodaMCMC = TRUE))
Bischof, Richard, Daniel Turek, Cyril Milleret, Torbjørn Ergon, Pierre Dupont, and de Valpine Perry. 2020. NimbleSCR: Spatial Capture-Recapture (Scr) Methods Using ’Nimble’.
de Valpine, Perry, Daniel Turek, Christopher J Paciorek, Clifford Anderson-Bergman, Duncan Temple Lang, and Rastislav Bodik. 2017. “Programming with Models: Writing Statistical Algorithms for General Model Structures with NIMBLE.” J. Comput. Graph. Stat. 26 (2): 403–13.
Milleret, Cyril, Pierre Dupont, Christophe Bonenfant, Henrik Brøseth, Øystein Flagstad, Chris Sutherland, and Richard Bischof. 2019. “A Local Evaluation of the Individual State-Space to Scale up Bayesian Spatial Capture–Recapture.” Ecology and Evolution 9 (1): 352–63.
NIMBLE Development Team. 2020. “NIMBLE: MCMC, Particle Filtering, and Programmable Hierarchical Modeling.” https://doi.org/10.5281/zenodo.1211190.
Royle, J Andrew, and Robert M Dorazio. 2012. “Parameter-Expanded Data Augmentation for Bayesian Analysis of Capture–Recapture Models.” Journal of Ornithology 152 (2): 521–37.
Turek, Daniel, Cyril Milleret, Torbjørn Ergon, Henrik Brøseth, Pierre Dupont, Richard Bischof, and Perry de Valpine. 2021. “Efficient Estimation of Large-Scale Spatial Capture–Recapture Models.” Ecosphere 12 (2): e03385.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.