---
title: "Forecasting Longitudinal Networks with lame"
author: "Cassy Dorff, Tosin Salau, Shahryar Minhas"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Forecasting Longitudinal Networks with lame}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
	collapse = TRUE,
	comment = "#>",
	fig.width = 7,
	fig.height = 5,
	eval = TRUE
)
```

## What this vignette covers

`lame()` with at least one dynamic component (`dynamic_beta`,
`dynamic_ab`, `dynamic_uv`) gives you a state-space model that is
exactly the right object for forecasting future periods. This vignette
walks through:

1. **Fit a model with dynamic coefficients.**
2. **Forecast `h` periods ahead** on the link scale and the response
   scale.
3. **Compute counterfactual forecasts** by feeding in alternative
   future covariates.
4. **Visualise** the forecast and the per-coefficient ribbon plot.
5. **Diagnose** whether the forecast variance has exploded (near-
   unit-root coefficients).

## Step 1: Fit

Forecasting needs a fit where the dynamic coefficient is actually
*identified* at each period, so we simulate a panel that delivers that:
30 actors over 6 periods, a single dyadic covariate `trade`, and a
directed binary outcome `cooperation`. The effect of `trade` on
cooperation is genuinely time-varying but **mean-reverting** -- it
follows an AR(1) around 0.5 with persistence 0.7, starting high (0.9)
and settling toward its long-run level. The network is moderately dense
(~37%), which is what makes per-period coefficients estimable; this is
the regime in which forecasting is well-posed.

A note on data choice: `dynamic_beta` needs both a reasonable number of
periods and enough ties per period to pin down a coefficient *for that
period*. Very sparse panels (say a rare-event sanctions network at ~2%
density over 3-4 years) do not carry that per-period information, the
AR(1) persistence runs up against a unit root, and the forecast variance
explodes (Step 5 shows exactly that failure mode). A simulated panel
lets us show forecasting working before we show it failing.

```{r fit}
library(lame)
set.seed(2026)

n_fc <- 30
T_fc <- 6
rho_true   <- 0.7
beta_bar   <- 0.5
beta_t_true <- numeric(T_fc)
beta_t_true[1] <- 0.9
for (t in 2:T_fc) {
	beta_t_true[t] <- beta_bar + rho_true * (beta_t_true[t - 1] - beta_bar) +
		rnorm(1, 0, 0.15)
}

Xdyad <- lapply(seq_len(T_fc), function(t) {
	x <- matrix(rnorm(n_fc * n_fc), n_fc, n_fc)
	array(x, dim = c(n_fc, n_fc, 1),
	      dimnames = list(NULL, NULL, "trade"))
})
a_fc <- rnorm(n_fc, 0, 0.4)
b_fc <- rnorm(n_fc, 0, 0.4)
Y <- lapply(seq_len(T_fc), function(t) {
	eta <- -0.5 + beta_t_true[t] * Xdyad[[t]][, , 1] + outer(a_fc, b_fc, "+")
	Yt  <- matrix(rbinom(n_fc * n_fc, 1, pnorm(eta)), n_fc, n_fc)
	diag(Yt) <- NA
	rownames(Yt) <- colnames(Yt) <- sprintf("a%02d", seq_len(n_fc))
	Yt
})
names(Y) <- paste0("t", seq_len(T_fc))

round(beta_t_true, 3)                                   # the truth we recover
sapply(Y, function(y) round(mean(y, na.rm = TRUE), 3))  # per-period density

fit <- lame(
	Y, Xdyad = Xdyad,
	family = "binary", R = 0,
	dynamic_beta = "dyad",      # the trade coefficient is AR(1)
	dynamic_beta_kind = "ar1",  # mean-reverting; switch to "rw1" for drift
	nscan = 150, burn = 30, odens = 5,
	verbose = FALSE
)

