Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Code for the post about SECR with upaired cameras at http://www.mikemeredith.net/blog/2017/SECR_unpaired_cameras.htm
# SECR with irregular block of habitat with parallel processing
# Several possible scenarios:
#
# A: single cameras, ca. 10 clusters of 8,
# analyse only data for flank with largest number of captures = toss out 1 flank
# B: single cameras, ca. 10 clusters of 8,
# analyse data for both flanks as separate sessions = two sessions
# C: paired cameras, ca. 10 clusters of 8 (so double the number of cameras) = many pairs
# D: paired cameras, ca. 10 clusters of 4 (so same number of cameras) = few pairs
#
library(secr)
library(spatstat)
library(beeswarm) # For the beeswarm plot
## PRELIMINARIES
## =============
# Create the habitat patch polygon:
patch <- data.frame(x = c(130455, 132789, 134476, 136032, 138497,
139794, 144204, 144594, 144983, 145761, 147577, 150820, 155749,
156398, 159122, 160289, 162494, 165218, 165478, 167553, 167812,
168980, 170537, 172482, 174687, 176763, 178060, 179617, 181562,
183638, 183767, 184935, 185194, 185065, 179876, 178968, 179876,
181303, 184416, 186881, 189086, 189215, 187659, 187270, 186362,
183897, 182730, 182211, 182470, 183249, 183508, 183249, 182859,
181173, 179227, 178190, 178319, 176114, 175855, 175336, 174428,
173650, 173779, 172482, 171185, 168331, 167683, 168331, 168980,
168202, 166645, 164829, 162883, 160159, 159640, 158732, 158603,
157824, 154971, 154063, 151987, 151209, 146669, 143686, 142648,
141221, 139924, 137848, 136162, 136292, 134995, 133697, 131881,
130844, 129676, 128509, 126044, 124488, 124358, 124747, 125525,
126174, 125525, 124617, 125007, 124488, 125655, 127731, 128249,
129287, 129287, 128898, 130065, 130325, 129936, 130455),
y = c(802999,
803388, 802220, 802350, 805852, 806371, 805593, 804426, 798588,
798588, 797940, 798329, 804815, 805334, 806242, 805982, 809874,
812209, 814154, 816878, 819991, 821937, 821548, 820381, 818824,
818175, 815841, 812857, 809744, 806631, 804426, 802350, 799626,
798070, 795605, 793011, 790935, 788990, 787303, 786395, 786136,
785358, 784190, 782504, 782115, 782504, 781855, 780688, 778353,
775240, 771737, 765641, 763306, 762398, 762917, 763825, 765771,
770570, 772775, 773683, 773553, 771867, 766289, 765641, 766030,
767976, 769403, 770051, 771867, 773683, 774980, 775629, 777704,
780039, 779780, 779131, 777575, 776018, 775110, 773553, 773553,
772905, 772775, 773683, 775759, 777185, 776667, 775369, 773294,
771997, 771478, 771608, 771348, 772775, 774202, 774461, 774202,
774721, 777185, 783671, 784450, 785358, 786266, 787433, 789249,
790416, 791195, 790935, 792103, 792751, 794178, 795475, 797032,
800794, 801702, 802999))
# A function to calculate the number of captures and the number of relocations,
# ie, captures of an animal in a new location
getCapsRelocs <- function(CH) {
# function to do 1 session
getCapsRelocs1 <- function(ch) {
nCap <- nrow(ch)
if(nCap > 0) {
tmp <- apply(ch, c(1,3), sum)
reLoc <- sum(rowSums(tmp > 0) > 1) # number of relocations
} else {
reLoc <- NA
}
return(c(nCap=nCap, reLoc=reLoc))
}
if(is.list(CH)) {
nSes <- length(CH)
out <- matrix(NA_real_, nSes, 2)
colnames(out) <- c("nCap", "reLoc")
for(ses in 1:nSes)
out[ses, ] <- getCapsRelocs1(CH[[ses]])
} else {
out <- getCapsRelocs1(CH)
}
return(out)
}
## WORK THROUGH ONE SIMULATION
## ===========================
# Check the habitat patch
head(patch)
MASS::eqscplot(patch, type='l') # Small "L", not "1" one
# Convert to 'owin' object
win <- owin(poly=patch[115:1, ])
plot(win) # check
( area <- area.owin(win) / 1e4 ) # ha
area / 100 # km2
# Data generation
# ---------------
# Biological model
# ''''''''''''''''
# Generate activity centres
N <- 20
( D <- N / area ) # animals / ha
( DD <- round(D *100*100, 3) ) # animals / 100 sq km
AC <- rSSI(3000, N, win) # ACs at least 3km apart
points(AC, pch=19, col='red') # add to the plot
## Observation model
## '''''''''''''''''
# Lots of trap clusters spread out
trapSpacing <- 2500
# for A, B, and C, each cluster 3 x 3 but hollow
cluster <- make.grid(3, 3, spacing=trapSpacing, hollow = TRUE,
detector='proximity')
traps <- make.systematic(c(5, 4), cluster, patch)
plot(traps, add=TRUE)
nrow(traps) # number of traps
# For scenario D
clusterD <- make.grid(2, 2, spacing=trapSpacing, detector='proximity')
trapsD <- make.systematic(c(5, 4), clusterD, patch)
points(trapsD, col= 'blue', pch=4)
# Generate capture histories
g0 <- 0.02
sigma <- 2000
nOcc <- 90 # number of trapping occasions
popn <- data.frame(x=AC$x, y=AC$y)
class(popn) <- c("popn", "data.frame")
# Both Left and Right captured (scenario C):
simCH_C <- sim.capthist(traps, popn=popn,
detectpar=list(g0=g0, sigma=sigma), noccasions=nOcc)
# Split randomly into 2
is.L <- array(rbinom(prod(dim(simCH_C)), 1, 0.5), dim=dim(simCH_C)) # random 1/0 array
simCHL <- simCH_C * is.L
simCHL <- subset(simCHL, dropnullCH = TRUE) # Keep only flanks captured
simCHR <- simCH_C * (1 - is.L)
simCHR <- subset(simCHR, dropnullCH = TRUE)
sum(simCHL, simCHR) == sum(simCH_C) # check total number of capture events
( nCapL <- nrow(simCHL) )
( nCapR <- nrow(simCHR) )
# For scenario A, use whichever flank has most captures:
simCH_A <- if (nCapR > nCapL) simCHR else simCHL
# For scenario B, combine into a multisession capture history:
simCH_B <- MS.capthist(L = simCHL, R = simCHR)
# Scenario D, use the smaller trap array
simCH_D <- sim.capthist(trapsD, popn=popn,
detectpar=list(g0=g0, sigma=sigma),
noccasions=nOcc)
nrow(simCH_D)
# Data analysis
# -------------
# Create a fairly sparse mask to reduce run time.
mask <- make.mask(traps, buffer=25000, nx=32, ny=32, poly=patch)
plot(mask, add=TRUE)
# Starting values may be necessary for some simulated data sets
secrStart <- c(D=log(D), g0=qlogis(g0), sigma=log(sigma))
fitA <- secr.fit(simCH_A, mask=mask, start=secrStart) # A = toss out one flank
predict(fitA)
predict(fitA)[1, -1]*100*100 # Convert to animals per 100km2
DD
fitB <- secr.fit(simCH_B, mask=mask, start=secrStart) # B = two sessions
predict(fitB)
predict(fitB)[[1]][1, -1]*100*100 # Convert to animals per 100km2
DD
fitC <- secr.fit(simCH_C, mask=mask, start=secrStart) # C = many pairs
predict(fitC)
predict(fitC)[1, -1]*100*100 # Convert to animals per 100km2
DD
fitD <- secr.fit(simCH_D, mask=mask, start=secrStart) # D = few pairs
predict(fitD)
predict(fitD)[1, -1]*100*100 # Convert to animals per 100km2
DD
# DO HUNDREDS
# ===========
# Set up the parallel stuff:
library(foreach)
library(parallel) # parallel is part of the R core distribution
library(doParallel)
( ncore <- detectCores() - 1 ) # I like to keep 1 thread free
( cl <- makeCluster(ncore) )
clusterSetRNGStream(cl = cl) # set up random number stream
registerDoParallel(cl) # register as backend for 'foreach'
# The bits needed to run the loop:
# be sure to run patch and getCapsRelocs
win <- owin(poly=patch[115:1, ])
area <- area.owin(win) / 1e4 # ha
trapSpacing <- 2500
cluster8 <- make.grid(3, 3, spacing=trapSpacing, hollow = TRUE, detector='proximity')
cluster4 <- make.grid(2, 2, spacing=trapSpacing, detector='proximity')
N <- 20
D <- N / area # animals / ha
g0 <- 0.02
sigma <- 2000
nOcc <- 90
trapsTmp <- make.systematic(c(5, 4), cluster8, patch) # only for mask
mask <- make.mask(trapsTmp, buffer=25000, nx=32, ny=32, poly=patch)
secrStart <- c(D=log(D), g0=qlogis(g0), sigma=log(sigma))
seeds <- 1:301 # I was using 7 cores, so 43 per core = 301
system.time(
result <- foreach(x=seeds, .combine=rbind, .packages=c("secr", "spatstat")) %dopar% {
set.seed(x)
# Trap layouts
traps <- make.systematic(c(5, 4), cluster8, patch) # scenarios A, B, C
trapsD <- make.systematic(c(5, 4), cluster4, patch) # scenario D
outTraps <- c(nTrapsABC = nrow(traps), nTrapsD = nrow(trapsD))
# Population, ie, activity centre locations (same for all scenarios)
AC <- rSSI(3000, N, win) # ACs at least 3km apart
popn <- data.frame(x=AC$x, y=AC$y)
class(popn) <- c("popn", "data.frame")
# List to hold capture histories (allows us to use a loop for the model fitting)
CHlist <- vector("list", 4)
names(CHlist) <- c("A", "B", "C", "D")
# Capture histories for scenario C, "many pairs"
CHlist$C <- simCH_C <- sim.capthist(traps, popn=popn,
detectpar=list(g0=g0, sigma=sigma), noccasions=nOcc)
# Split into 2 for use with scenarios A and B
is.L <- array(rbinom(prod(dim(simCH_C)), 1, 0.5), dim=dim(simCH_C)) # random 1/0 array
simCHL <- simCH_C * is.L
simCHL <- subset(simCHL, dropnullCH = TRUE) # Keep only flanks captured
simCHR <- simCH_C * (1 - is.L)
simCHR <- subset(simCHR, dropnullCH = TRUE)
nCapL <- nrow(simCHL)
nCapR <- nrow(simCHR)
CHlist$A <- if (nCapR > nCapL) simCHR else simCHL # scenario A, "toss out"
CHlist$B <- MS.capthist(simCHL, simCHR) # scenario B, "2 sessions"
# capture histories for scenario D, "few pairs"
CHlist$D <- sim.capthist(trapsD, popn=popn,
detectpar=list(g0=g0, sigma=sigma), noccasions=nOcc)
# Analyse each of the capture histories
outMat <- matrix(NA_real_, 6, 4)
# rownames(outMat) <- c("nCaps", "nRelocs", "Dest", "Dse", "Dlow", "Dupp")
# colnames(outMat) <- names(CHlist) # names are stripped out before returning anyway
for(ch in 1:4) {
# Check the number of animals captured, and relocations:
capRel <- getCapsRelocs(CHlist[[ch]])
if(is.matrix(capRel))
capRel <- colSums(capRel)
outMat[1:2, ch] <- capRel
if(capRel[1] > 0 && capRel[2] > 0) { # proceed with analysis
fit <- secr.fit(CHlist[[ch]], mask=mask, start=secrStart)
pred <- predict(fit)
if(!is.data.frame(pred)) # multisession fit gives a list of data frames
pred <- pred[[1]]
outMat[3:6, ch] <- as.matrix(pred[-1])[1, ]*100*100 # extract the density estimate
}
}
c(outTraps, as.vector(outMat))
} ) # Took 1 hr
# save(result, file="OneSide_sims_20170205.Rda")
dim(result) # check
# Give the columns informative names
colnames(result) <- c("nTrapsABC", "nTrapsD",
t(outer(c("A", "B", "C", "D"), c("nCaps", "nRelocs", "Dest", "Dse", "Dlow", "Dupp"),
paste, sep="_")))
colMeans(result, na.rm=TRUE)
( DD <- D*100*100) # true density, animals per 100 sq km
colSums(is.na(result)) # How many failed?
mean(result[, 1]) # average numbers of traps
mean(result[, 2])
density <- result[, c(5, 11,17,23)] # per 100km2
colMeans(density, na.rm=TRUE)
# Calculate bias and RMSE
( bias <- colMeans(density, na.rm=TRUE) - DD )
getRmse <- function(x, true)
sqrt(mean((x - true)^2, na.rm=TRUE))
( RMSE <- apply(density, 2, getRmse, true = DD) )
# Organise data for plotting
density0 <- density
density0[is.na(density0)] <- -0.1 # Replace NAs with small neg value for plotting
densityDF <- as.data.frame(density)
densityDF0 <- as.data.frame(density0)
beeswarm(densityDF0, method='hex', pch=19, cex=0.5, las=1,
pwcol=c('blue', 'red')[(density0 < 0) +1], ylim=c(-0.2, 5),
labels=c("A:Toss out", "B:Two sessions", "C:Many pairs", "D:Few pairs"),
main="Estimates of density", ylab="Density (animals per 100 sq km)")
boxplot(densityDF, yaxt='n', xaxt='n', pch=NA, add=TRUE)
abline(h=DD, col='red', lwd=2, lty=3)
abline(h=0)
text(1:4, 3.5, paste("bias:", round(bias, 2)), cex=0.8, pos=4)
text(1:4, 3.3, paste("RMSE:", round(RMSE, 2)), cex=0.8, pos=4)
# Check confidence interval coverage:
colnames(result)
mean(result[, 'A_Dlow'] < DD & result[, 'A_Dupp'] > DD, na.rm=TRUE)
mean(result[, 'B_Dlow'] < DD & result[, 'B_Dupp'] > DD, na.rm=TRUE)
mean(result[, 'C_Dlow'] < DD & result[, 'C_Dupp'] > DD, na.rm=TRUE)
mean(result[, 'D_Dlow'] < DD & result[, 'D_Dupp'] > DD, na.rm=TRUE)
@mikemeredith

This comment has been minimized.

Copy link
Owner Author

commented Mar 19, 2017

Line 194 should read # be sure to run patch and getCapsRelocs, patch not path. Now corrected.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.