Coordinate-spatial structured effects • drmTMB

Use spatial() when named sites have coordinates and nearby sites may have similar location deviations after fixed effects have been included. The fitted coordinate route uses a table supplied as coords = coords. A separate, fixed-kappa mesh/SPDE route is now available for one univariate Gaussian location intercept: its observation field is A_st %*% omega, not the dense coordinate route or a site-to-mesh-node lookup. If distance between sites is not what couples your observations, the structural-dependence overview compares this route against the relatedness- and tree-based ones.

What is fitted today

Question	Syntax	Status
Does one Gaussian response have smooth site-level location deviations?	`spatial(1 \| site, coords = coords)` in `mu`	Fitted first coordinate-spatial intercept slice.
Does one Gaussian response have a mesh/SPDE location field on projected coordinates?	`bf(y ~ spatial(1 \| site, mesh = mesh), sigma ~ 1)`	Fixed-kappa intercept at `point_fit_recovery` for the exact tested fixed-domain `n = 128, 256` designs. The retained `n = 64` rung failed; intervals, coverage, and range remain unclaimed.
Does one predictor have a spatially varying slope?	`spatial(1 + depth \| site, coords = coords)` in `mu`	Fitted one numeric-slope slice. The intercept and slope fields are independent and have separate SDs.
Do two Gaussian response means share a coordinate-spatial correlation?	matching `spatial(1 \| p \| site, coords = coords)` terms in `mu1` and `mu2`	Fitted q=2 location-location slice. `corpairs(level = "spatial")` reports the latent spatial row separately from residual `rho12`.
Do spatial location and scale deviations covary across two Gaussian responses?	matching `spatial(1 \| p \| site, coords = coords)` terms in `mu1`, `mu2`, `sigma1`, and `sigma2`	Fitted constant q=4 location-scale slice. Six latent spatial rows are reported through `corpairs(level = "spatial")`; q=4 correlations are derived and unavailable for intervals.

Start with the smallest useful model

For one response, start with a coordinate-spatial location intercept:

fit_spatial <- drmTMB(
  y ~ treatment + spatial(1 | site, coords = coords),
  data = dat,
  family = gaussian()
)

Geographic coordinates and the fixed-kappa mesh route

Longitude and latitude are not model coordinates: decimal degrees are not a metric distance system. Choose a projected CRS appropriate for the study area, transform explicitly, and then make the mesh. spatial_coords() never chooses a UTM zone for you. kappa has inverse projected-coordinate units and remains fixed configuration in this first slice; it is not a fitted range parameter.

mesh_dat <- data.frame(
  y = c(1.1, 1.7, 2.4, 2.0),
  longitude = c(-123.10, -123.05, -123.00, -123.07),
  latitude = c(49.20, 49.23, 49.21, 49.25),
  site = letters[1:4]
)
coords_xy <- spatial_coords(mesh_dat, longitude, latitude, crs_out = "EPSG:32610")
mesh <- make_mesh(coords_xy, kappa = 1 / 10000)
fit_mesh <- drmTMB(
  bf(y ~ spatial(1 | site, mesh = mesh), sigma ~ 1),
  data = mesh_dat, family = gaussian(),
  control = drm_control(se = FALSE)
)
c(
  vertices = ncol(mesh$A_st),
  observations = nrow(mesh$A_st),
  max_projection_row_error = max(abs(Matrix::rowSums(mesh$A_st) - 1))
)
#>                 vertices             observations max_projection_row_error 
#>                       27                        4                        0
mesh_parameters <- summary(fit_mesh)$parameters
mesh_parameters[mesh_parameters$parm == "sd:mu:spatial(1 | site)", ]
#>                                component dpar              term     estimate
#> sd:mu:spatial(1 | site) random-effect-sd   mu spatial(1 | site) 3.173802e-10
#>                         std_error minimum maximum    scale
#> sd:mu:spatial(1 | site)        NA      NA      NA response
#>                                            parm profile_ready
#> sd:mu:spatial(1 | site) sd:mu:spatial(1 | site)         FALSE
#>                                                   profile_note
#> sd:mu:spatial(1 | site) mesh_field_scale_intervals_unvalidated
ranef(fit_mesh, "spatial_mu")$projected
#>             1             2             3             4 
#> -1.068156e-12  1.366843e-13  1.266405e-12 -7.435916e-14
profile_targets(fit_mesh)[, c("parm", "profile_ready", "profile_note")]
#>                      parm profile_ready                           profile_note
#> 1    fixef:mu:(Intercept)          TRUE                                  ready
#> 2 fixef:sigma:(Intercept)          TRUE                                  ready
#> 3                   sigma          TRUE                                  ready
#> 4 sd:mu:spatial(1 | site)         FALSE mesh_field_scale_intervals_unvalidated

