---
title: "Unit-response two-stage continuous space-time model with spatial SVC"
subtitle: "Washington biomass prediction with row-level TCC, FIA-unit variances, and a positive-stage spatial SVC"
author: "Andrew O. Finley"
format:
  html:
    toc: true
    toc-location: right
    toc-depth: 2
    number-sections: true
    code-fold: true
    code-overflow: wrap
    html-math-method: katex
    css: stlmm-vignette.css
execute:
  warning: false
  message: false
bibliography: references.bib
---

```{r setup}
#| code-summary: "Show setup code"
library(stLMM)
library(RhpcBLASctl)
library(sf)
library(tidyverse)
library(knitr)
library(kableExtra)

source("article-utils.R")

set.seed(1)
theme_set(theme_bw(base_size = 12))

article_id <- "wa-unit-two-stage-continuous-space-time-tcc-fiaunit-var-svc"
data_dir <- "wa_data"
```

## Purpose

This article carries forward the two-stage unit-response target used in the previous articles and changes the latent process. The binary stage still models `agb_live_mg_ha > 0`, the positive stage still models square-root positive `agb_live_mg_ha`, row-level TCC remains the auxiliary covariate, and the county-year target remains the grid-supported mean of `agb_live_mg_ha`.

The new feature is the space-time process. The earlier unit-response articles used county-year areal effects, first through simple county and year effects and then through `car_time(county_fips, year)`. Here, the latent process is indexed by continuous projected coordinates and measurement time.

A continuous space-time process also makes prediction less tied to the administrative units used for fitting. Once the model is fit, predictions can be made on any user-defined spatial extent or grid with the required covariates, then summarized over counties, ownerships, ecological regions, management areas, or other supports. This is a central advantage of the unit-response space-time approach used in related forest-biomass applications [@Finley01032011; @MayFinley2025SpatioTemporal].

The positive Gaussian stage carries forward the FIA-unit-specific residual variance from the previous article. That variance refinement is paired here with a positive-stage spatially varying coefficient (SVC), so the TCC-biomass relationship can change smoothly over space rather than being forced to use one global positive-biomass slope.

The fitted model is intentionally rich. It combines a two-stage response, continuous space-time NNGP intercept processes, FIA-unit-specific residual variances, and a spatially varying TCC effect in the positive biomass stage. Before using this kind of model for a final analysis, check that the data identify the added parameters and compare predictive performance using WAIC or, preferably, a cross-validation study.

For plot row $j$ with coordinates $\mathbf{s}_j$, measurement time $t_j$, county $a_j$, FIA unit $u_j$, and standardized TCC $x_j$, the model is

$$
z_j \sim \operatorname{bernoulli}(p_j), \qquad
\operatorname{logit}(p_j) = \gamma_0 + \gamma_1 x_j + w^z(\mathbf{s}_j, t_j),
$$

and, for positive biomass rows,

$$
q_j = \delta_0 + \{\delta_1 + v^+(\mathbf{s}_j)\}x_j + w^+(\mathbf{s}_j, t_j) + r_j,\qquad
r_j \sim \operatorname{normal}(0,\tau^2_{u_j}),
$$

where $z_j = \mathbb{1}(y_j > 0)$ and $q_j = \sqrt{y_j^+}$. The processes $w^z(\mathbf{s}_j, t_j)$ and $w^+(\mathbf{s}_j, t_j)$ are occurrence-stage and positive-stage space-time intercept surfaces, respectively, and $v^+(\mathbf{s}_j)$ is the spatial departure from the global positive-stage TCC slope. Prediction combines draw-matched binary and positive-stage draws on grid row $g$,

$$
\tilde{y}^{(m)}_g = \tilde{z}^{(m)}_g\{\tilde{q}^{(m)}_g\}^2,
$$

then averages over grid rows within county-year.

::: {.stlmm-resource-box .stlmm-resource-source}
<div class="stlmm-resource-title">Source Code</div>

- [Quarto source](wa-unit-two-stage-continuous-space-time-tcc-fiaunit-var-svc.qmd)
- [Extracted R code](wa-unit-two-stage-continuous-space-time-tcc-fiaunit-var-svc.R)
:::

::: {.stlmm-resource-box .stlmm-resource-data}
<div class="stlmm-resource-title">Data Bundle</div>

- [Download Washington biomass data](wa_data.zip)
:::

## Data

