Forecasting Longitudinal Networks with lame

Step 1: Fit

We simulate 30 actors over 6 periods, a single dyadic covariate trade, and a directed binary outcome cooperation. The effect of trade on cooperation is time-varying and mean-reverting in the simulation: 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%). Even with many observed dyads, six periods provide limited information about temporal persistence. The chains below are short, to show the interface.

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. A sparse panel with only a few periods informs both the coefficient path and its persistence weakly. Step 5 examines the resulting uncertainty using the fitted posterior and forecast intervals.

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) # simulated coefficient path
#> [1] 0.900 0.858 0.589 0.583 0.545 0.432
sapply(Y, function(y) round(mean(y, na.rm = TRUE), 3)) # per-period density
#>    t1    t2    t3    t4    t5    t6 
#> 0.400 0.374 0.378 0.367 0.352 0.357

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
)

# estimated path of the trade coefficient
coef_path <- coef(fit)
trade_row <- grep("^trade[._]dyad$|trade", rownames(coef_path), value = TRUE)[1]
round(coef_path[trade_row, ], 3)
#>    t1    t2    t3    t4    t5    t6 
#> 0.896 0.812 0.593 0.507 0.511 0.490

Compare the estimated coefficients with the simulated path above. Step 5 also checks the upper posterior bound on $\rho_\beta$: a short panel carries limited information about the persistence even when the process mean-reverts.

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, examine the persistence draws, chain diagnostics and forecast uncertainty. Use dynamic_beta_kind = "rw1" when permanent drift is an appropriate model: RW1 fixes the persistence at one and changes the forecast law.

Step 2: Forecast `h` periods ahead

# forecasts draw the AR(1) state at random, so seed them
set.seed(1)
# 3-step-ahead forecast on the link scale
fc_link <- predict(fit, h = 3, type = "link")
length(fc_link) # 3, one matrix per future period
#> [1] 3
dim(fc_link[[1]]) # n x n (30 x 30 here)
#> [1] 30 30

# response scale: the inverse link is applied to each draw
fc_resp <- predict(fit, h = 3, type = "response")
# posterior-mean predicted probabilities
range(fc_resp[[1]], na.rm = TRUE)
#> [1] 0.009019194 0.931545467

# 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
#> [1] 30 30  3 60

# lower, median and upper matrices per period (95% by default)
fc_ci <- predict(fit, h = 3, type = "response", interval = "credible")
str(fc_ci[[1]]) # lower, median and upper, each n x n
#> List of 3
#>  $ lower : num [1:30, 1:30] 0.0538 0.034 0.0978 0.1101 0.1143 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:30] "a01" "a02" "a03" "a04" ...
#>   .. ..$ : chr [1:30] "a01" "a02" "a03" "a04" ...
#>  $ median: num [1:30, 1:30] 0.1671 0.0717 0.2443 0.1734 0.1839 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:30] "a01" "a02" "a03" "a04" ...
#>   .. ..$ : chr [1:30] "a01" "a02" "a03" "a04" ...
#>  $ upper : num [1:30, 1:30] 0.335 0.123 0.449 0.303 0.266 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:30] "a01" "a02" "a03" "a04" ...
#>   .. ..$ : chr [1:30] "a01" "a02" "a03" "a04" ...

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 (the posterior mean of rho_beta is 0.85; Step 5 reports its upper bound). 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.

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_distiller(palette = "Blues", direction = 1, 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.line = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "top",
        plot.title = element_text(hjust = 0),
        plot.subtitle = element_text(hjust = 0),
        plot.caption = element_text(hjust = 0)
    )

# right: across-dyad mean probability with the mean 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.line = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "top",
        plot.title = element_text(hjust = 0),
        plot.subtitle = element_text(hjust = 0),
        plot.caption = element_text(hjust = 0)
    )

p_heat + p_horizon

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.

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 covers the supported outcome families; the binary example above is illustrative. For family = "normal", type = "response" returns the forecasted conditional mean. For family = "poisson", each matched parameter draw has conditional mean e * exp(η + ve / 2), where e is the future-period exposure (set with newexposure). With by_draw = FALSE the point forecast is the mean count at the draw-averaged predictor and variance, e * exp(mean(η) + mean(ve) / 2), the form of the in-sample fitted() value; by_draw = TRUE returns each draw’s own conditional mean. Use interval = "prediction" for predictive count quantiles, or interval = "credible" for quantiles of the per-draw conditional means. These are distinct summaries.

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: conditional on the latent log-rate Z, the Poisson mean is period_exposure[t] * exp(Z_t). Integrating the Gaussian observation noise within a parameter draw gives period_exposure[t] * exp(η_t + ve / 2); at the posterior means of η_t and ve this is the mean count fitted() reports. Zero exposure gives a structural zero count. 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).

# seed the forecast draws
set.seed(1)
# baseline: hold trade at its last observed value
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: one more unit (one standard deviation) of trade for every dyad
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)

