---
title: "Direct-estimate CAR-time model"
subtitle: "Area-response smoothing over counties and years"
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)

source("article-utils.R")

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

article_id <- "wa-direct-car-time"
data_dir <- "wa_data"
```

## Purpose

This article fits the first model in the Washington county biomass series. The response is a county-year direct estimate of live aboveground biomass density. The model treats each direct estimate as a noisy observation of a latent county-year biomass mean and uses a CAR-time random effect to borrow strength across neighboring counties and adjacent years.

This is an **area-response** model: each row represents a county-year, not a plot. It follows the Fay-Herriot small-area model structure [@FayHerriot1979], with a CAR-time process replacing an independent area effect.

Source code and data links for reproducing this article are collected below.

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

- [Quarto source](wa-direct-car-time.qmd)
- [Extracted R code](wa-direct-car-time.R)
:::

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

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

## Model

Let $i$ index an observed county-year direct estimate, with county $a_i$ and year $t_i$. The direct estimate $\hat{y}_i$ is modeled as

$$
\hat{y}_i = \theta_i + e_i,\qquad e_i \sim N(0,\hat{v}_i),
$$

where $\theta_i = \theta_{a_i,t_i}$ is the latent county-year biomass mean and $\hat{v}_i$ is the direct-estimate variance. The latent mean model is

$$
\theta_i = \beta_0 + w_{a_i,t_i},
$$

where $w_{a,t}$ is a separable CAR-time process over counties and years. Its spatial component is defined by the county graph, and its temporal component is defined over the annual support.

## Data

The data chunk below reads the county polygons and county-year direct estimates. Rows with missing `direct_biomass_fh` are retained because they define county-year latent support for recovery and fitted values, but they do not contribute to the likelihood.

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

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

direct_estimates <- direct_estimates |>
  mutate(
    county_fips = as.character(county_fips),
    year = as.integer(year),
    direct_estimate_in_fh = direct_biomass_vhat_source == "direct_biomass_var",
    direct_biomass_fh = if_else(direct_estimate_in_fh, direct_biomass, NA_real_),
    direct_biomass_se_fh = if_else(direct_estimate_in_fh, direct_biomass_se, NA_real_),
    direct_biomass_vhat_fh = if_else(direct_estimate_in_fh, direct_biomass_var, NA_real_)
  ) |>
  arrange(county_fips, year)

observed_direct <- direct_estimates |>
  filter(!is.na(direct_biomass_fh))
```

The table summarizes the area-response data used in this article. County-years with a single FIA plot remain in the county-year support, but they are not used as direct-estimate likelihood rows because their design-based variance cannot be estimated from one plot. Multi-plot county-years with zero raw direct-estimate variance are also excluded here because this fixed-variance model treats each supplied variance as known.

```{r data-summary}
#| code-summary: "Show data summary code"
direct_summary <- tibble(
  quantity = c(
    "counties",
    "years",
    "county-years",
    "positive-variance direct estimates",
    "missing model responses",
    "single-plot rows retained but excluded",
    "zero-variance rows retained but excluded",
    "median plot count among modeled county-years",
    "median direct-estimate SE"
  ),
  value = c(
    n_distinct(direct_estimates$county_fips),
    n_distinct(direct_estimates$year),
    nrow(direct_estimates),
    sum(direct_estimates$direct_estimate_in_fh),
    sum(is.na(direct_estimates$direct_biomass_fh)),
    sum(direct_estimates$direct_estimate_status == "single_plot_no_variance"),
    sum(direct_estimates$direct_biomass_vhat_source == "floor", na.rm = TRUE),
    median(observed_direct$n, na.rm = TRUE),
    median(observed_direct$direct_biomass_se_fh, na.rm = TRUE)
  )
) |>
  mutate(value = fmt_num(as.numeric(value), 1))

show_table(direct_summary, caption = "County-year direct-estimate data used in the CAR-time model.")
```