```{r read-data}
#| code-summary: "Show data-reading code"
wa_counties <- read_rds(file.path(data_dir, "wa_counties.rds"))
unit_plots <- read_csv(file.path(data_dir, "wa_unit_plots.csv"), show_col_types = FALSE)
direct_estimates <- read_csv(file.path(data_dir, "wa_direct_estimates.csv"), show_col_types = FALSE)
prediction_grid <- read_csv(
  file.path(data_dir, "wa_unit_prediction_grid_5km.csv"),
  show_col_types = FALSE
)

wa_counties <- wa_counties |>
  mutate(
    county_fips = as.character(county_fips),
    fia_unit = factor(fia_unit)
  )

unit_plots <- unit_plots |>
  mutate(
    county_fips = as.character(county_fips),
    county_fips = factor(county_fips),
    fia_unit = factor(fia_unit),
    year = as.integer(year),
    time = as.numeric(time),
    x_km = x / 1000,
    y_km = y / 1000,
    biomass_positive = as.integer(agb_live_mg_ha > 0)
  ) |>
  arrange(county_fips, year, plot_id)

direct_estimates <- direct_estimates |>
  mutate(
    county_fips = as.character(county_fips),
    year = as.integer(year)
  )

plot_years <- sort(unique(unit_plots$year))
tcc_center <- mean(unit_plots$tcc_mean, na.rm = TRUE)
tcc_scale <- sd(unit_plots$tcc_mean, na.rm = TRUE)

unit_plots <- unit_plots |>
  mutate(tcc_mean_scaled = (tcc_mean - tcc_center) / tcc_scale)

county_levels <- levels(unit_plots$county_fips)
fia_unit_levels <- levels(unit_plots$fia_unit)
fia_unit_names <- unit_plots |>
  distinct(fia_unit, fia_unit_name) |>
  arrange(fia_unit)

prediction_grid <- prediction_grid |>
  mutate(
    county_fips = as.character(county_fips),
    county_fips = factor(county_fips, levels = county_levels),
    fia_unit = factor(fia_unit, levels = fia_unit_levels),
    year = as.integer(year),
    time = as.numeric(year),
    tcc_mean_scaled = (tcc_mean - tcc_center) / tcc_scale
  ) |>
  filter(year %in% plot_years) |>
  arrange(county_fips, year, x, y)

positive_plots <- unit_plots |>
  filter(biomass_positive == 1) |>
  mutate(biomass_sqrt = sqrt(agb_live_mg_ha))
```

The plot rows include `time`, a decimal measurement year used by the continuous space-time process. The prediction grid is annual, so its prediction time is set to the integer `year` for the county-year support being summarized. Plot-level and grid-level TCC values are annual point samples from the USDA Forest Service Tree Canopy Cover raster products [@USFSTreeCanopyCover2026]. The model uses the same plot-based TCC standardization for plot rows and prediction-grid rows.

```{r data-summary}
#| code-summary: "Show data summary code"
data_summary <- tibble(
  quantity = c(
    "plot rows",
    "positive-biomass plot rows",
    "FIA units",
    "5 km prediction-grid rows",
    "measurement years"
  ),
  value = c(
    fmt_int(nrow(unit_plots)),
    fmt_int(nrow(positive_plots)),
    fmt_int(n_distinct(unit_plots$fia_unit)),
    fmt_int(nrow(prediction_grid)),
    paste(range(plot_years), collapse = "-")
  )
)

show_table(data_summary, caption = "Unit-response data for the continuous space-time model with spatial SVC.")
```

```{r plot-map, fig.width=7.2, fig.height=5.8}
#| code-summary: "Show plot map code"
ggplot() +
  geom_sf(data = wa_counties, fill = "grey96", color = "grey70", linewidth = 0.2) +
  geom_point(
    data = unit_plots,
    aes(x, y, color = time),
    size = 0.35,
    alpha = 0.55
  ) +
  coord_sf(expand = FALSE) +
  scale_color_gradientn(colors = stlmm_palette(), name = "measurement time") +
  labs(x = NULL, y = NULL)
```

## Model Settings

The occurrence and positive-biomass intercept processes use `multi_res_sep_exp`, the two-scale separable exponential space-time covariance proposed by @MayFinley2025SpatioTemporal. The covariance is a weighted sum of broad and local separable exponential components. The broad component represents large-scale, slowly varying space-time structure, while the local component represents shorter-range departures. Coordinates are modeled in kilometers and time is modeled in years so the spatial and temporal decay parameters have direct range interpretations.

The positive-stage TCC SVC uses a simpler spatial exponential NNGP. That reflects the working assumption that the positive-biomass relationship with TCC can differ among regions but is not expected to change materially over the measurement period.

`stLMM` parameterizes exponential correlations as `exp(-phi * distance)`. If `range_05` is the distance where correlation reaches 0.05, then `phi = -log(0.05) / range_05`. The broad component uses lower spatial decay and the local component uses higher spatial decay through non-overlapping prior ranges.

