Feature: distance sampling for Poisson

[dmi3kno]

I am trying to implement the model, where I have the aggregated counts on a grid. Some counts are low (or even missing) because they are undersampled. I wanted to implement a distance sampling algorithm, but then I discovered that samplers argument is only defined for family="cp".

Literature:

• Royle et al (2004) "Modeling abundance effects in distance sampling"
• Eberhardt (1967). "Some Developments in 'Distance Sampling'"
• Miller et al (2013) "Spatial models for distance sampling data: recent developments and future directions", section "Count as response"

How do I adjust for "distance to road" if I have total counts per pixel?

Here's an example of implementation in INLA

[finnlindgren]

If you precompute the "distance to road" at the grid locations and add them as distance in the data data.frame, you can use them in the predictor expression in the same way as in the transect distance sampling case.

formula <- observed_counts ~ some + component + combination + log_detection_probability(dist, some, parameters)

The cp model uses domain and samplers to call fm_int(domain = domain, samplers = samplers), but in your case you already have all the required information, so there's no need to create any extra integration points.

[dmi3kno]

Sorry for the thesis panic. It will pass soon. Thanks for your help!

[dmi3kno]

library(inlabru)
library(terra)
library(fmesher)
library(INLA)
library(ggplot2)

mexdolphin <- mexdolphin_sf
mesh <- mexdolphin$mesh
mxcrs <- sf::st_crs(mexdolphin$points)$proj4string

I create rasters for dolphin encounters and the distance to transects

# aggregate the data
flat_rast <- rast(terra::ext(mexdolphin$depth), crs = crs(mxcrs), resolution=20)
mexdolphin_rast <- rasterize(vect(mexdolphin$points), flat_rast, fun = sum, background = 0) |>
  mask(vect(mexdolphin$ppoly))

plot(mexdolphin_rast)

# raster for distance to transects
dist_rast <- rasterize(vect(mexdolphin$samplers), flat_rast, fun=sum) |> terra::distance()
dist_rast |> plot()

Here's the aggregated count raster on the background of the original data

# plot the new data
ggplot() +
  gg(mexdolphin_rast) +
  gg(mexdolphin$ppoly, fill = NA) +
  gg(mexdolphin$samplers, color = "grey") +
  gg(mexdolphin$points, size = 0.2, alpha = 1, color = "firebrick") +
  theme(legend.key.width = unit(x = 0.2, "cm"), legend.key.height = unit(x = 0.3, "cm")) +
  theme(legend.text = element_text(size = 6))

I extract the coordinates for the cells and cast them to sf format

# coordinates data
mexdolphin_pxdf <- crds(mexdolphin_rast, df = TRUE, na.rm = TRUE) |>
  cbind(values(mexdolphin_rast, na.rm = TRUE)) |>
  dplyr::rename(count = sum) |>
  sf::st_as_sf(coords = c("x", "y"), crs = crs(mxcrs))

Here's a new (regularly spaced) mesh

# make new (regular) mesh
mexdolphin_pxmesh <- INLA::inla.mesh.2d(
  loc = sf::as_Spatial(mexdolphin_pxdf),
  boundary = sf::as_Spatial(mexdolphin$ppoly) |> inla.sp2segment(),
  max.edge = c(100, 200), crs = crs(mxcrs)
)

# make new prior
pxmatern <- inla.spde2.pcmatern(mexdolphin_pxmesh,
  prior.sigma = c(2, 0.01),
  prior.range = c(50, 0.01)
)

ggplot() + gg(mexdolphin_pxmesh)

The components and the formula

# distance sampling
hn <- function(distance, sigma) {
  exp(-0.5 * (distance / sigma)^2)
}

cmp <- count ~ mySPDE(main = geometry, model = pxmatern) +
  # to transform the prior of sigma to be exponentially distributed
  sigma(1, prec.linear = 1,  marginal = bru_mapper_marginal(qexp, pexp, dexp, rate = 1 / 100)) +
  # distance to transect cells
  distroad(dist_rast$layer, model = "linear") +
  Intercept(1)
frm <- count ~ mySPDE + log(hn(distroad, sigma)) + Intercept + log(2)

