Appendix 5. Model fitting for non central place foraging (NCPF) species

Modelling habitat selection using tracking data from central place foraging species: A practical guide for quantitative ecologists

Modelling Approaches

We compared three modelling approaches that differ in how they represent spatial intensity, account for uncertainty, and handle autocorrelation in telemetry data:

• Inhomogeneous Poisson Point Process (IPP), fitted using a down‑weighted Poisson regression approach (Renner et al. 2015; Matthiopoulos, Fieberg, and Aarts 2023).
  This model assumes independence among observations and estimates a log‑linear intensity surface as a function of spatial covariates.

• Log‑Gaussian Cox Process (LGCP), which extends the IPP by modelling the log‑intensity as a Gaussian random field (Dovers, Stoklosa, and Warton 2024).
  This formulation captures spatially-structured overdispersion and accounts for unobserved environmental variation.

• Space–Time Point Process model (STPP), which incorporates both spatial covariates and temporal dependence between observations (Johnson, Hooten, and Kuhn 2013; Hooten et al. 2017).
  By explicitly accounting for serial autocorrelation in telemetry data, the STPP provides a more realistic representation of movement-driven clustering. Two formulations were considered: a full STPP and a marginalised STPP.

These approaches span a progression from models that treat locations as independent, to models that incorporate spatial structure, through to models that explicitly account for movement processes. Comparing them allows us to assess how different assumptions about dependence and availability influence habitat-selection inference.

Below, we describe each approach in detail, corresponding to the models implemented in the analysis. We assessed each model performance by comparing estimated coefficients with the true values used in the simulation, and by evaluating how well each model reproduced the underlying intensity surface through visual comparison of predicted spatial patterns.

Settings and data processing

Settings

# ------------------------------------------------------------
# Load required packages
# ------------------------------------------------------------
library(dplyr)
library(sf)
library(ggplot2)
library(mgcv)
library(raster)
library(spatstat)
library(gsl)
library(mvtnorm)
library(tidyr)

Load data

The data corresponds to the data of 1000 individuals, each with one single movement track lasting 360 minutes, and reconstructed step‑selection sequences for all individuals.

# ------------------------------------------------------------
# Reproducibility settings
# ------------------------------------------------------------
seed <- 123
set.seed(seed)

# ------------------------------------------------------------
# Simulation scenario selector
# ------------------------------------------------------------
cov <- 1

# ------------------------------------------------------------
# Load simulation data and habitat raster
# ------------------------------------------------------------
load(paste0("NCPF/Simulations/NCPF_cov", cov, "_1000ind.RData"))
load("cov_raster.RData")

Observed relative selection

To obtain an empirical measure of habitat preference comparable across scenarios, we calculated the observed relative selection following the approach of Lauret et al. (2025). Specifically, we computed the log ratio between the number of steps ending in suitable habitat and the number of steps ending in unsuitable habitat.