In the model settings below, chains are kept intentionally short so the vignette renders quickly. Longer runs with at least 15,000 MCMC samples per chain produced satisfactory convergence diagnostics for all model parameters. The starting values below are based on those longer runs. For an actual analysis, use dispersed starting values and longer runs, with appropriate warmup and convergence checks.

```{r model-settings}
#| code-summary: "Show model settings code"
cor_at_range <- 0.05

decay_from_range <- function(range_05) {
  -log(cor_at_range) / range_05
}

## Broad process: approximately 25-100 km and 1-500 years.
phi_1_bounds <- decay_from_range(rev(c(25, 100)))
lambda_1_bounds <- decay_from_range(rev(c(1, 500)))

## Local process: approximately 1-25 km and 1-500 years.
phi_2_bounds <- decay_from_range(rev(c(1, 25)))
lambda_2_bounds <- decay_from_range(rev(c(1, 500)))

## Spatial SVC process: approximately 10-150 km.
phi_svc_bounds <- decay_from_range(rev(c(10, 150)))

n_samples <- 1000
burnin <- 1
chains <- 3
n_keep <- 25
posterior_thin <- max(1L, floor((n_samples - burnin) / n_keep))
posterior_sub_sample <- list(start = burnin + 1, thin = posterior_thin)
n_report <- 500

## Let stLMM use OpenMP threads while preventing BLAS from nesting threads
## inside the sampler and prediction calculations.
blas_set_num_threads(1)
n_omp_threads <- as.integer(Sys.getenv("STLMM_OMP_THREADS", unset = "25"))

fit_st_scale <- 2
pred_st_scale <- fit_st_scale

sqrt_response_sd <- sd(positive_plots$biomass_sqrt, na.rm = TRUE)
```

```{r cache-settings}
#| include: false
#| purl: false
fit_cache_tag <- paste(
  "two-stage-multi-res-nngp-tcc-fiaunit-var-spatial-svc-short-warmup-mcmc",
  paste(n_samples, burnin, chains, sep = "-"),
  sep = "-"
)

posterior_cache_tag <- paste(
  fit_cache_tag,
  paste0("keep-", n_keep),
  sep = "-"
)

fit_cache_name <- paste0("fits-", fit_cache_tag)
recovery_cache_name <- paste0("recovery-", posterior_cache_tag)
prediction_cache_name <- paste0("prediction-", posterior_cache_tag, "-svc-grid-one-year")
```

The `st_scale` value affects NNGP ordering and neighbor search, not the covariance function itself. Here, one year of separation is treated like about two kilometers for ordering and neighbor search. The covariance still evaluates spatial distance in kilometers and temporal distance in years.

```{r starting-values}
#| code-summary: "Show starting value code"
jitter_start <- function(value, amount = 0.01, bounds = NULL) {
  out <- rep(value, chains) * exp(rnorm(chains, 0, amount))

  if (!is.null(bounds)) {
    out <- pmin(pmax(out, bounds[1] + 1e-8), bounds[2] - 1e-8)
  }

  out
}


occurrence_start <- list(
  beta = c(
    `(Intercept)` = -0.18,
    tcc_mean_scaled = 10
  ),
  nngp_1 = list(
    sigma_sq = 80,
    alpha = 0.52,
    phi_1 = 0.031,
    lambda_1 = 0.039,
    phi_2 = 1.7,
    lambda_2 = 1.6
  )
)

biomass_start <- list(
  resid = list(
    tau_sq = c(
      `5` = 28,
      `6` = 25,
      `7` = 26,
      `8` = 12,
      `9` = 11
    )
  ),
  nngp_1 = list(
    sigma_sq = 2.3,
    alpha = 0.74,
    phi_1 = 0.037,
    lambda_1 = 0.059,
    phi_2 = 1.6,
    lambda_2 = 1.5
  ),
  nngp_2 = list(
    sigma_sq = 2.4,
    phi = 0.027
  )
)

occurrence_starting <- list(
  beta = occurrence_start$beta,
  nngp_1 = list(
    sigma_sq = jitter_start(occurrence_start$nngp_1$sigma_sq),
    alpha = jitter_start(occurrence_start$nngp_1$alpha, bounds = c(0, 1)),
    phi_1 = jitter_start(occurrence_start$nngp_1$phi_1, bounds = phi_1_bounds),
    lambda_1 = jitter_start(occurrence_start$nngp_1$lambda_1, bounds = lambda_1_bounds),
    phi_2 = jitter_start(occurrence_start$nngp_1$phi_2, bounds = phi_2_bounds),
    lambda_2 = jitter_start(occurrence_start$nngp_1$lambda_2, bounds = lambda_2_bounds)
  )
)

biomass_starting <- list(
  resid = biomass_start$resid,
  nngp_1 = list(
    sigma_sq = jitter_start(biomass_start$nngp_1$sigma_sq),
    alpha = jitter_start(biomass_start$nngp_1$alpha, bounds = c(0, 1)),
    phi_1 = jitter_start(biomass_start$nngp_1$phi_1, bounds = phi_1_bounds),
    lambda_1 = jitter_start(biomass_start$nngp_1$lambda_1, bounds = lambda_1_bounds),
    phi_2 = jitter_start(biomass_start$nngp_1$phi_2, bounds = phi_2_bounds),
    lambda_2 = jitter_start(biomass_start$nngp_1$lambda_2, bounds = lambda_2_bounds)
  ),
  nngp_2 = list(
    sigma_sq = jitter_start(biomass_start$nngp_2$sigma_sq),
    phi = jitter_start(biomass_start$nngp_2$phi, bounds = phi_svc_bounds)
  )
)
```