I want to use these in the poisson model. Perhaps my grid is too coarse, but the model would not fit.

# commented away because it would not run
# fit <- bru(
#   cmp,
#   like(
#     data = mexdolphin_pxdf, family = "poisson",
#     formula = frm,
#     options = list(safe=TRUE, inla.mode="experimental" )
#   )
# )

^{Created on 2023-09-14 with reprex v2.0.2}

[dmi3kno]

I think my precomputed distances are not fine enough to make a difference. it is mostly zero in the encounter locations and of course, there's no field to cast over plain zeros. I basically need some detection to occur outside of the transect pixels. I might need to simulate some more dolphins ;)

I wonder if things could be better if I move distance into the binomial portion of ZAP.

[finnlindgren]

Yes, the mexdolphins data has a much larger domain of interest than the width of the transects, which is why we implement the model by integrating along the lines, and treating distance as a separate dimension. The approach you're looking at above is only useful when the transects/plots are clearly wider/larger than the model resolution.

[finnlindgren]

But writing my last comment I realised there's a potential problem with the inlabru line transect implementation, when the point data is the original actual locations; to make the likelihood approximation stable, we should be using the "nearest point on the transect" instead. I haven't noticed any issues with this though, likely due to only applying it do data with narrow transect widths.

[dmi3kno]

we implement the model by integrating along the lines, and treating distance as a separate dimension

The original LGCP model

form <- geometry + distance ~ mySPDE + log(hn(distance, sigma)) + Intercept + log(2) 
fit <- lgcp(components = cmp, data =  mexdolphin$points,
  samplers = mexdolphin$samplers,
  domain = list( geometry = mesh,  distance = fm_mesh_1d(seq(0, 8, length.out = 30)) ),
  formula = form
  )

Do you mean that samplers is specified separate from distance? In my example, I dont have the lines to integrate along. It is all distance, unfortunately.

when the point data is the original actual locations; to make the likelihood approximation stable, we should be using the "nearest point on the transect" instead.

Do you mean that the distance is calculated to the actual observation points and not as the shortest (perpendicular) distance to the nearest transect line?

I have been thinking about the "observational polygons", such as national parks and reserves (or even urban locations in some cases). The animals within the boundaries are most likely to be detected by the visitors and those further out in the wilderness are less likely to be spotted. So here you have distance sampling with decaying probability away from a polygon (or perhaps lines+polygons taking into account the roads to/from parks).

[finnlindgren]

1. No, I meant that the point process domain for (narrow) line transects is implemented as "location along lines" x "distance from line".
2. Neither. I'm saying that it doesn't force the \sum_i \log(lambda(y_i)) part of the likelihood to use y_i = the nearest point on the line, but instead uses y_i = actual location in space. The potential problem with that is quite technical.
3. I'd avoid calling that distance sampling, as it doesn't necessarily involve the same kinds of mechanisms. The (more likely to detect inside the park) is rather one type of a "preferential sampling" situation. Or simply "varying detectability" in general.

[dmi3kno]

Another issue that I have is that I want a bounded detection function. My distance is on [0,1] with 1 meaning "spot on" and 0 means "furthest away". I guess i could just logit transform to unbounded space, but I am not sure how the math will work out and whether I will run into the non-log-linear transformation issue ( as the case is with hazard detection function). I will report here shortly.

[finnlindgren]

I think this is a duplicate of the comment I replied to above.

[dmi3kno]

Sorry, moved the comment and edited it for something else :)

[finnlindgren]

If you make the model "too non-parameteric" you will run into the same issue as in the zip/zap vignette draft of a non-identifiable model. Bounded variables that can attain their bounds are tricky when modelling them via transformed Gaussians.

[finnlindgren]

The hazard rate function works well, since it's parametrically bound to start at 1 at distance zero, and to remain between 0 and 1 elsewhere.

[finnlindgren]

It's probably essential to have a model that's "anchored" like that in some way.

[dmi3kno]

I am still not certain what is happening with the distance sampling when distance is provided via the raster

  library(tidyverse)
  library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
  library(inlabru)
  library(terra)
#> terra 1.7.46
  
