Exponential spatial range: metafor SPEXP rho vs brms gp lscale mismatch (other estimates match)

Interfaces brms
,
1

Hi all, I am comparing an exponential spatial meta-analytic model fitted in metafor and brms. The overall mean effect and variance components are very similar across packages, but the spatial “range” parameter looks very different at first glance.
I would like to confirm whether this is expected (due to parameterisation and weak identifiability) or whether I am missing something in the brms specification.

Data (brief description)

I am using a published global dataset on plant responses to intensified fire regimes (Grau-Andrés et al. 2024; https://doi.org/10.1111/geb.13858). The analysis dataset contains k = 2361 effect sizes (Hedges’ d), their corresponding sampling variances (var_Hedges), and geographic coordinates (latitude/longitude). I remove rows with missing longitude and treat each row as a unique effect size (effect_id).
The goal is to fit a spatially correlated structure using an exponential correlation function.

Spatial preprocessing (Lat/Long to km) and model overview

I convert WGS84 lat/long (EPSG:4326) to planar coordinates (EPSG:3857) and then to kilometers:

dat <- read.csv("Roger_etal_2024.csv") |>
  filter(!is.na(longitude)) |>
  mutate(const = 1,
         effect_id = seq_len(n()))
dat$effect_id <- factor(as.character(dat$effect_id),
                    levels = as.character(sort(unique(dat$effect_id))))
dat_sf  <- st_as_sf(dat, coords = c("longitude", "latitude"), crs = 4326)
dat_prj <- st_transform(dat_sf, crs = 3857)
xy_m <- st_coordinates(dat_prj)
dat$x_km <- xy_m[,1] / 1000
dat$y_km <- xy_m[,2] / 1000

And here are the models used:

## metafor ----
fit_metafor <- rma.mv(
  d_Hedges, var_Hedges,
  random = list(
    ~ 1 | effect_id,
    ~ x_km + y_km | const),
  struct = "SPEXP",
  data = dat,
  sparse = TRUE,
  method = "REML",
  test = "t")

## brms ----
vcv <- diag(dat$var_Hedges)
rownames(vcv) <- colnames(vcv) <- as.character(dat$effect_id)

bf_sp <- bf(
  d_Hedges ~ 1 +
    (1 | gr(effect_id, cov = vcv)) +
    gp(x_km, y_km, cov = "exponential", scale = FALSE))

prior <- default_prior(bf_sp, data = dat, data2 = list(vcv = vcv), family = gaussian())

# fix SD for the known sampling-variance term
prior$prior[5] <- "constant(1)"

fit_brms <- brm(
  bf_sp,
  data = dat,
  data2 = list(vcv = vcv),
  family = gaussian(),
  prior = prior,
  chains = 2, iter = 6000, warmup = 3000)

The table below provides a summary of the model estimates:

Estimate metafor brms
Overall mean -0.33 -0.32
Spatial variance \tau^2 = 1.23 \text{sdgp}^2 = 1.14^2 = 1.30
Effect size variance \sigma^2 = 0.79 \sigma^2 = 0.97^2 = 0.94
Spatial range \rho = 0.174 \text{lscale} = 128 km

What matches and what looks different?

  • Mean effect (Intercept) and residual variance are similar across packages.
  • Spatial variance is similar: metafor \tau^2 \approx 1.23 (metafor) versus \text{sdgp}^2 \approx 1.30 (brms).
  • The range parameter looks different: \rho \approx 0.174 (metafor) versus \text{lscale} \approx 128 km (brms).

My understanding

My understanding is that the exponential kernel is parameterised differently:

  • metafor - SPEXP: \mathrm{cor}(d) = \exp(- d \cdot \rho)
  • brms - exponential GP: \mathrm{cor}(d) = \exp(- d/\text{lscale})

So I expect \rho \approx 1/lscale (with d in km because I use x_km and y_km and set scale = FALSE).
With \text{lscale} \approx 127.84 km, 1/\text{lscale} \approx 0.0078 per km, which seems compatible with the very flat or weakly identified interval for \rho reported by metafor.

Questions

  1. In brms::gp(x_km, y_km, cov = "exponential", scale = FALSE), is the exponential GP correlation parameterised as \mathrm{cor}(d) = \exp(- d/\text{lscale}) with d measured in the same units as x_km and y_km (km in my case)?
    If so, is it appropriate to compare this to metafor’s struct = "SPEXP" parameterisation, which uses \mathrm{cor}(d) = \exp(- d \cdot \rho), through the mapping \rho \approx 1/\text{lscale}?
  2. Considering that metafor reports a very wide 95% CI for \rho (\rho < 0.0174 to 0.8572), is it reasonable to interpret the apparent difference in the point estimates of the range parameter as weak identifiability of the distance scale rather than an implementation mismatch between metafor and brms?

Thank you very much!

Session info

R 4.5.2 (2025-10-31), macOS Sequoia 15.7.3 (aarch64-apple-darwin20),
brms 2.23.0, metafor 4.8-0, sf 1.0-23
(I can paste the full sessionInfo() if needed)