library(CircStats)
library(Rcpp)
library(RcppArmadillo)
library(sp)
library(geosphere)
library(maps)
library(raster)
library(spatstat)
library(sf)
library(geoR)
library(ggplot2)
set.seed(1234)The simulation domain was defined as a 2,500 km\(^2\) area (50 \(\times\) 50 km), discretised into a regular grid with a spatial resolution of 100 m. Habitat suitability was represented as a binary chequerboard pattern consisting of alternating 2.5 \(\times\) 2.5 km cells, classified as suitable (value = 1) and unsuitable (value = 0).
# Define grid in meters: 100 m resolution over 50 km
x_grid <- seq(0, by = 100, length.out = 500) # 0 to 50,000 meters
y_grid <- seq(0, by = 100, length.out = 500)
# Create coordinate grid
coord_df <- expand.grid(x = x_grid, y = y_grid)
# Checkerboard pattern with 2.5 km squares
x_band <- floor(coord_df$x / 2500) %% 2
y_band <- floor(coord_df$y / 2500) %% 2
df <- data.frame(
x = coord_df$x,
y = coord_df$y,
value = (x_band == y_band) + 0
)
df$value <- (df$value - min(df$value)) / (max(df$value) - min(df$value))
cov_raster <- rasterFromXYZ(df)
crs(cov_raster) <- "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs"
# Save cov_raster
save(cov_raster, file = "cov_raster.RData")
ggplot(df, aes(x = x, y = y, fill = value)) +
geom_tile() +
labs(fill = "Value") +
xlab("X coordinates")+
ylab("Y coordinates") +
scale_fill_viridis_c()+
coord_fixed(expand = FALSE) +
theme_bw() Animal movement was modelled using a three-state behavioural framework designed to capture the main movement modes of non-central place foragers, which do not undertake directed trips to or from a central location. Specifically, State 1 (travelling) represents fast, directed movement; State 2 (search) corresponds to moderately fast movement with some directional persistence; and State 3 (foraging) reflects slower, more tortuous movement associated with intensive resource use. Transitions between behavioural states were influenced by habitat conditions, such that individuals were more likely to enter search (State 2) or foraging (State 3) when in suitable habitat, and more likely to leave the foraging state as habitat suitability decreased or as time increased.
Transitions between behavioural states were governed by a probabilistic matrix informed by Euclidean distance to the colony, time since departure from the nest, and habitat suitability. Transition probabilities were computed using a multinomial logit formulation, with parameters chosen to reflect biologically realistic behaviour—for example, an increased probability of switching to foraging in suitable habitat (Michelot et al. 2017).
\[ T = \begin{pmatrix} 1 - p_{12} & p_{12} & 0 \\ 0 & 1 - p_{23} & p_{23} \\ r_{31} & r_{32} & r_{33} \end{pmatrix} \]
where \(p_{12}\) corresponds to switching from travelling \(\rightarrow\) searching, \(p_{23}\) searching \(\rightarrow\) foraging and \(r_{31}, r_{32}, r_{33}\) the transitions from foraging to other states.
We defined:
\[p_{12} = \text{logit}^{-1}(\eta_{12})\] \[p_{23} = \text{logit}^{-1}(\eta_{23})\] \[r_{31} = \frac{e^{\eta_{31}}}{1 + e^{\eta_{31}} + e^{\eta_{32}}}\] \[r_{32} = \frac{e^{\eta_{32}}}{1 + e^{\eta_{31}} + e^{\eta_{32}}}\] \[r_{33} = 1 - r_{31} - r_{32}\]
The linear predictors were specified as:
\[ \eta_{12} = -2 + 1.5 \cdot \text{Habitat} \]
\[ \eta_{23} = -1.5 + 3 \cdot \text{Habitat} \]
\[ \eta_{31} = -3 - 3 \cdot \text{Habitat} + 0.05 \cdot \text{Time} \]
\[ \eta_{32} = -3 - 1.5 \cdot \text{Habitat} + 0.05 \cdot \text{Time} \]
# -------------------------------------------------------------
# Transition matrix
# ------------------------------------------------------------
addLogit_rcpp <- function(x1, x2) {
p1 <- exp(x1)
p2 <- exp(x2)
res <- c(p1 / (1 + p1 + p2), p2 / (1 + p1 + p2), 1 - (p1 / (1 + p1 + p2)) - (p2 / (1 + p1 + p2)))
return(res)
}
trMatrix_rcpp <- function(beta, time, cov) {
trMat <- matrix(ncol = 3, nrow = 3)
# From State 1
p12 <- plogis(beta[1] + beta[2] * cov)
trMat[1, 1] <- 1 - p12
trMat[1, 2] <- p12
trMat[1, 3] <- 0
# From State 2
p23 <- plogis(beta[3] + beta[4] * cov)
trMat[2, 1] <- 0
trMat[2, 2] <- 1 - p23
trMat[2, 3] <- p23
# From State 3
p31 <- beta[5] + beta[6] * cov + beta[7] * time
p32 <- beta[8] + beta[9] * cov + beta[10] * time
r3 <- addLogit_rcpp(p31, p32)
trMat[3, 1] <- r3[1]
trMat[3, 2] <- r3[2]
trMat[3, 3] <- r3[3]
return(trMat)
}
# HMM transition coefficients
beta <- c(
-2, 1.5, # State 1 → 2
-1.5, 3, # State 2 → 3
-3, -3, 0.05, # State 3 → 1 (less likely to leave)
-3, -1.5, 0.05 # State 3 → 2 (less likely to leave)
)Higher habitat suitability increased the probability of transitioning into search and foraging states, whereas the likelihood of leaving the foraging state increased with time spent within a patch. This transition structure captures realistic patch-use dynamics typical of non-central place foragers.
State-specific behaviours process were modelled as correlated random walks, with associated parameters varying among behavioural states. Step lengths and turning angles were modelled using gamma and von Mises distributions, respectively, with state-specific parameter values.
| State | Mean step length (m) | SD (m) | von Mises \(\kappa\) | Turning angle mean |
|---|---|---|---|---|
| 1 (Travel) | 571.8 | 195.4 | 14 | 0 |
| 2 (Search) | 350.0 | 190.0 | 5 | \(\pi\) |
| 3 (Foraging) | 108.0 | 131.2 | 0.84 | 0 |
New locations were then calculated from sampled step lengths and
turning angles using great-circle geometry implemented in the R package
geosphere (Hijmans 2026).
mu <- c(571.8, 350, 108)
sigma <- c(195.4, 190, 131.2)
am <- c(0, -0.00647766, pi)
kappa <- c(8, 5, 0.84)
nbStates <- 3At each time step, we drew 10,000 candidate steps from the movement kernel associated with the current behavioural state. We calculated a relative selection weight for each candidate by exponentiating the linear predictor and then selected one end point probabilistically as the animal’s next location. Each candidate step \(i\) was evaluated using a log‑linear model of the form:
\[ \log(w_i) = \beta_1 \times L_i + \beta_2 \times \log L_i + \beta_3 \times h_i \]
\[ w_i = \exp(\eta_i) \quad \mathrm{Pr}(i) = \frac{\exp(\eta_i)}{\sum_j \exp(\eta_j)} \]
where \(L_i\) is the step length, \(h_i\) is the habitat covariate and \(\beta_1, \beta_2\) and \(\beta_3\) are the corresponding selection coefficients.
The coefficients for step length and log step length were derived analytically from the gamma distribution parameters (Klappstein, Thomas, and Michelot 2023; Avgar et al. 2016) while the habitat coefficient \(\beta_3\) was varied systematically across simulations with values ranging from 0.5 to 1.5 (in 0.25 increments) as well as 0 (no habitat selection) to assess how models performed under different selection strengths.
Individuals were initialised at a random location within the simulation area at time \(t = 0\), with the first step generated by drawing a direction \(\phi_0 \sim \mathrm{von\ Mises}(0, 1)\) and a step length \(L_0 \sim \mathrm{Gamma}(k, \theta)\). The latter was parameterised in terms of the state-specific mean \(\mu\) and standard deviation \(\sigma\) for travel, such that \(k = \mu^2 / \sigma^2\) and \(\theta = \sigma^2 / \mu\). The resulting displacement defines the first location.
### SSF parameters
# Equations exist to determine the scale and shape parameters of a gamma distribution
# based on its mean and sd
scale <- sigma^2 / mu
shape <- mu^2 / sigma^2
# The equation presented in Appendix A gives the relationship between the
# gamma distribution parameters and the selection coefficients (referred to as β1 and β2 for L and log(L) respectively).
# β1 −1/scale and β2 = shape − 2.
beta1 <- -1 / scale
beta2 <- shape - 2
cov_beta <- 0 #0.5, 0.75, 1, 1.25, 1.5, 1.75,2We simulate a full population of 1,000 individuals and one track per individual with duration of 360 minutes.
nbTracks <- 1
nbID <- 1000
for(c in 1:length(cov_beta)){
cov_gamma <- rep(1, 3)
cov_coefs <- cov_beta[c] * cov_gamma
ssf_betas <- matrix(c(beta1, beta2, kappa, cov_coefs),
ncol = nbStates, byrow = TRUE)
ssf_formula <- ~ cov
n_zeros <- 10000
data <- NULL
for (id in 1:nbID) {
print(paste0("Beta: ", cov_beta[c],"; Individual: ",id ))
track_data <- NULL
state <- 1
phi <- rvm(1, 0, 1)
len <- rgamma(1,
shape = mu[1]^2 / sigma[1]^2,
scale = sigma[1]^2 / mu[1])
centre <- c(sample(seq(500,48500,100),1),
sample(seq(500,48500,100),1))
newLoc <- centre + c(len * cos(phi), len * sin(phi))
while(!newLoc[1] >= 100 & newLoc[1] <= 49900 &
newLoc[2] >= 100 & newLoc[2] <= 49900){
centre <- c(sample(seq(500,48500,100),1),
sample(seq(500,48500,100),1))
newLoc <- centre + c(len * cos(phi), len * sin(phi))
}
angle <- rvm(1, mean = am[state], k = kappa[state])
track_data <- rbind(track_data, data.frame(
x = centre[1], y = centre[2],
step = NA, angle = NA, bearing = NA,
state = 1, time = 0, ID = id))
time <- 0
while (time < 360) {
time <- time + 1
phi0 <- phi
cov <- raster::extract(cov_raster, matrix(newLoc, nrow = 1))
trMat <- trMatrix_rcpp(beta, time, cov)
state <- sample(1:nbStates, 1, prob = trMat[state, ])
angle_candidates <- rvm(n_zeros, mean = am[state], k = kappa[state])
phi_candidates <- phi + angle_candidates
len_candidates <- rgamma(n_zeros,
shape = mu[state]^2 / sigma[state]^2,
scale = sigma[state]^2 / mu[state])
dx <- len_candidates * cos(phi_candidates)
dy <- len_candidates * sin(phi_candidates)
endpoints <- cbind(newLoc[1] + dx, newLoc[2] + dy)
inside <- endpoints[,1] >= 100 & endpoints[,1] <= 49900 &
endpoints[,2] >= 100 & endpoints[,2] <= 49900
while (sum(inside) == 0) {
# sample new turning angles centred on reverse direction
angle_candidates <- rvm(
n_zeros,
mean = pi,
k = kappa[state] # or a separate kappa_back
)
phi_candidates <- phi + angle_candidates
dx <- len_candidates * cos(phi_candidates)
dy <- len_candidates * sin(phi_candidates)
endpoints <- cbind(
newLoc[1] + dx,
newLoc[2] + dy
)
inside <- endpoints[,1] >= 100 & endpoints[,1] <= 49900 &
endpoints[,2] >= 100 & endpoints[,2] <= 49900
}
endpoints <- endpoints[inside, ,drop = FALSE]
angle_candidates <- angle_candidates[inside]
phi_candidates <- phi_candidates[inside]
len_candidates <- len_candidates[inside]
cov_vals <- raster::extract(cov_raster, endpoints)
endpoints.df <- data.frame(
x = endpoints[,1],
y = endpoints[,2],
cov = cov_vals
)
mod_matrix <- model.matrix(ssf_formula, endpoints.df)
mod_matrix <- mod_matrix[, -1, drop = FALSE]
linear_pred <- mod_matrix %*% as.matrix(ssf_betas[4, state])
ssf_prob <- exp(linear_pred)
ssf_prob <- ssf_prob / sum(ssf_prob)
chosen <- sample(seq_len(nrow(endpoints)), 1, prob = ssf_prob)
newLoc <- endpoints[chosen, ]
phi <- phi_candidates[chosen]
len <- len_candidates[chosen]
track_data <- rbind(track_data, data.frame(
x = newLoc[1], y = newLoc[2],
step = len,
angle = phi - phi0,
bearing = phi,
state = state,
time = time,
ID = id
))
}
data <- rbind(data, track_data)
}
# Convert to sf
data_sf <- st_as_sf(data, coords = c("x", "y"), crs = 32630)
# Save only the data object
save(data_sf, file = paste0("NCPF_cov", cov_beta[c], "_", nbID, "ind.RData"))
}
# plot tracks
ggplot() +
geom_tile(data=df, aes(x,y,fill=factor(value)), alpha=0.1)+
geom_path(data= data[data$ID==100,], aes(x, y, col= factor(state), group = ID)) +
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("Travel", "Search", "Foraging"))+
theme_bw() +
coord_fixed()