# detection function
hn <- function(distance, sigma) {
  exp(-0.5 * (distance / sigma)^2)
}

# detection function
hzrt <- function(distance, sigma, bt=2) {
  1 - exp(-(distance / sigma)^(-bt))
}

counts_df <- read_csv("obs_df.csv") |> 
  st_as_sf(coords=c("x", "y"), crs=4326) 
#> Rows: 163 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> dbl (4): count, present, x, y
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

This is a very simple dataset with the count and coordinates. The present is set to 1 for all non-zero counts.

dist2road_rstr <- rast("dist_road.tiff")
dist2road_rstr |> plot()

My "distance to road" raster

px_mesh <- fm_mesh_2d_inla(
  loc=st_as_sfc(counts_df),
  boundary = st_convex_hull(st_union(counts_df)),
  max.edge = c(2, 4),
  crs=st_crs(counts_df)
) 

px_matern <- INLA::inla.spde2.pcmatern(
  px_mesh,
  prior.sigma = c(0.1, 0.01),
  prior.range = c(0.1, 0.01)
)

ggplot() +
  geom_fm(data = px_mesh) +
  # geom_sf(data=counts_df, aes(color=count), size=1, pch=4)+
  theme_minimal()

Here I am plotting all the counts including those I suppressed (made into zero or NA)

comps <- ~ dist(dist2road_rstr$max, model="linear") +
  sigma(1, prec.linear = 1, marginal=bru_mapper_marginal(qexp, pexp, dexp, rate=1/1e3)) +
  field(geometry, model = px_matern) + Intercept(1)

fit_poi <- bru(
  comps,
  like(
    family = "poisson", data = counts_df,
    formula = count ~ field + Intercept + 
      log(hn(dist, sigma))+ log(2)
  )
)
summary(fit_poi)
#> inlabru version: 2.9.0.9011
#> INLA version: 23.10.19-1
#> Components:
#> dist: main = linear(dist2road_rstr$max), group = exchangeable(1L), replicate = iid(1L)
#> sigma: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> field: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L)
#> Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'poisson'
#>     Data class: 'sf', 'tbl_df', 'tbl', 'data.frame'
#>     Predictor: count ~ field + Intercept + log(hn(dist, sigma)) + log(2)
#> Time used:
#>     Pre = 0.705, Running = 0.956, Post = 0.497, Total = 2.16 
#> Fixed effects:
#>           mean   sd 0.025quant 0.5quant 0.975quant mode kld
#> dist      0.00 0.00      0.000     0.00      0.001 0.00   0
#> sigma     0.00 1.00     -1.961     0.00      1.961 0.00   0
#> Intercept 0.93 0.04      0.851     0.93      1.009 0.93   0
#> 
#> Random effects:
#>   Name     Model
#>     field SPDE2 model
#> 
#> Model hyperparameters:
#>                  mean    sd 0.025quant 0.5quant 0.975quant  mode
#> Range for field 0.962 1.445      0.083     0.54       4.50 0.206
#> Stdev for field 0.031 0.038      0.003     0.02       0.13 0.007
#> 
#> Deviance Information Criterion (DIC) ...............: 725.78
#> Deviance Information Criterion (DIC, saturated) ....: 179.26
#> Effective number of parameters .....................: 2.28
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: 727.61
#> Effective number of parameters .................: 3.71
#> 
#> Marginal log-Likelihood:  -380.03 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

I am concerned by the zero coefficient by the distance and the 0,1 posterior on the scale (sigma). Something is wrong and it is not getting any better under a different model.

comps_zap <- ~ dist(dist2road_rstr$max, model="linear") +
  sigma(1, prec.linear = 1, marginal=bru_mapper_marginal(qexp, pexp, dexp, rate=1/1e3)) +
  field_count(geometry, model = px_matern) + Intercept_count(1)+
  field_present(geometry, model = px_matern) + Intercept_present(1)