## Fit

The binary model uses all plot rows and models whether `agb_live_mg_ha` is positive.

```{r fit-occurrence}
#| code-fold: false
#| eval: false
#| purl: true
occurrence_fit <- stLMM(
  biomass_positive ~
    tcc_mean_scaled +
    nngp(x_km, y_km, time,
         m = 15,
         cov_model = "multi_res_sep_exp",
         ordering = "maxmin",
         st_scale = fit_st_scale),
  data = unit_plots,
  family = "binomial",
  starting = occurrence_starting,
  priors = list(
    beta = normal(mean = 0, sd = 2.5),
    nngp_1 = list(
      sigma_sq = log_normal(meanlog = log(100), sdlog = 0.15),
      alpha = uniform(0, 1),
      phi_1 = uniform(phi_1_bounds[1], phi_1_bounds[2]),
      lambda_1 = uniform(lambda_1_bounds[1], lambda_1_bounds[2]),
      phi_2 = uniform(phi_2_bounds[1], phi_2_bounds[2]),
      lambda_2 = uniform(lambda_2_bounds[1], lambda_2_bounds[2])
    )
  ),
  n_samples = n_samples,
  chains = chains,
  save_process = posterior_sub_sample,
  n_omp_threads = n_omp_threads,
  verbose = TRUE,
  n_report = n_report
)
```

The positive stage models square-root positive biomass, carries forward FIA-unit-specific residual variances, and adds a positive-stage spatial TCC SVC. The term `tcc_mean_scaled:nngp(x_km, y_km, ...)` is a spatial process multiplied by the standardized TCC value, so it represents regional departures from the global positive-biomass TCC slope.

```{r fit-positive-biomass}
#| code-fold: false
#| eval: false
#| purl: true
biomass_fit <- stLMM(
  biomass_sqrt ~
    tcc_mean_scaled +
    nngp(x_km, y_km, time,
         m = 15,
         cov_model = "multi_res_sep_exp",
         ordering = "maxmin",
         st_scale = fit_st_scale) +
    tcc_mean_scaled:nngp(x_km, y_km,
         m = 15,
         cov_model = "exp",
         ordering = "maxmin") +
    resid(model = "group", group = fia_unit),
  data = positive_plots,
  starting = biomass_starting,
  priors = list(
    beta = normal(mean = 0, sd = 2 * sqrt_response_sd),
    resid = list(tau_sq = half_t(df = 3, scale = sqrt_response_sd)),
    nngp_1 = list(
      sigma_sq = half_t(df = 3, scale = sqrt_response_sd),
      alpha = uniform(0, 1),
      phi_1 = uniform(phi_1_bounds[1], phi_1_bounds[2]),
      lambda_1 = uniform(lambda_1_bounds[1], lambda_1_bounds[2]),
      phi_2 = uniform(phi_2_bounds[1], phi_2_bounds[2]),
      lambda_2 = uniform(lambda_2_bounds[1], lambda_2_bounds[2])
    ),
    nngp_2 = list(
      sigma_sq = half_t(df = 3, scale = 0.5 * sqrt_response_sd),
      phi = uniform(phi_svc_bounds[1], phi_svc_bounds[2])
    )
  ),
  n_samples = n_samples,
  chains = chains,
  n_omp_threads = n_omp_threads,
  verbose = TRUE,
  n_report = n_report
)
```