# the recovered per-period coefficient tracks the simulated path
coef_path <- coef(fit)
trade_row <- grep("^trade[._]dyad$|trade", rownames(coef_path), value = TRUE)[1]
round(coef_path[trade_row, ], 3)
```

The recovered trade coefficient path tracks the simulated truth, and (as we
confirm in Step 5) the posterior on $\rho_\beta$ sits comfortably below
1, so the forecast is well-posed.

A quick diagnostic: `summary(fit)` prints the per-block posterior-mean ρ_β
and fires a stationarity warning when, for any block, the 5th
percentile of the posterior on ρ_β is ≥ 0.97 **and** the IQR is < 0.1
(i.e. at least 95% of the mass sits near a unit root with little
spread). The forecast-time warning fires on a complementary trigger:
the upper 97.5% credible bound on ρ_β reaches 0.99. If either warning
fires, refit with `dynamic_beta_kind = "rw1"`. RW1 is unit-root by
construction and is the right prior for permanent-drift coefficients.

## Step 2: Forecast `h` periods ahead

```{r forecast}
set.seed(1)   # predict(h=) propagates the AR(1) state stochastically; seed for reproducibility
# 3-step-ahead forecast on the link scale (linear predictor)
fc_link <- predict(fit, h = 3, type = "link")
length(fc_link)         # 3 (one matrix per future period)
dim(fc_link[[1]])       # n x n  (30 x 30 here; unipartite)

# response-scale forecast applies the family inverse link per draw
fc_resp <- predict(fit, h = 3, type = "response")
# binary: each entry is a posterior-mean predicted probability
range(fc_resp[[1]], na.rm = TRUE)

# by_draw = TRUE returns the full [n, n, h, n_draws] array
fc_full <- predict(fit, h = 3, type = "response", by_draw = TRUE)
dim(fc_full)            # n x n x 3 x n_draws

# interval = "credible" returns a list of length-3 (lower / median / upper)
# matrices per period at the requested quantiles (default 95%)
fc_ci <- predict(fit, h = 3, type = "response", interval = "credible")
str(fc_ci[[1]])         # list($lower, $median, $upper), each n x n
```

The typical cell of `fc_ci[[1]]` carries an informative `lower` /
`median` / `upper` triple -- a genuine sub-interval of [0, 1] rather
than the whole range, because the coefficient path is mean-reverting
(`rho_beta` stays comfortably below 1). A small fraction of cells do
reach the boundary, but the vast majority stay well inside it. The
interval widens with the horizon `h` as forecast uncertainty
accumulates; on a stationary fit it rarely collapses toward the
degenerate [0, 1] you see pervasively near a unit root (Step 5 shows
that failure mode and the warning that catches it).

### Visualising the forecast

Two views make the forecast concrete. The left panel is the predicted
tie-probability matrix for the first future period (`fc_resp[[1]]`) --
the actual network the model expects next. The right panel collapses
each horizon to its across-dyad average predicted probability with the
average 95% credible band, so you can watch the band widen as the
horizon grows.

```{r forecast-plot, fig.width = 7, fig.height = 3.4, fig.alt="Left: heatmap of predicted tie probabilities among the 30 actors for the first forecast period, darker cells indicating higher predicted probability. Right: the across-dyad average predicted probability at horizons h = 1, 2, 3 with a shaded 95 percent credible band that widens with the horizon."}
library(ggplot2)
library(patchwork)

# left: forecast heatmap for the first future period
P1 <- fc_resp[[1]]
rn <- rownames(P1); if (is.null(rn)) rn <- sprintf("a%02d", seq_len(nrow(P1)))
cn <- colnames(P1); if (is.null(cn)) cn <- sprintf("a%02d", seq_len(ncol(P1)))
heat_df <- data.frame(
	sender   = factor(rn[row(P1)], levels = rev(rn)),
	receiver = factor(cn[col(P1)], levels = cn),
	prob     = as.vector(P1)
)
p_heat <- ggplot(heat_df, aes(receiver, sender, fill = prob)) +
	geom_tile() +
	scale_fill_viridis_c(limits = c(0, 1), name = "P(tie)") +
	labs(title = "Forecast for t7", x = "Receiver", y = "Sender") +
	theme_bw(base_size = 9) +
	theme(panel.border = element_blank(),
	      axis.text = element_blank(), axis.ticks = element_blank(),
	      panel.grid = element_blank(), legend.position = "top")

