\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}

Generating and plotting elliptical home ranges

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').

Simulating elliptical detection data

Generating detection histories

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.

Wrapper function to generate BVN data and fit SECR model

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.

Heterogeneous elliptical home ranges {#heterogeneity}

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') 

Anisotropic home ranges: a partial solution {#anisotropic}

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.

Distances in transformed space

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
}

Example

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)).

Function to fit anisotropic model

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

Fixed orientation

If the orientation is known then there are two solutions:

  1. Re-write anisodistfn for miscparm of length 1, with hard-coded value for psiA, or
  2. Fix the beta parameter corresponding to psiA.

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)

Inadequate data

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).

Package limitations

This package has the limited goal of determining how range elongation affects estimates of density in simple SECR models, and these specific limitations:

  1. Only binary proximity detectors are supported.
  2. The spatial distribution of activity centres is assumed to be homogeneous Poisson.
  3. Ellipses are specified using either 'sigmaX' and 'sigmaY' as separate arguments (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.

References

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

Appendix. Code for simulations of Efford (in prep) {#Appendix}

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)
}

Main simulations

10 x 10 grid

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

Straight line 36-detectors

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))

Hollow square 36 detectors

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))

Reduced spacing ($\sigma$ instead of $2\sigma$)

10 x 10 grid

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'))

Straight line

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'))

Hollow square

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'))

Spacing comparisons

10 x 10 grid

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')

Straight line

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')

Hollow square

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')

Variations

Low density (0.5/ha)

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)

Small array (6 x 6 grid)

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))

Common random orientation

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)

Grid with varying orientation

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'))

Anisotropic model

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).

Hollow square

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)

Straight line

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.

10 x 10 grid

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'))


MurrayEfford/secrBVN documentation built on Oct. 14, 2022, 5:20 a.m.