\vspace{16pt}
The small package secrBVN is used to evaluate the performance of SECR estimators when home ranges are bivariate normal (BVN) or uniform (flat-topped) ellipses. We assume detection hazard is directly proportional to home range utilisation (activity). Code to use secrBVN for the simulations of Efford (in prep.) is provided in the Appendix.
The key user-visible functions are simpopn.bvn
, plotpopn.bvn
, simcapt.bvn
, runEllipseSim
, simsum
, simplot
and anisotropic.fit
.
\vspace{12pt}
simpopn.bvn
is a wrapper for the secr function sim.popn
that adds attributes specifying a bivariate normal home range shape, size and orientation for each individual. By default, shape and size are the same for all individuals, but a mechanism is provided to vary them individually (see Heterogeneous elliptical home ranges).
First, load the package.
library(secrBVN) simfolder <- "d:/density communication/noncircularity/paper/simulations/" runall <- FALSE # skip lengthy simulations
tempgrid <- make.grid(nx = 10, ny = 10) par(mfrow = c(2,4), mar = c(2,2,2.6,2), xpd = TRUE) for (i in 1:4) { s2xy <- 25^2 * c(1/i, i) random.pop <- simpopn.bvn(s2xy = s2xy, core = tempgrid, buffer = 100, D = 1) plotpopn.bvn(random.pop, col = 'lightblue') mtext(side=3, line=1.5, i) } for (i in 1:4) { s2xy <- 25^2 * c(1/i, i) aligned.pop <- simpopn.bvn(s2xy = s2xy, core = tempgrid, buffer = 100, D = 1, theta = -1) plotpopn.bvn(aligned.pop, col = 'lightblue') }
Fig. 1. Elliptical home ranges with varying ratio of major to minor axes as shown. Upper row oriented randomly and independently, lower row with shared random orientation ('randomly aligned').
Normally in secr we use sim.capthist
to generate capture histories, but that is limited to circular detection functions. The secrBVN function simcapt.bvn
is a partial replacement for sim.capthist
that models detection with a bivariate normal. Specifically, the cumulative hazard of detection is a constant times the bivariate normal probability density at the detector. The constant is $\lambda_0 2 \pi \sigma_X \sigma_Y$. The user provides a 'popn' object that includes the BVN parameter values ($\sigma_x^2, \sigma_y^2, \theta$) for each animal (row), as generated above by simpopn.bvn
. The constant scales the BVN density so that the maximum hazard is $\lambda_0$.
The preceding paragraph describes the default behaviour of simcapt.bvn
. If type = uniform
is selected then a uniform (flat-topped) elliptical home range is simulated. The uniform probability of detection within the ellipse is controlled by g0
, not lambda0
.
The function runEllipseSim
is a wrapper for the preceding steps (simpopn.bvn
, simcapt.bvn
) that also fits a standard (circular) SECR model with secr.fit
[^footnote1]. The default population has a fixed number of individuals within the rectangular buffered area around the detectors (Ndist = 'fixed'
).
[^footnote1]: An elliptical SECR model also may be fitted - see Anisotropic home ranges.
For this example we use a $6 \times 6$ grid of binary proximity detectors with 50 metre spacing. The code in secrBVN does not allow for competition among detectors (secr detector type 'multi') or other other secr detector types. A density of 4/ha gives exactly 169 animals in the buffered area. Conditional likelihood is used for speed; the default extractfn = derived
is compatible with both CL = TRUE
and CL = FALSE
. The 200-m buffer allows for the longest ranges ($\sigma_y = 50$ m).
tr <- make.grid(6, 6, spacing = 50, detector = 'proximity') simrandomBVN36 <- vector('list', 4) for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 simrandomBVN36[[i]] <- runEllipseSim (nrepl = 500, sigmaX, sigmaY, buffer = 200, ncores = 20, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = TRUE, detectfn = 'HHN', details = list(distribution = 'binomial')) }
The secrBVN function simplot
is used to summarize the results
load(file = paste0(simfolder, 'simrandomBVN36.RData'))
par(xpd = FALSE, mar = c(4,4,4,4)) simplot(list(BVN = simrandomBVN36))
Fig. 2. Relative bias of density estimated by fitting circular SECR detection model to elliptical data, with 95\% confidence limit for simulated values.
Next print a summary of the simulation results. There is no apparent effect of range elongation itself on the bias of the estimates.
simsum(list(BVN = simrandomBVN36))
The Appendix has many other examples.
simpopn.bvn
by default generates a population with equal-sized home ranges. The size and elongation of each range may be varied by providing a user-defined function as the argument s2xy
:
tr <- make.grid(6,6, spacing = 50, detector = 'proximity') rs2xy <- function(N, scale = 25) { aspectratio <- 1 + runif(N) * 3 cbind(scale^2 / aspectratio, scale^2 * aspectratio) } pop <- simpopn.bvn(s2xy = rs2xy, core = tr, buffer = 100, D = 1) par(mfrow = c(1,1), mar = c(4,4,4,4), xpd = TRUE) plotpopn.bvn(pop, col = 'grey')
Heterogeneous populations may also be simulated in runEllipseSim
by passing a function as the argument 'sigmaX'.
sims <- runEllipseSim (10, sigmaX = rs2xy, buffer = 200, ncores = 2, traps = tr, g0 = 0.2, D = 4, type = 'uniform')
In principle, we can deal with elongated ranges by replacing Euclidean distances with distances in a space transformed to render home ranges circular (Murphy et al. 2016). The transformation compresses home ranges along their major axis. The transformation uses two parameters: $\psi_A$ for the shared orientation measured in radians and $\psi_R$ for the compression ratio ($\psi_R \ge 1$; $\psi_R = 1$ corresponds to no compression). The parameter $\psi_A$ (psiA) corresponds to the argument theta
used by simpopn.bvn
. Thanks to Ben Augustine for pointing out the geoR function coords.aniso
that lets us do this (Ribeiro and Diggle 2018)[^footnote].
[^footnote]:The package intamap (Pebesma et al. 2010) also offers anisotropic transformation in function rotateAnisotropicData
.
The method does not work if the detector array is strictly linear (2-D data are required to estimate the orientation and compression parameters). It is tested in the Appendix on simulated data from a hollow square array.
The code uses a user-defined distance function that computes a Euclidean distance in the transformed space. In this version both transformation parameters are estimated (cf Murphy et al. 2016). Transformation parameters are handled in the undocumented 'miscparm' feature of secr.fit
. This allows supernumerary unmodelled parameters to be passed to the distance function. The parameters are included in the vector of coefficients (beta parameters) over which the likelihood is maximised. 'Unmodelled' here means that the parameter takes a single value for all animals, times and places. The link function is 'identity' for $\psi_A$ and 'log' for $\psi_R-1$.
anisodistfn <- function (xy1, xy2, mask) { if (missing(xy1)) return(character(0)) xy1 <- as.matrix(xy1) xy2 <- as.matrix(xy2) miscparm <- attr(mask, 'miscparm') psiA <- miscparm[1] # anisotropy angle; identity link psiR <- 1 + exp(miscparm[2]) # anisotropy ratio; log link aniso.xy1 <- geoR::coords.aniso(xy1, aniso.pars = c(psiA, psiR)) aniso.xy2 <- geoR::coords.aniso(xy2, aniso.pars = c(psiA, psiR)) secr::edist(aniso.xy1, aniso.xy2) # nrow(xy1) x nrow(xy2) matrix }
Here is a simple application with oblique elliptical home ranges (theta = pi/4
).
tr <- make.grid(10, 10, spacing = 25, hollow = TRUE, detector = 'proximity') pop <- simpopn.bvn(s2xy = (25*c(0.5,2))^2, theta = pi/4, core = tr, buffer = 200, D = 4, Ndist = 'fixed') CH <- simcapt.bvn(tr, pop, type = 'BVN', lambda = 0.4, noccasions = 5)
par(mfrow = c(1,1)) plot(CH, tracks = TRUE) plotpopn.bvn(pop, border='grey', add = TRUE) plot(tr, add = TRUE)
First we fit a naive circular model.
fit0 <- secr.fit(CH, buffer = 200, detectfn = 'HHN', trace = FALSE, details = list(distribution = 'binomial')) predict(fit0)
Next, the model with transformation to isotropy.
details <- list(distribution = 'binomial', userdist = anisodistfn, miscparm = c(psiA = 0.5, psiR = 1.5)) # initial values fit1 <- secr.fit(CH, buffer = 200, detectfn = 'HHN', trace = FALSE, details = details) predict(fit1) coef(fit1)
The estimated bearing is r round(coef(fit1)['psiA','beta']*360/2/pi,2)
degrees. The estimated aspect ratio is r round(1 + exp(coef(fit1)['psiR', 'beta']),2)
, (95% CI r round(1 + exp(coef(fit1)['psiR', c('lcl','ucl')]), 2)
).
The secrBVN function anisotropic.fit
streamlines model fitting and does not require anisodistfn
to be defined externally. Use it as you would use secr.fit
. For example,
fit2 <- anisotropic.fit(CH, buffer = 200, detectfn = 'HHN', trace = FALSE, details = list(distribution = "binomial")) predictAniso(fit2) # estimates of transformation parameters psiA (in degrees) and psiR
If the orientation is known then there are two solutions:
Here we demonstrate (2), fixing psiA at pi/4 (45 degrees). The secr.fit
details argument fixedbeta is a vector of values, one per 'beta' parameter, using NA for each beta to be estimated. The first three beta parameters are for D, lambda0 and sigma; the next two are for psiA and psiR. If you model D, lambda0 or sigma there will be extra beta parameters, pushing psiA and psiR back.
fit3 <- anisotropic.fit(CH, buffer = 200, detectfn = 'HHN', trace = FALSE, details = list(distribution = "binomial", fixedbeta = c(NA,NA,NA,pi/4,NA))) predictAniso(fit3) # estimate of psiR only predict(fit3)
Detectors in a straight line provide very little information with which to estimate $\psi_A$ and $\psi_R$. Here is an example.
tr <- make.grid(100, 1, spacing = 25, detector = 'proximity') pop <- simpopn.bvn(s2xy = (25*c(0.5,2))^2, theta = pi/4, core = tr, buffer = 200, D = 4, Ndist = 'fixed') CH <- simcapt.bvn(tr, pop, type = 'BVN', lambda = 0.4, noccasions = 5) fit4 <- anisotropic.fit(CH, buffer = 200, detectfn = 'HHN', trace = FALSE, details = list(distribution = "binomial")) predictAniso(fit4) # estimates of transformation parameters psiA (in degrees) and psiR predict(fit4)
(Perhaps the parameters are not strictly nonidentifiable).
This package has the limited goal of determining how range elongation affects estimates of density in simple SECR models, and these specific limitations:
runEllipseSim
) or as a vector of the two values, squared ('s2xy' in simpopn.bvn
). This is confusing but it's better at this point not to change.Efford, M. G. (2004) Density estimation in live-trapping studies. Oikos 106, 598--610.
Efford, M. G. In prep. Non-circular home ranges and the estimation of population density.
Huggins, R. M. (1989) On the statistical analysis of capture experiments. Biometrika 76, 133--140.
Ivan, J. S., White, G. C. and Shenk, T. M. (2013) Using simulation to compare methods for estimating density from capture--recapture data. Ecology 94, 817--826.
Murphy, S. M., Cox, J. J., Augustine, B. C., Hast, J. T., Guthrie, J. M., Wright, J., McDermott, J., Maehr, S. C. and Plaxico, J. H. (2016) Characterizing recolonization by a reintroduced bear population using genetic spatial capture--recapture. Journal of Wildlife Management 80, 1390--1407.
Pebesma, E., Cornford, D., Dubois, G., Heuvelink, G.B.M., Hristopoulos, D., Pilz, J., Stoehlker, U., Morin, G. and Skoien, J.O. (2010) INTAMAP: the design and implementation of an interoperable automated interpolation web service. Computers & Geosciences 37, 343--352.
Ribeiro, P. J. Jr and Diggle, P. J. (2018) geoR: Analysis of Geostatistical Data. R package version 1.7-5.2.1. https://CRAN.R-project.org/package=geoR.
\pagebreak
library(secrBVN) simfolder <- "d:/density communication/noncircularity/paper/simulations/" nrepl <- 500 ncores <- 20 # for machine with at least 20 cores runsim <- FALSE # change to TRUE to execute simulations details <- list(distribution = 'binomial') # all simulations assume fixed N
Functions for heterogeneous aspect ratio.
# Uniform on 1-4 rs2xy <- function(N, scale = 25) { i <- 1 + runif(N) * 3 cbind(scale^2 /i, scale^2 * i) } # 2 classes 1,4 rs2xy2 <- function(N, scale = 25, prob = c(0.5,0.5)) { i <- sample(c(1,4), size = N, replace = TRUE, prob = prob) cbind(scale^2 /i, scale^2 * i) }
Elongated ranges are oriented at random with respect to the grid (default theta = NULL
).
tr <- make.grid(10,10, spacing = 50, detector = 'proximity') simrandom100.3 <- vector('list', 6) names(simrandom100.3) <- c(1:4,'1-4','1,4') simrandomBVN100.3 <- simrandom100.3 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, theta = NULL, D = 4, CL = TRUE, detectfn = 'HHN', details = details, seed = 347) for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) # uniform args$type <- 'uniform'; args$g0 <- 0.2 simrandom100.3[[i]] <- do.call(runEllipseSim, args) # bvn args$type <- 'BVN'; args$lambda0 <- 0.4 simrandomBVN100.3[[i]] <- do.call(runEllipseSim, args) message ('Completed aspect ratio', i) } # heterogeneous aspect ratio 1,4, uniform args$type <- 'uniform'; args$sigmaX <- rs2xy simrandom100.3[[5]] <- do.call(runEllipseSim, args) # heterogeneous aspect ratio 1,4, BVN args$type <- 'BVN'; args$sigmaX <- rs2xy simrandomBVN100.3[[5]] <- do.call(runEllipseSim, args) # heterogeneous aspect ratio 1,4, uniform args$type <- 'uniform'; args$sigmaX <- rs2xy2 simrandom100.3[[6]] <- do.call(runEllipseSim, args) #heterogeneous aspect ratio 1,4, BVN args$type <- 'BVN'; args$sigmaX <- rs2xy2; args$seed <- 347 simrandomBVN100.3[[6]] <- do.call(runEllipseSim, args) save(simrandom100.3, file = paste0(simfolder, 'simrandom100.3.RData')) save(simrandomBVN100.3, file = paste0(simfolder, 'simrandomBVN100.3.RData'))
load(file = paste0(simfolder, 'simrandom100.3.RData')) load(file = paste0(simfolder, 'simrandomBVN100.3.RData')) par(mfrow = c(1,1), mar = c(4,4,4,4), xpd = FALSE) # select only first 4 scenarios for plotting simplot(list(Uniform = simrandom100.3[1:4], BVN = simrandomBVN100.3[1:4]), legend = TRUE) # tabulate all simsum(list(Uniform = simrandom100.3, BVN = simrandomBVN100.3), compact = NULL, dec = 4) simsum(list(Uniform = simrandom100.3, BVN = simrandomBVN100.3), dec = 4)
# compact table for paper simsum(list(Uniform = simrandom100.3, BVN = simrandomBVN100.3)) # spatial scale parameter output <- simsum(list(Uniform = simrandom100.3, BVN = simrandomBVN100.3), component = 'pred', parm = 'sigma', compact = c("av.parmhat", "sd.parmhat")) lapply(output, '/', 50) # in units of detector spacing
tr <- make.grid(36, 1, spacing = 50, detector = 'proximity') simgridalignlinw1.1 <- vector('list', 6) names(simgridalignlinw1.1) <- c(1:4,'1-4','1,4') simgridalignlinw4.1 <- simgridalignlinw3.1 <- simgridalignlinw2.1 <- simgridalignlinw1.1 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = TRUE, detectfn = 'HHN', details = details) for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) # bvn perpendicular to traps args$theta <- 0; simgridalignlinw1.1[[i]] <- do.call(runEllipseSim, args) # bvn oblique to traps args$theta <- pi/4; simgridalignlinw2.1[[i]] <- do.call(runEllipseSim, args) # bvn parallel to traps args$theta <- pi/2; simgridalignlinw3.1[[i]] <- do.call(runEllipseSim, args) # bvn random orientation args$theta <- NULL; simgridalignlinw4.1[[i]] <- do.call(runEllipseSim, args) } args <- c(baseargs, list(sigmaX = rs2xy)) args$theta <- 0; simgridalignlinw1.1[[5]] <- do.call(runEllipseSim, args) args$theta <- pi/4; simgridalignlinw2.1[[5]] <- do.call(runEllipseSim, args) args$theta <- pi/2; simgridalignlinw3.1[[5]] <- do.call(runEllipseSim, args) args$theta <- NULL; simgridalignlinw4.1[[5]] <- do.call(runEllipseSim, args) args <- c(baseargs, list(sigmaX = rs2xy2)) args$theta <- 0; simgridalignlinw1.1[[6]] <- do.call(runEllipseSim, args) args$theta <- pi/4; simgridalignlinw2.1[[6]] <- do.call(runEllipseSim, args) args$theta <- pi/2; simgridalignlinw3.1[[6]] <- do.call(runEllipseSim, args) args$theta <- NULL; simgridalignlinw4.1[[6]] <- do.call(runEllipseSim, args) save(simgridalignlinw1.1, file = paste0(simfolder, 'simgridalignlinw1.1.RData')) save(simgridalignlinw2.1, file = paste0(simfolder, 'simgridalignlinw2.1.RData')) save(simgridalignlinw3.1, file = paste0(simfolder, 'simgridalignlinw3.1.RData')) save(simgridalignlinw4.1, file = paste0(simfolder, 'simgridalignlinw4.1.RData'))
# compact table linear load(file = paste0(simfolder, 'simgridalignlinw1.1.RData')) load(file = paste0(simfolder, 'simgridalignlinw2.1.RData')) load(file = paste0(simfolder, 'simgridalignlinw3.1.RData')) load(file = paste0(simfolder, 'simgridalignlinw4.1.RData')) simsum(list(Perpendicular = simgridalignlinw1.1, Oblique = simgridalignlinw2.1, Parallel = simgridalignlinw3.1, Random = simgridalignlinw4.1))
tr <- make.grid(10, 10, spacing = 50, detector = 'proximity', hollow = TRUE) simgridalignsqw1.1 <- vector('list', 4) # 'w' for wide names(simgridalignsqw1.1) <- 1:4 simgridalignsqw3.1 <- simgridalignsqw2.1 <- simgridalignsqw1.1 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = TRUE, detectfn = 'HHN', details = details) for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) # bvn aligned to traps args$theta <- 0; simgridalignsqw1.1[[i]] <- do.call(runEllipseSim, args) # bvn oblique to traps args$theta <- pi/4; simgridalignsqw2.1[[i]] <- do.call(runEllipseSim, args) # bvn random orientation args$theta <- NULL; simgridalignsqw3.1[[i]] <- do.call(runEllipseSim, args) } args <- c(baseargs, list(sigmaX = rs2xy)) args$theta <- 0; simgridalignsqw1.1[[5]] <- do.call(runEllipseSim, args) args$theta <- pi/4; simgridalignsqw2.1[[5]] <- do.call(runEllipseSim, args) args$theta <- NULL; simgridalignsqw3.1[[5]] <- do.call(runEllipseSim, args) args <- c(baseargs, list(sigmaX = rs2xy2)) args$theta <- 0; simgridalignsqw1.1[[6]] <- do.call(runEllipseSim, args) args$theta <- pi/4; simgridalignsqw2.1[[6]] <- do.call(runEllipseSim, args) args$theta <- NULL; simgridalignsqw3.1[[6]] <- do.call(runEllipseSim, args) save(simgridalignsqw1.1, file = paste0(simfolder, 'simgridalignsqw1.1.RData')) save(simgridalignsqw2.1, file = paste0(simfolder, 'simgridalignsqw2.1.RData')) save(simgridalignsqw3.1, file = paste0(simfolder, 'simgridalignsqw3.1.RData'))
# compact table square load(file = paste0(simfolder, 'simgridalignsqw1.1.RData')) load(file = paste0(simfolder, 'simgridalignsqw2.1.RData')) load(file = paste0(simfolder, 'simgridalignsqw3.1.RData')) simsum(list(Aligned = simgridalignsqw1.1, Oblique = simgridalignsqw2.1, Random = simgridalignsqw3.1))
tr <- make.grid(10,10, spacing = 25, detector = 'proximity') simrandomBVNs100.1 <- vector('list', 4) names(simrandomBVNs100.1) <- 1:4 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, theta = NULL, D = 4, CL = TRUE, detectfn = 'HHN', details = details, seed = 347, type = 'BVN', lambda = 0.4) for (i in 1:4) { args <- c(baseargs, list(sigmaX = 25/i^0.5, sigmaY = 25*i^0.5)) simrandomBVNs100.1[[i]] <- do.call(runEllipseSim, args) } save(simrandomBVNs100.1, file = paste0(simfolder, 'simrandomBVNs100.1.RData'))
tr <- make.grid(36, 1, spacing = 25, detector = 'proximity') simgridalignlinw1.1 <- vector('list', 6) names(simgridalignlinw1.1) <- c(1:4,'1-4','1,4') simgridalignlinw4.1 <- simgridalignlinw3.1 <- simgridalignlinw2.1 <- simgridalignlinw1.1 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', theta = 0, CL = TRUE, detectfn = 'HHN', details = details) for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) # bvn perpendicular to traps args$theta <- 0; simgridalignlin1.3[[i]] <- do.call(runEllipseSim, args) # bvn oblique to traps args$theta <- pi/4; simgridalignlin2.3[[i]] <- do.call(runEllipseSim, args) # bvn parallel to traps args$theta <- pi/2; simgridalignlin3.3[[i]] <- do.call(runEllipseSim, args) # bvn random orientation args$theta <- NULL; simgridalignlin4.3[[i]] <- do.call(runEllipseSim, args) } args <- c(baseargs, list(sigmaX = rs2xy)) args$theta <- 0; simgridalignlin1.3[[5]] <- do.call(runEllipseSim, args) args$theta <- pi/4; simgridalignlin2.3[[5]] <- do.call(runEllipseSim, args) args$theta <- pi/2; simgridalignlin3.3[[5]] <- do.call(runEllipseSim, args) args$theta <- NULL; simgridalignlin4.3[[5]] <- do.call(runEllipseSim, args) args <- c(baseargs, list(sigmaX = rs2xy2)) args$theta <- 0; simgridalignlin1.3[[6]] <- do.call(runEllipseSim, args) args$theta <- pi/4; simgridalignlin2.3[[6]] <- do.call(runEllipseSim, args) args$theta <- pi/2; simgridalignlin3.3[[6]] <- do.call(runEllipseSim, args) args$theta <- NULL; simgridalignlin4.3[[6]] <- do.call(runEllipseSim, args) save(simgridalignlin1.3, file = paste0(simfolder, 'simgridalignlin1.3.RData')) save(simgridalignlin2.3, file = paste0(simfolder, 'simgridalignlin2.3.RData')) save(simgridalignlin3.3, file = paste0(simfolder, 'simgridalignlin3.3.RData')) save(simgridalignlin4.3, file = paste0(simfolder, 'simgridalignlin4.3.RData'))
tr <- make.grid(10, 10, spacing = 25, detector = 'proximity', hollow = TRUE) simgridalignsq1.1 <- vector('list', 4) # 'w' for wide names(simgridalignsq1.1) <- 1:4 simgridalignsq3.1 <- simgridalignsq2.1 <- simgridalignsq1.1 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = TRUE, detectfn = 'HHN', details = details) for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) # bvn aligned to traps args$theta <- 0; simgridalignsq1.1[[i]] <- do.call(runEllipseSim, args) # bvn oblique to traps args$theta <- pi/4; simgridalignsq2.1[[i]] <- do.call(runEllipseSim, args) # bvn random orientation args$theta <- NULL; simgridalignsq3.1[[i]] <- do.call(runEllipseSim, args) } args <- c(baseargs, list(sigmaX = rs2xy)) args$theta <- 0; simgridalignsq1.1[[5]] <- do.call(runEllipseSim, args) args$theta <- pi/4; simgridalignsq2.1[[5]] <- do.call(runEllipseSim, args) args$theta <- NULL; simgridalignsq3.1[[5]] <- do.call(runEllipseSim, args) args <- c(baseargs, list(sigmaX = rs2xy2)) args$theta <- 0; simgridalignsq1.1[[6]] <- do.call(runEllipseSim, args) args$theta <- pi/4; simgridalignsq2.1[[6]] <- do.call(runEllipseSim, args) args$theta <- NULL; simgridalignsq3.1[[6]] <- do.call(runEllipseSim, args) save(simgridalignsq1.2, file = paste0(simfolder, 'simgridalignsq1.2.RData')) save(simgridalignsq2.2, file = paste0(simfolder, 'simgridalignsq2.2.RData')) save(simgridalignsq3.2, file = paste0(simfolder, 'simgridalignsq3.2.RData'))
load(file = paste0(simfolder, 'simrandomBVN100.3.RData')) load(file = paste0(simfolder, 'simrandomBVNs100.1.RData')) par(mfrow=c(1,2)) simplot(list(BVN = simrandomBVNs100.1), legend = FALSE) mtext(side=3, text = '25 m') simplot(list(BVN = simrandomBVN100.3[1:4]), legend = FALSE) mtext(side=3, text = '50 m')
load(file = paste0(simfolder, 'simgridalignlin1.3.RData')) load(file = paste0(simfolder, 'simgridalignlin2.3.RData')) load(file = paste0(simfolder, 'simgridalignlin3.3.RData')) load(file = paste0(simfolder, 'simgridalignlin4.3.RData')) load(file = paste0(simfolder, 'simgridalignlinw1.1.RData')) load(file = paste0(simfolder, 'simgridalignlinw2.1.RData')) load(file = paste0(simfolder, 'simgridalignlinw3.1.RData')) load(file = paste0(simfolder, 'simgridalignlinw4.1.RData')) par(mfrow=c(1,2)) simplot(list(Perpendicular = simgridalignlin1.3[1:4], Oblique = simgridalignlin2.3[1:4], Parallel = simgridalignlin3.3[1:4], Random = simgridalignlin4.3[1:4]), legend = TRUE, ylim = c(-1,3.5)) mtext(side=3, text = '25 m') simplot(list(Perpendicular = simgridalignlinw1.1[1:4], Oblique = simgridalignlinw2.1[1:4], Parallel = simgridalignlinw3.1[1:4], Random = simgridalignlinw4.1[1:4]), legend = TRUE, ylim = c(-1,3.5)) mtext(side=3, text = '50 m')
load(file = paste0(simfolder, 'simgridalignsq1.2.RData')) load(file = paste0(simfolder, 'simgridalignsq2.2.RData')) load(file = paste0(simfolder, 'simgridalignsq3.2.RData')) load(file = paste0(simfolder, 'simgridalignsqw1.1.RData')) load(file = paste0(simfolder, 'simgridalignsqw2.1.RData')) load(file = paste0(simfolder, 'simgridalignsqw3.1.RData')) par(mfrow =c(1,2)) simplot(list(Aligned = simgridalignsq1.2[1:4], Oblique = simgridalignsq2.2[1:4], Random = simgridalignsq3.2[1:4]), legend = TRUE, ylim = c(-1,3.5)) mtext(side=3, text = '25 m') simplot(list(Aligned = simgridalignsqw1.1[1:4], Oblique = simgridalignsqw2.1[1:4], Random = simgridalignsqw3.1[1:4]), legend = TRUE, ylim = c(-1,3.5)) mtext(side=3, text = '50 m')
tr <- make.grid(10, 10, spacing = 50, detector = 'proximity') simrandom100low.2 <- vector('list', 4) names(simrandom100low.2) <- 1:4 simrandomBVN100low.2 <- simrandom100low.2 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, D = 0.5, CL = TRUE, detectfn = 'HHN', details = details) for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) args$type <- 'uniform'; args$g0 <- 0.2 simrandom100low.2[[i]] <- do.call(runEllipseSim, args) args$type <- 'BVN'; args$lambda0 <- 0.4 simrandomBVN100low.2[[i]] <- do.call(runEllipseSim, args) } save(simrandom100low.2, file = paste0(simfolder, 'simrandom100low.2.RData')) save(simrandomBVN100low.2, file = paste0(simfolder, 'simrandomBVN100low.2.RData'))
load(file = paste0(simfolder, 'simrandom100low.2.RData')) load(file = paste0(simfolder, 'simrandomBVN100low.2.RData')) par(mfrow = c(1,1), mar = c(4,4,4,4), xpd = FALSE) simplot(list(Uniform = simrandom100low.2, BVN = simrandomBVN100low.2), trueval = 0.5, legend = TRUE) simsum(list(Uniform = simrandom100low.2, BVN = simrandomBVN100low.2), trueval = 0.5, compact = NULL)
tr <- make.grid(6, 6, spacing = 50, detector = 'proximity') simrandom36.2 <- vector('list', 4) names(simrandom36.2) <- 1:4 simrandomBVN36.2 <- simrandom36.2 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, D = 4, CL = TRUE, detectfn = 'HHN', details = details) for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) args$type <- 'uniform'; args$g0 <- 0.2 simrandom36.2[[i]] <- do.call(runEllipseSim, args) args$type <- 'BVN'; args$lambda0 <- 0.4 simrandomBVN36.2[[i]] <- do.call(runEllipseSim, args) } save(simrandom36.2, file = paste0(simfolder, 'simrandom36.2.RData')) save(simrandomBVN36.2, file = paste0(simfolder, 'simrandomBVN36.2.RData'))
load(file = paste0(simfolder, 'simrandom36.2.RData')) load(file = paste0(simfolder, 'simrandomBVN36.2.RData')) par(mfrow = c(1,1), mar = c(4,4,4,4), xpd = FALSE) simplot(list(Uniform = simrandom36.2[1:4], BVN = simrandomBVN36.2[1:4]), legend = TRUE) simsum(list(Uniform = simrandom36.2, BVN = simrandomBVN36.2))
tr <- make.grid(10, 10, spacing = 50, detector = 'proximity') simrandom100C.2 <- vector('list', 4) names(simrandom100C.2) <- 1:4 simrandomBVN100C.2 <- simrandom100C.2 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, theta = -1, D = 4, CL = TRUE, detectfn = 'HHN', details = details) for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) args$type <- 'uniform'; args$g0 <- 0.2 simrandom100C.2[[i]] <- do.call(runEllipseSim, args) args$type <- 'BVN'; args$lambda0 <- 0.4 simrandomBVN100C.2[[i]] <- do.call(runEllipseSim, args) } save(simrandom100C.2, file = paste0(simfolder, 'simrandom100C.2.RData')) save(simrandomBVN100C.2, file = paste0(simfolder, 'simrandomBVN100C.2.RData'))
load(file = paste0(simfolder, 'simrandom100C.2.RData')) load(file = paste0(simfolder, 'simrandomBVN100C.2.RData')) simplot(list(Uniform = simrandom100C.2, BVN = simrandomBVN100C.2), legend = TRUE)
Spacings narrow and wide (w); orientation parallel (p) and oblique(o). Only consider aspect ratio 4 for now.
tr <- make.grid(10,10, spacing = 25, detector = 'proximity') simrandomBVN100p.1 <- vector('list', 4) names(simrandomBVN100p.1) <- 1:4 simrandomBVN100o.1 <- simrandomBVN100p.1 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, type = 'BVN', lambda0 = 0.4, D = 4, CL = TRUE, detectfn = 'HHN', details = details) for (i in 4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY, theta = 0)) simrandomBVN100p.1[[i]] <- do.call(runEllipseSim, args) args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY, theta = pi/4)) simrandomBVN100o.1[[i]] <- do.call(runEllipseSim, args) } tr <- make.grid(10,10, spacing = 50, detector = 'proximity') simrandomBVN100pw.1 <- vector('list', 4) names(simrandomBVN100pw.1) <- 1:4 simrandomBVN100ow.1 <- simrandomBVN100pw.1 for (i in 4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY, theta = 0)) simrandomBVN100pw.1[[i]] <- do.call(runEllipseSim, args) args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY, theta = pi/4)) simrandomBVN100ow.1[[i]] <- do.call(runEllipseSim, args) } save(simrandomBVN100p.1, file = paste0(simfolder, 'simrandomBVN100p.1.RData')) save(simrandomBVN100o.1, file = paste0(simfolder, 'simrandomBVN100o.1.RData')) save(simrandomBVN100pw.1, file = paste0(simfolder, 'simrandomBVN100pw.1.RData')) save(simrandomBVN100ow.1, file = paste0(simfolder, 'simrandomBVN100ow.1.RData'))
Using function anisotropic.fit
. Extract the fitted coefficients corresponding to psiA
and psiR
with the coef
method for secr objects, and rely on direct estimation of density (CL = FALSE).
tr <- make.grid(10, 10, spacing = 25, detector = 'proximity', hollow = TRUE) simaniso1 <- vector('list', 4); names(simaniso1) <- 1:4 simaniso4 <- simaniso3 <- simaniso2 <- simaniso1 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = FALSE, detectfn = 'HHN', details = list(distribution = 'binomial'), extractfn = coef, secrfn = 'anisotropic.fit') for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) args$theta <- 0 # aligned to detectors simaniso1[[i]] <- do.call(runEllipseSim, args) args$theta <- pi/4 # oblique to detectors simaniso2[[i]] <- do.call(runEllipseSim, args) args$theta <- NULL # random orientation simaniso3[[i]] <- do.call(runEllipseSim, args) args$theta <- -1 # common random orientation simaniso4[[i]] <- do.call(runEllipseSim, args) } save(simaniso1, file = paste0(simfolder, 'simaniso1.RData')) save(simaniso2, file = paste0(simfolder, 'simaniso2.RData')) save(simaniso3, file = paste0(simfolder, 'simaniso3.RData')) save(simaniso4, file = paste0(simfolder, 'simaniso4.RData'))
load(file = paste0(simfolder, 'simgridalignsq1.2.RData')) load(file = paste0(simfolder, 'simgridalignsq2.2.RData')) load(file = paste0(simfolder, 'simgridalignsq3.2.RData')) load(file = paste0(simfolder, 'simaniso1.RData')) load(file = paste0(simfolder, 'simaniso2.RData')) load(file = paste0(simfolder, 'simaniso3.RData')) load(file = paste0(simfolder, 'simaniso4.RData')) par(mfrow = c(1,2)) simplot(list(Aligned = simgridalignsq1.2[1:4], Oblique = simgridalignsq2.2[1:4], Random = simgridalignsq3.2[1:4]), ylim = c(-0.4,0.45), legend = TRUE, pchi = c(23, 22, 16)) text(-0.4, 0.5, 'a.', cex = 1.45, xpd = TRUE) text(3, -0.35, 'isotropic') simplot(list(Aligned = simaniso1, Oblique = simaniso2, Random = simaniso3), component = "pred", ylim = c(-0.4,0.45), legend = TRUE, pchi = c(23, 22, 16)) text(-0.4, 0.5, 'b.', cex = 1.45, xpd = TRUE) text(3, -0.35, 'anisotropic')
Tabular summary for anisotropic model.
simsum(list(Aligned = simaniso1, Oblique = simaniso2, Random = simaniso3), component = "pred")
Tabular comparison of isotropic and anisotropic models for hollow grid.
tab1 <- simsum(list(Aligned = simgridalignsq1.2[1:4], Oblique = simgridalignsq2.2[1:4], Random = simgridalignsq3.2[1:4])) # find D-hat in the 'fit' component tab2 <- simsum(list(Aligned = simaniso1, Oblique = simaniso2, Random = simaniso3), component = "pred") # find D-hat in the 'pred' component fn <- function(t1,t2) cbind(t1[,-1],t2[,-1]) mapply(fn, tab1, tab2, SIMPLIFY = FALSE)
tr <- make.grid(36, 1, spacing = 25, detector = 'proximity') simanisolin1 <- vector('list', 4); names(simanisolin1) <- 1:4 simanisolin5 <- simanisolin4 <- simanisolin3 <- simanisolin2 <- simanisolin1 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = FALSE, detectfn = 'HHN', details = list(distribution = 'binomial'), extractfn = coef, secrfn = 'anisotropic.fit') for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) args$theta <- 0 # perpendicular to detectors simanisolin1[[i]] <- do.call(runEllipseSim, args) args$theta <- pi/4 # oblique to detectors simanisolin2[[i]] <- do.call(runEllipseSim, args) args$theta <- NULL # random orientation simanisolin3[[i]] <- do.call(runEllipseSim, args) args$theta <- -1 # common random orientation simanisolin4[[i]] <- do.call(runEllipseSim, args) args$theta <- pi/2 # parallel simanisolin5[[i]] <- do.call(runEllipseSim, args) # partial run 2018-08-23 } save(simanisolin1, file = paste0(simfolder, 'simanisolin1.RData')) save(simanisolin2, file = paste0(simfolder, 'simanisolin2.RData')) save(simanisolin3, file = paste0(simfolder, 'simanisolin3.RData')) save(simanisolin4, file = paste0(simfolder, 'simanisolin4.RData')) save(simanisolin5, file = paste0(simfolder, 'simanisolin5.RData'))
load(file = paste0(simfolder, 'simgridalignlin1.3.RData')) load(file = paste0(simfolder, 'simgridalignlin2.3.RData')) load(file = paste0(simfolder, 'simgridalignlin4.3.RData')) load(file = paste0(simfolder, 'simanisolin1.RData')) load(file = paste0(simfolder, 'simanisolin2.RData')) load(file = paste0(simfolder, 'simanisolin3.RData')) par(mfrow = c(1,2)) simplot(list(Aligned = simgridalignlin1.3[1:4], Oblique = simgridalignlin2.3[1:4], Random = simgridalignlin4.3[1:4]), ylim = c(-1,3.5), legend = TRUE, pchi = c(23, 22, 16)) text(-0.1, 3.5, 'a.', cex = 1.45, xpd = TRUE) text(1, -0.8, 'isotropic', adj=0) simplot(list(Perpendicular = simanisolin1, Oblique = simanisolin2, Random = simanisolin3), component = "pred", ylim = c(-1,3.5), legend = TRUE, pchi = c(23, 22, 16)) text(-0.1, 3.5, 'b.', cex = 1.45, xpd = TRUE) text(1, -0.8, 'anisotropic', adj=0)
Conclude that anisotropic transformation largely ineffective with linear array.
tr <- make.grid(10, 10, spacing = 25, detector = 'proximity') simanisogrid1 <- vector('list', 4); names(simanisogrid1) <- 1:4 simanisogrid4 <- simanisogrid3 <- simanisogrid2 <- simanisogrid1 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = FALSE, detectfn = 'HHN', details = list(distribution = 'binomial'), extractfn = coef, secrfn = 'anisotropic.fit') for (i in 1:4) { sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 args <- c(baseargs, list(sigmaX = sigmaX, sigmaY = sigmaY)) args$theta <- 0 # aligned to detectors simanisogrid1[[i]] <- do.call(runEllipseSim, args) args$theta <- pi/4 # oblique to detectors simanisogrid2[[i]] <- do.call(runEllipseSim, args) args$theta <- NULL # random orientation simanisogrid3[[i]] <- do.call(runEllipseSim, args) args$theta <- -1 # common random orientation simanisogrid4[[i]] <- do.call(runEllipseSim, args) } save(simanisogrid1, file = paste0(simfolder, 'simanisogrid1.RData')) save(simanisogrid2, file = paste0(simfolder, 'simanisogrid2.RData')) save(simanisogrid3, file = paste0(simfolder, 'simanisogrid3.RData')) save(simanisogrid4, file = paste0(simfolder, 'simanisogrid4.RData'))
load(file = paste0(simfolder, 'simanisogrid3.RData')) load(file = paste0(simfolder, 'simrandomBVNs100.1.RData')) simsum(list(Simple = simrandomBVNs100.1)) simsum(list(Anisotropic = simanisogrid3), component = "pred")
Performs well with full grid data (narrow spacing).
Now run all geometries with wide spacing (w) for aspect ratio 4.0.
tr <- make.grid(10, 10, spacing = 50, detector = 'proximity') simanisogridw1 <- vector('list', 4); names(simanisogridw1) <- 1:4 simanisogridw3 <- simanisogridw2 <- simanisogridw1 i <- 4 sigmaX <- 25/i^0.5; sigmaY <- 25*i^0.5 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = FALSE, detectfn = 'HHN', details = list(distribution = 'binomial'), extractfn = coef, secrfn = 'anisotropic.fit', sigmaX = sigmaX, sigmaY = sigmaY) args <- c(baseargs, list(theta = 0)) # aligned to detectors simanisogridw1[[i]] <- do.call(runEllipseSim, args) args$theta <- pi/4 # oblique to detectors simanisogridw2[[i]] <- do.call(runEllipseSim, args) args$theta <- NULL # random orientation simanisogridw3[[i]] <- do.call(runEllipseSim, args) save(simanisogridw1, file = paste0(simfolder, 'simanisogridw1.RData')) save(simanisogridw2, file = paste0(simfolder, 'simanisogridw2.RData')) save(simanisogridw3, file = paste0(simfolder, 'simanisogridw3.RData')) tr <- make.grid(36, 1, spacing = 50, detector = 'proximity') baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = FALSE, detectfn = 'HHN', details = list(distribution = 'binomial'), extractfn = coef, secrfn = 'anisotropic.fit', sigmaX = sigmaX, sigmaY = sigmaY) simanisolinw1 <- vector('list', 4); names(simanisolinw1) <- 1:4 simanisolinw4 <- simanisolinw3 <- simanisolinw2 <- simanisolinw1 args <- c(baseargs, list(theta = 0)) # perpendicular to detectors simanisolinw1[[i]] <- do.call(runEllipseSim, args) args$theta <- pi/4 # oblique to detectors simanisolinw2[[i]] <- do.call(runEllipseSim, args) args$theta <- pi/2 # parallel to detectors simanisolinw3[[i]] <- do.call(runEllipseSim, args) args$theta <- NULL # random orientation simanisolinw4[[i]] <- do.call(runEllipseSim, args) save(simanisolinw1, file = paste0(simfolder, 'simanisolinw1.RData')) save(simanisolinw2, file = paste0(simfolder, 'simanisolinw2.RData')) save(simanisolinw3, file = paste0(simfolder, 'simanisolinw3.RData')) save(simanisolinw4, file = paste0(simfolder, 'simanisolinw4.RData')) tr <- make.grid(10, 10, spacing = 50, detector = 'proximity', hollow = TRUE) simanisow1 <- vector('list', 4); names(simanisow1) <- 1:4 simanisow3 <- simanisow2 <- simanisow1 baseargs <- list(nrepl = nrepl, buffer = 200, ncores = ncores, traps = tr, lambda0 = 0.4, D = 4, type = 'BVN', CL = FALSE, detectfn = 'HHN', details = list(distribution = 'binomial'), extractfn = coef, secrfn = 'anisotropic.fit', sigmaX = sigmaX, sigmaY = sigmaY) args <- c(baseargs, list(theta = 0)) # aligned to detectors simanisow1[[i]] <- do.call(runEllipseSim, args) args$theta <- pi/4 # oblique to detectors simanisow2[[i]] <- do.call(runEllipseSim, args) args$theta <- NULL # random orientation simanisow3[[i]] <- do.call(runEllipseSim, args) save(simanisow1, file = paste0(simfolder, 'simanisow1.RData')) save(simanisow2, file = paste0(simfolder, 'simanisow2.RData')) save(simanisow3, file = paste0(simfolder, 'simanisow3.RData'))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.