---
title: "Unit-response two-stage CAR-time model with TCC"
subtitle: "Plot-level Washington biomass prediction with county-time structure"
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(sf)
library(tidyverse)
library(knitr)
library(kableExtra)
library(patchwork)

source("article-utils.R")

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

article_id <- "wa-unit-two-stage-car-time-tcc"
data_dir <- "wa_data"
```

## Purpose

This article continues the unit-response part of the Washington county biomass series. The opening unit-response article introduced BHF-style models with row-level TCC, county effects, and a shared AR(1) year effect. Here the model replaces that simple county-plus-year structure with a county-year effect, so each county can have its own temporal departure while still borrowing across the county graph.

The response is the plot-level `agb_live_mg_ha` variable used throughout the series. To respect the point mass at zero, the model is two-stage from the start:

1. a Bernoulli model for `agb_live_mg_ha > 0`;
2. a Gaussian model for square-root transformed positive `agb_live_mg_ha`.

Both stages include row-level tree canopy cover (TCC) and a separable county-year CAR-time process. This keeps the model aligned with the same series target used in the area-response articles, while giving the unit-response model a more useful spatial-temporal structure. The two-stage product follows the same general hurdle-style strategy used in related forest-biomass models [@Finley01032011; @MayFinley2025SpatioTemporal].

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

- [Quarto source](wa-unit-two-stage-car-time-tcc.qmd)
- [Extracted R code](wa-unit-two-stage-car-time-tcc.R)
:::

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

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

## Working Target

Let $j$ index an observed FIA plot row, with county $a_j$, measurement year $t_j$, and standardized row-level TCC covariate $x_j$. The plot-level response is $y_j =$ `agb_live_mg_ha`.

The county-year target is the mean of `agb_live_mg_ha` over grid rows assigned to county $a$ and year $t$. Let $g$ index those grid rows. The model predicts $\tilde{y}_g$ on the original Mg/ha scale and averages posterior predictive draws over grid rows within each county-year.

For equivalence with the area-response direct estimates in this series, the two-stage decomposition is based on `agb_live_mg_ha` itself,

$$
y_j = z_j y_j^+,\qquad z_j = \mathbb{1}(y_j > 0),
$$

where $y_j^+$ is the positive biomass value for rows with $y_j > 0$.

The binary stage is

$$
z_j \sim \operatorname{bernoulli}(p_j),\qquad
\operatorname{logit}(p_j) =
\gamma_0 + \gamma_1 x_j + w^z_{a_j,t_j},
$$

where $w^z_{a,t}$ is a separable CAR-time process over counties and years. The positive-biomass stage models the square root of positive `agb_live_mg_ha`,

$$
q_j = \sqrt{y_j^+},\qquad
q_j = \delta_0 + \delta_1 x_j + w^+_{a_j,t_j} + r_j,\qquad
r_j \sim \operatorname{normal}(0,\tau_+^2),
$$

where $w^+_{a,t}$ is a second CAR-time process for the positive biomass magnitude. The two stages are fit separately and recombined after prediction.

For posterior prediction, superscript $(m)$ indexes one retained posterior predictive draw. The grid-row response draw is

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

and the county-year draw is obtained by averaging these grid-row draws within county $a$ and year $t$.

## Data

The plot and grid TCC variables are annual point samples from the USDA Forest Service Tree Canopy Cover raster products [@USFSTreeCanopyCover2026]. The model uses the same standardization for plot rows and prediction-grid rows, with the center and scale estimated from the observed plot rows.

```{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))

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

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),
    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 5 km grid is the prediction support for this single-model unit-response article. The 10 km grid remains useful for fast local checks.