The model is y = X beta + A_st omega + epsilon, with omega ~ Normal(0, s^2 Q(kappa)^(-1)) and Q(kappa) = kappa^4 C0 + 2 kappa^2 C1 + C2. Thus the sd:mu:spatial(1 | site) row of summary(fit_mesh)$parameters reports the fitted GMRF field scale s; after projection, the marginal SD at an observation generally varies with its row of A_st. Do not treat it as a single uniform marginal field SD. The existing coords = coords route remains a distinct dense covariance model and is unchanged.

ranef(fit_mesh, "spatial_mu")$latent contains conditional values at mesh vertices. Use $projected for the corresponding observation-level conditional field values. The displayed profile_targets() row is deliberately not ready: this local-fit slice does not claim a field-scale interval.

The raw GMRF field scale has point-recovery evidence for the exact tested fixed-domain n = 128 and n = 256 designs. The retained n = 64 rung failed, so this is not a universal n >= 128 guarantee. The mesh slice rejects raw geographic degrees, mesh slopes or labels, sigma/shape mesh effects, non-Gaussian or bivariate models, mesh-plus-coords formulas, extrapolation beyond the mesh, range estimation, anisotropy, barriers, replicated fields, and spatiotemporal fields.

Add one numeric slope only when the scientific question is about spatial variation in that slope:

fit_spatial_slope <- drmTMB(
  y ~ treatment + depth + spatial(1 + depth | site, coords = coords),
  data = dat,
  family = gaussian()
)

For the exact Arc 1a REML route, keep sigma ~ 1, use an unlabelled intercept or independent intercept-plus-one-numeric-slope shape, and set REML = TRUE:

fit_spatial_reml <- drmTMB(
  bf(
    y ~ depth + spatial(1 + depth | site, coords = coords),
    sigma ~ 1
  ),
  data = dat,
  family = gaussian(),
  REML = TRUE
)

The multi-seed campaign used the coordinate representation shown here, with n_each = 20 and exactly M = {8, 16, 32} sites. This is not a continuous minimum-sample-size claim, and it does not admit estimated range, labelled, slope-only, multiple-slope, scale-side, other bivariate, or non-Gaussian REML routes. The exact bivariate exception is the matched labelled location-intercept cell below.

For two response means, use matching labelled terms and read the latent spatial correlation with corpairs():

fit_spatial_q2 <- drmTMB(
  mu1 = trait1 ~ treatment +
    spatial(1 | p | site, coords = coords),
  mu2 = trait2 ~ treatment +
    spatial(1 | p | site, coords = coords),
  data = dat,
  family = biv_gaussian()
)
corpairs(fit_spatial_q2, level = "spatial")
rho12(fit_spatial_q2)

corpairs() reports the fitted latent coordinate-spatial correlation among site-level location deviations. rho12() reports the residual correlation between paired responses after fixed effects and random effects have been included.

That exact q2 location-intercept model also admits native REML when the three residual parameters are intercept-only:

fit_spatial_q2_reml <- drmTMB(
  bf(
    mu1 = trait1 ~ treatment +
      spatial(1 | p | site, coords = coords),
    mu2 = trait2 ~ treatment +
      spatial(1 | p | site, coords = coords),
    sigma1 = ~ 1,
    sigma2 = ~ 1,
    rho12 = ~ 1
  ),
  data = dat,
  family = biv_gaussian(),
  REML = TRUE
)

Here the coordinates define a fixed spatial covariance matrix. This cell has dense-oracle and retained-denominator point-recovery evidence only; it does not authorize interval, coverage, range-estimation, slope, scale-side, or q4 claims.

For a constant location-scale spatial block, use the same labelled spatial() term in all four bivariate Gaussian endpoints:

fit_spatial_q4 <- drmTMB(
  mu1 = trait1 ~ treatment +
    spatial(1 | p | site, coords = coords),
  mu2 = trait2 ~ treatment +
    spatial(1 | p | site, coords = coords),
  sigma1 = ~ treatment +
    spatial(1 | p | site, coords = coords),
  sigma2 = ~ treatment +
    spatial(1 | p | site, coords = coords),
  rho12 = ~ 1,
  data = dat,
  family = biv_gaussian()
)
corpairs(fit_spatial_q4, level = "spatial")

