Nothing
R/sim.capthist.R
In secr: Spatially Explicit Capture-Recapture
## package 'secr'
## sim.capthist.R
## simulate capture histories
## 2023-11-05 update availability
###############################################################################
expands <- function (x, s, default = 1) {
if (is.null(x))
rep(default, s)
else {
y <- numeric(s)
y[] <- x
y
}
}
expandsk <- function (x, s, k, default = 1) {
if (is.null(x))
matrix(default, nrow = s, ncol = k)
else
matrix(x, nrow = s, ncol = k)
}
sim.capthist <- function (
traps,
popn = list(D = 5, buffer = 100, Ndist = 'poisson'),
detectfn = 0,
detectpar = list(),
noccasions = 5,
nsessions = 1,
binomN = NULL,
exactN = NULL,
p.available = 1,
renumber = TRUE,
seed = NULL,
maxperpoly = 100,
chulltol = 0.001,
userdist = NULL,
savepopn = FALSE
)
#
# Simulate detection of known population
#
## A note on sort order of records 2009 12 01
## Simulation routines return the primary detection data in 'ski' order,
## meaning that occasion (s) changes fastest and individual (i) slowest -
## this allows efficient sequential output as new animals are detected.
## In R the order is changed to 'isk' as this is pictorially more natural.
## Secondary data (xy locations, signal strength) are generated in the order
## 'kis' (detector (k) changing fastest), because the simulation routines use -
## for (s=0; s<*ss; s++)
## for (i=0; i<*N; i++)
## for (k=0; k<*kk; k++) {
## ...
## }
## or similar loops.
## For consistency, all data are sorted to 'isk' order before returning to R.
## Secondary simulated data are re-sorted into 'ksi' order by creating the
## index 'start' within C as required for prwipolygon & prwitransect.
## The functions 'animalID', 'occasion', and 'trap' are the safest way to
## retrieve values for detections in isk order.
{
poplist <- inherits(popn,'popn') && inherits(popn,'list')
# moved here 2020-12-06
if (is.character(detectfn)) {
detectfn <- detectionfunctionnumber(detectfn)
}
if (detectfn %in% c(0:7, 14:20)) {
df0name <- if (detectfn %in% (0:7)) 'g0' else 'lambda0'
dfsname <- 'sigma'
dfzname <- if (detectfn %in% (5:6)) 'w' else 'z'
}
else {
df0name <- dfsname <- dfzname <- ""
}
##------------------------------------------------------------------------
## multi-session loop
if (poplist | (nsessions > 1)) {
if (poplist && (nsessions>1) && (length(popn) != nsessions))
stop ("incompatible use of popn list and nsessions>1")
nsessions <- ifelse (poplist, length(popn), nsessions)
if (ms(traps) & (length(traps) != nsessions))
stop ("incompatible use of traps list and nsessions>1")
output <- vector(nsessions, mode='list')
nocc <- numeric(nsessions)
nocc[] <- noccasions
##################
## set random seed
## copied from simulate.lm
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
#####################################################################
## generate population list if not provided
if (!inherits(popn,'popn')) {
popn <- replacedefaults(list(D = 5, buffer = 100,
Ndist = 'poisson'), popn)
## will fail with multiple traps
popn <- sim.popn (popn$D, core = traps, buffer = popn$buffer,
covariates = NULL, Ndist = popn$Ndist, nsessions = nsessions)
poplist <- TRUE
}
#####################################################################
## availability preliminaries
if (poplist) {
if (any(p.available != 1) && length(unique(sapply(popn, nrow))) != 1) {
warning ("incomplete availability not implemented ",
"for population lists with varying population size")
}
if(renumber)
warning ("typically use renumber = FALSE for multiple sessions")
}
if (!(length(p.available) %in% 1:2))
stop ("p.available must be vector of 1 or 2 probabilities")
if (all(p.available == 1)) {
availability <- 'complete'
}
else if (length(p.available) == 1) {
## random temporary emigration
availability <- 'random'
}
else {
## Markovian temporary emigration
availability <- 'Markov'
}
#####################################################################
## session-specific detection 2015-04-01; extended to HPX 2021-03-25
if (detectfn %in% c(0:7, 14:20)) {
if (!is.null(detectpar$individual) && detectpar$individual) {
df0 <- dfs <- dfz <- NA
}
else {
detectpar$individual <- FALSE
df0 <- expands(detectpar[[df0name]], nsessions, default = NULL)
dfs <- expands(detectpar[[dfsname]], nsessions, default = NULL)
dfz <- expands(detectpar[[dfzname]], nsessions, default = NULL)
}
}
else {
df0 <- dfs <- dfz <- NA
}
# detectfn, 9, c('g0', 'sigma'), c('b0', 'b1'))
# detectfn, 10:13, c('g0', 'sigma'), c('beta0','beta1'))
########################################################################
## loop over sessions
for (t in 1:nsessions) {
## if (R > 1) ## modified 2015-10-29
if (poplist)
temppop <- popn[[t]]
else
temppop <- popn
if (t == 1) {
if (availability == 'complete') {
available <- rep(TRUE,nrow(temppop))
}
else if (availability == 'random') {
available <- runif(nrow(temppop)) < p.available
}
else {
## Markovian temporary emigration
equilibrium.p <- (p.available[2] / (1-p.available[1]+p.available[2]))
available <- runif(nrow(temppop)) < equilibrium.p
}
}
if (availability != 'complete') {
##-------------------------------------------------------------
## update availability in preparation for next session
if (availability == 'random') {
available <- runif(nrow(temppop)) < p.available
}
else {
p.vect <- ifelse(available, p.available[1], p.available[2])
available <- runif(nrow(temppop)) < p.vect
}
#--------------------------------------------------------------
}
temppop <- subset(temppop, available)
## select session-specific traps if necessary
trps <- if (ms(traps)) traps[[t]] else traps
## session-specific detection parameters; extended to HPX 2021-03-25
if (detectfn %in% c(0:7, 14:20) && !detectpar$individual) {
detectpar[[df0name]] <- df0[t]
detectpar[[dfsname]] <- dfs[t]
detectpar[[dfzname]] <- dfz[t]
}
# detectfn, 9, c('g0', 'sigma'), c('b0', 'b1'))
# detectfn, 10:13, c('g0', 'sigma'), c('beta0','beta1'))
## recursive call for each session
## we have set the seed above
output[[t]] <- sim.capthist(trps, temppop, detectfn, detectpar,
noccasions = nocc[t], nsessions = 1,
binomN = binomN, exactN = exactN, p.available = 1,
renumber = renumber, seed = NULL, maxperpoly)
session( output[[t]] ) <- t ## added 2011-09-09
}
########################################################################
class(output) <- c('capthist', 'list')
names(output) <- 1:nsessions
output
} ## end of multi-session
else {
##################
## set random seed
## copied from simulate.lm
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
##################
if (is.null(detector(traps)))
stop ("'traps' lacks detector type")
if (length(detector(traps)) != noccasions)
detector(traps) <- rep(detector(traps)[1], noccasions)
usge <- usage(traps, noccasions)
validatepar <- function (x, xrange) {
xname <- deparse(substitute(x))
if (is.null(x))
stop ("no value for ", xname)
if (any(is.na(x)))
stop ("NA is not a valid value for ", xname)
if (any(x < xrange[1]))
warning ("value for ", xname, " is less than minimum ",
xrange[1])
if (any(x > xrange[2]))
warning ("value for ", xname, " is greater than maximum ",
xrange[2])
}
if (all(detector(traps) %in% c('signal'))) {
if (!(detectfn %in% c(10,11))) {
warning ("forcing detection function = 10 for signal detectors")
detectfn <- 10
}
}
else if (all(detector(traps) %in% c('signalnoise'))) {
if (!(detectfn %in% c(12,13))) {
warning ("forcing detection function = 12 for signalnoise detectors")
detectfn <- 12
}
}
if (any(detector(traps) %in% .localstuff$polydetectors)) {
if (!all(detectfn %in% 14:20)) {
stop("polygon and transect detectors use hazard detection functions 14:20 (HHN etc.)")
}
}
## Detection function parameters
## 0 halfnormal
## 1 hazard rate
## 2 exponential
## 3 compound halfnormal
## 4 uniform
## 5 w-exponential
## 6 annular normal
## 7 cumulative lognormal
## 8 cumulative gamma??
## 9 binary signal strength (b0 = (beta0-c)/sdS, b1 = beta1/sdS)
## 10 signal strength (signal detectors only)
## 11 signal strength with spherical spreading (signal detectors only)
## 12 signal-noise (signalnoise detectors only)
## 13 signal-noise with spherical spreading (signalnoise detectors only)
## 14 hazard halfnormal
## 15 hazard hazard-rate
## 16 hazard exponential
## 17 hazard annular normal
## 18 hazard cumulative gamma
## 19 hazard variable power
## 20 hazard pixelar 2021-03-25
if (!is.null(usage(traps))) {
if (!is.null(noccasions)) {
if (noccasions != ncol(usage(traps)))
warning ("noccasions differs from ncol(usage); using latter")
}
noccasions <- ncol(usage(traps))
}
if (detectfn %in% c(0:4)) defaults <- list(g0 = 0.2, sigma = 25, z = 1)
if (detectfn %in% c(5)) defaults <- list(g0 = 0.2, sigma = 25, w = 10)
if (detectfn %in% c(6)) defaults <- list(g0 = 0.2, sigma = 25, w = 10)
if (detectfn %in% c(7,8)) defaults <- list(g0 = 0.2, sigma = 25, z = 5)
if (detectfn %in% c(9)) defaults <- list(b0 = 1, b1=-0.1, cutval = 60,
tx = 'identity')
if (detectfn %in% c(10,11)) defaults <- list(beta0 = 90, beta1=-0.2,
sdS = 2, cutval = 60, sdM = 0, tx = 'identity')
if (detectfn %in% c(14:16,19)) defaults <- list(lambda0 = 0.2, sigma = 25, z = 1)
if (detectfn %in% c(17)) defaults <- list(lambda0 = 0.2, sigma = 25, w =10)
if (detectfn %in% c(18)) defaults <- list(lambda0 = 0.2, sigma = 25, z = 5)
# 2021-03-25
if (detectfn %in% c(20)) defaults <- list(lambda0 = 0.2, sigma = spacing(traps)/2, z = 1)
if (detectfn %in% c(12,13)) defaults <- list(beta0 = 90, beta1=-0.2,
sdS = 2, cutval = 10, muN = 40, sdN = 2, sdM = 0, tx = 'identity')
else defaults <- c(defaults, list(truncate = 1e+10, recapfactor = 1.0, bk = FALSE))
if (is.null(detectpar$individual)) {
detectpar$individual <- FALSE
}
if (is.null(binomN)) {
detectpar$binomN <- 0 ## default Poisson counts
}
else {
detectpar$binomN <- if (tolower(binomN[1]) == 'usage')
1 ## code for 'binomial size from usage' 2012-12-22
else
binomN ## as input
}
detectpar <- replacedefaults(defaults, detectpar)
# extended to HPX 2021-03-25
if (detectfn %in% c(0:7, 14:20)) {
if (detectfn %in% (0:7)) {
g0 <- expandsk(detectpar$g0, s = noccasions, k = ndetector(traps))
if (any(detector(traps) %in% .localstuff$countdetectors)) {
if (detectpar$binomN == 0)
validatepar(g0, c(0,Inf)) ## Poisson lambda
else
validatepar(g0, c(0,1)) ## Binomial p
}
df0 <- g0
}
else {
lambda0 <- expandsk(detectpar$lambda0, s = noccasions, k = ndetector(traps))
validatepar(lambda0, c(0,Inf)) ## Poisson lambda
df0 <- lambda0
}
dfs <- matrix(expands(detectpar$sigma, noccasions), ncol = 1)
dfz <- expands(ifelse(detectfn %in% c(5,6), detectpar$w, detectpar$z), noccasions)
validatepar(dfs[,1], c(1e-10,Inf))
validatepar(dfz, c(0,Inf))
}
# detectfn, 9, c('g0', 'sigma'), c('b0', 'b1'))
# detectfn, 10:13, c('g0', 'sigma'), c('beta0','beta1'))
# Acoustic detection function parameters
if (detectfn %in% c(10,11,12,13)) {
df0 <- 1 ## dummy value - not used
tx <- detectpar$tx
cutval <- detectpar$cutval
sdM <- detectpar$sdM
beta0 <- expands(detectpar$beta0, noccasions)
beta1 <- expands(detectpar$beta1, noccasions)
sdS <- expands(detectpar$sdS, noccasions)
muN <- expands(detectpar$muN, noccasions)
sdN <- expands(detectpar$sdN, noccasions)
validatepar(beta0, c(-Inf,Inf))
validatepar(beta1, c(-Inf,Inf))
validatepar(sdS, c(0,Inf))
validatepar(sdN, c(0,Inf))
}
else if (detectfn %in% c(9)) {
df0 <- expandsk(detectpar$b0, s = noccasions, k = ndetector(traps))
dfs <- expands(detectpar$b1, noccasions)
dfz <- 0
cutval <- detectpar$cutval
}
else {
cutval <- NULL
truncate <- ifelse(is.null(detectpar$truncate),
1e+10, detectpar$truncate)
validatepar(truncate, c(1e-10, Inf)) ## must be positive
}
# Allow general learned response for traps
if (any(detector(traps) %in% c('single','multi')))
{
rfl <- length(detectpar$recapfactor)
if (!(rfl %in% 1:2)) stop ("invalid recapfactor")
detectpar$recapfactor <- c(detectpar$recapfactor,1)
df0 <- cbind(df0, df0 * detectpar$recapfactor[1])
dfs <- cbind(dfs, dfs * detectpar$recapfactor[2])
}
else {
if (any(detectpar$recapfactor != 1.0))
stop("learned response available only for 'single', 'multi' detectors")
}
#-----------------------------------------------------------------------------#
if (!inherits(popn,'popn')) # generate if not provided
{
popn <- replacedefaults(list(D = 5, buffer = 100,
Ndist = 'poisson'), popn)
popn <- sim.popn (popn$D, core = traps, buffer = popn$buffer,
covariates = NULL, Ndist = popn$Ndist)
}
#-----------------------------------------------------------------------------#
if (detectpar$individual) {
# check covariates present
if(!all(c(df0name, dfsname) %in% names(covariates(popn))))
stop ("individual detection parameters missing from covariates(popn)")
}
## user-provided distances
if (is.null(userdist)) {
distmat2 <- getdistmat2(traps, popn, NULL)
}
else {
## move towards IHP/k simulations
## mask is NULL unless IHP or linear
mask <- attr(popn, 'mask')
## ASSUME MASK HAS NONEUC REQUIREMENTS 2014-10-30
## fails if any(detector %in% .localstuff$polydetectors)
distmat2 <- valid.userdist(userdist,
detector(traps),
xy1 = traps,
xy2 = popn, # animals, 2017-04-06
mask = mask)^2
}
#-----------------------------------------------------------------
stopiferror <- function(resultcode, fnname) {
if (resultcode[1] != 0)
stop ("call to '", fnname[1], "' failed")
}
##-----------------------------------------------------------------
N <- nrow(popn)
animals <- as.matrix(popn)
K <- ndetector(traps)
dettype <- detectorcode(traps, noccasions = noccasions)
simfunctionname <- paste0('trapping', detector(traps))
if (length(simfunctionname)==1 & noccasions>1)
simfunctionname <- rep(simfunctionname, noccasions)
if (length(simfunctionname)>1 & length(simfunctionname)!=noccasions)
stop("provide one detector name for each occasion")
## vector binomN
detectpar$binomN <- expandbinomN(detectpar$binomN, dettype)
if (detectfn %in% c(0:7, 14:20)) {
trappingargs <- list(
simfunctionname[1], # 1
as.double(df0), # 2
as.double(dfs), # 3
as.double(dfz), # 4
as.matrix(distmat2), # 5 assume ntraps x nanimals
as.matrix(usge), # 6
as.integer(detectfn), # 7
as.double(truncate^2), # 8
as.integer(detectpar$binomN), # 9
as.logical(detectpar$bk) # 10
)
}
## 'count' includes presence
##-----------------------------------------------------------------
## remember: output from trappingxxxx is sorted so captured animals
## precede all others
## caught[] contains capture sequence for each of N animals
if ((all(detector(traps) %in% c('single','multi','proximity','count', 'capped')))
|| detectfn==20) {
w <- array(0, dim=c(N, noccasions, K), dimnames = list(rownames(popn), 1:noccasions, 1:K))
if (length(unique(detector(traps)))==1 && detectpar$individual) {
## Individual detectpar code added 2024-03-15
sm <- all(dettype %in% c(-1,0)) # single, multi
dmat2 <- trappingargs[[5]] # traps x animals
# one individual at a time...
for (i in 1:N) {
df0 <- matrix(covariates(popn)[i, df0name], nrow = noccasions, ncol = K)
dfs <- matrix(covariates(popn)[i, dfsname], nrow = noccasions, ncol = K)
if (any(detector(traps) %in% c('single','multi')))
{
df0 <- cbind(df0, df0 * detectpar$recapfactor[1])
dfs <- cbind(dfs, dfs * detectpar$recapfactor[2])
}
trappingargs[[2]] <- df0
trappingargs[[3]] <- dfs
trappingargs[[5]] <- dmat2[,i,drop = FALSE]
temp <- do.call(simfunctionname[1], trappingargs[-1])
stopiferror(temp$resultcode, simfunctionname[1])
if (temp$n > 0) {
if (sm) {
# maximum 1 capture/animal/occasion
occ <- which(temp$value>0)
trp <- temp$value[temp$value>0]
id <- rep(i, length(trp))
w[cbind(id, occ, trp)] <- 1
}
else {
occ <- rep(1:noccasions, K)
trp <- rep(1:K, each = noccasions)
id <- rep(i, noccasions*K)
w[cbind(id, occ, trp)] <- temp$value
}
}
}
}
else if (all(dettype %in% c(-1,0))) {
# only for single- and multi-catch
temp <- do.call(simfunctionname[1], trappingargs[-1])
stopiferror(temp$resultcode, simfunctionname[1])
# if any animals caught...
if (temp$n > 0) {
# maximum 1 capture/animal/occasion
occ <- rep(1:noccasions, N)[temp$value>0]
trp <- temp$value[temp$value>0]
id <- rep(1:N, each = noccasions)[temp$value>0]
id <- which(temp$caught>0)[id] ## restore original id 2024-05-29
w[cbind(id, occ, trp)] <- 1
}
}
else {
## 2017-12-02
## external tracking of previous captures is needed to allow single-occasion calls to trappingxxx
## functions; I hope this also works with multi-occasion calls (currently disabled)
lastcapt <- rep(0, N) # NOT USED 2024-03-18
for(s in 1:noccasions) ## about 25% slower with one occasion at a time
{
sm <- dettype[s] %in% c(-1,0) # single, multi
# K1 <- if (sm) 1 else K unused?
trappingargs[[1]] <- simfunctionname[s]
trappingargs[[2]] <- as.double(df0[s,])
trappingargs[[3]] <- as.double(dfs[s,])
trappingargs[[6]] <- as.matrix(usge[,s, drop = FALSE]) ## usage on occasion s
trappingargs[[9]] <- as.integer(detectpar$binomN[s])
temp <- do.call(simfunctionname[s], trappingargs[-1])
stopiferror(temp$resultcode, simfunctionname[s])
# if any animals caught...
if (temp$n > 0) {
if (sm) {
# maximum 1 capture/animal/occasion
occ <- rep(s, sum(temp$caught>0))
trp <- temp$value[temp$caught]
id <- (1:N)[temp$caught>0]
w[cbind(id, occ, trp)] <- 1
lastcapt[temp$caught>0] <- s ## for recap
}
else {
occ <- rep(s, N*K)
trp <- rep(1:K, each = N)
id <- rep(1:N, K)
w[cbind(id, occ, trp)] <- temp$value
# lastcapt not updated because recapfactor does not apply
}
}
}
}
## drop empty histories
w <- w[apply(w,1,sum)>0,,, drop = FALSE]
class(w) <- 'capthist'
traps(w) <- traps
xy(w) <- NULL
}
##-----------------------------------------------------------------------
else if (detector(traps)[1] %in% c('polygonX','transectX')) {
if (detector(traps)[1] == 'polygonX') {
polynames <- levels(polyID(traps))
nk <- length(polynames)
K <- table(polyID(traps))
temp <- trappingpolygonX (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.integer(nk),
as.integer(K),
as.matrix(animals),
as.matrix(traps),
as.matrix(usge),
as.integer(detectfn),
as.double(truncate^2),
as.integer(detectpar$binomN)
)
}
else {
polynames <- levels(transectID(traps))
nk <- length(polynames)
if (nk>1) stop("simulation not working for multiple transects in secr 4.6")
K <- table(transectID(traps))
temp <- trappingtransectX (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.integer(nk),
as.integer(K),
as.matrix(animals),
as.matrix(traps),
as.matrix(usge),
as.integer(detectfn),
as.double(truncate^2)
)
}
if (temp$resultcode != 0)
stop ("call to '", simfunctionname, "' failed")
w <- array(0, dim=c(N, noccasions, nk),
dimnames = list(rownames(popn), 1:noccasions, polynames))
if (temp$n > 0) {
# maximum 1 capture/animal/occasion
# sorted by animal, occasion
occ <- rep(1:noccasions, temp$n)
trp <- temp$value[1:(temp$n*noccasions)]
id <- rep(1:temp$n, each = noccasions)
w[cbind(id, occ, trp)] <- 1
}
## drop empty histories
w <- w[apply(w,1,sum)>0,,, drop = FALSE]
class(w) <- 'capthist'
traps(w) <- traps
if (temp$n > 0) {
## put XY coordinates in attribute
nd <- sum(abs(w) > 0)
detectedXY <- data.frame(matrix(ncol = 2,
temp$detectedXY[1:(2*nd)]))
names(detectedXY) <- c('x','y')
xy(w) <- detectedXY
}
else
xy(w) <- NULL
}
##-----------------------------------------------------------------------
else if (detector(traps)[1] %in% c('polygon','transect','telemetry')) {
if (detector(traps)[1] == 'polygon') {
nk <- length(levels(polyID(traps)))
K <- table(polyID(traps))
temp <- trappingpolygon (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.integer(nk),
as.integer(K),
as.matrix(animals),
as.matrix(traps),
as.matrix(usge),
as.integer(detectfn),
as.double(truncate^2),
as.integer(detectpar$binomN),
as.integer(maxperpoly)
)
polynames <- levels(polyID(traps))
}
##-----------------------------------------------------------------------
if (detector(traps)[1] == 'transect') {
nk <- length(levels(transectID(traps)))
if (nk>1) stop("simulation not working for multiple transects in secr 4.6")
K <- table(transectID(traps))
temp <- trappingtransect (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.integer(nk),
as.integer(K),
as.matrix(animals),
as.matrix(traps),
as.matrix(usge),
as.integer(detectfn),
as.double(truncate^2),
as.integer(detectpar$binomN),
as.integer(maxperpoly)
)
polynames <- levels(polyID(traps))
}
##-----------------------------------------------------------------------
if (detector(traps)[1] == 'telemetry') {
nk <- 1
polynames <- '1'
maxperpoly <- 2000
if (is.null(exactN)) exactN <- 0
temp <- trappingtelemetry (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.matrix(animals),
as.integer(noccasions),
as.integer(detectfn),
as.double(truncate^2),
as.integer(detectpar$binomN),
as.integer(exactN),
as.integer(maxperpoly)
)
}
if (temp$resultcode != 0) stop("result code", temp$resultcode, "from external simulation code")
w <- array(dim=c(noccasions, nk, temp$n),
dimnames = list(1:noccasions, polynames, rownames(animals)[temp$caught]))
if (temp$n > 0) {
w[,,] <- temp$value[1:prod(dim(w))]
}
w <- aperm(w, c(3,1,2))
class(w) <- 'capthist' ## moved here so we can use animalID(w) 2017-01-11
traps(w) <- traps
if (temp$resultcode != 0) {
if (temp$resultcode == 2)
stop ("more than ", maxperpoly, " detections per animal",
" per polygon per occasion")
else
stop ("call to ",simfunctionname, " failed")
}
detectedXY <- NULL
if (temp$n > 0) {
## put XY coordinates in attribute
nd <- sum(abs(w))
detectedXY <- data.frame(matrix(ncol = 2,
temp$detectedXY[1:(2*nd)]))
names(detectedXY) <- c('x','y')
}
if (detector(traps)[1] == 'telemetry') {
ID <- animalID(w)
xyl <- split(detectedXY, ID)
telemetryxy(w) <- xyl
}
else {
xy(w) <- detectedXY
}
}
##-----------------------------------------------------------------------
else if (detector(traps)[1] %in% c('signal','signalnoise')) {
if (any(detectpar$binomN != 1))
stop ("binomN != 1 not yet implemented for signal detectors")
temp <- trappingsignal (
as.double(beta0),
as.double(beta1),
as.double(sdS),
as.double(cutval),
as.double(muN), # used only if signalnoise detector type
as.double(sdN), # used only if signalnoise detector type
as.double(sdM),
as.matrix(animals),
as.matrix(traps),
as.matrix(distmat2),
as.matrix(usge),
as.integer(detectfn)
)
# n = integer(1),
# caught = integer(N),
# signal = double(N*noccasions*K),
# noise = double(N*noccasions*K),
# value = integer(N*noccasions*K), # detected/not detected
# resultcode = integer(1) # 0,1,2
if (temp$resultcode != 0)
stop ("call to 'trappingsignal' failed with resultcode ", temp$resultcode)
w <- array(dim=c(noccasions, K, temp$n), dimnames =
list(1:noccasions,NULL,NULL))
if (temp$n>0) {
w[,,] <- temp$value[1:(temp$n * noccasions * K)]
}
w <- aperm(w, c(3,1,2))
class(w) <- 'capthist'
traps(w) <- traps
}
else stop ('Unrecognised detector type')
##-----------------------------------------------------------------------
if (!is.null(covariates(popn))) {
if (all(detector(traps) %in% c('single','multi','proximity','count','capped')))
covariates(w) <- covariates(popn)[row.names(w),, drop=F]
else
covariates(w) <- covariates(popn)[as.logical(temp$caught),, drop=F]
}
attr(w, 'cutval') <- cutval
attr(w, 'seed') <- RNGstate ## save random seed
attr(w, 'detectpar') <- detectpar
session(w) <- '1' ## dummy session values for now
## optionally save population, whether simulated or input 2014-04-27
if (savepopn)
attr(w, 'popn') <- popn
if (detector(traps)[1] %in% c('signal','signalnoise')) {
if (temp$n>0) {
signal(w) <- temp$signal[1:sum(w)]
if (detector(traps)[1] %in% c('signalnoise'))
noise(w) <- temp$noise[1:sum(w)]
}
}
if (renumber & (nrow(w)>0))
rownames(w) <- 1:nrow(w)
else {
if (!all(detector(traps) %in% c('single','multi','proximity','count','capped'))) {
rown <- rownames(popn)[temp$caught > 0]
caught <- temp$caught[temp$caught>0]
rownames(w) <- rown[order(caught)]
}
}
w
}
}
############################################################################################
sim.resight <- function (traps, popn = list(D = 5, buffer = 100, Ndist = 'poisson'), ...,
pID = 1, unmarked = TRUE, nonID = TRUE, unresolved = FALSE,
unsighted = TRUE, pmark = 0.5, Nmark = NULL, markingmask = NULL) {
## markocc is vector length = noccasions; 0 = resighting, 1 = marking
## pmark is probability an individual is marked, across the entire population,
## if all occasions are sighting occasions
############################################################################
## checks & setup
if (is.null(markocc(traps))) {
warning ("using default markocc 1 0 0 0 0")
markocc <- c(1,0,0,0,0)
}
else
markocc <- markocc(traps)
markocc(traps) <- NULL ## discard!
S <- length(markocc) ## number of sampling occasions
K <- ndetector(traps)
allsighting <- !any(markocc>0)
unres <- markocc == -1
proximityocc <- any(detector(traps) %in% 'proximity')
dettype <- detectorcode(traps, MLonly = FALSE, noccasions = S)
if (!all(dettype[markocc<1] %in% c(1,2,6,7)))
stop ("only for sightings at binary or count proximity detectors",
" or polygon or transect searches")
## special case: Nmark animals distributed according to mask covariate
distributedmarking <- FALSE
if (!is.null(markingmask) & !is.null(Nmark)) {
if (!is.null(covariates(markingmask)))
distributedmarking <- 'marking' %in% names(covariates(markingmask))
}
############################################################################
## make complete, identified capthist including all detections
## marking and sighting
dots <- list(...)
dots$traps <- traps
#################################################################
## 2017-07-26
seed <- dots$seed
dots$seed <- NULL
## set random seed
## copied from simulate.lm
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
#################################################################
dots$noccasions <- S ## override noccasions
dots$renumber <- FALSE ## to match animalID
if (!inherits(popn,'popn')) ## generate popn if not provided
{
popn <- replacedefaults(list(D = 5, buffer = 100,
Ndist = 'poisson'), popn)
if (!(distributedmarking & allsighting))
popn <- sim.popn (popn$D, core = traps, buffer = popn$buffer,
covariates = NULL, Ndist = popn$Ndist)
}
#---------------------------------------------------------------------------
if (allsighting) {
if (distributedmarking) {
## 2017-05-26
## Nunmark <- if (popn$Ndist == 'fixed') {
Ndist <- attr(popn, 'Ndist')
if (is.null(Ndist)) {
warning ("popn does not have Ndist attribute; assuming Poisson")
Ndist <- "poisson"
}
Nbuffer <- attr(popn, 'Nbuffer')
D <- attr(popn, 'D')
if (is.null(Nbuffer) & (is.null(D) | inherits(D, 'mask')))
stop ("sim.resight option not available for missing Nbuffer,D or D as mask")
Nunmark <- if (Ndist == 'fixed') {
if (!is.null(Nbuffer)) Nbuffer - Nmark
else D * masksize(markingmask) - Nmark
}
else {
rpois(1, D * masksize(markingmask)) - Nmark
}
if (Nunmark<0) {
warning ("Nunmark < 0; setting to 0")
Nunmark <- 0
}
covariates(markingmask)$pi.marking <- covariates(markingmask)$marking / sum(covariates(markingmask)$marking)
covariates(markingmask)$pi.unmarking <- 1 - covariates(markingmask)$marking
popnM <- sim.popn(D = 'pi.marking', core = markingmask, Ndist = 'fixed', Nbuffer = Nmark, model2D='IHP')
popnU <- sim.popn(D = 'pi.unmarking', core = markingmask, Ndist = 'fixed', Nbuffer = Nunmark, model2D='IHP')
}
else {
gotten <- rep(TRUE, nrow(popn))
## circumscribe marked animals
if (!is.null(markingmask)) {
gotten <- pointsInPolygon(popn, markingmask)
if (!is.null(covariates(markingmask)))
if ('marking' %in% names(covariates(markingmask)))
warning("uniform distribution over marking mask;",
" covariate 'marking' ignored")
}
if (!is.null(Nmark)) {
if (Nmark > sum(gotten)) {
warning("Nmark exceeds population size, marking all")
}
else {
gottenN <- sample.int(sum(gotten), size = Nmark, replace = FALSE)
gotten[gotten] <- (1:sum(gotten)) %in% gottenN
}
}
else {
gotten <- gotten & (runif(length(gotten)) < pmark)
}
popnM <- subset(popn, gotten) ## inside A0 and marked
popnU <- subset(popn, !gotten) ## outside A0 or not marked
}
dots$popn <- popnM
} ## end of allsighting block
else {
dots$popn <- popn
}
#---------------------------------------------------------------------------
if (ms(popn) & length(popn)==2) {
## experimental combination of 'marking' and 'sighting' populations
## these have distinct x-y, but rows corresp to the same individuals
lambda0 <- dots$detectpar$lambda0
usge <- usage(dots$traps)
dots$popn <- popn[[1]]
detector(dots$traps) <- detector(traps)[markocc==1]
markocc(dots$traps) <- NULL
dots$noccasions <- length(detector(dots$traps))
usage(dots$traps) <- usage(dots$traps)[,markocc==1]
dots$detectpar$lambda0 <- lambda0[markocc==1]
CH1 <- do.call('sim.capthist', dots)
dots$popn <- popn[[2]]
detector(dots$traps) <- detector(traps)[markocc==0]
dots$noccasions <- length(detector(dots$traps))
usage(dots$traps) <- usge[,markocc==0]
dots$detectpar$lambda0 <- lambda0[markocc==0]
CH2 <- do.call('sim.capthist', dots)
## coerce detectors purely for the sake of join()
detector(traps(CH1)) <- 'count'
detector(traps(CH2)) <- 'count'
capthist <- join(list(CH1, CH2))
detector(traps(capthist)) <- detector(traps)
}
else {
capthist <- do.call('sim.capthist', dots)
}
############################################################################
## transform simulated capthist object into resighting data
## 2021-05-19 using ksn as safe universal order regardless of polygon/signal
df <- data.frame(session = rep(1,sum(abs(capthist))),
animalID = animalID(capthist, sortorder = 'ksn'),
occasion = occasion(capthist, sortorder = 'ksn'),
trapID = trap(capthist, sortorder = 'ksn'),
stringsAsFactors = FALSE)
if (all(detector(traps) %in% c('polygon', 'transect'))) {
df <- cbind(df[,-4], xy(capthist))
df$trapID <- xyinpoly(df[,4:5], traps(capthist))
}
df <- df[order(df$animalID, df$occasion, df$trapID),]
if (nrow(df) == 0) {
dfID <- data.frame (session = 1, animalID = "NONE", occasion = 1, trapID = 1)
dfnonID <- dfID[-1,]
}
else {
# data frame of detections on marking occasions
dfmarking <- df[markocc[df$occasion]>0,,drop=FALSE]
# when was each animal first marked?
df$markingoccasion <- dfmarking$occasion[match(df$animalID, dfmarking$animalID)]
df$markingoccasion[is.na(df$markingoccasion)] <- S+1 ## never marked
if (allsighting) df$markingoccasion <- rep(0, nrow(df)) ## all marked
# separate data frames for detections of marked and unmarked animals
dfmarked <- df[df$occasion >= df$markingoccasion,, drop = FALSE]
dfunmarked <- df[df$occasion < df$markingoccasion,, drop = FALSE]
# which of the marked detections were identified?
markingoccasions <- (1:S)[markocc>0]
dfmarked$ID <- (dfmarked$occasion %in% markingoccasions) |
(runif(nrow(dfmarked)) < pID)
# separate data frames for detections of ID and nonID marked animals
dfnonID <- dfmarked[!dfmarked$ID, , drop = FALSE]
dfID <- dfmarked[dfmarked$ID, , drop = FALSE]
}
############################################################################
## Now build the new capthist object from dfID
if (all(detector(traps) %in% c('polygon', 'transect'))) {
CH <- make.capthist(dfID[,1:5, drop=FALSE], traps, fmt = 'XY', noccasions = S)
}
else {
CH <- make.capthist(dfID[,1:4, drop=FALSE], traps, fmt = 'trapID', noccasions = S)
}
markocc(traps(CH)) <- markocc
if (all(detector(traps(capthist)) %in% c('polygon', 'transect')))
trapnames <- levels(polyID(traps(CH)))
else
trapnames <- row.names(traps(CH))
## 2017-09-01
CH[,unres,] <- 0
############################################################################
## add all-zero histories for animals marked but never sighted (allsighting only)
if (allsighting & unsighted) {
covariates(CH) <- NULL # beats problem in addzeroCH 2016-11-19
nzero <- nrow(popnM) - nrow(CH)
CH <- addzeroCH(CH, nzero)
}
############################################################################
## add sightings of unmarked animals
## apply same sampling to unmarked fraction of population
if (allsighting) {
dots$popn <- popnU
capthistU <- do.call('sim.capthist', dots)
}
##------------------------------------------------------------------------------
## sightings of definitely unmarked animals
if (unmarked) {
countfn <- function(x) {
x <- x * col(x) ## x 0/1
tabulate(x[x>0], nbins = K)
}
if (allsighting) {
Tu <- as.matrix(table(
factor(trap(capthistU, sortorder = 'snk'), levels = trapnames),
factor(occasion(capthistU, sortorder = 'snk'), levels = 1:S)))
}
else {
Tu <- as.matrix(table(factor(dfunmarked$trapID, levels = trapnames),
factor(dfunmarked$occasion, levels = 1:S)))
}
## 2017-03-16, presence only if detector == proximity
Tu[, proximityocc] <- pmin (1, Tu[, proximityocc])
Tu[,markocc>0] <- 0 ## just in case...
Tu[,unres] <- 0
Tu(CH) <- Tu
}
##------------------------------------------------------------------------------
## sightings of marked animals that were not individuated
if (nonID) {
Tm <- as.matrix(table(factor(dfnonID$trapID, levels = trapnames),
factor(dfnonID$occasion, levels = 1:S)))
Tm[, proximityocc] <- pmin (1, Tm[, proximityocc])
Tm[,unres] <- 0
Tm(CH) <- Tm
}
##------------------------------------------------------------------------------
## sightings on occasions when mark status not recorded
if (unresolved) {
if (any(unres)) {
Tn <- as.matrix(table(factor(df$trapID, levels = trapnames),
factor(df$occasion, levels = 1:S)))
+ as.matrix(table(factor(dfunmarked$trapID, levels = trapnames),
factor(dfunmarked$occasion, levels = 1:S)))
Tn[,!unres] <- 0
Tn(CH) <- Tn
## optionally drop all-zero histories
if (allsighting & !unsighted) {
OK <- apply(CH,1,sum)>0
CH <- subset(CH,OK)
}
}
}
##------------------------------------------------------------------------------
## if savepopn = TRUE...
if (!is.null(attr(capthist, 'popn'))) {
if (allsighting) {
popn <- rbind(popnM, popnU)
covariates(popn) <- data.frame(marked = rep(c(TRUE, FALSE), c(nrow(popnM), nrow(popnU))))
attr(CH, 'popn') <- popn
}
else {
## 2015-12-18; this could be moved to sim.capthist
covariates(popn) <- data.frame(marked = rownames(popn) %in% dfmarked$animalID)
attr(CH, 'popn') <- popn ## NULL unless
}
}
attr(CH, 'seed') <- RNGstate ## save random seed
CH
}
############################################################################################
Try the secr package in your browser
Any scripts or data that you put into this service are public.
secr documentation built on Nov. 4, 2024, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## package 'secr'
## sim.capthist.R
## simulate capture histories
## 2023-11-05 update availability
###############################################################################
expands <- function (x, s, default = 1) {
if (is.null(x))
rep(default, s)
else {
y <- numeric(s)
y[] <- x
y
}
}
expandsk <- function (x, s, k, default = 1) {
if (is.null(x))
matrix(default, nrow = s, ncol = k)
else
matrix(x, nrow = s, ncol = k)
}
sim.capthist <- function (
traps,
popn = list(D = 5, buffer = 100, Ndist = 'poisson'),
detectfn = 0,
detectpar = list(),
noccasions = 5,
nsessions = 1,
binomN = NULL,
exactN = NULL,
p.available = 1,
renumber = TRUE,
seed = NULL,
maxperpoly = 100,
chulltol = 0.001,
userdist = NULL,
savepopn = FALSE
)
#
# Simulate detection of known population
#
## A note on sort order of records 2009 12 01
## Simulation routines return the primary detection data in 'ski' order,
## meaning that occasion (s) changes fastest and individual (i) slowest -
## this allows efficient sequential output as new animals are detected.
## In R the order is changed to 'isk' as this is pictorially more natural.
## Secondary data (xy locations, signal strength) are generated in the order
## 'kis' (detector (k) changing fastest), because the simulation routines use -
## for (s=0; s<*ss; s++)
## for (i=0; i<*N; i++)
## for (k=0; k<*kk; k++) {
## ...
## }
## or similar loops.
## For consistency, all data are sorted to 'isk' order before returning to R.
## Secondary simulated data are re-sorted into 'ksi' order by creating the
## index 'start' within C as required for prwipolygon & prwitransect.
## The functions 'animalID', 'occasion', and 'trap' are the safest way to
## retrieve values for detections in isk order.
{
poplist <- inherits(popn,'popn') && inherits(popn,'list')
# moved here 2020-12-06
if (is.character(detectfn)) {
detectfn <- detectionfunctionnumber(detectfn)
}
if (detectfn %in% c(0:7, 14:20)) {
df0name <- if (detectfn %in% (0:7)) 'g0' else 'lambda0'
dfsname <- 'sigma'
dfzname <- if (detectfn %in% (5:6)) 'w' else 'z'
}
else {
df0name <- dfsname <- dfzname <- ""
}
##------------------------------------------------------------------------
## multi-session loop
if (poplist | (nsessions > 1)) {
if (poplist && (nsessions>1) && (length(popn) != nsessions))
stop ("incompatible use of popn list and nsessions>1")
nsessions <- ifelse (poplist, length(popn), nsessions)
if (ms(traps) & (length(traps) != nsessions))
stop ("incompatible use of traps list and nsessions>1")
output <- vector(nsessions, mode='list')
nocc <- numeric(nsessions)
nocc[] <- noccasions
##################
## set random seed
## copied from simulate.lm
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
#####################################################################
## generate population list if not provided
if (!inherits(popn,'popn')) {
popn <- replacedefaults(list(D = 5, buffer = 100,
Ndist = 'poisson'), popn)
## will fail with multiple traps
popn <- sim.popn (popn$D, core = traps, buffer = popn$buffer,
covariates = NULL, Ndist = popn$Ndist, nsessions = nsessions)
poplist <- TRUE
}
#####################################################################
## availability preliminaries
if (poplist) {
if (any(p.available != 1) && length(unique(sapply(popn, nrow))) != 1) {
warning ("incomplete availability not implemented ",
"for population lists with varying population size")
}
if(renumber)
warning ("typically use renumber = FALSE for multiple sessions")
}
if (!(length(p.available) %in% 1:2))
stop ("p.available must be vector of 1 or 2 probabilities")
if (all(p.available == 1)) {
availability <- 'complete'
}
else if (length(p.available) == 1) {
## random temporary emigration
availability <- 'random'
}
else {
## Markovian temporary emigration
availability <- 'Markov'
}
#####################################################################
## session-specific detection 2015-04-01; extended to HPX 2021-03-25
if (detectfn %in% c(0:7, 14:20)) {
if (!is.null(detectpar$individual) && detectpar$individual) {
df0 <- dfs <- dfz <- NA
}
else {
detectpar$individual <- FALSE
df0 <- expands(detectpar[[df0name]], nsessions, default = NULL)
dfs <- expands(detectpar[[dfsname]], nsessions, default = NULL)
dfz <- expands(detectpar[[dfzname]], nsessions, default = NULL)
}
}
else {
df0 <- dfs <- dfz <- NA
}
# detectfn, 9, c('g0', 'sigma'), c('b0', 'b1'))
# detectfn, 10:13, c('g0', 'sigma'), c('beta0','beta1'))
########################################################################
## loop over sessions
for (t in 1:nsessions) {
## if (R > 1) ## modified 2015-10-29
if (poplist)
temppop <- popn[[t]]
else
temppop <- popn
if (t == 1) {
if (availability == 'complete') {
available <- rep(TRUE,nrow(temppop))
}
else if (availability == 'random') {
available <- runif(nrow(temppop)) < p.available
}
else {
## Markovian temporary emigration
equilibrium.p <- (p.available[2] / (1-p.available[1]+p.available[2]))
available <- runif(nrow(temppop)) < equilibrium.p
}
}
if (availability != 'complete') {
##-------------------------------------------------------------
## update availability in preparation for next session
if (availability == 'random') {
available <- runif(nrow(temppop)) < p.available
}
else {
p.vect <- ifelse(available, p.available[1], p.available[2])
available <- runif(nrow(temppop)) < p.vect
}
#--------------------------------------------------------------
}
temppop <- subset(temppop, available)
## select session-specific traps if necessary
trps <- if (ms(traps)) traps[[t]] else traps
## session-specific detection parameters; extended to HPX 2021-03-25
if (detectfn %in% c(0:7, 14:20) && !detectpar$individual) {
detectpar[[df0name]] <- df0[t]
detectpar[[dfsname]] <- dfs[t]
detectpar[[dfzname]] <- dfz[t]
}
# detectfn, 9, c('g0', 'sigma'), c('b0', 'b1'))
# detectfn, 10:13, c('g0', 'sigma'), c('beta0','beta1'))
## recursive call for each session
## we have set the seed above
output[[t]] <- sim.capthist(trps, temppop, detectfn, detectpar,
noccasions = nocc[t], nsessions = 1,
binomN = binomN, exactN = exactN, p.available = 1,
renumber = renumber, seed = NULL, maxperpoly)
session( output[[t]] ) <- t ## added 2011-09-09
}
########################################################################
class(output) <- c('capthist', 'list')
names(output) <- 1:nsessions
output
} ## end of multi-session
else {
##################
## set random seed
## copied from simulate.lm
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
##################
if (is.null(detector(traps)))
stop ("'traps' lacks detector type")
if (length(detector(traps)) != noccasions)
detector(traps) <- rep(detector(traps)[1], noccasions)
usge <- usage(traps, noccasions)
validatepar <- function (x, xrange) {
xname <- deparse(substitute(x))
if (is.null(x))
stop ("no value for ", xname)
if (any(is.na(x)))
stop ("NA is not a valid value for ", xname)
if (any(x < xrange[1]))
warning ("value for ", xname, " is less than minimum ",
xrange[1])
if (any(x > xrange[2]))
warning ("value for ", xname, " is greater than maximum ",
xrange[2])
}
if (all(detector(traps) %in% c('signal'))) {
if (!(detectfn %in% c(10,11))) {
warning ("forcing detection function = 10 for signal detectors")
detectfn <- 10
}
}
else if (all(detector(traps) %in% c('signalnoise'))) {
if (!(detectfn %in% c(12,13))) {
warning ("forcing detection function = 12 for signalnoise detectors")
detectfn <- 12
}
}
if (any(detector(traps) %in% .localstuff$polydetectors)) {
if (!all(detectfn %in% 14:20)) {
stop("polygon and transect detectors use hazard detection functions 14:20 (HHN etc.)")
}
}
## Detection function parameters
## 0 halfnormal
## 1 hazard rate
## 2 exponential
## 3 compound halfnormal
## 4 uniform
## 5 w-exponential
## 6 annular normal
## 7 cumulative lognormal
## 8 cumulative gamma??
## 9 binary signal strength (b0 = (beta0-c)/sdS, b1 = beta1/sdS)
## 10 signal strength (signal detectors only)
## 11 signal strength with spherical spreading (signal detectors only)
## 12 signal-noise (signalnoise detectors only)
## 13 signal-noise with spherical spreading (signalnoise detectors only)
## 14 hazard halfnormal
## 15 hazard hazard-rate
## 16 hazard exponential
## 17 hazard annular normal
## 18 hazard cumulative gamma
## 19 hazard variable power
## 20 hazard pixelar 2021-03-25
if (!is.null(usage(traps))) {
if (!is.null(noccasions)) {
if (noccasions != ncol(usage(traps)))
warning ("noccasions differs from ncol(usage); using latter")
}
noccasions <- ncol(usage(traps))
}
if (detectfn %in% c(0:4)) defaults <- list(g0 = 0.2, sigma = 25, z = 1)
if (detectfn %in% c(5)) defaults <- list(g0 = 0.2, sigma = 25, w = 10)
if (detectfn %in% c(6)) defaults <- list(g0 = 0.2, sigma = 25, w = 10)
if (detectfn %in% c(7,8)) defaults <- list(g0 = 0.2, sigma = 25, z = 5)
if (detectfn %in% c(9)) defaults <- list(b0 = 1, b1=-0.1, cutval = 60,
tx = 'identity')
if (detectfn %in% c(10,11)) defaults <- list(beta0 = 90, beta1=-0.2,
sdS = 2, cutval = 60, sdM = 0, tx = 'identity')
if (detectfn %in% c(14:16,19)) defaults <- list(lambda0 = 0.2, sigma = 25, z = 1)
if (detectfn %in% c(17)) defaults <- list(lambda0 = 0.2, sigma = 25, w =10)
if (detectfn %in% c(18)) defaults <- list(lambda0 = 0.2, sigma = 25, z = 5)
# 2021-03-25
if (detectfn %in% c(20)) defaults <- list(lambda0 = 0.2, sigma = spacing(traps)/2, z = 1)
if (detectfn %in% c(12,13)) defaults <- list(beta0 = 90, beta1=-0.2,
sdS = 2, cutval = 10, muN = 40, sdN = 2, sdM = 0, tx = 'identity')
else defaults <- c(defaults, list(truncate = 1e+10, recapfactor = 1.0, bk = FALSE))
if (is.null(detectpar$individual)) {
detectpar$individual <- FALSE
}
if (is.null(binomN)) {
detectpar$binomN <- 0 ## default Poisson counts
}
else {
detectpar$binomN <- if (tolower(binomN[1]) == 'usage')
1 ## code for 'binomial size from usage' 2012-12-22
else
binomN ## as input
}
detectpar <- replacedefaults(defaults, detectpar)
# extended to HPX 2021-03-25
if (detectfn %in% c(0:7, 14:20)) {
if (detectfn %in% (0:7)) {
g0 <- expandsk(detectpar$g0, s = noccasions, k = ndetector(traps))
if (any(detector(traps) %in% .localstuff$countdetectors)) {
if (detectpar$binomN == 0)
validatepar(g0, c(0,Inf)) ## Poisson lambda
else
validatepar(g0, c(0,1)) ## Binomial p
}
df0 <- g0
}
else {
lambda0 <- expandsk(detectpar$lambda0, s = noccasions, k = ndetector(traps))
validatepar(lambda0, c(0,Inf)) ## Poisson lambda
df0 <- lambda0
}
dfs <- matrix(expands(detectpar$sigma, noccasions), ncol = 1)
dfz <- expands(ifelse(detectfn %in% c(5,6), detectpar$w, detectpar$z), noccasions)
validatepar(dfs[,1], c(1e-10,Inf))
validatepar(dfz, c(0,Inf))
}
# detectfn, 9, c('g0', 'sigma'), c('b0', 'b1'))
# detectfn, 10:13, c('g0', 'sigma'), c('beta0','beta1'))
# Acoustic detection function parameters
if (detectfn %in% c(10,11,12,13)) {
df0 <- 1 ## dummy value - not used
tx <- detectpar$tx
cutval <- detectpar$cutval
sdM <- detectpar$sdM
beta0 <- expands(detectpar$beta0, noccasions)
beta1 <- expands(detectpar$beta1, noccasions)
sdS <- expands(detectpar$sdS, noccasions)
muN <- expands(detectpar$muN, noccasions)
sdN <- expands(detectpar$sdN, noccasions)
validatepar(beta0, c(-Inf,Inf))
validatepar(beta1, c(-Inf,Inf))
validatepar(sdS, c(0,Inf))
validatepar(sdN, c(0,Inf))
}
else if (detectfn %in% c(9)) {
df0 <- expandsk(detectpar$b0, s = noccasions, k = ndetector(traps))
dfs <- expands(detectpar$b1, noccasions)
dfz <- 0
cutval <- detectpar$cutval
}
else {
cutval <- NULL
truncate <- ifelse(is.null(detectpar$truncate),
1e+10, detectpar$truncate)
validatepar(truncate, c(1e-10, Inf)) ## must be positive
}
# Allow general learned response for traps
if (any(detector(traps) %in% c('single','multi')))
{
rfl <- length(detectpar$recapfactor)
if (!(rfl %in% 1:2)) stop ("invalid recapfactor")
detectpar$recapfactor <- c(detectpar$recapfactor,1)
df0 <- cbind(df0, df0 * detectpar$recapfactor[1])
dfs <- cbind(dfs, dfs * detectpar$recapfactor[2])
}
else {
if (any(detectpar$recapfactor != 1.0))
stop("learned response available only for 'single', 'multi' detectors")
}
#-----------------------------------------------------------------------------#
if (!inherits(popn,'popn')) # generate if not provided
{
popn <- replacedefaults(list(D = 5, buffer = 100,
Ndist = 'poisson'), popn)
popn <- sim.popn (popn$D, core = traps, buffer = popn$buffer,
covariates = NULL, Ndist = popn$Ndist)
}
#-----------------------------------------------------------------------------#
if (detectpar$individual) {
# check covariates present
if(!all(c(df0name, dfsname) %in% names(covariates(popn))))
stop ("individual detection parameters missing from covariates(popn)")
}
## user-provided distances
if (is.null(userdist)) {
distmat2 <- getdistmat2(traps, popn, NULL)
}
else {
## move towards IHP/k simulations
## mask is NULL unless IHP or linear
mask <- attr(popn, 'mask')
## ASSUME MASK HAS NONEUC REQUIREMENTS 2014-10-30
## fails if any(detector %in% .localstuff$polydetectors)
distmat2 <- valid.userdist(userdist,
detector(traps),
xy1 = traps,
xy2 = popn, # animals, 2017-04-06
mask = mask)^2
}
#-----------------------------------------------------------------
stopiferror <- function(resultcode, fnname) {
if (resultcode[1] != 0)
stop ("call to '", fnname[1], "' failed")
}
##-----------------------------------------------------------------
N <- nrow(popn)
animals <- as.matrix(popn)
K <- ndetector(traps)
dettype <- detectorcode(traps, noccasions = noccasions)
simfunctionname <- paste0('trapping', detector(traps))
if (length(simfunctionname)==1 & noccasions>1)
simfunctionname <- rep(simfunctionname, noccasions)
if (length(simfunctionname)>1 & length(simfunctionname)!=noccasions)
stop("provide one detector name for each occasion")
## vector binomN
detectpar$binomN <- expandbinomN(detectpar$binomN, dettype)
if (detectfn %in% c(0:7, 14:20)) {
trappingargs <- list(
simfunctionname[1], # 1
as.double(df0), # 2
as.double(dfs), # 3
as.double(dfz), # 4
as.matrix(distmat2), # 5 assume ntraps x nanimals
as.matrix(usge), # 6
as.integer(detectfn), # 7
as.double(truncate^2), # 8
as.integer(detectpar$binomN), # 9
as.logical(detectpar$bk) # 10
)
}
## 'count' includes presence
##-----------------------------------------------------------------
## remember: output from trappingxxxx is sorted so captured animals
## precede all others
## caught[] contains capture sequence for each of N animals
if ((all(detector(traps) %in% c('single','multi','proximity','count', 'capped')))
|| detectfn==20) {
w <- array(0, dim=c(N, noccasions, K), dimnames = list(rownames(popn), 1:noccasions, 1:K))
if (length(unique(detector(traps)))==1 && detectpar$individual) {
## Individual detectpar code added 2024-03-15
sm <- all(dettype %in% c(-1,0)) # single, multi
dmat2 <- trappingargs[[5]] # traps x animals
# one individual at a time...
for (i in 1:N) {
df0 <- matrix(covariates(popn)[i, df0name], nrow = noccasions, ncol = K)
dfs <- matrix(covariates(popn)[i, dfsname], nrow = noccasions, ncol = K)
if (any(detector(traps) %in% c('single','multi')))
{
df0 <- cbind(df0, df0 * detectpar$recapfactor[1])
dfs <- cbind(dfs, dfs * detectpar$recapfactor[2])
}
trappingargs[[2]] <- df0
trappingargs[[3]] <- dfs
trappingargs[[5]] <- dmat2[,i,drop = FALSE]
temp <- do.call(simfunctionname[1], trappingargs[-1])
stopiferror(temp$resultcode, simfunctionname[1])
if (temp$n > 0) {
if (sm) {
# maximum 1 capture/animal/occasion
occ <- which(temp$value>0)
trp <- temp$value[temp$value>0]
id <- rep(i, length(trp))
w[cbind(id, occ, trp)] <- 1
}
else {
occ <- rep(1:noccasions, K)
trp <- rep(1:K, each = noccasions)
id <- rep(i, noccasions*K)
w[cbind(id, occ, trp)] <- temp$value
}
}
}
}
else if (all(dettype %in% c(-1,0))) {
# only for single- and multi-catch
temp <- do.call(simfunctionname[1], trappingargs[-1])
stopiferror(temp$resultcode, simfunctionname[1])
# if any animals caught...
if (temp$n > 0) {
# maximum 1 capture/animal/occasion
occ <- rep(1:noccasions, N)[temp$value>0]
trp <- temp$value[temp$value>0]
id <- rep(1:N, each = noccasions)[temp$value>0]
id <- which(temp$caught>0)[id] ## restore original id 2024-05-29
w[cbind(id, occ, trp)] <- 1
}
}
else {
## 2017-12-02
## external tracking of previous captures is needed to allow single-occasion calls to trappingxxx
## functions; I hope this also works with multi-occasion calls (currently disabled)
lastcapt <- rep(0, N) # NOT USED 2024-03-18
for(s in 1:noccasions) ## about 25% slower with one occasion at a time
{
sm <- dettype[s] %in% c(-1,0) # single, multi
# K1 <- if (sm) 1 else K unused?
trappingargs[[1]] <- simfunctionname[s]
trappingargs[[2]] <- as.double(df0[s,])
trappingargs[[3]] <- as.double(dfs[s,])
trappingargs[[6]] <- as.matrix(usge[,s, drop = FALSE]) ## usage on occasion s
trappingargs[[9]] <- as.integer(detectpar$binomN[s])
temp <- do.call(simfunctionname[s], trappingargs[-1])
stopiferror(temp$resultcode, simfunctionname[s])
# if any animals caught...
if (temp$n > 0) {
if (sm) {
# maximum 1 capture/animal/occasion
occ <- rep(s, sum(temp$caught>0))
trp <- temp$value[temp$caught]
id <- (1:N)[temp$caught>0]
w[cbind(id, occ, trp)] <- 1
lastcapt[temp$caught>0] <- s ## for recap
}
else {
occ <- rep(s, N*K)
trp <- rep(1:K, each = N)
id <- rep(1:N, K)
w[cbind(id, occ, trp)] <- temp$value
# lastcapt not updated because recapfactor does not apply
}
}
}
}
## drop empty histories
w <- w[apply(w,1,sum)>0,,, drop = FALSE]
class(w) <- 'capthist'
traps(w) <- traps
xy(w) <- NULL
}
##-----------------------------------------------------------------------
else if (detector(traps)[1] %in% c('polygonX','transectX')) {
if (detector(traps)[1] == 'polygonX') {
polynames <- levels(polyID(traps))
nk <- length(polynames)
K <- table(polyID(traps))
temp <- trappingpolygonX (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.integer(nk),
as.integer(K),
as.matrix(animals),
as.matrix(traps),
as.matrix(usge),
as.integer(detectfn),
as.double(truncate^2),
as.integer(detectpar$binomN)
)
}
else {
polynames <- levels(transectID(traps))
nk <- length(polynames)
if (nk>1) stop("simulation not working for multiple transects in secr 4.6")
K <- table(transectID(traps))
temp <- trappingtransectX (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.integer(nk),
as.integer(K),
as.matrix(animals),
as.matrix(traps),
as.matrix(usge),
as.integer(detectfn),
as.double(truncate^2)
)
}
if (temp$resultcode != 0)
stop ("call to '", simfunctionname, "' failed")
w <- array(0, dim=c(N, noccasions, nk),
dimnames = list(rownames(popn), 1:noccasions, polynames))
if (temp$n > 0) {
# maximum 1 capture/animal/occasion
# sorted by animal, occasion
occ <- rep(1:noccasions, temp$n)
trp <- temp$value[1:(temp$n*noccasions)]
id <- rep(1:temp$n, each = noccasions)
w[cbind(id, occ, trp)] <- 1
}
## drop empty histories
w <- w[apply(w,1,sum)>0,,, drop = FALSE]
class(w) <- 'capthist'
traps(w) <- traps
if (temp$n > 0) {
## put XY coordinates in attribute
nd <- sum(abs(w) > 0)
detectedXY <- data.frame(matrix(ncol = 2,
temp$detectedXY[1:(2*nd)]))
names(detectedXY) <- c('x','y')
xy(w) <- detectedXY
}
else
xy(w) <- NULL
}
##-----------------------------------------------------------------------
else if (detector(traps)[1] %in% c('polygon','transect','telemetry')) {
if (detector(traps)[1] == 'polygon') {
nk <- length(levels(polyID(traps)))
K <- table(polyID(traps))
temp <- trappingpolygon (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.integer(nk),
as.integer(K),
as.matrix(animals),
as.matrix(traps),
as.matrix(usge),
as.integer(detectfn),
as.double(truncate^2),
as.integer(detectpar$binomN),
as.integer(maxperpoly)
)
polynames <- levels(polyID(traps))
}
##-----------------------------------------------------------------------
if (detector(traps)[1] == 'transect') {
nk <- length(levels(transectID(traps)))
if (nk>1) stop("simulation not working for multiple transects in secr 4.6")
K <- table(transectID(traps))
temp <- trappingtransect (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.integer(nk),
as.integer(K),
as.matrix(animals),
as.matrix(traps),
as.matrix(usge),
as.integer(detectfn),
as.double(truncate^2),
as.integer(detectpar$binomN),
as.integer(maxperpoly)
)
polynames <- levels(polyID(traps))
}
##-----------------------------------------------------------------------
if (detector(traps)[1] == 'telemetry') {
nk <- 1
polynames <- '1'
maxperpoly <- 2000
if (is.null(exactN)) exactN <- 0
temp <- trappingtelemetry (
as.double(df0),
as.double(dfs),
as.double(dfz),
as.matrix(animals),
as.integer(noccasions),
as.integer(detectfn),
as.double(truncate^2),
as.integer(detectpar$binomN),
as.integer(exactN),
as.integer(maxperpoly)
)
}
if (temp$resultcode != 0) stop("result code", temp$resultcode, "from external simulation code")
w <- array(dim=c(noccasions, nk, temp$n),
dimnames = list(1:noccasions, polynames, rownames(animals)[temp$caught]))
if (temp$n > 0) {
w[,,] <- temp$value[1:prod(dim(w))]
}
w <- aperm(w, c(3,1,2))
class(w) <- 'capthist' ## moved here so we can use animalID(w) 2017-01-11
traps(w) <- traps
if (temp$resultcode != 0) {
if (temp$resultcode == 2)
stop ("more than ", maxperpoly, " detections per animal",
" per polygon per occasion")
else
stop ("call to ",simfunctionname, " failed")
}
detectedXY <- NULL
if (temp$n > 0) {
## put XY coordinates in attribute
nd <- sum(abs(w))
detectedXY <- data.frame(matrix(ncol = 2,
temp$detectedXY[1:(2*nd)]))
names(detectedXY) <- c('x','y')
}
if (detector(traps)[1] == 'telemetry') {
ID <- animalID(w)
xyl <- split(detectedXY, ID)
telemetryxy(w) <- xyl
}
else {
xy(w) <- detectedXY
}
}
##-----------------------------------------------------------------------
else if (detector(traps)[1] %in% c('signal','signalnoise')) {
if (any(detectpar$binomN != 1))
stop ("binomN != 1 not yet implemented for signal detectors")
temp <- trappingsignal (
as.double(beta0),
as.double(beta1),
as.double(sdS),
as.double(cutval),
as.double(muN), # used only if signalnoise detector type
as.double(sdN), # used only if signalnoise detector type
as.double(sdM),
as.matrix(animals),
as.matrix(traps),
as.matrix(distmat2),
as.matrix(usge),
as.integer(detectfn)
)
# n = integer(1),
# caught = integer(N),
# signal = double(N*noccasions*K),
# noise = double(N*noccasions*K),
# value = integer(N*noccasions*K), # detected/not detected
# resultcode = integer(1) # 0,1,2
if (temp$resultcode != 0)
stop ("call to 'trappingsignal' failed with resultcode ", temp$resultcode)
w <- array(dim=c(noccasions, K, temp$n), dimnames =
list(1:noccasions,NULL,NULL))
if (temp$n>0) {
w[,,] <- temp$value[1:(temp$n * noccasions * K)]
}
w <- aperm(w, c(3,1,2))
class(w) <- 'capthist'
traps(w) <- traps
}
else stop ('Unrecognised detector type')
##-----------------------------------------------------------------------
if (!is.null(covariates(popn))) {
if (all(detector(traps) %in% c('single','multi','proximity','count','capped')))
covariates(w) <- covariates(popn)[row.names(w),, drop=F]
else
covariates(w) <- covariates(popn)[as.logical(temp$caught),, drop=F]
}
attr(w, 'cutval') <- cutval
attr(w, 'seed') <- RNGstate ## save random seed
attr(w, 'detectpar') <- detectpar
session(w) <- '1' ## dummy session values for now
## optionally save population, whether simulated or input 2014-04-27
if (savepopn)
attr(w, 'popn') <- popn
if (detector(traps)[1] %in% c('signal','signalnoise')) {
if (temp$n>0) {
signal(w) <- temp$signal[1:sum(w)]
if (detector(traps)[1] %in% c('signalnoise'))
noise(w) <- temp$noise[1:sum(w)]
}
}
if (renumber & (nrow(w)>0))
rownames(w) <- 1:nrow(w)
else {
if (!all(detector(traps) %in% c('single','multi','proximity','count','capped'))) {
rown <- rownames(popn)[temp$caught > 0]
caught <- temp$caught[temp$caught>0]
rownames(w) <- rown[order(caught)]
}
}
w
}
}
############################################################################################
sim.resight <- function (traps, popn = list(D = 5, buffer = 100, Ndist = 'poisson'), ...,
pID = 1, unmarked = TRUE, nonID = TRUE, unresolved = FALSE,
unsighted = TRUE, pmark = 0.5, Nmark = NULL, markingmask = NULL) {
## markocc is vector length = noccasions; 0 = resighting, 1 = marking
## pmark is probability an individual is marked, across the entire population,
## if all occasions are sighting occasions
############################################################################
## checks & setup
if (is.null(markocc(traps))) {
warning ("using default markocc 1 0 0 0 0")
markocc <- c(1,0,0,0,0)
}
else
markocc <- markocc(traps)
markocc(traps) <- NULL ## discard!
S <- length(markocc) ## number of sampling occasions
K <- ndetector(traps)
allsighting <- !any(markocc>0)
unres <- markocc == -1
proximityocc <- any(detector(traps) %in% 'proximity')
dettype <- detectorcode(traps, MLonly = FALSE, noccasions = S)
if (!all(dettype[markocc<1] %in% c(1,2,6,7)))
stop ("only for sightings at binary or count proximity detectors",
" or polygon or transect searches")
## special case: Nmark animals distributed according to mask covariate
distributedmarking <- FALSE
if (!is.null(markingmask) & !is.null(Nmark)) {
if (!is.null(covariates(markingmask)))
distributedmarking <- 'marking' %in% names(covariates(markingmask))
}
############################################################################
## make complete, identified capthist including all detections
## marking and sighting
dots <- list(...)
dots$traps <- traps
#################################################################
## 2017-07-26
seed <- dots$seed
dots$seed <- NULL
## set random seed
## copied from simulate.lm
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
#################################################################
dots$noccasions <- S ## override noccasions
dots$renumber <- FALSE ## to match animalID
if (!inherits(popn,'popn')) ## generate popn if not provided
{
popn <- replacedefaults(list(D = 5, buffer = 100,
Ndist = 'poisson'), popn)
if (!(distributedmarking & allsighting))
popn <- sim.popn (popn$D, core = traps, buffer = popn$buffer,
covariates = NULL, Ndist = popn$Ndist)
}
#---------------------------------------------------------------------------
if (allsighting) {
if (distributedmarking) {
## 2017-05-26
## Nunmark <- if (popn$Ndist == 'fixed') {
Ndist <- attr(popn, 'Ndist')
if (is.null(Ndist)) {
warning ("popn does not have Ndist attribute; assuming Poisson")
Ndist <- "poisson"
}
Nbuffer <- attr(popn, 'Nbuffer')
D <- attr(popn, 'D')
if (is.null(Nbuffer) & (is.null(D) | inherits(D, 'mask')))
stop ("sim.resight option not available for missing Nbuffer,D or D as mask")
Nunmark <- if (Ndist == 'fixed') {
if (!is.null(Nbuffer)) Nbuffer - Nmark
else D * masksize(markingmask) - Nmark
}
else {
rpois(1, D * masksize(markingmask)) - Nmark
}
if (Nunmark<0) {
warning ("Nunmark < 0; setting to 0")
Nunmark <- 0
}
covariates(markingmask)$pi.marking <- covariates(markingmask)$marking / sum(covariates(markingmask)$marking)
covariates(markingmask)$pi.unmarking <- 1 - covariates(markingmask)$marking
popnM <- sim.popn(D = 'pi.marking', core = markingmask, Ndist = 'fixed', Nbuffer = Nmark, model2D='IHP')
popnU <- sim.popn(D = 'pi.unmarking', core = markingmask, Ndist = 'fixed', Nbuffer = Nunmark, model2D='IHP')
}
else {
gotten <- rep(TRUE, nrow(popn))
## circumscribe marked animals
if (!is.null(markingmask)) {
gotten <- pointsInPolygon(popn, markingmask)
if (!is.null(covariates(markingmask)))
if ('marking' %in% names(covariates(markingmask)))
warning("uniform distribution over marking mask;",
" covariate 'marking' ignored")
}
if (!is.null(Nmark)) {
if (Nmark > sum(gotten)) {
warning("Nmark exceeds population size, marking all")
}
else {
gottenN <- sample.int(sum(gotten), size = Nmark, replace = FALSE)
gotten[gotten] <- (1:sum(gotten)) %in% gottenN
}
}
else {
gotten <- gotten & (runif(length(gotten)) < pmark)
}
popnM <- subset(popn, gotten) ## inside A0 and marked
popnU <- subset(popn, !gotten) ## outside A0 or not marked
}
dots$popn <- popnM
} ## end of allsighting block
else {
dots$popn <- popn
}
#---------------------------------------------------------------------------
if (ms(popn) & length(popn)==2) {
## experimental combination of 'marking' and 'sighting' populations
## these have distinct x-y, but rows corresp to the same individuals
lambda0 <- dots$detectpar$lambda0
usge <- usage(dots$traps)
dots$popn <- popn[[1]]
detector(dots$traps) <- detector(traps)[markocc==1]
markocc(dots$traps) <- NULL
dots$noccasions <- length(detector(dots$traps))
usage(dots$traps) <- usage(dots$traps)[,markocc==1]
dots$detectpar$lambda0 <- lambda0[markocc==1]
CH1 <- do.call('sim.capthist', dots)
dots$popn <- popn[[2]]
detector(dots$traps) <- detector(traps)[markocc==0]
dots$noccasions <- length(detector(dots$traps))
usage(dots$traps) <- usge[,markocc==0]
dots$detectpar$lambda0 <- lambda0[markocc==0]
CH2 <- do.call('sim.capthist', dots)
## coerce detectors purely for the sake of join()
detector(traps(CH1)) <- 'count'
detector(traps(CH2)) <- 'count'
capthist <- join(list(CH1, CH2))
detector(traps(capthist)) <- detector(traps)
}
else {
capthist <- do.call('sim.capthist', dots)
}
############################################################################
## transform simulated capthist object into resighting data
## 2021-05-19 using ksn as safe universal order regardless of polygon/signal
df <- data.frame(session = rep(1,sum(abs(capthist))),
animalID = animalID(capthist, sortorder = 'ksn'),
occasion = occasion(capthist, sortorder = 'ksn'),
trapID = trap(capthist, sortorder = 'ksn'),
stringsAsFactors = FALSE)
if (all(detector(traps) %in% c('polygon', 'transect'))) {
df <- cbind(df[,-4], xy(capthist))
df$trapID <- xyinpoly(df[,4:5], traps(capthist))
}
df <- df[order(df$animalID, df$occasion, df$trapID),]
if (nrow(df) == 0) {
dfID <- data.frame (session = 1, animalID = "NONE", occasion = 1, trapID = 1)
dfnonID <- dfID[-1,]
}
else {
# data frame of detections on marking occasions
dfmarking <- df[markocc[df$occasion]>0,,drop=FALSE]
# when was each animal first marked?
df$markingoccasion <- dfmarking$occasion[match(df$animalID, dfmarking$animalID)]
df$markingoccasion[is.na(df$markingoccasion)] <- S+1 ## never marked
if (allsighting) df$markingoccasion <- rep(0, nrow(df)) ## all marked
# separate data frames for detections of marked and unmarked animals
dfmarked <- df[df$occasion >= df$markingoccasion,, drop = FALSE]
dfunmarked <- df[df$occasion < df$markingoccasion,, drop = FALSE]
# which of the marked detections were identified?
markingoccasions <- (1:S)[markocc>0]
dfmarked$ID <- (dfmarked$occasion %in% markingoccasions) |
(runif(nrow(dfmarked)) < pID)
# separate data frames for detections of ID and nonID marked animals
dfnonID <- dfmarked[!dfmarked$ID, , drop = FALSE]
dfID <- dfmarked[dfmarked$ID, , drop = FALSE]
}
############################################################################
## Now build the new capthist object from dfID
if (all(detector(traps) %in% c('polygon', 'transect'))) {
CH <- make.capthist(dfID[,1:5, drop=FALSE], traps, fmt = 'XY', noccasions = S)
}
else {
CH <- make.capthist(dfID[,1:4, drop=FALSE], traps, fmt = 'trapID', noccasions = S)
}
markocc(traps(CH)) <- markocc
if (all(detector(traps(capthist)) %in% c('polygon', 'transect')))
trapnames <- levels(polyID(traps(CH)))
else
trapnames <- row.names(traps(CH))
## 2017-09-01
CH[,unres,] <- 0
############################################################################
## add all-zero histories for animals marked but never sighted (allsighting only)
if (allsighting & unsighted) {
covariates(CH) <- NULL # beats problem in addzeroCH 2016-11-19
nzero <- nrow(popnM) - nrow(CH)
CH <- addzeroCH(CH, nzero)
}
############################################################################
## add sightings of unmarked animals
## apply same sampling to unmarked fraction of population
if (allsighting) {
dots$popn <- popnU
capthistU <- do.call('sim.capthist', dots)
}
##------------------------------------------------------------------------------
## sightings of definitely unmarked animals
if (unmarked) {
countfn <- function(x) {
x <- x * col(x) ## x 0/1
tabulate(x[x>0], nbins = K)
}
if (allsighting) {
Tu <- as.matrix(table(
factor(trap(capthistU, sortorder = 'snk'), levels = trapnames),
factor(occasion(capthistU, sortorder = 'snk'), levels = 1:S)))
}
else {
Tu <- as.matrix(table(factor(dfunmarked$trapID, levels = trapnames),
factor(dfunmarked$occasion, levels = 1:S)))
}
## 2017-03-16, presence only if detector == proximity
Tu[, proximityocc] <- pmin (1, Tu[, proximityocc])
Tu[,markocc>0] <- 0 ## just in case...
Tu[,unres] <- 0
Tu(CH) <- Tu
}
##------------------------------------------------------------------------------
## sightings of marked animals that were not individuated
if (nonID) {
Tm <- as.matrix(table(factor(dfnonID$trapID, levels = trapnames),
factor(dfnonID$occasion, levels = 1:S)))
Tm[, proximityocc] <- pmin (1, Tm[, proximityocc])
Tm[,unres] <- 0
Tm(CH) <- Tm
}
##------------------------------------------------------------------------------
## sightings on occasions when mark status not recorded
if (unresolved) {
if (any(unres)) {
Tn <- as.matrix(table(factor(df$trapID, levels = trapnames),
factor(df$occasion, levels = 1:S)))
+ as.matrix(table(factor(dfunmarked$trapID, levels = trapnames),
factor(dfunmarked$occasion, levels = 1:S)))
Tn[,!unres] <- 0
Tn(CH) <- Tn
## optionally drop all-zero histories
if (allsighting & !unsighted) {
OK <- apply(CH,1,sum)>0
CH <- subset(CH,OK)
}
}
}
##------------------------------------------------------------------------------
## if savepopn = TRUE...
if (!is.null(attr(capthist, 'popn'))) {
if (allsighting) {
popn <- rbind(popnM, popnU)
covariates(popn) <- data.frame(marked = rep(c(TRUE, FALSE), c(nrow(popnM), nrow(popnU))))
attr(CH, 'popn') <- popn
}
else {
## 2015-12-18; this could be moved to sim.capthist
covariates(popn) <- data.frame(marked = rownames(popn) %in% dfmarked$animalID)
attr(CH, 'popn') <- popn ## NULL unless
}
}
attr(CH, 'seed') <- RNGstate ## save random seed
CH
}
############################################################################################
Try the secr package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.