sqpad: Single-bin QPAD (SQPAD)

In detect: Analyzing Wildlife Data with Detection Error

Single-bin QPAD (SQPAD)

Description

Single bin QPAD (SQPAD) approach by Lele and Solymos (2025).

Usage

sqpad(formula, data, dis, dur,
    type = c("full", "approx"), det = c("joint", "pq"),
    Nmax = NULL, K = NULL, A = NULL,
    montecarlo = FALSE, np = 1000, distcorr = 2/3, ...)

sqpad.fit(Y, dis, dur, X = NULL, Z = NULL, 
    type = c("full", "approx", "conv"), det = c("joint", "pq"),
    init = NULL, method = "Nelder-Mead", hessian = TRUE,
    tt1 = NULL, A = NULL, Nmax = NULL, K = NULL, 
    montecarlo = FALSE, np = 1000, distcorr = 2/3, dislist = NULL, ...)

## S3 method for class 'sqpad'
print(x, digits, ...)
## S3 method for class 'sqpad'
vcov(object, ...)
## S3 method for class 'sqpad'
fitted(object, ...)
## S3 method for class 'sqpad'
logLik(object, ...)
## S3 method for class 'sqpad'
summary(object, ...)
## S3 method for class 'summary.sqpad'
print(x, digits, ...)
## S3 method for class 'sqpad'
predict(object, newdata = NULL, 
    type = c("link", "response"), ...)

rtriang(n, r = 1)

Arguments

formula	formula, LHS takes a vector, RHS is either 1 or some covariates, see Examples.
data	data.
dis, dur	a vector with distance and duration values.
A	a vector with area values. Area is calculated from dis, use A if distance is not directly related to area (e.g. for ARU data).
type	character, one of "full" (full likelihood), "approx" (approximate method), and "conv" (convolution likelihood).
det	character, the detection model: "joint" or independence ("pq").
Y	response vector, non-negative integers.
X, Z	design matrices, X is the matrix with covariates for the density parameters. Z is the matrix with covariates for the cue rate model. NULL means constant density and detection.
K	truncation value, K=1 implements the occupancy version of SQPAD. K>1 implements the truncated count models.
Nmax	maximum abundance value for numerical integration under the full and convolution likelihoods.
montecarlo, np, distcorr	should mean probability (TRUE) using np points or approximation be used with distance correction distcorr in the detection model.
dislist	distance list for the convolution likelihood.
tt1	vector with time to 1st detection values, same units as for dur.
init	optional initial values for optim.
method	method for optim.
hessian	logical, should the Hessian matrix be computed by optim.
object, x	fitted model object.
newdata	optionally, a data frame in which to look for variables with which to predict. If omitted, the fitted linear predictors are used.
digits	digits to use when providing summaries.
...	additional options for optim or the methods.
n	number of random variates to generate from the triangular distribution.
r	maximum value for the triangular distribution.

Details

Single bin QPAD (SQPAD) approach for robust analysis of point count data with detection error by Lele and Solymos (2025).

Value

An object of class 'sqpad'.

Author(s)

Peter Solymos and Subhash Lele

References

Solymos, P., Lele, S. R., 2025. Single bin QPAD (SQPAD) approach for robust analysis of point count data with detection error. Ornithological Applications, xx, xx–xx. <doi:10.1093/ornithapp/duaf078>

Supporting info for the SQPAD method: https://github.com/psolymos/sqpad-paper, <doi:10.5281/zenodo.16172209>

Examples

set.seed(0)
n <- 100
x <- rnorm(n)
D <- exp(-2 + 0.5 * x)
phi <- 0.25
tau <- 1
dur <- sample(1:10, n, replace=TRUE)
dis <- sample(seq(0.5, 2, 0.25), n, replace=TRUE)
A <- dis^2 * pi
dcorr <- 2/3
p <- 1 - exp(-dur * phi * exp(-(dis*dcorr)^2/tau^2))
N <- rpois(n, D*A)
Y <- rbinom(n, N, p)