```{r fit-models-execute}
#| include: false
#| purl: false
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
fit_cache <- use_article_cache(article_id, fit_cache_name, {
  occurrence_fit <- stLMM(
    biomass_positive ~
      tcc_mean_scaled +
      nngp(x_km, y_km, time,
           m = 15,
           cov_model = "multi_res_sep_exp",
           ordering = "maxmin",
           st_scale = fit_st_scale),
    data = unit_plots,
    family = "binomial",
    starting = occurrence_starting,
    priors = list(
      beta = normal(mean = 0, sd = 2.5),
      nngp_1 = list(
        sigma_sq = log_normal(meanlog = log(100), sdlog = 0.15),
        alpha = uniform(0, 1),
        phi_1 = uniform(phi_1_bounds[1], phi_1_bounds[2]),
        lambda_1 = uniform(lambda_1_bounds[1], lambda_1_bounds[2]),
        phi_2 = uniform(phi_2_bounds[1], phi_2_bounds[2]),
        lambda_2 = uniform(lambda_2_bounds[1], lambda_2_bounds[2])
      )
    ),
    n_samples = n_samples,
    chains = chains,
    save_process = posterior_sub_sample,
    n_omp_threads = n_omp_threads,
    verbose = TRUE,
    n_report = n_report
  )

  biomass_fit <- stLMM(
    biomass_sqrt ~
      tcc_mean_scaled +
      nngp(x_km, y_km, time,
           m = 15,
           cov_model = "multi_res_sep_exp",
           ordering = "maxmin",
           st_scale = fit_st_scale) +
      tcc_mean_scaled:nngp(x_km, y_km,
           m = 15,
           cov_model = "exp",
           ordering = "maxmin") +
      resid(model = "group", group = fia_unit),
    data = positive_plots,
    starting = biomass_starting,
    priors = list(
      beta = normal(mean = 0, sd = 2 * sqrt_response_sd),
      resid = list(tau_sq = half_t(df = 3, scale = sqrt_response_sd)),
      nngp_1 = list(
        sigma_sq = half_t(df = 3, scale = sqrt_response_sd),
        alpha = uniform(0, 1),
        phi_1 = uniform(phi_1_bounds[1], phi_1_bounds[2]),
        lambda_1 = uniform(lambda_1_bounds[1], lambda_1_bounds[2]),
        phi_2 = uniform(phi_2_bounds[1], phi_2_bounds[2]),
        lambda_2 = uniform(lambda_2_bounds[1], lambda_2_bounds[2])
      ),
      nngp_2 = list(
        sigma_sq = half_t(df = 3, scale = 0.5 * sqrt_response_sd),
        phi = uniform(phi_svc_bounds[1], phi_svc_bounds[2])
      )
    ),
    n_samples = n_samples,
    chains = chains,
    n_omp_threads = n_omp_threads,
    verbose = TRUE,
    n_report = n_report
  )

  list(
    occurrence_fit = occurrence_fit,
    occurrence_summary = summary(occurrence_fit, burn = burnin),
    biomass_fit = biomass_fit,
    biomass_summary = summary(biomass_fit, burn = burnin)
  )
})

list2env(fit_cache, envir = environment())
```

```{r fit-summaries}
#| code-summary: "Show fit summaries"
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
occurrence_summary$parameters
biomass_summary$parameters
```

The `tcc_mean_scaled` row in the positive-stage summary is the global TCC slope on the square-root positive-biomass scale. The `nngp_2` rows describe the spatial process that modifies that slope. Thus, at a given location, the positive-stage TCC coefficient is the global slope plus the local draw from the SVC process. This SVC is used only in the positive biomass stage; the occurrence stage keeps a single global TCC slope.

## Fitted Space-Time Ranges

The table below reports posterior median practical ranges, defined as the spatial or temporal lag where correlation is approximately 0.05. The positive-stage SVC has its own spatial covariance parameter because it is a second NNGP process.

```{r fitted-covariance-summary}
#| code-summary: "Show fitted covariance summary code"
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
posterior_medians <- function(fit) {
  draws <- as_samples(fit, burn = burnin, metadata = FALSE)
  apply(draws, 2, stats::median)
}

param_value <- function(x, name) {
  unname(x[name])
}

occurrence_medians <- posterior_medians(occurrence_fit)
biomass_medians <- posterior_medians(biomass_fit)

practical_ranges <- function(p, model, process) {
  tibble(
    model = model,
    component = c("broad", "local"),
    spatial_range_km = -log(cor_at_range) / c(
      param_value(p, paste0(process, "_phi_1")),
      param_value(p, paste0(process, "_phi_2"))
    ),
    temporal_range_years = -log(cor_at_range) / c(
      param_value(p, paste0(process, "_lambda_1")),
      param_value(p, paste0(process, "_lambda_2"))
    ),
    variance_weight = c(
      param_value(p, paste0(process, "_alpha")),
      1 - param_value(p, paste0(process, "_alpha"))
    )
  )
}

range_summary <- bind_rows(
  practical_ranges(occurrence_medians, "positive-biomass occurrence", "nngp_1"),
  practical_ranges(biomass_medians, "sqrt positive biomass intercept", "nngp_1"),
  tibble(
    model = "sqrt positive biomass TCC SVC",
    component = "spatial",
    spatial_range_km = -log(cor_at_range) / param_value(biomass_medians, "nngp_2_phi"),
    temporal_range_years = NA_real_,
    variance_weight = NA_real_
  )
) |>
  mutate(across(where(is.numeric), ~ round(.x, 2)))

show_table(
  range_summary,
  caption = "Posterior median practical ranges for the intercept and spatial SVC covariance components."
)
```

