Nothing
############################################################################################
## package 'secr'
## fewstervarn.R
## adapted in part from RMF
## simplified from derivedSystematic 2018-12-24-27
## 2019-12-22 lacework option
## 2019-12-23 single-detector cluster option
## 2019-12-29 ncores
############################################################################################
Fewstervarn <- function (nj, xy, design, esa, detectfn, detectpar, nocc,
basenx, df, maskspacing, keep, ncores, extrapolate = TRUE)
{
# design is a list with components --
# cluster -- traps object for one cluster
# region -- sfc
# spacing -- distance between clusters)
## or --
# lacework -- logical TRUE
# region -- sfc
# spacing -- 2-vector (distance between lines, between detectors along lines)
# and possibly --
# exclude
# edgemethod
# exclmethod
design$region <- boundarytoSF(design$region)
minxy <- st_bbox(design$region)[1:2]
lacework <- !is.null(design$lacework) && design$lacework
setNumThreads(ncores) ## for multithreading in pdot
#########################################################################
## Discrete sample space (bases for systematic samples)
rx <- design$spacing[1] ## spacing between cluster centres or lacework lines
drx <- rx / basenx ## spacing between base points
bx <- seq(drx/2, rx, drx)
b <- expand.grid(x = bx, y = bx) ## assume square grid
B <- basenx^2 ## number of base points
## 2019-12-26
b <- b-rx
#########################################################################
## Boxlet centres clipped to region
boxlets <- make.mask(
type = 'polygon',
poly = design$region,
spacing = maskspacing)
## optionally clip out 'lakes'
if (!is.null(design$exclude)) {
boxlets <- subset(boxlets, !pointsInPolygon(boxlets, design$exclude))
}
#########################################################################
## Model trend, using esa(j) to represent 'size' of each sample
njxy <- data.frame (nj = nj, xy) ## dataframe with one row per cluster
if (is.null(df)) {
lgam <- gam (nj ~ s(x, y), data = njxy,
family = poisson(link = log), offset = log(esa))
}
else {
lgam <- gam (nj ~ s(x, y, fx = TRUE, k = df+1), data = njxy,
family = poisson(link = log), offset = log(esa))
}
#########################################################################
## Sum over boxlets for each basis point (starting position b)
## bxy is 2-vector of x- and y- displacement corresp to b
## pdot (X, trps, ...) gives p.(x) = 1 - (1-p(x))^nocc for animal at X
## where p(x) is aggregate over all detectors
## msk is subset of region near detectors
sumoverboxlets <- function (bxy) {
## place grid for base 'bxy' and clip to region
design$origin <- bxy + minxy
design$lacework <- NULL ## not an argument of make.lacework
makefn <- if (lacework) 'make.lacework' else 'make.systematic'
trps.b <- do.call(makefn, design)
## proceed cluster by cluster for speed, and to reduce memory demand
onecluster <- function (msk, trps) {
pxy <- covariates(msk)$pxy
a <- pdot (msk, trps, detectfn, detectpar, nocc, ncores = ncores)
c(sum(pxy * a), sum(a) * attr(msk, 'area'))
}
splittrps <- split(trps.b, clusterID(trps.b))
splitbox <- split(boxlets, clusters = splittrps)
bycluster <- mapply(onecluster, splitbox, splittrps, SIMPLIFY = TRUE)
## return 2-vector of sum(pxy.a), sum(a) over clusters
apply(bycluster, 1, sum)
}
#########################################################################
## predicted p, saved as covariate of boxlet mask
pxy <- predict.gam(lgam, newdata = boxlets, type = "response")
covariates(boxlets) <- data.frame(pxy = pxy)
if (!extrapolate) {
inner <- xy[chull(xy),]
inner <- pointsInPolygon(boxlets, inner)
innerboxlets <- subset(boxlets, inner)
outerboxlets <- subset(boxlets, !inner)
nearest <- nearesttrap(outerboxlets, innerboxlets)
covariates(boxlets)$pxy[!inner] <- covariates(boxlets)$pxy[nearest]
}
covariates(boxlets)$pxy <- covariates(boxlets)$pxy /
sum(covariates(boxlets)$pxy, na.rm=TRUE)
#########################################################################
## Qb including pdot weight
detectfn <- valid.detectfn(detectfn) ## converts character code to numeric
tmp <- apply(b, 1, sumoverboxlets) ## 2 x b matrix
Qb <- tmp[1,]
Ab <- tmp[2,]
#########################################################################
## Finally...
A <- sum(esa)
regionarea <- masksize(boxlets) ## hectares or km
N <- sum(nj) / A * regionarea ## H-T estimate
out <- 1/B * sum ((N * Qb * (1-Qb) + N^2 * Qb^2) / Ab^2) -
sum((N * Qb / Ab) / B)^2
out <- out * A^2 ## scale var(n/A) to var(n)
#########################################################################
## And optionally save lots of intermediate values
if (keep) {
attr(out, 'xy') <- xy
attr(out, 'N') <- N
attr(out, 'design') <- design
attr(out, 'b') <- b
attr(out, 'boxlets') <- boxlets # pxy is covariate of mask
attr(out, 'Qb') <- Qb
attr(out, 'Ab') <- Ab
}
out # var (n)
}
############################################################################
derivedSystematic <- function ( object, xy, design = list(),
basenx = 10, df = 9, extrapolate = TRUE, alpha = 0.05, loginterval = TRUE,
independent.esa = FALSE, keep = FALSE, ncores = NULL) {
## 'design' should comprise, e.g.,
## list(cluster = hollowsquare, spacing = 600, region = possumarea)
## or
## list(lacework = TRUE, spacing = c(500,50), region = possumarea)
warning("derivedSystematic is experimental in secr 4.6", call. = FALSE)
if (!inherits (object, 'secr') | !is.null(object$groups))
stop ("requires fitted secr model without groups")
else if (ms(object)) {
## multisession, cluster = session
if (('session' %in% object$vars) | ('Session' %in% object$vars) |
(!is.null(object$sessioncov))) {
stop ("derivedSystematic assumes detection model constant over sessions")
}
der <- derived(object, se.esa = TRUE, se.D = FALSE, bycluster = FALSE)
nj <- sapply(object$capthist, nrow) ## animals per session
maskspacing <- spacing(object$mask[[1]])
detectpar <- detectpar(object)[[1]]
nocc <- dim(object$capthist[[1]])[2]
if (is.null(design$region)) {
design$region <- attr(object$mask[[1]], 'polygon')
}
}
else {
## single session, native clusters
## except that lacework treats each detector as a cluster
if (is.null(clusterID(object$capthist))) {
warning("clusters not defined; individual detectors as clusters")
clusterID(traps(object$capthist)) <- 1:nrow(traps(object$capthist))
}
der <- derived(object, se.esa = TRUE, se.D = FALSE, bycluster = TRUE)
nj <- cluster.counts(object$capthist) ## animals per cluster
maskspacing <- spacing(object$mask)
detectpar <- detectpar(object)
nocc <- dim(object$capthist)[2]
if (is.null(design$region)) {
design$region <- attr(object$mask, 'polygon')
}
}
der <- do.call(rbind, der)
esa <- der[grepl('esa', rownames(der)),1]
se.esa <- der[grepl('esa', rownames(der)),2]
J <- length(nj) ## number of clusters
if (is.null(design$region))
stop ("region not found")
## convert matrix to sfc if not already
design$region <- boundarytoSF(design$region) ## see utility.R
design$exclude <- boundarytoSF(design$exclude) ## may be NULL
#########
## var(n)
varn <- Fewstervarn (nj, xy, design, esa, object$detectfn, detectpar, nocc,
basenx, df, maskspacing, keep, ncores)
##########
## var (A)
A <- sum(esa)
if (independent.esa)
varA <- sum(se.esa^2)
else
varA <- J^2 * sum(esa/A * se.esa^2)
## unitary A - experimental
lacework <- !is.null(design$lacework) && design$lacework
if (lacework) {
der <- derived(object, se.esa = TRUE, se.D = FALSE, bycluster = FALSE)
A <- der['esa','estimate']
varA <- der['esa','SE.estimate']^2
}
##################
## assemble output
n <- sum (nj) ## total animals
if (lacework)
D <- der['D','estimate']
else
D <- n / A
varD <- D^2 * (varn/n^2 + varA/A^2)
temp <- data.frame(row.names = c('esa','D'),
estimate = c(A,D),
SE.estimate = c(varA,varD)^0.5)
temp <- add.cl(temp, alpha, loginterval)
temp$CVn <- c(NA, varn^0.5/n)
temp$CVa <- c(NA, sqrt(varA)/A)
temp$CVD <- c(NA, varD^0.5/D)
attr(temp, 'nj') <- nj
attr(temp, 'esa') <- esa
attr(temp, 'se.esa') <- se.esa
######################
## optional attributes
if (keep) {
attr(temp, 'xy') <- xy
attr(temp, 'design') <- attr(varn, 'design')
attr(temp, 'N') <- attr(varn, 'N')
attr(temp, 'b') <- attr(varn, 'b')
attr(temp, 'boxlets') <- attr(varn, 'boxlets') # pxy is covariate of mask
attr(temp, 'Qb') <- attr(varn, 'Qb')
attr(temp, 'Ab') <- attr(varn, 'Ab')
}
temp
}
plotSystematic <- function (out, dec = 0, legend = TRUE, textcex = 0.6) {
nj <- attr(out, 'nj')
esa <- attr(out, 'esa')
b <- attr(out, 'b')
xy <- attr(out, 'xy')
N <- attr(out, 'N')
design <- attr(out, 'design')
boxlets <- attr(out, 'boxlets')
Qb <- attr(out, 'Qb')
par(mfrow = c(2,3), mar = c(2,1,2,1), cex = 0.7)
##################################################
## nj
plot(as(design$region, "Spatial")) # avoid sf::plot
text (xy[,1], xy[,2], as.character(nj), cex = textcex)
mtext(side = 3, line = -0.5, 'nj')
##################################################
## nj/esa
plot(as(design$region, "Spatial")) # avoid sf::plot
text (xy[,1], xy[,2], as.character(round(nj/esa,1)), cex = textcex)
mtext(side = 3, line = -0.5, 'nj/esa')
##################################################
## base points b Not reliable for lacework
opar <- par(mar=c(5,5,1,1))
plot(as(design$region, "Spatial")) # avoid sf::plot
bbox <- st_bbox(design$region)
points(b[,1] + bbox[1],
b[,2] + bbox[2], xpd = TRUE)
mtext(side = 3, line = -0.5, 'base points b')
par(opar)
##################################################
## grid for b[1]
design$origin <- b[1,] + bbox[1:2]
lacework <- !is.null(design$lacework) && design$lacework
design$lacework <- NULL
if (lacework) {
trps <- do.call(make.lacework, design)
}
else
trps <- do.call(make.systematic, design)
plot(as(design$region, "Spatial")) # avoid sf::plot
plot(trps, detpar = list(cex=0.7), add = TRUE)
mtext(side = 3, line = -0.5, 'grid for b[1]')
##################################################
## pxy
covariates(boxlets)$pxy4 <- covariates(boxlets)$pxy * 1e4
plot(boxlets)
plot(boxlets, cov = 'pxy4', breaks = 5, add = TRUE, legend = legend)
mtext(side = 3, line = -0.5, 'pxy * 1E4')
plot(as(design$region, "Spatial"), add = TRUE) # avoid sf::plot
##################################################
## E(n|b)
opar <- par(pty = 's', mar = c(4,4,2,1))
plot(b, type = 'n', xlab = '', ylab = '', axes = FALSE,
xlim = c(-design$spacing[1],0), ylim=c(-design$spacing[1],0))
axis(1, at = unique(b$x))
axis(2, at = unique(b$y))
text (b[,1], b[,2], round(N*Qb, dec+1), cex = 0.7)
mtext(side = 3, line = 1, 'E(n|b)')
par(opar)
##################################################
}
# plotSystematic(tmp)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.