I've added extra functions to the 'makeJAGSmask'
package since it first appeared in August 2016  see the old
post
here for details of the original functions. It's now possible to create a mask from a raster
object and to use the values of the raster as a spatial
covariate. You can also define a 'core' area smaller than the
habitat mask and get an estimate of the number of activity
centres in the core. You can download and install the latest version
(currently 0.1.1) of the package with
devtools::install_github("mikemeredith/makeJAGSmask")
Note that you will need to have packages abind ,
plotrix and raster
installed to run makeJAGSmask .
The idea of a 'core' area
Often suitable habitat extends outside the study area, and we
can't ignore the possibility that animals with activity centres
(ACs) outside the study area will appear in our traps. We
usually reckon that animals with ACs less than 4 \(\sigma\) from
any trap have some probability of capture, and work with a "state space"
which includes all suitable habitat within that distance. But
the abundance or density of animals in this larger state space
is not the figure we are looking for, and it involves
extrapolating estimates a long way from the traps, into areas
where we have little information other than habitat suitability.
The solution is to produce an estimate of the number or
density of ACs in the smaller area we are interested in, the
core area. The same methods can be used to estimate populations
in regions within the main study area, as in the example below.
Implementation
We'll demonstrate this using example data from the
makeJAGSmask package.
library(makeJAGSmask)
data(simSCR)
str(simSCR)
attach(simSCR) # remember to detach later
plot(JAGSmask)
Take a look at the posterior distributions of the AC
locations:
# Convert AC locations to the original units, remove phantoms:
AC < convertOutput(sims.list$S, JAGSmask)
AC[sims.list$w == 0] < NA
plotACs(which=1:30, AC=AC, traps=traps, Y=Y, hab=patchDF, show.labels=FALSE)
There is a river flowing through the study area, and the
portion north of the river lies in the municipality of Arkadia.
polygon(Arkadia, border='red')
Here Arkadia is a 2column data frame with the vertices of
the polygon forming the boundary, but a SpatialPolygons or
SpatialPolygonsDataFrame
object imported from a GIS file could also be used.
We add Arkadia to the mask as a core area (yellow in the plot
below), then use getACinCore to obtain an MCMC chain for the number,
N , of ACs in Arkadia.
coremask < addCore(JAGSmask, type="poly", poly=Arkadia)
plot(coremask)
N < getACinCore(sims.list$S, sims.list$w, coremask)
plot(table(N)) # posterior probability distribution
mean(N) # posterior mean
quantile(N, probs=c(0.025, 0.975)) # 95% credible interval
Of course it's possible to make more than one JAGSmask object
with different core areas and get population estimates for
different subsets of the state space.
Using rasters
The new version of makeJAGSmask has a function to generate a
1/0 habitat mask from an object of class RasterLayer
or RasterStack . This
ensures that the pixels of the 1/0 mask align exactly with those
of the raster. The raster itself is converted to a matrix which
can be used as a covariate in the analysis in JAGS.
We'll use the example data set above, where simSCR$patchRS is
a RasterStack :
library(raster)
plot(patchRS)
The "rand" layer is a random field while "edge" is the
distance in metres from the
edge of the patch. Nonhabitat pixels have NA values.
We will use the distancefromedge values as a spatial
covariate, but first we need to standardise to mean 0, SD 1. To
keep the code readable, we'll pull out the "edge" layer:
edge < patchRS$edge
patchRS$edgeS < (edge  mean(values(edge), na.rm=TRUE)) / sd(values(edge), na.rm=TRUE)
patchRS
Now we create our JAGSmask object which will
include the raster information as matrices:
mymask < convertRaster(patchRS, traps)
str(mymask)
Note that mymask has three covariate matrices
with the original names from the raster stack, but only the
first of these appears in the plot. If we had used a single
raster layer, the covariate would have been named covMat .
With only 5 animals captured (out of 6 total) we can't expect to
be able to estimate coefficients properly, but this example will
run quickly. We'll only use the edgeS covariate for
the analysis below.
In the JAGS model, we prepare a matrix, probs , with the
probability that an AC will occur in the corresponding pixel,
which is a function of the distancefromedge covariate and the
coefficient beta . Unfortunately, the link function exp
is not vectorised, so we have to do that step pixel by pixel. We ensure that the probability of landing in
nonhabitat is zero, and that the probabilities sum to one.
Candidate AC locations are generated as usual from a uniform
prior, we look up the probability of landing in that pixel in
the probs matrix, then use the 'zeros trick' (Lunn
et al, 2013, §9.5.2,
Kéry
& Royle, 2016, §5.8) to incorporate
this into the likelihood.
model {
sigma ~ dunif(0, 10) # set up the priors for the parameters
alpha < 1/(2*sigma^2)
p0 ~ dbeta(1,1)
omega ~ dbeta(0.001, 1)
beta ~ dunif(5, 5)
# spatial covariate
for(i in 1:(upperLimit[1]1)) {
for(j in 1:(upperLimit[2]1)) {
lam[i, j] < exp(beta * edgeS[i, j]) # 'exp' ensures 'lam' is nonnegative
}
}
lam0 < lam * habMat # convert 'lam' to 0 for nonhabitat
probs < lam0 / sum(lam0) # 'probs' must sum to 1
for (i in 1:M){ # loop through the augmented population
w[i] ~ dbern(omega) # state of individual 'i' (real or imaginary)
S[i, 1] ~ dunif(1, upperLimit[1]) # uniform priors for the activity centres for each individual
S[i, 2] ~ dunif(1, upperLimit[2])
# the zeros trick:
negLogDen[i] < log(probs[trunc(S[i,1]), trunc(S[i,2])])
zeros[i] ~ dpois(negLogDen[i]) # 'zeros' is a vector of zeros in the data
# loop through the camera trap locations
for(j in 1:nTraps) {
Dsq[i,j] < (S[i,1]trapMat[j,1])^2 + (S[i,2]trapMat[j,2])^2
p[i,j] < p0 * exp(Dsq[i,j] * alpha)
y[i,j] ~ dbin(p[i,j] * w[i], nOcc)
}
}
N < sum(w) # derive number (check that N << M)
}
We save the model in the file "patch_covar.jags", then
prepare the other inputs and do a short run:
M < 25 # size of the augmented data set
nTraps < nrow(traps)
yAug < rbind(Y, matrix(0, M  nrow(Y), nTraps))
# organise data for JAGS
dataList < c(mymask, list(y = yAug, M = M, nTraps = nTraps, zeros = rep(0, M), nOcc = 90))
str(dataList)
# Initial values
inits < function() list(w = rep(1, M), S = matrix(15, M, 2))
wanted < c("N", "omega", "p0", "beta", "sigma", "S", "w")
library(jagsUI)
result < jags(dataList, inits, wanted, "patch_covar.jags", DIC=FALSE,
n.chains=3, n.adapt=1000, n.iter=5500, n.burnin=500, n.thin=10, parallel=TRUE)
This takes about 4 mins to run and convergence is good,
though n.eff 's are small. The randomlyplaced ACs did tend to be near the edge
as can be seen in the plot in the previous section,
so an apparent negative effect of distancefromedge is expected.
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
N 7.049 2.079 5.000 6.000 7.000 8.000 12.000 1.004 670
omega 0.271 0.115 0.098 0.185 0.258 0.337 0.544 1.002 728
p0 0.023 0.019 0.006 0.012 0.018 0.027 0.080 1.050 220
beta 0.742 1.032 3.474 1.214 0.544 0.036 0.712 1.010 1500
sigma 2.922 0.664 1.996 2.432 2.809 3.266 4.586 1.001 1500
...