# change in each dyad's predicted probability at h = 3
delta_h3 <- fc_up[[3]] - fc_cf[[3]]
summary(as.vector(delta_h3))
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> 0.006484 0.095982 0.114325 0.107554 0.126730 0.141251

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 properties of 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

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

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.

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). For the fit above, that bound is 0.98; it is below the warning threshold. The simulated process mean-reverts, and the bound reports what this panel says about its persistence. The two transition models have different conditional forecast variances:

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 supports eventual mean reversion, assess AR(1) forecasts over the horizons supported by held-out validation. The warning printed by predict() states whichever of these two variance behaviours applies to your dynamic_beta_kind.

A sparse, short, drifting panel

We now simulate 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 in these four periods. The fitted persistence and forecast intervals show how the AR(1) model represents that uncertainty:

set.seed(2026)
n_sp <- 20
T_sp <- 4
beta_drift <- c(0.4, 0.9, 1.4, 1.9) # drifts upward without mean reversion
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
#>    t1    t2    t3    t4 
#> 0.034 0.055 0.118 0.161

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
)
#> ℹ This panel has 6 mutual dyads in total, so the data inform the dyadic
#>   correlation `rho` mainly through the absence of mutual ties.
#> ℹ Its posterior can stay near the prior, proportional to `1 - rho^2` on (-1,
#>   1), or move well away from it, toward -1 when the fitted means predict more
#>   mutual ties than the panel has; compare `fit$VC[, "rho"]` with that prior.
#>   `dcor = FALSE` fits a model without reciprocal residual dependence.
round(quantile(fit_sparse$RHO_BETA, c(0.5, 0.975)), 3)
#>   50% 97.5% 
#> 0.897 0.994

The upper 97.5% bound on $\rho_\beta$ is now 0.994, against the 0.99 at which predict() raises its warning, and the median sits at 0.897. The forecasts are already behaving like a unit root:

set.seed(1)
fc_sparse <- predict(
    fit_sparse,
    h = 6,
    type = "response",
    interval = "credible"
)
#> Warning: ! Posterior 95% interval for `rho_beta` reaches 0.994.
#> ℹ Under AR(1) the `h`-step forecast variance saturates at the stationary level
#>   `sigma_beta^2 / (1 - rho_beta^2)`, which itself explodes as `rho_beta -> 1`;
#>   horizons beyond "3" become uninformative in this regime.
#> ℹ Consider refitting with `dynamic_beta_kind = "rw1"` if the underlying process
#>   is genuinely unit-root.

The upper bound reaches the warning threshold. Compare how the average 95% forecast-interval width grows with the horizon on this fit versus the mean-reverting Step 1 fit:

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)
#>                 [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
#> mean_reverting 0.228 0.252 0.277 0.291 0.309 0.318
#> near_unit_root 0.240 0.274 0.332 0.298 0.329 0.327

# share of high-|sanction_risk| dyads whose h = 3 interval spans nearly [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)
#> [1] 0.45

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.line = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "top",
        plot.title = element_text(hjust = 0),
        plot.subtitle = element_text(hjust = 0),
        plot.caption = element_text(hjust = 0),
        legend.title = element_blank()
    )

Line plot comparing the across-dyad mean 95 percent forecast interval width at horizons one through six for the Step 1 fit and the sparse drifting-panel fit. These intervals are on the bounded probability scale.

The ratio of mean interval width at h = 6 to that at h = 1 is 1.4 for the Step 1 fit and 1.4 for the sparse fit. These are probability intervals, so their widths also reflect the nonlinear probit transformation and the [0, 1] bounds, so a flat probability-width curve shows the bounds as well as the coefficient uncertainty.

Among dyads in the top decile of |sanction_risk|, 45% of the sparse fit’s h = 3 intervals extend below 0.05 and above 0.95. The maximum observed period density is 16%. Those broad intervals provide little precision about individual tie probabilities. The leave-future-out scores below compare reporting horizons, and RW1 is the choice when permanent drift is the assumption.

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

# refit before each of the last two periods and score that period
# seed the forecast draws inside lfo()
set.seed(1)
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
)
#> ℹ 1 of the 2 window refits had a `rho_beta` posterior at the unit root, as
#>   short training windows often do; their one-step forecasts carry the unit-root
#>   variance.
print(lfo_res)
#> 
#> ── Exact rolling-origin LFO ──
#> 
#> Periods evaluated: 5 and 6
#> Refit per leave-out: TRUE
#> Total elpd_lfo: -960.3
#>  period      elpd n_obs
#>       5 -486.3009   870
#>       6 -473.9959   870

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). This demonstrates the output format; an analysis gives every training fit a longer chain.

Compare the two reported leave-out scores directly. A sharp decline at the last origin can suggest a late change in the data-generating process, but short chains and predictive variability can also move a score. After checking those sources of uncertainty, detect_change_point(fit) can identify a coefficient path and period for closer inspection. Its tail ratio is a heuristic reference, not a formal test or a Bayes factor.

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? The PIT uses the posterior predictive CDF, averaging the conditional distributions over matched posterior draws. For binary / Poisson / ordinal outcomes it randomises between the CDF just below and at the observation, using the Czado-Gneiting-Held PIT (Czado, Gneiting & Held 2009). A well-calibrated forecast produces PIT values that are Uniform(0, 1).