# right: across-dyad average predicted probability + average 95% band
horizon_df <- data.frame(
	period = factor(paste0("t", T_fc + seq_along(fc_ci)),
	                levels = paste0("t", T_fc + seq_along(fc_ci))),
	median = sapply(fc_ci, function(m) mean(m$median, na.rm = TRUE)),
	lower  = sapply(fc_ci, function(m) mean(m$lower,  na.rm = TRUE)),
	upper  = sapply(fc_ci, function(m) mean(m$upper,  na.rm = TRUE))
)
p_horizon <- ggplot(horizon_df, aes(period, median, group = 1)) +
	geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey85") +
	geom_line() + geom_point(size = 2) +
	labs(title = "Average forecast by horizon",
	     x = "Forecast period", y = "Mean P(tie)") +
	theme_bw(base_size = 9) +
	theme(panel.border = element_blank(), axis.ticks = element_blank(),
	      legend.position = "top")

p_heat + p_horizon
```

The heatmap shows the model does not forecast a uniform network: some
sender/receiver pairs carry a much higher predicted tie probability than
others, inherited from their estimated additive effects and the recovered
trade coefficient. The horizon panel shows the average predicted
probability staying roughly level (the coefficient is mean-reverting) while
its credible band widens with `h` -- the visual signature of accumulating
forecast uncertainty.

The forecast propagates `(β_t, a_t, b_t, U_t, V_t)` forward by drawing
from the joint posterior of the dynamic state-space hyperparameters
`(ρ_β, σ_β, ρ_ab, σ_ab, ρ_uv, σ_uv)` for each posterior draw. The same
forecasting machinery works for any family; the binary example above
is just for illustration. For `family = "normal"`, `type = "response"`
returns the forecasted continuous outcome; for `family = "poisson"`,
the response-scale forecast multiplies `exp(η)` by the future-period
exposure (pass via `newexposure = ...`; see the `period_exposure`
section below).

### Poisson exposure offsets via `period_exposure`

For `family = "poisson"`, pass `period_exposure = e` to `lame()` (a
non-negative numeric vector of length `T`) when the per-period rate
should be scaled by a known exposure (e.g. weeks observed, population
size, number of edge-formation opportunities). The sampler treats
`log(period_exposure[t])` as a fixed offset on the linear predictor:
the Poisson mean for period `t` becomes
`period_exposure[t] * exp(η_t)` instead of `exp(η_t)`. Setting a
non-trivial (any value not equal to 1) `period_exposure` for a
non-Poisson family is rejected with an error.
When you call `predict(fit, h = K)` later without `newexposure`, the
response-scale forecast reuses the last observed exposure for every
future period; pass `newexposure = rep(e_new, K)` to forecast under a
different exposure path.

## Step 3: Counterfactual forecasts

Pass `newdata = list_of_h_future_X_arrays` to combine alternative
future covariates with the forecast. Each array is `n x n x p` and the
slice order along the 3rd dimension must match the original
`dimnames(Xdyad[[1]])[[3]]` from the fit (here a single slice,
`trade`).

```{r counterfactual}
set.seed(1)   # forecast draws are stochastic; seed so the printed delta summary reproduces
# baseline scenario: freeze the trade covariate at the last observed period
last_X <- Xdyad[[length(Xdyad)]]
X_future <- list(last_X, last_X, last_X)
fc_cf <- predict(fit, h = 3, type = "response", newdata = X_future)

# counterfactual: shift every dyad's trade covariate up by one unit
# going forward. Because `trade` is a standardised (mean-zero) covariate,
# an additive +1 shift -- "one standard deviation more trade for every
# pair" -- is the interpretable counterfactual; a multiplicative scaling
# would be meaningless on a centred covariate (see caution 2 below).
X_future_up <- lapply(X_future, function(x) {
	x[, , 1] <- x[, , 1] + 1
	x
})
fc_up <- predict(fit, h = 3, type = "response", newdata = X_future_up)

