Nothing
R/utility.R
In secr: Spatially Explicit Capture-Recapture
Defines functions
gamsetup
getknownclass
clean.mlogit
getMeanSD
fixpmix
var.in.model
expanddet
valid.model
Documented in
getMeanSD
################################################################################
## package 'secr'
## utility.R
################################################################################
## 2019-07-27 secr 4.0.0
## 2019-07-27 makelookupcpp replaces makelookup
## 2019-08-12 individualcovariates
## 2019-08-14 setNumThreads
## 2019-11-29 secr 4.1.0
## 2020-01-08 distancetotrap revised for polygon detectors
## 2020-01-26 selectCHsession
## 2020-02-21 secr 4.2.0
## 2020-05-15 stringsAsFactors function
## 2020-07-14 secr 4.3.0 distmat
## 2020-09-05 getknownclass factor bug fixed
## 2021-04-25 4.4.0
## 2021-12-16 tidy up transformations, allow arbitrary link X(), invX(), se.invX()
## 2022-01-04 4.5.0
## 2022-01-18 4.5.1
## 2022-01-23 uniformusage() function to assign all-ones usage matrix
## 2022-02-10 4.5.2 shift to sf where possible
## 2022-04-22 4.5.4
## 2022-05-31 4.5.5
## 2022-10-08 4.5.7
## 2022-11-29 4.5.8
## 2023-03-09 4.5.9
## 2023-03-10 4.5.11
## 2023-03-10 distancetotrap and nearesttrap moved to separate file
## 2023-03-10 setNumThreads moved to separate file
## 2023-05-21 4.6.0
## 2023-09-17 im2mask converts spatstat im object
################################################################################
# Global variables in namespace
#
## define a local environment for temporary variables e.g. iter
## e.g. Roger Peng https://stat.ethz.ch/pipermail/r-devel/2009-March/052883.html
.localstuff <- new.env()
## .localstuff$packageType <- ' pre-release'
.localstuff$packageType <- ''
.localstuff$validdetectors <- c('single','multi','proximity','count',
'polygonX', 'transectX', 'signal', 'polygon', 'transect',
'capped', 'null','null','null','null', 'telemetry', 'signalnoise')
.localstuff$simpledetectors <- c('single','multi','proximity','count', 'capped')
.localstuff$individualdetectors <- c('single','multi','proximity','count',
'polygonX', 'transectX', 'signal', 'signalnoise', 'polygon', 'transect',
'telemetry', 'capped')
.localstuff$pointdetectors <- c('single','multi','proximity','count',
'signal', 'signalnoise', 'unmarked','presence','capped')
.localstuff$polydetectors <- c('polygon','transect','polygonX','transectX')
.localstuff$exclusivedetectors <- c('single','multi','polygonX','transectX')
.localstuff$countdetectors <- c('count','polygon','transect','unmarked','telemetry')
.localstuff$iter <- 0 ## counter 1
.localstuff$iter2 <- 0 ## counter 2
.localstuff$detectionfunctions <-
c('halfnormal',
'hazard rate',
'exponential',
'compound halfnormal',
'uniform',
'w exponential',
'annular normal',
'cumulative lognormal',
'cumulative gamma',
'binary signal strength',
'signal strength',
'signal strength spherical',
'signal-noise',
'signal-noise spherical',
'hazard halfnormal',
'hazard hazard rate',
'hazard exponential',
'hazard annular normal',
'hazard cumulative gamma',
'hazard variable power',
'hazard pixelar')
.localstuff$DFN <- c('HN', 'HR', 'EX', 'CHN', 'UN', 'WEX', 'ANN', 'CLN', 'CG',
'BSS', 'SS', 'SSS', 'SN', 'SNS',
'HHN', 'HHR', 'HEX', 'HAN', 'HCG', 'HVP','HPX')
.localstuff$learnedresponses <- c('b', 'bk', 'B', 'k', 'Bk') ## Bk added 2020-02-26
#-------------------------------------------------------------------------------
detectionfunctionname <- function (fn) {
.localstuff$detectionfunctions[fn+1]
}
#-------------------------------------------------------------------------------
detectionfunctionnumber <- function (detname) {
dfn <- match (toupper(detname), .localstuff$DFN)
if (is.na(dfn))
dfn <- match (tolower(detname), .localstuff$detectionfunctions)
if (is.na(dfn))
stop ("unrecognised detection function ", detname)
dfn-1
}
#-------------------------------------------------------------------------------
parnames <- function (detectfn) {
parnames <- switch (detectfn+1,
c('g0','sigma'), ## 0
c('g0','sigma','z'),
c('g0','sigma'),
c('g0','sigma','z'),
c('g0','sigma'),
c('g0','sigma','w'),
c('g0','sigma','w'),
c('g0','sigma','z'),
c('g0','sigma','z'),
c('b0','b1'),
c('beta0','beta1', 'sdS'), ## include cutval?
c('beta0','beta1', 'sdS'), ## include cutval?
c('beta0','beta1', 'sdS','muN','sdN'),
c('beta0','beta1', 'sdS','muN','sdN'),
c('lambda0','sigma'),
c('lambda0','sigma','z'),
c('lambda0','sigma'),
c('lambda0','sigma','w'),
c('lambda0','sigma','z'),
c('lambda0','sigma','z'),
c('lambda0') ## 20
)
}
#-------------------------------------------------------------------------------
getdfn <- function (detectfn) {
switch (detectfn+1, HN, HR, EX, CHN, UN, WEX, ANN, CLN, CG, BSS, SS, SSS,
SN, SNS, HHN, HHR, HEX, HAN, HCG, HVP, HPX)
}
#-------------------------------------------------------------------------------
valid.detectfn <- function (detectfn, valid = c(0:3,5:19, 20)) {
# exclude 4 uniform: too numerically flakey
if (is.null(detectfn))
stop ("requires 'detectfn'")
if (is.character(detectfn))
detectfn <- detectionfunctionnumber(detectfn)
if (!(detectfn %in% valid))
stop ("invalid detection function")
detectfn
}
#-------------------------------------------------------------------------------
valid.detectpar <- function (detectpar, detectfn) {
if (is.null(detectpar) | is.null(detectfn))
stop ("requires valid 'detectpar' and 'detectfn'")
## 2013-07-19, 2013-10-22
## replace a0 with g0 or lambda0 as appropriate to detectfn
if ('a0' %in% names(detectpar)) {
aname <- if (detectfn %in% 0:8) 'g0' else 'lambda0'
lambda0 <- detectpar[['a0']] / 2 / pi / detectpar[[2]]^2 * 10000
detectpar[[aname]] <- if (detectfn %in% 0:8) 1-exp(-lambda0) else lambda0
}
if (!all(parnames(detectfn) %in% names(detectpar)))
stop ("requires 'detectpar' ", paste(parnames(detectfn), collapse=','),
" for ", detectionfunctionname(detectfn), " detectfn")
detectpar[parnames(detectfn)]
}
#-------------------------------------------------------------------------------
valid.model <- function(model, CL, detectfn, hcov, userdist, sessioncovnames) {
## 2014-08-25
if (any(sapply(model, badsmooths)))
warning("smooth term may be unsuitable for secr: does not specify k or fx where required")
}
#-------------------------------------------------------------------------------
getuserdistnames <- function (userdist) {
## return the names of any supplementary arguments of user-provided function
## for non-euclidean distance computations
if (is.function(userdist)) {
distnames <- try(userdist(), silent = TRUE)
if (!is.character(distnames))
stop("invalid userdist function - ",
"should return parameter names when called with no arguments")
distnames
}
else
character(0)
}
#-------------------------------------------------------------------------------
valid.pnames <- function (details, CL, detectfn, alltelem, sighting, nmix) {
## modelled parameters
pnames <- switch (detectfn+1,
c('g0','sigma'), # 0 halfnormal
c('g0','sigma','z'), # 1 hazard rate
c('g0','sigma'), # 2 exponential
c('g0','sigma','z'), # 3
c('g0','sigma'), # 4
c('g0','sigma','w'), # 5
c('g0','sigma','w'), # 6
c('g0','sigma','z'), # 7
c('g0','sigma','z'), # 8
c('b0','b1'), # 9
c('beta0','beta1','sdS'), # 10
c('beta0','beta1','sdS'), # 11
c('beta0','beta1','sdS'), # 12 cf parnames() in utility.R: muN, sdN?
c('beta0','beta1','sdS'), # 13 cf parnames() in utility.R: muN, sdN?
c('lambda0','sigma'), # 14 hazard halfnormal
c('lambda0','sigma','z'), # 15 hazard hazard rate
c('lambda0','sigma'), # 16 hazard exponential
c('lambda0','sigma','w'), # 17
c('lambda0','sigma','z'), # 18
c('lambda0','sigma','z'), # 19
c('lambda0','sigma')) # 20 hazard pixelar 2021-03-25
if (details$param %in% c(2,6))
pnames[1] <- 'esa'
if (details$param %in% c(3,5))
pnames[1] <- 'a0'
if (details$param %in% 4:6) {
pnames[2] <- 'sigmak'
pnames <- c(pnames, 'c')
pnames <- c(pnames, 'd')
}
if (!CL)
pnames <- c('D', pnames)
if ('noneuc' %in% getuserdistnames(details$userdist))
pnames <- c(pnames, 'noneuc')
if (sighting)
pnames <- c(pnames, 'pID')
# if (alltelem) {
# rnum <- match(c('D','lambda0','a0','esa','g0'), pnames)
# rnum[is.na(rnum)] <- 0
# pnames <- pnames[-rnum]
# }
if (nmix>1)
pnames <- c(pnames, 'pmix')
pnames
}
#-------------------------------------------------------------------------------
valid.userdist <- function (userdist, detector, xy1, xy2, mask, sessnum) {
if (is.null(userdist)) {
## default to Euclidean distance
result <- edist(xy1, xy2)
}
else {
if (any(detector %in% .localstuff$polydetectors)) {
stop ("userdist cannot be used with polygon detector types;")
}
if (is.function(userdist))
{
OK <- getuserdistnames(userdist) %in% names(covariates(mask))
if ((length(OK)>0) & !all(OK))
stop ("covariates required by userdist function not in mask : ",
paste(getuserdistnames(userdist)[!OK], collapse=','))
# 2023-02-06 selected columns 1:2 only (mask passes miscparm)
result <- do.call(userdist, c(list(xy1[,1:2], xy2[,1:2], mask)))
}
else {
if (is.character(userdist)) {
userdist <- get(userdist, pos=-1)
}
if (is.list(userdist) & !is.data.frame(userdist)) {
if (missing(sessnum) || is.na(sessnum))
stop("This use does not yet allow for session-specific userdist")
else
result <- userdist[[sessnum]]
}
else {
result <- userdist
}
}
if (!all(dim(result) == c(nrow(xy1), nrow(xy2))))
stop ("invalid distance matrix dim = ", dim(result)[1], ',', dim(result)[2])
baddist <- (!is.finite(result)) | (result<0) | is.na(result)
if (any(baddist)) {
warning ("replacing infinite, negative and NA userdist values with 1e10")
result[baddist] <- 1e10
}
}
result
}
#-------------------------------------------------------------------------------
new.param <- function (details, model, CL) {
esa <- 'esa' %in% names(model)
a0 <- 'a0' %in% names(model)
sigmak <- 'sigmak' %in% names(model)
newparam <- details$param
if (esa & !sigmak) {
newparam <- 2
}
if (a0 & !sigmak) {
newparam <- 3
}
if (sigmak) {
if (esa) {
newparam <- 6
}
else {
if (CL)
stop ("sigmak parameterization requires full model, not CL, unless also 'esa'")
newparam <- ifelse(a0, 5, 4)
}
}
if (newparam != details$param)
warning ("Using parameterization details$param = ", newparam)
newparam
}
#-------------------------------------------------------------------------------
# secr 3.0 2016-10-05
## NEEDS MULTI-SESSION FORM
detectorcode <- function (object, MLonly = TRUE, noccasions = NULL) {
## numeric detector code from a traps object
detcode <- sapply(detector(object), switch,
single = -1,
multi = 0,
proximity = 1,
count = 2,
polygonX = 3,
transectX = 4,
signal = 5,
polygon = 6,
transect = 7,
capped = 8,
unmarked = 10,
presence = 11,
signalnoise = 12,
telemetry = 13,
-2)
if (MLonly) {
detcode <- ifelse (detcode==-1, rep(0,length(detcode)), detcode)
if (any(detcode<0))
stop ("Unrecognised detector type")
}
if (!is.null(noccasions) & (length(detcode)==1))
detcode <- rep(detcode, noccasions)
detcode
}
#-------------------------------------------------------------------------------
expanddet <- function(CH) {
trps <- traps(CH)
if (is.null(trps))
return ('nonspatial')
else {
det <- detector(trps)
if (length(det)<ncol(CH))
rep(det[1], ncol(CH))
else det
}
}
#-------------------------------------------------------------------------------
ndetectpar <- function (detectfn) {
length(parnames(detectfn))
}
#-------------------------------------------------------------------------------
replacedefaults <- function (default, user) replace(default, names(user), user)
#-------------------------------------------------------------------------------
discreteN <- function (n, N) {
tN <- trunc(N)
if (N != tN) tN + sample (x = c(1,0), prob = c(N-tN, 1-(N-tN)),
replace = T, size = n)
else rep(tN,n)
}
#-------------------------------------------------------------------------------
ndetector <- function (traps) {
if (is.null(traps))
return(1)
else if (all(detector(traps) %in% .localstuff$polydetectors))
length(levels(polyID(traps)))
else
nrow(traps)
}
#-------------------------------------------------------------------------------
memo <- function (text, trace) {
## could use message(text), but does not immediately flush console
if (trace) { cat (text, '\n')
flush.console() }
}
#-------------------------------------------------------------------------------
insertdim <- function (x, dimx, dims) {
## make vector of values
## using x repeated so as to fill array
## with dim = dims and the x values occupying dimension(s) dimx
olddim <- 1:length(dims)
olddim <- c(olddim[dimx], olddim[-dimx])
temp <- array (dim=c(dims[dimx], dims[-dimx]))
tempval <- array(dim=dims[dimx])
if (length(x) > length(tempval))
tempval[] <- x[1:length(tempval)]
else
tempval[] <- x ## repeat as needed
temp[] <- tempval ## repeat as needed
if (is.factor(x))
factor(levels(x), levels=levels(x))[aperm(temp, order(olddim))] ## 2010 02 25
else
as.vector(aperm(temp, order(olddim)))
}
#-------------------------------------------------------------------------------
pad1 <- function (x, n) {
## pad x to length n with dummy (first value)
if (is.factor(x)) {
xc <- as.character(x)
xNA <- c(xc, rep(xc[1], n-length(xc)))
out <- factor(xNA, levels=levels(x))
}
else out <- c(x, rep(x[1], n-length(x)))
out
}
#-------------------------------------------------------------------------------
padarray <- function (x, dims) {
temp <- array(dim=dims)
dimx <- dim(x)
if (all(dimx>0)) {
if (length(dimx)<2 | length(dimx)>3)
stop ("invalid array")
if (length(dimx)>2) temp[1:dimx[1], 1:dimx[2], 1:dimx[3]] <- x
else temp[1:dimx[1], 1:dimx[2]] <- x
}
temp
}
#-------------------------------------------------------------------------------
## regularize a list of formulae
stdform <- function (flist) {
LHS <- function (form) {
trms <- as.character (form)
if (length(trms)==2) '' else trms[2]
}
RHS <- function (form) {
trms <- as.character (form)
## 2020-05-14 for compatibility with R 4.0
if (length(trms)==3) as.formula(paste(trms[c(1,3)], collapse = " ")) else form
}
lhs <- sapply(flist, LHS)
temp <- lapply(flist, RHS)
if (is.null(names(flist))) names(temp) <- lhs
else names(temp) <- ifelse(names(flist) == '', lhs, names(flist))
temp
}
#-------------------------------------------------------------------------------
## miscellaneous functions
invlogit <- function (y) 1/(1+exp(-y)) # plogis(y)
logit <- function (x) log(x/(1-x)) # qlogis(x), except for invalid argument
sine <- function (x) asin (x*2-1)
invsine <- function (y) (sin(y)+1) / 2
odds <- function (x) x / (1-x)
invodds <- function (y) y / (1+y)
#-------------------------------------------------------------------------------
lnbinomial <- function (x,size,prob) {
# dbinom allowing non-integer x, forcing log = TRUE
if (x <= size) {
lgamma (size+1) - lgamma (size-x+1) - lgamma (x+1) +
x * log(prob) + (size-x) * log (1-prob)
}
else {
-Inf
}
}
#-------------------------------------------------------------------------------
model.string <- function (model, userDfn) {
# 2023-04-16 Note: model should be a list
if (!is.null(userDfn)) {
if (!is.null(model$D))
model$D <- paste('~userD', userDfn('name'), sep='.')
}
temp <- paste (names(model), as.character(model), collapse=' ', sep='')
temp
}
#-------------------------------------------------------------------------------
fixed.string <- function (fixed) {
if (is.null(fixed) | length(fixed)==0) 'none'
else paste (names(fixed), as.character(fixed), collapse=', ', sep=' = ')
}
#-------------------------------------------------------------------------------
var.in.model <- function(v,m) v %in% unlist(lapply(m, all.vars))
#-------------------------------------------------------------------------------
get.nmix <- function (model, capthist, hcov) {
model$D <- NULL ## ignore density model
model$pmix <- NULL ## pmix alone cannot make this a mixture model
nmix <- 1
if (any(var.in.model('h2', model))) {
nmix <- 2
if (any(var.in.model('h3', model)))
stop ("do not combine h2 and h3")
}
if (any(var.in.model('h3', model))) {
nmix <- 3
}
if ((nmix == 1) & (!is.null(hcov))) {
if (ms(capthist))
capthist <- capthist[[1]]
if (is.factor(covariates(capthist)[,hcov]))
lev <- levels(covariates(capthist)[,hcov])
else
lev <- levels(factor(covariates(capthist)[,hcov]))
if (all(is.na(covariates(capthist)[,hcov])))
stop ("hcov missing for all individuals, but detection model invariant")
if (length(lev) < 2)
stop ("hcov covariate not found or has fewer than 2 levels")
if (length(lev) > 2)
warning ("hcov covariate has more than 2 levels; using only first two")
nmix <- 2
}
nmix
}
#-------------------------------------------------------------------------------
fixpmix <- function(x, nmix) {
## x is a list with component pmix that is a matrix (dataframe)
## with columns 'estimate' and 'se' (and possibly others)
## referring to the linear predictor of pmix (i.e. on mlogit
## scale) and rows corresponding to rows in newdata
## (i.e. arbitrary combinations of predictors, including mixture
## class h2 or h3)
####################################################
## It is necessary that newdata include all levels
## of the mixture class.
####################################################
## 2013-10-29
## assuming mixture is always last dimension...
## previously used in collate, model.average and predict.secr
## 2015-09-30 incorporated in secr.lpredictor
temp <- matrix(x$pmix[,'estimate'], ncol = nmix)
if (nmix==2) temp[,x$pmix[,'h2']] <- x$pmix[,'estimate']
if (nmix==3) temp[,x$pmix[,'h3']] <- x$pmix[,'estimate']
temp2 <- apply(temp, 1, clean.mlogit)
x$pmix[,'estimate'] <- as.numeric(t(temp2))
if (nmix==2)
x$pmix[as.numeric(x$pmix$h2)==1,'se'] <- x$pmix[as.numeric(x$pmix$h2)==2,'se']
else
x$pmix[,'se'] <- rep(NA, nrow(x$pmix)) ## don't know how
x
}
#-------------------------------------------------------------------------------
add.cl <- function (df, alpha, loginterval, lowerbound = 0) {
## add lognormal or standard Wald intervals to dataframe with columns
## 'estimate' and 'SE.estimate'
## lowerbound added 2011-07-15
z <- abs(qnorm(1-alpha/2))
if (loginterval) {
delta <- df$estimate - lowerbound
df$lcl <- delta / exp(z * sqrt(log(1 + (df$SE.estimate /
delta)^2))) + lowerbound
df$ucl <- delta * exp(z * sqrt(log(1 + (df$SE.estimate /
delta)^2))) + lowerbound
}
else {
df$lcl <- pmax(lowerbound, df$estimate - z * df$SE.estimate)
df$ucl <- df$estimate + z * df$SE.estimate
}
df
}
#-------------------------------------------------------------------------------
spatialscale <- function (object, detectfn, session = '') {
if (inherits(object, 'secr')) {
if (ms(object))
detpar <- detectpar(object)[[session]]
else
detpar <- detectpar(object)
cutval <- object$details$cutval
}
else {
detpar <- object
cutval <- object$cutval
}
if (!is.null(detpar$sigma)) detpar$sigma
else if (detectfn == 10) {
(cutval - detpar$beta0) / detpar$beta1
}
else if (detectfn == 11) {
d11 <- function(d, beta0, beta1, c) beta0 +
beta1 * (d-1) - 10 * log10(d^2) - c
interval <- c(0,10 * (cutval - detpar$beta0) / detpar$beta1)
uniroot (d11, interval, detpar$beta0, detpar$beta1, cutval)$root
}
else if (detectfn == 9) {
- 1 / detpar$b1
}
else stop ("unrecognised detectfn")
}
#-------------------------------------------------------------------------------
## logical for whether object specifies userDfn
userD <- function (object) {
if (!inherits(object, c('secr','ipsecr')))
stop ("requires fitted model")
!is.null(object$details$userDfn)
}
#-------------------------------------------------------------------------------
## mean and SD if x numeric
getMeanSD <- function(xy) {
MeanSD <- function (x) {
if (is.numeric(x))
c(mean(x, na.rm = TRUE), sd(x, na.rm = TRUE))
else
c(NA,NA)
}
as.data.frame (apply(xy, 2, MeanSD))
}
#-------------------------------------------------------------------------------
nclusters <- function (capthist) {
if (ms(capthist)) {
lapply(capthist, nclusters)
}
else {
nmash <- attr(capthist, 'n.mash')
ifelse (is.null(nmash), 1, length(nmash))
}
}
#-------------------------------------------------------------------------------
## clunky but effective re-write 2012-09-04, improved 2016-02-20, 2016-05-10
leadingzero <- function (x) {
xc <- as.character(x)
w <- max(nchar(xc))
n0 <- function(n) paste(rep('0',n), collapse='')
paste(sapply(w-nchar(xc), n0), x, sep='')
## or, 2016-01-15, 2016-02-20 BUT DOESN'T HANDLE NON-INTEGER 2016-05-10
# if (is.character(x)) x <- as.numeric(x)
# sprintf(paste("%0", w, "d", sep = ""), x)
}
#-------------------------------------------------------------------------------
## Detection functions
HN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]
g0 * exp (-r^2 / 2 / sigma^2)
}
HR <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
g0 * (1 - exp (-(r / sigma)^-z))
}
EX <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]
g0 * exp (-r / sigma)
}
UN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]
ifelse (r<=sigma, g0, 0)
}
CHN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
g0 * ( 1 - (1 - exp (-r^2 / 2 / sigma^2)) ^ z )
}
WEX <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; w <- pars[3]
ifelse(r<=w, g0, g0*exp (-(r-w) / sigma))
}
ANN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; w <- pars[3]
g0 * exp (-(r-w)^2 / 2 / sigma^2)
}
CLN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
CV2 <- (z/sigma)^2
sdlog <- log(1 + CV2)^0.5
meanlog <- log(sigma) - sdlog^2/2
g0 * plnorm(r, meanlog, sdlog, lower.tail = FALSE)
}
CG <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
g0 * pgamma(r, shape=z, scale=sigma/z, lower.tail = FALSE)
}
CN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
x <- z * (r - sigma)
g0 * (1 + (1 - exp(x)) / (1 + exp(x)))/2
}
BSS <- function (r, pars, cutval) {
b0 <- pars[1]; b1 <- pars[2]
gam <- -(b0 + b1 * r);
pnorm (gam, mean=0, sd=1, lower.tail=FALSE)
}
SS <- function (r, pars, cutval) {
beta0 <- pars[1]; beta1 <- pars[2]; sdS <- pars[3]
if (is.null(cutval))
stop ("require 'details$cutval' for signal strength plot")
mu <- beta0 + beta1 * r
1 - pnorm (q=cutval, mean=mu, sd=sdS)
}
SSS <- function (r, pars, cutval) {
beta0 <- pars[1]; beta1 <- pars[2]; sdS <- pars[3]
if (is.null(cutval))
stop ("require 'details$cutval' for signal strength plot")
## spherical so assume distance r measured from 1 m
mu <- beta0 - 10 * log ( r^2 ) / 2.302585 + beta1 * (r-1)
mu[r<1] <- beta0
1 - pnorm (q=cutval, mean=mu, sd=sdS)
}
SN <- function (r, pars, cutval) {
beta0 <- pars[1]; beta1 <- pars[2]; sdS <- pars[3];
muN <- pars[4]; sdN <- pars[5]
muS <- beta0 + beta1 * r
1 - pnorm (q=cutval, mean=muS-muN, sd=sqrt(sdS^2+sdN^2))
}
SNS <- function (r, pars, cutval) {
beta0 <- pars[1]; beta1 <- pars[2]; sdS <- pars[3];
muN <- pars[4]; sdN <- pars[5]
## spherical so assume distance r measured from 1 m
muS <- beta0 - 10 * log ( r^2 ) / 2.302585 + beta1 * (r-1)
muS[r<1] <- beta0
1 - pnorm (q=cutval, mean=muS-muN, sd=sqrt(sdS^2+sdN^2))
}
HHN <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]
1 - exp(-lambda0 * exp (-r^2 / 2 / sigma^2))
}
HHR <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
1 - exp(-lambda0 * ( 1 - exp (-(r / sigma)^-z)))
}
HEX <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]
1 - exp(-lambda0 * exp (-r / sigma))
}
HAN <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]; w <- pars[3]
1 - exp(-lambda0 * exp (-(r-w)^2 / 2 / sigma^2))
}
HCG <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
lambda0 * pgamma(r, shape=z, scale=sigma/z, lower.tail = FALSE)
}
HVP <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
1 - exp(-lambda0 * exp(-(r/sigma)^z))
}
HPX <- function (r, pars, cutval) {
g0 <- 1-exp(-pars[1])
radius <- pars[2]
ifelse (r<=radius, g0, 0) # circular, not square! crude approx
}
#-------------------------------------------------------------------------------
gradient <- function (pars, fun, eps=0.001, ...)
## quick & dirty 2009 09 14
## used by plot.secr for delta method limits
{
est <- pars
g <- pars
for (i in 1:length(est))
{
temp <- est[i]
if (temp != 0.0) delta <- eps * abs(temp)
else delta <- eps
est[i] <- temp - delta
fminus <- fun (est, ...)
est[i] <- temp + delta
fplus <- fun (est, ...)
g[i] <- (fplus - fminus) / (2.0 * delta)
est[i] <- temp;
}
g
}
#-------------------------------------------------------------------------------
# transformation tidy up 2021-12-16
# arbitrary link function specified with functions X, invX, se.invX
transform <- function (x, link) {
switch (link,
identity = x,
i1000 = x * 1000,
log = log(x),
neglog = log(-x),
logit = logit(x),
odds = odds(x),
sin = sine(x),
do.call(link, list(x))
)
}
#-------------------------------------------------------------------------------
# used only in model.average, modelAverage
se.transform <- function (real, sereal, link) {
switch (link,
identity = sereal,
i1000 = sereal / 1000,
log = log((sereal/real)^2 + 1)^0.5,
neglog = log((sereal/-real)^2 + 1)^0.5,
logit = sereal / real / (1 - real),
sin = NA,
do.call(paste0('se.',link), list(real, sereal) )
)
}
#-------------------------------------------------------------------------------
untransform <- function (beta, link) {
switch (link,
identity = beta,
i1000 = beta / 1000,
log = exp(beta),
neglog = -exp(beta),
logit = invlogit(beta),
odds = invodds(beta),
sin = invsine(beta),
do.call(paste0('inv',link), list(beta))
)
}
#-------------------------------------------------------------------------------
se.untransform <- function (beta, sebeta, link) {
# Approximate translation of SE to untransformed scale
# Delta method cf Lebreton et al 1992 p 77
switch (link,
identity = sebeta,
i1000 = sebeta / 1000,
log = exp(beta) * sqrt(exp(sebeta^2)-1),
neglog = exp(beta) * sqrt(exp(sebeta^2)-1),
logit = invlogit(beta) * (1-invlogit(beta)) * sebeta,
sin = NA, ####!!!!
do.call(paste0('se.inv', link), list(beta=beta, sebeta=sebeta))
)
}
#-------------------------------------------------------------------------------
# vectorized transformations
Xtransform <- function (real, linkfn, varnames) {
mapply(transform, real, linkfn[varnames])
}
se.Xtransform <- function (real, sereal, linkfn, varnames) {
mapply(se.transform, real, sereal, linkfn[varnames])
}
Xuntransform <- function (beta, linkfn, varnames) {
mapply(untransform, beta, linkfn[varnames])
}
se.Xuntransform <- function (beta, sebeta, linkfn, varnames)
{
if (length(beta)!=length(sebeta))
stop ("'beta' and 'sebeta' do not match")
if (!all(varnames %in% names(linkfn)))
stop ("'linkfn' component missing for at least one real variable")
mapply(se.untransform, beta, sebeta, linkfn[varnames])
}
#-------------------------------------------------------------------------------
mlogit.untransform <- function (beta, latentmodel) {
if (!missing(latentmodel)) {
for (i in unique(latentmodel))
beta[latentmodel==i] <- mlogit.untransform(beta[latentmodel==i])
beta
}
else {
## beta should include values for all classes (mixture components)
nmix <- length(beta)
if (sum(is.na(beta)) != 1) {
## require NA for a single reference class
rep(NA, length(beta))
}
else {
nonreference <- !is.na(beta) # not reference class
b <- beta[nonreference]
pmix <- numeric(nmix)
pmix[nonreference] <- exp(b) / (1+sum(exp(b)))
pmix[!nonreference] <- 1 - sum(pmix[nonreference])
pmix
}
}
}
#-------------------------------------------------------------------------------
clean.mlogit <- function(x) {
## 2014-08-19 for robustness...
if (is.na(x[2])) x[2] <- 1-x[1]
x[1] <- NA ## assumed reference class
logit(mlogit.untransform(x))
}
#-------------------------------------------------------------------------------
mlogit <- function (x) {
## return the mlogit of an unscaled vector of positive values
## 2013-04-14
logit(x/sum(x))
}
## End of miscellaneous functions
#-------------------------------------------------------------------------------
group.levels <- function (capthist, groups, sep='.') {
# 2016-06-05 use also for trap strata
if (inherits(capthist, 'list')) {
temp <- lapply(capthist, group.levels, groups, sep)
sort(unique(unlist(temp))) ## vector of global levels
}
else {
if (is.null(groups)) 0
else {
if (!all(groups %in% names(covariates(capthist))))
stop ("one or more grouping variables is missing ",
"from covariates")
temp <- as.data.frame(covariates(capthist)[,groups])
# omit null combinations, sort as with default of factor levels
sort(levels(interaction(temp, drop=T, sep=sep)))
}
}
}
#-------------------------------------------------------------------------------
h.levels <- function (capthist, hcov, nmix) {
## determine the first nmix levels of a factor individual covariate
if (is.null(hcov))
as.character(1:nmix)
else {
if (ms(capthist)) {
## take first session as we can assume factor covariates have same levels in
## all sessions
capthist <- capthist[[1]]
}
hcov <- covariates(capthist)[,hcov]
if (!is.factor(hcov)) {
warning("hcov was coerced to a factor")
hcov <- factor(hcov)
}
levels(hcov)[1:nmix]
}
}
#-------------------------------------------------------------------------------
n.occasion <- function (capthist) {
## return the number of sampling occasions for each session in capthist
if (inherits(capthist, 'list')) {
sapply(capthist, n.occasion)
}
else {
ncol(capthist)
}
}
#-------------------------------------------------------------------------------
group.factor <- function (capthist, groups, sep='.')
## convert a set of grouping factors to a single factor (g)
## levels common to all sessions
{
if (inherits(capthist, 'list')) {
temp <- lapply(capthist, group.factor, groups) ## recursive call
grouplevels <- group.levels(capthist, groups)
if (length(grouplevels)<2)
temp
else
# list; force shared factor levels on each component
lapply (temp, factor, levels=grouplevels)
}
else {
if (is.null(groups) | (length(groups)==0) )
return (factor(rep(1, nrow(capthist)), levels = 1)) # added levels 2017-04-18
temp <- as.data.frame(covariates(capthist)[,groups])
if (ncol(temp) != length(groups))
stop ("one or more grouping variables is missing from ",
"covariates(capthist)")
temp <- interaction(temp, drop=T, sep=sep) # omit null combinations
temp
}
}
#-------------------------------------------------------------------------------
getgrpnum <- function (capthist, groups) {
if (is.null(groups))
rep(1, nrow(capthist))
else
match(group.factor(capthist, groups), group.levels(capthist, groups))
}
#-------------------------------------------------------------------------------
make.lookup <- function (tempmat) {
## should add something to protect make.lookup from bad data...
nrw <- nrow(tempmat)
ncl <- ncol(tempmat)
nam <- colnames(tempmat)
df <- is.data.frame(tempmat)
if (df) {
lev <- lapply(tempmat, levels)
tempmat[] <- sapply(tempmat, as.numeric)
tempmat <- as.matrix(tempmat)
}
dimnames(tempmat) <- NULL
temp <- makelookupcpp(tempmat)
lookup <- temp$lookup
colnames(lookup) <- nam
if (df) {
lookup <- as.data.frame(lookup)
## restore factors
for (i in 1: length(lev))
if (!is.null(lev[[i]]))
lookup[,i] <- factor(lev[[i]][lookup[,i]], levels = lev[[i]])
}
list (lookup=lookup, index=temp$index)
}
#-------------------------------------------------------------------------------
## Return an integer vector of class membership defined by a categorical
## individual covariate in a capthist object. Individuals of unknown
## class (including those with class exceeding nmix) are coded 1,
## others as (class number + 1). When no mixture is specified (nmix == 1)
## all are coded as unknown.
## knownclass 1 'unknown'
## knownclass 2 'latent class 1'
## knownclass 3 'latent class 2'
getknownclass <- function(capthist, nmix, hcov) {
if (ms(capthist)) {
lapply(capthist, getknownclass, nmix = nmix, hcov = hcov)
}
else {
if ((nmix>1) & (!is.null(hcov))) {
## 2020-09-05 use as.factor() instead of factor() to coerce
## (if already factor, coercing with factor() loses old levels)
var <- as.factor(covariates(capthist)[,hcov])
tmp <- as.numeric(var) + 1
tmp[is.na(tmp) | (tmp>(nmix+1))] <- 1
attr(tmp,'levels') <- levels(factor(covariates(capthist)
[,hcov]))[1:nmix]
tmp
}
else
rep(1,nrow(capthist))
}
}
#-------------------------------------------------------------------------------
getnmix <- function (details) {
if (is.null(details$nmix))
1
else
details$nmix
}
#-------------------------------------------------------------------------------
## expand beta parameter vector using template of 'fixed beta'
## fixed beta fb input is missing (NA) for estimated beta parameters
fullbeta <- function (beta, fb) {
if (!is.null(fb)) {
fb[is.na(fb)] <- beta ## partial beta (varying only)
beta <- fb ## complete beta
}
beta
}
#-------------------------------------------------------------------------------
## inflate a convex outline along all radii by linear factor 'rmult'
inflate <- function (xy, rmult = 1) {
xy <- as.matrix(xy)
centre <- apply(xy, 2, mean)
xy <- sweep(xy, MARGIN = 2, STATS = centre, FUN = '-')
r <- apply(xy, 1, function(z) sqrt(sum(z^2)))
theta <- atan2 (xy[,2], xy[,1])
r <- r * rmult
xy <- cbind(r * cos(theta), r * sin(theta))
sweep(xy, MARGIN = 2, STATS = centre, FUN = '+')
}
#-------------------------------------------------------------------------------
## moved from pdot.R 2013-11-09
## scalar 2016-10-14
getbinomN <- function (binomN, detectr) {
if (any(detectr %in% .localstuff$countdetectors)) {
if (is.null(binomN))
return(0)
else if (binomN == 'usage')
return(1)
else
return(binomN)
}
else
return(1)
}
#-------------------------------------------------------------------------------
## convert telemetryxy attribute of a combined dataset into a standalone capthist
## TO DO: option of telemetry or polygon output
xy2CH <- function (CH, inflation = 1e-8) {
if (ms(CH)) {
out <- lapply(CH, xy2CH, inflation)
class(out) <- c('capthist', 'list')
out
}
else {
xylist <- telemetryxy(CH)
if (is.null(xylist))
stop ("requires 'telemetryxy' attribute")
n <- length(xylist)
neach <- sapply(xylist, nrow)
allxy <- do.call(rbind, xylist)
trps <- allxy[chull(allxy),]
trps <- rbind(trps, trps[1,,drop=F])
trps <- inflate(trps, 1 + inflation) ## see also telemetry.R
trps <- as.data.frame(trps)
dimnames(trps) <- list(1:nrow(trps), c('x','y'))
class(trps) <- c("traps","data.frame")
detector(trps) <- "polygon"
polyID(trps) <- factor(rep(1,nrow(trps)))
rown <- rep(names(xylist), neach)
newCH <- array(neach, dim = c(n, 1, 1))
attr(newCH, "detectedXY") <- allxy
if (!is.null(covariates(CH))) {
rowlookup <- match(names(xylist), rownames(CH))
covariates(newCH) <- covariates(CH)[rowlookup,, drop=FALSE]
}
class(newCH) <- "capthist"
traps(newCH) <- trps
newCH
}
}
#-------------------------------------------------------------------------------
inflatechull <- function (poly, r, ntheta = 60) {
theta <- (2*pi) * (1:ntheta) / ntheta
## add supernumerary vertices
temp <- data.frame(x = apply(expand.grid(poly$x, r * cos(theta)),1,sum),
y = apply(expand.grid(poly$y, r * sin(theta)),1,sum))
hull <- chull(temp)
temp[c(hull,hull[1]), ]
}
#-------------------------------------------------------------------------------
maskarea <- function (mask, sessnum = 1) {
if (!ms(mask)) nrow(mask) * attr(mask,'area')
else nrow(mask[[sessnum]]) * attr(mask[[sessnum]],'area')
}
#-------------------------------------------------------------------------------
masklength <- function (mask, sessnum = 1) {
if (!ms(mask)) nrow(mask) * attr(mask,'spacing')/1000
else nrow(mask[[sessnum]]) * attr(mask[[sessnum]],'spacing')/1000
}
#-------------------------------------------------------------------------------
masksize <- function (mask, sessnum = 1) {
if (inherits(mask, 'linearmask'))
masklength(mask, sessnum)
else
maskarea(mask, sessnum)
}
#-------------------------------------------------------------------------------
complete.beta <- function (object) {
fb <- object$details$fixedbeta
# modified 2022-04-02 for consistency with ipsecr
beta <- if (inherits(object, 'secr')) object$fit$par else object$beta
if (!is.null(fb)) {
nbeta <- length(fb)
fb[is.na(fb)] <- beta
beta <- fb
}
beta
}
#-------------------------------------------------------------------------------
complete.beta.vcv <- function (object) {
fb <- object$details$fixedbeta
if (!is.null(fb)) {
nbeta <- length(fb)
beta.vcv <- matrix(NA, nrow = nbeta, ncol = nbeta)
beta.vcv[is.na(fb[row(beta.vcv)]) & is.na(fb[col(beta.vcv)])] <- object$beta.vcv
}
else {
beta.vcv <- object$beta.vcv
}
beta.vcv
}
#-------------------------------------------------------------------------------
smooths <- function (formula) {
## which terms in formula are smooths?
## returns logical vector
labels <- attr(terms(formula), "term.labels")
if (length(labels) > 0)
sapply(labels, function (x) any(sapply(c("s\\(", "te\\(", "poly\\("), grepl, x)))
else
logical(0)
}
#-------------------------------------------------------------------------------
polys <- function (formula) {
## which terms in formula are orthogonal polynomials?
## returns logical vector
labels <- attr(terms(formula), "term.labels")
if (length(labels) > 0)
sapply(labels, grepl, pattern = "poly\\(")
else
logical(0)
}
#-------------------------------------------------------------------------------
badsmooths <- function (formula) {
## does smooth specification conform to secr requirements?
## returns TRUE/FALSE
labels <- attr(terms(formula), "term.labels")
if (length(labels) > 0) {
smoothterms <- sapply(labels, function (x)
any(sapply(c("s\\(", "te\\("), grepl, x)))
labels <- labels[smoothterms]
any(sapply(labels, function(x)
grepl("s\\(", x) & !grepl("k =", x))) |
any(sapply(labels, function(x)
grepl("te\\(", x) & (!grepl("fx = TRUE", x) | !grepl("k =", x))))
}
else
FALSE
}
#-------------------------------------------------------------------------------
gamsetup <- function(formula, data, ...) {
## use 'session' column as dummy LHS so gam does not gag
## (cf secrgam:::make.density.design.matrix)
## session is always present in detection data, must be added for D
if (is.null(data$session)) data$session <- rep(1,nrow(data))
formula <- update.formula(formula, session ~ .)
setup <- gam(formula, data = data, fit = FALSE, ...)
colnames(setup$X) <- setup$term.names
setup
}
#-------------------------------------------------------------------------------
general.model.matrix <- function (formula, data, gamsmth = NULL,
contrasts = NULL, ...) {
## A function to compute the design matrix for the model in
## 'formula' given the data in 'data'. This is merely the result
## of model.matrix() unless 'formula' includes smooth terms -- s()
## or te() as described in mgcv ?formula.gam.
## If smooth terms are present then the matrix may be based on a
## previous gam setup (provided in the argument 'gamsmth') or
## computed de novo with gam(..., fit = FALSE)
## note 2014-08-24
## orthogonal polynomials e.g. poly(x,2) are handled by model.matrix,
## but the information needed for prediction at new data is not
## saved by secr.fit, so predict.secr generally fails with message
## "'degree' must be less than number of unique points"
## head(eval(parse(text = attr(terms(~ poly(x,y, degree=2)),
## 'term.labels')[1]), env=possummask))
## 2014-08-24, 2014-09-09, 2017-11-30
## 2019-10-12 drop row names
## 2021-12-09 f optional argument
dots <- list(...)
if (any(polys(formula)))
stop ("orthogonal polynomials are temporarily blocked") ## 2014-09-12
if (any(smooths(formula))) {
if (is.null(gamsmth)) {
## setup knots etc from scratch
mat <- gamsetup(formula, data, ...)$X
}
else {
## fool predict.gam into generating the necessary
## predictor matrix from previous setup
class (gamsmth) <- 'gam'
gamsmth$coefficients <- rep(NA, ncol(gamsmth$X))
mat <- mgcv::predict.gam(gamsmth, newdata = data, type = 'lpmatrix')
colnames(mat) <- colnames(gamsmth$X)
}
}
else {
## model.matrix(formula, data, ...)
mat <- model.matrix(formula, data = data, contrasts.arg = contrasts)
}
rownames (mat) <- NULL
mat
}
#-------------------------------------------------------------------------------
makerealparameters <- function (design, beta, parindx, link, fixed) {
modelfn <- function(i) {
## linear predictor for real parameter i
Yp <- design$designMatrices[[i]] %*% beta[parindx[[i]]]
if (names(link)[i] == 'pmix') {
## 2013-04-14 index of class groups (pmix sum to 1.0 within latentmodel)
cols <- dimnames(design$designMatrices[[i]])[[2]]
h2 <- grep('.h2', cols, fixed=T)
h3 <- grep('.h3', cols, fixed=T)
h2c <- grep(':h2', cols, fixed=T)
h3c <- grep(':h3', cols, fixed=T)
h.cols <- c(h2,h3,h2c,h3c)
tmp <- design$designMatrices[[i]][,-h.cols, drop = FALSE]
tmph <- design$designMatrices[[i]][,h.cols, drop = FALSE]
## 2018-02-23 why as.numeric()?
latentmodel <- as.numeric(factor(apply(tmp,1,paste, collapse='')))
refclass <- apply(tmph,1,sum) == 0
Yp[refclass] <- NA
Yp <- mlogit.untransform(Yp, latentmodel)
Yp[design$parameterTable[,i]]
}
else {
Yp <- untransform(Yp, link[[i]])
Yp[design$parameterTable[,i]] ## replicate as required
}
}
## construct matrix of detection parameters
nrealpar <- length(design$designMatrices)
parindx$D <- NULL ## detection parameters only
link$D <- NULL ## detection parameters only
parindx$noneuc <- NULL ## detection parameters only
link$noneuc <- NULL ## detection parameters only
detectionparameters <- names(link)
fixed.dp <- fixed[detectionparameters[detectionparameters %in% names(fixed)]]
if (length(fixed.dp)>0)
for (a in names(fixed.dp)) ## bug fixed by adding this line 2011-09-28
link[[a]] <- NULL
if (length(link) != nrealpar)
stop ("number of links does not match design matrices")
if (nrealpar == 0) {
return(matrix(unlist(fixed.dp),nrow = 1))
}
temp <- sapply (1:nrealpar, modelfn)
if (nrow(design$parameterTable)==1) temp <- t(temp)
nrw <- nrow(temp)
## make new matrix and insert columns in right place
temp2 <- as.data.frame(matrix(nrow = nrw, ncol = length(detectionparameters)))
names(temp2) <- detectionparameters
temp2[ , names(design$designMatrices)] <- temp ## modelled
if (!is.null(fixed.dp) & length(fixed.dp)>0)
temp2[ , names(fixed.dp)] <- sapply(fixed.dp, rep, nrw) ## fixed
as.matrix(temp2)
}
#-------------------------------------------------------------------------------
secr.lpredictor <- function (formula, newdata, indx, beta, field, beta.vcv=NULL,
smoothsetup = NULL, contrasts = NULL, Dfn = NULL) {
## form linear predictor for a single 'real' parameter
## smoothsetup should be provided whenever newdata differs from
## data used to fit model and the model includes smooths from gam
vars <- all.vars(formula)
OK <- vars %in% names(newdata)
if (any(!OK)) {
missingvars <- paste(vars[!OK], collapse = ', ')
if (sum(!OK) == 1)
stop ("model covariate ", missingvars, " not found in 'newdata'")
else
stop ("model covariates ", missingvars, " not found in 'newdata'")
}
newdata <- as.data.frame(newdata)
lpred <- matrix(ncol = 2, nrow = nrow(newdata), dimnames = list(NULL,c('estimate','se')))
if (!is.null(Dfn) && field == 'D') {
warning("secr.lpredictor is not ready for D as function - do not use estimates")
nsess <- length(unique(newdata$session))
Yp <- Dfn(newdata[,vars[1]], beta = beta[indx], dimD = c(nrow(newdata)/nsess,1,nsess))
mat <- as.matrix(newdata[,vars[1], drop = FALSE])
}
else {
mat <- general.model.matrix(formula, data = newdata, gamsmth = smoothsetup,
contrasts = contrasts)
if (nrow(mat) < nrow(newdata))
warning ("missing values in predictors?")
nmix <- 1
if (field=='pmix') {
## drop pmix beta0 column from design matrix (always zero)
mat <- mat[,-1,drop=FALSE]
if ('h2' %in% names(newdata)) nmix <- 2
if ('h3' %in% names(newdata)) nmix <- 3
mixfield <- c('h2','h3')[nmix-1]
}
###############################
Yp <- mat %*% beta[indx]
###############################
## A latent model comprises one row for each latent class.
## Back transformation of pmix in mlogit.untransform() requires all rows of
## each latent model. That function splits vector Yp by latent model.
if (field == 'pmix') {
nonh <- newdata[, names(newdata) != mixfield, drop = FALSE]
latentmodel <- factor(apply(nonh, 1, paste, collapse = ''))
refclass <- as.numeric(newdata[, mixfield]) == 1
Yp[refclass] <- NA
Yp <- mlogit.untransform(Yp, latentmodel)
Yp <- logit(Yp) # return to logit scale for later untransform!
if (nmix==2) {
h2.1 <- as.numeric(newdata$h2)==1
h2.2 <- as.numeric(newdata$h2)==2
}
}
}
lpred[,1] <- Yp
if (is.null(beta.vcv) || (any(is.na(beta[indx])))) return ( cbind(newdata,lpred) )
else {
if (is.null(Dfn) || field != 'D') {
vcv <- beta.vcv[indx,indx, drop = FALSE]
vcv[is.na(vcv)] <- 0
nrw <- nrow(mat)
vcv <- apply(expand.grid(1:nrw, 1:nrw), 1, function(ij)
mat[ij[1],, drop=F] %*% vcv %*% t(mat[ij[2],, drop=F]))
vcv <- matrix (vcv, nrow = nrw)
if (field=='pmix') {
if (nmix==2)
vcv[h2.1,h2.1] <- vcv[h2.2,h2.2]
else
vcv[,] <- NA
}
lpred[,2] <- diag(vcv)^0.5
}
else {
vcv <- NULL
}
temp <- cbind(newdata,lpred)
attr(temp, 'vcv') <- vcv
return(temp)
}
}
#-------------------------------------------------------------------------------
edist <- function (xy1, xy2) {
nr <- nrow(xy1)
nc <- nrow(xy2)
x1 <- matrix(xy1[,1], nr, nc)
x2 <- matrix(xy2[,1], nr, nc, byrow=T)
y1 <- matrix(xy1[,2], nr, nc)
y2 <- matrix(xy2[,2], nr, nc, byrow=T)
sqrt((x1-x2)^2 + (y1-y2)^2)
}
#-------------------------------------------------------------------------------
## least cost paths from mask including barriers to movement
## use edist for equivalent Euclidean distances
nedist <- function (xy1, xy2, mask, inf = Inf, ...) {
newargs <- list(...)
if (missing(mask)) mask <- xy2
noneuc <- covariates(mask)$noneuc
if (is.null(noneuc)) noneuc <- rep(1, nrow(mask))
defaultargs <- list(transitionFunction = mean, directions = 16)
args <- replace(defaultargs, names(newargs), newargs)
args$x <- raster(mask, values = noneuc)
if (requireNamespace('gdistance', quietly = TRUE)) { ## 2015-01-23
tr <- do.call(gdistance::transition, args)
tr <- gdistance::geoCorrection(tr, type = "c", multpl = FALSE)
out <- gdistance::costDistance(tr, as.matrix(xy1), as.matrix(xy2))
}
else stop ("package gdistance is required for nedist")
if (is.finite(inf)) out[!is.finite(out)] <- inf
out
}
#-------------------------------------------------------------------------------
getcellsize <- function (mask) {
if (inherits(mask, 'linearmask'))
cell <- attr(mask, 'spacing') / 1000 ## per km
else
cell <- attr(mask, 'area') ## per ha
if (is.null(cell))
stop ("mask lacks valid cell size (area or spacing)")
cell
}
#-------------------------------------------------------------------------------
## intercept and fix certain models with bad defaults
updatemodel <- function (model, detectfn, detectfns, oldvar, newvar, warn = FALSE) {
if (detectfn %in% detectfns) {
for (i in 1:length(oldvar)) {
if (oldvar[i] %in% names(model)) {
names(model)[names(model) == oldvar[i]] <- newvar[i]
if (warn)
warning ("replacing ", oldvar[i], " by ", newvar[i],
" in model for detectfn ", detectfn)
}
}
}
model
}
#-------------------------------------------------------------------------------
## Manually remove some mask points
# simplified 2022-02-03
deleteMaskPoints <- function (mask, onebyone = TRUE, add = FALSE, poly = NULL,
poly.habitat = FALSE, ...) {
## interface does not work properly in RStudio
if (ms(mask)) { ## a list of mask objects
if (inherits(poly, 'list') & (!is.data.frame(poly)))
stop ("lists of polygons not implemented in 'make.mask'")
temp <- lapply (mask, deleteMaskPoints, onebyone = onebyone, add = add,
poly = poly, poly.habitat = poly.habitat, ...)
class (temp) <- c('mask', 'list')
temp
}
else {
plot(mask, add = add, ...)
if (!is.null(poly)) {
if (poly.habitat)
pointstodrop <- (1:nrow(mask))[!pointsInPolygon(mask, poly)]
else
pointstodrop <- (1:nrow(mask))[pointsInPolygon(mask, poly)]
}
else if (onebyone) {
cat ('Click to select points; right-click to stop\n')
flush.console()
xy <- locator(type = 'p', pch=1, col='red')
pointstodrop <- if (length(xy$x)==0)
numeric(0)
else
nearesttrap(xy, mask)
}
else {
cat ('Click to select polygon vertices; right-click to stop\n')
flush.console()
xy <- locator(type = 'l', col='red')
xy <- as.data.frame(xy)
xy <- rbind(xy, xy[1,])
if (poly.habitat)
pointstodrop <- (1:nrow(mask))[!pointsInPolygon(mask, xy)]
else
pointstodrop <- (1:nrow(mask))[pointsInPolygon(mask, xy)]
}
npts <- length(pointstodrop)
if (npts>0) {
points(mask[pointstodrop,], pch = 16, col = 'red')
if(.Platform$OS.type == "windows") {
pl <- if (npts>1) 's' else ''
msg <- paste ('Delete ', npts, ' red point',pl, '?', sep='')
response <- utils::winDialog(type = "okcancel", msg)
} else {
response <- 'OK'
}
if (response == 'OK') {
mask <- subset(mask, -pointstodrop)
if (npts==1)
message("1 point deleted")
else
message(npts, " points deleted")
}
else
message ("point(s) not deleted")
}
else
message ("no points to delete")
plot(mask, col='green')
mask
}
}
#-------------------------------------------------------------------------------
nparameters <- function (object) {
Npar <- max(unlist(object$parindx))
Npar <- Npar + length(object$details$miscparm)
## allow for fixed beta parameters
if (!is.null(object$details$fixedbeta))
Npar <- Npar - sum(!is.na(object$details$fixedbeta))
Npar
}
#-------------------------------------------------------------------------------
mapbeta <- function (parindx0, parindx1, beta0, betaindex)
## Extend beta vector from simple model (beta0) to a more complex (i.e. general)
## model, inserting neutral values (zero) as required.
## For each real parameter, a 1:1 match is assumed between
## beta values until all beta values from the simpler model are
## used up. THIS ASSUMPTION MAY NOT BE JUSTIFIED.
## betaindex is a user-controlled alternative.
{
## list of zeroed vectors, one per real parameter
beta1 <- lapply(parindx1, function (x) {x[]<-0; x})
if (!is.null(betaindex)) {
beta1 <- unlist(beta1)
if (sum(betaindex>0) != length(beta0))
stop ("invalid 'betaindex'")
beta1[betaindex] <- beta0
beta1
}
else {
## indx is within-parameter rather than absolute index
## for each _original_ real parameter
indx <- lapply(parindx0, function(x) x-x[1]+1)
## for (j in 1:length(beta1))
## improved replace by name2015-11-17
for (j in names(beta1)) {
if (j %in% names(beta0))
beta1[[j]][indx[[j]]] <- beta0[parindx0[[j]]]
}
unlist(beta1)
}
}
#-------------------------------------------------------------------------------
xyinpoly <- function (xy, trps) {
ptinside <- function (i,k) {
## is point i inside poly k?
polyxy <- as.matrix(lxy[[k]])
polyxy <- rbind(polyxy, polyxy[1,]) ## close 2014-08-28
nr <- nrow(polyxy)
temp <- insidecpp(unlist(xy[i,]), 0, nr-1, as.matrix(polyxy))
}
lxy <- split (trps, polyID(trps))
firstinside <- function (i) {
frstk <- 0
for (k in 1:length(lxy)) {
if (ptinside(i,k)) {
frstk <- k
break
}
}
frstk
}
sapply(1:nrow(xy), firstinside)
}
#-------------------------------------------------------------------------------
addzerodf <- function (df, oldCH, sess) {
## add dummy detection records to dataframe for 'all-zero' case
## that arises in sighting-only mark-resight with known marks
allzero <- apply(oldCH,1,sum)==0
naz <- sum(allzero)
if (naz > 0) {
df0 <- expand.grid(
newID = rownames(oldCH)[allzero],
newocc = NA,
newtrap = trap(oldCH)[1],
alive = TRUE,
sess = sess,
stringsAsFactors = FALSE)
df$x <- NULL; df$y <- NULL ## 2021-04-08
df <- rbind(df,df0)
if (!is.null(xy(oldCH))) {
df$x <- c(xy(oldCH)$x, rep(NA, naz))
df$y <- c(xy(oldCH)$y, rep(NA, naz))
}
if (!is.null(signal(oldCH))) {
df$signal <- c(signal(oldCH), rep(NA, naz))
}
}
df
}
#-------------------------------------------------------------------------------
## including pre-marked animals never sighted
## cov is optional dataframe of covariates
addzeroCH <- function (CH, nzero, cov = NULL) {
if (nzero == 0)
return(CH)
else {
nc <- nrow(CH)
chdim <- dim(CH)
chdim[1] <- nzero
extra <- array(0, dim=chdim)
dimnames(extra) <- c(list(paste('Z', 1:nzero, sep='')), dimnames(CH)[2:3])
CH2 <- abind(CH, extra, along = 1)
class(CH2) <- 'capthist'
traps(CH2) <- traps(CH)
xy(CH2) <- xy(CH) ## order is not affected by adding zero histories
if (!is.null(covariates(CH)) & (nrow(CH)>0)) {
if (is.null(cov)) {
cov <- covariates(CH)[rep(1,nzero),]
cov[,] <- NA ## covariates are unknown
}
covariates(CH2) <- rbind(covariates(CH), cov[1:nzero,])
}
## ... and other essential attributes?
CH2
}
}
#-------------------------------------------------------------------------------
expandbinomN <- function (binomN, detectorcodes) {
# assumes detectorcodes is a vector of length = noccasions
binomN <- ifelse (detectorcodes %in% c(2,6,7), binomN, 1)
if (any(is.na(binomN))) stop ("NA value in binomN")
binomN
}
#-------------------------------------------------------------------------------
check3D <- function (object) {
if (ms(object)) {
out <- lapply(object, check3D)
class(out) <- class(object)
out
}
else {
if (is.matrix(object)) {
warning("secr >= 3.0 requires 3-D capthist; using updateCH() to convert")
updateCH(object)
}
else {
object
}
}
}
#-------------------------------------------------------------------------------
updateCH <- function(object) {
if (!inherits(object, 'capthist'))
stop ("requires capthist object")
# following replaces this old code 2020-08-29
# reduce(object, dropunused = FALSE)
if (ms(object)) {
out <- lapply(object, updateCH)
class (out) <- c("capthist", "list")
out
}
else {
if (length(dim(object)) == 3) {
return(object)
}
else {
K <- ndetector(traps(object))
ch <- array(0, dim = c(dim(object), K), dimnames =
list(rownames(object), colnames(object), 1:K))
OK <- as.logical(object!=0)
animal <- row(object)[OK]
occ <- col(object)[OK]
detn <- object[OK]
ch[cbind(animal, occ, detn)] <- 1
traps(ch) <- traps(object)
class (ch) <- "capthist"
session(ch) <- session(object)
ch
}
}
}
#-------------------------------------------------------------------------------
newstr <-function (strings) {
## compress a character vector
## use run length encoding function
rl <- rle(strings)
st <- rl$values
le <- paste0(' (',as.character(rl$lengths), ')')
le[le==' (1)'] <- ''
paste(paste0(st, le), collapse = ', ')
}
# newstr(c("single", rep("proximity",4)))
#-------------------------------------------------------------------------------
outsidemask <- function(CH, mask, threshold = spacing(mask) / sqrt(2)) {
xylist <- telemetryxy(CH)
dfun <- function(xy) {
centres <- matrix(apply(xy, 2, mean), ncol = 2)
distancetotrap(centres, mask)
}
sapply(xylist, dfun) > threshold
}
#-------------------------------------------------------------------------------
shareFactorLevels <- function (object, columns = NULL, stringsAsFactors = TRUE) {
## stringsAsFactors added 2020-05-16
if (ms(object)) {
if (!is.null(covariates(object))) {
df <- do.call(rbind, covariates(object))
if (is.null(columns)) {
columns <- 1:ncol(df)
}
if (stringsAsFactors) {
df[,columns] <- stringsAsFactors(df[,columns, drop = FALSE])
}
for (i in columns) {
if (is.factor(df[,i])) {
levelsi <- levels(df[,i])
for (sess in 1:length(object)) {
covariates(object[[sess]])[,i] <-
factor(covariates(object[[sess]])[,i],
levels = levelsi)
}
}
}
}
}
else {
# modified 2021-04-27 to apply to covariates, not object itself
if (!is.null(covariates(object))) {
if (stringsAsFactors) {
df <- covariates(object)
if (is.null(columns)) {
columns <- 1:ncol(df)
}
df[,columns] <- stringsAsFactors(df[,columns, drop = FALSE])
covariates(object) <- df
}
}
}
object
}
#-------------------------------------------------------------------------------
allzero <- function (object) {
if (!inherits(object, 'capthist'))
stop ("requires 'capthist' object")
if (ms(object)) {
lapply(object, allzero)
}
else {
telemocc <- detector(traps(object))=='telemetry'
apply(object[,!telemocc,,drop=FALSE],1,sum)==0
}
}
#-------------------------------------------------------------------------------
primarysessions <- function(intervals) {
primarysession <- cumsum(c(0,intervals))
match(primarysession, unique(primarysession))
}
#-------------------------------------------------------------------------------
secondarysessions <- function(intervals) {
primary <- primarysessions(intervals)
unname(unlist(sapply(table(primary), seq_len)))
}
#-------------------------------------------------------------------------------
boundarytoSF <- function (poly) {
if (is.null(poly)) {
NULL
}
else if(inherits(poly, c('sf','sfc'))) {
poly <- st_geometry(poly) # extract sfc if not already sfc
geomtype <- st_geometry_type(poly, by_geometry = FALSE)
if (geomtype == 'GEOMETRY') { # 2023-06-02
geomtype <- st_geometry_type(poly, by_geometry = TRUE)
}
if (!all(geomtype %in% c("POLYGON", "MULTIPOLYGON"))) {
stop ("poly sf/sfc should be of type POLYGON or MULTIPOLYGON")
}
poly
}
else if (inherits(poly, 'SpatialPolygons')) { # also SPDF?
st_as_sfc(poly)
}
else if (inherits(poly, 'SpatVector')) {
st_as_sfc(as(poly,"Spatial"))
}
else if (inherits(poly, c('matrix', 'data.frame'))) {
## input is 2-column matrix for a single polygon
poly <- matrix(unlist(poly), ncol = 2)
poly <- rbind (poly, poly[1,]) ## force closure of polygon
st_sfc(st_polygon(list(poly)))
}
else stop (class(poly), " not valid input to boundarytoSF")
}
#-------------------------------------------------------------------------------
pointsInPolygon <- function (xy, poly, logical = TRUE) {
# xy is 2-column matrix or data.frame of coordinates
if (inherits(poly, 'mask')) {
if (ms(poly))
stop ("multi-session masks not supported")
sp <- spacing(poly)
minx <- min(poly$x, na.rm = TRUE)
miny <- min(poly$y, na.rm = TRUE)
mask <- sweep(poly, MARGIN = 2, FUN = '+', STATS = c(-minx, -miny))
mask <- round(mask/sp) + 1
xy <- matrix(unlist(xy), ncol = 2) ## in case dataframe
xy <- sweep(xy, MARGIN = 2, FUN = '+', STATS = c(-minx, -miny))
xy <- round(xy/sp) + 1
xy[xy<=0] <- NA
xy[,1][xy[,1]>max(mask$x, na.rm = TRUE)] <- NA
xy[,2][xy[,2]>max(mask$y, na.rm = TRUE)] <- NA
maskmatrix <- matrix(0, ncol = max(mask$y, na.rm = TRUE), nrow = max(mask$x, na.rm = TRUE))
maskmatrix[as.matrix(mask)] <- 1:nrow(mask)
inside <- maskmatrix[as.matrix(xy)]
inside[is.na(inside)] <- 0
if (logical)
inside <- inside > 0
inside
}
else {
poly <- boundarytoSF(poly)
if (inherits(poly, c('sf','sfc'))) {
xy <- st_as_sf(data.frame(xy), coords = 1:2)
st_crs(xy) <- st_crs(poly)
apply(st_within(xy, poly, sparse = FALSE), 1, any)
}
else {
stop ("unknown input to pointsInPolygon")
}
}
}
#-------------------------------------------------------------------------------
## return indices of first occasion and detector for which PIAx is non-zero
firstsk <- function (PIAx) {
## PIAx dim n,s,k
wh <- function(d2) {
match(TRUE, d2>0)
}
apply(PIAx,1,wh)
}
#-------------------------------------------------------------------------------
maskboolean <- function (ch, mask, threshold) {
if (ms(ch)) {
if (!ms(mask)) stop ("masklookup: multisession ch requires multisession mask")
outlist <- mapply(maskboolean, ch, mask, MoreArgs = list(threshold = threshold), SIMPLIFY = FALSE)
outlist
}
else {
id <- animalID(ch, names = FALSE, sortorder = 'snk')
tr <- trap(ch, names = FALSE, sortorder = 'snk')
trps <- traps(ch)
m <- nrow(mask)
if (!is.null(threshold) && all(detector(trps) %in% .localstuff$pointdetectors)) {
df <- data.frame(id = id, x = trps$x[tr], y = trps$y[tr])
x <- tapply(df$x, df$id, mean, na.rm=T)
y <- tapply(df$y, df$id, mean, na.rm=T)
xy <- data.frame(x=x,y=y)
d2 <- edist2cpp(as.matrix(xy), as.matrix(mask))
out <- (d2 <= threshold^2)
}
else {
## NULL option
out <- matrix(TRUE, nrow = nrow(ch), ncol = m)
}
out
}
}
#-------------------------------------------------------------------------------
telemcode <- function(object, ...) {
if (inherits(object, 'traps') && !ms(object))
switch (telemetrytype(object), none = 0,
independent = 1, dependent = 2, concurrent = 3, 0)
else
NA
}
#-------------------------------------------------------------------------------
uniquerownames <- function (capthist) {
if (!ms(capthist)) {
return(capthist)
}
else {
last <- 0
for (i in 1:length(capthist)) {
nr <- nrow(capthist[[i]])
if (nr > 0) {
rownames(capthist[[i]]) <- last + (1:nr)
last <- last+nr
}
}
capthist
}
}
#-------------------------------------------------------------------------------
selectCHsession <- function(capthist, sessnum) {
if (ms(capthist))
capthist[[sessnum]]
else
capthist
}
#-------------------------------------------------------------------------------
stringsAsFactors <- function (DF) {
# convert any character columns of a data.frame (or list) to factor
if (is.list(DF) && length(DF)>0) { ## bug fix 2020-08-14
chr <- sapply(DF, is.character)
DF[chr] <- lapply(DF[chr], as.factor)
}
DF
}
#-------------------------------------------------------------------------------
# see also getuserdist in loglikhelperfn.R
# 2021-03-30 moved from preparedata.R
# 2021-03-30 integrate HPX
getdistmat2 <- function (traps, mask, userdist, HPX = FALSE) {
## Static distance matrix
if (is.function(userdist)) {
NULL ## compute dynamically later
}
else {
if (HPX) {
if (any(detector(traps) %in% .localstuff$polydetectors)) {
trps <- split(traps, polyID(traps))
inside <- t(sapply(trps, pointsInPolygon, xy = mask))
d2 <- 1-inside # 0 inside, 1 outside
d2[d2>0] <- 1e10 # 0 inside, 1e10 outside
d2
}
else {
# maximum of squared distance in x- or y- directions
xydist2cpp(as.matrix(traps), as.matrix(mask))
}
}
else {
if (any(detector(traps) %in% .localstuff$polydetectors)) {
## do not use result if detector is one of
## polygonX, polygon, transectX, transect, OR telemetry?
matrix(0, nrow = nrow(traps), ncol = nrow(mask))
}
else {
# Euclidean distance
## edist(as.matrix(traps), as.matrix(mask))^2
edist2cpp(as.matrix(traps), as.matrix(mask))
}
}
}
}
#-------------------------------------------------------------------------------
## function to assign all-ones usage matrix
uniformusage <- function(object, noccasions) {
if (inherits(object, 'capthist')) {
if (ms(object)) {
for (r in 1:length(object)) {
ndet <- dim(object[[r]])[3]
noccasions <- dim(object[[r]])[2]
usage(traps(object[[r]])) <- matrix(1, ndet, noccasions)
}
}
else {
ndet <- dim(object)[3]
noccasions <- dim(object)[2]
usage(traps(object)) <- matrix(1, ndet, noccasions)
}
}
else if (inherits(object, 'traps')) {
if (missing(noccasions)) {
stop ('noccasions should be specified for traps input')
}
if (ms(object)) {
for (r in 1:length(object)) {
ndet <- ndetector(object[[r]])
usage(object[[r]]) <- matrix(1, ndet, noccasions)
}
}
else {
ndet <- ndetector(object)
usage(object) <- matrix(1, ndet, noccasions)
}
}
object
}
#-------------------------------------------------------------------------------
sfrotate <- function (x, degrees, centrexy = NULL, usecentroid = FALSE) {
rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
gx <- st_geometry(x)
if (is.null(centrexy)) {
if (usecentroid) {
centrexy <- st_centroid(gx)[1,] # unique centre
}
else {
centrexy <- st_centroid(st_as_sfc(st_bbox(x)))
}
}
else {
centrexy <- st_sfc(st_point(centrexy) )
}
(gx - centrexy) * rot(degrees/360*2*pi) + centrexy
}
#-------------------------------------------------------------------------------
# Based on Tim Salabim stackoverflow Jul 12 2018
# https://stackoverflow.com/questions/51292952/snap-a-point-to-the-closest-point-on-a-line-segment-using-sf
snap_points <- function(x, y, max_dist = 1000) {
if (inherits(x, "sf")) n = nrow(x)
if (inherits(x, "sfc")) n = length(x)
out = do.call(c,
lapply(seq(n), function(i) {
nrst = st_nearest_points(st_geometry(x)[i], y)
nrst_len = st_length(nrst)
nrst_mn = which.min(nrst_len)
if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
return(st_cast(nrst[nrst_mn], "POINT")[2])
})
)
return(out)
}
#-------------------------------------------------------------------------------
# random truncated Poisson
rtpois <- function(n, lambda) {
qpois(runif(n, dpois(0, lambda), 1), lambda)
}
#-------------------------------------------------------------------------------
im2mask <- function(im) {
# spatstat im object to mask
df <- as.data.frame(im)
names(df) <- c('x','y','Lambda')
df$Lambda <- df$Lambda * 1e4 # per hectare
read.mask(data = df, spacing = im$xstep)
}
#-------------------------------------------------------------------------------
# LIST 2023-03-10
# detectionfunctionname (fn)
# detectionfunctionnumber (detname)
# parnames (detectfn)
# getdfn (detectfn)
# valid.detectfn (detectfn, valid = c(0:3,5:19, 20))
# valid.detectpar (detectpar, detectfn)
# valid.model(model, CL, detectfn, hcov, userdist, sessioncovnames)
# getuserdistnames (userdist)
# valid.pnames (details, CL, detectfn, alltelem, sighting, nmix)
# valid.userdist (userdist, detector, xy1, xy2, mask, sessnum)
# new.param (details, model, CL)
# detectorcode (object, MLonly = TRUE, noccasions = NULL)
# expanddet(CH)
# ndetectpar (detectfn)
# replacedefaults (default, user)
# discreteN (n, N)
# ndetector (traps)
# memo (text, trace)
# insertdim (x, dimx, dims)
# pad1 (x, n)
# padarray (x, dims)
# stdform (flist)
# invlogit (y) 1/(1+exp(-y)) # plogis(y)
# logit (x) log(x/(1-x)) # qlogis(x), except for invalid argument
# sine (x) asin (x*2-1)
# invsine (y) (sin(y)+1) / 2
# odds (x) x / (1-x)
# invodds (y) y / (1+y)
# lnbinomial (x,size,prob)
# model.string (model, userDfn)
# fixed.string (fixed)
# var.in.model(v,m)
# get.nmix (model, capthist, hcov)
# fixpmix(x, nmix)
# add.cl (df, alpha, loginterval, lowerbound = 0)
# spatialscale (object, detectfn, session = '')
#
# # logical for whether object specifies userDfn
# userD (object)
# # mean and SD if x numeric
# getMeanSD(xy)
# nclusters (capthist)
# leadingzero (x)
# HN (r, pars, cutval)
# HR (r, pars, cutval)
# EX (r, pars, cutval)
# UN (r, pars, cutval)
# CHN (r, pars, cutval)
# WEX (r, pars, cutval)
# ANN (r, pars, cutval)
# CLN (r, pars, cutval)
# CG (r, pars, cutval)
# CN (r, pars, cutval)
# BSS (r, pars, cutval)
# SS (r, pars, cutval)
# SSS (r, pars, cutval)
# SN (r, pars, cutval)
# SNS (r, pars, cutval)
# HHN (r, pars, cutval)
# HHR (r, pars, cutval)
# HEX (r, pars, cutval)
# HAN (r, pars, cutval)
# HCG (r, pars, cutval)
# HVP (r, pars, cutval)
# HPX (r, pars, cutval)
# gradient (pars, fun, eps=0.001, ...)
# transform (x, link)
# se.transform (real, sereal, link)
# untransform (beta, link)
# se.untransform (beta, sebeta, link)
# Xtransform (real, linkfn, varnames)
# se.Xtransform (real, sereal, linkfn, varnames)
# Xuntransform (beta, linkfn, varnames)
# se.Xuntransform (beta, sebeta, linkfn, varnames)
# mlogit.untransform (beta, latentmodel)
# clean.mlogit(x)
# mlogit (x)
# group.levels (capthist, groups, sep='.')
# h.levels (capthist, hcov, nmix)
# n.occasion (capthist)
# group.factor (capthist, groups, sep='.')
# getgrpnum (capthist, groups)
# make.lookup (tempmat)
# getknownclass(capthist, nmix, hcov)
# getnmix (details)
# fullbeta (beta, fb)
# inflate (xy, rmult = 1)
# getbinomN (binomN, detectr)
# xy2CH (CH, inflation = 1e-8)
# inflatechull (poly, r, ntheta = 60)
# maskarea (mask, sessnum = 1)
# masklength (mask, sessnum = 1)
# masksize (mask, sessnum = 1)
# complete.beta (object)
# complete.beta.vcv (object)
# smooths (formula)
# polys (formula)
# badsmooths (formula)
# gamsetup(formula, data, ...)
# general.model.matrix (formula, data, gamsmth = NULL,
# makerealparameters (design, beta, parindx, link, fixed)
# secr.lpredictor (formula, newdata, indx, beta, field, beta.vcv=NULL,
# edist (xy1, xy2)
# nedist (xy1, xy2, mask, inf = Inf, ...)
# getcellsize (mask)
# updatemodel (model, detectfn, detectfns, oldvar, newvar, warn = FALSE)
# deleteMaskPoints (mask, onebyone = TRUE, add = FALSE, poly = NULL,
# nparameters (object)
# mapbeta (parindx0, parindx1, beta0, betaindex)
# xyinpoly (xy, trps)
# addzerodf (df, oldCH, sess)
# addzeroCH (CH, nzero, cov = NULL)
# expandbinomN (binomN, detectorcodes)
# check3D (object)
# updateCH(object)
# newstr <-function (strings)
# outsidemask(CH, mask, threshold = spacing(mask) / sqrt(2))
# shareFactorLevels (object, columns = NULL, stringsAsFactors = TRUE)
# allzero (object)
# primarysessions(intervals)
# secondarysessions(intervals)
# boundarytoSF (poly)
# pointsInPolygon (xy, poly, logical = TRUE)
# firstsk (PIAx)
# maskboolean (ch, mask, threshold)
# telemcode(object, ...)
# uniquerownames (capthist)
# selectCHsession(capthist, sessnum)
# stringsAsFactors (DF)
# getdistmat2 (traps, mask, userdist, HPX = FALSE)
# uniformusage(object, noccasions)
# sfrotate (x, degrees, centrexy = NULL, usecentroid = FALSE)
# snap_points(x, y, max_dist = 1000)
# rtpois(n, lambda)
Try the secr package in your browser
Any scripts or data that you put into this service are public.
secr documentation built on Oct. 18, 2023, 1:07 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.
Defines functions gamsetup getknownclass clean.mlogit getMeanSD fixpmix var.in.model expanddet valid.model
Documented in getMeanSD
################################################################################
## package 'secr'
## utility.R
################################################################################
## 2019-07-27 secr 4.0.0
## 2019-07-27 makelookupcpp replaces makelookup
## 2019-08-12 individualcovariates
## 2019-08-14 setNumThreads
## 2019-11-29 secr 4.1.0
## 2020-01-08 distancetotrap revised for polygon detectors
## 2020-01-26 selectCHsession
## 2020-02-21 secr 4.2.0
## 2020-05-15 stringsAsFactors function
## 2020-07-14 secr 4.3.0 distmat
## 2020-09-05 getknownclass factor bug fixed
## 2021-04-25 4.4.0
## 2021-12-16 tidy up transformations, allow arbitrary link X(), invX(), se.invX()
## 2022-01-04 4.5.0
## 2022-01-18 4.5.1
## 2022-01-23 uniformusage() function to assign all-ones usage matrix
## 2022-02-10 4.5.2 shift to sf where possible
## 2022-04-22 4.5.4
## 2022-05-31 4.5.5
## 2022-10-08 4.5.7
## 2022-11-29 4.5.8
## 2023-03-09 4.5.9
## 2023-03-10 4.5.11
## 2023-03-10 distancetotrap and nearesttrap moved to separate file
## 2023-03-10 setNumThreads moved to separate file
## 2023-05-21 4.6.0
## 2023-09-17 im2mask converts spatstat im object
################################################################################
# Global variables in namespace
#
## define a local environment for temporary variables e.g. iter
## e.g. Roger Peng https://stat.ethz.ch/pipermail/r-devel/2009-March/052883.html
.localstuff <- new.env()
## .localstuff$packageType <- ' pre-release'
.localstuff$packageType <- ''
.localstuff$validdetectors <- c('single','multi','proximity','count',
'polygonX', 'transectX', 'signal', 'polygon', 'transect',
'capped', 'null','null','null','null', 'telemetry', 'signalnoise')
.localstuff$simpledetectors <- c('single','multi','proximity','count', 'capped')
.localstuff$individualdetectors <- c('single','multi','proximity','count',
'polygonX', 'transectX', 'signal', 'signalnoise', 'polygon', 'transect',
'telemetry', 'capped')
.localstuff$pointdetectors <- c('single','multi','proximity','count',
'signal', 'signalnoise', 'unmarked','presence','capped')
.localstuff$polydetectors <- c('polygon','transect','polygonX','transectX')
.localstuff$exclusivedetectors <- c('single','multi','polygonX','transectX')
.localstuff$countdetectors <- c('count','polygon','transect','unmarked','telemetry')
.localstuff$iter <- 0 ## counter 1
.localstuff$iter2 <- 0 ## counter 2
.localstuff$detectionfunctions <-
c('halfnormal',
'hazard rate',
'exponential',
'compound halfnormal',
'uniform',
'w exponential',
'annular normal',
'cumulative lognormal',
'cumulative gamma',
'binary signal strength',
'signal strength',
'signal strength spherical',
'signal-noise',
'signal-noise spherical',
'hazard halfnormal',
'hazard hazard rate',
'hazard exponential',
'hazard annular normal',
'hazard cumulative gamma',
'hazard variable power',
'hazard pixelar')
.localstuff$DFN <- c('HN', 'HR', 'EX', 'CHN', 'UN', 'WEX', 'ANN', 'CLN', 'CG',
'BSS', 'SS', 'SSS', 'SN', 'SNS',
'HHN', 'HHR', 'HEX', 'HAN', 'HCG', 'HVP','HPX')
.localstuff$learnedresponses <- c('b', 'bk', 'B', 'k', 'Bk') ## Bk added 2020-02-26
#-------------------------------------------------------------------------------
detectionfunctionname <- function (fn) {
.localstuff$detectionfunctions[fn+1]
}
#-------------------------------------------------------------------------------
detectionfunctionnumber <- function (detname) {
dfn <- match (toupper(detname), .localstuff$DFN)
if (is.na(dfn))
dfn <- match (tolower(detname), .localstuff$detectionfunctions)
if (is.na(dfn))
stop ("unrecognised detection function ", detname)
dfn-1
}
#-------------------------------------------------------------------------------
parnames <- function (detectfn) {
parnames <- switch (detectfn+1,
c('g0','sigma'), ## 0
c('g0','sigma','z'),
c('g0','sigma'),
c('g0','sigma','z'),
c('g0','sigma'),
c('g0','sigma','w'),
c('g0','sigma','w'),
c('g0','sigma','z'),
c('g0','sigma','z'),
c('b0','b1'),
c('beta0','beta1', 'sdS'), ## include cutval?
c('beta0','beta1', 'sdS'), ## include cutval?
c('beta0','beta1', 'sdS','muN','sdN'),
c('beta0','beta1', 'sdS','muN','sdN'),
c('lambda0','sigma'),
c('lambda0','sigma','z'),
c('lambda0','sigma'),
c('lambda0','sigma','w'),
c('lambda0','sigma','z'),
c('lambda0','sigma','z'),
c('lambda0') ## 20
)
}
#-------------------------------------------------------------------------------
getdfn <- function (detectfn) {
switch (detectfn+1, HN, HR, EX, CHN, UN, WEX, ANN, CLN, CG, BSS, SS, SSS,
SN, SNS, HHN, HHR, HEX, HAN, HCG, HVP, HPX)
}
#-------------------------------------------------------------------------------
valid.detectfn <- function (detectfn, valid = c(0:3,5:19, 20)) {
# exclude 4 uniform: too numerically flakey
if (is.null(detectfn))
stop ("requires 'detectfn'")
if (is.character(detectfn))
detectfn <- detectionfunctionnumber(detectfn)
if (!(detectfn %in% valid))
stop ("invalid detection function")
detectfn
}
#-------------------------------------------------------------------------------
valid.detectpar <- function (detectpar, detectfn) {
if (is.null(detectpar) | is.null(detectfn))
stop ("requires valid 'detectpar' and 'detectfn'")
## 2013-07-19, 2013-10-22
## replace a0 with g0 or lambda0 as appropriate to detectfn
if ('a0' %in% names(detectpar)) {
aname <- if (detectfn %in% 0:8) 'g0' else 'lambda0'
lambda0 <- detectpar[['a0']] / 2 / pi / detectpar[[2]]^2 * 10000
detectpar[[aname]] <- if (detectfn %in% 0:8) 1-exp(-lambda0) else lambda0
}
if (!all(parnames(detectfn) %in% names(detectpar)))
stop ("requires 'detectpar' ", paste(parnames(detectfn), collapse=','),
" for ", detectionfunctionname(detectfn), " detectfn")
detectpar[parnames(detectfn)]
}
#-------------------------------------------------------------------------------
valid.model <- function(model, CL, detectfn, hcov, userdist, sessioncovnames) {
## 2014-08-25
if (any(sapply(model, badsmooths)))
warning("smooth term may be unsuitable for secr: does not specify k or fx where required")
}
#-------------------------------------------------------------------------------
getuserdistnames <- function (userdist) {
## return the names of any supplementary arguments of user-provided function
## for non-euclidean distance computations
if (is.function(userdist)) {
distnames <- try(userdist(), silent = TRUE)
if (!is.character(distnames))
stop("invalid userdist function - ",
"should return parameter names when called with no arguments")
distnames
}
else
character(0)
}
#-------------------------------------------------------------------------------
valid.pnames <- function (details, CL, detectfn, alltelem, sighting, nmix) {
## modelled parameters
pnames <- switch (detectfn+1,
c('g0','sigma'), # 0 halfnormal
c('g0','sigma','z'), # 1 hazard rate
c('g0','sigma'), # 2 exponential
c('g0','sigma','z'), # 3
c('g0','sigma'), # 4
c('g0','sigma','w'), # 5
c('g0','sigma','w'), # 6
c('g0','sigma','z'), # 7
c('g0','sigma','z'), # 8
c('b0','b1'), # 9
c('beta0','beta1','sdS'), # 10
c('beta0','beta1','sdS'), # 11
c('beta0','beta1','sdS'), # 12 cf parnames() in utility.R: muN, sdN?
c('beta0','beta1','sdS'), # 13 cf parnames() in utility.R: muN, sdN?
c('lambda0','sigma'), # 14 hazard halfnormal
c('lambda0','sigma','z'), # 15 hazard hazard rate
c('lambda0','sigma'), # 16 hazard exponential
c('lambda0','sigma','w'), # 17
c('lambda0','sigma','z'), # 18
c('lambda0','sigma','z'), # 19
c('lambda0','sigma')) # 20 hazard pixelar 2021-03-25
if (details$param %in% c(2,6))
pnames[1] <- 'esa'
if (details$param %in% c(3,5))
pnames[1] <- 'a0'
if (details$param %in% 4:6) {
pnames[2] <- 'sigmak'
pnames <- c(pnames, 'c')
pnames <- c(pnames, 'd')
}
if (!CL)
pnames <- c('D', pnames)
if ('noneuc' %in% getuserdistnames(details$userdist))
pnames <- c(pnames, 'noneuc')
if (sighting)
pnames <- c(pnames, 'pID')
# if (alltelem) {
# rnum <- match(c('D','lambda0','a0','esa','g0'), pnames)
# rnum[is.na(rnum)] <- 0
# pnames <- pnames[-rnum]
# }
if (nmix>1)
pnames <- c(pnames, 'pmix')
pnames
}
#-------------------------------------------------------------------------------
valid.userdist <- function (userdist, detector, xy1, xy2, mask, sessnum) {
if (is.null(userdist)) {
## default to Euclidean distance
result <- edist(xy1, xy2)
}
else {
if (any(detector %in% .localstuff$polydetectors)) {
stop ("userdist cannot be used with polygon detector types;")
}
if (is.function(userdist))
{
OK <- getuserdistnames(userdist) %in% names(covariates(mask))
if ((length(OK)>0) & !all(OK))
stop ("covariates required by userdist function not in mask : ",
paste(getuserdistnames(userdist)[!OK], collapse=','))
# 2023-02-06 selected columns 1:2 only (mask passes miscparm)
result <- do.call(userdist, c(list(xy1[,1:2], xy2[,1:2], mask)))
}
else {
if (is.character(userdist)) {
userdist <- get(userdist, pos=-1)
}
if (is.list(userdist) & !is.data.frame(userdist)) {
if (missing(sessnum) || is.na(sessnum))
stop("This use does not yet allow for session-specific userdist")
else
result <- userdist[[sessnum]]
}
else {
result <- userdist
}
}
if (!all(dim(result) == c(nrow(xy1), nrow(xy2))))
stop ("invalid distance matrix dim = ", dim(result)[1], ',', dim(result)[2])
baddist <- (!is.finite(result)) | (result<0) | is.na(result)
if (any(baddist)) {
warning ("replacing infinite, negative and NA userdist values with 1e10")
result[baddist] <- 1e10
}
}
result
}
#-------------------------------------------------------------------------------
new.param <- function (details, model, CL) {
esa <- 'esa' %in% names(model)
a0 <- 'a0' %in% names(model)
sigmak <- 'sigmak' %in% names(model)
newparam <- details$param
if (esa & !sigmak) {
newparam <- 2
}
if (a0 & !sigmak) {
newparam <- 3
}
if (sigmak) {
if (esa) {
newparam <- 6
}
else {
if (CL)
stop ("sigmak parameterization requires full model, not CL, unless also 'esa'")
newparam <- ifelse(a0, 5, 4)
}
}
if (newparam != details$param)
warning ("Using parameterization details$param = ", newparam)
newparam
}
#-------------------------------------------------------------------------------
# secr 3.0 2016-10-05
## NEEDS MULTI-SESSION FORM
detectorcode <- function (object, MLonly = TRUE, noccasions = NULL) {
## numeric detector code from a traps object
detcode <- sapply(detector(object), switch,
single = -1,
multi = 0,
proximity = 1,
count = 2,
polygonX = 3,
transectX = 4,
signal = 5,
polygon = 6,
transect = 7,
capped = 8,
unmarked = 10,
presence = 11,
signalnoise = 12,
telemetry = 13,
-2)
if (MLonly) {
detcode <- ifelse (detcode==-1, rep(0,length(detcode)), detcode)
if (any(detcode<0))
stop ("Unrecognised detector type")
}
if (!is.null(noccasions) & (length(detcode)==1))
detcode <- rep(detcode, noccasions)
detcode
}
#-------------------------------------------------------------------------------
expanddet <- function(CH) {
trps <- traps(CH)
if (is.null(trps))
return ('nonspatial')
else {
det <- detector(trps)
if (length(det)<ncol(CH))
rep(det[1], ncol(CH))
else det
}
}
#-------------------------------------------------------------------------------
ndetectpar <- function (detectfn) {
length(parnames(detectfn))
}
#-------------------------------------------------------------------------------
replacedefaults <- function (default, user) replace(default, names(user), user)
#-------------------------------------------------------------------------------
discreteN <- function (n, N) {
tN <- trunc(N)
if (N != tN) tN + sample (x = c(1,0), prob = c(N-tN, 1-(N-tN)),
replace = T, size = n)
else rep(tN,n)
}
#-------------------------------------------------------------------------------
ndetector <- function (traps) {
if (is.null(traps))
return(1)
else if (all(detector(traps) %in% .localstuff$polydetectors))
length(levels(polyID(traps)))
else
nrow(traps)
}
#-------------------------------------------------------------------------------
memo <- function (text, trace) {
## could use message(text), but does not immediately flush console
if (trace) { cat (text, '\n')
flush.console() }
}
#-------------------------------------------------------------------------------
insertdim <- function (x, dimx, dims) {
## make vector of values
## using x repeated so as to fill array
## with dim = dims and the x values occupying dimension(s) dimx
olddim <- 1:length(dims)
olddim <- c(olddim[dimx], olddim[-dimx])
temp <- array (dim=c(dims[dimx], dims[-dimx]))
tempval <- array(dim=dims[dimx])
if (length(x) > length(tempval))
tempval[] <- x[1:length(tempval)]
else
tempval[] <- x ## repeat as needed
temp[] <- tempval ## repeat as needed
if (is.factor(x))
factor(levels(x), levels=levels(x))[aperm(temp, order(olddim))] ## 2010 02 25
else
as.vector(aperm(temp, order(olddim)))
}
#-------------------------------------------------------------------------------
pad1 <- function (x, n) {
## pad x to length n with dummy (first value)
if (is.factor(x)) {
xc <- as.character(x)
xNA <- c(xc, rep(xc[1], n-length(xc)))
out <- factor(xNA, levels=levels(x))
}
else out <- c(x, rep(x[1], n-length(x)))
out
}
#-------------------------------------------------------------------------------
padarray <- function (x, dims) {
temp <- array(dim=dims)
dimx <- dim(x)
if (all(dimx>0)) {
if (length(dimx)<2 | length(dimx)>3)
stop ("invalid array")
if (length(dimx)>2) temp[1:dimx[1], 1:dimx[2], 1:dimx[3]] <- x
else temp[1:dimx[1], 1:dimx[2]] <- x
}
temp
}
#-------------------------------------------------------------------------------
## regularize a list of formulae
stdform <- function (flist) {
LHS <- function (form) {
trms <- as.character (form)
if (length(trms)==2) '' else trms[2]
}
RHS <- function (form) {
trms <- as.character (form)
## 2020-05-14 for compatibility with R 4.0
if (length(trms)==3) as.formula(paste(trms[c(1,3)], collapse = " ")) else form
}
lhs <- sapply(flist, LHS)
temp <- lapply(flist, RHS)
if (is.null(names(flist))) names(temp) <- lhs
else names(temp) <- ifelse(names(flist) == '', lhs, names(flist))
temp
}
#-------------------------------------------------------------------------------
## miscellaneous functions
invlogit <- function (y) 1/(1+exp(-y)) # plogis(y)
logit <- function (x) log(x/(1-x)) # qlogis(x), except for invalid argument
sine <- function (x) asin (x*2-1)
invsine <- function (y) (sin(y)+1) / 2
odds <- function (x) x / (1-x)
invodds <- function (y) y / (1+y)
#-------------------------------------------------------------------------------
lnbinomial <- function (x,size,prob) {
# dbinom allowing non-integer x, forcing log = TRUE
if (x <= size) {
lgamma (size+1) - lgamma (size-x+1) - lgamma (x+1) +
x * log(prob) + (size-x) * log (1-prob)
}
else {
-Inf
}
}
#-------------------------------------------------------------------------------
model.string <- function (model, userDfn) {
# 2023-04-16 Note: model should be a list
if (!is.null(userDfn)) {
if (!is.null(model$D))
model$D <- paste('~userD', userDfn('name'), sep='.')
}
temp <- paste (names(model), as.character(model), collapse=' ', sep='')
temp
}
#-------------------------------------------------------------------------------
fixed.string <- function (fixed) {
if (is.null(fixed) | length(fixed)==0) 'none'
else paste (names(fixed), as.character(fixed), collapse=', ', sep=' = ')
}
#-------------------------------------------------------------------------------
var.in.model <- function(v,m) v %in% unlist(lapply(m, all.vars))
#-------------------------------------------------------------------------------
get.nmix <- function (model, capthist, hcov) {
model$D <- NULL ## ignore density model
model$pmix <- NULL ## pmix alone cannot make this a mixture model
nmix <- 1
if (any(var.in.model('h2', model))) {
nmix <- 2
if (any(var.in.model('h3', model)))
stop ("do not combine h2 and h3")
}
if (any(var.in.model('h3', model))) {
nmix <- 3
}
if ((nmix == 1) & (!is.null(hcov))) {
if (ms(capthist))
capthist <- capthist[[1]]
if (is.factor(covariates(capthist)[,hcov]))
lev <- levels(covariates(capthist)[,hcov])
else
lev <- levels(factor(covariates(capthist)[,hcov]))
if (all(is.na(covariates(capthist)[,hcov])))
stop ("hcov missing for all individuals, but detection model invariant")
if (length(lev) < 2)
stop ("hcov covariate not found or has fewer than 2 levels")
if (length(lev) > 2)
warning ("hcov covariate has more than 2 levels; using only first two")
nmix <- 2
}
nmix
}
#-------------------------------------------------------------------------------
fixpmix <- function(x, nmix) {
## x is a list with component pmix that is a matrix (dataframe)
## with columns 'estimate' and 'se' (and possibly others)
## referring to the linear predictor of pmix (i.e. on mlogit
## scale) and rows corresponding to rows in newdata
## (i.e. arbitrary combinations of predictors, including mixture
## class h2 or h3)
####################################################
## It is necessary that newdata include all levels
## of the mixture class.
####################################################
## 2013-10-29
## assuming mixture is always last dimension...
## previously used in collate, model.average and predict.secr
## 2015-09-30 incorporated in secr.lpredictor
temp <- matrix(x$pmix[,'estimate'], ncol = nmix)
if (nmix==2) temp[,x$pmix[,'h2']] <- x$pmix[,'estimate']
if (nmix==3) temp[,x$pmix[,'h3']] <- x$pmix[,'estimate']
temp2 <- apply(temp, 1, clean.mlogit)
x$pmix[,'estimate'] <- as.numeric(t(temp2))
if (nmix==2)
x$pmix[as.numeric(x$pmix$h2)==1,'se'] <- x$pmix[as.numeric(x$pmix$h2)==2,'se']
else
x$pmix[,'se'] <- rep(NA, nrow(x$pmix)) ## don't know how
x
}
#-------------------------------------------------------------------------------
add.cl <- function (df, alpha, loginterval, lowerbound = 0) {
## add lognormal or standard Wald intervals to dataframe with columns
## 'estimate' and 'SE.estimate'
## lowerbound added 2011-07-15
z <- abs(qnorm(1-alpha/2))
if (loginterval) {
delta <- df$estimate - lowerbound
df$lcl <- delta / exp(z * sqrt(log(1 + (df$SE.estimate /
delta)^2))) + lowerbound
df$ucl <- delta * exp(z * sqrt(log(1 + (df$SE.estimate /
delta)^2))) + lowerbound
}
else {
df$lcl <- pmax(lowerbound, df$estimate - z * df$SE.estimate)
df$ucl <- df$estimate + z * df$SE.estimate
}
df
}
#-------------------------------------------------------------------------------
spatialscale <- function (object, detectfn, session = '') {
if (inherits(object, 'secr')) {
if (ms(object))
detpar <- detectpar(object)[[session]]
else
detpar <- detectpar(object)
cutval <- object$details$cutval
}
else {
detpar <- object
cutval <- object$cutval
}
if (!is.null(detpar$sigma)) detpar$sigma
else if (detectfn == 10) {
(cutval - detpar$beta0) / detpar$beta1
}
else if (detectfn == 11) {
d11 <- function(d, beta0, beta1, c) beta0 +
beta1 * (d-1) - 10 * log10(d^2) - c
interval <- c(0,10 * (cutval - detpar$beta0) / detpar$beta1)
uniroot (d11, interval, detpar$beta0, detpar$beta1, cutval)$root
}
else if (detectfn == 9) {
- 1 / detpar$b1
}
else stop ("unrecognised detectfn")
}
#-------------------------------------------------------------------------------
## logical for whether object specifies userDfn
userD <- function (object) {
if (!inherits(object, c('secr','ipsecr')))
stop ("requires fitted model")
!is.null(object$details$userDfn)
}
#-------------------------------------------------------------------------------
## mean and SD if x numeric
getMeanSD <- function(xy) {
MeanSD <- function (x) {
if (is.numeric(x))
c(mean(x, na.rm = TRUE), sd(x, na.rm = TRUE))
else
c(NA,NA)
}
as.data.frame (apply(xy, 2, MeanSD))
}
#-------------------------------------------------------------------------------
nclusters <- function (capthist) {
if (ms(capthist)) {
lapply(capthist, nclusters)
}
else {
nmash <- attr(capthist, 'n.mash')
ifelse (is.null(nmash), 1, length(nmash))
}
}
#-------------------------------------------------------------------------------
## clunky but effective re-write 2012-09-04, improved 2016-02-20, 2016-05-10
leadingzero <- function (x) {
xc <- as.character(x)
w <- max(nchar(xc))
n0 <- function(n) paste(rep('0',n), collapse='')
paste(sapply(w-nchar(xc), n0), x, sep='')
## or, 2016-01-15, 2016-02-20 BUT DOESN'T HANDLE NON-INTEGER 2016-05-10
# if (is.character(x)) x <- as.numeric(x)
# sprintf(paste("%0", w, "d", sep = ""), x)
}
#-------------------------------------------------------------------------------
## Detection functions
HN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]
g0 * exp (-r^2 / 2 / sigma^2)
}
HR <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
g0 * (1 - exp (-(r / sigma)^-z))
}
EX <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]
g0 * exp (-r / sigma)
}
UN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]
ifelse (r<=sigma, g0, 0)
}
CHN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
g0 * ( 1 - (1 - exp (-r^2 / 2 / sigma^2)) ^ z )
}
WEX <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; w <- pars[3]
ifelse(r<=w, g0, g0*exp (-(r-w) / sigma))
}
ANN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; w <- pars[3]
g0 * exp (-(r-w)^2 / 2 / sigma^2)
}
CLN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
CV2 <- (z/sigma)^2
sdlog <- log(1 + CV2)^0.5
meanlog <- log(sigma) - sdlog^2/2
g0 * plnorm(r, meanlog, sdlog, lower.tail = FALSE)
}
CG <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
g0 * pgamma(r, shape=z, scale=sigma/z, lower.tail = FALSE)
}
CN <- function (r, pars, cutval) {
g0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
x <- z * (r - sigma)
g0 * (1 + (1 - exp(x)) / (1 + exp(x)))/2
}
BSS <- function (r, pars, cutval) {
b0 <- pars[1]; b1 <- pars[2]
gam <- -(b0 + b1 * r);
pnorm (gam, mean=0, sd=1, lower.tail=FALSE)
}
SS <- function (r, pars, cutval) {
beta0 <- pars[1]; beta1 <- pars[2]; sdS <- pars[3]
if (is.null(cutval))
stop ("require 'details$cutval' for signal strength plot")
mu <- beta0 + beta1 * r
1 - pnorm (q=cutval, mean=mu, sd=sdS)
}
SSS <- function (r, pars, cutval) {
beta0 <- pars[1]; beta1 <- pars[2]; sdS <- pars[3]
if (is.null(cutval))
stop ("require 'details$cutval' for signal strength plot")
## spherical so assume distance r measured from 1 m
mu <- beta0 - 10 * log ( r^2 ) / 2.302585 + beta1 * (r-1)
mu[r<1] <- beta0
1 - pnorm (q=cutval, mean=mu, sd=sdS)
}
SN <- function (r, pars, cutval) {
beta0 <- pars[1]; beta1 <- pars[2]; sdS <- pars[3];
muN <- pars[4]; sdN <- pars[5]
muS <- beta0 + beta1 * r
1 - pnorm (q=cutval, mean=muS-muN, sd=sqrt(sdS^2+sdN^2))
}
SNS <- function (r, pars, cutval) {
beta0 <- pars[1]; beta1 <- pars[2]; sdS <- pars[3];
muN <- pars[4]; sdN <- pars[5]
## spherical so assume distance r measured from 1 m
muS <- beta0 - 10 * log ( r^2 ) / 2.302585 + beta1 * (r-1)
muS[r<1] <- beta0
1 - pnorm (q=cutval, mean=muS-muN, sd=sqrt(sdS^2+sdN^2))
}
HHN <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]
1 - exp(-lambda0 * exp (-r^2 / 2 / sigma^2))
}
HHR <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
1 - exp(-lambda0 * ( 1 - exp (-(r / sigma)^-z)))
}
HEX <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]
1 - exp(-lambda0 * exp (-r / sigma))
}
HAN <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]; w <- pars[3]
1 - exp(-lambda0 * exp (-(r-w)^2 / 2 / sigma^2))
}
HCG <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
lambda0 * pgamma(r, shape=z, scale=sigma/z, lower.tail = FALSE)
}
HVP <- function (r, pars, cutval) {
lambda0 <- pars[1]; sigma <- pars[2]; z <- pars[3]
1 - exp(-lambda0 * exp(-(r/sigma)^z))
}
HPX <- function (r, pars, cutval) {
g0 <- 1-exp(-pars[1])
radius <- pars[2]
ifelse (r<=radius, g0, 0) # circular, not square! crude approx
}
#-------------------------------------------------------------------------------
gradient <- function (pars, fun, eps=0.001, ...)
## quick & dirty 2009 09 14
## used by plot.secr for delta method limits
{
est <- pars
g <- pars
for (i in 1:length(est))
{
temp <- est[i]
if (temp != 0.0) delta <- eps * abs(temp)
else delta <- eps
est[i] <- temp - delta
fminus <- fun (est, ...)
est[i] <- temp + delta
fplus <- fun (est, ...)
g[i] <- (fplus - fminus) / (2.0 * delta)
est[i] <- temp;
}
g
}
#-------------------------------------------------------------------------------
# transformation tidy up 2021-12-16
# arbitrary link function specified with functions X, invX, se.invX
transform <- function (x, link) {
switch (link,
identity = x,
i1000 = x * 1000,
log = log(x),
neglog = log(-x),
logit = logit(x),
odds = odds(x),
sin = sine(x),
do.call(link, list(x))
)
}
#-------------------------------------------------------------------------------
# used only in model.average, modelAverage
se.transform <- function (real, sereal, link) {
switch (link,
identity = sereal,
i1000 = sereal / 1000,
log = log((sereal/real)^2 + 1)^0.5,
neglog = log((sereal/-real)^2 + 1)^0.5,
logit = sereal / real / (1 - real),
sin = NA,
do.call(paste0('se.',link), list(real, sereal) )
)
}
#-------------------------------------------------------------------------------
untransform <- function (beta, link) {
switch (link,
identity = beta,
i1000 = beta / 1000,
log = exp(beta),
neglog = -exp(beta),
logit = invlogit(beta),
odds = invodds(beta),
sin = invsine(beta),
do.call(paste0('inv',link), list(beta))
)
}
#-------------------------------------------------------------------------------
se.untransform <- function (beta, sebeta, link) {
# Approximate translation of SE to untransformed scale
# Delta method cf Lebreton et al 1992 p 77
switch (link,
identity = sebeta,
i1000 = sebeta / 1000,
log = exp(beta) * sqrt(exp(sebeta^2)-1),
neglog = exp(beta) * sqrt(exp(sebeta^2)-1),
logit = invlogit(beta) * (1-invlogit(beta)) * sebeta,
sin = NA, ####!!!!
do.call(paste0('se.inv', link), list(beta=beta, sebeta=sebeta))
)
}
#-------------------------------------------------------------------------------
# vectorized transformations
Xtransform <- function (real, linkfn, varnames) {
mapply(transform, real, linkfn[varnames])
}
se.Xtransform <- function (real, sereal, linkfn, varnames) {
mapply(se.transform, real, sereal, linkfn[varnames])
}
Xuntransform <- function (beta, linkfn, varnames) {
mapply(untransform, beta, linkfn[varnames])
}
se.Xuntransform <- function (beta, sebeta, linkfn, varnames)
{
if (length(beta)!=length(sebeta))
stop ("'beta' and 'sebeta' do not match")
if (!all(varnames %in% names(linkfn)))
stop ("'linkfn' component missing for at least one real variable")
mapply(se.untransform, beta, sebeta, linkfn[varnames])
}
#-------------------------------------------------------------------------------
mlogit.untransform <- function (beta, latentmodel) {
if (!missing(latentmodel)) {
for (i in unique(latentmodel))
beta[latentmodel==i] <- mlogit.untransform(beta[latentmodel==i])
beta
}
else {
## beta should include values for all classes (mixture components)
nmix <- length(beta)
if (sum(is.na(beta)) != 1) {
## require NA for a single reference class
rep(NA, length(beta))
}
else {
nonreference <- !is.na(beta) # not reference class
b <- beta[nonreference]
pmix <- numeric(nmix)
pmix[nonreference] <- exp(b) / (1+sum(exp(b)))
pmix[!nonreference] <- 1 - sum(pmix[nonreference])
pmix
}
}
}
#-------------------------------------------------------------------------------
clean.mlogit <- function(x) {
## 2014-08-19 for robustness...
if (is.na(x[2])) x[2] <- 1-x[1]
x[1] <- NA ## assumed reference class
logit(mlogit.untransform(x))
}
#-------------------------------------------------------------------------------
mlogit <- function (x) {
## return the mlogit of an unscaled vector of positive values
## 2013-04-14
logit(x/sum(x))
}
## End of miscellaneous functions
#-------------------------------------------------------------------------------
group.levels <- function (capthist, groups, sep='.') {
# 2016-06-05 use also for trap strata
if (inherits(capthist, 'list')) {
temp <- lapply(capthist, group.levels, groups, sep)
sort(unique(unlist(temp))) ## vector of global levels
}
else {
if (is.null(groups)) 0
else {
if (!all(groups %in% names(covariates(capthist))))
stop ("one or more grouping variables is missing ",
"from covariates")
temp <- as.data.frame(covariates(capthist)[,groups])
# omit null combinations, sort as with default of factor levels
sort(levels(interaction(temp, drop=T, sep=sep)))
}
}
}
#-------------------------------------------------------------------------------
h.levels <- function (capthist, hcov, nmix) {
## determine the first nmix levels of a factor individual covariate
if (is.null(hcov))
as.character(1:nmix)
else {
if (ms(capthist)) {
## take first session as we can assume factor covariates have same levels in
## all sessions
capthist <- capthist[[1]]
}
hcov <- covariates(capthist)[,hcov]
if (!is.factor(hcov)) {
warning("hcov was coerced to a factor")
hcov <- factor(hcov)
}
levels(hcov)[1:nmix]
}
}
#-------------------------------------------------------------------------------
n.occasion <- function (capthist) {
## return the number of sampling occasions for each session in capthist
if (inherits(capthist, 'list')) {
sapply(capthist, n.occasion)
}
else {
ncol(capthist)
}
}
#-------------------------------------------------------------------------------
group.factor <- function (capthist, groups, sep='.')
## convert a set of grouping factors to a single factor (g)
## levels common to all sessions
{
if (inherits(capthist, 'list')) {
temp <- lapply(capthist, group.factor, groups) ## recursive call
grouplevels <- group.levels(capthist, groups)
if (length(grouplevels)<2)
temp
else
# list; force shared factor levels on each component
lapply (temp, factor, levels=grouplevels)
}
else {
if (is.null(groups) | (length(groups)==0) )
return (factor(rep(1, nrow(capthist)), levels = 1)) # added levels 2017-04-18
temp <- as.data.frame(covariates(capthist)[,groups])
if (ncol(temp) != length(groups))
stop ("one or more grouping variables is missing from ",
"covariates(capthist)")
temp <- interaction(temp, drop=T, sep=sep) # omit null combinations
temp
}
}
#-------------------------------------------------------------------------------
getgrpnum <- function (capthist, groups) {
if (is.null(groups))
rep(1, nrow(capthist))
else
match(group.factor(capthist, groups), group.levels(capthist, groups))
}
#-------------------------------------------------------------------------------
make.lookup <- function (tempmat) {
## should add something to protect make.lookup from bad data...
nrw <- nrow(tempmat)
ncl <- ncol(tempmat)
nam <- colnames(tempmat)
df <- is.data.frame(tempmat)
if (df) {
lev <- lapply(tempmat, levels)
tempmat[] <- sapply(tempmat, as.numeric)
tempmat <- as.matrix(tempmat)
}
dimnames(tempmat) <- NULL
temp <- makelookupcpp(tempmat)
lookup <- temp$lookup
colnames(lookup) <- nam
if (df) {
lookup <- as.data.frame(lookup)
## restore factors
for (i in 1: length(lev))
if (!is.null(lev[[i]]))
lookup[,i] <- factor(lev[[i]][lookup[,i]], levels = lev[[i]])
}
list (lookup=lookup, index=temp$index)
}
#-------------------------------------------------------------------------------
## Return an integer vector of class membership defined by a categorical
## individual covariate in a capthist object. Individuals of unknown
## class (including those with class exceeding nmix) are coded 1,
## others as (class number + 1). When no mixture is specified (nmix == 1)
## all are coded as unknown.
## knownclass 1 'unknown'
## knownclass 2 'latent class 1'
## knownclass 3 'latent class 2'
getknownclass <- function(capthist, nmix, hcov) {
if (ms(capthist)) {
lapply(capthist, getknownclass, nmix = nmix, hcov = hcov)
}
else {
if ((nmix>1) & (!is.null(hcov))) {
## 2020-09-05 use as.factor() instead of factor() to coerce
## (if already factor, coercing with factor() loses old levels)
var <- as.factor(covariates(capthist)[,hcov])
tmp <- as.numeric(var) + 1
tmp[is.na(tmp) | (tmp>(nmix+1))] <- 1
attr(tmp,'levels') <- levels(factor(covariates(capthist)
[,hcov]))[1:nmix]
tmp
}
else
rep(1,nrow(capthist))
}
}
#-------------------------------------------------------------------------------
getnmix <- function (details) {
if (is.null(details$nmix))
1
else
details$nmix
}
#-------------------------------------------------------------------------------
## expand beta parameter vector using template of 'fixed beta'
## fixed beta fb input is missing (NA) for estimated beta parameters
fullbeta <- function (beta, fb) {
if (!is.null(fb)) {
fb[is.na(fb)] <- beta ## partial beta (varying only)
beta <- fb ## complete beta
}
beta
}
#-------------------------------------------------------------------------------
## inflate a convex outline along all radii by linear factor 'rmult'
inflate <- function (xy, rmult = 1) {
xy <- as.matrix(xy)
centre <- apply(xy, 2, mean)
xy <- sweep(xy, MARGIN = 2, STATS = centre, FUN = '-')
r <- apply(xy, 1, function(z) sqrt(sum(z^2)))
theta <- atan2 (xy[,2], xy[,1])
r <- r * rmult
xy <- cbind(r * cos(theta), r * sin(theta))
sweep(xy, MARGIN = 2, STATS = centre, FUN = '+')
}
#-------------------------------------------------------------------------------
## moved from pdot.R 2013-11-09
## scalar 2016-10-14
getbinomN <- function (binomN, detectr) {
if (any(detectr %in% .localstuff$countdetectors)) {
if (is.null(binomN))
return(0)
else if (binomN == 'usage')
return(1)
else
return(binomN)
}
else
return(1)
}
#-------------------------------------------------------------------------------
## convert telemetryxy attribute of a combined dataset into a standalone capthist
## TO DO: option of telemetry or polygon output
xy2CH <- function (CH, inflation = 1e-8) {
if (ms(CH)) {
out <- lapply(CH, xy2CH, inflation)
class(out) <- c('capthist', 'list')
out
}
else {
xylist <- telemetryxy(CH)
if (is.null(xylist))
stop ("requires 'telemetryxy' attribute")
n <- length(xylist)
neach <- sapply(xylist, nrow)
allxy <- do.call(rbind, xylist)
trps <- allxy[chull(allxy),]
trps <- rbind(trps, trps[1,,drop=F])
trps <- inflate(trps, 1 + inflation) ## see also telemetry.R
trps <- as.data.frame(trps)
dimnames(trps) <- list(1:nrow(trps), c('x','y'))
class(trps) <- c("traps","data.frame")
detector(trps) <- "polygon"
polyID(trps) <- factor(rep(1,nrow(trps)))
rown <- rep(names(xylist), neach)
newCH <- array(neach, dim = c(n, 1, 1))
attr(newCH, "detectedXY") <- allxy
if (!is.null(covariates(CH))) {
rowlookup <- match(names(xylist), rownames(CH))
covariates(newCH) <- covariates(CH)[rowlookup,, drop=FALSE]
}
class(newCH) <- "capthist"
traps(newCH) <- trps
newCH
}
}
#-------------------------------------------------------------------------------
inflatechull <- function (poly, r, ntheta = 60) {
theta <- (2*pi) * (1:ntheta) / ntheta
## add supernumerary vertices
temp <- data.frame(x = apply(expand.grid(poly$x, r * cos(theta)),1,sum),
y = apply(expand.grid(poly$y, r * sin(theta)),1,sum))
hull <- chull(temp)
temp[c(hull,hull[1]), ]
}
#-------------------------------------------------------------------------------
maskarea <- function (mask, sessnum = 1) {
if (!ms(mask)) nrow(mask) * attr(mask,'area')
else nrow(mask[[sessnum]]) * attr(mask[[sessnum]],'area')
}
#-------------------------------------------------------------------------------
masklength <- function (mask, sessnum = 1) {
if (!ms(mask)) nrow(mask) * attr(mask,'spacing')/1000
else nrow(mask[[sessnum]]) * attr(mask[[sessnum]],'spacing')/1000
}
#-------------------------------------------------------------------------------
masksize <- function (mask, sessnum = 1) {
if (inherits(mask, 'linearmask'))
masklength(mask, sessnum)
else
maskarea(mask, sessnum)
}
#-------------------------------------------------------------------------------
complete.beta <- function (object) {
fb <- object$details$fixedbeta
# modified 2022-04-02 for consistency with ipsecr
beta <- if (inherits(object, 'secr')) object$fit$par else object$beta
if (!is.null(fb)) {
nbeta <- length(fb)
fb[is.na(fb)] <- beta
beta <- fb
}
beta
}
#-------------------------------------------------------------------------------
complete.beta.vcv <- function (object) {
fb <- object$details$fixedbeta
if (!is.null(fb)) {
nbeta <- length(fb)
beta.vcv <- matrix(NA, nrow = nbeta, ncol = nbeta)
beta.vcv[is.na(fb[row(beta.vcv)]) & is.na(fb[col(beta.vcv)])] <- object$beta.vcv
}
else {
beta.vcv <- object$beta.vcv
}
beta.vcv
}
#-------------------------------------------------------------------------------
smooths <- function (formula) {
## which terms in formula are smooths?
## returns logical vector
labels <- attr(terms(formula), "term.labels")
if (length(labels) > 0)
sapply(labels, function (x) any(sapply(c("s\\(", "te\\(", "poly\\("), grepl, x)))
else
logical(0)
}
#-------------------------------------------------------------------------------
polys <- function (formula) {
## which terms in formula are orthogonal polynomials?
## returns logical vector
labels <- attr(terms(formula), "term.labels")
if (length(labels) > 0)
sapply(labels, grepl, pattern = "poly\\(")
else
logical(0)
}
#-------------------------------------------------------------------------------
badsmooths <- function (formula) {
## does smooth specification conform to secr requirements?
## returns TRUE/FALSE
labels <- attr(terms(formula), "term.labels")
if (length(labels) > 0) {
smoothterms <- sapply(labels, function (x)
any(sapply(c("s\\(", "te\\("), grepl, x)))
labels <- labels[smoothterms]
any(sapply(labels, function(x)
grepl("s\\(", x) & !grepl("k =", x))) |
any(sapply(labels, function(x)
grepl("te\\(", x) & (!grepl("fx = TRUE", x) | !grepl("k =", x))))
}
else
FALSE
}
#-------------------------------------------------------------------------------
gamsetup <- function(formula, data, ...) {
## use 'session' column as dummy LHS so gam does not gag
## (cf secrgam:::make.density.design.matrix)
## session is always present in detection data, must be added for D
if (is.null(data$session)) data$session <- rep(1,nrow(data))
formula <- update.formula(formula, session ~ .)
setup <- gam(formula, data = data, fit = FALSE, ...)
colnames(setup$X) <- setup$term.names
setup
}
#-------------------------------------------------------------------------------
general.model.matrix <- function (formula, data, gamsmth = NULL,
contrasts = NULL, ...) {
## A function to compute the design matrix for the model in
## 'formula' given the data in 'data'. This is merely the result
## of model.matrix() unless 'formula' includes smooth terms -- s()
## or te() as described in mgcv ?formula.gam.
## If smooth terms are present then the matrix may be based on a
## previous gam setup (provided in the argument 'gamsmth') or
## computed de novo with gam(..., fit = FALSE)
## note 2014-08-24
## orthogonal polynomials e.g. poly(x,2) are handled by model.matrix,
## but the information needed for prediction at new data is not
## saved by secr.fit, so predict.secr generally fails with message
## "'degree' must be less than number of unique points"
## head(eval(parse(text = attr(terms(~ poly(x,y, degree=2)),
## 'term.labels')[1]), env=possummask))
## 2014-08-24, 2014-09-09, 2017-11-30
## 2019-10-12 drop row names
## 2021-12-09 f optional argument
dots <- list(...)
if (any(polys(formula)))
stop ("orthogonal polynomials are temporarily blocked") ## 2014-09-12
if (any(smooths(formula))) {
if (is.null(gamsmth)) {
## setup knots etc from scratch
mat <- gamsetup(formula, data, ...)$X
}
else {
## fool predict.gam into generating the necessary
## predictor matrix from previous setup
class (gamsmth) <- 'gam'
gamsmth$coefficients <- rep(NA, ncol(gamsmth$X))
mat <- mgcv::predict.gam(gamsmth, newdata = data, type = 'lpmatrix')
colnames(mat) <- colnames(gamsmth$X)
}
}
else {
## model.matrix(formula, data, ...)
mat <- model.matrix(formula, data = data, contrasts.arg = contrasts)
}
rownames (mat) <- NULL
mat
}
#-------------------------------------------------------------------------------
makerealparameters <- function (design, beta, parindx, link, fixed) {
modelfn <- function(i) {
## linear predictor for real parameter i
Yp <- design$designMatrices[[i]] %*% beta[parindx[[i]]]
if (names(link)[i] == 'pmix') {
## 2013-04-14 index of class groups (pmix sum to 1.0 within latentmodel)
cols <- dimnames(design$designMatrices[[i]])[[2]]
h2 <- grep('.h2', cols, fixed=T)
h3 <- grep('.h3', cols, fixed=T)
h2c <- grep(':h2', cols, fixed=T)
h3c <- grep(':h3', cols, fixed=T)
h.cols <- c(h2,h3,h2c,h3c)
tmp <- design$designMatrices[[i]][,-h.cols, drop = FALSE]
tmph <- design$designMatrices[[i]][,h.cols, drop = FALSE]
## 2018-02-23 why as.numeric()?
latentmodel <- as.numeric(factor(apply(tmp,1,paste, collapse='')))
refclass <- apply(tmph,1,sum) == 0
Yp[refclass] <- NA
Yp <- mlogit.untransform(Yp, latentmodel)
Yp[design$parameterTable[,i]]
}
else {
Yp <- untransform(Yp, link[[i]])
Yp[design$parameterTable[,i]] ## replicate as required
}
}
## construct matrix of detection parameters
nrealpar <- length(design$designMatrices)
parindx$D <- NULL ## detection parameters only
link$D <- NULL ## detection parameters only
parindx$noneuc <- NULL ## detection parameters only
link$noneuc <- NULL ## detection parameters only
detectionparameters <- names(link)
fixed.dp <- fixed[detectionparameters[detectionparameters %in% names(fixed)]]
if (length(fixed.dp)>0)
for (a in names(fixed.dp)) ## bug fixed by adding this line 2011-09-28
link[[a]] <- NULL
if (length(link) != nrealpar)
stop ("number of links does not match design matrices")
if (nrealpar == 0) {
return(matrix(unlist(fixed.dp),nrow = 1))
}
temp <- sapply (1:nrealpar, modelfn)
if (nrow(design$parameterTable)==1) temp <- t(temp)
nrw <- nrow(temp)
## make new matrix and insert columns in right place
temp2 <- as.data.frame(matrix(nrow = nrw, ncol = length(detectionparameters)))
names(temp2) <- detectionparameters
temp2[ , names(design$designMatrices)] <- temp ## modelled
if (!is.null(fixed.dp) & length(fixed.dp)>0)
temp2[ , names(fixed.dp)] <- sapply(fixed.dp, rep, nrw) ## fixed
as.matrix(temp2)
}
#-------------------------------------------------------------------------------
secr.lpredictor <- function (formula, newdata, indx, beta, field, beta.vcv=NULL,
smoothsetup = NULL, contrasts = NULL, Dfn = NULL) {
## form linear predictor for a single 'real' parameter
## smoothsetup should be provided whenever newdata differs from
## data used to fit model and the model includes smooths from gam
vars <- all.vars(formula)
OK <- vars %in% names(newdata)
if (any(!OK)) {
missingvars <- paste(vars[!OK], collapse = ', ')
if (sum(!OK) == 1)
stop ("model covariate ", missingvars, " not found in 'newdata'")
else
stop ("model covariates ", missingvars, " not found in 'newdata'")
}
newdata <- as.data.frame(newdata)
lpred <- matrix(ncol = 2, nrow = nrow(newdata), dimnames = list(NULL,c('estimate','se')))
if (!is.null(Dfn) && field == 'D') {
warning("secr.lpredictor is not ready for D as function - do not use estimates")
nsess <- length(unique(newdata$session))
Yp <- Dfn(newdata[,vars[1]], beta = beta[indx], dimD = c(nrow(newdata)/nsess,1,nsess))
mat <- as.matrix(newdata[,vars[1], drop = FALSE])
}
else {
mat <- general.model.matrix(formula, data = newdata, gamsmth = smoothsetup,
contrasts = contrasts)
if (nrow(mat) < nrow(newdata))
warning ("missing values in predictors?")
nmix <- 1
if (field=='pmix') {
## drop pmix beta0 column from design matrix (always zero)
mat <- mat[,-1,drop=FALSE]
if ('h2' %in% names(newdata)) nmix <- 2
if ('h3' %in% names(newdata)) nmix <- 3
mixfield <- c('h2','h3')[nmix-1]
}
###############################
Yp <- mat %*% beta[indx]
###############################
## A latent model comprises one row for each latent class.
## Back transformation of pmix in mlogit.untransform() requires all rows of
## each latent model. That function splits vector Yp by latent model.
if (field == 'pmix') {
nonh <- newdata[, names(newdata) != mixfield, drop = FALSE]
latentmodel <- factor(apply(nonh, 1, paste, collapse = ''))
refclass <- as.numeric(newdata[, mixfield]) == 1
Yp[refclass] <- NA
Yp <- mlogit.untransform(Yp, latentmodel)
Yp <- logit(Yp) # return to logit scale for later untransform!
if (nmix==2) {
h2.1 <- as.numeric(newdata$h2)==1
h2.2 <- as.numeric(newdata$h2)==2
}
}
}
lpred[,1] <- Yp
if (is.null(beta.vcv) || (any(is.na(beta[indx])))) return ( cbind(newdata,lpred) )
else {
if (is.null(Dfn) || field != 'D') {
vcv <- beta.vcv[indx,indx, drop = FALSE]
vcv[is.na(vcv)] <- 0
nrw <- nrow(mat)
vcv <- apply(expand.grid(1:nrw, 1:nrw), 1, function(ij)
mat[ij[1],, drop=F] %*% vcv %*% t(mat[ij[2],, drop=F]))
vcv <- matrix (vcv, nrow = nrw)
if (field=='pmix') {
if (nmix==2)
vcv[h2.1,h2.1] <- vcv[h2.2,h2.2]
else
vcv[,] <- NA
}
lpred[,2] <- diag(vcv)^0.5
}
else {
vcv <- NULL
}
temp <- cbind(newdata,lpred)
attr(temp, 'vcv') <- vcv
return(temp)
}
}
#-------------------------------------------------------------------------------
edist <- function (xy1, xy2) {
nr <- nrow(xy1)
nc <- nrow(xy2)
x1 <- matrix(xy1[,1], nr, nc)
x2 <- matrix(xy2[,1], nr, nc, byrow=T)
y1 <- matrix(xy1[,2], nr, nc)
y2 <- matrix(xy2[,2], nr, nc, byrow=T)
sqrt((x1-x2)^2 + (y1-y2)^2)
}
#-------------------------------------------------------------------------------
## least cost paths from mask including barriers to movement
## use edist for equivalent Euclidean distances
nedist <- function (xy1, xy2, mask, inf = Inf, ...) {
newargs <- list(...)
if (missing(mask)) mask <- xy2
noneuc <- covariates(mask)$noneuc
if (is.null(noneuc)) noneuc <- rep(1, nrow(mask))
defaultargs <- list(transitionFunction = mean, directions = 16)
args <- replace(defaultargs, names(newargs), newargs)
args$x <- raster(mask, values = noneuc)
if (requireNamespace('gdistance', quietly = TRUE)) { ## 2015-01-23
tr <- do.call(gdistance::transition, args)
tr <- gdistance::geoCorrection(tr, type = "c", multpl = FALSE)
out <- gdistance::costDistance(tr, as.matrix(xy1), as.matrix(xy2))
}
else stop ("package gdistance is required for nedist")
if (is.finite(inf)) out[!is.finite(out)] <- inf
out
}
#-------------------------------------------------------------------------------
getcellsize <- function (mask) {
if (inherits(mask, 'linearmask'))
cell <- attr(mask, 'spacing') / 1000 ## per km
else
cell <- attr(mask, 'area') ## per ha
if (is.null(cell))
stop ("mask lacks valid cell size (area or spacing)")
cell
}
#-------------------------------------------------------------------------------
## intercept and fix certain models with bad defaults
updatemodel <- function (model, detectfn, detectfns, oldvar, newvar, warn = FALSE) {
if (detectfn %in% detectfns) {
for (i in 1:length(oldvar)) {
if (oldvar[i] %in% names(model)) {
names(model)[names(model) == oldvar[i]] <- newvar[i]
if (warn)
warning ("replacing ", oldvar[i], " by ", newvar[i],
" in model for detectfn ", detectfn)
}
}
}
model
}
#-------------------------------------------------------------------------------
## Manually remove some mask points
# simplified 2022-02-03
deleteMaskPoints <- function (mask, onebyone = TRUE, add = FALSE, poly = NULL,
poly.habitat = FALSE, ...) {
## interface does not work properly in RStudio
if (ms(mask)) { ## a list of mask objects
if (inherits(poly, 'list') & (!is.data.frame(poly)))
stop ("lists of polygons not implemented in 'make.mask'")
temp <- lapply (mask, deleteMaskPoints, onebyone = onebyone, add = add,
poly = poly, poly.habitat = poly.habitat, ...)
class (temp) <- c('mask', 'list')
temp
}
else {
plot(mask, add = add, ...)
if (!is.null(poly)) {
if (poly.habitat)
pointstodrop <- (1:nrow(mask))[!pointsInPolygon(mask, poly)]
else
pointstodrop <- (1:nrow(mask))[pointsInPolygon(mask, poly)]
}
else if (onebyone) {
cat ('Click to select points; right-click to stop\n')
flush.console()
xy <- locator(type = 'p', pch=1, col='red')
pointstodrop <- if (length(xy$x)==0)
numeric(0)
else
nearesttrap(xy, mask)
}
else {
cat ('Click to select polygon vertices; right-click to stop\n')
flush.console()
xy <- locator(type = 'l', col='red')
xy <- as.data.frame(xy)
xy <- rbind(xy, xy[1,])
if (poly.habitat)
pointstodrop <- (1:nrow(mask))[!pointsInPolygon(mask, xy)]
else
pointstodrop <- (1:nrow(mask))[pointsInPolygon(mask, xy)]
}
npts <- length(pointstodrop)
if (npts>0) {
points(mask[pointstodrop,], pch = 16, col = 'red')
if(.Platform$OS.type == "windows") {
pl <- if (npts>1) 's' else ''
msg <- paste ('Delete ', npts, ' red point',pl, '?', sep='')
response <- utils::winDialog(type = "okcancel", msg)
} else {
response <- 'OK'
}
if (response == 'OK') {
mask <- subset(mask, -pointstodrop)
if (npts==1)
message("1 point deleted")
else
message(npts, " points deleted")
}
else
message ("point(s) not deleted")
}
else
message ("no points to delete")
plot(mask, col='green')
mask
}
}
#-------------------------------------------------------------------------------
nparameters <- function (object) {
Npar <- max(unlist(object$parindx))
Npar <- Npar + length(object$details$miscparm)
## allow for fixed beta parameters
if (!is.null(object$details$fixedbeta))
Npar <- Npar - sum(!is.na(object$details$fixedbeta))
Npar
}
#-------------------------------------------------------------------------------
mapbeta <- function (parindx0, parindx1, beta0, betaindex)
## Extend beta vector from simple model (beta0) to a more complex (i.e. general)
## model, inserting neutral values (zero) as required.
## For each real parameter, a 1:1 match is assumed between
## beta values until all beta values from the simpler model are
## used up. THIS ASSUMPTION MAY NOT BE JUSTIFIED.
## betaindex is a user-controlled alternative.
{
## list of zeroed vectors, one per real parameter
beta1 <- lapply(parindx1, function (x) {x[]<-0; x})
if (!is.null(betaindex)) {
beta1 <- unlist(beta1)
if (sum(betaindex>0) != length(beta0))
stop ("invalid 'betaindex'")
beta1[betaindex] <- beta0
beta1
}
else {
## indx is within-parameter rather than absolute index
## for each _original_ real parameter
indx <- lapply(parindx0, function(x) x-x[1]+1)
## for (j in 1:length(beta1))
## improved replace by name2015-11-17
for (j in names(beta1)) {
if (j %in% names(beta0))
beta1[[j]][indx[[j]]] <- beta0[parindx0[[j]]]
}
unlist(beta1)
}
}
#-------------------------------------------------------------------------------
xyinpoly <- function (xy, trps) {
ptinside <- function (i,k) {
## is point i inside poly k?
polyxy <- as.matrix(lxy[[k]])
polyxy <- rbind(polyxy, polyxy[1,]) ## close 2014-08-28
nr <- nrow(polyxy)
temp <- insidecpp(unlist(xy[i,]), 0, nr-1, as.matrix(polyxy))
}
lxy <- split (trps, polyID(trps))
firstinside <- function (i) {
frstk <- 0
for (k in 1:length(lxy)) {
if (ptinside(i,k)) {
frstk <- k
break
}
}
frstk
}
sapply(1:nrow(xy), firstinside)
}
#-------------------------------------------------------------------------------
addzerodf <- function (df, oldCH, sess) {
## add dummy detection records to dataframe for 'all-zero' case
## that arises in sighting-only mark-resight with known marks
allzero <- apply(oldCH,1,sum)==0
naz <- sum(allzero)
if (naz > 0) {
df0 <- expand.grid(
newID = rownames(oldCH)[allzero],
newocc = NA,
newtrap = trap(oldCH)[1],
alive = TRUE,
sess = sess,
stringsAsFactors = FALSE)
df$x <- NULL; df$y <- NULL ## 2021-04-08
df <- rbind(df,df0)
if (!is.null(xy(oldCH))) {
df$x <- c(xy(oldCH)$x, rep(NA, naz))
df$y <- c(xy(oldCH)$y, rep(NA, naz))
}
if (!is.null(signal(oldCH))) {
df$signal <- c(signal(oldCH), rep(NA, naz))
}
}
df
}
#-------------------------------------------------------------------------------
## including pre-marked animals never sighted
## cov is optional dataframe of covariates
addzeroCH <- function (CH, nzero, cov = NULL) {
if (nzero == 0)
return(CH)
else {
nc <- nrow(CH)
chdim <- dim(CH)
chdim[1] <- nzero
extra <- array(0, dim=chdim)
dimnames(extra) <- c(list(paste('Z', 1:nzero, sep='')), dimnames(CH)[2:3])
CH2 <- abind(CH, extra, along = 1)
class(CH2) <- 'capthist'
traps(CH2) <- traps(CH)
xy(CH2) <- xy(CH) ## order is not affected by adding zero histories
if (!is.null(covariates(CH)) & (nrow(CH)>0)) {
if (is.null(cov)) {
cov <- covariates(CH)[rep(1,nzero),]
cov[,] <- NA ## covariates are unknown
}
covariates(CH2) <- rbind(covariates(CH), cov[1:nzero,])
}
## ... and other essential attributes?
CH2
}
}
#-------------------------------------------------------------------------------
expandbinomN <- function (binomN, detectorcodes) {
# assumes detectorcodes is a vector of length = noccasions
binomN <- ifelse (detectorcodes %in% c(2,6,7), binomN, 1)
if (any(is.na(binomN))) stop ("NA value in binomN")
binomN
}
#-------------------------------------------------------------------------------
check3D <- function (object) {
if (ms(object)) {
out <- lapply(object, check3D)
class(out) <- class(object)
out
}
else {
if (is.matrix(object)) {
warning("secr >= 3.0 requires 3-D capthist; using updateCH() to convert")
updateCH(object)
}
else {
object
}
}
}
#-------------------------------------------------------------------------------
updateCH <- function(object) {
if (!inherits(object, 'capthist'))
stop ("requires capthist object")
# following replaces this old code 2020-08-29
# reduce(object, dropunused = FALSE)
if (ms(object)) {
out <- lapply(object, updateCH)
class (out) <- c("capthist", "list")
out
}
else {
if (length(dim(object)) == 3) {
return(object)
}
else {
K <- ndetector(traps(object))
ch <- array(0, dim = c(dim(object), K), dimnames =
list(rownames(object), colnames(object), 1:K))
OK <- as.logical(object!=0)
animal <- row(object)[OK]
occ <- col(object)[OK]
detn <- object[OK]
ch[cbind(animal, occ, detn)] <- 1
traps(ch) <- traps(object)
class (ch) <- "capthist"
session(ch) <- session(object)
ch
}
}
}
#-------------------------------------------------------------------------------
newstr <-function (strings) {
## compress a character vector
## use run length encoding function
rl <- rle(strings)
st <- rl$values
le <- paste0(' (',as.character(rl$lengths), ')')
le[le==' (1)'] <- ''
paste(paste0(st, le), collapse = ', ')
}
# newstr(c("single", rep("proximity",4)))
#-------------------------------------------------------------------------------
outsidemask <- function(CH, mask, threshold = spacing(mask) / sqrt(2)) {
xylist <- telemetryxy(CH)
dfun <- function(xy) {
centres <- matrix(apply(xy, 2, mean), ncol = 2)
distancetotrap(centres, mask)
}
sapply(xylist, dfun) > threshold
}
#-------------------------------------------------------------------------------
shareFactorLevels <- function (object, columns = NULL, stringsAsFactors = TRUE) {
## stringsAsFactors added 2020-05-16
if (ms(object)) {
if (!is.null(covariates(object))) {
df <- do.call(rbind, covariates(object))
if (is.null(columns)) {
columns <- 1:ncol(df)
}
if (stringsAsFactors) {
df[,columns] <- stringsAsFactors(df[,columns, drop = FALSE])
}
for (i in columns) {
if (is.factor(df[,i])) {
levelsi <- levels(df[,i])
for (sess in 1:length(object)) {
covariates(object[[sess]])[,i] <-
factor(covariates(object[[sess]])[,i],
levels = levelsi)
}
}
}
}
}
else {
# modified 2021-04-27 to apply to covariates, not object itself
if (!is.null(covariates(object))) {
if (stringsAsFactors) {
df <- covariates(object)
if (is.null(columns)) {
columns <- 1:ncol(df)
}
df[,columns] <- stringsAsFactors(df[,columns, drop = FALSE])
covariates(object) <- df
}
}
}
object
}
#-------------------------------------------------------------------------------
allzero <- function (object) {
if (!inherits(object, 'capthist'))
stop ("requires 'capthist' object")
if (ms(object)) {
lapply(object, allzero)
}
else {
telemocc <- detector(traps(object))=='telemetry'
apply(object[,!telemocc,,drop=FALSE],1,sum)==0
}
}
#-------------------------------------------------------------------------------
primarysessions <- function(intervals) {
primarysession <- cumsum(c(0,intervals))
match(primarysession, unique(primarysession))
}
#-------------------------------------------------------------------------------
secondarysessions <- function(intervals) {
primary <- primarysessions(intervals)
unname(unlist(sapply(table(primary), seq_len)))
}
#-------------------------------------------------------------------------------
boundarytoSF <- function (poly) {
if (is.null(poly)) {
NULL
}
else if(inherits(poly, c('sf','sfc'))) {
poly <- st_geometry(poly) # extract sfc if not already sfc
geomtype <- st_geometry_type(poly, by_geometry = FALSE)
if (geomtype == 'GEOMETRY') { # 2023-06-02
geomtype <- st_geometry_type(poly, by_geometry = TRUE)
}
if (!all(geomtype %in% c("POLYGON", "MULTIPOLYGON"))) {
stop ("poly sf/sfc should be of type POLYGON or MULTIPOLYGON")
}
poly
}
else if (inherits(poly, 'SpatialPolygons')) { # also SPDF?
st_as_sfc(poly)
}
else if (inherits(poly, 'SpatVector')) {
st_as_sfc(as(poly,"Spatial"))
}
else if (inherits(poly, c('matrix', 'data.frame'))) {
## input is 2-column matrix for a single polygon
poly <- matrix(unlist(poly), ncol = 2)
poly <- rbind (poly, poly[1,]) ## force closure of polygon
st_sfc(st_polygon(list(poly)))
}
else stop (class(poly), " not valid input to boundarytoSF")
}
#-------------------------------------------------------------------------------
pointsInPolygon <- function (xy, poly, logical = TRUE) {
# xy is 2-column matrix or data.frame of coordinates
if (inherits(poly, 'mask')) {
if (ms(poly))
stop ("multi-session masks not supported")
sp <- spacing(poly)
minx <- min(poly$x, na.rm = TRUE)
miny <- min(poly$y, na.rm = TRUE)
mask <- sweep(poly, MARGIN = 2, FUN = '+', STATS = c(-minx, -miny))
mask <- round(mask/sp) + 1
xy <- matrix(unlist(xy), ncol = 2) ## in case dataframe
xy <- sweep(xy, MARGIN = 2, FUN = '+', STATS = c(-minx, -miny))
xy <- round(xy/sp) + 1
xy[xy<=0] <- NA
xy[,1][xy[,1]>max(mask$x, na.rm = TRUE)] <- NA
xy[,2][xy[,2]>max(mask$y, na.rm = TRUE)] <- NA
maskmatrix <- matrix(0, ncol = max(mask$y, na.rm = TRUE), nrow = max(mask$x, na.rm = TRUE))
maskmatrix[as.matrix(mask)] <- 1:nrow(mask)
inside <- maskmatrix[as.matrix(xy)]
inside[is.na(inside)] <- 0
if (logical)
inside <- inside > 0
inside
}
else {
poly <- boundarytoSF(poly)
if (inherits(poly, c('sf','sfc'))) {
xy <- st_as_sf(data.frame(xy), coords = 1:2)
st_crs(xy) <- st_crs(poly)
apply(st_within(xy, poly, sparse = FALSE), 1, any)
}
else {
stop ("unknown input to pointsInPolygon")
}
}
}
#-------------------------------------------------------------------------------
## return indices of first occasion and detector for which PIAx is non-zero
firstsk <- function (PIAx) {
## PIAx dim n,s,k
wh <- function(d2) {
match(TRUE, d2>0)
}
apply(PIAx,1,wh)
}
#-------------------------------------------------------------------------------
maskboolean <- function (ch, mask, threshold) {
if (ms(ch)) {
if (!ms(mask)) stop ("masklookup: multisession ch requires multisession mask")
outlist <- mapply(maskboolean, ch, mask, MoreArgs = list(threshold = threshold), SIMPLIFY = FALSE)
outlist
}
else {
id <- animalID(ch, names = FALSE, sortorder = 'snk')
tr <- trap(ch, names = FALSE, sortorder = 'snk')
trps <- traps(ch)
m <- nrow(mask)
if (!is.null(threshold) && all(detector(trps) %in% .localstuff$pointdetectors)) {
df <- data.frame(id = id, x = trps$x[tr], y = trps$y[tr])
x <- tapply(df$x, df$id, mean, na.rm=T)
y <- tapply(df$y, df$id, mean, na.rm=T)
xy <- data.frame(x=x,y=y)
d2 <- edist2cpp(as.matrix(xy), as.matrix(mask))
out <- (d2 <= threshold^2)
}
else {
## NULL option
out <- matrix(TRUE, nrow = nrow(ch), ncol = m)
}
out
}
}
#-------------------------------------------------------------------------------
telemcode <- function(object, ...) {
if (inherits(object, 'traps') && !ms(object))
switch (telemetrytype(object), none = 0,
independent = 1, dependent = 2, concurrent = 3, 0)
else
NA
}
#-------------------------------------------------------------------------------
uniquerownames <- function (capthist) {
if (!ms(capthist)) {
return(capthist)
}
else {
last <- 0
for (i in 1:length(capthist)) {
nr <- nrow(capthist[[i]])
if (nr > 0) {
rownames(capthist[[i]]) <- last + (1:nr)
last <- last+nr
}
}
capthist
}
}
#-------------------------------------------------------------------------------
selectCHsession <- function(capthist, sessnum) {
if (ms(capthist))
capthist[[sessnum]]
else
capthist
}
#-------------------------------------------------------------------------------
stringsAsFactors <- function (DF) {
# convert any character columns of a data.frame (or list) to factor
if (is.list(DF) && length(DF)>0) { ## bug fix 2020-08-14
chr <- sapply(DF, is.character)
DF[chr] <- lapply(DF[chr], as.factor)
}
DF
}
#-------------------------------------------------------------------------------
# see also getuserdist in loglikhelperfn.R
# 2021-03-30 moved from preparedata.R
# 2021-03-30 integrate HPX
getdistmat2 <- function (traps, mask, userdist, HPX = FALSE) {
## Static distance matrix
if (is.function(userdist)) {
NULL ## compute dynamically later
}
else {
if (HPX) {
if (any(detector(traps) %in% .localstuff$polydetectors)) {
trps <- split(traps, polyID(traps))
inside <- t(sapply(trps, pointsInPolygon, xy = mask))
d2 <- 1-inside # 0 inside, 1 outside
d2[d2>0] <- 1e10 # 0 inside, 1e10 outside
d2
}
else {
# maximum of squared distance in x- or y- directions
xydist2cpp(as.matrix(traps), as.matrix(mask))
}
}
else {
if (any(detector(traps) %in% .localstuff$polydetectors)) {
## do not use result if detector is one of
## polygonX, polygon, transectX, transect, OR telemetry?
matrix(0, nrow = nrow(traps), ncol = nrow(mask))
}
else {
# Euclidean distance
## edist(as.matrix(traps), as.matrix(mask))^2
edist2cpp(as.matrix(traps), as.matrix(mask))
}
}
}
}
#-------------------------------------------------------------------------------
## function to assign all-ones usage matrix
uniformusage <- function(object, noccasions) {
if (inherits(object, 'capthist')) {
if (ms(object)) {
for (r in 1:length(object)) {
ndet <- dim(object[[r]])[3]
noccasions <- dim(object[[r]])[2]
usage(traps(object[[r]])) <- matrix(1, ndet, noccasions)
}
}
else {
ndet <- dim(object)[3]
noccasions <- dim(object)[2]
usage(traps(object)) <- matrix(1, ndet, noccasions)
}
}
else if (inherits(object, 'traps')) {
if (missing(noccasions)) {
stop ('noccasions should be specified for traps input')
}
if (ms(object)) {
for (r in 1:length(object)) {
ndet <- ndetector(object[[r]])
usage(object[[r]]) <- matrix(1, ndet, noccasions)
}
}
else {
ndet <- ndetector(object)
usage(object) <- matrix(1, ndet, noccasions)
}
}
object
}
#-------------------------------------------------------------------------------
sfrotate <- function (x, degrees, centrexy = NULL, usecentroid = FALSE) {
rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
gx <- st_geometry(x)
if (is.null(centrexy)) {
if (usecentroid) {
centrexy <- st_centroid(gx)[1,] # unique centre
}
else {
centrexy <- st_centroid(st_as_sfc(st_bbox(x)))
}
}
else {
centrexy <- st_sfc(st_point(centrexy) )
}
(gx - centrexy) * rot(degrees/360*2*pi) + centrexy
}
#-------------------------------------------------------------------------------
# Based on Tim Salabim stackoverflow Jul 12 2018
# https://stackoverflow.com/questions/51292952/snap-a-point-to-the-closest-point-on-a-line-segment-using-sf
snap_points <- function(x, y, max_dist = 1000) {
if (inherits(x, "sf")) n = nrow(x)
if (inherits(x, "sfc")) n = length(x)
out = do.call(c,
lapply(seq(n), function(i) {
nrst = st_nearest_points(st_geometry(x)[i], y)
nrst_len = st_length(nrst)
nrst_mn = which.min(nrst_len)
if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
return(st_cast(nrst[nrst_mn], "POINT")[2])
})
)
return(out)
}
#-------------------------------------------------------------------------------
# random truncated Poisson
rtpois <- function(n, lambda) {
qpois(runif(n, dpois(0, lambda), 1), lambda)
}
#-------------------------------------------------------------------------------
im2mask <- function(im) {
# spatstat im object to mask
df <- as.data.frame(im)
names(df) <- c('x','y','Lambda')
df$Lambda <- df$Lambda * 1e4 # per hectare
read.mask(data = df, spacing = im$xstep)
}
#-------------------------------------------------------------------------------
# LIST 2023-03-10
# detectionfunctionname (fn)
# detectionfunctionnumber (detname)
# parnames (detectfn)
# getdfn (detectfn)
# valid.detectfn (detectfn, valid = c(0:3,5:19, 20))
# valid.detectpar (detectpar, detectfn)
# valid.model(model, CL, detectfn, hcov, userdist, sessioncovnames)
# getuserdistnames (userdist)
# valid.pnames (details, CL, detectfn, alltelem, sighting, nmix)
# valid.userdist (userdist, detector, xy1, xy2, mask, sessnum)
# new.param (details, model, CL)
# detectorcode (object, MLonly = TRUE, noccasions = NULL)
# expanddet(CH)
# ndetectpar (detectfn)
# replacedefaults (default, user)
# discreteN (n, N)
# ndetector (traps)
# memo (text, trace)
# insertdim (x, dimx, dims)
# pad1 (x, n)
# padarray (x, dims)
# stdform (flist)
# invlogit (y) 1/(1+exp(-y)) # plogis(y)
# logit (x) log(x/(1-x)) # qlogis(x), except for invalid argument
# sine (x) asin (x*2-1)
# invsine (y) (sin(y)+1) / 2
# odds (x) x / (1-x)
# invodds (y) y / (1+y)
# lnbinomial (x,size,prob)
# model.string (model, userDfn)
# fixed.string (fixed)
# var.in.model(v,m)
# get.nmix (model, capthist, hcov)
# fixpmix(x, nmix)
# add.cl (df, alpha, loginterval, lowerbound = 0)
# spatialscale (object, detectfn, session = '')
#
# # logical for whether object specifies userDfn
# userD (object)
# # mean and SD if x numeric
# getMeanSD(xy)
# nclusters (capthist)
# leadingzero (x)
# HN (r, pars, cutval)
# HR (r, pars, cutval)
# EX (r, pars, cutval)
# UN (r, pars, cutval)
# CHN (r, pars, cutval)
# WEX (r, pars, cutval)
# ANN (r, pars, cutval)
# CLN (r, pars, cutval)
# CG (r, pars, cutval)
# CN (r, pars, cutval)
# BSS (r, pars, cutval)
# SS (r, pars, cutval)
# SSS (r, pars, cutval)
# SN (r, pars, cutval)
# SNS (r, pars, cutval)
# HHN (r, pars, cutval)
# HHR (r, pars, cutval)
# HEX (r, pars, cutval)
# HAN (r, pars, cutval)
# HCG (r, pars, cutval)
# HVP (r, pars, cutval)
# HPX (r, pars, cutval)
# gradient (pars, fun, eps=0.001, ...)
# transform (x, link)
# se.transform (real, sereal, link)
# untransform (beta, link)
# se.untransform (beta, sebeta, link)
# Xtransform (real, linkfn, varnames)
# se.Xtransform (real, sereal, linkfn, varnames)
# Xuntransform (beta, linkfn, varnames)
# se.Xuntransform (beta, sebeta, linkfn, varnames)
# mlogit.untransform (beta, latentmodel)
# clean.mlogit(x)
# mlogit (x)
# group.levels (capthist, groups, sep='.')
# h.levels (capthist, hcov, nmix)
# n.occasion (capthist)
# group.factor (capthist, groups, sep='.')
# getgrpnum (capthist, groups)
# make.lookup (tempmat)
# getknownclass(capthist, nmix, hcov)
# getnmix (details)
# fullbeta (beta, fb)
# inflate (xy, rmult = 1)
# getbinomN (binomN, detectr)
# xy2CH (CH, inflation = 1e-8)
# inflatechull (poly, r, ntheta = 60)
# maskarea (mask, sessnum = 1)
# masklength (mask, sessnum = 1)
# masksize (mask, sessnum = 1)
# complete.beta (object)
# complete.beta.vcv (object)
# smooths (formula)
# polys (formula)
# badsmooths (formula)
# gamsetup(formula, data, ...)
# general.model.matrix (formula, data, gamsmth = NULL,
# makerealparameters (design, beta, parindx, link, fixed)
# secr.lpredictor (formula, newdata, indx, beta, field, beta.vcv=NULL,
# edist (xy1, xy2)
# nedist (xy1, xy2, mask, inf = Inf, ...)
# getcellsize (mask)
# updatemodel (model, detectfn, detectfns, oldvar, newvar, warn = FALSE)
# deleteMaskPoints (mask, onebyone = TRUE, add = FALSE, poly = NULL,
# nparameters (object)
# mapbeta (parindx0, parindx1, beta0, betaindex)
# xyinpoly (xy, trps)
# addzerodf (df, oldCH, sess)
# addzeroCH (CH, nzero, cov = NULL)
# expandbinomN (binomN, detectorcodes)
# check3D (object)
# updateCH(object)
# newstr <-function (strings)
# outsidemask(CH, mask, threshold = spacing(mask) / sqrt(2))
# shareFactorLevels (object, columns = NULL, stringsAsFactors = TRUE)
# allzero (object)
# primarysessions(intervals)
# secondarysessions(intervals)
# boundarytoSF (poly)
# pointsInPolygon (xy, poly, logical = TRUE)
# firstsk (PIAx)
# maskboolean (ch, mask, threshold)
# telemcode(object, ...)
# uniquerownames (capthist)
# selectCHsession(capthist, sessnum)
# stringsAsFactors (DF)
# getdistmat2 (traps, mask, userdist, HPX = FALSE)
# uniformusage(object, noccasions)
# sfrotate (x, degrees, centrexy = NULL, usecentroid = FALSE)
# snap_points(x, y, max_dist = 1000)
# rtpois(n, lambda)
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.