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.1726824599322867

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.030356344856249

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. Poisson also accepts coords and estimates the spatial range on its Laplace route. NB2, Gamma, and Beta accept spatial(1 | site) with a supplied covariance matrix K; they do not accept the coordinate workflow shown here. Binomial does not accept the spatial marker. Check the family-specific restrictions in Detailed capabilities and limits.

See also ​