# delta is the per-dyad change in posterior-mean cooperation probability
# at h = 3, attributable solely to the +1 trade shift
delta_h3 <- fc_up[[3]] - fc_cf[[3]]
summary(as.vector(delta_h3))
```

The trade coefficient is positive and mean-reverting, so a one-unit
increase in `trade` raises the posterior-mean cooperation probability at
`h = 3` on average (the mean and median of `delta_h3` are both positive).
The per-dyad shift varies in size, and a minority of dyads even show a
small *negative* shift: the three-step-ahead forecast of the coefficient
is itself uncertain (its posterior at `h = 3` is centred well above zero
but has nonzero mass below it), and the marginal probit effect
$\phi(\eta)\,\beta$ is smallest for dyads whose baseline linear predictor
$\eta = a_i + b_j + \beta x_{ij}$ already sits far from the 0.5
probability boundary. For substantive interpretation, summarise
`delta_h3` for the dyads that matter rather than reporting the aggregate
`summary()`.

Two cautions on counterfactual covariates. (1) Because `dynamic_beta`
draws the coefficient path forward from its posterior, the same
covariate change can produce different magnitudes at different horizons,
because the multiplier on `trade` is itself uncertain and drifting. (2)
Multiplicative counterfactuals on covariates that take negative or zero
values (such as this centred `trade` covariate) do not make substantive
sense; prefer additive shifts there, which is why we used `+ 1` rather
than `* 2` above.

## Step 4: Visualise the coefficient path

```{r autoplot, fig.width = 7, fig.height = 4, fig.alt="In-sample posterior coefficient paths from autoplot.lame: each panel shows the posterior median line and a 95 percent credible interval ribbon across the fitted training periods (t1-t6) for one regression coefficient. This plots the estimated in-sample path, not a forecast."}
library(ggplot2)
autoplot(fit, probs = c(0.025, 0.5, 0.975)) +
	labs(title = "Posterior coefficient paths",
	     subtitle = "median + 95% interval",
	     x = "Time period", y = "Coefficient value")
