DRModels.jl

Appearance

Coordinate-spatial structured effects ​

Status: Stable (Gaussian mean, coordinate-direct)

Mirrors drmTMB's Coordinate-spatial structured effects. In DRModels.jl today: spatial(1 | site) on the Gaussian mean with site coordinates; an exponential spatial correlation K(ρ) = exp(-d / ρ) whose range ρ is estimated jointly. Closed-form GLS (coordinate-direct; mesh/SPDE is planned).

Not the same model as drmTMB's spatial()

The same spatial(1 | site) term means different things in the two packages. DRModels.jl estimates the range ρ. drmTMB fixes the range at the median positive pairwise distance between sites and estimates only the spatial SD. On the same data and formula the two packages therefore report different log-likelihoods and spatial SDs; the DRModels.jl fit has one more parameter. Compare the two only after accounting for this.

Nearby sites tend to be similar. spatial(1 | site) adds a random intercept with a distance-based correlation built from the site coordinates, u ~ N(0, σ_s² K(ρ)), and estimates the spatial range ρ along with the spatial SD. Pass one coordinate row per site and label the rows with ids. The labels must match the values in the grouping column. A Dict from site labels to coordinate vectors is another option. Unlabelled matrices retain first-seen site ordering and emit a warning; their rows are not sorted by label. Extra labelled sites are dropped, while a missing site label raises an error. If the model also uses animal(), ids labels the pedigree matrix A; use a Dict to label the spatial coordinates in that case.

julia
using DRModels, Random, LinearAlgebra
Random.seed!(8)

G = 50
coords = rand(G, 2) .* 10.0                  # site locations in [0,10]²
Ddist = [sqrt(sum(abs2, coords[k, :] .- coords[l, :])) for k in 1:G, l in 1:G]
K = exp.(-Ddist ./ 2.0) + 1e-8 * I           # true correlation, range = 2
u = 0.9 .* (cholesky(Symmetric(K)).L * randn(G))

m = 4; n = G * m
site = repeat(1:G, inner = m)
x = randn(n)
y = 0.3 .+ 0.5 .* x .+ u[site] .+ 0.4 .* randn(n)

fit = drm(bf(@formula(y ~ x + spatial(1 | site)), @formula(sigma ~ 1)),
          Gaussian(); data = (; y, x, site), coords = coords, ids = collect(1:G))

re_sd(fit)[:site]      # spatial SD (square it for the spatial variance)
1.172682459930717

re_sd(fit)[:site] is the spatial SD σ_s; vc(fit) is the matching random-effect (co)variance summary for the structured block.

julia
exp(coef(fit, :range)[1])     # estimated spatial range ρ
4.03035634483914

This is one simulated realization; it does not establish general parameter recovery or interval calibration. Inspect uncertainty for both the spatial SD and range. For Gaussian mean models, spatial, phylo, animal, and relmat share the closed-form structured GLS engine; the source of the correlation differs (coordinates, tree, pedigree, or supplied matrix).

Check the fitted model before interpreting the spatial SD or range. Start with check_drm and read its convergence, gradient, and Hessian messages; then use Improving convergence if any of those checks warn. The Detailed capabilities & limits page states what this Gaussian spatial route does and does not establish about uncertainty.

Spatial routes depend on the family

This tutorial uses the exact Gaussian coordinate-spatial route, which estimates the range. The non-Gaussian families that take coords, Poisson included, fix the range as drmTMB does: the kernel is exp(-d / m) + 1e-6 I, with m the median positive pairwise distance between sites.

  • Gamma accepts spatial(1 | site) with coords (native kernel) or with a supplied K, not both. A 1-D vector of coordinates is accepted.

  • NB2 accepts spatial(1 | site) with a supplied K. With coords it accepts only a slope-only spatial(0 + x | site) term with its sd(...) formula.

  • Beta accepts spatial(1 | site) with a supplied K when sigma ~ 1 and no SD formula is used. With sd(site, spatial) ~ z it takes coords (native kernel) and refuses K.

  • Binomial refuses a bare spatial(1 | site) on the mean. It accepts spatial with an sd(site, spatial) formula and coords. That route keeps the older fixed range (the mean pairwise distance, 1e-8 ridge) because drmTMB has no Binomial spatial route to match.

Check the family-specific restrictions in Detailed capabilities and limits.

See also ​