The next table shows selected-year data coverage. The model includes all counties in every year, but years with fewer positive-variance direct estimates rely more heavily on the model structure.

```{r year-coverage}
#| code-summary: "Show year coverage code"
year_coverage <- direct_estimates |>
  group_by(year) |>
  summarise(
    counties_with_direct_estimates = sum(direct_estimate_in_fh),
    median_plot_n = median(n[n > 0], na.rm = TRUE),
    median_direct_biomass = median(direct_biomass_fh, na.rm = TRUE),
    median_direct_se = median(direct_biomass_se_fh, na.rm = TRUE),
    .groups = "drop"
  )

show_table(
  year_coverage |>
    filter(year %in% c(2016, 2018, 2020, 2022, 2024, 2025)) |>
    mutate(
      median_plot_n = fmt_num(median_plot_n, 0),
      median_direct_biomass = fmt_num(median_direct_biomass, 1),
      median_direct_se = fmt_num(median_direct_se, 1)
    ),
  caption = "Selected-year direct-estimate coverage and precision."
)
```

## County Graph

The county graph defines spatial neighbors for the CAR component. Washington has island components, so the graph construction explicitly adds nearest-neighbor bridge edges. The island rule is part of the model definition.

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

## Fit

The sampler settings below are chosen to give stable posterior summaries for the example while keeping the workflow practical to rerun. For a final analysis, increase the run length as needed and check convergence diagnostics before interpreting results.

```{r sampler-settings}
#| code-fold: false
n_samples <- 12000
burnin <- 6000
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 = 1, dispersion = 1.5)
warmup_control <- list(batch_length = 25, min_batches = 10)
```

The `n_keep` value controls how many post-burn-in draws per chain are retained for recovery and fitted-value summaries. The fit summary uses all post-burn-in MCMC samples.