This q=4 route estimates four coordinate-spatial endpoint SDs and six latent correlations: one location-location, four location-scale, and one scale-scale row. It is still a constant intercept block, not a spatial slope or predictor-dependent spatial correlation model.

What to inspect

After fitting, inspect the spatial layer before interpreting it:

Output	Use
`check_drm(fit)`	Confirm the spatial layer was recognized. Mesh fits report vertex count, fixed `kappa`, the exact tested point-recovery designs, and the remaining interval/range boundary.
`summary(fit)$parameters`	Read fitted spatial location SDs, including separate intercept and slope SDs when a one-slope model is used.
`ranef(fit, "spatial_mu")`	For `coords =`, inspect conditional site deviations. For a mesh, use `$projected` for observation-level field values and `$latent` for mesh vertices.
`summary(fit)$covariance`	Check how spatial SDs and q=2 or q=4 spatial correlations are reported beside other covariance layers.
`profile_targets(fit)`	See which spatial SD or constant q=2 correlation targets can be profiled directly, and which q=4 correlation rows are derived-unavailable for intervals.
`corpairs(fit, level = "spatial")`	Read fitted constant spatial correlation rows.

The one-slope route is deliberately narrow. The formula term spatial(1 + depth | site, coords = coords) fits an intercept field and one numeric slope field with the same coordinate precision and separate SDs. It does not estimate an intercept-slope correlation. Spatial sigma is supported through a separate route – a standalone sigma ~ spatial(1 | site) or sigma ~ spatial(1 + depth | site) field, or the matched location-scale block – which fits at recovery grade (trust the point estimate, not the interval); this mu one-slope term is not that route.

Rendered checks

The small example below is only a guide to the output grain. The coordinate surface is fitted in the location predictor mu; raw response values remain on the response scale and should not be plotted as if they were spatial SDs or correlations.

spatial_example <- simulate_spatial_guide_data()
spatial_dat <- spatial_example$data
coords <- spatial_example$coords

fit_spatial <- drmTMB(
  drm_formula(
    y ~ depth + temp + spatial(1 | site, coords = coords),
    sigma ~ depth
  ),
  family = gaussian(),
  data = spatial_dat
)

fit_spatial_slope <- drmTMB(
  drm_formula(
    y ~ depth + temp + spatial(1 + depth | site, coords = coords),
    sigma ~ depth
  ),
  family = gaussian(),
  data = spatial_dat
)

if (requireNamespace("ggplot2", quietly = TRUE)) {
  spatial_effect <- ranef(fit_spatial, "spatial_mu")$terms[[1]]
  spatial_field <- data.frame(
    site = names(spatial_effect),
    fitted_spatial_deviation = unname(spatial_effect),
    coords[names(spatial_effect), , drop = FALSE],
    row.names = NULL
  )
  field_limit <- max(abs(spatial_field$fitted_spatial_deviation))
  if (!is.finite(field_limit) || field_limit == 0) field_limit <- 1

  ggplot2::ggplot(
    spatial_field,
    ggplot2::aes(
      x = x,
      y = y,
      fill = fitted_spatial_deviation
    )
  ) +
    ggplot2::geom_hline(yintercept = 0, colour = "grey90", linewidth = 0.4) +
    ggplot2::geom_vline(xintercept = 0, colour = "grey90", linewidth = 0.4) +
    ggplot2::geom_point(
      shape = 21,
      size = 7,
      colour = "grey20",
      stroke = 0.35
    ) +
    ggplot2::scale_fill_gradient2(
      low = "#D55E00",
      mid = "white",
      high = "#009E73",
      midpoint = 0,
      limits = c(-field_limit, field_limit),
      name = "Fitted\nspatial deviation"
    ) +
    ggplot2::coord_equal() +
    spatial_guide_theme() +
    ggplot2::labs(
      title = "Fitted spatial location field",
      subtitle = "Conditional fitted deviations; uncertainty not shown",
      x = "Coordinate x",
      y = "Coordinate y"
    )
}

Map of twelve simulated sampled sites. Each point is positioned at the site coordinates and coloured by the fitted conditional spatial location deviation; positive deviations are teal and negative deviations are orange.

Simulated example of coordinate-spatial fitted site deviations from ranef(fit_spatial, "spatial_mu"). Points are conditional location-effect estimates; uncertainty is not shown.

