--- 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} # the models run when the tarball is built; R CMD check re-knits only the text run_vignette <- !nzchar(Sys.getenv("_R_CHECK_PACKAGE_NAME_")) if (!run_vignette) knitr::knit_hooks$set(evaluate.inline = function(code, envir) "") knitr::opts_chunk$set( fig.retina = 1, dpi = 72, collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, eval = run_vignette ) ``` ## 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 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. ```{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) # simulated coefficient path 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 ) # 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) ``` 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 ```{r forecast} # 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 dim(fc_link[[1]]) # n x n (30 x 30 here) # 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) # 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 # 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 ``` 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 `r round(mean(fit$RHO_BETA), 2)`; 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. ```{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_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 ``` 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`). ```{r counterfactual} # 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)) ``` 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 ```{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). For the fit above, that bound is `r round(quantile(fit$RHO_BETA, 0.975), 3)`; `r if (quantile(fit$RHO_BETA, 0.975) >= 0.99) "it reaches the warning threshold" else "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: ```{r sparse-fail} 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 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)`, against the 0.99 at which `predict()` raises its warning, and the median sits at `r round(quantile(fit_sparse$RHO_BETA, 0.5), 3)`. The forecasts are already behaving like a unit root: ```{r sparse-fail-warn} set.seed(1) fc_sparse <- predict( fit_sparse, h = 6, type = "response", interval = "credible" ) ``` The upper bound `r if (quantile(fit_sparse$RHO_BETA, 0.975) >= 0.99) "reaches the warning threshold" else "is below 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: ```{r sparse-fail-width, fig.width = 6, fig.height = 3.2, fig.alt="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."} 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) # 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) 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() ) ``` The ratio of mean interval width at `h = 6` to that at `h = 1` is `r round(w_healthy[6] / w_healthy[1], 1)` for the Step 1 fit and `r round(w_sparse[6] / w_sparse[1], 1)` 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|`, `r round(100 * deg_h3)`% of the sparse fit's `h = 3` intervals extend below 0.05 and above 0.95. The maximum observed period density is `r round(100 * max(sapply(Y_sp, mean, na.rm = TRUE)))`%. 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: ```{r lfo} # 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 ) 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`). 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). ```{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 ) # seed the randomised PIT set.seed(1) 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 training periods. Its upper 97.5% posterior bound on $\rho_\beta$ is `r round(quantile(fit_train$RHO_BETA, 0.975), 3)`, which `r if (quantile(fit_train$RHO_BETA, 0.975) >= 0.99) "reaches" else "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 `r round(pit$cover_95, 3)`, 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`: ```{r loo-compare, eval = FALSE} # 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 ``` ```{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. 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: ```{r chunked, eval = FALSE} 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 ``` ## See also - The [dynamic effects article](https://netify-dev.github.io/lame/articles/dynamic_effects.html): 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.