```{r data-summary}
#| code-summary: "Show data summary code"
data_summary <- tibble(
  quantity = c(
    "plot rows",
    "counties",
    "measurement years",
    "positive-biomass plot rows",
    "zero-biomass plot rows",
    "5 km prediction-grid rows",
    "plot TCC mean",
    "plot TCC SD"
  ),
  value = c(
    fmt_int(nrow(unit_plots)),
    fmt_int(n_distinct(unit_plots$county_fips)),
    paste(range(plot_years), collapse = "-"),
    fmt_int(sum(unit_plots$biomass_positive == 1)),
    fmt_int(sum(unit_plots$biomass_positive == 0)),
    fmt_int(nrow(prediction_grid)),
    fmt_num(tcc_center, 1),
    fmt_num(tcc_scale, 1)
  )
)

show_table(data_summary, caption = "Unit-response data for the two-stage CAR-time TCC model.")
```

## County Graph

The county graph defines the spatial component of both CAR-time processes. Washington has island components, so the graph construction adds nearest-neighbor bridge edges.

```{r county-graph}
#| code-fold: false
g <- car_graph(
  wa_counties,
  id = "county_fips",
  island = "nearest",
  island_k = 4
)

g$island_added_edges
```

The figure checks the graph visually. Thin lines connect neighboring county centroids; thicker colored lines show bridge edges added for island components.

```{r county-graph-map, fig.width=7, fig.height=5.2}
#| code-summary: "Show graph plotting code"
coord_dat <- data.frame(
  county_fips = wa_counties$county_fips,
  sf::st_coordinates(sf::st_point_on_surface(sf::st_geometry(wa_counties)))
)

edge_index <- which(as.matrix(g$adjacency) != 0, arr.ind = TRUE)
edge_index <- edge_index[edge_index[, "row"] < edge_index[, "col"], , drop = FALSE]
edge_dat <- tibble(
  from = g$ids[edge_index[, "row"]],
  to = g$ids[edge_index[, "col"]]
) |>
  mutate(
    key = paste(pmin(from, to), pmax(from, to), sep = "--"),
    x = coord_dat$X[match(from, coord_dat$county_fips)],
    y = coord_dat$Y[match(from, coord_dat$county_fips)],
    xend = coord_dat$X[match(to, coord_dat$county_fips)],
    yend = coord_dat$Y[match(to, coord_dat$county_fips)]
  )

island_key <- paste(
  pmin(g$island_added_edges$from, g$island_added_edges$to),
  pmax(g$island_added_edges$from, g$island_added_edges$to),
  sep = "--"
)
edge_dat$island <- edge_dat$key %in% island_key

ggplot(wa_counties) +
  geom_sf(fill = "grey96", color = "white", linewidth = 0.15) +
  geom_segment(
    data = edge_dat,
    aes(x = x, y = y, xend = xend, yend = yend, color = island, linewidth = island)
  ) +
  geom_point(data = coord_dat, aes(X, Y), color = stlmm_color("primary"), size = 1.3) +
  coord_sf(expand = FALSE) +
  scale_color_manual(
    values = c("FALSE" = "grey45", "TRUE" = stlmm_color("secondary")),
    guide = "none"
  ) +
  scale_linewidth_manual(values = c("FALSE" = 0.25, "TRUE" = 1), guide = "none") +
  theme(plot.margin = margin(0, 0, 0, 0))
```

## Fit

This article fits one unit-response model, but the model has two likelihood stages. The binary stage and positive-biomass stage have separate fixed effects and separate CAR-time processes.