```

The ribbon shows the posterior 95% credible interval around the
posterior median per coefficient and period; faceted by coefficient.

## Step 5: Forecast diagnostics

`predict(fit, h = ...)` fires a one-per-fit warning if the upper
97.5% credible bound on ρ_β reaches 0.99 (i.e. the posterior puts
non-trivial mass near the unit root). On the mean-reverting fit above
that bound sits at `r round(quantile(fit$RHO_BETA, 0.975), 3)`, so
**the warning does not fire** and the forecasts are well-posed. The two regimes still have qualitatively
different forecast-variance behavior, and getting this distinction
right matters for picking your horizon:

- **AR(1), `dynamic_beta_kind = "ar1"`** with $|\rho| < 1$: the
  conditional $h$-step variance of $\beta_{T+h}$ given the training
  window is $\sigma_\beta^2 \cdot (1 - \rho^{2h}) / (1 - \rho^2)$. It
  *saturates* (the increment from $h$ to $h+1$ shrinks geometrically
  to zero) and the limit is the stationary variance
  $\sigma_\beta^2 / (1 - \rho^2)$. That stationary level itself blows
  up as $\rho \to 1$, so AR(1) is still informative at long horizons
  only when $\rho$ is comfortably bounded away from 1.
- **RW1, `dynamic_beta_kind = "rw1"`**: $\rho = 1$ by construction and
  the $h$-step variance is exactly $\sigma_\beta^2 \cdot h$, *linear*
  growth in $h$, never saturating.

In words: AR(1) eventually forgets the training window and falls back
to a stationary distribution; RW1 keeps adding innovation variance
every period. AR(1) forecasts can stay tight for moderate $h$ when
$\rho$ is small; RW1 forecasts get linearly wider every step
regardless. When the AR(1) stationarity warning fires and the
underlying process is genuinely unit-root, refit with
`dynamic_beta_kind = "rw1"`; if the AR(1) posterior on $\rho$ is
near-but-not-at 1 and the science says "should mean-revert eventually,"
stick with AR(1) but cap your reported horizon at $h \approx 3$. The
warning printed by `predict()` states whichever of these two variance
behaviours applies to your `dynamic_beta_kind`.

### Watching it fail: a sparse, short, drifting panel

Everything above describes the failure mode; here it is actually
happening. We simulate the kind of panel the Step 1 note warned
about: fewer actors, only four periods, single-digit tie density in
the early periods, and a covariate effect that *drifts upward*
instead of mean-reverting. There is no long-run level for the AR(1)
to find, so the posterior on $\rho_\beta$ runs up against the unit
root:

```{r sparse-fail}
set.seed(2026)
n_sp <- 20
T_sp <- 4
beta_drift <- c(0.4, 0.9, 1.4, 1.9)   # drifts up; never mean-reverts
X_sp <- lapply(seq_len(T_sp), function(t) {
	x <- matrix(rnorm(n_sp * n_sp), n_sp, n_sp)
	array(x, dim = c(n_sp, n_sp, 1),
	      dimnames = list(NULL, NULL, "sanction_risk"))
})
Y_sp <- lapply(seq_len(T_sp), function(t) {
	eta <- -2 + beta_drift[t] * X_sp[[t]][, , 1]
	Yt  <- matrix(rbinom(n_sp * n_sp, 1, pnorm(eta)), n_sp, n_sp)
	diag(Yt) <- NA
	rownames(Yt) <- colnames(Yt) <- sprintf("s%02d", seq_len(n_sp))
	Yt
})
names(Y_sp) <- paste0("t", seq_len(T_sp))
sapply(Y_sp, function(y) round(mean(y, na.rm = TRUE), 3))  # sparse early on

fit_sparse <- lame(
	Y_sp, Xdyad = X_sp,
	family = "binary", R = 0,
	dynamic_beta = "dyad", dynamic_beta_kind = "ar1",
	nscan = 400, burn = 100, odens = 5,
	verbose = FALSE
)
round(quantile(fit_sparse$RHO_BETA, c(0.5, 0.975)), 3)
```

The upper 97.5% bound on $\rho_\beta$ is now
`r round(quantile(fit_sparse$RHO_BETA, 0.975), 3)` -- past the 0.99
threshold -- so the forecast call itself raises the warning:

```{r sparse-fail-warn}
set.seed(1)
fc_sparse <- predict(fit_sparse, h = 6, type = "response",
                     interval = "credible")
```

That is the diagnostic doing its job at forecast time. To see what the
near-unit-root posterior does to the intervals, compare how the average
95% forecast-interval width grows with the horizon on this fit versus
the mean-reverting Step 1 fit:

```{r sparse-fail-width, fig.width = 6, fig.height = 3.2, fig.alt="Line plot of the across-dyad mean 95 percent forecast interval width at horizons one through six for two fits. The mean-reverting Step 1 fit's width is flat from roughly horizon four on, the saturation an AR(1) below the unit root predicts, while the near-unit-root sparse fit's width is still rising at horizon six, roughly double its one-step value."}
set.seed(1)
fc_healthy <- predict(fit, h = 6, type = "response",
                      interval = "credible")
width_by_h <- function(fc) {
	sapply(fc, function(m) mean(m$upper - m$lower, na.rm = TRUE))
}
w_healthy <- width_by_h(fc_healthy)
w_sparse  <- width_by_h(fc_sparse)
round(rbind(mean_reverting = w_healthy, near_unit_root = w_sparse), 3)