Spatial intercept and slope SDs have different units, so placing them on one quantitative axis would imply a comparison that is not meaningful. The compact display below reports each estimate in its own unit instead.

spatial_sd <- data.frame(
  Component = c("Spatial intercept SD", "Spatial depth-slope SD"),
  Estimate = formatC(
    summary(fit_spatial_slope)$parameters[
      match(
        c("sd:mu:spatial(1 | site)", "sd:mu:spatial(0 + depth | site)"),
        summary(fit_spatial_slope)$parameters$parm
      ),
      "estimate"
    ],
    digits = 4,
    format = "g"
  ),
  Unit = c("Response units", "Response units per depth unit"),
  Status = c(
    "Point estimate; interval not validated",
    "Near-zero boundary; interval not validated"
  ),
  check.names = FALSE
)

knitr::kable(spatial_sd, align = c("l", "r", "l", "l"))

Component	Estimate	Unit	Status
Spatial intercept SD	0.3808	Response units	Point estimate; interval not validated
Spatial depth-slope SD	2.667e-06	Response units per depth unit	Near-zero boundary; interval not validated

For the exact fixed-kappa bivariate Gaussian location model, the calibrated M rung has 36 sites with three complete response pairs per site and the baseline ring geometry. Native REML for this cell requires unit weights, intercept-only sigma1, sigma2, and rho12, no known meta_V() covariance, and no additional ordinary random effect, direct-SD formula, or corpair() regression.

spatial_q2_example <- simulate_spatial_q2_guide_data()
spatial_q2_dat <- spatial_q2_example$data
spatial_q2_coords <- spatial_q2_example$coords

fit_spatial_q2_example <- drmTMB(
  drm_formula(
    mu1 = y1 ~ x1 + spatial(1 | p | site, coords = spatial_q2_coords),
    mu2 = y2 ~ x2 + spatial(1 | p | site, coords = spatial_q2_coords),
    sigma1 = ~ 1,
    sigma2 = ~ 1,
    rho12 = ~ 1
  ),
  family = c(gaussian(), gaussian()),
  data = spatial_q2_dat,
  REML = TRUE
)

The prospective campaign retained every attempted dataset. At M, all-attempt coverage was 0.938, 0.932, and 0.938 for the first spatial SD, second spatial SD, and latent spatial correlation; finite-profile rates were 1.000, 1.000, and 0.986. The higher H rung (36 sites x 8 observations) also passed jointly. The smaller L rung (12 x 3) failed and is not part of the interval claim.

spatial_q2_targets <- c(
  "sd:mu:mu1:spatial(1 | p | site)",
  "sd:mu:mu2:spatial(1 | p | site)",
  "cor:spatial:cor(mu1:(Intercept),mu2:(Intercept) | p | site)"
)

spatial_q2_profile <- stats::confint(
  fit_spatial_q2_example,
  parm = spatial_q2_targets,
  method = "profile",
  profile_engine = "endpoint"
)
spatial_q2_target_table <- profile_targets(fit_spatial_q2_example)
spatial_q2_estimate <- spatial_q2_target_table$estimate[
  match(spatial_q2_targets, spatial_q2_target_table$parm)
]
spatial_q2_eye <- data.frame(
  target = factor(
    c(
      "Spatial SD: response 1",
      "Spatial SD: response 2",
      "Latent spatial correlation"
    ),
    levels = c(
      "Spatial SD: response 1",
      "Spatial SD: response 2",
      "Latent spatial correlation"
    )
  ),
  estimate = unname(spatial_q2_estimate),
  lower = spatial_q2_profile$lower,
  upper = spatial_q2_profile$upper
)
spatial_q2_eye_region <- do.call(
  rbind,
  lapply(seq_len(nrow(spatial_q2_eye)), function(i) {
    eye_x <- seq(
      spatial_q2_eye$lower[i],
      spatial_q2_eye$upper[i],
      length.out = 101
    )
    left_width <- max(
      spatial_q2_eye$estimate[i] - spatial_q2_eye$lower[i],
      .Machine$double.eps
    )
    right_width <- max(
      spatial_q2_eye$upper[i] - spatial_q2_eye$estimate[i],
      .Machine$double.eps
    )
    taper <- ifelse(
      eye_x <= spatial_q2_eye$estimate[i],
      (eye_x - spatial_q2_eye$lower[i]) / left_width,
      (spatial_q2_eye$upper[i] - eye_x) / right_width
    )
    half_height <- 0.10 * sqrt(pmax(taper, 0))
    data.frame(
      target = spatial_q2_eye$target[i],
      eye_x = c(eye_x, rev(eye_x)),
      eye_y = c(half_height, rev(-half_height))
    )
  })
)