```{r mcmc-settings}
#| code-summary: "Show MCMC settings"
n_samples <- 10000
burnin <- 5000
chains <- 3
n_keep <- 100
posterior_thin <- max(1L, floor((n_samples - burnin) / n_keep))
posterior_sub_sample <- list(start = burnin + 1, thin = posterior_thin)
chain_control <- list(seed = 13, dispersion = 1.5)
warmup_control <- list(batch_length = 25, min_batches = 10)

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

fit_summary_parameters <- list(
  occurrence = c(
    "(Intercept)", "tcc_mean_scaled",
    "car_time_1_sigma_sq", "car_time_1_rho", "car_time_1_phi"
  ),
  biomass = c(
    "(Intercept)", "tcc_mean_scaled", "tau_sq",
    "car_time_1_sigma_sq", "car_time_1_rho", "car_time_1_phi"
  )
)

select_summary_parameters <- function(x, parameters) {
  keep <- intersect(parameters, rownames(x$parameters))
  x$parameters <- x$parameters[keep, , drop = FALSE]
  if (!is.null(x$diagnostics)) {
    x$diagnostics <- x$diagnostics[x$diagnostics$parameter %in% keep, , drop = FALSE]
  }
  x
}
```

```{r cache-settings}
#| include: false
#| purl: false
fit_cache_tag <- paste(
  "two-stage-car-time-tcc-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)
```

### Positive-Biomass Occurrence

The Bernoulli model uses all plot rows. The `family = "binomial"` argument gives a Bernoulli likelihood because no trial-count column is supplied.

```{r fit-occurrence}
#| code-fold: false
#| eval: false
#| purl: true
occurrence_fit <- stLMM(
  biomass_positive ~
    tcc_mean_scaled +
    car_time(county_fips, year, graph = g, car_model = "leroux"),
  data = unit_plots,
  family = "binomial",
  priors = list(
    beta = normal(mean = 0, sd = 2.5),
    car_time_1 = list(
      sigma_sq = half_t(df = 3, scale = 1),
      rho = uniform(0.01, 0.99),
      phi = uniform(-0.99, 0.99)
    )
  ),
  n_samples = n_samples,
  chains = chains,
  chain_control = chain_control,
  warmup = warmup_control,
  verbose = TRUE,
  n_report = 500
)

occurrence_summary <- summary(
  occurrence_fit,
  burn = burnin
) |>
  select_summary_parameters(fit_summary_parameters$occurrence)
```

### Positive Biomass

The positive-biomass stage is fit only to rows where `agb_live_mg_ha` is positive. The square-root transformation keeps the Gaussian stage closer to symmetric, and prediction draws are squared when returning to Mg/ha.

```{r fit-positive-biomass}
#| code-fold: false
#| eval: false
#| purl: true
biomass_fit <- stLMM(
  biomass_sqrt ~
    tcc_mean_scaled +
    car_time(county_fips, year, graph = g, car_model = "leroux"),
  data = positive_plots,
  priors = list(
    beta = normal(mean = 0, sd = 2 * sqrt_response_sd),
    resid = list(tau_sq = half_t(df = 3, scale = sqrt_response_sd)),
    car_time_1 = list(
      sigma_sq = half_t(df = 3, scale = sqrt_response_sd),
      rho = uniform(0.01, 0.99),
      phi = uniform(-0.99, 0.99)
    )
  ),
  n_samples = n_samples,
  chains = chains,
  chain_control = chain_control,
  warmup = warmup_control,
  verbose = TRUE,
  n_report = 500
)

biomass_summary <- summary(
  biomass_fit,
  burn = burnin
) |>
  select_summary_parameters(fit_summary_parameters$biomass)
```

