Question about sampling units in MHQ samples but assessed as absent

[florisvdh]

Related to #49 and #54. Using current master (a6d53d2). Needed for code to run is n2khab_strata.rds.zip; after unzipping you can run below code to load it. I've also included an update of the read function in order to capture the point_code column so that your join can be made.

Code to generate mhq_terr_samples_grts

n2khab_strata <- readRDS("n2khab_strata.rds")

#' Read and tidy csv file with MHQ samples
#'
#' @param path File path.
#' @param grts_var Column name to be used as GRTS address.
#' @param single_type Optional string to set a single type that represents all
#'   rows.
read_csv_mhq_samples <- function(path,
                                 grts_var = "grts_ranking_draw",
                                 single_type = NULL) {
  (
    if (is.null(single_type)) {
      read_delim(
        file = path,
        delim = ";",
        col_types = cols_only(
          {{ grts_var }} := col_integer(),
          habitattype = col_character(),
          point_code = col_character()
        )
      ) %>%
        select(
          stratum = habitattype,
          grts_address = {{ grts_var }},
          point_code
        )
    } else {
      read_delim(
        file = path,
        delim = ";",
        col_types = cols_only(
          {{ grts_var }} := col_integer(),
          point_code = col_character()
        )
      ) %>%
        mutate(stratum = single_type) %>%
        select(
          stratum,
          grts_address = {{ grts_var }},
          point_code
        )
    }
  ) %>%
    mutate(
      # shortcut a complication for forests, where > 1 type is sometimes noted
      stratum = str_extract(stratum, "^\\w+(\\+$)?"),
      stratum = parse_factor(stratum, levels = levels(n2khab_strata$stratum))
    ) %>%
    arrange(stratum, grts_address)
}
>
> # mhq_samples_textdatapath <- "/substitute/with/correct/filepath"
> mhq_samples_heath_haymeadow <-
    read_csv_mhq_samples(
      file.path(mhq_samples_textdatapath, "mhq_terr_cyclus2_anb.csv")
    )
> mhq_samples_mire_grassland <-                                                                                    
    read_csv_mhq_samples(
      file.path(mhq_samples_textdatapath, "mhq_terr_cyclus2_inbo.csv")
    )
> mhq_samples_dunes <-                                                                                             
    read_csv_mhq_samples(file.path(
      mhq_samples_textdatapath,
      "mhq_terr_cyclus2_duinen_inbo_2023-05-11.csv"
    ))
> mhq_samples_forest <-                                                                                            
    read_csv_mhq_samples(
      file.path(mhq_samples_textdatapath, "mhq_forests_cycle1.csv"),
      grts_var = "grts_ranking"
    )
> mhq_terr_samples_grts <-                                                                                         
    bind_rows(
      mhq_samples_heath_haymeadow,
      mhq_samples_mire_grassland,
      mhq_samples_dunes,
      mhq_samples_forest,
    ) %>% 
    # reset some peculiarities in read_csv_mhq_samples()
    rename(type = stratum, grts_ranking_draw = grts_address) %>% 
    arrange(type, grts_ranking_draw)

We can see that about 1/5 of the recent (2020-2024) 'type absence' assessments is a member of the MHQ samples (for the type x point_code combination). This corresponds to the second line of the result below (in_sample is TRUE).

Question: is this because the MHQ samples lag behind the type absence information? Or do I misunderstand something?

> mhq_terr_datapath <- file.path(dirname(gitroot), "n2khab-sample-admin/data/mhq_terr/rapportage2025")
> mhq_terr_assessments <-
    read_vc("mhq_terr_assessments", root = mhq_terr_datapath) %>%
    as_tibble()
>
> mhq_terr_assessments %>% 
    filter(
      !is.na(type),
      !is_present,
      year(assessment_date) >= year(today()) - 5
    ) %>% 
    select(point_code, type, is_present, assessment_date) %>% 
    left_join(
      mhq_terr_samples_grts %>% 
        mutate(in_sample = TRUE),
      join_by(point_code, type)
    ) %>% 
    mutate(in_sample = ifelse(is.na(in_sample), FALSE, in_sample)) %>% 
    count(in_sample)
# A tibble: 2 × 2
  in_sample     n
  <lgl>     <int>
1 FALSE       462
2 TRUE        110