```{r fitted-residual-sd-table}
#| code-summary: "Show FIA-unit residual SD table code"
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
biomass_draws <- as_samples(biomass_fit, burn = burnin, metadata = FALSE)
resid_var_names <- paste0("tau_sq_", fia_unit_levels)
if (!all(resid_var_names %in% colnames(biomass_draws))) {
  resid_var_names <- fia_unit_levels
}
resid_sd_draws <- sqrt(biomass_draws[, resid_var_names, drop = FALSE])

fia_resid_sd <- tibble(
  fia_unit = factor(fia_unit_levels, levels = fia_unit_levels),
  sd_mean = colMeans(resid_sd_draws),
  sd_lower = apply(resid_sd_draws, 2, quantile, probs = 0.025),
  sd_upper = apply(resid_sd_draws, 2, quantile, probs = 0.975)
) |>
  left_join(fia_unit_names, by = "fia_unit")

show_table(
  fia_resid_sd |>
    transmute(
      `FIA unit` = fia_unit_name,
      `posterior mean SD` = fmt_num(sd_mean, 2),
      `lower` = fmt_num(sd_lower, 2),
      `upper` = fmt_num(sd_upper, 2)
    ),
  caption = "FIA-unit residual standard deviations for square-root positive biomass."
)
```

## Prediction and Aggregation

The model predicts at all annual 5 km grid rows. For the full grid aggregation, this article uses independent NNGP prediction for computational tractability. Joint Vecchia prediction is the heavier alternative when the primary focus is aggregate posterior uncertainty over a dense prediction lattice. The Gaussian positive-stage predictions use posterior predictive draws on the square-root scale before squaring, so the recombined biomass draws include positive-stage residual uncertainty. A separate one-year prediction keeps draws of the positive-stage SVC process so the grid-level TCC coefficient can be summarized without storing process draws for every annual grid row.

```{r aggregation-helpers}
#| code-summary: "Show aggregation helper code"
prediction_sample_matrix <- function(pred, sample, include_w = FALSE) {
  as.matrix(as_samples(
    pred,
    sample = sample,
    include_w = include_w,
    metadata = FALSE
  ))
}

prediction_w_matrix <- function(pred, term) {
  draws <- prediction_sample_matrix(pred, "mu", include_w = TRUE)
  w_cols <- grep(paste0("^w_", term, "_"), colnames(draws), value = TRUE)

  if (!length(w_cols)) {
    stop("No predicted process samples found for ", term, call. = FALSE)
  }

  draws[, w_cols, drop = FALSE]
}

aggregate_grid_draws <- function(draws, grid, prefix = "") {
  county_year <- grid |>
    distinct(county_fips, county, year) |>
    arrange(county_fips, year)

  county_samples <- matrix(
    NA_real_,
    nrow = nrow(draws),
    ncol = nrow(county_year)
  )

  for (j in seq_len(nrow(county_year))) {
    ii <- which(
      grid$county_fips == county_year$county_fips[j] &
        grid$year == county_year$year[j]
    )
    county_samples[, j] <- rowMeans(draws[, ii, drop = FALSE])
  }

  county_year |>
    bind_cols(
      summarize_draw_matrix(county_samples, prefix = prefix) |>
        select(-prediction_row)
    )
}
```

```{r recover-models}
#| code-fold: false
#| eval: false
#| purl: true
occurrence_rec <- recover(
  occurrence_fit,
  sub_sample = posterior_sub_sample
)

biomass_rec <- recover(
  biomass_fit,
  sub_sample = posterior_sub_sample
)
```

```{r recover-models-execute}
#| include: false
#| purl: false
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
recovery_cache <- use_article_cache(article_id, recovery_cache_name, {
  occurrence_rec <- recover(
    occurrence_fit,
    sub_sample = posterior_sub_sample
  )

  biomass_rec <- recover(
    biomass_fit,
    sub_sample = posterior_sub_sample
  )

  list(
    occurrence_rec = occurrence_rec,
    biomass_rec = biomass_rec
  )
})

list2env(recovery_cache, envir = environment())
```