# at h = 3, how many of the dyads the covariate actually moves (top
# decile of |sanction_risk|) get an interval spanning nearly all of
# [0, 1]?
x_last <- X_sp[[T_sp]][, , 1]
hi_x   <- abs(x_last) >= quantile(abs(x_last), 0.9)
deg_h3 <- mean(fc_sparse[[3]]$lower[hi_x] < 0.05 &
               fc_sparse[[3]]$upper[hi_x] > 0.95, na.rm = TRUE)
round(deg_h3, 2)

width_df <- data.frame(
	h     = rep(seq_along(w_healthy), 2),
	fit   = rep(c("mean-reverting (Step 1 fit)",
	              "near-unit-root (sparse fit)"), each = length(w_healthy)),
	width = c(w_healthy, w_sparse)
)
ggplot(width_df, aes(h, width, linetype = fit)) +
	geom_line() + geom_point(size = 1.8) +
	labs(title = "Forecast interval width by horizon",
	     x = "Forecast horizon h", y = "Mean 95% interval width") +
	theme_bw(base_size = 9) +
	theme(panel.border = element_blank(), axis.ticks = element_blank(),
	      legend.position = "top", legend.title = element_blank())
```

Read the growth, not the levels: the sparse fit's absolute widths are
smaller only because most of its predicted probabilities sit near zero
(the network is sparse). The mean-reverting fit's width grows by a
factor of just `r round(w_healthy[6] / w_healthy[1], 1)` across the
whole horizon and is flat from about `h = 4` on -- the saturation the
AR(1) bullet above predicts. The near-unit-root fit's width has grown
by a factor of `r round(w_sparse[6] / w_sparse[1], 1)` by `h = 6` and
is still climbing -- the growth you get when $\rho_\beta$ sits at the
unit root, damped here only by the bounded [0, 1] probability scale. And the degeneracy lands exactly where it hurts: among
the dyads the covariate actually moves (the top decile of
`|sanction_risk|`), `r round(100 * deg_h3)`% of the `h = 3` intervals
already span essentially the whole unit interval, on a network whose
observed density never exceeded
`r round(100 * max(sapply(Y_sp, mean, na.rm = TRUE)))`%. A forecast
interval of [0, 1] for a rare-event dyad says nothing at all -- which is
precisely what the warning is telling you. The remedy is the one stated
above: refit with `dynamic_beta_kind = "rw1"` if the drift is real, and
in either case do not report horizons past $h \approx 3$.

For a longer-horizon diagnostic, use the exact rolling-origin
leave-future-out CV:

```{r lfo}
# refit on Y[1:(t-1)] for each t in periods and score Y[[t]].
# Skip the very first leave-out (would leave a 1-period training window
# which is too short for dynamic_beta); use the last 2 origins.
set.seed(1)  # the internal h=1 forecast draws inside lfo() are stochastic; seed for reproducibility
T_fit <- length(fit$YPM)
lfo_periods <- tail(seq_len(T_fit), max(1L, min(2L, T_fit - 2L)))
lfo_res <- lfo(fit, periods = lfo_periods, refit = TRUE,
			   nscan = 100, burn = 25, odens = 5, verbose = FALSE)
