## ----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
  )
}

## -----------------------------------------------------------------------------
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")]

## ----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
)

## ----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-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"))

## ----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"
  )
)

## ----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))
    )
  })
)

## ----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"
    )
}