```{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 +
      car_time(county_fips, year, graph = g, car_model = "leroux"),
    data = unit_plots,
    family = "binomial",
    priors = list(
      beta = normal(mean = 0, sd = 2.5),
      car_time_1 = list(
        sigma_sq = half_t(df = 3, scale = 1),
        rho = uniform(0.01, 0.99),
        phi = uniform(-0.99, 0.99)
      )
    ),
    n_samples = n_samples,
    chains = chains,
    chain_control = chain_control,
    warmup = warmup_control,
    verbose = TRUE,
    n_report = 500
  )

  biomass_fit <- stLMM(
    biomass_sqrt ~
      tcc_mean_scaled +
      car_time(county_fips, year, graph = g, car_model = "leroux"),
    data = positive_plots,
    priors = list(
      beta = normal(mean = 0, sd = 2 * sqrt_response_sd),
      resid = list(tau_sq = half_t(df = 3, scale = sqrt_response_sd)),
      car_time_1 = list(
        sigma_sq = half_t(df = 3, scale = sqrt_response_sd),
        rho = uniform(0.01, 0.99),
        phi = uniform(-0.99, 0.99)
      )
    ),
    n_samples = n_samples,
    chains = chains,
    chain_control = chain_control,
    warmup = warmup_control,
    verbose = TRUE,
    n_report = 500
  )

  list(
    occurrence_fit = occurrence_fit,
    occurrence_summary = summary(
      occurrence_fit,
      burn = burnin
    ) |>
      select_summary_parameters(fit_summary_parameters$occurrence),
    biomass_fit = biomass_fit,
    biomass_summary = summary(
      biomass_fit,
      burn = burnin
    ) |>
      select_summary_parameters(fit_summary_parameters$biomass)
  )
})

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
biomass_summary
```

The fit summaries focus on fixed effects and variance parameters. The many county-year process coefficients are recovered and summarized through predictions instead of printed in the main table.

## Prediction and Aggregation

The 5 km grid rows define the county-year prediction support. For each posterior predictive draw, the binary and positive stages are combined on each grid row and then averaged within county-year. The Gaussian positive-stage predictions use `y_samples = TRUE`, so the squared positive-biomass draws include residual predictive uncertainty before recombination.

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

Structured CAR-time process terms are collapsed during fitting. Recovery adds the county-year process draws needed for prediction.

```{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())
```

Prediction uses the same plot-based TCC standardization that was used for fitting.

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

biomass_pred <- predict(
  biomass_rec,
  newdata = prediction_grid,
  y_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, {
  occurrence_pred <- predict(
    occurrence_rec,
    newdata = prediction_grid,
    y_samples = TRUE
  )

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

  occurrence_draws <- prediction_sample_matrix(occurrence_pred, "y")
  occurrence_prob_draws <- prediction_sample_matrix(occurrence_pred, "mu")
  biomass_sqrt_draws <- prediction_sample_matrix(biomass_pred, "y")

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

  positive_biomass_draws <- biomass_sqrt_draws^2
  response_draws <- occurrence_draws * positive_biomass_draws

  occurrence_county_summary <- aggregate_grid_draws(
    occurrence_prob_draws,
    prediction_grid,
    prefix = "p_positive_"
  )

  positive_biomass_county_summary <- aggregate_grid_draws(
    positive_biomass_draws,
    prediction_grid,
    prefix = "positive_"
  )

  county_summary <- aggregate_grid_draws(
    response_draws,
    prediction_grid,
    prefix = "theta_"
  ) |>
    left_join(
      occurrence_county_summary,
      by = c("county_fips", "county", "year")
    ) |>
    left_join(
      positive_biomass_county_summary,
      by = c("county_fips", "county", "year")
    )

  rm(
    occurrence_pred, biomass_pred,
    occurrence_draws, occurrence_prob_draws, biomass_sqrt_draws,
    positive_biomass_draws, response_draws
  )

  list(
    county_summary = county_summary,
    occurrence_county_summary = occurrence_county_summary,
    positive_biomass_county_summary = positive_biomass_county_summary
  )
})

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

The resulting county-year summary has the common shape used throughout the series: county, year, direct-estimate information where available, and posterior summaries for the model-based county-year mean.

```{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 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,
      p_positive_mean, positive_mean
    ) |>
    slice_head(n = 10) |>
    mutate(
      across(
        c(
          direct_biomass_model,
          theta_mean, theta_lower, theta_upper,
          p_positive_mean, positive_mean
        ),
        ~ fmt_num(.x, 1)
      )
    ),
  caption = "Highest 2024 county-year posterior predictive means from the two-stage CAR-time TCC model."
)
```

## Maps