```{r predict-grid}
#| code-fold: false
#| eval: false
#| purl: true
occurrence_pred <- predict(
  occurrence_rec,
  newdata = prediction_grid,
  y_samples = TRUE,
  st_scale = pred_st_scale,
  return_w_samples = FALSE
)

biomass_pred <- predict(
  biomass_rec,
  newdata = prediction_grid,
  y_samples = TRUE,
  return_w_samples = FALSE
)

svc_map_year <- 2024
svc_prediction_grid <- prediction_grid |>
  filter(year == svc_map_year)

biomass_svc_pred <- predict(
  biomass_rec,
  newdata = svc_prediction_grid,
  y_samples = FALSE,
  return_w_samples = TRUE
)
```

```{r predict-grid-execute}
#| include: false
#| purl: false
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
prediction_cache <- use_article_cache(article_id, prediction_cache_name, {
  map_years <- c(2016, 2020, 2024)
  svc_map_year <- max(map_years)
  svc_prediction_grid <- prediction_grid |>
    filter(year == svc_map_year)

  occurrence_pred <- predict(
    occurrence_rec,
    newdata = prediction_grid,
    y_samples = TRUE,
    st_scale = pred_st_scale,
    return_w_samples = FALSE
  )

  biomass_pred <- predict(
    biomass_rec,
    newdata = prediction_grid,
    y_samples = TRUE,
    return_w_samples = FALSE
  )

  biomass_svc_pred <- predict(
    biomass_rec,
    newdata = svc_prediction_grid,
    y_samples = FALSE,
    return_w_samples = TRUE
  )

  occurrence_prob_draws <- prediction_sample_matrix(occurrence_pred, "mu")
  occurrence_draws <- prediction_sample_matrix(occurrence_pred, "y")
  biomass_sqrt_draws <- prediction_sample_matrix(biomass_pred, "y")
  biomass_svc_draws <- prediction_w_matrix(biomass_svc_pred, "nngp_2")

  n_draw <- min(nrow(occurrence_draws), nrow(biomass_sqrt_draws))
  occurrence_draws <- occurrence_draws[seq_len(n_draw), , drop = FALSE]
  biomass_sqrt_draws <- biomass_sqrt_draws[seq_len(n_draw), , drop = FALSE]

  response_draws <- occurrence_draws * biomass_sqrt_draws^2

  biomass_rec_draws <- as_samples(biomass_rec, metadata = FALSE)
  n_svc_draw <- min(nrow(biomass_svc_draws), nrow(biomass_rec_draws))
  tcc_coef_draws <- sweep(
    biomass_svc_draws[seq_len(n_svc_draw), , drop = FALSE],
    1,
    biomass_rec_draws$tcc_mean_scaled[seq_len(n_svc_draw)],
    `+`
  )
  tcc_coef_summary <- summarize_draw_matrix(tcc_coef_draws, prefix = "tcc_coef_") |>
    select(-prediction_row)

  county_summary <- aggregate_grid_draws(
    response_draws,
    prediction_grid,
    prefix = "theta_"
  )

  map_rows <- which(prediction_grid$year %in% map_years)
  prediction_map <- prediction_grid[map_rows, ] |>
    mutate(
      occurrence_prob_median = apply(
        occurrence_prob_draws[, map_rows, drop = FALSE],
        2,
        median
      ),
      theta_median = apply(
        response_draws[, map_rows, drop = FALSE],
        2,
        median
      )
    )

  tcc_coef_map <- svc_prediction_grid |>
    bind_cols(tcc_coef_summary) |>
    filter(year == svc_map_year)

  rm(
    occurrence_pred, biomass_pred, biomass_svc_pred,
    occurrence_prob_draws, occurrence_draws,
    biomass_sqrt_draws, biomass_svc_draws,
    biomass_rec_draws, tcc_coef_draws,
    response_draws
  )

  list(
    county_summary = county_summary,
    prediction_map = prediction_map,
    tcc_coef_map = tcc_coef_map,
    svc_map_year = svc_map_year
  )
})

list2env(prediction_cache, envir = environment())
```

```{r combine-direct-estimates}
#| code-summary: "Show direct-estimate join code"
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
direct_for_compare <- direct_estimates |>
  filter(year %in% plot_years) |>
  select(
    county_fips, county, year, n,
    direct_biomass_model, direct_biomass_se_model
  )

county_summary <- county_summary |>
  mutate(county_fips = as.character(county_fips)) |>
  left_join(
    direct_for_compare,
    by = c("county_fips", "county", "year")
  )
```