The Confidence Eye treats each interval as a small pale tapered region and marks the estimate with a hollow circle. The eye’s horizontal span is the interval; there is no separate interval bar. Separate facet scales keep standard deviations and correlation on their own units.

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(spatial_q2_eye) +
    ggplot2::geom_vline(
      data = data.frame(
        target = factor(
          "Latent spatial correlation",
          levels = levels(spatial_q2_eye$target)
        ),
        zero = 0
      ),
      ggplot2::aes(xintercept = zero),
      inherit.aes = FALSE,
      linetype = "dotted",
      linewidth = 0.5,
      colour = "grey55"
    ) +
    ggplot2::geom_polygon(
      data = spatial_q2_eye_region,
      ggplot2::aes(
        x = eye_x,
        y = eye_y,
        group = target,
        fill = target
      ),
      inherit.aes = FALSE,
      alpha = 0.24,
      colour = NA
    ) +
    ggplot2::geom_point(
      ggplot2::aes(
        x = estimate,
        y = 0,
        colour = target
      ),
      shape = 21,
      fill = "white",
      size = 4.2,
      stroke = 1.2
    ) +
    ggplot2::facet_wrap(~target, scales = "free_x", nrow = 1) +
    ggplot2::scale_x_continuous(
      expand = ggplot2::expansion(mult = c(0.20, 0.20))
    ) +
    ggplot2::scale_y_continuous(
      NULL,
      breaks = NULL,
      limits = c(-0.22, 0.32),
      expand = c(0, 0)
    ) +
    ggplot2::scale_fill_manual(
      values = c(
        "Spatial SD: response 1" = "#0072B2",
        "Spatial SD: response 2" = "#D55E00",
        "Latent spatial correlation" = "#009E73"
      ),
      guide = "none"
    ) +
    ggplot2::scale_colour_manual(
      values = c(
        "Spatial SD: response 1" = "#0072B2",
        "Spatial SD: response 2" = "#D55E00",
        "Latent spatial correlation" = "#009E73"
      ),
      guide = "none"
    ) +
    ggplot2::labs(
      x = "Target value (facet-specific scale)",
      title = "Profile uncertainty for the calibrated spatial q2 target set",
      subtitle = "Each eye is a 95% endpoint profile interval; hollow circle marks the estimate"
    ) +
    ggplot2::theme_minimal(base_size = 12) +
    ggplot2::theme(
      panel.grid.major.x = ggplot2::element_line(
        colour = "grey90",
        linewidth = 0.35
      ),
      panel.grid.major.y = ggplot2::element_blank(),
      panel.grid.minor = ggplot2::element_blank(),
      panel.spacing.x = grid::unit(1.3, "lines"),
      strip.text = ggplot2::element_text(face = "bold"),
      plot.title.position = "plot"
    )
}

Confidence Eye for the three direct fixed-kappa Gaussian q2 spatial targets at the tested M rung (36 sites x 3 observations, baseline ring geometry). Each coloured pale eye spans a 95% endpoint profile-likelihood interval; the larger hollow circle is the point estimate. Calibration passed jointly at the exact M and H rungs and failed at L.

This result supports inference_ready_with_caveats only for the exact tested M/H fixed-kappa ring configurations. It does not establish mesh intervals, estimated range, spatial slopes, q4+, non-Gaussian spatial models, spatial scale models, derived observed correlations, geometry robustness, or the supported tier.

Boundaries

The following spatial routes remain deferred:

multiple spatial slopes and spatial slope correlations;
partial spatial terms in sigma, plus spatial terms in nu, zero-inflation, or rho12;
direct spatial SD surfaces;
predictor-dependent spatial corpair() regressions;
simultaneous phylo() plus spatial() layers in the same formula;
non-Gaussian spatial structured effects outside the exact ordinary Poisson/NB2 q1 spatial mu intercept-plus-one-slope, recovery-grade NB2 q1 spatial sigma, Student-t spatial mu, Poisson spatial zi, fixed-zi Poisson spatial mu, and fixed-zi NB2 spatial mu gates.

Use the structural-dependence overview when you are choosing among animal(), phylo(), spatial(), and relmat(). Use the detailed structural-dependence tutorial when you need the current worked examples, equations, and broader parity ladder.