---
title: "Coordinate-spatial structured effects"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Coordinate-spatial structured effects}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.1,
  dpi = 144
)
library(drmTMB)

spatial_guide_theme <- function() {
  ggplot2::theme_minimal(base_size = 11) +
    ggplot2::theme(
      panel.grid.minor = ggplot2::element_blank(),
      axis.title = ggplot2::element_text(colour = "grey15"),
      axis.text = ggplot2::element_text(colour = "grey25"),
      plot.title = ggplot2::element_text(
        face = "bold",
        colour = "grey10",
        margin = ggplot2::margin(b = 4)
      ),
      plot.subtitle = ggplot2::element_text(
        colour = "grey30",
        margin = ggplot2::margin(b = 8)
      ),
      legend.position = "bottom"
    )
}

spatial_eye_theme <- function() {
  spatial_guide_theme() +
    ggplot2::theme(
      panel.grid.major.y = ggplot2::element_blank(),
      legend.position = "none"
    )
}

simulate_spatial_guide_data <- function(seed = 20260534) {
  set.seed(seed)
  sites <- paste0("site_", seq_len(12))
  theta <- seq(0, 1.8 * pi, length.out = length(sites))
  coords <- data.frame(
    x = cos(theta) + seq_along(sites) / 30,
    y = sin(theta)
  )
  rownames(coords) <- sites

  distance <- as.matrix(dist(coords))
  range <- stats::median(distance[distance > 0])
  spatial_cov <- exp(-distance / range)
  diag(spatial_cov) <- diag(spatial_cov) + 1e-6
  spatial_intercept <- as.vector(
    t(chol(spatial_cov)) %*% rnorm(length(sites), sd = 0.45)
  )
  spatial_slope <- as.vector(
    t(chol(spatial_cov)) %*% rnorm(length(sites), sd = 0.20)
  )
  names(spatial_intercept) <- sites
  names(spatial_slope) <- sites

  site <- rep(sites, each = 6)
  temp <- rnorm(length(site))
  depth <- rnorm(length(site))
  sigma <- exp(-1.25 + 0.15 * depth)
  y <- 0.4 + 0.35 * temp + 0.20 * depth +
    spatial_intercept[site] + spatial_slope[site] * depth +
    rnorm(length(site), sd = sigma)

  list(
    data = data.frame(y = y, temp = temp, depth = depth, site = site),
    coords = coords
  )
}

simulate_spatial_q2_guide_data <- function(seed = 260805001L) {
  set.seed(seed)
  sites <- paste0("site_", seq_len(36))
  theta <- seq(0, 1.5 * pi, length.out = length(sites))
  coords <- data.frame(
    x = cos(theta) + seq_along(sites) / (3 * length(sites)),
    y = sin(theta)
  )
  rownames(coords) <- sites

  distance <- as.matrix(dist(coords))
  range <- stats::median(distance[distance > 0])
  spatial_cov <- exp(-distance / range)
  diag(spatial_cov) <- diag(spatial_cov) + 1e-6

  z1 <- rnorm(length(sites))
  z2 <- 0.45 * z1 + sqrt(1 - 0.45^2) * rnorm(length(sites))
  u1 <- as.vector(t(chol(spatial_cov)) %*% z1) * 0.55
  u2 <- as.vector(t(chol(spatial_cov)) %*% z2) * 0.55
  names(u1) <- sites
  names(u2) <- sites

  site <- rep(sites, each = 3)
  x1 <- rnorm(length(site))
  x2 <- rnorm(length(site))
  e1 <- rnorm(length(site))
  e2 <- -0.10 * e1 + sqrt(1 - (-0.10)^2) * rnorm(length(site))

  list(
    data = data.frame(
      site = site,
      x1 = x1,
      x2 = x2,
      y1 = 0.35 + 0.25 * x1 + u1[site] + 0.18 * e1,
      y2 = -0.20 - 0.30 * x2 + u2[site] + 0.20 * e2
    ),
    coords = coords
  )
}
```

<style>
@media (max-width: 575.98px) {
  main h1 {
    hyphens: none;
    overflow-wrap: normal;
    word-break: normal;
  }

  main table {
    display: block;
    min-width: 720px;
    max-width: 100%;
    overflow-x: auto;
  }
}
</style>

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](structural-dependence.html) 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:

```r
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.

```{r}
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))
)
mesh_parameters <- summary(fit_mesh)$parameters
mesh_parameters[mesh_parameters$parm == "sd:mu:spatial(1 | site)", ]
ranef(fit_mesh, "spatial_mu")$projected
profile_targets(fit_mesh)[, c("parm", "profile_ready", "profile_note")]
```

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:

```r
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`:

```r
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()`:

```r
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:

```r
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:

```r
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.

```{r spatial-guide-fit}
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
)
```

```{r spatial-site-field-figure, fig.width = 6.6, fig.height = 4.4, fig.cap = "Simulated example of coordinate-spatial fitted site deviations from `ranef(fit_spatial, \"spatial_mu\")`. Points are conditional location-effect estimates; uncertainty is not shown.", fig.alt = "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."}
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"
    )
}
```

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.

```{r spatial-sd-summary}
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"))
```

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.

```{r spatial-q2-correlation-fit}
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 = biv_gaussian(),
  data = spatial_q2_dat,
  REML = TRUE,
  control = drm_control(
    optimizer = list(eval.max = 1000L, iter.max = 1000L),
    fallback_optimizer = "BFGS"
  )
)
```

The fixed seed is the first retained M-rung baseline-ring smoke dataset. It
keeps this tutorial fit reproducible without re-running the coverage campaign.

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.

```{r spatial-q2-profile-intervals}
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.

```{r spatial-q2-confidence-eye, fig.width = 9.2, fig.height = 3.1, fig.cap = "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.", fig.alt = "Three side-by-side facets show coloured tapered confidence eyes for two spatial standard deviations and one latent spatial correlation. Each eye's horizontal width is its 95 percent profile interval, and a large hollow circle marks the estimate. Each facet uses its own horizontal scale."}
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"
    )
}
```

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](structural-dependence.html) when you
are choosing among `animal()`, `phylo()`, `spatial()`, and `relmat()`. Use the
detailed [structural-dependence tutorial](https://itchyshin.github.io/drmTMB/articles/phylogenetic-spatial.html) when you
need the current worked examples, equations, and broader parity ladder.