```{r prediction-map, fig.width=9, fig.height=4.8}
#| code-summary: "Show prediction map code"
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
ggplot() +
  geom_raster(
    data = prediction_map,
    aes(x, y, fill = theta_median)
  ) +
  geom_sf(data = wa_counties, fill = NA, color = "grey35", linewidth = 0.15) +
  facet_wrap(~ year, nrow = 1) +
  coord_sf(expand = FALSE) +
  scale_fill_gradientn(
    colors = stlmm_palette(),
    name = "Mg/ha",
    na.value = "grey92"
  ) +
  labs(
    title = "posterior predictive median biomass",
    x = NULL,
    y = NULL
  ) +
  theme(
    plot.title = element_text(size = 10, face = "bold", margin = margin(0, 0, 2, 0)),
    plot.margin = margin(0, 0, 0, 0)
  )
```

The SVC term can also be viewed directly on the prediction grid. The map below summarizes the positive-stage TCC coefficient, which is the global `tcc_mean_scaled` slope plus the predicted local SVC draw. The coefficient is on the square-root positive-biomass scale and is measured per one standard deviation increase in TCC. Because this SVC is spatial-only, the coefficient surface does not change by year; the displayed year is just the grid slice used for mapping.

```{r tcc-svc-grid-map, fig.width=7, fig.height=5.2}
#| code-summary: "Show TCC SVC grid map code"
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
ggplot() +
  geom_raster(
    data = tcc_coef_map,
    aes(x, y, fill = tcc_coef_mean)
  ) +
  geom_sf(data = wa_counties, fill = NA, color = "grey35", linewidth = 0.15) +
  coord_sf(expand = FALSE) +
  scale_fill_gradientn(
    colors = stlmm_continuous_palette(palette = "navia"),
    name = "sqrt(Mg/ha)\nper SD TCC",
    na.value = "grey92"
  ) +
  labs(
    title = "posterior mean positive-stage TCC coefficient",
    subtitle = paste("prediction grid,", svc_map_year),
    x = NULL,
    y = NULL
  ) +
  theme(
    plot.title = element_text(size = 10, face = "bold", margin = margin(0, 0, 2, 0)),
    plot.subtitle = element_text(size = 9, margin = margin(0, 0, 2, 0)),
    plot.margin = margin(0, 0, 0, 0)
  )
```

```{r selected-summary-table}
#| code-summary: "Show selected county-year summary code"
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
show_table(
  county_summary |>
    filter(year == 2024) |>
    arrange(desc(theta_mean)) |>
    select(county, n, direct_biomass_model, theta_mean, theta_lower, theta_upper) |>
    slice_head(n = 10) |>
    mutate(
      across(
        c(direct_biomass_model, theta_mean, theta_lower, theta_upper),
        ~ fmt_num(.x, 1)
      )
    ),
  caption = "Highest 2024 county-year posterior predictive means from the continuous space-time model with spatial SVC."
)
```

## County Time Series

```{r county-time-series, fig.width=9, fig.height=5.8}
#| code-summary: "Show county time-series code"
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
profile_counties <- c("Clallam", "King", "Okanogan", "Yakima")

profile_dat <- county_summary |>
  filter(county %in% profile_counties) |>
  mutate(county = factor(county, levels = profile_counties))

ggplot(profile_dat, aes(year, theta_mean)) +
  geom_ribbon(
    aes(ymin = theta_lower, ymax = theta_upper),
    fill = stlmm_color("primary"),
    alpha = 0.15,
    color = NA
  ) +
  geom_line(color = stlmm_color("primary"), linewidth = 0.75) +
  geom_point(
    aes(y = direct_biomass_model),
    color = "grey20",
    size = 1.4,
    alpha = 0.75,
    na.rm = TRUE
  ) +
  facet_wrap(~ county, scales = "free_y", ncol = 2) +
  labs(
    x = "year",
    y = "mean agb_live_mg_ha (Mg/ha)"
  )
```

The direct estimates remain points of comparison rather than the model target. The main comparison with the previous unit-response articles is how the continuous space-time process changes the fitted county-year trajectories after prediction and aggregation over the same grid-supported response.

## Takeaways

This article replaces county-year areal smoothing with continuous space-time NNGP intercept processes. The response definition, TCC covariate, two-stage product, FIA-unit residual variance in the positive stage, and county-year aggregation target are carried forward from the previous unit-response articles. The added SVC makes the positive-biomass TCC effect vary over continuous space.

The continuous process is useful when the model should borrow over distance and time without forcing all variation through county-year areal effects. The tradeoff is computational and inferential: dense grid prediction is more expensive than the areal CAR-time workflow, and the richer positive-stage mean model should be checked carefully before treating it as the preferred analysis model.