\[ \log\left(\frac{\#\text{ steps ending in suitable habitat}}{\#\text{ steps ending in unsuitable habitat}}\right) \] Because the simulation design ensured that suitable and unsuitable habitats were equally available, this log‑ratio provides a direct estimate of the relative intensity of space use in the two habitat categories.

# ------------------------------------------------------------
# Extract habitat category for each observed location
# ------------------------------------------------------------

data_sf$habitat <- raster::extract(cov_raster, st_coordinates(data_sf))

# ------------------------------------------------------------
# Remove rows with missing bearings
# (these correspond to the first step of each track)
# ------------------------------------------------------------
data_sf <- data_sf[!is.na(data_sf$bearing), ]

# ------------------------------------------------------------
# Laurent et al. (2025) relative selection index
# log(# steps ending in suitable habitat / # steps ending in unsuitable habitat)
# ------------------------------------------------------------
p.est_laurent <- log(
  table(data_sf$habitat)[2] /
  table(data_sf$habitat)[1]
)

Sub-sampling individuals

To simulate data typically available in GPS‑tracking studies, where only a small number of individuals are monitored relative to the true population, we randomly sub-sampled 50 individuals from the pool of 1000 simulated individuals for each parameter combination. All model‑fitting procedures described below were conducted using tracks from this subset. This sub-sampling ensures that the evaluation of model performance reflects the level of information generally available in empirical telemetry studies rather than the unrealistically large amount of data used to generate the simulated population-level movement process.

# ------------------------------------------------------------
# Select a subset of individuals to mimic realistic tracking
# ------------------------------------------------------------
n <- 50
ind <- sample(unique(data_sf$ID), n, replace = FALSE)

tracked_birds <- data_sf[data_sf$ID %in% ind, ]

# Extract coordinates into separate columns
tracked_birds$x <- st_coordinates(tracked_birds)[, 1]
tracked_birds$y <- st_coordinates(tracked_birds)[, 2]

# Drop geometry to simplify later processing
tracked_birds <- st_drop_geometry(tracked_birds)

# ------------------------------------------------------------
# Convert habitat raster to a data frame for plotting
# ------------------------------------------------------------
df <- as.data.frame(cov_raster, xy = TRUE)

# ------------------------------------------------------------
# Plot habitat map with movement paths
# ------------------------------------------------------------
ncpf_plot <-ggplot() +
  geom_tile(
    data = df,
    aes(x, y, fill = factor(value)),
    alpha = 0.1
  ) +
  geom_path(
    data = tracked_birds,
    aes(x, y, col = factor(state), group = ID),
    linewidth = 0.4
  ) +
  scale_fill_manual(
    values = c("lightgrey", "darkgreen"),
    name   = "Habitat",
    labels = c("Unsuitable", "Suitable")
  ) +
  scale_colour_manual(
    values = c("black", "darkorange", "darkred"),
    name   = "Behaviours",
    labels = c("Outbound", "Search", "Foraging")
  ) +
  coord_fixed() +
  theme_bw()+ 
  labs(x = "", y = "") +
  ggplot2::theme(
    axis.text = ggplot2::element_text(size = 15, color = "black"),
    axis.text.y = ggplot2::element_text(
      margin = ggplot2::margin(t = 0, r = 5, b = 0, l = 0)
    ),
    axis.text.x = ggplot2::element_text(
      margin = ggplot2::margin(t = 5, r = 0, b = 0, l = 0)
    ),
    axis.title = ggplot2::element_text(size = 12),
    axis.title.y = ggplot2::element_text(
      margin = ggplot2::margin(t = 0, r = 20, b = 0, l = 0),
      angle = 90
    ),
    axis.title.x = ggplot2::element_text(
      margin = ggplot2::margin(t = 20, r = 0, b = 0, l = 0)
    )
  ) +
  theme(
    legend.text = element_text(size = 15),
    legend.spacing.y = unit(1.5, "cm"),
    legend.title = element_text(
      size = 18,
      face = 'bold',
      margin = margin(b = 15)
    )
  )
 
ncpf_plot

Prediction grid

A regular grid covering the study area is created to allow spatial prediction of intensity surfaces. For each grid cell centroid, we extract:

Habitat type

Cell area

Number of points per cell

This grid is later used for model predictions and visualisation.

# ------------------------------------------------------------
# Define study area polygon (50 km × 50 km square)
# ------------------------------------------------------------
pol <- st_polygon(list(cbind(
  c(0, 50000, 50000, 0, 0),
  c(50000, 50000, 0, 0, 50000)
)))

pol <- st_sfc(pol, crs = 32630)

# ------------------------------------------------------------
# Create regular grid for quadrature / prediction
# ------------------------------------------------------------
ncells <- 100   # number of grid cells per dimension

grid <- st_make_grid(pol, n = ncells)

# ------------------------------------------------------------
# Convert grid to data frame of centroid coordinates
# ------------------------------------------------------------
grid_df <- grid %>%
  st_centroid() %>%
  st_coordinates() %>%
  as.data.frame()

names(grid_df) <- c("x", "y")

# Area of each grid cell (constant)
grid_df$area <- as.numeric(st_area(grid))
grid_df$larea <- log(grid_df$area)

# ------------------------------------------------------------
# Extract covariates for each grid cell
# ------------------------------------------------------------

# Habitat category from raster
grid_df$habitat <- raster::extract(cov_raster, grid_df[, c("x", "y")])

# ------------------------------------------------------------
# Count points per grid cell
# ------------------------------------------------------------
grid_sf <- st_sf(geometry = grid)
grid_sf$cell_id <- seq_len(nrow(grid_sf))

coords <- st_as_sf(tracked_birds, coords = c("x", "y"), crs = 32630)

join_result <- st_join(coords, grid_sf, join = st_within)

counts_per_cell <- join_result %>%
  st_drop_geometry() %>%
  count(cell_id, name = "count")

grid_sf <- grid_sf %>%
  left_join(counts_per_cell, by = "cell_id") %>%
  mutate(count = replace_na(count, 0))

grid_df$count <- grid_sf$count

Inhomogeneous Poisson Point Process (IPP)

The inhomogeneous Poisson point process (IPP) provides a flexible framework for modelling spatial variation in animal space use. Under this model, observed locations are assumed to occur independently according to an intensity function that varies across the landscape.

Fitting the IPP requires evaluating the log-likelihood, which includes summing over observed locations and integrating the intensity over the study region. This integral is usually not tractable when the intensity depends on spatially varying covariates.

This limitation motivates the use of practical approximation strategies, which allow the IPP model to be fitted using standard regression tools. Two approaches are common (Aarts, Fieberg, and Matthiopoulos 2012). The first is a grid‑based approximation, in which the study area is partitioned into small cells and the number of observed points in each cell is recorded. Each cell count is treated as a Poisson response with a log‑linear mean that incorporates covariates evaluated at the cell centroid and an offset for the area of the cell. This approach converts the IPP into a familiar Poisson regression model while preserving most of the information contained in the original point pattern.

gridded_model<-glm(count ~ habitat  + offset(larea),
  family = poisson(),
  data   = grid_df)

summary(gridded_model)

## 
## Call:
## glm(formula = count ~ habitat + offset(larea), family = poisson(), 
##     data = grid_df)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -12.39327    0.01389  -892.2   <2e-16 ***
## habitat       0.90539    0.01646    55.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 39738  on 9999  degrees of freedom
## Residual deviance: 36396  on 9998  degrees of freedom
## AIC: 50359
## 
## Number of Fisher Scoring iterations: 6

co      <- coef(summary(gridded_model))
beta_hat <- co["habitat", "Estimate"]
se_hat   <- co["habitat", "Std. Error"]
ci_low   <- beta_hat - 1.96 * se_hat
ci_up    <- beta_hat + 1.96 * se_hat

gridded_model_est <- data.frame(
  model    = "Poisson Regression",
  estimate = beta_hat,
  low      = ci_low,
  up       = ci_up)

# ------------------------------------------------------------
# Predict intensity surface
# ------------------------------------------------------------
grid_df$predict_gridded_model <- exp(predict(gridded_model, newdata = grid_df))

ggplot(grid_df, aes(x, y, fill = predict_gridded_model)) +
  geom_tile() +
  scale_fill_viridis_c(name = "Predictions") +
  coord_equal() +
  theme_minimal()

The second approach is numerical quadrature, in which the spatial integral in the likelihood is approximated using a set of quadrature points and associated weights. This approximation forms the foundation of the down‑weighted Poisson regression method, where telemetry locations contribute only to the log‑intensity component of the model and quadrature points contribute to the integral component. By assigning extremely small weights to observed locations and appropriate area‑based weights to quadrature points, the resulting Poisson regression yields parameter estimates that are equivalent to those obtained under the full IPP likelihood. Let

\[ S_n = \{s_i\}_{i=1}^n \]

denote the set of observed locations within a region \(A\). The IPP assumes that points arise independently with an intensity function \(\lambda(s)\) that varies across space. In our case, the intensity is modelled as:

\[ \lambda(s_i) = \exp\left(\beta_0 + \beta_1 \text{Habitat}(s_i)\right) \] where Habitat\((s_i)\) is a binary indicator of habitat suitability.

The total number of points in region \(A\) follows:

\[ N(A) \sim \text{Poisson}\left(\int_A \lambda(s)\, ds\right) \]

Because the integral \(\int_A \lambda(s)\, ds\) rarely has a closed‑form solution when covariates vary spatially, we approximate it using numerical quadrature. A regular grid of quadrature points \(\{s_j\}_{j=1}^m\) is generated across the study area, and the integral is approximated as:

\[ \int_A \lambda(s)\, ds \approx \sum_{j=1}^m w_j \lambda(s_j) \] where \(w_j\) is the area represented by each quadrature point.

We fitted the IPP using the down‑weighted Poisson regression method (Renner et al. 2015). In this approach:

• Observed telemetry locations are assigned extremely small weights (1e‑6), ensuring they contribute only to the log‑intensity component of the likelihood.
• Quadrature points are assigned weights equal to the area of each grid cell, ensuring they correctly approximate the integral term.

# ------------------------------------------------------------
# Prepare presence and covariates for IPP (full data)
# ------------------------------------------------------------

# Quadrature points (grid) represent zeros
grid_df$presence <- 0

# Observed telemetry points represent presences

tracked_birds$presence <- 1
tracked_birds$area     <- 1e-6   # down‑weight observed points

# Extract covariates for observed points
tracked_birds$habitat <- raster::extract(cov_raster,
                                         tracked_birds[, c("x", "y")])

full_data <- rbind(
  tracked_birds[, c("x", "y", "area", "habitat", "presence")],
  grid_df[,      c("x", "y", "area", "habitat", "presence")])

# Fit IPP using down‑weighted Poisson regression
start_ipp_fulldata <- Sys.time()

dwp_fulldata <- glm(
  presence / area ~ habitat ,
  weights = area,
  data    = full_data,
  family  = poisson())

summary(dwp_fulldata)

## 
## Call:
## glm(formula = presence/area ~ habitat, family = poisson(), data = full_data, 
##     weights = area)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -12.46837    0.01442 -864.58   <2e-16 ***
## habitat       1.00933    0.01685   59.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 923650  on 27999  degrees of freedom
## Residual deviance: 919590  on 27998  degrees of freedom
## AIC: 919594
## 
## Number of Fisher Scoring iterations: 14

end_ipp_fulldata <- Sys.time()

# ------------------------------------------------------------
# Predict intensity surface across grid
# ------------------------------------------------------------
grid_df$predictions_dwp_fulldata <- exp(
  predict(dwp_fulldata, newdata = grid_df)
)
ggplot(grid_df, aes(x = x, y = y)) +
  geom_tile(aes(fill = predictions_dwp_fulldata)) +
  scale_fill_viridis_c(name = "Predictions") +
  coord_equal() +
  theme_minimal()

# ------------------------------------------------------------
# Extract habitat‑selection coefficient and CI
# ------------------------------------------------------------

co      <- coef(summary(dwp_fulldata))
beta_hat <- co["habitat", "Estimate"]
se_hat   <- co["habitat", "Std. Error"]
ci_low   <- beta_hat - 1.96 * se_hat
ci_up    <- beta_hat + 1.96 * se_hat

dwp_full_est <- data.frame(
  model    = "IPP",
  estimate = beta_hat,
  low      = ci_low,
  up       = ci_up)

Autocorrelation and thinning

A key assumption of the IPP is that points occur independently. High‑resolution telemetry data violate this assumption because an animal’s location at time \(t\) is strongly dependent on its location at time \(t-1\). Ignoring this serial autocorrelation can lead to underestimated uncertainty and spurious inference about habitat selection. To reduce temporal dependence and better satisfy the independence assumption, we applied random thinning, retaining each location with probability of 10%. This reduces autocorrelation while preserving broad spatial patterns of space use.

# ------------------------------------------------------------
# Randomly thin telemetry data to reduce temporal autocorrelation
# ------------------------------------------------------------

# Convert observed points to spatstat ppp object
data_ppm <- as.ppp(
  tracked_birds[, c("x", "y")],
  W = c(0, 50000, 0, 50000)
)

# Thin to 10% of points
data_t10 <- rthin(data_ppm, P = 0.10)

# ------------------------------------------------------------
# Build thinned dataset with covariates
# ------------------------------------------------------------
t10 <- data.frame(
  x        = data_t10$x,
  y        = data_t10$y,
  presence = 1,
  area     = 1e-6
)

# Extract habitat covariates
t10$habitat <- raster::extract(cov_raster, t10[, c("x", "y")])

Down weighted Poisson model

# ------------------------------------------------------------
# Prepare thinned dataset + quadrature points
# ------------------------------------------------------------

# Quadrature points (grid) represent zeros
grid_df$presence <- 0

# Combine thinned presences with quadrature points
ppm_data10 <- rbind(
  t10[,     c("x", "y", "area", "habitat","presence")],
  grid_df[, c("x", "y", "area", "habitat", "presence")]
)

# ------------------------------------------------------------
# Fit IPP using down‑weighted Poisson regression
# ------------------------------------------------------------
start_ipp <- Sys.time()

dwp<- glm(
  presence / area ~ habitat ,
  weight = area,
  data   = ppm_data10,
  family = poisson()
)

summary(dwp)

## 
## Call:
## glm(formula = presence/area ~ habitat, family = poisson(), data = ppm_data10, 
##     weights = area)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -14.77889    0.04578 -322.85   <2e-16 ***
## habitat       0.99952    0.05355   18.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 99198  on 11772  degrees of freedom
## Residual deviance: 98805  on 11771  degrees of freedom
## AIC: 98809
## 
## Number of Fisher Scoring iterations: 16

end_ipp <- Sys.time()

# ------------------------------------------------------------
# Extract habitat‑selection coefficient and CI
# ------------------------------------------------------------
co      <- coef(summary(dwp))
beta_hat <- co["habitat", "Estimate"]
se_hat   <- co["habitat", "Std. Error"]
ci_low   <- beta_hat - 1.96 * se_hat
ci_up    <- beta_hat + 1.96 * se_hat

dwp_est <- data.frame(
  model    = "IPP*",
  estimate = beta_hat,
  low      = ci_low,
  up       = ci_up)

# ------------------------------------------------------------
# Predict intensity surface across grid
# ------------------------------------------------------------
grid_df$predictions_dwp <- exp(
  predict(dwp, newdata = grid_df)
)

ggplot(grid_df, aes(x = x, y = y)) +
  geom_tile(aes(fill = predictions_dwp)) +
  scale_fill_viridis_c(name = "Predictions") +
  coord_equal() +
  theme_minimal()

Log‑Gaussian Cox Process (LGCP)

Presence‑only telemetry data can exhibit spatial autocorrelation in addition to serial (temporal) autocorrelation. Spatial dependence may arise from clustering in the point pattern and/or from unmeasured environmental variables that vary smoothly over space. When such residual spatial structure is present, the assumption of conditional independence is violated (e.g. when covariates do not explain the spatial structure in the data) and a simple IPP can result in biased inference. The log‑Gaussian Cox process (LGCP) addresses this by augmenting the log‑intensity with a Gaussian random field (GRF) that acts as a spatially correlated error term. In this formulation, the intensity depends on environmental covariates and on a latent spatial process: \[ \lambda(s_i) \;=\; \exp\{\beta_0 \;+\; \beta_1\,\text{Habitat}(s_i)\ +\;\xi(s_i)\} \] where \(\xi(s_i)\) is a mean‑zero spatial Gaussian process whose covariance decreases with increasing separation between locations. This structure allows the model to capture residual clustering not explained by the observed covariates, improving calibration and reducing spatial confounding.

An efficient way to fit an LGCP is to approximate the GRF with a basis‑function expansion, implemented as a smooth over coordinates (e.g., a Gaussian‑process or thin‑plate spline basis) using mgcv (Dovers, Stoklosa, and Warton 2024). In practice, this corresponds to adding a term such as s(x, y, bs = "gp") to the linear predictor to represent the latent spatial field, while retaining the fixed‑effect covariates (e.g., Habitat, Distance).

To ensure comparability with the IPP, we used the same quadrature scheme and weighting when fitting the LGCP. As before, the data was thinned prior to model fitting to reduce the strength of spatial and temporal dependence, helping to ensure that the remaining clustering can be attributed to ecological processes of interest rather than artefacts in movement.

# ------------------------------------------------------------
# Fit LGCP (Gaussian Process Spatial Random Field)
# ------------------------------------------------------------
start_lgcp <- Sys.time()

lgcp <- gam(
  presence / area ~ habitat  + s(x, y, bs = "gp", k = 100),
  weight = area,
  data   = ppm_data10,
  method = "REML",
  family = poisson()
)

summary(lgcp)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## presence/area ~ habitat + s(x, y, bs = "gp", k = 100)
## 
## Parametric coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -14.97401    0.04952 -302.40   <2e-16 ***
## habitat       0.98572    0.05363   18.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df Chi.sq p-value    
## s(x,y) 71.71  84.97    779  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  -0.00622   Deviance explained = 1.39%
## -REML =  49084  Scale est. = 1         n = 11773

end_lgcp <- Sys.time()

# ------------------------------------------------------------
# Extract habitat‑selection coefficient and CI
# ------------------------------------------------------------
coef_est <- summary(lgcp)

lgcp_est <- data.frame(
  model    = "LGCP*",
  estimate = coef_est$p.coeff["habitat"],
  up       = coef_est$p.coeff["habitat"] + 1.96 * coef_est$se["habitat"],
  low      = coef_est$p.coeff["habitat"] - 1.96 * coef_est$se["habitat"]
)

# ------------------------------------------------------------
# Predict intensity surface across grid
# ------------------------------------------------------------
grid_df$predict_lgcp <- exp(
  predict(lgcp, newdata = grid_df)
)

ggplot(grid_df, aes(x = x, y = y)) +
  geom_tile(aes(fill = predict_lgcp)) +
  scale_fill_viridis_c(name = "Predictions") +
  coord_equal() +
  theme_minimal()

Space–Time Point Process Model (STPP)

While thinning telemetry data is easy to apply, it has the drawback of discarding sometimes large amounts of data. To avoid this, space-time point process model accounts for serial dependence without discarding data (Johnson, Hooten, and Kuhn 2013; Hooten et al. 2017). We model animal movement and habitat selection using a spatio‑temporal point process (STPP) framework. Two model formulations were considered: a full (conditional) STPP and a marginalised STPP.

Movement kernel choice

As mentioned above, the STPP requires a movement kernel \(f(s_t \mid s_{t-1})\) that defines the relative probability of moving between locations over a single time-step.

The Brownian motion (BM) and Ornstein–Uhlenbeck (OU) movement kernels used in the STPP were chosen due to their wide usage in movement ecology and practicality of implementation. They have well-defined transition densities, can be evaluated easily at arbitrary quadrature points, and provide a clear separation between movement (availability) and habitat selection. However, they are not intended to exactly reproduce the movement process used to simulate the tracks, which were generated from a complex, state-dependent movement process, incorporating gamma-distributed step lengths, directional persistence and habitat-dependent switching between behavioural states. Instead, their role within the STPP is to determine which locations are available at a given step, while the habitat covariates estimate selection among those available locations. The BM and OU STPPs therefore provide a framework for estimating habitat selection conditional on simplified movement constraints.

The BM model represents unbiased, diffusive movement, where transitions depend only on distance from the previous location. The OU model extends this by introducing attraction towards a central location, providing a simple approximation to central-place behaviour. Further details on each movement model can be found in the sections below.

Full spatio-temporal point process model

The full spatio-temporal point process model is formulated at the scale of individual movement steps. At each time step \(t\), the animal’s location is modelled as a realisation from a conditional intensity function that depends on the previous observed location \(s_{t-1}\).

The conditional intensity is written as

\[ \lambda(s_t \mid s_{t-1}) = f(s_t \mid s_{t-1}) \cdot \exp\{x(s_t)^\top \beta\}, \]

where:

• \(f(s_t \mid s_{t-1})\) is the movement kernel, defining the set of locations available from \(s_{t-1}\),
• \(x(s_t)\) are habitat covariates evaluated at location \(s_t\),
• \(\beta\) are habitat-selection coefficients.

Taking logs gives:

\[ \log \lambda(s_t \mid s_{t-1}) = \log f(s_t \mid s_{t-1}) + x(s_t)^\top\beta \]

This formulation separates movement and habitat selection. The movement kernel defines which locations are available at each step, while the covariates estimate selection among those available locations.

Because the intensity cannot be evaluated over continuous space, the likelihood is approximated numerically using a quadrature scheme. For each observed step, the continuous spatial domain is replaced by a discrete set of quadrature points, and the model is fitted as a weighted regression over these points, where the movement kernel enters as an offset term.

Spatial-temporal quadrature grid

For each observed step, a spatial grid is constructed around the previous observed location \(s_{t-1}\). This grid represents the set of locations that are available to the animal over the interval \(\Delta t\).

The extent of the grid is defined using the empirical movement data. Specifically, the search radius is determined from the maximum observed displacement per unit time, ensuring that the quadrature domain includes only the region where the movement kernel has a meaningful influence. This avoids evaluating the model over an unnecessarily large spatial domain while capturing all plausible movements.

Within this search radius, space is discretised into a regular grid with cell dimensions defined by tile.dim. Each quadrature point corresponds to the centroid of a grid cell. The grid resolution therefore controls the balance between computational cost and approximation accuracy: finer grids provide a closer approximation to continuous space, while coarser grids reduce computational burden.

The observed location at time \(t\) is included explicitly as the first quadrature point and assigned \(\text{presence} = 1\). All other grid points are treated as pseudo-absence points with \(\text{presence} = 0\).

Each quadrature point has a spatio-temporal weighting, defined as the product of the grid cell area and the time interval \(\Delta t\). This ensures that the discrete likelihood provides a consistent approximation to the continuous space-time point process.

The total number of quadrature points varies across steps and depends on both the grid resolution and the movement-based search radius. Steps with larger potential displacement generate larger quadrature domains and therefore more quadrature points.

At each step, the movement kernel \(f(s_t \mid s_{t-1})\) is evaluated at all quadrature points to define their relative availability. Habitat covariates are extracted at the same locations, and the resulting quadrature data form the input to the STPP likelihood.

Under this construction, the observed location is treated as a realisation from the conditional intensity, while the pseudo-absence points represent locations that were available but not used.

SpatTempQuad <- function(track.data, time.col,
                         tile.dim,
                         mu_hat = NULL,
                         Lambda_hat = NULL,
                         phi_hat = NULL,
                         movement = c("OU", "BM")) {
  
  movement <- match.arg(movement)
  jitter   <- 1e-6
  
  # distance scaling: metres -> km
  scale_dist <- 1000
  
  # ------------------------------------------------------------
  # 1. Time
  # ------------------------------------------------------------
  times   <- track.data@data[, time.col, drop = TRUE]
  delta.t <- unique(diff(times))
  
  if (length(delta.t) != 1) stop("Time interval is not constant.")
  delta.t <- delta.t[1]
  
  # ------------------------------------------------------------
  # 2. Coordinates (metres, UTM)
  # ------------------------------------------------------------
  coords <- as.matrix(track.data@coords)
  
  # ------------------------------------------------------------
  # 3. Step distances (for grid extent only)
  # ------------------------------------------------------------
  step.dist <- sqrt(rowSums((coords[-1, ] - coords[-nrow(coords), ])^2))
  max.disp  <- max(step.dist, na.rm = TRUE)
  
  disp.rad.lookup <- list(
    x = max.disp / tile.dim[1],
    y = max.disp / tile.dim[2]
  )
  
  # ------------------------------------------------------------
  # 4. Quadrature loop
  # ------------------------------------------------------------
  df <- NULL
  
  # previous observed location
  mu <- coords[1, ]
  
  for (i in 2:nrow(coords)) {
    
    x_obs <- coords[i, 1]
    y_obs <- coords[i, 2]
    
    rad.x <- ceiling(disp.rad.lookup$x)
    rad.y <- ceiling(disp.rad.lookup$y)
    
    x.grid <- seq(
      mu[1] - rad.x * tile.dim[1],
      mu[1] + rad.x * tile.dim[1],
      tile.dim[1]
    )
    
    y.grid <- seq(
      mu[2] - rad.y * tile.dim[2],
      mu[2] + rad.y * tile.dim[2],
      tile.dim[2]
    )
    
    spqt <- SpatialPoints(
      expand.grid(x = x.grid, y = y.grid),
      proj4string = CRS(SRS_string = "EPSG:32630")
    )
    
    obs_pt <- SpatialPoints(
      cbind(x_obs, y_obs),
      proj4string = CRS(SRS_string = "EPSG:32630")
    )
    
    spqt <- rbind(obs_pt, spqt)
    
    out <- as(spqt, "data.frame")
    names(out)[1:2] <- c("x", "y")
    
    out$time     <- times[i]
    out$presence <- c(1, rep(0, nrow(out) - 1))
    out$area     <- prod(tile.dim)
    out$volume   <- out$area * delta.t
    out$response <- out$presence / out$volume
    
    # ------------------------------------------------------------
    # 5. Movement kernel
    # ------------------------------------------------------------
    if (movement == "BM") {
      
      # distances in metres
      d_m <- spDistsN1(
        as.matrix(out[, c("x", "y")]),
        mu,
        longlat = FALSE
      )
      
      # convert to km
      d_km <- d_m / scale_dist
      
      # Brownian kernel with D = 0.5
      out$moveKern <- exp(-0.5 * d_km^2 / delta.t)
      
    } else if (movement == "OU") {
      
      if (is.null(mu_hat) || is.null(Lambda_hat) || is.null(phi_hat)) {
        stop("mu_hat, Lambda_hat, phi_hat required for OU")
      }
      
      r <- exp(-delta.t / phi_hat)
      
      mu_t <- mu_hat + r * (mu - mu_hat)
      
      Sigma_t <- (1 - r^2) * Lambda_hat + diag(jitter, 2)
      
      out$moveKern <- mvtnorm::dmvnorm(
        x     = as.matrix(out[, c("x", "y")]),
        mean  = mu_t,
        sigma = Sigma_t,
        log   = FALSE
      )
    }
    
    # clean numerical issues
    out$moveKern[!is.finite(out$moveKern)] <- 0
    
    # ------------------------------------------------------------
    # 6. Update previous location
    # ------------------------------------------------------------
    mu <- c(x_obs, y_obs)
    
    df <- rbind(df, out)
  }
  
  return(df)
}

The function returns a quadrature data frame containing:

• quadrature coordinates,
• presence/absence indicator,
• spatial and temporal weights,
• kernel values.

Brownian Movement Kernel

Under a Brownian motion model, the location at time \(t\) depends on the previous location \(s_{t-1}\) through a Gaussian transition density:

\[ s_t \mid s_{t-1} \sim \mathcal N(s_{t-1}, 2D\Delta t I) \]

where \(s_t\) is the animal location at time \(t\), \(s_{t-1}\) is the previous observed location, \(\Delta t\) is the time interval between observations, \(D\) is the Brownian diffusion coefficient, and \(I\) is the \(2 \times 2\) identity matrix.

This gives the Brownian movement kernel:

\[ f_{\text{BM}}(s_t \mid s_{t-1}) \propto \exp\left( -\frac{\|s_t - s_{t-1}\|^2}{4D\Delta t} \right) \]

Using the Gaussian distribution, and following Johnson, Hooten, and Kuhn (2013), we fixed the diffusion coefficient at \(D = 0.5\). Substituting this value gives:

\[ f_{\text{BM}}(s_t \mid s_{t-1}) \propto \exp\left( -0.5 \frac{\|s_t - s_{t-1}\|^2}{\Delta t} \right) \]

As telemetry coordinates were recorded in metres, direct use of squared metre-scale distances resulted in an overly concentrated movement kernel relative to the spatial resolution of the quadrature grid (here 500 m). This happened even when the nearest available grid points received negligible movement weight, leading to numerical instability in the model.

To address this, distances were rescaled from metres to kilometres when evaluating the Brownian kernel. Let

\[ d_{\mathrm{km}}(s_t, s_{t-1}) = \frac{\|s_t - s_{t-1}\|}{1000} \]

The movement kernel used in the full STPP was therefore:

\[ f_{\text{BM}}(s_t \mid s_{t-1}) \propto \exp\left( -0.5 \frac{d_{\mathrm{km}}(s_t, s_{t-1})^2}{\Delta t} \right) \]

This rescaling affects only the numerical evaluation of the movement kernel and does not alter the spatial coordinates, covariate values, or quadrature construction. Instead, it ensures that the Brownian kernel is compatible with the discretised spatial grid, yielding a smooth and well-resolved availability surface.

The resulting movement kernel was included in the full STPP as an offset:

\[ \log \lambda(s_t \mid s_{t-1}) = \log f_{\text{BM}}(s_t \mid s_{t-1}) + x(s_t)^\top \beta + g(x_t, y_t) \]

where \(\lambda(s_t \mid s_{t-1})\) is the conditional intensity, \(x(s_t)\) are habitat covariates evaluated at quadrature locations, \(\beta\) are habitat-selection coefficients, and \(g(x_t, y_t)\) is a spatial smooth term.

Thus, the Brownian movement kernel defines step-level availability, while the habitat covariates estimate selection among locations available under this movement process.

full_bm_start<-Sys.time()
# ------------------------------------------------------------
# Prepare track combinations
# ------------------------------------------------------------

ids<-unique(tracked_birds$ID)

# ------------------------------------------------------------
# Compute Brownian movement kernel 
# ------------------------------------------------------------
quad_BM<-NULL

for (j in 1:length(ids)) {
  
  bird_id  <- ids[j]
  
  track_data <- tracked_birds %>%
    filter(ID == bird_id)
  
  #transform to spatial data
  coordinates(track_data) <- ~x+y
  
  temp_BM <- SpatTempQuad(track.data = track_data, 
                          time.col = "time", 
                          tile.dim = c(500, 500), # previous grid resolution 
                          movement   = "BM")
  
  temp_BM$ID<- bird_id
  
  quad_BM<-rbind(quad_BM, temp_BM)
}

quad_BM$log_moveKern<-log(quad_BM$moveKern +  1e-10)

ggplot(quad_BM) +
  geom_point(aes(x = x, y = y, colour = log_moveKern)) +
  scale_color_viridis_c() +
  theme_minimal()

# ------------------------------------------------------------
# Extract habitat value
# ------------------------------------------------------------
quad_BM$cov <- raster::extract(cov_raster, cbind(quad_BM$x, quad_BM$y))

# ------------------------------------------------------------
# Run model
# ------------------------------------------------------------

STPP_BM<- gam(response ~ cov+ s(x,y, fx=T), offset= log_moveKern, weights =volume, data=quad_BM, family=quasi(link=log, variance=mu), method="REML")

summary(STPP_BM)

## 
## Family: quasi 
## Link function: log 
## 
## Formula:
## response ~ cov + s(x, y, fx = T)
## 
## Parametric coefficients:
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) -15.98696    0.01416 -1128.69   <2e-16 ***
## cov           0.74619    0.01653    45.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##        edf Ref.df     F  p-value    
## s(x,y)  29     29 3.011 9.38e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0458   Deviance explained = 2.01%
## -REML = -5.8389e+06  Scale est. = 0.95724   n = 865777

STPP_est_BM<-data.frame(model="STPP BM",
                        estimate= summary(STPP_BM)[["p.coeff"]][["cov"]],
                        up = summary(STPP_BM)[["p.coeff"]][["cov"]] + 1.96*summary(STPP_BM)[["se"]][["cov"]],
                        low = summary(STPP_BM)[["p.coeff"]][["cov"]] - 1.96*summary(STPP_BM)[["se"]][["cov"]])

full_bm_end<-Sys.time()

# ------------------------------------------------------------
# Plot intensity
# ------------------------------------------------------------
#Window and point pattern
win <- owin(range(quad_BM$x), range(quad_BM$y))

pp <- ppp(
  x = quad_BM$x,
  y = quad_BM$y,
  window = win,
  marks = quad_BM$moveKern
)

# Kernel surface on a grid
dens_im <- density.ppp(
  pp,
  weights = quad_BM$moveKern,
  sigma = 0.1,   # set a fixed value for increase speed
  dimyx = c(100, 100)      # grid resolution
)

# Convert to data frame
dens_df <- as.data.frame(dens_im)

# Set as grid
grid_stpp <- dens_df
names(grid_stpp) <- c("x", "y", "moveKern_bm")

# Log version
grid_stpp$log_moveKern_bm<- log(grid_stpp$moveKern_bm + 1e-10)

# Plot final intensity surface
ggplot(grid_stpp, aes(x = x, y = y)) +
  geom_tile(aes(fill = log_moveKern_bm)) +
  scale_fill_viridis_c(name = "Predictions") +
  coord_equal() +
  theme_minimal()

# Extract habitat covariate values for each grid cell
grid_stpp$cov <- raster::extract(cov_raster, cbind(grid_stpp$x, grid_stpp$y))

# Get linear predictor from GAM (excludes offset)
eta_bm <- predict(STPP_BM, newdata = grid_stpp, type = "link")

# Combine habitat + movement, then convert to intensity scale
grid_stpp$pred_bm <- exp(eta_bm + grid_stpp$log_moveKern_bm)

# Plot final intensity surface
ggplot(grid_stpp, aes(x = x, y = y)) +
  geom_tile(aes(fill = pred_bm)) +
  scale_fill_viridis_c(name = "Predictions") +
  coord_equal() +
  theme_minimal()

Ornstein–Uhlenbeck (OU) movement kernel

The Ornstein–Uhlenbeck (OU) process (Johnson et al. 2008) provides a continuous-time model of movement in which locations evolve with temporal correlation and are attracted toward a central point. For two successive observations at times \(t-1\) and \(t\), the OU model specifies

\[ s_t \mid s_{t-1} \sim \mathcal{N}(\mu_t,\ \Lambda_t), \]

with conditional mean

\[ \mu_t = \mu + B_t (s_{t-1} - \mu), \]

and temporal decay matrix

\[ B_t = \begin{pmatrix} e^{-\Delta t / \phi} & 0 \\ 0 & e^{-\Delta t / \phi} \end{pmatrix}, \]

where \(\Delta t\) is the time interval and \(\phi\) is the autocorrelation timescale. The conditional covariance is

\[ \Lambda_t = \Lambda - B_t \Lambda B_t. \]

Here, \(\mu\) is defined as the mean \(x\) and \(y\) coordinates and \(\Lambda\) is the stationary covariance describing the spatial spread of movement around \(\mu\).

Estimation of \(\phi\) and \(\Lambda\)

To characterise overall spatial variability, a population-level covariance matrix \(\widehat{\Lambda}\) was calculated using all locations combined after centering on \(\mu\). This represents the overall spatial extent of movement across all tracks.

To estimate the autocorrelation timescale \(\phi\), parameters were fitted at the level of individual tracks. For each track, locations were centred on \(\mu\), and a track-specific covariance matrix \(\Lambda_i\) was computed. This track-level covariance was used when evaluating the OU transition likelihood, ensuring that the spatial scale of movement was correctly represented for each individual.

The autocorrelation parameter \(\phi_i\) was estimated separately for each track by maximising the OU transition likelihood based on consecutive observations. The resulting track-level estimates were then summarised using the median to obtain a population-level value \(\widehat{\phi}\). The final OU movement kernel used in the STPP therefore combines the population-level covariance matrix \(\widehat{\Lambda}\) with the population-level autocorrelation parameter \(\widehat{\phi}\).

A small diagonal jitter term was added to step-wise covariance matrices to ensure numerical stability and maintain positive definiteness.

The OU movement kernel was incorporated into the STPP as an offset term in the model formulation. For each quadrature point, the transition density \(k_{\text{OU}}(s_t|s_{t-1})\) was evaluated and included in log form as

\[ k_{\text{OU}}(s_t|s_{t-1}) = \mathcal N(s_t,\mu_t,\Lambda_t) \]

full_ou_start<-Sys.time()
# ------------------------------------------------------------
# Split tracks
# ------------------------------------------------------------

track_list <- split(tracked_birds, tracked_birds$ID)

# ------------------------------------------------------------
# Calculate mu_hat
# ------------------------------------------------------------
mu_hat <-c(mean(tracked_birds$x), mean(tracked_birds$y))

# ------------------------------------------------------------
# Calculate phi_hat
# ------------------------------------------------------------
estimate_phi <- function(td) {
  jitter <- 1e-6
  
  # Need at least 3 points
  if (nrow(td) < 3) return(NA_real_)
  # Ensure increasing time
  td <- td[order(td$time), ]
  
  # Center positions around the attraction point
  centered   <- cbind(td$x - mu_hat[1], td$y - mu_hat[2])
  Lambda_hat <- stats::cov(centered)
  
  # Negative log-likelihood as a function of log(phi_min)
  ou_nll_logphi <- function(log_phi) {
    phi_min <- exp(log_phi)   # minutes
    nll <- 0
    for (i in 2:nrow(td)) {
      delta_min <- td$time[i] - td$time[i - 1]       # minutes
      if (!is.finite(delta_min) || delta_min <= 0) next
      
      r <- exp(- delta_min / phi_min)                # unit-consistent
      mu_t <- mu_hat + r * (c(td$x[i - 1], td$y[i - 1]) - mu_hat)
      
      Sigma_t <- (1 - r^2) * Lambda_hat + diag(jitter, 2)
      
      nll <- nll - mvtnorm::dmvnorm(
        x     = c(td$x[i], td$y[i]),
        mean  = mu_t,
        sigma = Sigma_t,
        log   = TRUE
      )
    }
    nll
  }
  
  # Optimise log(phi_min) with sensible bounds (minutes)
  opt <- try(
    optim(par = log(60), fn = ou_nll_logphi, method = "L-BFGS-B",
          lower = log(1e-1), upper = log(1e5)),
    silent = TRUE
  )
  
  if (inherits(opt, "try-error")) return(NA_real_)
  return(exp(opt$par))# phi_hat in minutes
}


phi_vals <- sapply(track_list, estimate_phi)

phi_tbl <- data.frame(
  Track = names(phi_vals),
  phi_hat = phi_vals
)

phi_tbl <- phi_tbl[is.finite(phi_tbl$phi_hat), ]
phi_hat<-median(phi_tbl$phi_hat)

# ------------------------------------------------------------
# Calculate lambda_hat
# ------------------------------------------------------------
# Center locations on the central place
dx <- tracked_birds$x - mean(tracked_birds$x)
dy <- tracked_birds$y -  mean(tracked_birds$y)

# Empirical stationary covariance
Lambda_hat <- cov(cbind(dx, dy))

# ------------------------------------------------------------
# Prepare track combinations
# ------------------------------------------------------------
track_combos <- unique("tracked_birds$ID")

# ------------------------------------------------------------
# Compute OU movement kernel 
# ------------------------------------------------------------
quad_OU <- NULL

for (j in 1:length(ids)) {
  
  bird_id  <- ids[j]
  
  track_data <- tracked_birds %>%
    filter(ID == bird_id)
  
  if (nrow(track_data) < 2) next
  
  coordinates(track_data) <- ~x+y
  
  temp_OU <- SpatTempQuad(
    track.data = track_data,
    time.col   = "time",
    tile.dim   = c(500, 500),
    movement   = "OU",
    mu_hat     = mu_hat,
    Lambda_hat = Lambda_hat,
    phi_hat    = phi_hat
  )
  
  temp_OU$ID       <- bird_id

  quad_OU <- rbind(quad_OU, temp_OU)
}

quad_OU$log_moveKern<-log(quad_OU$moveKern +  1e-10)

ggplot(quad_OU) +
  geom_point(aes(x = x, y = y, colour = log_moveKern)) +
  scale_color_viridis_c() +
  theme_minimal()

# ------------------------------------------------------------
# Extract habitat value
# ------------------------------------------------------------
quad_OU$cov <- raster::extract(cov_raster, cbind(quad_OU$x, quad_OU$y))

# ------------------------------------------------------------
# Run model
# ------------------------------------------------------------
STPP_OU<- gam(response ~ cov + s(x,y, fx=T), offset = log_moveKern, weights =volume, data=quad_OU, family=quasi(link=log, variance=mu), method="REML")

summary(STPP_OU)

## 
## Family: quasi 
## Link function: log 
## 
## Formula:
## response ~ cov + s(x, y, fx = T)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.005436   0.014253  -0.381    0.703    
## cov          0.819626   0.016633  49.276   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##        edf Ref.df     F p-value    
## s(x,y)  29     29 5.521  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0239   Deviance explained = 2.25%
## -REML = -5.7948e+06  Scale est. = 0.96947   n = 865777

STPP_OU_est<-data.frame(model="STPP OU",
                        estimate= summary(STPP_OU)[["p.coeff"]][["cov"]],
                        up = summary(STPP_OU)[["p.coeff"]][["cov"]] + 1.96*summary(STPP_OU)[["se"]][["cov"]],
                        low = summary(STPP_OU)[["p.coeff"]][["cov"]] - 1.96*summary(STPP_OU)[["se"]][["cov"]])
full_ou_end<-Sys.time()

# ------------------------------------------------------------
# Plot full OU intensity
# ------------------------------------------------------------

# Compute OU movement kernel on grid

grid_stpp$moveKern_ou <- dmvnorm(
  cbind(grid_stpp$x, grid_stpp$y),
  mean  = mu_hat,
  sigma = Lambda_hat
)

# Avoid log(0)
grid_stpp$log_moveKern <- log(grid_stpp$moveKern_ou + 1e-10)

# Predict without offset
eta_ou <- predict(STPP_OU, newdata = grid_stpp, type = "link")

# Manually add offset (movement kernel)
grid_stpp$pred_ou<- exp(eta_ou + grid_stpp$log_moveKern)

#plot predictions
ggplot(grid_stpp, aes(x, y)) +
  geom_tile(aes(fill = pred_ou)) +
  scale_fill_viridis_c(name = "Predictions") +
  coord_equal() +
  theme_minimal()

Marginalised Space–Time Point Process Model (mSTPP)

In this framework, the telemetry data are treated as a space–time point process, and inference is performed after integrating out the time component. This marginalisation reduces the model to a purely spatial point process, allowing it to be fitted using standard Poisson regression methods. To account for serial dependence in movement, a spatial movement kernel is used to captures the expected spatial distribution of locations generated by the movement process alone, effectively representing long-term availability across space.

The model is implemented using the following steps:

1. Divide the study area into a grid of spatial cells, \(s\).
2. Evaluate the movement kernel at each cell centroid to obtain a spatial availability surface.
3. Aggregate telemetry locations into counts per grid cell.
4. Fit a Poisson GAM with spatial smooths, including the movement kernel as an offset.

The resulting model can be written as

\[ \log \lambda(s) = \log k(s) + x(s)^\top\beta \]

where \(k(s)\) is the marginal movement kernel, \(x(s)\) are habitat covariates and \(\beta\) are selection coefficients.

Including the movement kernel as an offset ensures that the expected spatial pattern arising from movement is accounted for explicitly, rather than being estimated from the data. This avoids confounding between movement behaviour and habitat selection, allowing the covariate effects to be interpreted as selection relative to the availability generated by the movement process.

Brownian Movement Kernel

For the marginalised spatio-temporal point process model, movement was incorporated through a spatial movement kernel that describes expected space use under Brownian motion. Brownian motion does not have a stationary spatial distribution and so the marginal Brownian movement kernel must be constructed by summing movement contributions from previous observed locations through time.

For a Brownian motion process, the conditional movement density from a previous location \(s_{t-1}\) to a location \(s\) depends on the squared Euclidean distance between the two locations and the elapsed time \(\Delta_t\). In principle, the marginal Brownian kernel can be obtained by integrating this transition density over time. This integration yields a closed-form expression involving the incomplete gamma function.

However, evaluation of the exact expression can be computationally expensive and numerically unstable for telemetry data, especially when distances are very small, time intervals are short, or the argument of the special function becomes extreme. Therefore, we used a stable exponential approximation to the marginal Brownian kernel. This approximation preserves the key behaviour of Brownian movement: locations closer to the previous observed location receive higher movement weight, while locations farther away receive lower movement weight.

For each observed movement step, we calculate the movement kernel

\[ \gamma_t(s) = \exp\left[ -\frac{\alpha d(s, s_{t-1})^2}{\Delta_t} \right] \]

where \(d(s, s_{t-1})^2\) is the squared Euclidean distance between the quadrature location \(s\) and the observed location \(s_{t}\), and \(\Delta_t\) is the time interval between two consecutive observations.

The parameter \(\alpha\) controls the spatial decay of the Brownian movement kernel. Larger values of \(\alpha\) produce a narrower movement kernel, meaning that availability declines rapidly with distance from the previous location. Smaller values of \(\alpha\) produce a broader kernel, meaning that locations farther from the previous location retain non-negligible movement weight.

Because Brownian motion has no stationary distribution, the marginal Brownian model requires an externally defined movement scale. Following the approach of Hooten et al. (2017), this scale was derived from observed movement velocities, \(\sigma=\sigma_{1}...\sigma_t\). Specifically, Hooten et al. (2017) used:

\[ \alpha = \frac{1}{(\bar\sigma)^2} \]

In our dataset, however, mean velocity was strongly influenced by numerous small step lengths associated with fine-scale foraging behaviour. This produced a large value of \(\alpha\), which made the movement kernel overly restrictive and concentrated availability too tightly around previous observed locations.

To avoid the marginal movement surface being dominated by slow, tortuous movement behaviour, we instead defined:

\[ \alpha = \frac{1}{\max({\sigma})^2} \]

This choice defines the marginal Brownian movement kernel using the upper bound of observed movement capacity. It therefore gives a more conservative and biologically realistic representation of the spatial area that could be available to the animal over the observed sampling interval.

The marginal Brownian movement kernel was then obtained by summing all movement kernels, \(\gamma_t(s)\) :

\[ k_{\text{BM}}(s) = \sum_{t=2}^{T} \gamma_t(s). \] The movement kernel was included in the marginalised STPP as an offset in the Poisson model:

\[ \log \lambda(s) = \log k_{\text{BM}}(s) + x(s)^\top\beta \]

where \(\lambda(s)\) is the spatial intensity at location \(s\), \(x(s)\) are habitat covariates, and \(\beta\) are habitat-selection coefficients.

Under this formulation, the Brownian movement kernel defines the background availability surface, while the habitat coefficients estimate selection relative to that movement-informed availability.

bm_start<-Sys.time()

# ------------------------------------------------------------
# Compute velocity and step characteristics
# ------------------------------------------------------------
tracked_birds <- tracked_birds %>%
  arrange(ID, time) %>%
  group_by(ID) %>%
  mutate(
    dx  = lead(x) - x,
    dy  = lead(y) - y,
    dt  = lead(time) - time,
    vel = sqrt(dx^2 + dy^2) / dt
  ) %>%
  ungroup()

# ------------------------------------------------------------
# Brownian kernel scaling parameter
# ------------------------------------------------------------
alpha <- 1 / (max(tracked_birds$vel, na.rm = TRUE)^2)

# ------------------------------------------------------------
# Prepare track combinations
# ------------------------------------------------------------
track_combos <- unique("tracked_birds$ID")

# ------------------------------------------------------------
# Compute Brownian movement kernel across grid
# ------------------------------------------------------------
z <- numeric(nrow(grid_df))

for (j in 1:length(ids)) {
  
  bird_id  <- ids[j]
  
  track_data <- tracked_birds %>%
    filter(ID == bird_id)
  
  for (i in 2:nrow(track_data)) {
    
    obs   <- track_data[i - 1, ]
    delta <- track_data$time[i] - track_data$time[i - 1]
    
    d2 <- spDistsN1(
      cbind(grid_df$x, grid_df$y),
      c(obs$x, obs$y),
      longlat = FALSE
    )^2
    
    # Brownian kernel
    kernel <- delta * exp(-alpha * d2 / delta) #-alpha * d2 * gamma_inc(0, alpha * d2 / delta)
    
    kernel[is.nan(kernel)] <- 0
    z <- z + kernel
  }
}

# ------------------------------------------------------------
# Add kernel to grid and visualise
# ------------------------------------------------------------
grid_df$kB    <- z
grid_df$ln_kB <- log(grid_df$kB + 1e-10)

ggplot(grid_df) +
  geom_point(aes(x = x, y = y, colour = ln_kB)) +
  scale_color_viridis_c() +
  theme_minimal()

# ------------------------------------------------------------
# Fit marginalised STPP (Brownian kernel)
# ------------------------------------------------------------
mstpp_brownian <- gam(
  count ~ habitat + offset(ln_kB),
  family = poisson(),
  method = "REML",
  data   = grid_df
)

summary(mstpp_brownian)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## count ~ habitat + offset(ln_kB)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.89510    0.01389 -280.42   <2e-16 ***
## habitat      0.79927    0.01646   48.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.365   Deviance explained = 10.5%
## -REML =  17944  Scale est. = 1         n = 10000

bm_end<-Sys.time()

# ------------------------------------------------------------
# Extract coefficient estimate
# ------------------------------------------------------------

coef_est <- summary(mstpp_brownian)

mstpp_brownian_est <- data.frame(
  model    = "mSTPP (BM)",
  estimate = coef_est$p.coeff["habitat"],
  up       = coef_est$p.coeff["habitat"] + 1.96 * coef_est$se["habitat"],
  low      = coef_est$p.coeff["habitat"] - 1.96 * coef_est$se["habitat"]
)

# ------------------------------------------------------------
# Predict intensity surface
# ------------------------------------------------------------
grid_df$predict_mSTPP <- exp(predict(mstpp_brownian, newdata = grid_df))

ggplot(grid_df, aes(x = x, y = y)) +
  geom_tile(aes(fill = predict_mSTPP)) +
  scale_fill_viridis_c(name = "Predictions") +
  coord_equal() +
  theme_minimal()

Ornstein–Uhlenbeck (OU)

For the Ornstein–Uhlenbeck (OU) process, the marginal distribution of locations has a closed-form stationary distribution:

\[ s \sim \mathcal N(\mu,\; \Lambda), \]

where \(\mu\) is the central location, calculated as mean \(x\) and \(y\) coordinates, and \(\Lambda\) is the stationary covariance matrix. This represents the long-term spatial distribution of the process after temporal dependence has been integrated out.

Because the OU process is mean-reverting and stationary, marginalising over time leads directly to a spatial availability kernel given by

\[ k_{\text{OU}}(s) = \mathcal N(s,\mu,\Lambda) \]

Thus, unlike the Brownian motion case, no explicit integration over time is required, and the autocorrelation timescale \(\phi\) does not appear in the marginal kernel as it does not affect the stationary spatial distribution.

In the marginal STPP, this kernel defines the baseline spatial availability across the study area. Locations closer to \(\mu\) have higher availability, while more distant locations are downweighted according to the covariance structure \(\Lambda\).

In practice, the kernel is evaluated at grid cell centroids, and its logarithm is included as an offset in the Poisson GAM:

\[ \log \lambda(s) = \log k_{\text{OU}}(s) + x(s)^\top\beta \]

where \(x(s)\) are habitat covariates. This formulation separates the long-term spatial pattern generated by movement from habitat selection effects, with the OU kernel defining availability and the covariates explaining habitat selection.

ou_start <- Sys.time()

# ------------------------------------------------------------
# Marginal OU movement kernel (stationary distribution)
# ------------------------------------------------------------

grid_df$kOU <- mvtnorm::dmvnorm(
  x     = as.matrix(grid_df[, c("x", "y")]),
  mean  = mu_hat,
  sigma = Lambda_hat,
  log   = FALSE
)

grid_df$ln_kOU <- log(grid_df$kOU + 1e-10)


# ------------------------------------------------------------
# Plot OU kernel
# ------------------------------------------------------------
ggplot(grid_df, aes(x, y, colour = ln_kOU)) +
  geom_point(size = 1) +
  scale_color_viridis_c() +
  coord_equal() +
  labs(color = "log kOU") +
  theme_minimal()

# ------------------------------------------------------------
# Fit marginalised STPP (OU kernel)
# ------------------------------------------------------------
mstpp_OU <- gam(
  count ~ habitat + offset(ln_kOU),
  family = poisson(),
  method = "REML",
  data   = grid_df
)

summary(mstpp_OU)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## count ~ habitat + offset(ln_kOU)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 21.58185    0.01389 1553.74   <2e-16 ***
## habitat      0.90526    0.01646   54.99   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.00414   Deviance explained = 7.88%
## -REML =  26506  Scale est. = 1         n = 10000

ou_end <- Sys.time()

# ------------------------------------------------------------
# Extract coefficient estimate
# ------------------------------------------------------------
coef_est <- summary(mstpp_OU)

mstpp_OU_est <- data.frame(
  model    = "mSTPP (OU)",
  estimate = coef_est$p.coeff["habitat"],
  up       = coef_est$p.coeff["habitat"] + 1.96 * coef_est$se["habitat"],
  low      = coef_est$p.coeff["habitat"] - 1.96 * coef_est$se["habitat"]
)

# ------------------------------------------------------------
# Predict intensity surface
# ------------------------------------------------------------
grid_df$predict_mSTPP_ou <- exp(
  predict(mstpp_OU, newdata = grid_df)
)

ggplot(grid_df, aes(x = x, y = y)) +
  geom_tile(aes(fill = predict_mSTPP_ou)) +
  scale_fill_viridis_c(name = "Predictions") +
  coord_equal() +
  theme_minimal()

For non‑central place foraging, clear differences were observed when plotting the model predictions for the BM and OU mSTPP. The OU-based mSTPP produced a smooth, symmetric prediction surface centred on \(\mu\), reflecting the stationary Gaussian structure imposed by the model. However, this spatial pattern did not match the underlying movement, which lacked a central attraction. In contrast, the BM-based mSTPP produced a more irregular and spatially heterogeneous prediction surface, better reflecting the diffuse and unconstrained movement patterns in the simulated data. Importantly, it is worth noticing that both approaches yield similar coefficient estimates despite producing very different spatial intensity surfaces. This highlights the importance of specifying an appropriate movement model: a misspecified movement model can impose incorrect spatial structure and therefore misrepresent patterns of space use.

References

Aarts, Geert, John Fieberg, and Jason Matthiopoulos. 2012. “Comparative Interpretation of Count, Presence–Absence and Point Methods for Species Distribution Models.” Methods in Ecology and Evolution 3 (1): 177–87.

Dovers, Elliot, Jakub Stoklosa, and David I Warton. 2024. “Fitting Log-Gaussian Cox Processes Using Generalized Additive Model Software.” The American Statistician 78 (4): 418–25.

Hooten, Mevin B, Devin S Johnson, Brett T McClintock, and Juan M Morales. 2017. Animal Movement: Statistical Models for Telemetry Data. CRC press.

Johnson, Devin S, Mevin B Hooten, and Carey E Kuhn. 2013. “Estimating Animal Resource Selection from Telemetry Data Using Point Process Models.” Journal of Animal Ecology 82 (6): 1155–64.

Johnson, Devin S, Dana L Thomas, Jay M Ver Hoef, and Aaron Christ. 2008. “A General Framework for the Analysis of Animal Resource Selection from Telemetry Data.” Biometrics 64 (3): 968–76.

Lauret, Valentin, Nicolas Courbin, Olivier Scher, and Aurélien Besnard. 2025. “Integrating Individual Tracking Data and Spatial Surveys to Improve Estimation of Animal Spatial Distribution.” Ecosphere 16 (5): e70283.

Matthiopoulos, Jason, John R Fieberg, and Geert Aarts. 2023. “Species-Habitat Associations: Spatial Data, Predictive Models, and Ecological Insights.”

Renner, Ian W, Jane Elith, Adrian Baddeley, William Fithian, Trevor Hastie, Steven J Phillips, Gordana Popovic, and David I Warton. 2015. “Point Process Models for Presence-Only Analysis.” Methods in Ecology and Evolution 6 (4): 366–79.