# 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
)
# seed the randomised PIT
set.seed(1)
pit <- forecast_pit(fit_train, y_future = list(Y[[6]]))
print(pit)
#> 
#> ── Forecast PIT calibration check ──
#> 
#> • Family: "binary"
#> • Forecast horizon: 1
#> • Held-out dyads: 870
#> • Kolmogorov-Smirnov statistic against Uniform(0, 1): 0.039
#> • Global calibration p-value: unavailable for dependent network cells
#> • Fraction in [0.025, 0.975]: 0.946 (expect 0.95)
#> These are descriptive diagnostics. Network dependence requires a calibrated
#> joint reference for global inference.

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

plot(pit)

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.

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 training periods. Its upper 97.5% posterior bound on $\rho_\beta$ is 0.955, which is below the forecast warning threshold; a bound at the threshold means the forecast carries the unit-root variance described in Step 5.

pit$cover_95 is the fraction of PIT values in [0.025, 0.975]; here it is 0.946, compared with 0.95 under marginal calibration. For discrete outcomes this uses the randomised PIT and is not the coverage of a discrete outcome interval. The KS statistic describes the distance between the empirical PIT distribution and the uniform reference.

pit$ks_p is NA: dependent network cells require a calibrated joint reference for a global significance test. Neither the PIT fraction nor the descriptive KS statistic supplies such a test. For binary outcomes, randomisation adds variability as well, and this vignette shows one held-out period; several held-out periods and the forecast scores above give the comparison. The seed above makes this displayed randomisation reproducible.

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 when their saved scores are valid observed likelihood factors. The following models explicitly set dcor = FALSE; estimated residual dyad correlation would require joint dyad factors. Both save the matched observed scores with save_log_lik = TRUE:

# two candidates with independent residuals and saved log-likelihoods
fit_static <- lame(
    Y,
    Xdyad = Xdyad,
    family = "binary",
    R = 0,
    dcor = FALSE,
    nscan = 5000,
    burn = 1000,
    odens = 25,
    save_log_lik = TRUE,
    verbose = FALSE
)
fit_dyn <- lame(
    Y,
    Xdyad = Xdyad,
    family = "binary",
    R = 0,
    dcor = FALSE,
    dynamic_beta = "dyad",
    nscan = 5000,
    burn = 1000,
    odens = 25,
    save_log_lik = TRUE,
    verbose = FALSE
)

# the best model has elpd_diff = 0; a named list labels the rows
cmp <- loo::loo_compare(list(
    fit_static = loo::loo(fit_static),
    fit_dyn = loo::loo(fit_dyn)
))
cmp

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. The Pareto-$k$ diagnostics in the loo output report which pointwise estimates are reliable.

These binary models store probit observed-cell likelihood factors and use independent residuals. normal, binary, poisson, and explicit-cutpoint ordinal models support observed marginal scores; Poisson scores use adaptive numerical integration with an accuracy check. A correct marginal score becomes a LOO factor only when the saved factorization contract supports it. Symmetric fits use each unordered dyad once; directed fits with estimated or nonzero residual correlation are refused by LOO/WAIC until joint dyad factors are supplied.

Capped binary, fixed-rank nomination and data-induced ordinal models do not supply an exact observed score through the current dispatcher. An explicit log_lik_method = "observed_ghk" request provides the documented legacy approximate diagnostic for those models; "augmented" provides a latent-Z diagnostic. Neither is accepted by LOO/WAIC. Use held-out evaluation appropriate to the fitted rank/cap model, rather than treating a diagnostic as its likelihood.

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 read by loo() when the same observed-factor contract passes. read_log_lik(fit) also reads marginal or diagnostic scores that cannot be used for LOO/WAIC:

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

Forecasting Longitudinal Networks with lame

Cassy Dorff, Tosin Salau, Shahryar Minhas

2026-11-04

What this vignette covers

Step 1: Fit

Step 2: Forecast `h` periods ahead

Visualising the forecast

Poisson exposure offsets via `period_exposure`

Step 3: Counterfactual forecasts

Step 4: Visualise the coefficient path

Step 5: Forecast diagnostics

A sparse, short, drifting panel

Probability-integral-transform (PIT) calibration

Comparing two fits with `loo_compare()`

Memory-conscious long runs

See also

Forecasting Longitudinal Networks with lame

Cassy Dorff, Tosin Salau, Shahryar Minhas

2026-11-04

What this vignette covers

Step 1: Fit

Step 2: Forecast h periods ahead

Visualising the forecast

Poisson exposure offsets via period_exposure

Step 3: Counterfactual forecasts

Step 4: Visualise the coefficient path

Step 5: Forecast diagnostics

A sparse, short, drifting panel

Probability-integral-transform (PIT) calibration

Comparing two fits with loo_compare()

Memory-conscious long runs

See also

Step 2: Forecast `h` periods ahead

Poisson exposure offsets via `period_exposure`

Comparing two fits with `loo_compare()`