df <- data.frame(x = x, y = Y)

m <- sqpad(y ~ x | 1, data = df, dis = dis, dur = dur, type = "full", det = "joint", K = NULL)

print(m)
summary(m)

coef(m)
nobs(m)
vcov(m)
confint(m)

logLik(m)
AIC(m)
BIC(m)

fitted(m)
predict(m, type = "link")
predict(m, newdata = df[1:10,], type = "response")
predict(m, newdata = df[1:10,], type = "link")

## Not run: 
m0 <- sqpad(y ~ 1 | 1, data = df, dis = dis, dur = dur, type = "full", det = "joint", K = NULL)
m1 <- sqpad(y ~ 1 | x, data = df, dis = dis, dur = dur, type = "full", det = "joint", K = NULL)
m2 <- sqpad(y ~ x | x, data = df, dis = dis, dur = dur, type = "full", det = "joint", K = NULL)

AIC(m, m0, m1, m2) # m2 is best
BIC(m, m0, m1, m2) # this is needed!

# average probability

set.seed(5)
n <- 1000
x <- rnorm(n)
D <- exp(-2 + 0.5 * x)
phi <- 0.25
tau <- 1
dur <- sample(1:10, n, replace=TRUE)
dis <- sample(seq(0.5, 2, 0.25), n, replace=TRUE)
A <- dis^2 * pi

dcorr <- 2/3
p <- 1 - exp(-dur * phi * exp(-(dis*dcorr)^2/tau^2))
N <- rpois(n, D*A)
Y <- rbinom(n, N, p)
df <- data.frame(x = x, y = Y)

np <- 1000
dij <- rtriang(as.integer(np), 1)
dij <- matrix(dij, nrow = n, ncol = np, byrow = TRUE)
p2 <- rowMeans(1 - exp(-dur * phi * exp(-(dis*dij)^2/tau^2)))
Y2 <- rbinom(n, N, p2)
df2 <- data.frame(x = x, y = Y2)

plot(p, p2)
abline(0, 1)

mAP <- sqpad(y ~ x | 1, data = df, dis = dis, dur = dur, type = "full", 
    det = "joint", K = NULL, montecarlo = FALSE, hessian = FALSE)
mMC <- sqpad(y ~ x | 1, data = df, dis = dis, dur = dur, type = "full", 
    det = "joint", K = NULL, montecarlo = TRUE, hessian = FALSE)
m2AP <- sqpad(y ~ x | 1, data = df2, dis = dis, dur = dur, type = "full", 
    det = "joint", K = NULL, montecarlo = FALSE, hessian = FALSE)
m2MC <- sqpad(y ~ x | 1, data = df2, dis = dis, dur = dur, type = "full", 
    det = "joint", K = NULL, montecarlo = TRUE, hessian = FALSE)

cbind(mAP=coef(mAP), mMC=coef(mMC), m2AP=coef(m2AP), m2MC=coef(m2MC))

# convolution likelihood

D <- 1
phi <- 0.25
tau <- 1.2
n <- 200

dur <- sample(1:10, n, replace=TRUE)
dis <- sample(seq(0.25, 2, 0.25), n, replace=TRUE)
A <- dis^2 * pi
N <- rpois(n, D*A)
dislist <- lapply(seq_len(n), function(i) {
    dij <- rtriang(N[i], dis[i])
    # pij is based on time of 1st detection and distance at 1st detection
    pij <- 1 - exp(-dur[i] * phi * exp(-(dij)^2/tau^2))
    Yij <- rbinom(length(pij), 1, pij)
    dij[Yij > 0]
})
Y <- sapply(dis, length)

m4 <- sqpad.fit(Y = Y, dis = dis, dur = dur, dislist = dislist, 
    type = "conv", det = "joint", Nmax = 20)
cbind(parameters = c(D = D, phi = phi, tau = tau), estimates = exp(m4$coef))

## End(Not run)

detect documentation built on Jan. 8, 2026, 5:07 p.m.