fit_zap <- bru(
  comps_zap,
  like(
    family = "binomial", data = counts_df,
    formula = present ~ field_present + Intercept_present + 
      log(hn(dist, sigma))+ log(2)
  ),
  like(
    family="nzpoisson", data = counts_df[counts_df$present>0,],
    formula = count ~ field_count + Intercept_count
  )
)
summary(fit_zap)
#> inlabru version: 2.9.0.9011
#> INLA version: 23.10.19-1
#> Components:
#> dist: main = linear(dist2road_rstr$max), group = exchangeable(1L), replicate = iid(1L)
#> sigma: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> field_count: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L)
#> Intercept_count: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> field_present: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L)
#> Intercept_present: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'binomial'
#>     Data class: 'sf', 'tbl_df', 'tbl', 'data.frame'
#>     Predictor: 
#>         present ~ field_present + Intercept_present + log(hn(dist, sigma)) + 
#>             log(2)
#>   Family: 'nzpoisson'
#>     Data class: 'sf', 'tbl_df', 'tbl', 'data.frame'
#>     Predictor: count ~ field_count + Intercept_count
#> Time used:
#>     Pre = 0.981, Running = 1.28, Post = 0.282, Total = 2.54 
#> Fixed effects:
#>                    mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> dist              0.001 0.002     -0.003    0.001      0.004 0.001   0
#> sigma             0.000 1.000     -1.961    0.000      1.961 0.000   0
#> Intercept_present 4.674 1.138      2.443    4.674      6.905 4.674   0
#> Intercept_count   1.604 0.036      1.533    1.604      1.674 1.604   0
#> 
#> Random effects:
#>   Name     Model
#>     field_present SPDE2 model
#>    field_count SPDE2 model
#> 
#> Model hyperparameters:
#>                          mean    sd 0.025quant 0.5quant 0.975quant  mode
#> Range for field_present 1.103 1.784      0.099    0.596      5.342 0.237
#> Stdev for field_present 0.025 0.028      0.002    0.016      0.098 0.005
#> Range for field_count   0.967 1.457      0.084    0.542      4.529 0.208
#> Stdev for field_count   0.031 0.037      0.003    0.020      0.128 0.007
#> 
#> Deviance Information Criterion (DIC) ...............: 732.49
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: 5.32
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: 736.66
#> Effective number of parameters .................: 6.75
#> 
#> Marginal log-Likelihood:  -381.82 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

^{Created on 2023-11-06 with reprex v2.0.2}

The files are attached here

obs_df.csv
dist_road.tiff.zip

[dmi3kno]

I thought that covariates are evaluated for all mesh points, not just observations. So I would need to supply a dataframe that contains covariates for every node of the mesh.

[finnlindgren]

There are two issues here; 1. Whether distance should have a component (it shouldn't) and 2. Where covariates given as raster objects as inputs to components are evaluated.

For 2., raster component inputs are evaluated wherever there is a "data point". For ordinary observation model, those points are precisely the ones provided by the user as having observed values associated with them. For family="cp" (that lgcp() is a simple wrapper for) it's different. There, the "data points" are constructed by inlabru to include all locations it needs in order To evaluate/approximate the point process likelihood; currently that by default becomes the union of the observed point locations, and (a subset of) the mesh nodes (but more generally can include other numerical integration scheme locations).

Since your model treats the data as georeferenced count values, the family="cp" aspects are not involved.

[dmi3kno]

Understood, thank you. I guess evaluating covariates on the mesh nodes may be important if you are trying to predict on mesh-derived pixels (like in this Poisson example)

[dmi3kno]

One more question: I understand that if we remove distance from components, then there will be no coefficient in front the distance term in the likelihood, other than the sigma coefficient, which scales the detection function. Scale is, of course non-negative. Does that mean I have restricted my coefficient to positive values only?

If I want the distance to reduce (not increase) the probability of zero in a ZAP model, I can't just write, can I? The meaning of "minus" sign in formula is quite different, isnt it?

formula = present ~ field_present + Intercept_present - (  log(hn(dist, sigma)) +  log(2) )

Or perhaps I could "fix" the sigma to 1 and reintroduce the distance to components? Then the coefficient in front of distance will be equivalent to transformed sigma (wrong scale, yes, but it can then be negative). Tell me, if I am missing something.

[finnlindgren]

I think you need to spell out the model you want in terms of mathematical formulae. Without that, it's not really possible to tell if the code implements the intended model or not.