pcount.spHDS.RdFunction fits an N-mixture model for a discrete state space with raster covariates, and a detection function which decreases with distance from the observer, assumed to be at the centre. See Kery & Royle (2016) Section 9.8.4 for details.
pcount.spHDS(formula, data, K, mixture = c("P", "NB", "ZIP"), starts,
method = "BFGS", se = TRUE, ...)Double right-hand side formula describing covariates of detection and abundance, in that order.
Detection model should be specified without an intercept, for example:
~ -1 + I(dist^2), where dist is a covariate giving the distance of each cell of the raster from the observer.
Internally this forces the intercept p(0) = 1, conventional for distance
sampling models (see Kery & Royle (2016) for explanation). More general models
work but may not honor that constraint. e.g.,
~ 1,
~ dist,
~ I(dist^2),
~ dist + I(dist^2)
an unmarkedFramePCount object supplying data to the model.
Integer upper index of integration for N-mixture. This should be set high enough so that it does not affect the parameter estimates. Note that computation time will increase with K.
character specifying mixture: Poisson (P), Negative-Binomial (NB), or Zero Inflated Poisson (ZIP).
vector of starting values
Optimization method used by optim.
logical specifying whether or not to compute standard errors.
Additional arguments to optim, such as lower and upper bounds
unmarkedFit object describing the model fit.
Kery & Royle (2016) Applied Hierarachical Modeling in Ecology Section 9.8.4
## Simulate some data to analyse
# This is based on Kery and Royle (2016) section 9.8.3
# See AHMbook::sim.spatialDS for more simulation options.
# We will simulate distance data for a logit detection function with sigma = 1,
# for a 6x6 square, divided into a 30 x 30 grid of pixels (900 in all), with the
# observer in the centre.
set.seed(2017)
## 1. Create coordinates for 30 x 30 grid
grx <- seq(0.1, 5.9, 0.2) # mid-point coordinates
gr <- expand.grid(grx, grx) # data frame with coordinates of pixel centres
## 2a. Simulate spatially correlated Habitat covariate
# Get the pair-wise distances between pixel centres
tmp <- as.matrix(dist(gr)) # a 900 x 900 matrix
# Correlation is a negative exponential function of distance, with scale parameter = 1
V <- exp(-tmp/1)
Habitat <- crossprod(t(chol(V)), rnorm(900))
## 2b. Do a detection covariate: the distance of each pixel centre from the observer
dist <- sqrt((gr[,1]-3)^2 + (gr[,2]-3)^2)
## 3. Simulate the true population
# Probability that an animal is in a pixel depends on the Habitat covariate, with
# coefficient beta:
beta <- 1
probs <- exp(beta*Habitat) / sum(exp(beta*Habitat))
# Allocate 600 animals to the 900 pixels, get the pixel ID for each animal
pixel.id <- sample(1:900, 600, replace=TRUE, prob=probs)
## 4. Simulate the detection process
# Get the distance of each animal from the observer
# (As an approximation, we'll treat animals as if they are at the pixel centre.)
d <- dist[pixel.id]
# Calculate probability of detection with logit detection function with
sigma <- 1
p <- 2*plogis(-d^2/(2*sigma^2))
# Simulate the 1/0 detection/nondetection vector
y <- rbinom(600, 1, p)
# Check the number of animals detected
sum(y)
#> [1] 112
# Select the pixel IDs for the animals detected and count the number in each pixel
detected.pixel.id <- pixel.id[y == 1]
pixel.count <- tabulate(detected.pixel.id, nbins=900)
## 5. Prepare the data for unmarked
# Centre the Habitat covariate
Habitat <- Habitat - mean(Habitat)
# Construct the unmarkedFramePCount object
umf <- unmarkedFramePCount(y=cbind(pixel.count), # y needs to be a 1-column matrix
siteCovs=data.frame(dist=dist, Habitat=Habitat))
summary(umf)
#> unmarkedFrame Object
#>
#> 900 sites
#> Maximum number of observations per site: 1
#> Mean number of observations per site: 1
#> Sites with at least one detection: 87
#>
#> Tabulation of y observations:
#> 0 1 2 3 5
#> 813 67 17 2 1
#>
#> Site-level covariates:
#> dist Habitat
#> Min. :0.1414 Min. :-3.084964
#> 1st Qu.:1.7029 1st Qu.:-0.621886
#> Median :2.4042 Median :-0.000016
#> Mean :2.2946 Mean : 0.000000
#> 3rd Qu.:2.9155 3rd Qu.: 0.570115
#> Max. :4.1012 Max. : 2.783497
## 6. Fit some models
(fm0 <- pcount.spHDS(~ -1 + I(dist^2) ~ 1, umf, K = 20))
#>
#> Call:
#> pcount.spHDS(formula = ~-1 + I(dist^2) ~ 1, data = umf, K = 20)
#>
#> Abundance:
#> Estimate SE z P(>|z|)
#> -0.569 0.131 -4.35 1.37e-05
#>
#> Detection:
#> Estimate SE z P(>|z|)
#> -0.602 0.0929 -6.48 9.14e-11
#>
#> AIC: 578.9329
(fm1 <- pcount.spHDS(~ -1 + I(dist^2) ~ Habitat, umf, K = 20))
#>
#> Call:
#> pcount.spHDS(formula = ~-1 + I(dist^2) ~ Habitat, data = umf,
#> K = 20)
#>
#> Abundance:
#> Estimate SE z P(>|z|)
#> (Intercept) -0.897 0.139 -6.47 1.01e-10
#> Habitat 1.134 0.134 8.46 2.66e-17
#>
#> Detection:
#> Estimate SE z P(>|z|)
#> -0.616 0.0821 -7.5 6.59e-14
#>
#> AIC: 505.6613
# The true Habitat coefficient (beta above) = 1
# fm1 has much lower AIC; look at the population estimate
sum(predict(fm1, type="state")[, 1])
#> [1] 631.1605