Implementation of point process models in gllvm

[wlangera]

I performed an analysis with presence-only data on the 2025 summer school.
It is an analysis based on citizen science from Flanders (BE) from waarnemingen.be.
@BertvanderVeen stated that this is close to a point process model that could be implemented in the gllvm package by weighting the likelihood (e.g. Renner et al., 2015).

I will summarise my analysis here, but you can also find my code on this repository.

step 1: data preparation

1. Temporal aggregation
   • We select the data from 2015-2024 (10 year span)
2. Validated data
   • We select data that has been validated (evidence, AI, expert knowledge ...)
3. Precision of coordinates
   • We select data with coordinate precision < 100 m
4. Quality filters (Van Eupen et al., 2021)
   • DETAIL: Filters records that include additional information beyond the default fields (e.g., behaviour, photos, comments). Used as a proxy for observer effort in unstructured data, where more detailed records are considered higher quality.
     • We keep all data “Validated based on evidence” (= picture). We keep the the rest of the data (e.g. validated by AI or expert knowledge) if they have more than 1 additional details
5. Taxonomic considerations
   • We aggregate subspecies to species level and only keep species level observations

Since we have a lot of observations across the whole of Flanders, we formulate a specific research question.
We look at the observations in heather habitat and want to know which species are more observed in wet versus dry heather habitat.

step 2: data selection

• We create a map of heather habitat with three categories
  1. dry heather
  2. mixed heather
  3. wet heather
• We only keep polygons of at least 10,000 m^2

• We select observations within these polygons
• Count the number of observations (ignoring individual count per observation) per polygon and fill NAs with (pseudo-)zeroes
  • If a species was not seen in a polygon, it is considered a pseudo-absence
• We only keep species that occur in at least 10 polygons

step 3: calculate sampling effort

We account for two biases in sampling effort:

1. Polygon size: some polygons are larger than others
   • Calculate polygon size
2. Number of visits: certain polygons could be more popular than others
   • Count total number of observations over all (original?) datasets per polygon

step 4: model specification

Quick data exploration shows we have zero inflation.
Also from a theoretical perspecive, it would be logical to fit a zero inflated model:

• in a zero inflated distribution, e.g. the zero inflated Poisson, zeroes are split in true zeroes (belonging to the Poisson distribution for example) and other zeroes
• in our case the pseudo-absences are also expected to be true absences (because the species are not there) and untrue absences (because they are not seen/recorded by citizen scientist)

We fit a zero-inflated negative binomial distribution where we use an offset for sampling effort.

Let $y_{ij}$ be the count of species $j$ at site or observation unit $i$.
Assume:

$$ y_{ij} \sim \text{NegBin}(\mu_{ij}, \theta_{j}) $$

with

$$ \log(\mu_{ij}) = \log(\text{E}(y_{ij})) = \beta_{0j} + \beta_{\text{mixed},j} + \beta_{\text{wet},j} + \log(X_{\text{area},i}) + \log(X_{\text{effort},i}) $$

or equivalently, by offsetting:

$$ \log\left( \frac{\text{E}(y_{ij})}{X_{\text{area},i} \cdot X_{\text{effort},i}} \right) = \beta_{0j} + \beta_{\text{mixed},j} + \beta_{\text{wet},j} $$

where:

• $\mu_{ij} = \mathbb{E}(y_{ij})$ is the expected count,
• $\theta_j$ is the dispersion parameter for species $j$,
• $X_{\text{area},i}$ and $X_{\text{effort},i}$ are offset terms (e.g. surveyed area and effort),
• All $\beta$'s are species-specific coefficients.

plot_data_wide$offset_var <- log(
  plot_data_wide$habAreaM2 * plot_data_wide$samplingEffort
)

We use the offset argument of the gllvm() function.

  model <- gllvm(
    y = butterflies_data_wide[, -1],
    X = plot_data_wide,
    formula = ~habClass,
    offset = plot_data_wide$offset_var,
    num.lv = 2,
    family = "ZINB"
  )

However, from the coefplot, we notice we want shrinkage of some effects.
We therefore fit the model with random intercepts instead of fixed before.

model2 <- gllvm(
  y = butterflies_data_wide[, -1],
  X = plot_data_wide,
  formula = ~(habClass | 1),
  offset = plot_data_wide$offset_var,
  num.lv = 2,
  family = "ZINB"
)

We visualise the parameters of interest.
We see that some species occur significantly more in wet and mixed heather than in dry heather.

estimates <- model2$params$Br[-1, ] %>%
  as_tibble(rownames = "param") %>%
  pivot_longer(-param, names_to = "species", values_to = "estimate")

ses_model2 <- getPredictErr(model2)$Br[-1, ] %>%
  as_tibble(rownames = "param") %>%
  pivot_longer(-param, names_to = "species", values_to = "se")

full_join(estimates, ses_model2, by = join_by(param, species)) %>%
  mutate(order = mean(estimate), .by = "species") %>%
  mutate(
    lcl = estimate + qnorm(0.025) * se,
    ucl = estimate + qnorm(0.975) * se,
    significant = !(0 > lcl & 0 < ucl),
    species = reorder(species, order)
  ) %>%
  ggplot(aes(x = estimate, y = species, colour = param, alpha =  significant)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "darkgrey") +
  geom_point(position = position_dodge(width = 0.8)) +
  geom_errorbar(aes(xmin = lcl, xmax = ucl),
                position = position_dodge(width = 0.8)) +
  labs(x = "compared to dry heather", y = "", colour = "Legend:") +
  theme_minimal() +
  guides(alpha = "none") +
  scale_alpha_manual(values = c(0.2, 1)) +
  theme(legend.position = "inside",
        legend.position.inside = c(0.8, 0.15),
        legend.background = element_rect(fill = "white"),
        axis.text.y = element_text(size = 6))

[wlangera]

The variance partitioning plot shows a strange result however.

As discussed, this is probably due to the pseudo absences.
The pseudo absences were created per polygon and therefore we expect more pseudo absences in smaller polygons than in larger polygons.

[BertvanderVeen]

Suspect similar problem as outlined in #237.