print(lfo_res)
```

The `nscan = 100, burn = 25, odens = 5` settings inside `lfo()` are
chosen so the vignette builds in well under a minute on a laptop. They
give `100 / 5 = 20` stored draws per leave-out refit (burn-in is run and discarded, not subtracted from `nscan`), enough to
return a sensible point estimate of per-period elpd for *demonstration*,
but tighter than you should use when the LFO output drives an inference.
For real analyses, pass `nscan = 5000, burn = 1000, odens = 25` (or
whatever matches your `lame()` baseline) and verify the per-period elpd
is stable across two independent runs before relying on it.

On this mean-reverting fit the two leave-out folds return per-period
elpd of roughly -560 to -620 over 870 scored dyads (about -0.65 to -0.7
nats/dyad), and the two folds are close to each other -- the model
forecasts the held-out period about as well at origin 5 as at origin 6,
which is what you want to see when there is no late regime change. If `elpd` instead
dropped sharply at the last fold, that would be the empirical signal of
a regime change near the end of your training window; reach for
`detect_change_point(fit)` (a heuristic Bayes-factor diagnostic) to
localise it.

### Probability-integral-transform (PIT) calibration

`forecast_pit()` complements `lfo()` by asking a sharper question:
*given the held-out outcomes, is the h-step posterior-predictive
distribution itself well calibrated?* For continuous families the PIT
is the Gaussian CDF evaluated at the observation; for binary /
Poisson / ordinal it is the Czado-Gneiting-Held randomised PIT
(Czado, Gneiting & Held 2009). A well-calibrated forecast produces
PIT values that are Uniform(0, 1).

```{r pit}
# fit on the first five periods, hold out the sixth
fit_train <- lame(
	Y[1:5], Xdyad = Xdyad[1:5],
	family = "binary", R = 0,
	dynamic_beta = "dyad",
	nscan = 100, burn = 30, odens = 5,
	verbose = FALSE
)
set.seed(1)   # the randomised PIT draws are stochastic; seed for reproducibility
pit <- forecast_pit(fit_train, y_future = list(Y[[6]]))
print(pit)
```

`plot(pit)` renders the PIT values against the Uniform(0, 1) reference
that a perfectly calibrated forecast would follow:

```{r pit-plot, fig.height = 3.6, fig.alt="PIT calibration plot for the held-out period: the empirical distribution of probability-integral-transform values compared with the Uniform(0,1) reference. Bars or a curve tracking the reference indicate calibration; a U-shape signals an over-confident (too-narrow) forecast and a central hump signals an under-confident (too-wide) one."}
plot(pit)
```

Bars that hug the uniform reference indicate a well-calibrated forecast.
A U-shape (mass piling up in the tails) is the signature of an
over-confident forecast whose intervals are too narrow; a central hump
is the opposite, an under-confident forecast with intervals that are too
wide. On this single binary held-out period the randomised PIT is noisy,
so read the plot for gross departures rather than small wiggles, and
corroborate it with `cover_95` below.

(`fit_train` has only five periods to pin down $\rho_\beta$, so its
posterior is wider than that of the six-period `fit` above; on very
short or sparse panels that extra uncertainty can spill past the
forecast-time unit-root threshold, exactly as the drifting sparse
panel earlier in Step 5 demonstrated. Here it stays clear of the
threshold -- the unit-root warning does not fire for `fit_train` and
the `print(pit)` output above is clean. If you *do* see that warning
on your own shorter-window fit, read it as a caution about the
training fit's forecast variance rather than a defect in the PIT
check -- the calibration summaries below remain meaningful either
way.)

Read the two summaries together. `pit$cover_95` -- the fraction of
held-out cells whose observed outcome falls in the central 95%
posterior-predictive interval, target 0.95 -- is the **stable** summary;
here it lands at about 0.94, close to nominal, so the forecast intervals
have roughly the right width. `pit$ks_p` tests whether the full PIT
distribution is Uniform(0, 1); a value below 0.05 flags
mis-calibration. For a **binary** outcome the PIT is the randomised
Czado-Gneiting-Held version, so on a single held-out period `ks_p`
carries randomisation noise on top of ordinary sampling noise. Here it
is `r round(pit$ks_p, 3)` -- above 0.05, so this run passes -- but
redrawing the randomisation alone can move that p-value from below 0.05
to far above it, so a single KS p-value on one held-out period neither
confirms nor condemns the forecast. (We seed the chunk only so the
printed value reproduces.) Treat `cover_95` as the headline and `ks_p`
as corroborating evidence. For a sharper test, hold out two or
more periods and pair `plot(pit)` with `pit$pit` to see whether any
mis-calibration is in the tails (over- vs under-dispersion) or the
centre (location bias). On a sparse, short, near-unit-root panel both
summaries degrade together -- `cover_95` drops well below 0.95 and
`ks_p` collapses toward 0 -- the calibration-side signature of the
forecast-variance explosion the Step 5 warning flags.

## Comparing two fits with `loo_compare()`

When you have two candidate specifications, say a static-`beta`
baseline and a `dynamic_beta = "dyad"` variant, `loo::loo_compare()`
is the one-liner that ranks them on out-of-sample predictive
performance. Both fits need `save_log_lik = TRUE` so `loo()` can read
the stored pointwise log-likelihood matrix off `fit$log_lik`:

```{r loo-compare, eval = FALSE}
# run this after fitting both candidates with converged chains
fit_static <- lame(
	Y, Xdyad = Xdyad, family = "binary", R = 0,
	nscan = 5000, burn = 1000, odens = 25,
	save_log_lik = TRUE,
	verbose = FALSE
)
fit_dyn <- lame(
	Y, Xdyad = Xdyad, family = "binary", R = 0,
	dynamic_beta = "dyad",
	nscan = 5000, burn = 1000, odens = 25,
	save_log_lik = TRUE,
	verbose = FALSE
)

