Manual creation of integration points for distance sampling model

[Josh Cullen]

Hi,

I'm attempting to fit a distance sampling model in {inlabru} as a thinned inhomogeneous Poisson point process, but using points where environmental values have already been extracted for both the animal sightings and the transect segments. I'm trying to use this approach instead of the method of supplying the transects as {sf} LINESTRINGS and a raster stack since I would need to provide 30 years of daily rasters, which I expect would prove computationally challenging. So it seems much simpler to just extract all environmental values prior to model fitting.

While I've built and fitted a model where the marginal effects for the environmental covariates make sense, i'm not sure if the way I'm generating the integration points is correct because I wasn't able to find a similar example available from this Google Group, Github repos, or the {inlabru} website. Any feedback as to whether the approach I'm using is correct would be much appreciated. (Relevant sample code below).

Thanks,

     Josh

```

# Define half-Normal detection function
log_hn_ls <- function(distance, lsig) {
  # lsig is the linear predictor for log(sigma)
  -0.5 * (distance / exp(lsig))^2
}

# Reformat transect data and scale covars
surv_sf2 <- surv_sf |>
  st_transform(proj) |>
  mutate(across(.cols = c(mld, sss, sst, eke, bathym),
                .fns = ~{scale(.x) |>
                    c()},
                .names = "{.col}_s"))

# Scale covars for distance sampling model
dat_sights3 <- dat_sights2 |>
  mutate(bathym_s = (bathym - mean(surv_sf2$bathym)) / sd(surv_sf2$bathym),
         mld_s = (mld - mean(surv_sf2$mld)) / sd(surv_sf2$mld),
         sss_s = (sss - mean(surv_sf2$sss)) / sd(surv_sf2$sss),
         sst_s = (sst - mean(surv_sf2$sst)) / sd(surv_sf2$sst),
         eke_s = (eke - mean(surv_sf2$eke)) / sd(surv_sf2$eke))

# Define 1D distance mesh

W <- max(dat_sights3$distance)  #truncation distance

mesh_dist <- fm_mesh_1d(
  seq(0, W, length.out = 5),
  degree = 2
)

### Generate integration points on subdivided mesh ###
mesh2 <- mesh
mesh2$loc <- cbind(mesh2$loc, 0)  #needs 3rd column to work
mesh_ips <- fm_subdivide(mesh2, n = 1)

surv_ips <- fm_int(domain = mesh_ips,
                   samplers = surv_sf2 |>
                     dplyr::select(segnum, transect, date, dist, eke_s, mld_s, sss_s,
                                   sst_s, bathym_s),
                   int.args = list(method = "stable", nsub = NULL)
                                   ) |>
  rename(weight_spatial = weight)

# Add back in environ covars by associated row ID
surv_ips2 <- surv_ips |>
  left_join(surv_sf2 |>
              st_drop_geometry() |>
              mutate(row_id = row_number()),
            by = c(".block" = "row_id"))

### Create cross-product of IPs across the detection distance domain ###

# Create detection distance ips
dist_ips <- fmesher::fm_int(mesh_dist, name = "distance") |>
  rename(weight_dist = weight) |>
  dplyr::select(distance, weight_dist)

# Perform a "Cross Join" with spatial points (cross-product of all rows from surv_ips w/ dist_ips)
surv_ips_full <- surv_ips2 |>
  st_drop_geometry() |>
  crossing(dist_ips) |> # This multiplies rows: (92k spatial) * (9 distance)
  st_as_sf(geometry = rep(st_geometry(surv_ips2), each = nrow(dist_ips))) |>  #add spatial component back

  mutate(weight = weight_spatial * weight_dist)  #(Weight of Transect Segment) * (Weight in 'distance' dimension)

### Define model components
cmp6 <- ~
  ## Intensity model
  Int_lambda(1) +

  # Functional GP over SST (1-D SPDE on covariate value)
  gp_sst(
    main   = sst_s,                  # numeric covariate at IPs
    model  = spde1d[["sst_s"]],
    mapper = bru_mapper_fmesher(covmesh1d[["sst_s"]])
  ) +
  # Functional GP over bathymetry
  gp_bathym(
    main   = bathym_s,
    model  = spde1d[["bathym_s"]],
    mapper = bru_mapper_fmesher(covmesh1d[["bathym_s"]])
  ) +
  # Functional GP over MLD
  gp_mld(
    main   = mld_s,
    model  = spde1d[["mld_s"]],
    mapper = bru_mapper_fmesher(covmesh1d[["mld_s"]])
  ) +
  # Functional GP over SSS
  gp_sss(
    main   = sss_s,
    model  = spde1d[["sss_s"]],
    mapper = bru_mapper_fmesher(covmesh1d[["sss_s"]])
  ) +
  # Functional GP over EKE
  gp_eke(
    main   = eke_s,
    model  = spde1d[["eke_s"]],
    mapper = bru_mapper_fmesher(covmesh1d[["eke_s"]])
  ) +
 
  ## Distance sampling model
  lsig0(1, model = "linear", mean.linear = 0, prec.linear = 1) +
  GRF(main = geometry, model = matern)

# Define likelihood for distance sampling model
lik_ds <- bru_obs(
  family = "cp",
  data   = dat_sights3 |>
    st_as_sf(coords = c("lon","lat"), crs = 4326) |>
    st_transform(proj),
  ips = surv_ips_full,
  domain = list(distance = mesh_dist,
                geometry = mesh),
  formula = distance + geometry ~ Int_lambda + gp_sst + gp_bathym + gp_mld + gp_sss + gp_eke +
    GRF +
    log_hn_ls(distance,  lsig0) +

   log(2)
    )

# Fit model
fit5 <- bru(
  components = cmp6,
  lik_ds,
  options = list(
    control.inla  = list(int.strategy = "eb", strategy = "adaptive"),
    control.compute = list(openmp.strategy = "huge"),
    inla.mode = "experimental",
    verbose = TRUE,
    bru_verbose = 4
  )
)

```