LGCPs for Distance Sampling
48 Hours of Estimating Animal Abundance — CEAUL @ FCUL, June 2026
Moritz Klaassen
2026-06-09
0. Before we start
The big picture
Today’s session is about doing distance sampling and spatial modelling simultaneously in a single coherent model.
we will instead write down a log-Gaussian Cox process that:
- models the spatial distribution of dolphins as a Log-Gaussian Cox Process (LGCP),
- models imperfect detection via a detection function,
- and estimates everything at once: getting a full posterior for the density surface and for N.
The key insight is that a line transect survey produces a thinned point process: we only see animals that happen to be detected. The thinning probability at perpendicular distance \(y\) is \(g(y)\), the detection function. So the observed intensity is
\[ \lambda_\text{obs}(s, y) = \lambda(s) \times g(y) \]
Installing the packages
Most packages install as usual from CRAN, except INLA — it is hosted on its own repository and must be installed separately. You should have downloaded INLA beforehand.
If you did not install INLA, run this once before the session (it can take a few minutes):
Package requirements
inlabruwrapsINLAfor point process likelihoods. You write the model in a high-level formula language and inlabru handles the integration domains, the thinning, and the SPDE bookkeeping. Full documentation and more examples at https://inlabru-org.github.io/inlabru/.
1. The data: Gulf of Mexico dolphins
Loading and inspecting
The mexdolphin_sf dataset is a line transect survey of
pantropical spotted dolphins (Stenella attenuata) in the Gulf
of Mexico. It ships with dsm.
# Load the sf-format version of the dataset
mexdolphin <- mexdolphin_sf
# What's in there?
names(mexdolphin)## [1] "points" "samplers" "mesh" "ppoly" "lambda" "simulated"
## [7] "depth"The dataset is a list with four main elements:
| Element | What it is |
|---|---|
$ppoly |
The study area polygon (prediction domain) |
$samplers |
The transect lines that were actually surveyed (the “effort”) |
$points |
The detected animal groups, with perpendicular distances |
$mesh |
A triangulation mesh over the study area — needed for the SPDE |
How many detections?
## [1] 4747 groups detected. The key column is distance — the perpendicular distance from the transect line to the detection, in kilometers
Plotting the survey
Let’s see where the transects went and where dolphins were seen.
ggplot() +
# Study area polygon
gg(mexdolphin$ppoly, fill = "steelblue", alpha = 0.08, colour = "steelblue") +
# Transect lines
gg(mexdolphin$samplers, colour = "grey50", linewidth = 0.5) +
# Detection locations
gg(mexdolphin$points, size = 1.5, colour = "firebrick", alpha = 0.9) +
coord_sf(datum = fm_crs(mexdolphin$ppoly)) +
labs(
title = "Gulf of Mexico dolphin survey",
subtitle = "Transects (grey) and detections (red)",
x = NULL, y = NULL
) +
theme_minimal(base_size = 12)
Survey design and detections. Grey lines = transects; dots = dolphin groups.
The triangulation mesh
The SPDE approach represents the Gaussian random field on a mesh of triangles. Let’s see what ours looks like.
ggplot() +
# gg(mesh) fails without a field to evaluate, so we convert to sf first
geom_sf(data = fm_as_sfc(mesh), colour = "grey80", linewidth = 0.2, fill = NA) +
gg(mexdolphin$ppoly, fill = NA, colour = "steelblue", linewidth = 1) +
gg(mexdolphin$samplers, colour = "grey50", linewidth = 0.4) +
coord_sf(datum = fm_crs(mexdolphin$ppoly)) +
labs(
title = "SPDE triangulation mesh",
subtitle = "The Gaussian random field is represented at the mesh nodes",
x = NULL, y = NULL
) +
theme_minimal(base_size = 12)
Triangulation mesh overlaid on the study area.
2. Exploratory look at distances
Before fitting any model, always look at your distance data.
W <- 8 # strip half-width in kilometers
ggplot(mexdolphin$points) +
geom_histogram(
aes(x = distance),
breaks = seq(0, W, length.out = 9), # 8 equal-width bins
boundary = 0,
fill = "steelblue",
colour = "white",
alpha = 0.8
) +
# Rug plot: each tick is one detection
geom_point(aes(x = distance), y = 0.5, shape = "|", size = 3, colour = "firebrick") +
labs(
title = "Distribution of perpendicular distances",
subtitle = paste0("W = ", W, " km | n = ", nrow(mexdolphin$points), " detections"),
x = "Perpendicular distance (km)",
y = "Count"
) +
theme_minimal(base_size = 12)
Histogram of perpendicular distances. The strip half-width is W = 8 km
What to look for
- In priciple :) : Counts should be highest near zero (most detections happen close to the transect line) and decline with distance: this is the hallmark of imperfect detection.
- The shape of the decline tells you something about the detection function: a smooth bell-shaped decline suggests a half-normal; a flatter shoulder near zero with a sharp drop suggests a hazard-rate.
3. The detection function
Defining g(y)
We’ll use the half-normal detection function:
\[g(y) = \exp\!\left(-\frac{y^2}{2\sigma^2}\right)\]
This says: at distance \(y=0\) (right at the transect), detection is certain (\(g(0)=1\)). Detection falls off smoothly and symmetrically as a Gaussian. The parameter \(\sigma\) controls how quickly detection drops with distance.
# The half-normal function.
# distance = the perpendicular distance
# sigma = the scale parameter (controls how fast detection drops)
#
# Note: sigma here is on the NATURAL scale (kilometer),
# but inlabru will estimate it on the log scale internally.
hn <- function(distance, sigma) {
exp(-0.5 * (distance / sigma)^2)
}Setting a prior on sigma
We need to tell the model something sensible about \(\sigma\) before seeing the data. A reasonable choice: animals very close to the transect are almost certainly detected, and detection should fall off over the scale of the strip width \(W\).
We use bm_marginal() to map a standard Normal latent
variable (what INLA uses internally) to an Exponential distribution with
rate \(1/W\), giving a prior with
mean \(\sigma = W =
8\).
# bm_marginal(qfun, pfun, dfun, ...) defines a bijective transformation from
# N(0,1) to some target distribution. Here: Exponential with rate 1/8.
# This gives a prior for sigma with:
# E[sigma] = 8 (the strip half-width)
# P(sigma < 1) ≈ 0.12 (very small sigma is unlikely)
# P(sigma > 20) ≈ 0.08 (very large sigma is possible but not favoured)
sigma_prior <- bm_marginal(qexp, pexp, dexp, rate = 1 / W)4. The spatial model: SPDE + PC priors
The Matérn covariance
The spatial random field \(W(s)\) is modelled with a Matérn covariance function. It has two key parameters:
- Range \(\rho\): the distance beyond which spatial correlation becomes negligible. Think of it as “how far apart two locations can be before knowing the density at one tells you almost nothing about the density at the other.”
- Marginal standard deviation \(\sigma_W\): how much the log-intensity varies spatially (on the log scale).
We set penalised complexity (PC) priors on both parameters. PC priors penalise departure from a simpler base model — here, independence (no spatial structure). They are specified as tail probability statements:
# inla.spde2.pcmatern() creates a Matérn SPDE model on the mesh.
#
# prior.range = c(rho0, alpha1): P(range < rho0) = alpha1
# -> P(range < 50 nmi) = 0.01 (range is very unlikely to be less than 50 nmi)
#
# prior.sigma = c(sigma0, alpha2): P(sd > sigma0) = alpha2
# -> P(sigma_W > 2) = 0.01 (spatial SD > 2 on log scale is very unlikely)
#
# These are weakly informative priors that keep the model from going wild.
# You should always think about what values are plausible for YOUR species and survey.
matern <- inla.spde2.pcmatern(
mesh,
prior.range = c(50, 0.01),
prior.sigma = c(2, 0.01)
)5. Model components
In inlabru, a model is built from
components named effects that together form the linear
predictor. We have three:
# cmp defines the model components (the "ingredients"):
#
# mySPDE(main = geometry, model = matern)
# -> A spatial Gaussian random field W(s).
# main = geometry tells inlabru to use the spatial coordinates of
# each observation as the input to the SPDE.
#
# sigma(1, prec.linear = 1, marginal = sigma_prior)
# -> The detection function scale parameter sigma.
# It's modelled as a single scalar (1 = "intercept-like"),
# with a unit precision on the latent (Normal) scale,
# transformed to Exponential(1/8) via sigma_prior.
#
# Intercept(1)
# -> The overall mean log-intensity mu. Think of it as the baseline
# log-density of animals, before adding the spatial field.
cmp <- ~ mySPDE(main = geometry, model = matern) +
sigma(1, prec.linear = 1, marginal = sigma_prior) +
Intercept(1)6. The model formula
The formula describes how the components combine into the linear predictor for the log-intensity of the observed point process:
\[ \log \lambda_\text{obs}(s, y) = \underbrace{W(s)}_{\text{mySPDE}} + \underbrace{\log g(y|\sigma)}_{\text{log(hn(...))}} + \underbrace{\mu}_{\text{Intercept}} + \underbrace{\log(2)}_{\text{both sides}} \]
# form is the model formula.
#
# Left-hand side: geometry + distance
# -> We observe BOTH spatial location (geometry) and perpendicular distance.
# Both dimensions get their own integration domain.
#
# Right-hand side: the linear predictor
# mySPDE -> spatial random field W(s)
# log(hn(...)) -> log of detection probability at each distance
# Intercept -> overall log-intensity intercept mu
form <- geometry + distance ~ mySPDE +
log(hn(distance, sigma)) +
Intercept + log(2)7. Fitting the model
Now we fit! lgcp() is the main fitting function in
inlabru for Log-Gaussian Cox Processes.
We specify two integration domains:
geometry = mesh: integrate the spatial intensity over the meshdistance = fm_mesh_1d(...): integrate the detection function over distances from 0 to W, using 30 quadrature nodes
# This will take a minute or two to run — INLA is doing approximate
# Bayesian inference over the joint posterior of all parameters.
#
# components = cmp -> the model components defined above
# mexdolphin$points -> the observed detections (response data)
# samplers -> the transect lines (the effort)
# domain -> the integration domains
# formula -> the linear predictor formula
fit <- lgcp(
components = cmp,
mexdolphin$points,
samplers = mexdolphin$samplers,
domain = list(
geometry = mesh,
distance = fm_mesh_1d(seq(0, W, length.out = 30))
),
formula = form
)
# A quick look at the fitted object
summary(fit)## inlabru version: 2.12.0.9022
## INLA version: 24.06.27
## Components:
## mySPDE: main = spde(geometry), group = exchangeable(1L), replicate = iid(1L), NULL
## sigma: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## Observation models:
## Family: 'cp'
## Tag: <No tag>
## Data class: 'sf', 'data.frame'
## Response class: 'numeric'
## Predictor:
## geometry + distance ~ mySPDE + log(hn(distance, sigma)) + Intercept +
## log(2)
## Additive/Linear: FALSE/FALSE
## Used components: effects[mySPDE, sigma, Intercept], latent[]
## Time used:
## Pre = 0.476, Running = 25.2, Post = 2.66, Total = 28.4
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## sigma -0.050 0.209 -0.461 -0.050 0.360 -0.050 0.000
## Intercept -8.156 0.494 -9.287 -8.107 -7.339 -8.034 0.004
##
## Random effects:
## Name Model
## mySPDE SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for mySPDE 295.50 147.305 110.592 262.948 673.84 209.846
## Stdev for mySPDE 0.81 0.251 0.416 0.777 1.39 0.717
##
## Deviance Information Criterion (DIC) ...............: -804.68
## Deviance Information Criterion (DIC, saturated) ....: -923.92
## Effective number of parameters .....................: -904.68
##
## Watanabe-Akaike information criterion (WAIC) ...: 193.11
## Effective number of parameters .................: 10.55
##
## Marginal log-Likelihood: -527.15
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')8. SPDE parameter posteriors
The SPDE has two parameters: the spatial range and the log-variance of the field. Let’s look at their posterior distributions.
# spde.posterior() extracts marginal posteriors for SPDE hyperparameters.
# what = "range" gives the posterior for rho (spatial range).
spde.range <- spde.posterior(fit, "mySPDE", what = "range")
plot(spde.range) +
labs(
title = "Posterior: spatial range",
subtitle = "How far (km) spatial correlation in dolphin density extends",
x = "Range (km)",
y = "Density"
) +
theme_minimal(base_size = 12)
Posterior distribution of the spatial range. Values represent nautical miles.
# what = "log.variance" gives the posterior for log(sigma_W^2).
spde.logvar <- spde.posterior(fit, "mySPDE", what = "log.variance")
plot(spde.logvar) +
labs(
title = "Posterior: log spatial variance",
subtitle = "How much the log-intensity varies spatially",
x = "Log variance",
y = "Density"
) +
theme_minimal(base_size = 12)
Posterior distribution of the log spatial variance.
How to read these plots:
- The range posterior tells you the scale of spatial clustering. A range of 250 means dolphin density is correlated up to about 250km. If the posterior is wide, the data don’t pin down the range very precisely.
- The log variance tells you the amplitude of spatial variation. A large variance means strong spatial heterogeneity (patchy distribution); a small variance means density is relatively uniform.
9. Predicted spatial intensity
Now let’s predict the dolphin density surface. We create a pixel grid
(fm_pixels) over the study area and predict the expected
intensity \(\exp(W(s) + \mu)\) at each
pixel.
# fm_pixels() creates a grid of prediction points inside ppoly.
# dims = c(200, 100) controls the pixel resolution (width x height).
# mask = ppoly ensures we only predict inside the study area.
pxl <- fm_pixels(mesh, dims = c(200, 100), mask = mexdolphin$ppoly)
# predict() evaluates a function of the model components at the new locations.
# ~ exp(mySPDE + Intercept) is the expected intensity:
# lambda(s) = exp(W(s) + mu)
# Note: we DON'T include the detection function here — we want the TRUE
# density surface, not the observed (thinned) surface.
pr.int <- predict(fit, pxl, ~ exp(mySPDE + Intercept))
# Plot the predicted mean intensity
ggplot() +
gg(pr.int, geom = "tile", aes(fill = mean)) +
scale_fill_viridis_c(
name = "Intensity\n(groups/km²)",
option = "plasma",
trans = "sqrt" # sqrt transform helps visualise skewed intensity surfaces
) +
gg(mexdolphin$ppoly, linewidth = 1, alpha = 0, colour = "white") +
gg(mexdolphin$samplers, colour = "grey80", linewidth = 0.4) +
gg(mexdolphin$points, size = 0.8, colour = "white", alpha = 0.7) +
coord_sf(datum = fm_crs(mexdolphin$ppoly)) +
labs(
title = "Predicted dolphin density surface",
subtitle = "exp(W(s) + μ) — the LGCP intensity (mean posterior)",
x = NULL, y = NULL
) +
theme_minimal(base_size = 12) +
theme(legend.position = "right")
Predicted mean dolphin density (groups per nmi²). The brighter the colour, the higher the predicted density.
What this map shows:
The colour scale represents predicted group density (groups per square km). Hot spots (bright colours) are where the model thinks dolphins are most likely to be found. The spatial field \(W(s)\) has pulled density towards areas where multiple detections cluster, while transect coverage patterns influence how much information we have in each area.
10. The detection function
Let’s extract and plot the posterior for the detection function \(g(y|\sigma)\).
# Create a data frame of distances at which to evaluate g(y)
distdf <- data.frame(distance = seq(0, W, length.out = 100))
# predict() evaluates hn(distance, sigma) at each distance value.
# inlabru automatically propagates posterior uncertainty in sigma
# through to uncertainty in g(y). So we get a posterior for the
# entire detection function curve.
dfun <- predict(fit, distdf, ~ hn(distance, sigma))
# Plot
plot(dfun) +
labs(
title = "Fitted detection function: Half-normal",
subtitle = "g(y) = exp(-0.5 (y/σ)²) | shaded = 95% credible interval",
x = "Perpendicular distance y (km)",
y = "Detection probability g(y)"
) +
theme_minimal(base_size = 12)
Posterior detection function. Shaded band = 95% credible interval.
11. Abundance estimation
Expected number of groups
The expected total number of groups in the study area is:
\[ E[N] = \int_A \lambda(s)\, ds \approx \sum_i w_i \exp(W(s_i) + \mu) \]
where \(w_i\) are the quadrature weights from the mesh integration.
# fm_int() creates integration points over the study area polygon using the mesh.
# These give us the quadrature nodes and weights for approximating the area integral.
predpts <- fm_int(mexdolphin$mesh, mexdolphin$ppoly)
# predict() with the formula ~ sum(weight * exp(...)) sums up intensity x area
# across all integration points — giving the posterior for E[N].
Lambda <- predict(
fit,
predpts,
~ sum(weight * exp(mySPDE + Intercept))
)
Lambda## Posterior mean N̂ : 244## Posterior SD : 47## 95% CI : [ 158 , 338 ]Note: \(\Lambda\) here is the expected number of groups given the random field. The actual number \(N\) follows a Poisson distribution with this mean. The posterior for \(N\) itself is more spread out (see below).
13. Exercise: the hazard-rate detection function
Why bother?
The half-normal forces detection to drop smoothly from 1 at \(y=0\). The hazard-rate model allows for a “shoulder” — near-certain detection close to the transect, then a sharper drop at larger distances. We fix b to 1 for computational demands, but can estimate it as well!:
\[ g(y) = 1 - \exp\!\left(-\left(\frac{y}{\sigma}\right)^{-1}\right) \]
This is often more realistic for aerial or vessel surveys where animals very close to the platform are almost certainly seen.
Your task: fit the hazard-rate LGCP, produce the same set of plots, and compare the results to the half-normal model.
Step 1 — Define the hazard-rate function
Step 2 — Update the formula
# Replace hn() with hr() in the formula.
# Everything else (components, priors, domain) stays the same!
#
# YOUR TURN: update the formula.
formula_hr <- geometry + distance ~ mySPDE +
log(hr(distance, sigma)) +
Intercept + log(2)Hint: The
cmpobject (components) can be reused exactly as-is. Only the formula changes. See the inlabru vignette for reference: https://inlabru-org.github.io/inlabru/articles/2d_lgcp_distancesampling.html
Step 3 — Fit the model
# Fit the hazard-rate model using lgcp().
# YOUR TURN: copy the lgcp() call from above and change `formula = form`
# to `formula = formula_hr`.
fit_hr <- lgcp(
components = cmp,
mexdolphin$points,
samplers = mexdolphin$samplers,
domain = list(
geometry = mesh,
distance = fm_mesh_1d(seq(0, W, length.out = 30))
),
formula = formula_hr
)Step 4 — SPDE parameter posteriors
# YOUR TURN: extract and plot the range and log-variance posteriors
# for fit_hr. Compare the range to what you got for the half-normal model.
# Do the spatial correlation structures look similar?
spde.range_hr <- spde.posterior(fit_hr, "mySPDE", what = "range")
plot(spde.range_hr) +
labs(title = "Hazard-rate: posterior spatial range", x = "Range (nmi)", y = "Density") +
theme_minimal(base_size = 12)
Hazard-rate model: SPDE posterior for spatial range.
Step 5 — Predicted spatial intensity
# YOUR TURN: predict the spatial intensity for fit_hr on the same pixel grid.
# Use the same pxl object from above.
pr.int_hr <- predict(fit_hr, pxl, ~ exp(mySPDE + Intercept))
ggplot() +
gg(pr.int_hr, geom = "tile", aes(fill = mean)) +
scale_fill_viridis_c(
name = "Intensity\n(groups/km²)",
option = "plasma",
trans = "sqrt"
) +
gg(mexdolphin$ppoly, linewidth = 1, alpha = 0, colour = "white") +
gg(mexdolphin$samplers, colour = "grey80", linewidth = 0.4) +
gg(mexdolphin$points, size = 0.8, colour = "white", alpha = 0.7) +
coord_sf(datum = fm_crs(mexdolphin$ppoly)) +
labs(
title = "Hazard-rate: predicted density surface",
subtitle = "exp(W(s) + μ) — mean posterior",
x = NULL, y = NULL
) +
theme_minimal(base_size = 12)
Hazard-rate model: predicted dolphin density surface.
Step 6 — Detection function
# YOUR TURN: predict the detection function for fit_hr.
# Use the same distdf object (distances from 0 to W).
# Notice the shoulder near y=0 compared to the half-normal!
dfun_hr <- predict(fit_hr, distdf, ~ hr(distance, sigma))
plot(dfun_hr) +
labs(
title = "Fitted detection function: Hazard-rate",
subtitle = "g(y) = 1 - exp(-(y/σ)⁻¹) | shaded = 95% credible interval",
x = "Perpendicular distance y (nmi)",
y = "Detection probability g(y)"
) +
theme_minimal(base_size = 12)
Hazard-rate detection function with 95% credible interval.
Step 7 — Abundance N
# YOUR TURN: estimate E[N] for the hazard-rate model,
# then compute the full posterior PMF for N.
# Compare the mean and CI to the half-normal result.
Lambda_hr <- predict(
fit_hr,
predpts,
~ sum(weight * exp(mySPDE + Intercept))
)
cat("Hazard-rate — Posterior mean N̂:", round(Lambda_hr$mean), "\n")## Hazard-rate — Posterior mean N̂: 305## Hazard-rate — 95% CI : [ 164 , 524 ]Abundance comparison
# Summary table comparing the two models
comparison <- data.frame(
Model = c("Half-normal", "Hazard-rate"),
Mean_N = c(round(Lambda$mean), round(Lambda_hr$mean)),
SD_N = c(round(Lambda$sd), round(Lambda_hr$sd)),
CI_lower = c(round(Lambda$q0.025), round(Lambda_hr$q0.025)),
CI_upper = c(round(Lambda$q0.975), round(Lambda_hr$q0.975))
)
comparisonDiscussion questions:
- Which model gives a higher N̂? Why might the hazard-rate lead to a larger estimate?
- The credible intervals are wide. What are the main sources of this uncertainty?
- Would you expect both models to give similar maps even if their N̂ values differ? Why?
14. Further reading and getting help
If you want to go deeper:
inlabru website: https://inlabru-org.github.io/inlabru/ — includes vignettes for 1D and 2D LGCPs, mark models, data integration, and more.
This vignette (the mexdolphin example): https://inlabru-org.github.io/inlabru/articles/2d_lgcp_distancesampling.html
INLA book (Krainski et al., 2019): Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA — available free at https://becarioprecario.bitbucket.io/spde-gitbook/
Lindgren et al. (2011): “An explicit link between Gaussian fields and Gaussian Markov random fields” — the SPDE paper. Journal of the Royal Statistical Society B, 73(4).
Miller et al. (2013): “Spatial models for distance sampling data” — the DSM approach (if you want to compare inlabru to a more classical framework). Methods in Ecology and Evolution, 4(11).
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Europe/London
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dsm_2.3.4 numDeriv_2016.8-1.1 mrds_3.0.1
## [4] mgcv_1.9-1 nlme_3.1-166 INLA_24.06.27
## [7] sp_2.2-0 Matrix_1.7-1 inlabru_2.12.0.9022
## [10] fmesher_0.7.0 dplyr_1.2.1 ggplot2_4.0.3
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.49 bslib_0.8.0
## [4] lattice_0.22-6 vctrs_0.7.3 tools_4.4.2
## [7] Rdpack_2.6.1 generics_0.1.4 stats4_4.4.2
## [10] parallel_4.4.2 tibble_3.3.0 proxy_0.4-27
## [13] pkgconfig_2.0.3 KernSmooth_2.23-24 RColorBrewer_1.1-3
## [16] S7_0.2.1 optimx_2025-4.9 lifecycle_1.0.5
## [19] truncnorm_1.0-9 compiler_4.4.2 farver_2.1.2
## [22] MatrixModels_0.5-4 mnormt_2.1.1 statmod_1.5.0
## [25] codetools_0.2-20 htmltools_0.5.8.1 class_7.3-22
## [28] sass_0.4.9 Rsolnp_2.0.1 yaml_2.3.10
## [31] pracma_2.4.4 pillar_1.10.2 nloptr_2.1.1
## [34] jquerylib_0.1.4 classInt_0.4-11 cachem_1.1.0
## [37] parallelly_1.41.0 Deriv_4.1.6 tidyselect_1.2.1
## [40] digest_0.6.37 future_1.34.0 sf_1.0-21
## [43] listenv_0.9.1 labeling_0.4.3 splines_4.4.2
## [46] fastmap_1.2.0 grid_4.4.2 cli_3.6.5
## [49] magrittr_2.0.3 future.apply_1.11.3 e1071_1.7-16
## [52] sn_2.1.1 withr_3.0.2 scales_1.4.0
## [55] rmarkdown_2.29 globals_0.16.3 evaluate_1.0.1
## [58] knitr_1.49 rbibutils_2.3 viridisLite_0.4.2
## [61] rlang_1.2.0 Rcpp_1.0.14 glue_1.8.0
## [64] DBI_1.2.3 rstudioapi_0.17.1 jsonlite_1.8.9
## [67] plyr_1.8.9 R6_2.6.1 units_0.8-7