# rank the two fits: the row with elpd_diff = 0 is the best,
# subsequent rows show elpd_diff and its SE relative to it. Pass a NAMED
# list so the rows are labelled by model name rather than "model1"/"model2".
cmp <- loo::loo_compare(list(fit_static = loo::loo(fit_static),
                             fit_dyn    = loo::loo(fit_dyn)))
cmp
```

```{r loo-compare-read, eval = FALSE}
best     <- rownames(cmp)[1]
runnerup <- rownames(cmp)[2]
diff_val <- abs(round(cmp[2, "elpd_diff"], 1))
se_val   <- round(cmp[2, "se_diff"], 1)
ratio    <- round(diff_val / se_val, 1)
```

**How to read this output.** The top row is the preferred model and has
`elpd_diff = 0`. Each subsequent row gives the difference in expected log
pointwise predictive density relative to that model and the standard error of
the difference. A difference that is large relative to its standard error
supports the preferred model; when the two are similar, the specifications are
not distinguishable on this predictive criterion. Read the Pareto-$k$
diagnostics before interpreting the ranking.

Because `family = "binary"` is one of the five families with an exact
closed-form log-likelihood (`normal`, `binary`, `cbin`, `poisson`,
`ordinal`; see the **log-lik scale** discussion in the
[Dynamic Effects vignette](dynamic_effects.html#loo-waic-via-save_log_lik-true)),
the two `elpd_loo` values here are on the response scale and the
absolute numbers are directly comparable to a `loo()` result from a
probit-link `brms` fit on the same data, modulo the `lame` Gibbs
sampler's R-hat / ESS being satisfactory. For the rank likelihood
`frn`, the default `observed_exact` falls back to an
augmented-Z normal approximation and is only valid for *relative*
comparison within `lame`; opt in to `log_lik_method = "observed_ghk"`
for the exact marginal on that family.

### Memory-conscious long runs

If the in-memory log-likelihood matrix from `save_log_lik = TRUE`
would exceed RAM, pass `save_log_lik = "chunked"` instead. The
on-disk chunks are recovered transparently by `loo()`, or can be
read back directly as a matrix with `read_log_lik(fit)`:

```{r chunked, eval = FALSE}
fit_chk <- lame(
	Y, Xdyad = Xdyad, family = "binary", R = 0,
	dynamic_beta = "dyad",
	nscan = 2000, burn = 500, odens = 5,
	save_log_lik = "chunked",          # per-column-chunk binary files
	log_lik_chunk_size = 5000L,
	verbose = FALSE
)
loo::loo(fit_chk)                       # reads the chunks transparently
```

## See also

- `vignettes/dynamic_effects.Rmd`: the full reference on
  `dynamic_beta` / `dynamic_ab` / `dynamic_uv` and the decision tree
  for picking `dynamic_beta_kind`.
- `?predict.lame`: argument reference for `h`, `newdata`, `type`,
  `by_draw`.
- `?autoplot.lame`: ggplot ribbon plot for dynamic coefficients.
- `?lfo`: exact rolling-origin LFO CV.
