LGCPs - Distance sampling
David Borchers and Finn Lindgren
Generated on 2023-12-03
Source:vignettes/articles/2d_lgcp_distancesampling.Rmd
2d_lgcp_distancesampling.Rmd
Introduction
We’re going to estimate distribution and abundance from a line
transect survey of dolphins in the Gulf of Mexico. These data are also
available in the R
package dsm
(where they go
under the name mexdolphins
). In inlabru
the
data are called mexdolphin
for sp
format, and
mexdolphin_sf
for sf
format.
Get the data
We’ll start by loading the data, renaming it, and extracting the mesh (for convenience).
mexdolphin <- mexdolphin_sf
mesh <- mexdolphin$mesh
Plot the data (the initial code below is just to get rid of tick marks, if desired)
noyticks <- theme(
axis.text.y = element_blank(),
axis.ticks = element_blank()
)
noxticks <- theme(
axis.text.x = element_blank(),
axis.ticks = element_blank()
)
ggplot() +
gg(mexdolphin$ppoly) +
gg(mexdolphin$samplers, color = "grey") +
gg(mexdolphin$points, size = 0.2, alpha = 1) +
theme(legend.key.width = unit(x = 0.2, "cm"), legend.key.height = unit(x = 0.3, "cm")) +
theme(legend.text = element_text(size = 6))
Spatial model with a half-normal detection function
The samplers
in this dataset are lines, not polygons, so
we need to tell inlabru
about the strip half-width,
W
, which in the case of these data is 8. We start by
plotting the distances and histogram of frequencies in distance
intervals:
W <- 8
ggplot(mexdolphin$points) +
geom_histogram(aes(x = distance),
breaks = seq(0, W, length.out = 9),
boundary = 0, fill = NA, color = "black"
) +
geom_point(aes(x = distance), y = 0, pch = "|", cex = 4)
We need to define a half-normal detection probability function. This must take distance as its first argument and the linear predictor of the sigma parameter as its second:
hn <- function(distance, sigma) {
exp(-0.5 * (distance / sigma)^2)
}
To control the prior distribution for the sigma
parameter, we define a transformation function that converts a \(N(0, 1)\) latent variable into an
exponentially distributed variable with expectation 8 (this is a
somewhat arbitrary value, but motivated by the maximum observation
distance W
):
sigma <- function(x) {
bru_forward_transformation(qexp, x, rate = 1 / 8)
}
Specify and fit an SPDE model to these data using a half-normal detection function form. We need to define a (Matérn) covariance function for the SPDE
matern <- inla.spde2.pcmatern(mexdolphin$mesh,
prior.sigma = c(2, 0.01),
prior.range = c(50, 0.01)
)
We need to now separately define the components of the model (the
SPDE, the Intercept and the detection function parameter
sigma_theta
)
cmp <- ~ mySPDE(main = geometry, model = matern) +
sigma_theta(1, prec.linear = 1) + Intercept(1)
… and the formula, which describes how these components are combined to form the linear predictor (remembering that we need an offset due to the unknown direction of the detections!):
From version 2.9.0.9004
, a more compact option is
available, by allowing the component itself do the transformation
between N(0,1) and Exponential(1/8):
# # An alternative approach available from 2.9.0.9004; requires matching
# # adjustments to the later predict() calls, etc.
# cmp <- ~ mySPDE(main = geometry, model = matern) +
# sigma(1,
# prec.linear = 1,
# marginal = bru_mapper_marginal(qexp, pexp, dexp, rate = 1 / 8)
# ) +
# Intercept(1)
# form <- geometry + distance ~ mySPDE +
# log(hn(distance, sigma)) +
# Intercept + log(2)
Then we fit the model, passing both the components and the formula
(previously the formula was constructed invisibly by
inlabru
), and specify integration domains for the spatial
and distance dimensions:
fit <- lgcp(
components = cmp,
mexdolphin$points,
samplers = mexdolphin$samplers,
domain = list(
geometry = mesh,
distance = fm_mesh_1d(seq(0, 8, length.out = 30))
),
formula = form
)
Look at the SPDE parameter posteriors
spde.range <- spde.posterior(fit, "mySPDE", what = "range")
plot(spde.range)
spde.logvar <- spde.posterior(fit, "mySPDE", what = "log.variance")
plot(spde.logvar)
Predict spatial intensity, and plot it:
pxl <- fm_pixels(mesh, dims = c(200, 100), mask = mexdolphin$ppoly)
pr.int <- predict(fit, pxl, ~ exp(mySPDE + Intercept))
ggplot() +
gg(pr.int, geom = "tile") +
gg(mexdolphin$ppoly, linewidth = 1, alpha = 0) +
gg(mexdolphin$samplers, color = "grey") +
gg(mexdolphin$points, size = 0.2, alpha = 1) +
theme(legend.key.width = unit(x = 0.2, "cm"), legend.key.height = unit(x = 0.3, "cm")) +
theme(legend.text = element_text(size = 6))
Predict the detection function and plot it, to generate a plot like
the one below. Here, we should make sure that it doesn’t try to evaluate
the effects of components that can’t be evaluated using the given input
data. From version 2.8.0, inlabru
automatically detects
which components are involved. See ?predict.bru
for more
information.
distdf <- data.frame(distance = seq(0, 8, length.out = 100))
dfun <- predict(fit, distdf, ~ hn(distance, sigma(sigma_theta)))
plot(dfun)
The average detection probability within the maximum detection distance
is estimated to be 0.7060268.
We can look at the posterior for expected number of dolphins as usual:
predpts <- fm_int(mexdolphin$mesh, mexdolphin$ppoly)
Lambda <- predict(fit, predpts, ~ sum(weight * exp(mySPDE + Intercept)))
Lambda
#> mean sd q0.025 q0.5 q0.975 median mean.mc_std_err
#> 1 248.4061 60.10638 163.3518 239.5366 373.7499 239.5366 6.010638
#> sd.mc_std_err
#> 1 6.205917
and including the randomness about the expected number. In this case, it turns out that you need lots of posterior samples, e.g. 2,000 to smooth out the Monte Carlo error in the posterior, and this takes a little while to compute:
Ns <- seq(50, 450, by = 1)
Nest <- predict(fit, predpts,
~ data.frame(
N = Ns,
density = dpois(Ns,
lambda = sum(weight * exp(mySPDE + Intercept))
)
),
n.samples = 2000
)
Nest$plugin_estimate <- dpois(Nest$N, lambda = Lambda$mean)
ggplot(data = Nest) +
geom_line(aes(x = N, y = mean, colour = "Posterior")) +
geom_ribbon(
aes(
x = N,
ymin = mean - 2 * mean.mc_std_err,
ymax = mean + 2 * mean.mc_std_err,
colour = NULL, fill = "Posterior"
),
alpha = 0.2
) +
geom_line(aes(x = N, y = plugin_estimate, colour = "Plugin", fill = "Plugin"))
Hazard-rate Detection Function
Try doing this all again, but use this hazard-rate detection function
model, with the same prior for the sigma
parameter as for
the half-Normal model (such parameters aren’t always comparable, but in
this example it’s a reasonable choice):
hr <- function(distance, sigma) {
1 - exp(-(distance / sigma)^-1)
}
Solution:
formula1 <- geometry + distance ~ mySPDE +
log(hr(distance, sigma(sigma_theta))) +
Intercept + log(2)
fit1 <- lgcp(
components = cmp,
mexdolphin$points,
samplers = mexdolphin$samplers,
domain = list(
geometry = mesh,
distance = fm_mesh_1d(seq(0, 8, length.out = 30))
),
formula = formula1
)
Plots:
spde.range <- spde.posterior(fit1, "mySPDE", what = "range")
plot(spde.range)
spde.logvar <- spde.posterior(fit1, "mySPDE", what = "log.variance")
plot(spde.logvar)
pr.int1 <- predict(fit1, pxl, ~ exp(mySPDE + Intercept))
ggplot() +
gg(pr.int1, geom = "tile") +
gg(mexdolphin$ppoly, linewidth = 1, alpha = 0) +
gg(mexdolphin$samplers, color = "grey") +
gg(mexdolphin$points, size = 0.2, alpha = 1) +
theme(legend.key.width = unit(x = 0.2, "cm"), legend.key.height = unit(x = 0.3, "cm")) +
theme(legend.text = element_text(size = 6))
distdf <- data.frame(distance = seq(0, 8, length.out = 100))
dfun1 <- predict(fit1, distdf, ~ hr(distance, sigma(sigma_theta)))
plot(dfun1)
predpts <- fm_int(mexdolphin$mesh, mexdolphin$ppoly)
Lambda1 <- predict(fit1, predpts, ~ sum(weight * exp(mySPDE + Intercept)))
Lambda1
#> mean sd q0.025 q0.5 q0.975 median mean.mc_std_err
#> 1 301.8794 84.14475 189.6172 286.8156 488.374 286.8156 8.414475
#> sd.mc_std_err
#> 1 6.265036
Ns <- seq(50, 650, by = 1)
Nest1 <- predict(
fit1,
predpts,
~ data.frame(
N = Ns,
density = dpois(Ns,
lambda = sum(weight * exp(mySPDE + Intercept))
)
),
n.samples = 2000
)
Nest1$plugin_estimate <- dpois(Nest1$N, lambda = Lambda1$mean)
ggplot(data = Nest1) +
geom_line(aes(x = N, y = mean, colour = "Posterior")) +
geom_ribbon(
aes(
x = N,
ymin = mean - 2 * mean.mc_std_err,
ymax = mean + 2 * mean.mc_std_err,
colour = NULL, fill = "Posterior"
),
alpha = 0.2
) +
geom_line(aes(x = N, y = plugin_estimate, colour = "Plugin", fill = "Plugin"))
Comparing the models
Look at the goodness-of-fit of the two models in the distance dimension
bc <- bincount(
result = fit,
observations = mexdolphin$points$distance,
breaks = seq(0, max(mexdolphin$points$distance), length.out = 9),
predictor = distance ~ hn(distance, sigma(sigma_theta))
)
attributes(bc)$ggp
bc1 <- bincount(
result = fit1,
observations = mexdolphin$points$distance,
breaks = seq(0, max(mexdolphin$points$distance), length.out = 9),
predictor = distance ~ hn(distance, sigma(sigma_theta))
)
attributes(bc1)$ggp
Fit Models only to the distance sampling data
Half-normal first
formula <- distance ~ log(hn(distance, sigma(sigma_theta))) + Intercept
cmp <- ~ sigma_theta(1, prec.linear = 1) + Intercept(1)
dfit <- lgcp(
components = cmp,
mexdolphin$points,
domain = list(distance = fm_mesh_1d(seq(0, 8, length.out = 30))),
formula = formula,
options = list(bru_initial = list(sigma_theta = 1, Intercept = 3))
)
detfun <- predict(dfit, distdf, ~ hn(distance, sigma(sigma_theta)))
Hazard-rate next
formula1 <- distance ~ log(hr(distance, sigma(sigma_theta))) + Intercept
cmp <- ~ sigma_theta(1, prec.linear = 1) + Intercept(1)
dfit1 <- lgcp(
components = cmp,
mexdolphin$points,
domain = list(distance = fm_mesh_1d(seq(0, 8, length.out = 30))),
formula = formula1
)
detfun1 <- predict(dfit1, distdf, ~ hr(distance, sigma(sigma_theta)))
Plot both lines on histogram of observations. First scale lines to have same area as that of histogram.
Half-normal:
hnline <- data.frame(distance = detfun$distance, p = detfun$mean, lower = detfun$q0.025, upper = detfun$q0.975)
wts <- diff(hnline$distance)
wts[1] <- wts[1] / 2
wts <- c(wts, wts[1])
hnarea <- sum(wts * hnline$p)
n <- length(mexdolphin$points$distance)
scale <- n / hnarea
hnline$En <- hnline$p * scale
hnline$En.lower <- hnline$lower * scale
hnline$En.upper <- hnline$upper * scale
Hazard-rate:
hrline <- data.frame(distance = detfun1$distance, p = detfun1$mean, lower = detfun1$q0.025, upper = detfun1$q0.975)
wts <- diff(hrline$distance)
wts[1] <- wts[1] / 2
wts <- c(wts, wts[1])
hrarea <- sum(wts * hrline$p)
n <- length(mexdolphin$points$distance)
scale <- n / hrarea
hrline$En <- hrline$p * scale
hrline$En.lower <- hrline$lower * scale
hrline$En.upper <- hrline$upper * scale
Combine lines in a single object for plotting
Plot without the 95% credible intervals
ggplot(data.frame(mexdolphin$points)) +
geom_histogram(aes(x = distance), breaks = seq(0, 8, length.out = 9), alpha = 0.3) +
geom_point(aes(x = distance), y = 0.2, shape = "|", size = 3) +
geom_line(data = dlines, aes(x = distance, y = En, group = model, col = model))
Plot with the 95% credible intervals (without taking the count rescaling into account)
ggplot(data.frame(mexdolphin$points)) +
geom_histogram(aes(x = distance), breaks = seq(0, 8, length.out = 9), alpha = 0.3) +
geom_point(aes(x = distance), y = 0.2, shape = "|", size = 3) +
geom_line(data = dlines, aes(x = distance, y = En, group = model, col = model)) +
geom_ribbon(
data = dlines, aes(x = distance, ymin = En.lower, ymax = En.upper, group = model, col = model, fill = model),
alpha = 0.2, lty = 2
)