The maps separate the two model stages before showing the recombined series target. The first row is the average positive-biomass probability, the second row is the positive-biomass magnitude from the Gaussian stage, and the third row is the recombined county-year mean. The two Mg/ha rows use the same color scale, anchored by the positive-biomass magnitude, so the recombined mean can be read as a damped product of occurrence and positive magnitude.

```{r prediction-map, fig.width=10, fig.height=7.2}
#| code-summary: "Show prediction map code"
#| eval: !expr identical(tolower(Sys.getenv("STLMM_RUN_LONG_ARTICLES", unset = "false")), "true")
panel_years <- c(2016, 2018, 2020, 2022, 2024)

map_dat <- wa_counties |>
  left_join(
    county_summary |> filter(year %in% panel_years),
    by = "county_fips"
  )

biomass_map_limits <- c(0, max(map_dat$positive_mean, na.rm = TRUE))

component_map <- function(fill_var, legend_title, title, limits = NULL) {
  ggplot(map_dat) +
    geom_sf(aes(fill = .data[[fill_var]]), color = "grey80", linewidth = 0.12) +
    coord_sf(expand = FALSE) +
    facet_wrap(~ year, nrow = 1) +
    scale_fill_gradientn(
      colors = stlmm_palette(),
      name = legend_title,
      limits = limits,
      na.value = "grey92"
    ) +
    labs(title = title) +
    theme(
      panel.spacing = grid::unit(0.03, "lines"),
      strip.text = element_text(margin = margin(1, 1, 1, 1)),
      plot.title = element_text(size = 10, face = "bold", margin = margin(0, 0, 2, 0)),
      plot.margin = margin(0, 0, 0, 0)
    )
}

component_map(
  "p_positive_mean",
  "prob.",
  "positive-biomass probability"
) /
  component_map(
    "positive_mean",
    "Mg/ha",
    "positive-biomass magnitude",
    limits = biomass_map_limits
  ) /
  component_map(
    "theta_mean",
    "Mg/ha",
    "recombined county-year mean",
    limits = biomass_map_limits
  )
```

## County Time Series

The time-series figure uses the same four counties as the opening unit-response article. Direct estimates are shown where county-year direct estimates are model-ready.

```{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)"
  )
```

As in the previous BHF analysis, the Yakima predictions are consistently below many of the direct estimates. Importantly, the direct estimates are not the target for these models; we include them to show where the unit-response summaries agree with the direct estimates and where they differ. Yakima is useful for seeing the difference between fitting the observed plot rows and aggregating over a prediction support. In the data used here, Yakima's observed plot rows have higher mean TCC than the prediction grid. Averaged over the modeled years, the observed Yakima plot rows have mean TCC about 31, while the 5 km prediction grid has mean TCC about 25. Because TCC is positively related to both positive-biomass occurrence and positive-biomass magnitude, applying the model over the lower-TCC grid support shifts the county-year prediction downward.

This is not a failure of the fitted model. If the biomass-TCC relationship is estimated well and the grid is the intended support, the unit-response model may be compensating for a sampled plot support that is not representative of the county-wide TCC distribution. The important caveats are that this interpretation depends on the model: TCC must be an adequate adjustment variable, the fitted relationship must transport from sampled plots to grid cells, and the grid must represent the target support we intend to estimate. Under those assumptions, disagreement with the simple direct estimate can be a useful feature of the unit-response approach rather than a symptom that the model missed Yakima.

## Takeaways

This model is the first structured unit-response analogue of the area-response CAR-time plus TCC articles. It keeps the same series target, uses `agb_live_mg_ha > 0` for the two-stage split, and adds a county-year process so counties can have their own temporal patterns.

The next unit-response article can keep this same two-stage CAR-time TCC mean structure and add FIA-unit-specific residual variance in the positive Gaussian stage. The binary stage has no Gaussian residual variance term, so the heteroskedastic extension belongs only in the positive-biomass stage.