Comparing offsets to offsets calculated by BAM

[see24]

First of all it would be great if the README included the calculation of the offset that should actually be used in modelling as this last step is missing and although I think it is simply log(Apq) based on the Solymos paper it would be nice if this was laid out in the language of the napops package.

When I tried to calculate the offsets using this method I got very different offset values than the offsets I have from BAM calculated using the QPAD package. The code to do this is here https://github.com/LandSciTech/ROFBirds/blob/ef289a176b8aab30c019ea1ec9af9968e8d652f0/analysis/scripts/prep_BBS_data.R
(Apologies that this won't be entirely reproducible but you should be able to do most of it if you skip the bits with local data, and at least if gives you an idea what I am doing.)

And this pdf shows a graph comparing mean offset for each species.
compare_offsets_BAM_NAPOP.pdf

I actually found after messing around a bit that if I change the radius used to calculate area from 400 m to 4 I get much more comparable offsets but that is weird units so I am wondering if you have any thoughts?

[see24]

It looks like the distance in m is divided by 100 in the qpad-offsets repository functions https://github.com/borealbirds/qpad-offsets/blob/b5660ab80f1266de2c44d8ecb7d53df91ce210d9/functions.R#L108

So maybe that is needed here too?

[see24]

Here is a more reproducible work flow that compares napops with wildRtrax which implements QPAD offsets. They definitely don't match but I am not sure where or what the problem might be.

library(tidyverse)
library(wildRtrax)
library(napops)

wt_auth()

projects <- wt_get_download_summary("PC")

dat_pc <- projects %>%
  # pick a project in boreal and with few tasks
  filter(project == "CWS-Ontario Atlas Digital Point Counts Boreal FMUs 2022") %>%
  pull(project_id) %>%
  wt_download_report(sensor_id = "PC", reports = "main", weather_cols = FALSE)

dat_pc <- dat_pc %>% wt_tidy_species(sensor = "PC", zerofill = TRUE)

dat_pc_clean <- dat_pc %>%
  wt_make_wide(sensor = "PC")

# QPAD offsets
dat_pc_off <- wt_qpad_offsets(dat_pc_clean, sensor = "PC")

# napops offsets
# fetch_data()

# Need to get covariates

# Use NALCMS data same as napops paper. I am using 2020 data could in theory use
# version closest to actual year of data collection
url_nalcms_2020_Can <- "http://www.cec.org/files/atlas_layers/1_terrestrial_ecosystems/1_01_0_land_cover_2020_30m/can_land_cover_2020_30m_tif.zip"
nalcms_2020_Can <- reproducible::prepInputs(url = url_nalcms_2020_Can,
                                            destinationPath = tempDir,
                                            writeTo = NULL)

lcc_lu <- tibble::tribble(
  ~from, ~to, ~becomes,
  1,       2,        1,
  5,       6,        1)

# get boundary of bird data + ~2000m so window doesn't include edges
SA <- dat_pc %>% ungroup() %>%
  summarise(across(c(longitude, latitude), lst(min = \(x) min(x, na.rm = TRUE),
                                               max = \(x) max(x, na.rm = TRUE)))) %>%
  select(xmin = longitude_min, ymin = latitude_min, xmax = longitude_max,
         ymax = latitude_max) %>% unlist() %>%
  sf::st_bbox() %>%
  {. + 0.02} %>%
  sf::st_as_sfc() %>%
  sf::st_set_crs(4326) %>%
  sf::st_transform(sf::st_crs(nalcms_2020_Can))

# Crop and convert to prop forest in 5x5 window
for_cov <- reproducible::Cache(nalcms_2020_Can %>%
                   terra::crop(terra::vect(SA)) %>%
                   terra::classify(rcl = lcc_lu, others = 0, right = NA))
for_cov150 <- reproducible::Cache(terra::focal(for_cov, w = 5, fun = "mean"))

obs_for_cov <- dat_pc %>% filter(!is.na(latitude)) %>%
  select(project_id, location_id, survey_id, species_code, latitude, longitude) %>%
  sf::st_as_sf(coords = c("longitude", "latitude"),
               crs = 4326) %>%
  sf::st_transform(sf::st_crs(nalcms_2020_Can)) %>%
  {terra::extract(for_cov150, ., bind = TRUE)} %>%
  sf::st_as_sf() %>%
  sf::st_drop_geometry() %>%
  rename(for_cov = focal_mean)

dat_pc_cov <- dat_pc %>%
  filter(!is.na(latitude)) %>%
  left_join(obs_for_cov, by = join_by(project_id, location_id, survey_id, species_code)) %>%
  # time zone says UTC but is really unknown
  mutate(timez = lutz::tz_lookup_coords(latitude, longitude, method = "accurate")) %>%
  mutate(survey_date = lubridate::force_tzs(survey_date, timez, tzone_out = "UTC")) %>%
  mutate(sunrise = suncalc::getSunlightTimes(
  data = data.frame(date = lubridate::as_date(survey_date), lon = longitude, lat = latitude),
  keep = "sunrise")$sunrise) %>%
  mutate(tssr = difftime(survey_date, sunrise, units = "hours") %>% as.numeric(),
         od = lubridate::yday(survey_date))

nap_out <- dat_pc_cov %>%
  select(species_code, od, tssr, for_cov, detection_distance,
         survey_duration_method, survey_date) %>%
  slice(2) %>%
  mutate(
    on_road = FALSE,
    avail_off = avail(species_code, "best", od = od, tssr = tssr,
                      time = 5, pairwise = TRUE)$p,
    precep_off = percept(species_code, "best", forest = for_cov, road = on_road,
                         distance = 400, pairwise = TRUE)$q,
    cr_off = cue_rate(species_code, "best", od = od, tssr = tssr,
                      pairwise = TRUE)$CR_est[,1],
    edr_off = edr(species_code, "best", forest = for_cov, road = on_road,
                  pairwise = TRUE)$EDR_est[,1],
    off = log(pi*400^2)+log(avail_off)+log(precep_off))


wt_qpad_offsets(dat_pc_clean %>% slice(2), species = "REVI", sensor = "PC")

# when distance is INF napops gives 0 for q

# wt_qpad_offsets forces dist to be INF co compare use debug and after make_x save x to global env with <<-
# try with maxdis = 400
debug(wt_qpad_offsets)

# fix duration that was being cut off and set distance
qpad_off_x$MAXDUR <- 300
qpad_off_x$MAXDIS <- 400

wildRtrax:::.make_off("REVI", qpad_off_x)

nap_out %>% select(contains("off"))
## # A tibble: 1 × 5
##  avail_off precep_off cr_off edr_off   off
##      <dbl>      <dbl>  <dbl>   <dbl> <dbl>
## 1     0.847     0.0360  0.376    75.9  9.64
qpad_off
##  p            q        A correction    offset       tau       phi
## 1 1 4.414925e-06 502654.8   2.219183 0.7971392 0.8404689 0.4069309