Nothing
R/admbsecr.r
In admbsecr: SECR models with supplementary location information
#' Fitting SECR models in ADMB
#'
#' Fits an SECR model, with our without supplementary information
#' relevant to animal location. Parameter estimation is done by
#' maximum likelihood through an AD Model Builder (ADMB) executable.
#'
#' ADMB uses a quasi-Newton method to find maximum likelihood
#' estimates for the model parameters. Standard errors are calculated
#' by taking the inverse of the negative of the
#' Hessian. Alternatively, \link{boot.admbsecr} can be used to carry
#' out a parametric bootstrap procedure, from which parameter
#' uncertainty can also be inferred.
#'
#' The class of model fitted by this function (and, indeed, around
#' which this package is based) was first proposed by Borchers et
#' al. (in press); this reference is a good starting point for
#' practitioners looking to implement these methods.
#'
#' If the data are from an acoustic survey where individuals call more
#' than once (i.e., the argument \code{call.freqs} contains values
#' that are not 1), then standard errors calculated from the inverse
#' of the negative Hessian are not correct. They are therefore not
#' provided in this case. The method used by the function
#' \link{boot.admbsecr} is currently the only way to calculate these
#' reliably (see Stevenson et al., in prep., for details).
#'
#' @section The \code{capt} argument:
#' The \code{capt} argument is a list with named components. Each
#' component must be an \eqn{n} by \eqn{k} matrix, where \eqn{n} is
#' the number of detections made, and \eqn{k} is the number of traps
#' (or detectors) deployed. A component named \code{bincapt} is
#' compulsory.
#'
#' Further optional component names each refer to a type of
#' information which is informative on animal location collected on
#' each detection. Possible choices are: \code{bearing}, \code{dist},
#' \code{ss}, \code{toa}, and \code{mrds}.
#'
#' If the \eqn{i}th individual evaded the \eqn{j}th trap (or
#' detector), then the \eqn{j}th element in the \eqn{i}th row should
#' be 0 for all components. Otherwise, if the \eqn{i}th individual was
#' trapped (or detected) by the \eqn{j}th trap (or detector), then:
#' \itemize{
#' \item For the \code{bincapt} component, the element should be 1.
#' \item For the \code{bearing} component, the element should be the
#' estimated bearing from which the detector detected the
#' individual.
#' \item For the \code{dist} component, the element should be the
#' estimated distance between the individual and the detector
#' at the time of the detection.
#' \item For the \code{ss} component, the element should be the
#' measured signal strength of an acoustic signal detected by
#' the detector (only possible when the detectors are
#' microphones).
#' \item For the \code{toa} component, the element should be the
#' measured time of arrival (in seconds) since the start of
#' the survey (or some other reference time) of an acoustic
#' signal detected by the detector (only possible when the
#' detectors are microphones).
#' \item For the \code{mrds} component, the element should be the
#' \emph{known} (not estimated) distance between the individual
#' and the detector at the time of the detection.
#' }
#'
#' @section Fitted parameters:
#'
#' The parameter \code{D}, the density of individuals (or, in an
#' acoustic survey, the density of calls) is always fitted. The
#' effective survey area, \code{esa}, (see Borchers, 2012, for
#' details) is always provided as a derived parameter, with a standard
#' error calculated using the delta method.
#'
#' Further parameters to be fitted depend on the choice of the
#' detection function (i.e., the \code{detfn} argument), and the types
#' of additional information collected (i.e., the components in the
#' \code{capt}).
#'
#' Details of the detection functions are as follows:
#'
#' For \code{detfn = "hn"}:
#' \itemize{
#' \item Estimated paramters are \code{g0} and \code{sigma}.
#' \item \eqn{g(d) = g_0\ exp(-d^2/(2\sigma^2))}{g(d) = g0 * exp( -d^2 / (2 * sigma^2 ))}
#' }
#'
#' For \code{detfn = "hr"}:
#' \itemize{
#' \item Estimated parameters are \code{g0}, \code{sigma}, and
#' \code{z}.
#' \item \eqn{g(d) = g_0\ (1 - exp(-(d/\sigma)^{-z}))}{g(d) = g0 * ( 1 - exp( -(d/sigma)^{-z} ) )}
#' }
#'
#' For \code{detfn = "lth"}:
#' \itemize{
#' \item Estimated parameters are \code{shape.1}
#' \ifelse{latex}{(\eqn{\kappa})}{}, \code{shape.2}
#' \ifelse{latex}{(\eqn{\nu})}{}, and \code{scale}
#' \ifelse{latex}{(\eqn{\tau})}{}.
#' \item \eqn{g(d) = 0.5 - 0.5\ erf(\kappa - exp(\nu - \tau d))}{g(d) = 0.5 - 0.5 * erf( shape.1 - exp( shape.2 - scale * d ) )}
#' }
#'
#' For \code{detfn = "th"}:
#' \itemize{
#' \item Estimated parameters are \code{shape}
#' \ifelse{latex}{(\eqn{\kappa})}{} and \code{scale}
#' \ifelse{latex}{(\eqn{\tau})}{}.
#' \item \eqn{g(d) = 0.5 - 0.5\ erf(d/\tau - \kappa)}{g(d) = 0.5 - 0.5 * erf( d/scale - shape )}
#' }
#'
#' For \code{detfn = "ss"}:
#' \itemize{
#' \item The signal strength detection function is special in that
#' it requires signal strength information to be collected in
#' order for all parameters to be estimated.
#' \item Estimated parameters are \code{b0.ss}, \code{b1.ss}, and
#' \code{sigma.ss}.
#' \item The expected signal strength is modelled as:
#' \eqn{E(SS) = h^{-1}(\beta_0 - \beta_1d)}{E(SS) = h^{-1}(b0.ss - b1.ss*d)},
#' where \eqn{h} is specified by the argument \code{ss.link}.
#' }
#'
#' Details of the parameters associated with different additional data
#' types are as follows:
#'
#' For data type \code{"bearing"}, \code{kappa} is estimated. This is
#' the concerntration parameter of the von-Mises distribution used for
#' measurement error in estimated bearings.
#'
#' For data type \code{"dist"}, \code{alpha} is estimated. This is the
#' shape parameter of the gamma distribution used for measurement
#' error in estimated distances.
#'
#' For data type \code{"toa"}, \code{sigma.toa} is estimated. This is
#' the standard deviation parameter of the normal distribution used
#' for measurement error in recorded times of arrival.
#'
#' For data type \code{"mrds"}, no extra parameters are
#' estimated. Animal location is assumed to be known.
#'
#' @section Convergence:
#'
#' If maximum likelihood estimates could not be found during
#' optimisation, then \code{admbsecr} will usually show a warning that
#' the maximum gradient component is large, or possibly throw an error
#' reporting that a \code{.par} file is missing.
#'
#' The best approach to fixing convergence issues is to re-run the
#' \code{admbsecr} function with the argument \code{trace} set to
#' \code{TRUE}. Parameter values will be printed out for each step of
#' the optimisation algorithm.
#'
#' First, look for a large jump in a parameter to a value far from
#' what is feasible. This issue can be fixed by using the
#' \code{bounds} argument to restrict the parameter space over which
#' ADMB searches for the maximum likelihood estimate.
#'
#' Alternatively, try a different set of start values using the
#' argument \code{sv}; by default \code{admbsecr} will choose some
#' start values, but these are not necessarily sensible. The start
#' values that were used appear as the first line of text when
#' \code{trace} is \code{TRUE}.
#'
#' Sometimes the algorithm appears to converge, but nevertheless
#' perseveres reporting the same parameter values again and again for
#' a while (prior to the calculation of the Hessian). This is because
#' ADMB has failed to detect convergence as at least one of the
#' gradient components is still larger than the convergence criterion
#' (by default, 0.0001). It is possible to speed things up and help
#' ADMB detect convergence earlier by either tightening parameter
#' bounds (as above), or by setting appropriate scalefactors (using
#' the argument \code{sf}). To do this, first identify which
#' parameters have large gradient components from the "final
#' statistics" section of the \code{trace} output. Next, find the
#' default settings of the scalefactors by printing the object
#' \code{fit$args$sf}, where \code{fit} is the original object returned
#' by \code{admbsecr}. Finally, rerun \code{admbsecr} again, but this
#' time set the argument \code{sf} manually. Set scalefactors for any
#' parameters with small gradient components to the same as the
#' defaults ascertained above, and increase those associated with
#' large gradient components by a factor of 10. If the problem
#' persists, repeat this process (e.g., if the same parameters still
#' have large gradient components, increase the associated
#' scalefactors by another factor of 10).
#'
#' @section Local integration:
#'
#' For SECR models, the likelihood is calculated by integrating over
#' the unobserved animal activity centres (see Borchers & Efford,
#' 2008). Here, the integral is approximated numerically by taking a
#' finite sum over the mask points. The integrand is negligible in
#' size for mask points far from detectors that detected a particular
#' individual, and so to increase computational efficiency the region
#' over which this sum takes place can be reduced.
#'
#' Setting \code{local} to \code{TRUE} will only carry out this sum
#' across mask points that are within the mask buffer distance of
#' \emph{all} detectors that made a detection. So long as the buffer
#' suitably represents a distance beyond which detection is
#' practically impossible, the effect this has on parameter estimates
#' is negligible, but processing time can be substantially reduced.
#'
#' Note that this increases the parameter estimates' sensitivity to
#' the buffer. A buffer that is too small will lead to inaccurate
#' results.
#'
#' @references Borchers, D. L., and Efford, M. G. (2008) Spatially
#' explicit maximum likelihood methods for capture-recapture
#' studies. \emph{Biometrics}, \strong{64}: 377--385.
#'
#' @references Borchers, D. L. (2012) A non-technical overview of
#' spatially explicit capture-recapture models. \emph{Journal of
#' Ornithology}, \strong{152}: 435--444.
#'
#' @references Borchers, D. L., Stevenson, B. C., Kidney, D., Thomas,
#' L., and Marques, T. A. (in press) A unifying model for
#' capture-recapture and distance sampling surveys of wildlife
#' populations. \emph{Journal of the American Statistical
#' Association}.
#'
#' @references Stevenson, B. C., Borchers, D. L., Altwegg, R., Swift,
#' R. J., Gillespie, D. M., and Measey, G. J. (submitted) A general
#' framework for animal density estimation from acoustic detections
#' across a fixed microphone array.
#'
#' @return A list of class \code{"admbsecr"}. Components contain
#' information such as estimated parameters and standard errors. The
#' best way to access such information, however, is through the
#' variety of helper functions provided by the admbsecr package.
#'
#' @param capt A list with named components, containing the capture
#' history and supplementary information. The function
#' \link{create.capt} will return a suitable object. See 'Details'
#' below.
#' @param traps A matrix with two columns. Each row provides Cartesian
#' coordinates for the location of a trap (or detector).
#' @param mask A matrix with two columns. Each row provides Cartesian
#' coordinates for the location of a mask point. The function
#' \link{create.mask} will return a suitable object.
#' @param detfn A character string specifying the detection function
#' to be used. Options are "hn" (halfnormal), "hr" (hazard rate), "th"
#' (threshold), "lth" (log-link threshold), or "ss" (signal
#' strength). If the latter is used, signal strength information must
#' be provided in \code{capt}.
#' @param sv A named list. Component names are parameter names, and
#' each component is a start value for the associated parameter. See
#' 'Details' for further information on the parameters to be fitted.
#' @param bounds A named list. Component names are parameter names,
#' and each components is a vector of length two, specifying the
#' bounds for the associated parameter.
#' @param fix A named list. Component names are parameter names to be
#' fixed, and each component is the fixed value for the associated
#' parameter.
#' @param sf A named list. Component names are parameter names, and
#' each component is a scalefactor for the associated parameter. The
#' default behaviour is to automatically select scalefactors based on
#' parameter start values. See the section on convergence below.
#' @param ss.link A character string, either \code{"identity"} or
#' \code{"log"}, which specifies the link function for the signal
#' strength detection function. Only required when \code{detfn} is
#' \code{"ss"} (i.e., when there is signal strength information in
#' \code{capt}).
#' @param cutoff The signal strength threshold, above which sounds are
#' identified as detections. Only required when \code{detfn} is
#' \code{"ss"}.
#' @param call.freqs A vector of call frequencies collected
#' independently to an acoustic survey.
#' @param sound.speed The speed of sound in metres per second,
#' defaults to 330 (the speed of sound in air). Only used when
#' \code{"toa"} is a component name of \code{capt}.
#' @param local Logical, if \code{TRUE} integration over unobserved
#' animal activity centres is only carried out in a region local to
#' detectors that detected individuals. See 'Details'.
#' @param hess Logical, if \code{TRUE} the Hessian is estimated,
#' allowing for calculation of standard errors, the
#' variance-covariance matrix, and the correlation matrix, at the
#' expense of a little processing time. If \code{FALSE}, the Hessian
#' is not estimated. Note that if individuals are detectable more than
#' once (e.g., by calling more than once on an acoustic survey) then
#' parameter uncertainty is not properly represented by these
#' calculations.
#' @param trace Logical, if \code{TRUE} parameter values at each step
#' of the optimisation algorithm are printed to the R console.
#' @param clean Logical, if \code{TRUE} ADMB output files are
#' removed. Otherwise, ADMB output file will remain in a directory,
#' the location of which is reported after the model is fitted.
#' @param cbs The CMPDIF_BUFFER_SIZE, set using the \code{-cbs} option
#' of the executable created by ADMB. This can be increased to speed
#' up optimisation if \code{cmpdiff.tmp} gets too large (please
#' ignore, unless you are familiar with ADMB and know what you are
#' doing).
#' @param gbs The GRADSTACK_BUFFER_SIZE, set using the \code{-gbs}
#' option of the executable created by ADMB. This can be increased to
#' speed up optimisation if \code{gradfil1.tmp} gets too large (please
#' ignore, unless you are familiar with ADMB and know what you are
#' doing).
#' @param exe.type Character string, either \code{"old"} or
#' \code{"new"}, depending on which executable is to be used (for
#' development purposes only; please ignore).
#'
#' @seealso \link{boot.admbsecr} to calculate standard errors and
#' estimate bias using a parametric bootstrap.
#' @seealso \link{coef.admbsecr}, \link{stdEr.admbsecr}, and
#' \link{vcov.admbsecr} to extract estimated parameters, standard
#' errors, and the variance-covariance matrix, respectively.
#' @seealso \link{confint.admbsecr} to calculate confidence intervals.
#' @seealso \link{summary.admbsecr} to get a summary of estimates and
#' standard errors.
#' @seealso \link{show.detfn} to plot the estimated detection
#' function.
#' @seealso \link{locations} to plot estimated locations of particular
#' individuals or calls.
#'
#' @examples
#' \dontrun{
#' simple.capt <- example$capt["bincapt"]
#' simple.hn.fit <- admbsecr(capt = simple.capt, traps = example$traps,
#' mask = example$mask, fix = list(g0 = 1))
#' simple.hr.fit <- admbsecr(capt = simple.capt, traps = example$traps,
#' mask = example$mask, detfn = "hr")
#' bearing.capt <- example.capt[c("bincapt", "bearing")]
#' bearing.hn.fit <- admbsecr(capt = bearing.capt, traps = example$traps,
#' mask = example$mask, fix = list(g0 = 1))
#' }
#'
#' @export
#'
admbsecr <- function(capt, traps, mask, detfn = "hn", sv = NULL, bounds = NULL,
fix = NULL, sf = NULL, ss.link = "identity",
cutoff = NULL, call.freqs = NULL, sound.speed = 330,
local = FALSE, hess = !any(call.freqs > 1), trace = FALSE,
clean = TRUE, cbs = NULL, gbs = NULL, exe.type = "old"){
arg.names <- names(as.list(environment()))
capt.bin <- capt$bincapt
## Checking for bincapt.
if (is.null(capt.bin)){
stop("The binary capture history must be provided as a component of 'capt'.")
}
## Checking for correct number of trap locations.
if (ncol(capt.bin) != nrow(traps)){
stop("There must be a trap location for each column in the components of 'capt'.")
}
## Checking that each component of 'capt' is a matrix.
if (any(!laply(capt, is.matrix))){
stop("At least one component of 'capt' is not a matrix.")
}
## Checking for agreement in matrix dimensions.
if (length(capt) > 1){
all.dims <- laply(capt, dim)
if (any(aaply(all.dims, 2, function(x) diff(range(x))) != 0)){
stop("Components of 'capt' object have different dimensions.")
}
}
## Various checks for other arguments.
if (!is.list(sv) & !is.null(sv)){
stop("The 'sv' argument must be 'NULL' or a list.")
}
if (!is.list(bounds) & !is.null(bounds)){
stop("The 'bounds' argument must be 'NULL' or a list.")
}
if (is.list(bounds)){
if (any(laply(bounds, length) != 2)){
stop("Each component of 'bounds' must be a vector of length 2.")
}
}
if (!is.list(fix) & !is.null(fix)){
stop("The 'fix' argument must be 'NULL' or a list.")
}
n <- nrow(capt.bin)
n.traps <- nrow(traps)
n.mask <- nrow(mask)
A <- attr(mask, "area")
buffer <- attr(mask, "buffer")
## Removing attributes from mask.
mask <- as.matrix(mask)
attr(mask, "area") <- A
attr(mask, "buffer") <- buffer
## TODO: Sort out how to determine supplementary parameter names.
supp.types <- c("bearing", "dist", "ss", "toa", "mrds")
fit.types <- supp.types %in% names(capt)
names(fit.types) <- supp.types
## Logical indicators for additional information types.
fit.bearings <- fit.types["bearing"]
fit.dists <- fit.types["dist"]
fit.ss <- fit.types["ss"]
fit.toas <- fit.types["toa"]
fit.mrds <- fit.types["mrds"]
## Generating ordered binary capture history.
capt.bin.order <- do.call(order, as.data.frame(capt.bin))
capt.bin.unique <- capt.bin[capt.bin.order, ]
capt.bin.freqs <- as.vector(table(apply(capt.bin.unique, 1, paste, collapse = "")))
names(capt.bin.freqs) <- NULL
capt.bin.unique <- capt.bin.unique[!duplicated(as.data.frame(capt.bin.unique)), ]
n.unique <- nrow(capt.bin.unique)
unique.changes <- cumsum(c(0, capt.bin.freqs[-n.unique])) + 1
## Reordering all capture history components.
capt.ord <- capt
for (i in 1:length(capt)){
capt.ord[[i]] <- capt[[i]][capt.bin.order, ]
}
## Capture histories for additional information types (if they exist)
capt.bearing <- if (fit.bearings) capt.ord$bearing else 0
capt.dist <- if (fit.dists) capt.ord$dist else 0
capt.ss <- if (fit.ss) capt.ord$ss else 0
capt.toa <- if (fit.toas) capt.ord$toa else 0
mrds.dist <- if (fit.mrds) capt.ord$mrds else 0
suppar.names <- c("kappa", "alpha", "sigma.toa")[fit.types[c("bearing", "dist", "toa")]]
if (fit.ss){
## Warning for failure to provide 'cutoff'.
if (missing(cutoff)){
stop("Argument 'cutoff' is missing.")
}
if (!missing(detfn) & detfn != "ss"){
warning("Argument 'detfn' is being ignored as signal strength information is provided in 'capt'. A signal strength detection function has been fitted instead.")
}
if (ss.link == "identity"){
detfn <- "ss"
linkfn.id <- 1
} else if (ss.link == "log"){
detfn <- "log.ss"
linkfn.id <- 2
} else {
stop("ss.link must be either \"identity\" or \"log\"")
}
} else {
## Not sure what a linkfn.id of 3 means? Probably throws an error in ADMB.
linkfn.id <- 3
}
detfns <- c("hn", "hr", "th", "lth", "ss", "log.ss")
## Sets detection function ID number for use in ADMB:
## 1 = Half normal
## 2 = Hazard rate
## 3 = Threshold
## 4 = Log-link threshold
## 5 = Identity-link signal strength
## 6 = Log-link signal strength.
detfn.id <- which(detfn == detfns)
detpar.names <- switch(detfn,
hn = c("g0", "sigma"),
hr = c("g0", "sigma", "z"),
th = c("shape", "scale"),
lth = c("shape.1", "shape.2", "scale"),
ss = c("b0.ss", "b1.ss", "sigma.ss"),
log.ss = c("b0.ss", "b1.ss", "sigma.ss"))
par.names <- c("D", detpar.names, suppar.names)
n.detpars <- length(detpar.names)
n.suppars <- length(suppar.names)
any.suppars <- n.suppars > 0
n.pars <- length(par.names)
## Checking par.names against names of sv, fix, bounds, and sf.
for (i in c("sv", "fix", "bounds", "sf")){
obj <- get(i)
if (!is.null(obj)){
obj.fitted <- names(obj) %in% par.names
if(!all(obj.fitted)){
msg <- paste("Some parameters listed in '", i, "' are not being used. These are being removed.",
sep = "")
warning(msg)
assign(i, obj[obj.fitted])
}
}
}
## Sets link function ID number for use in ADMB:
## 1 = identity
## 2 = log
## 3 = logit
links <- list(D = 2,
g0 = 3,
sigma = 2,
shape = 1,
shape.1 = 2,
shape.2 = 1,
scale = 2,
b0.ss = 2,
b1.ss = 2,
sigma.ss = 2,
z = 2,
sigma.toa = 2,
kappa = 2,
alpha = 2)[par.names]
link.list <- list(identity, log.link, logit.link)
unlink.list <- list(identity, exp, inv.logit)
par.links <- llply(links, function(x, link.list) link.list[[x]], link.list)
par.unlinks <- llply(links, function(x, unlink.list) unlink.list[[x]], unlink.list)
## Sorting out start values. Start values are set to those provided,
## or else are determined automatically from functions in
## autofuns.r.
sv.link <- vector("list", length = n.pars)
names(sv.link) <- par.names
sv.link[names(sv)] <- sv
sv.link[names(fix)] <- fix
auto.names <- par.names[sapply(sv.link, is.null)]
sv.funs <- paste("auto", auto.names, sep = "")
## Done in reverse so that D is calculated last (requires detfn parameters).
## D not moved to front as it should appear as the first parameter in any output.
for (i in rev(seq(1, length(auto.names), length.out = length(auto.names)))){
sv.link[auto.names[i]] <- eval(call(sv.funs[i],
list(capt = capt, detfn = detfn,
detpar.names = detpar.names,
mask = mask, traps = traps,
sv = sv.link, cutoff = cutoff,
ss.link = ss.link, A = A)))
}
## Converting start values to link scale.
sv <- sv.link
for (i in names(sv.link)){
sv.link[[i]] <- link.list[[links[[i]]]](sv.link[[i]])
}
## Sorting out phases.
## TODO: Add phases parameter so that these can be controlled by user.
phases <- vector("list", length = n.pars)
names(phases) <- par.names
for (i in par.names){
if (any(i == names(fix))){
## Phase of -1 in ADMB fixes parameter at starting value.
phases[[i]] <- -1
} else {
phases[[i]] <- 0
}
}
D.phase <- phases[["D"]]
detpars.phase <- c(phases[detpar.names], recursive = TRUE)
if (any.suppars){
suppars.phase <- c(phases[suppar.names], recursive = TRUE)
} else {
suppars.phase <- -1
}
## Sorting out bounds.
## Below bounds are the defaults.
default.bounds <- list(D = c(n/(A*n.mask), 1e8),
g0 = c(0, 1),
sigma = c(0, 1e8),
shape = c(-100, 100),
shape.1 = c(0, 1e8),
shape.2 = c(-100, 100),
scale = c(0, 1e8),
b0.ss = c(0, 1e8),
b1.ss = c(0, 1e8),
sigma.ss = c(0, 1e8),
z = c(0, 1e8),
sigma.toa = c(0, 1e8),
kappa = c(0, 700),
alpha = c(0, 1e8))[par.names]
bound.changes <- bounds
bounds <- default.bounds
for (i in names(default.bounds)){
if (i %in% names(bound.changes)){
bounds[[i]] <- bound.changes[[i]]
}
}
## Converting bounds to link scale.
bounds.link <- bounds
for (i in names(bounds)){
bounds.link[[i]] <- link.list[[links[[i]]]](bounds[[i]])
}
D.bounds <- bounds.link[["D"]]
D.lb <- D.bounds[1]
D.ub <- D.bounds[2]
detpar.bounds <- bounds.link[detpar.names]
detpars.lb <- sapply(detpar.bounds, function(x) x[1])
detpars.ub <- sapply(detpar.bounds, function(x) x[2])
if (any.suppars){
suppar.bounds <- bounds.link[suppar.names]
suppars.lb <- sapply(suppar.bounds, function(x) x[1])
suppars.ub <- sapply(suppar.bounds, function(x) x[2])
} else {
suppars.lb <- 0
suppars.ub <- 0
}
## Sorting out scalefactors.
if (is.null(sf)){
sv.vec <- c(sv.link, recursive = TRUE)
## Currently, by default, the scalefactors are the inverse
## fraction of each starting value to the largest starting
## value. Not sure how sensible this is.
##bound.ranges <- laply(bounds.link, function(x) diff(range(x)))
##sf <- max(bound.ranges)/bound.ranges
sf <- abs(max(sv.vec)/sv.vec)
names(sf) <- par.names
} else if (is.list(sf)){
sf <- numeric(n.pars)
names(sf) <- par.names
for (i in par.names){
sf[i] <- ifelse(i %in% names(sf), sf[[i]], 1)
}
} else if (is.vector(sf) & length(sf) == 1){
sf <- rep(sf, length(par.names))
names(sf) <- par.names
}
## Replacing infinite scalefactors.
sf[!is.finite(sf)] <- 1
D.sf <- sf[["D"]]
detpars.sf <- c(sf[detpar.names], recursive = TRUE)
if (any.suppars){
suppars.sf <- c(sf[suppar.names], recursive = TRUE)
} else {
suppars.sf <- 1
}
## Creating link objects to pass to ADMB.
detpars.link <- c(links[detpar.names], recursive = TRUE)
if (any.suppars){
suppars.link <- c(links[suppar.names], recursive = TRUE)
} else {
suppars.link <- 1
}
## Setting small number so that numerical under/overflow in ADMB
## does not affect estimation.
dbl.min <- 1e-150
## Calculating distances and angles.
dists <- distances(traps, mask)
if (fit.bearings){
bearings <- bearings(traps, mask)
} else {
bearings <- 0
}
if (fit.toas){
toa.ssq <- make_toa_ssq(capt.ord$toa, dists, sound.speed)
} else {
toa.ssq <- 0
}
if (is.null(cutoff)){
cutoff <- 0
}
## Kludge to fix number of parameters for no supplementary
## information.
if (!any.suppars){
n.suppars <- max(c(n.suppars, 1))
sv.link$dummy <- 0
}
## Sorting out which mask points are local to each detection.
if (local){
all.which.local <- find_local(capt.bin.unique, dists, buffer)
all.n.local <- laply(all.which.local, length)
all.which.local <- c(all.which.local, recursive = TRUE)
} else {
all.n.local <- rep(1, n.unique)
all.which.local <- rep(0, n.unique)
}
## Stuff for the .dat file.
data.list <- list(
n_unique = n.unique, local = as.numeric(local), all_n_local = all.n.local,
all_which_local = all.which.local, D_lb = D.lb, D_ub = D.ub, D_phase =
D.phase, D_sf = D.sf, n_detpars = n.detpars, detpars_lb = detpars.lb,
detpars_ub = detpars.ub, detpars_phase = detpars.phase, detpars_sf =
detpars.sf, detpars_linkfns = detpars.link, n_suppars = n.suppars,
suppars_lb = suppars.lb, suppars_ub = suppars.ub, suppars_phase =
suppars.phase, suppars_sf = suppars.sf, suppars_linkfns =
suppars.link, detfn_id = detfn.id, buffer = buffer, trace =
as.numeric(trace), DBL_MIN = dbl.min, n = n, n_traps = n.traps, n_mask
= n.mask, A = A, capt_bin_unique = capt.bin.unique, capt_bin_freqs =
capt.bin.freqs, fit_angs = as.numeric(fit.bearings), capt_ang =
capt.bearing, fit_dists = as.numeric(fit.dists), capt_dist =
capt.dist, fit_ss = as.numeric(fit.ss), cutoff = cutoff, linkfn_id =
linkfn.id, capt_ss = capt.ss, fit_toas = as.numeric(fit.toas),
capt_toa = capt.toa, fit_mrds = as.numeric(fit.mrds), mrds_dist =
mrds.dist, dists = dists, angs = bearings, toa_ssq = toa.ssq)
## Determining whether or not standard errors should be calculated.
if (!is.null(call.freqs)){
fit.freqs <- any(call.freqs != 1)
} else {
fit.freqs <- FALSE
}
## Idea of running executable as below taken from glmmADMB.
## Working out correct command to run from command line.
if (exe.type == "new"){
exe.name <- "secr_new"
} else if (exe.type == "old"){
exe.name <- "secr"
} else if (exe.type == "test"){
exe.name <- "secr_test"
} else {
stop("Argument 'exe.type' must be \"old\" or \"new\".")
}
prefix.name <- exe.name
if (.Platform$OS == "windows"){
os.type <- "windows"
exe.name <- paste(prefix.name, ".exe", sep = "")
} else if (.Platform$OS == "unix"){
if (Sys.info()["sysname"] == "Linux"){
os.type <- "linux"
} else if (Sys.info()["sysname"] == "Darwin"){
os.type <- "mac"
} else {
stop("Unknown OS type.")
}
} else {
stop("Unknown OS type.")
}
## Finding executable folder (possible permission problems?).
exe.dir <- paste(system.file(package = "admbsecr"), "ADMB", "bin", os.type, sep = "/")
exe.loc <- paste(exe.dir, exe.name, sep = "/")
## Creating command to run using system().
curr.dir <- getwd()
## Creating temporary directory.
temp.dir <- tempfile("admbsecr", curr.dir)
dir.create(temp.dir)
setwd(temp.dir)
## Creating .pin and .dat files.
write_pin("secr", sv.link)
write_dat("secr", data.list)
## Creating link to executable.
if (os.type == "windows"){
file.copy(exe.loc, exe.name)
} else {
file.symlink(exe.loc, exe.name)
}
## Sorting out -cbs and -gbs.
if (!is.null(cbs)){
cbs.cmd <- paste(" -cbs", format(cbs, scientific = FALSE))
} else {
cbs.cmd <- NULL
}
if (!is.null(gbs)){
gbs.cmd <- paste(" -gbs", format(gbs, scientific = FALSE))
} else {
gbs.cmd <- NULL
}
## Running ADMB executable.
cmd <- paste("./"[os.type != "windows"], exe.name,
" -ind secr.dat -ainp secr.pin",
" -nohess"[!hess], cbs.cmd, gbs.cmd, sep = "")
if (os.type == "windows"){
system(cmd, ignore.stdout = !trace, show.output.on.console = trace)
} else {
system(cmd, ignore.stdout = !trace)
}
## Reading in model results.
options(warn = -1)
if (exe.type == "test"){
prefix.name <- strsplit(list.files(), "\\.")[[which(substr(list.files(),
nchar(list.files()) - 3,
nchar(list.files())) == ".par")]][1]
}
out <- read.admbsecr(prefix.name)
options(warn = 0)
setwd(curr.dir)
## Cleaning up files.
if (clean){
unlink(temp.dir, recursive = TRUE)
} else {
cat("ADMB files found in:", "\n", temp.dir, "\n")
}
## Warning for non-convergence.
if (out$maxgrad < -0.1){
warning("Failed convergence -- maximum gradient component is large.")
}
## Moving back to original working directory.
setwd(curr.dir)
## Removing fixed coefficients from list.
if (!hess){
out$coeflist[c(D.phase, detpars.phase, suppars.phase) == -1] <- NULL
}
## Creating coefficients vector.
est.pars <- c("D", detpar.names, suppar.names)[c(D.phase, detpars.phase, suppars.phase) > -1]
n.est.pars <- length(est.pars)
out$coefficients <- numeric(2*n.est.pars + 1)
names(out$coefficients) <- c(paste(est.pars, "_link", sep = ""), est.pars, "esa")
for (i in 1:n.est.pars){
out$coefficients[i] <- out$coeflist[[i]]
}
for (i in 1:n.est.pars){
out$coefficients[n.est.pars + i] <-
unlink.list[[links[[est.pars[i]]]]](out$coeflist[[i]])
}
## Adding extra components to list.
if (detfn == "log.ss") detfn <- "ss"
## Putting in updated argument names.
args <- vector(mode = "list", length = length(arg.names))
names(args) <- arg.names
for (i in arg.names){
if (!is.null(get(i))){
args[[i]] <- get(i)
}
}
out$args <- args
out$fit.types <- fit.types
out$infotypes <- names(fit.types)[fit.types]
out$detpars <- detpar.names
out$suppars <- suppar.names
out$phases <- phases
out$par.links <- par.links
out$par.unlinks <- par.unlinks
## Putting in esa estimate.
out$coefficients[2*n.est.pars + 1] <- p.dot(out, esa = TRUE)
## Putting in call frequency information and correct parameter names.
if (fit.freqs){
mu.freqs <- mean(call.freqs)
Da <- get.par(out, "D")/mu.freqs
names.vec <- c(names(out[["coefficients"]]), "Da", "mu.freqs")
coefs.updated <- c(out[["coefficients"]], Da, mu.freqs)
names(coefs.updated) <- names.vec
out[["coefficients"]] <- coefs.updated
## Removing ses, cor, vcov matrices.
cor.updated <- matrix(NA, nrow = length(names.vec),
ncol = length(names.vec))
dimnames(cor.updated) <- list(names.vec, names.vec)
vcov.updated <- matrix(NA, nrow = length(names.vec),
ncol = length(names.vec))
dimnames(vcov.updated) <- list(names.vec, names.vec)
if (hess){
ses.updated <- c(out[["se"]], rep(NA, 2))
max.ind <- length(names.vec) - 2
cor.updated[1:max.ind, 1:max.ind] <- out[["cor"]]
vcov.updated[1:max.ind, 1:max.ind] <- out[["vcov"]]
} else {
ses.updated <- rep(NA, length(names.vec))
}
names(ses.updated) <- names.vec
out[["se"]] <- ses.updated
out[["cor"]] <- cor.updated
out[["vcov"]] <- vcov.updated
if (trace){
if (!hess){
cat("NOTE: Standard errors not calculated; use boot.admbsecr().", "\n")
} else {
cat("NOTE: Standard errors are probably not correct; use boot.admbsecr().", "\n")
}
}
} else {
if (hess){
## Putting correct parameter names into se, cor, vcov.
replace <- substr(names(out$se), 1, 8) == "par_ests"
names(out$se)[replace] <- rownames(out$vcov)[replace] <-
colnames(out$vcov)[replace] <- rownames(out$cor)[replace] <-
colnames(out$cor)[replace] <- est.pars
replace <- 1:length(est.pars)
names(out$se)[replace] <- rownames(out$vcov)[replace] <-
colnames(out$vcov)[replace] <- rownames(out$cor)[replace] <-
colnames(out$cor)[replace] <- paste(est.pars, "_link", sep = "")
} else {
## Filling se, cor, vcov with NAs.
names.vec <- names(out[["coefficients"]])
ses.updated <- rep(NA, length(names.vec))
names(ses.updated) <- names.vec
cor.updated <- matrix(NA, nrow = length(names.vec),
ncol = length(names.vec))
dimnames(cor.updated) <- list(names.vec, names.vec)
vcov.updated <- matrix(NA, nrow = length(names.vec),
ncol = length(names.vec))
dimnames(vcov.updated) <- list(names.vec, names.vec)
out[["se"]] <- ses.updated
out[["cor"]] <- cor.updated
out[["vcov"]] <- vcov.updated
}
}
out$fit.freqs <- fit.freqs
if (out$maxgrad < -0.01){
warning("Maximum gradient component is large.")
}
class(out) <- c("admbsecr", "admb")
out
}
## Roxygen code for NAMESPACE and datasets.
## Package imports for roxygenise to pass to NAMESPACE.
#' @import plyr Rcpp R2admb
#' @importFrom CircStats dvm rvm
#' @importFrom lattice wireframe
#' @importFrom matrixStats colProds
#' @importFrom secr make.capthist make.mask read.mask read.traps sim.popn
#' @useDynLib admbsecr
NULL
## Data documentation.
#' \emph{Arthroleptella lightfooti} survey data
#'
#' Data from an acoustic survey of the Western Cape moss frog
#' \emph{Arthroleptella lightfooti}. These data are from a 25 s subset
#' of the original recording, taken on 16 May 2012 at Silvermine,
#' Table Mountain National Park, Cape Town, South Africa. Acoustic
#' signal strengths and times of arrival were measured, and this
#' information is contained in the capture history object.
#'
#' This object is a list which contains components:
#' \itemize{
#' \item \code{capt}: A capture history object.
#' \item \code{traps}: A traps object.
#' \item \code{mask}: A suitable mask object.
#' \item \code{cutoff}: The microphone cutoff value.
#' \item \code{freqs}: A vector of call frequencies measured
#' independently to the acoustic survey.
#' }
#'
#' @name lightfooti
#' @format A list.
#' @usage lightfooti
#' @docType data
#' @keywords datasets
NULL
#' Example data
#'
#' This object contains simulated data with all types of supplementary
#' information, corresponding trap locations, and a suitable mask
#' object. Also included are some example model fits, which were
#' generated from these data using the \link{admbsecr} function.
#'
#' This object is a list which contains components:
#' \itemize{
#' \item \code{capt}: A capture history object.
#' \item \code{traps}: A traps object.
#' \item \code{mask}: A suitable mask object.
#' \item \code{cutoff}: The cutoff value used to simluate these data.
#' \item \code{fits}: Some example model fits.
#' }
#'
#' @name example
#' @format A list.
#' @usage example
#' @docType data
#' @keywords datasets
NULL
Try the admbsecr package in your browser
Any scripts or data that you put into this service are public.
admbsecr documentation built on May 2, 2019, 5:21 p.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.
#' Fitting SECR models in ADMB
#'
#' Fits an SECR model, with our without supplementary information
#' relevant to animal location. Parameter estimation is done by
#' maximum likelihood through an AD Model Builder (ADMB) executable.
#'
#' ADMB uses a quasi-Newton method to find maximum likelihood
#' estimates for the model parameters. Standard errors are calculated
#' by taking the inverse of the negative of the
#' Hessian. Alternatively, \link{boot.admbsecr} can be used to carry
#' out a parametric bootstrap procedure, from which parameter
#' uncertainty can also be inferred.
#'
#' The class of model fitted by this function (and, indeed, around
#' which this package is based) was first proposed by Borchers et
#' al. (in press); this reference is a good starting point for
#' practitioners looking to implement these methods.
#'
#' If the data are from an acoustic survey where individuals call more
#' than once (i.e., the argument \code{call.freqs} contains values
#' that are not 1), then standard errors calculated from the inverse
#' of the negative Hessian are not correct. They are therefore not
#' provided in this case. The method used by the function
#' \link{boot.admbsecr} is currently the only way to calculate these
#' reliably (see Stevenson et al., in prep., for details).
#'
#' @section The \code{capt} argument:
#' The \code{capt} argument is a list with named components. Each
#' component must be an \eqn{n} by \eqn{k} matrix, where \eqn{n} is
#' the number of detections made, and \eqn{k} is the number of traps
#' (or detectors) deployed. A component named \code{bincapt} is
#' compulsory.
#'
#' Further optional component names each refer to a type of
#' information which is informative on animal location collected on
#' each detection. Possible choices are: \code{bearing}, \code{dist},
#' \code{ss}, \code{toa}, and \code{mrds}.
#'
#' If the \eqn{i}th individual evaded the \eqn{j}th trap (or
#' detector), then the \eqn{j}th element in the \eqn{i}th row should
#' be 0 for all components. Otherwise, if the \eqn{i}th individual was
#' trapped (or detected) by the \eqn{j}th trap (or detector), then:
#' \itemize{
#' \item For the \code{bincapt} component, the element should be 1.
#' \item For the \code{bearing} component, the element should be the
#' estimated bearing from which the detector detected the
#' individual.
#' \item For the \code{dist} component, the element should be the
#' estimated distance between the individual and the detector
#' at the time of the detection.
#' \item For the \code{ss} component, the element should be the
#' measured signal strength of an acoustic signal detected by
#' the detector (only possible when the detectors are
#' microphones).
#' \item For the \code{toa} component, the element should be the
#' measured time of arrival (in seconds) since the start of
#' the survey (or some other reference time) of an acoustic
#' signal detected by the detector (only possible when the
#' detectors are microphones).
#' \item For the \code{mrds} component, the element should be the
#' \emph{known} (not estimated) distance between the individual
#' and the detector at the time of the detection.
#' }
#'
#' @section Fitted parameters:
#'
#' The parameter \code{D}, the density of individuals (or, in an
#' acoustic survey, the density of calls) is always fitted. The
#' effective survey area, \code{esa}, (see Borchers, 2012, for
#' details) is always provided as a derived parameter, with a standard
#' error calculated using the delta method.
#'
#' Further parameters to be fitted depend on the choice of the
#' detection function (i.e., the \code{detfn} argument), and the types
#' of additional information collected (i.e., the components in the
#' \code{capt}).
#'
#' Details of the detection functions are as follows:
#'
#' For \code{detfn = "hn"}:
#' \itemize{
#' \item Estimated paramters are \code{g0} and \code{sigma}.
#' \item \eqn{g(d) = g_0\ exp(-d^2/(2\sigma^2))}{g(d) = g0 * exp( -d^2 / (2 * sigma^2 ))}
#' }
#'
#' For \code{detfn = "hr"}:
#' \itemize{
#' \item Estimated parameters are \code{g0}, \code{sigma}, and
#' \code{z}.
#' \item \eqn{g(d) = g_0\ (1 - exp(-(d/\sigma)^{-z}))}{g(d) = g0 * ( 1 - exp( -(d/sigma)^{-z} ) )}
#' }
#'
#' For \code{detfn = "lth"}:
#' \itemize{
#' \item Estimated parameters are \code{shape.1}
#' \ifelse{latex}{(\eqn{\kappa})}{}, \code{shape.2}
#' \ifelse{latex}{(\eqn{\nu})}{}, and \code{scale}
#' \ifelse{latex}{(\eqn{\tau})}{}.
#' \item \eqn{g(d) = 0.5 - 0.5\ erf(\kappa - exp(\nu - \tau d))}{g(d) = 0.5 - 0.5 * erf( shape.1 - exp( shape.2 - scale * d ) )}
#' }
#'
#' For \code{detfn = "th"}:
#' \itemize{
#' \item Estimated parameters are \code{shape}
#' \ifelse{latex}{(\eqn{\kappa})}{} and \code{scale}
#' \ifelse{latex}{(\eqn{\tau})}{}.
#' \item \eqn{g(d) = 0.5 - 0.5\ erf(d/\tau - \kappa)}{g(d) = 0.5 - 0.5 * erf( d/scale - shape )}
#' }
#'
#' For \code{detfn = "ss"}:
#' \itemize{
#' \item The signal strength detection function is special in that
#' it requires signal strength information to be collected in
#' order for all parameters to be estimated.
#' \item Estimated parameters are \code{b0.ss}, \code{b1.ss}, and
#' \code{sigma.ss}.
#' \item The expected signal strength is modelled as:
#' \eqn{E(SS) = h^{-1}(\beta_0 - \beta_1d)}{E(SS) = h^{-1}(b0.ss - b1.ss*d)},
#' where \eqn{h} is specified by the argument \code{ss.link}.
#' }
#'
#' Details of the parameters associated with different additional data
#' types are as follows:
#'
#' For data type \code{"bearing"}, \code{kappa} is estimated. This is
#' the concerntration parameter of the von-Mises distribution used for
#' measurement error in estimated bearings.
#'
#' For data type \code{"dist"}, \code{alpha} is estimated. This is the
#' shape parameter of the gamma distribution used for measurement
#' error in estimated distances.
#'
#' For data type \code{"toa"}, \code{sigma.toa} is estimated. This is
#' the standard deviation parameter of the normal distribution used
#' for measurement error in recorded times of arrival.
#'
#' For data type \code{"mrds"}, no extra parameters are
#' estimated. Animal location is assumed to be known.
#'
#' @section Convergence:
#'
#' If maximum likelihood estimates could not be found during
#' optimisation, then \code{admbsecr} will usually show a warning that
#' the maximum gradient component is large, or possibly throw an error
#' reporting that a \code{.par} file is missing.
#'
#' The best approach to fixing convergence issues is to re-run the
#' \code{admbsecr} function with the argument \code{trace} set to
#' \code{TRUE}. Parameter values will be printed out for each step of
#' the optimisation algorithm.
#'
#' First, look for a large jump in a parameter to a value far from
#' what is feasible. This issue can be fixed by using the
#' \code{bounds} argument to restrict the parameter space over which
#' ADMB searches for the maximum likelihood estimate.
#'
#' Alternatively, try a different set of start values using the
#' argument \code{sv}; by default \code{admbsecr} will choose some
#' start values, but these are not necessarily sensible. The start
#' values that were used appear as the first line of text when
#' \code{trace} is \code{TRUE}.
#'
#' Sometimes the algorithm appears to converge, but nevertheless
#' perseveres reporting the same parameter values again and again for
#' a while (prior to the calculation of the Hessian). This is because
#' ADMB has failed to detect convergence as at least one of the
#' gradient components is still larger than the convergence criterion
#' (by default, 0.0001). It is possible to speed things up and help
#' ADMB detect convergence earlier by either tightening parameter
#' bounds (as above), or by setting appropriate scalefactors (using
#' the argument \code{sf}). To do this, first identify which
#' parameters have large gradient components from the "final
#' statistics" section of the \code{trace} output. Next, find the
#' default settings of the scalefactors by printing the object
#' \code{fit$args$sf}, where \code{fit} is the original object returned
#' by \code{admbsecr}. Finally, rerun \code{admbsecr} again, but this
#' time set the argument \code{sf} manually. Set scalefactors for any
#' parameters with small gradient components to the same as the
#' defaults ascertained above, and increase those associated with
#' large gradient components by a factor of 10. If the problem
#' persists, repeat this process (e.g., if the same parameters still
#' have large gradient components, increase the associated
#' scalefactors by another factor of 10).
#'
#' @section Local integration:
#'
#' For SECR models, the likelihood is calculated by integrating over
#' the unobserved animal activity centres (see Borchers & Efford,
#' 2008). Here, the integral is approximated numerically by taking a
#' finite sum over the mask points. The integrand is negligible in
#' size for mask points far from detectors that detected a particular
#' individual, and so to increase computational efficiency the region
#' over which this sum takes place can be reduced.
#'
#' Setting \code{local} to \code{TRUE} will only carry out this sum
#' across mask points that are within the mask buffer distance of
#' \emph{all} detectors that made a detection. So long as the buffer
#' suitably represents a distance beyond which detection is
#' practically impossible, the effect this has on parameter estimates
#' is negligible, but processing time can be substantially reduced.
#'
#' Note that this increases the parameter estimates' sensitivity to
#' the buffer. A buffer that is too small will lead to inaccurate
#' results.
#'
#' @references Borchers, D. L., and Efford, M. G. (2008) Spatially
#' explicit maximum likelihood methods for capture-recapture
#' studies. \emph{Biometrics}, \strong{64}: 377--385.
#'
#' @references Borchers, D. L. (2012) A non-technical overview of
#' spatially explicit capture-recapture models. \emph{Journal of
#' Ornithology}, \strong{152}: 435--444.
#'
#' @references Borchers, D. L., Stevenson, B. C., Kidney, D., Thomas,
#' L., and Marques, T. A. (in press) A unifying model for
#' capture-recapture and distance sampling surveys of wildlife
#' populations. \emph{Journal of the American Statistical
#' Association}.
#'
#' @references Stevenson, B. C., Borchers, D. L., Altwegg, R., Swift,
#' R. J., Gillespie, D. M., and Measey, G. J. (submitted) A general
#' framework for animal density estimation from acoustic detections
#' across a fixed microphone array.
#'
#' @return A list of class \code{"admbsecr"}. Components contain
#' information such as estimated parameters and standard errors. The
#' best way to access such information, however, is through the
#' variety of helper functions provided by the admbsecr package.
#'
#' @param capt A list with named components, containing the capture
#' history and supplementary information. The function
#' \link{create.capt} will return a suitable object. See 'Details'
#' below.
#' @param traps A matrix with two columns. Each row provides Cartesian
#' coordinates for the location of a trap (or detector).
#' @param mask A matrix with two columns. Each row provides Cartesian
#' coordinates for the location of a mask point. The function
#' \link{create.mask} will return a suitable object.
#' @param detfn A character string specifying the detection function
#' to be used. Options are "hn" (halfnormal), "hr" (hazard rate), "th"
#' (threshold), "lth" (log-link threshold), or "ss" (signal
#' strength). If the latter is used, signal strength information must
#' be provided in \code{capt}.
#' @param sv A named list. Component names are parameter names, and
#' each component is a start value for the associated parameter. See
#' 'Details' for further information on the parameters to be fitted.
#' @param bounds A named list. Component names are parameter names,
#' and each components is a vector of length two, specifying the
#' bounds for the associated parameter.
#' @param fix A named list. Component names are parameter names to be
#' fixed, and each component is the fixed value for the associated
#' parameter.
#' @param sf A named list. Component names are parameter names, and
#' each component is a scalefactor for the associated parameter. The
#' default behaviour is to automatically select scalefactors based on
#' parameter start values. See the section on convergence below.
#' @param ss.link A character string, either \code{"identity"} or
#' \code{"log"}, which specifies the link function for the signal
#' strength detection function. Only required when \code{detfn} is
#' \code{"ss"} (i.e., when there is signal strength information in
#' \code{capt}).
#' @param cutoff The signal strength threshold, above which sounds are
#' identified as detections. Only required when \code{detfn} is
#' \code{"ss"}.
#' @param call.freqs A vector of call frequencies collected
#' independently to an acoustic survey.
#' @param sound.speed The speed of sound in metres per second,
#' defaults to 330 (the speed of sound in air). Only used when
#' \code{"toa"} is a component name of \code{capt}.
#' @param local Logical, if \code{TRUE} integration over unobserved
#' animal activity centres is only carried out in a region local to
#' detectors that detected individuals. See 'Details'.
#' @param hess Logical, if \code{TRUE} the Hessian is estimated,
#' allowing for calculation of standard errors, the
#' variance-covariance matrix, and the correlation matrix, at the
#' expense of a little processing time. If \code{FALSE}, the Hessian
#' is not estimated. Note that if individuals are detectable more than
#' once (e.g., by calling more than once on an acoustic survey) then
#' parameter uncertainty is not properly represented by these
#' calculations.
#' @param trace Logical, if \code{TRUE} parameter values at each step
#' of the optimisation algorithm are printed to the R console.
#' @param clean Logical, if \code{TRUE} ADMB output files are
#' removed. Otherwise, ADMB output file will remain in a directory,
#' the location of which is reported after the model is fitted.
#' @param cbs The CMPDIF_BUFFER_SIZE, set using the \code{-cbs} option
#' of the executable created by ADMB. This can be increased to speed
#' up optimisation if \code{cmpdiff.tmp} gets too large (please
#' ignore, unless you are familiar with ADMB and know what you are
#' doing).
#' @param gbs The GRADSTACK_BUFFER_SIZE, set using the \code{-gbs}
#' option of the executable created by ADMB. This can be increased to
#' speed up optimisation if \code{gradfil1.tmp} gets too large (please
#' ignore, unless you are familiar with ADMB and know what you are
#' doing).
#' @param exe.type Character string, either \code{"old"} or
#' \code{"new"}, depending on which executable is to be used (for
#' development purposes only; please ignore).
#'
#' @seealso \link{boot.admbsecr} to calculate standard errors and
#' estimate bias using a parametric bootstrap.
#' @seealso \link{coef.admbsecr}, \link{stdEr.admbsecr}, and
#' \link{vcov.admbsecr} to extract estimated parameters, standard
#' errors, and the variance-covariance matrix, respectively.
#' @seealso \link{confint.admbsecr} to calculate confidence intervals.
#' @seealso \link{summary.admbsecr} to get a summary of estimates and
#' standard errors.
#' @seealso \link{show.detfn} to plot the estimated detection
#' function.
#' @seealso \link{locations} to plot estimated locations of particular
#' individuals or calls.
#'
#' @examples
#' \dontrun{
#' simple.capt <- example$capt["bincapt"]
#' simple.hn.fit <- admbsecr(capt = simple.capt, traps = example$traps,
#' mask = example$mask, fix = list(g0 = 1))
#' simple.hr.fit <- admbsecr(capt = simple.capt, traps = example$traps,
#' mask = example$mask, detfn = "hr")
#' bearing.capt <- example.capt[c("bincapt", "bearing")]
#' bearing.hn.fit <- admbsecr(capt = bearing.capt, traps = example$traps,
#' mask = example$mask, fix = list(g0 = 1))
#' }
#'
#' @export
#'
admbsecr <- function(capt, traps, mask, detfn = "hn", sv = NULL, bounds = NULL,
fix = NULL, sf = NULL, ss.link = "identity",
cutoff = NULL, call.freqs = NULL, sound.speed = 330,
local = FALSE, hess = !any(call.freqs > 1), trace = FALSE,
clean = TRUE, cbs = NULL, gbs = NULL, exe.type = "old"){
arg.names <- names(as.list(environment()))
capt.bin <- capt$bincapt
## Checking for bincapt.
if (is.null(capt.bin)){
stop("The binary capture history must be provided as a component of 'capt'.")
}
## Checking for correct number of trap locations.
if (ncol(capt.bin) != nrow(traps)){
stop("There must be a trap location for each column in the components of 'capt'.")
}
## Checking that each component of 'capt' is a matrix.
if (any(!laply(capt, is.matrix))){
stop("At least one component of 'capt' is not a matrix.")
}
## Checking for agreement in matrix dimensions.
if (length(capt) > 1){
all.dims <- laply(capt, dim)
if (any(aaply(all.dims, 2, function(x) diff(range(x))) != 0)){
stop("Components of 'capt' object have different dimensions.")
}
}
## Various checks for other arguments.
if (!is.list(sv) & !is.null(sv)){
stop("The 'sv' argument must be 'NULL' or a list.")
}
if (!is.list(bounds) & !is.null(bounds)){
stop("The 'bounds' argument must be 'NULL' or a list.")
}
if (is.list(bounds)){
if (any(laply(bounds, length) != 2)){
stop("Each component of 'bounds' must be a vector of length 2.")
}
}
if (!is.list(fix) & !is.null(fix)){
stop("The 'fix' argument must be 'NULL' or a list.")
}
n <- nrow(capt.bin)
n.traps <- nrow(traps)
n.mask <- nrow(mask)
A <- attr(mask, "area")
buffer <- attr(mask, "buffer")
## Removing attributes from mask.
mask <- as.matrix(mask)
attr(mask, "area") <- A
attr(mask, "buffer") <- buffer
## TODO: Sort out how to determine supplementary parameter names.
supp.types <- c("bearing", "dist", "ss", "toa", "mrds")
fit.types <- supp.types %in% names(capt)
names(fit.types) <- supp.types
## Logical indicators for additional information types.
fit.bearings <- fit.types["bearing"]
fit.dists <- fit.types["dist"]
fit.ss <- fit.types["ss"]
fit.toas <- fit.types["toa"]
fit.mrds <- fit.types["mrds"]
## Generating ordered binary capture history.
capt.bin.order <- do.call(order, as.data.frame(capt.bin))
capt.bin.unique <- capt.bin[capt.bin.order, ]
capt.bin.freqs <- as.vector(table(apply(capt.bin.unique, 1, paste, collapse = "")))
names(capt.bin.freqs) <- NULL
capt.bin.unique <- capt.bin.unique[!duplicated(as.data.frame(capt.bin.unique)), ]
n.unique <- nrow(capt.bin.unique)
unique.changes <- cumsum(c(0, capt.bin.freqs[-n.unique])) + 1
## Reordering all capture history components.
capt.ord <- capt
for (i in 1:length(capt)){
capt.ord[[i]] <- capt[[i]][capt.bin.order, ]
}
## Capture histories for additional information types (if they exist)
capt.bearing <- if (fit.bearings) capt.ord$bearing else 0
capt.dist <- if (fit.dists) capt.ord$dist else 0
capt.ss <- if (fit.ss) capt.ord$ss else 0
capt.toa <- if (fit.toas) capt.ord$toa else 0
mrds.dist <- if (fit.mrds) capt.ord$mrds else 0
suppar.names <- c("kappa", "alpha", "sigma.toa")[fit.types[c("bearing", "dist", "toa")]]
if (fit.ss){
## Warning for failure to provide 'cutoff'.
if (missing(cutoff)){
stop("Argument 'cutoff' is missing.")
}
if (!missing(detfn) & detfn != "ss"){
warning("Argument 'detfn' is being ignored as signal strength information is provided in 'capt'. A signal strength detection function has been fitted instead.")
}
if (ss.link == "identity"){
detfn <- "ss"
linkfn.id <- 1
} else if (ss.link == "log"){
detfn <- "log.ss"
linkfn.id <- 2
} else {
stop("ss.link must be either \"identity\" or \"log\"")
}
} else {
## Not sure what a linkfn.id of 3 means? Probably throws an error in ADMB.
linkfn.id <- 3
}
detfns <- c("hn", "hr", "th", "lth", "ss", "log.ss")
## Sets detection function ID number for use in ADMB:
## 1 = Half normal
## 2 = Hazard rate
## 3 = Threshold
## 4 = Log-link threshold
## 5 = Identity-link signal strength
## 6 = Log-link signal strength.
detfn.id <- which(detfn == detfns)
detpar.names <- switch(detfn,
hn = c("g0", "sigma"),
hr = c("g0", "sigma", "z"),
th = c("shape", "scale"),
lth = c("shape.1", "shape.2", "scale"),
ss = c("b0.ss", "b1.ss", "sigma.ss"),
log.ss = c("b0.ss", "b1.ss", "sigma.ss"))
par.names <- c("D", detpar.names, suppar.names)
n.detpars <- length(detpar.names)
n.suppars <- length(suppar.names)
any.suppars <- n.suppars > 0
n.pars <- length(par.names)
## Checking par.names against names of sv, fix, bounds, and sf.
for (i in c("sv", "fix", "bounds", "sf")){
obj <- get(i)
if (!is.null(obj)){
obj.fitted <- names(obj) %in% par.names
if(!all(obj.fitted)){
msg <- paste("Some parameters listed in '", i, "' are not being used. These are being removed.",
sep = "")
warning(msg)
assign(i, obj[obj.fitted])
}
}
}
## Sets link function ID number for use in ADMB:
## 1 = identity
## 2 = log
## 3 = logit
links <- list(D = 2,
g0 = 3,
sigma = 2,
shape = 1,
shape.1 = 2,
shape.2 = 1,
scale = 2,
b0.ss = 2,
b1.ss = 2,
sigma.ss = 2,
z = 2,
sigma.toa = 2,
kappa = 2,
alpha = 2)[par.names]
link.list <- list(identity, log.link, logit.link)
unlink.list <- list(identity, exp, inv.logit)
par.links <- llply(links, function(x, link.list) link.list[[x]], link.list)
par.unlinks <- llply(links, function(x, unlink.list) unlink.list[[x]], unlink.list)
## Sorting out start values. Start values are set to those provided,
## or else are determined automatically from functions in
## autofuns.r.
sv.link <- vector("list", length = n.pars)
names(sv.link) <- par.names
sv.link[names(sv)] <- sv
sv.link[names(fix)] <- fix
auto.names <- par.names[sapply(sv.link, is.null)]
sv.funs <- paste("auto", auto.names, sep = "")
## Done in reverse so that D is calculated last (requires detfn parameters).
## D not moved to front as it should appear as the first parameter in any output.
for (i in rev(seq(1, length(auto.names), length.out = length(auto.names)))){
sv.link[auto.names[i]] <- eval(call(sv.funs[i],
list(capt = capt, detfn = detfn,
detpar.names = detpar.names,
mask = mask, traps = traps,
sv = sv.link, cutoff = cutoff,
ss.link = ss.link, A = A)))
}
## Converting start values to link scale.
sv <- sv.link
for (i in names(sv.link)){
sv.link[[i]] <- link.list[[links[[i]]]](sv.link[[i]])
}
## Sorting out phases.
## TODO: Add phases parameter so that these can be controlled by user.
phases <- vector("list", length = n.pars)
names(phases) <- par.names
for (i in par.names){
if (any(i == names(fix))){
## Phase of -1 in ADMB fixes parameter at starting value.
phases[[i]] <- -1
} else {
phases[[i]] <- 0
}
}
D.phase <- phases[["D"]]
detpars.phase <- c(phases[detpar.names], recursive = TRUE)
if (any.suppars){
suppars.phase <- c(phases[suppar.names], recursive = TRUE)
} else {
suppars.phase <- -1
}
## Sorting out bounds.
## Below bounds are the defaults.
default.bounds <- list(D = c(n/(A*n.mask), 1e8),
g0 = c(0, 1),
sigma = c(0, 1e8),
shape = c(-100, 100),
shape.1 = c(0, 1e8),
shape.2 = c(-100, 100),
scale = c(0, 1e8),
b0.ss = c(0, 1e8),
b1.ss = c(0, 1e8),
sigma.ss = c(0, 1e8),
z = c(0, 1e8),
sigma.toa = c(0, 1e8),
kappa = c(0, 700),
alpha = c(0, 1e8))[par.names]
bound.changes <- bounds
bounds <- default.bounds
for (i in names(default.bounds)){
if (i %in% names(bound.changes)){
bounds[[i]] <- bound.changes[[i]]
}
}
## Converting bounds to link scale.
bounds.link <- bounds
for (i in names(bounds)){
bounds.link[[i]] <- link.list[[links[[i]]]](bounds[[i]])
}
D.bounds <- bounds.link[["D"]]
D.lb <- D.bounds[1]
D.ub <- D.bounds[2]
detpar.bounds <- bounds.link[detpar.names]
detpars.lb <- sapply(detpar.bounds, function(x) x[1])
detpars.ub <- sapply(detpar.bounds, function(x) x[2])
if (any.suppars){
suppar.bounds <- bounds.link[suppar.names]
suppars.lb <- sapply(suppar.bounds, function(x) x[1])
suppars.ub <- sapply(suppar.bounds, function(x) x[2])
} else {
suppars.lb <- 0
suppars.ub <- 0
}
## Sorting out scalefactors.
if (is.null(sf)){
sv.vec <- c(sv.link, recursive = TRUE)
## Currently, by default, the scalefactors are the inverse
## fraction of each starting value to the largest starting
## value. Not sure how sensible this is.
##bound.ranges <- laply(bounds.link, function(x) diff(range(x)))
##sf <- max(bound.ranges)/bound.ranges
sf <- abs(max(sv.vec)/sv.vec)
names(sf) <- par.names
} else if (is.list(sf)){
sf <- numeric(n.pars)
names(sf) <- par.names
for (i in par.names){
sf[i] <- ifelse(i %in% names(sf), sf[[i]], 1)
}
} else if (is.vector(sf) & length(sf) == 1){
sf <- rep(sf, length(par.names))
names(sf) <- par.names
}
## Replacing infinite scalefactors.
sf[!is.finite(sf)] <- 1
D.sf <- sf[["D"]]
detpars.sf <- c(sf[detpar.names], recursive = TRUE)
if (any.suppars){
suppars.sf <- c(sf[suppar.names], recursive = TRUE)
} else {
suppars.sf <- 1
}
## Creating link objects to pass to ADMB.
detpars.link <- c(links[detpar.names], recursive = TRUE)
if (any.suppars){
suppars.link <- c(links[suppar.names], recursive = TRUE)
} else {
suppars.link <- 1
}
## Setting small number so that numerical under/overflow in ADMB
## does not affect estimation.
dbl.min <- 1e-150
## Calculating distances and angles.
dists <- distances(traps, mask)
if (fit.bearings){
bearings <- bearings(traps, mask)
} else {
bearings <- 0
}
if (fit.toas){
toa.ssq <- make_toa_ssq(capt.ord$toa, dists, sound.speed)
} else {
toa.ssq <- 0
}
if (is.null(cutoff)){
cutoff <- 0
}
## Kludge to fix number of parameters for no supplementary
## information.
if (!any.suppars){
n.suppars <- max(c(n.suppars, 1))
sv.link$dummy <- 0
}
## Sorting out which mask points are local to each detection.
if (local){
all.which.local <- find_local(capt.bin.unique, dists, buffer)
all.n.local <- laply(all.which.local, length)
all.which.local <- c(all.which.local, recursive = TRUE)
} else {
all.n.local <- rep(1, n.unique)
all.which.local <- rep(0, n.unique)
}
## Stuff for the .dat file.
data.list <- list(
n_unique = n.unique, local = as.numeric(local), all_n_local = all.n.local,
all_which_local = all.which.local, D_lb = D.lb, D_ub = D.ub, D_phase =
D.phase, D_sf = D.sf, n_detpars = n.detpars, detpars_lb = detpars.lb,
detpars_ub = detpars.ub, detpars_phase = detpars.phase, detpars_sf =
detpars.sf, detpars_linkfns = detpars.link, n_suppars = n.suppars,
suppars_lb = suppars.lb, suppars_ub = suppars.ub, suppars_phase =
suppars.phase, suppars_sf = suppars.sf, suppars_linkfns =
suppars.link, detfn_id = detfn.id, buffer = buffer, trace =
as.numeric(trace), DBL_MIN = dbl.min, n = n, n_traps = n.traps, n_mask
= n.mask, A = A, capt_bin_unique = capt.bin.unique, capt_bin_freqs =
capt.bin.freqs, fit_angs = as.numeric(fit.bearings), capt_ang =
capt.bearing, fit_dists = as.numeric(fit.dists), capt_dist =
capt.dist, fit_ss = as.numeric(fit.ss), cutoff = cutoff, linkfn_id =
linkfn.id, capt_ss = capt.ss, fit_toas = as.numeric(fit.toas),
capt_toa = capt.toa, fit_mrds = as.numeric(fit.mrds), mrds_dist =
mrds.dist, dists = dists, angs = bearings, toa_ssq = toa.ssq)
## Determining whether or not standard errors should be calculated.
if (!is.null(call.freqs)){
fit.freqs <- any(call.freqs != 1)
} else {
fit.freqs <- FALSE
}
## Idea of running executable as below taken from glmmADMB.
## Working out correct command to run from command line.
if (exe.type == "new"){
exe.name <- "secr_new"
} else if (exe.type == "old"){
exe.name <- "secr"
} else if (exe.type == "test"){
exe.name <- "secr_test"
} else {
stop("Argument 'exe.type' must be \"old\" or \"new\".")
}
prefix.name <- exe.name
if (.Platform$OS == "windows"){
os.type <- "windows"
exe.name <- paste(prefix.name, ".exe", sep = "")
} else if (.Platform$OS == "unix"){
if (Sys.info()["sysname"] == "Linux"){
os.type <- "linux"
} else if (Sys.info()["sysname"] == "Darwin"){
os.type <- "mac"
} else {
stop("Unknown OS type.")
}
} else {
stop("Unknown OS type.")
}
## Finding executable folder (possible permission problems?).
exe.dir <- paste(system.file(package = "admbsecr"), "ADMB", "bin", os.type, sep = "/")
exe.loc <- paste(exe.dir, exe.name, sep = "/")
## Creating command to run using system().
curr.dir <- getwd()
## Creating temporary directory.
temp.dir <- tempfile("admbsecr", curr.dir)
dir.create(temp.dir)
setwd(temp.dir)
## Creating .pin and .dat files.
write_pin("secr", sv.link)
write_dat("secr", data.list)
## Creating link to executable.
if (os.type == "windows"){
file.copy(exe.loc, exe.name)
} else {
file.symlink(exe.loc, exe.name)
}
## Sorting out -cbs and -gbs.
if (!is.null(cbs)){
cbs.cmd <- paste(" -cbs", format(cbs, scientific = FALSE))
} else {
cbs.cmd <- NULL
}
if (!is.null(gbs)){
gbs.cmd <- paste(" -gbs", format(gbs, scientific = FALSE))
} else {
gbs.cmd <- NULL
}
## Running ADMB executable.
cmd <- paste("./"[os.type != "windows"], exe.name,
" -ind secr.dat -ainp secr.pin",
" -nohess"[!hess], cbs.cmd, gbs.cmd, sep = "")
if (os.type == "windows"){
system(cmd, ignore.stdout = !trace, show.output.on.console = trace)
} else {
system(cmd, ignore.stdout = !trace)
}
## Reading in model results.
options(warn = -1)
if (exe.type == "test"){
prefix.name <- strsplit(list.files(), "\\.")[[which(substr(list.files(),
nchar(list.files()) - 3,
nchar(list.files())) == ".par")]][1]
}
out <- read.admbsecr(prefix.name)
options(warn = 0)
setwd(curr.dir)
## Cleaning up files.
if (clean){
unlink(temp.dir, recursive = TRUE)
} else {
cat("ADMB files found in:", "\n", temp.dir, "\n")
}
## Warning for non-convergence.
if (out$maxgrad < -0.1){
warning("Failed convergence -- maximum gradient component is large.")
}
## Moving back to original working directory.
setwd(curr.dir)
## Removing fixed coefficients from list.
if (!hess){
out$coeflist[c(D.phase, detpars.phase, suppars.phase) == -1] <- NULL
}
## Creating coefficients vector.
est.pars <- c("D", detpar.names, suppar.names)[c(D.phase, detpars.phase, suppars.phase) > -1]
n.est.pars <- length(est.pars)
out$coefficients <- numeric(2*n.est.pars + 1)
names(out$coefficients) <- c(paste(est.pars, "_link", sep = ""), est.pars, "esa")
for (i in 1:n.est.pars){
out$coefficients[i] <- out$coeflist[[i]]
}
for (i in 1:n.est.pars){
out$coefficients[n.est.pars + i] <-
unlink.list[[links[[est.pars[i]]]]](out$coeflist[[i]])
}
## Adding extra components to list.
if (detfn == "log.ss") detfn <- "ss"
## Putting in updated argument names.
args <- vector(mode = "list", length = length(arg.names))
names(args) <- arg.names
for (i in arg.names){
if (!is.null(get(i))){
args[[i]] <- get(i)
}
}
out$args <- args
out$fit.types <- fit.types
out$infotypes <- names(fit.types)[fit.types]
out$detpars <- detpar.names
out$suppars <- suppar.names
out$phases <- phases
out$par.links <- par.links
out$par.unlinks <- par.unlinks
## Putting in esa estimate.
out$coefficients[2*n.est.pars + 1] <- p.dot(out, esa = TRUE)
## Putting in call frequency information and correct parameter names.
if (fit.freqs){
mu.freqs <- mean(call.freqs)
Da <- get.par(out, "D")/mu.freqs
names.vec <- c(names(out[["coefficients"]]), "Da", "mu.freqs")
coefs.updated <- c(out[["coefficients"]], Da, mu.freqs)
names(coefs.updated) <- names.vec
out[["coefficients"]] <- coefs.updated
## Removing ses, cor, vcov matrices.
cor.updated <- matrix(NA, nrow = length(names.vec),
ncol = length(names.vec))
dimnames(cor.updated) <- list(names.vec, names.vec)
vcov.updated <- matrix(NA, nrow = length(names.vec),
ncol = length(names.vec))
dimnames(vcov.updated) <- list(names.vec, names.vec)
if (hess){
ses.updated <- c(out[["se"]], rep(NA, 2))
max.ind <- length(names.vec) - 2
cor.updated[1:max.ind, 1:max.ind] <- out[["cor"]]
vcov.updated[1:max.ind, 1:max.ind] <- out[["vcov"]]
} else {
ses.updated <- rep(NA, length(names.vec))
}
names(ses.updated) <- names.vec
out[["se"]] <- ses.updated
out[["cor"]] <- cor.updated
out[["vcov"]] <- vcov.updated
if (trace){
if (!hess){
cat("NOTE: Standard errors not calculated; use boot.admbsecr().", "\n")
} else {
cat("NOTE: Standard errors are probably not correct; use boot.admbsecr().", "\n")
}
}
} else {
if (hess){
## Putting correct parameter names into se, cor, vcov.
replace <- substr(names(out$se), 1, 8) == "par_ests"
names(out$se)[replace] <- rownames(out$vcov)[replace] <-
colnames(out$vcov)[replace] <- rownames(out$cor)[replace] <-
colnames(out$cor)[replace] <- est.pars
replace <- 1:length(est.pars)
names(out$se)[replace] <- rownames(out$vcov)[replace] <-
colnames(out$vcov)[replace] <- rownames(out$cor)[replace] <-
colnames(out$cor)[replace] <- paste(est.pars, "_link", sep = "")
} else {
## Filling se, cor, vcov with NAs.
names.vec <- names(out[["coefficients"]])
ses.updated <- rep(NA, length(names.vec))
names(ses.updated) <- names.vec
cor.updated <- matrix(NA, nrow = length(names.vec),
ncol = length(names.vec))
dimnames(cor.updated) <- list(names.vec, names.vec)
vcov.updated <- matrix(NA, nrow = length(names.vec),
ncol = length(names.vec))
dimnames(vcov.updated) <- list(names.vec, names.vec)
out[["se"]] <- ses.updated
out[["cor"]] <- cor.updated
out[["vcov"]] <- vcov.updated
}
}
out$fit.freqs <- fit.freqs
if (out$maxgrad < -0.01){
warning("Maximum gradient component is large.")
}
class(out) <- c("admbsecr", "admb")
out
}
## Roxygen code for NAMESPACE and datasets.
## Package imports for roxygenise to pass to NAMESPACE.
#' @import plyr Rcpp R2admb
#' @importFrom CircStats dvm rvm
#' @importFrom lattice wireframe
#' @importFrom matrixStats colProds
#' @importFrom secr make.capthist make.mask read.mask read.traps sim.popn
#' @useDynLib admbsecr
NULL
## Data documentation.
#' \emph{Arthroleptella lightfooti} survey data
#'
#' Data from an acoustic survey of the Western Cape moss frog
#' \emph{Arthroleptella lightfooti}. These data are from a 25 s subset
#' of the original recording, taken on 16 May 2012 at Silvermine,
#' Table Mountain National Park, Cape Town, South Africa. Acoustic
#' signal strengths and times of arrival were measured, and this
#' information is contained in the capture history object.
#'
#' This object is a list which contains components:
#' \itemize{
#' \item \code{capt}: A capture history object.
#' \item \code{traps}: A traps object.
#' \item \code{mask}: A suitable mask object.
#' \item \code{cutoff}: The microphone cutoff value.
#' \item \code{freqs}: A vector of call frequencies measured
#' independently to the acoustic survey.
#' }
#'
#' @name lightfooti
#' @format A list.
#' @usage lightfooti
#' @docType data
#' @keywords datasets
NULL
#' Example data
#'
#' This object contains simulated data with all types of supplementary
#' information, corresponding trap locations, and a suitable mask
#' object. Also included are some example model fits, which were
#' generated from these data using the \link{admbsecr} function.
#'
#' This object is a list which contains components:
#' \itemize{
#' \item \code{capt}: A capture history object.
#' \item \code{traps}: A traps object.
#' \item \code{mask}: A suitable mask object.
#' \item \code{cutoff}: The cutoff value used to simluate these data.
#' \item \code{fits}: Some example model fits.
#' }
#'
#' @name example
#' @format A list.
#' @usage example
#' @docType data
#' @keywords datasets
NULL
Try the admbsecr 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.