```{r cache-settings}
#| include: false
#| purl: false
fit_cache_tag <- paste(
  "positive-raw-var-mcmc",
  paste(n_samples, burnin, chains, sep = "-"),
  sep = "-"
)

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

In the model call below, `direct_biomass_fh` is $\hat{y}_i$, and `direct_biomass_vhat_fh` is $\hat{v}_i$. Fixed direct-estimate variances enter through `resid(model = "fixed", variance = direct_biomass_vhat_fh)`.

```{r fit-model}
#| code-fold: false
#| eval: false
#| purl: true
fit <- stLMM(
  direct_biomass_fh ~
    car_time(county_fips, year, graph = g, car_model = "leroux") +
    resid(model = "fixed", variance = direct_biomass_vhat_fh),
  data = direct_estimates,
  priors = list(
    car_time_1 = list(
      sigma_sq = half_t(
        df = 3,
        scale = sd(observed_direct$direct_biomass_fh, na.rm = TRUE)
      ),
      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
)

fit_summary <- summary(fit, burn = burnin)
```

```{r fit-model-execute}
#| include: false
#| purl: false
fit_cache <- use_article_cache(article_id, paste0("fit-", fit_cache_tag), {
  fit <- stLMM(
    direct_biomass_fh ~
      car_time(county_fips, year, graph = g, car_model = "leroux") +
      resid(model = "fixed", variance = direct_biomass_vhat_fh),
    data = direct_estimates,
    priors = list(
      car_time_1 = list(
        sigma_sq = half_t(
          df = 3,
          scale = sd(observed_direct$direct_biomass_fh, na.rm = TRUE)
        ),
        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(
    fit = fit,
    fit_summary = summary(fit, burn = burnin)
  )
})

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

```{r fit-summary}
#| code-summary: "Show fit summary code"
#| purl: true
fit_summary
```

## Recovery and Fitted Means

Structured process terms are collapsed during Gaussian model fitting. The `recover()` call draws the county-year CAR-time effects from retained posterior parameter draws so fitted values can include the latent process contribution.

```{r recover-model}
#| code-fold: false
#| eval: false
#| purl: true
rec <- recover(
  fit,
  sub_sample = posterior_sub_sample
)
```

```{r recover-model-execute}
#| include: false
#| purl: false
recovery_cache <- use_article_cache(article_id, paste0("recovery-", posterior_cache_tag), {
  rec <- recover(
    fit,
    sub_sample = posterior_sub_sample
  )

  list(rec = rec)
})

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

The fitted values are posterior draws of the latent county-year mean for the same county-year support used in the fit. This support includes county-years with missing direct estimates. Those rows did not contribute to the likelihood, but their CAR-time effects can still be recovered because their county and year define latent support nodes.

```{r fitted-county-years}
#| code-summary: "Show fitted-value summary code"
#| code-fold: true
#| eval: false
#| purl: true
fitted_draws <- as.matrix(as_samples(fitted(rec, summary = FALSE), metadata = FALSE))

county_year_summary <- direct_estimates |>
  bind_cols(
    summarize_draw_matrix(fitted_draws, prefix = "theta_") |>
      select(-prediction_row)
  )
```

```{r fitted-county-years-execute}
#| include: false
#| purl: false
fitted_cache <- use_article_cache(article_id, paste0("fitted-", posterior_cache_tag), {
  fitted_draws <- as.matrix(as_samples(fitted(rec, summary = FALSE), metadata = FALSE))

  county_year_summary <- direct_estimates |>
    bind_cols(
      summarize_draw_matrix(fitted_draws, prefix = "theta_") |>
        select(-prediction_row)
    )

  list(
    fitted_draws = fitted_draws,
    county_year_summary = county_year_summary
  )
})

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

The resulting summary table has the common output shape used throughout the series: one row per county-year, direct-estimate information where available, and posterior means and 95% credible intervals for the model-based county-year biomass mean.

```{r selected-summary-table}
#| code-summary: "Show selected county-year summary code"
show_table(
  county_year_summary |>
    filter(year == 2024) |>
    arrange(desc(theta_mean)) |>
    select(
      county, n, direct_biomass_fh, direct_biomass_se_fh,
      theta_mean, theta_lower, theta_upper
    ) |>
    slice_head(n = 10) |>
    mutate(
      across(
        c(direct_biomass_fh, direct_biomass_se_fh, theta_mean, theta_lower, theta_upper),
        ~ fmt_num(.x, 1)
      )
    ),
  caption = "Highest 2024 model-smoothed county-year biomass means, with 95% posterior credible intervals."
)
```

## Direct Estimates and Smoothed Means

The next figure compares positive-variance direct estimates with model-smoothed county-year means. Counties with no positive-variance direct estimate in a year can still receive a model-based estimate because their county-year latent effect is part of the CAR-time support.

```{r direct-vs-smoothed-map, fig.width=10, fig.height=3.8}
#| code-summary: "Show direct-versus-smoothed map code"
panel_years <- c(2016, 2018, 2020, 2022, 2024)

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

map_long <- bind_rows(
  map_dat |> mutate(quantity = "direct estimate", value = direct_biomass_fh),
  map_dat |> mutate(quantity = "model-smoothed\nmean", value = theta_mean)
) |>
  mutate(quantity = factor(quantity, levels = c("direct estimate", "model-smoothed\nmean")))

ggplot(map_long) +
  geom_sf(aes(fill = value), color = "grey80", linewidth = 0.12) +
  coord_sf(expand = FALSE) +
  facet_grid(quantity ~ year) +
  scale_fill_gradientn(
    colors = stlmm_palette(),
    name = "Mg/ha",
    na.value = "grey92"
  ) +
  theme(
    panel.spacing = grid::unit(0.02, "lines"),
    strip.text = element_text(margin = margin(1, 1, 1, 1)),
    plot.margin = margin(0, 0, 0, 0)
  )
```

The scatter plot shows the same comparison only for county-years used in the direct-estimate likelihood. The model is not trying to reproduce every noisy direct estimate; it smooths them according to the supplied sampling variances and the county-time dependence structure.

```{r direct-vs-smoothed-scatter, fig.width=6.5, fig.height=5.5}
#| code-summary: "Show direct-versus-smoothed scatter code"
scatter_dat <- county_year_summary |>
  filter(!is.na(direct_biomass_fh))

fit_axis_limits <- range(
  c(scatter_dat$direct_biomass_fh, scatter_dat$theta_mean, 0),
  na.rm = TRUE
)

ggplot(scatter_dat, aes(x = direct_biomass_fh, y = theta_mean)) +
  geom_abline(slope = 1, intercept = 0, color = "grey45", linewidth = 0.45) +
  geom_point(
    aes(size = 1 / direct_biomass_vhat_fh, color = year),
    alpha = 0.55
  ) +
  scale_size_continuous(range = c(0.5, 3.5), guide = "none") +
  scale_color_viridis_c(option = "C") +
  coord_equal(xlim = fit_axis_limits, ylim = fit_axis_limits) +
  labs(
    x = "Direct estimate Mg/ha",
    y = "CAR-time posterior mean Mg/ha",
    color = "Year"
  )
```

## County Time Series

Time series show how the model smooths through years with sparse or missing positive-variance direct estimates. The example counties below are selected to include different FIA units and biomass ranges.

```{r county-time-series, fig.width=10, fig.height=6.5}
#| code-summary: "Show county time-series code"
selected_counties <- c("Clallam", "King", "Okanogan", "Yakima")

series_dat <- county_year_summary |>
  filter(county %in% selected_counties) |>
  mutate(county = factor(county, levels = selected_counties))

ggplot(series_dat, aes(x = year)) +
  geom_ribbon(
    aes(ymin = theta_lower, ymax = theta_upper),
    fill = "#9ecae1",
    alpha = 0.45
  ) +
  geom_line(aes(y = theta_mean), color = stlmm_color("primary"), linewidth = 0.7) +
  geom_point(
    aes(y = direct_biomass_fh, size = n),
    color = stlmm_color("secondary"),
    alpha = 0.75,
    na.rm = TRUE
  ) +
  facet_wrap(~ county, ncol = 2, scales = "free_y") +
  scale_x_continuous(breaks = seq(min(series_dat$year), max(series_dat$year), by = 5)) +
  scale_size_continuous(range = c(1, 3.5), name = "Plot n") +
  labs(
    x = "Year",
    y = "Live aboveground biomass Mg/ha"
  )
```

The ribbons are 95% posterior credible intervals for the latent county-year mean. Points are positive-variance direct estimates, with point size proportional to the number of FIA plot rows used in that county-year direct estimate. These intervals do not include posterior predictive error for a new direct estimate; the direct-estimate observation variance is represented by the noisy points, not added to the ribbon.

## Fixed Direct-Estimate Variances

This article follows the classical Fay-Herriot convention of treating the supplied direct-estimate variances as fixed. In that formulation, `direct_biomass_vhat_fh` describes uncertainty in the direct estimate $\hat{y}_i$, not residual variation in the latent county-year mean $\theta_i$ after the CAR-time process is included.

This distinction matters because the direct-estimate variances are also estimated from the sample. For county-years with few FIA plots, those variance estimates can be noisy. Treating them as fixed can therefore make the model more certain about some direct estimates than the data support.

The next direct-estimate article relaxes this assumption by using the supplied variance estimates to inform, rather than fix, the observation variances.

## Summary

This first article establishes the area-response workflow:

- the response rows are positive-variance county-year direct estimates;
- fixed direct-estimate variances enter through `resid(model = "fixed", variance = ...)`;
- the county graph defines spatial borrowing;
- `car_time()` adds separable county-time smoothing;
- `recover()` is needed before computing fitted values with Gaussian structured-process models;
- the common output is a county-year table of posterior means and 95% credible intervals.

The next area-response article relaxes the fixed-variance assumption by treating the supplied direct-estimate variances as a template that the model can recalibrate.