--- title: "Dynamic Effects in Longitudinal AME Models" author: "Cassy Dorff, Tosin Salau, Shahryar Minhas" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Dynamic Effects in Longitudinal AME Models} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) ``` ## Introduction Networks change. Allies drift apart, tastes evolve. The standard AME model pools across time, assuming each actor has a fixed latent position and fixed additive effects in every period. The `lame` package provides three mechanisms for relaxing that assumption: `dynamic_uv` for time-varying latent positions, `dynamic_ab` for time-varying additive effects, and `dynamic_beta` for time-varying regression coefficients. This vignette explains what each option does, when to use it, and how to interpret the results. ## What Are Dynamic Effects? ### The Static Baseline In the standard AME model for longitudinal data, the tie between actors $i$ and $j$ at time $t$ is: $$y_{ij,t} = \beta'x_{ij,t} + a_i + b_j + u_i'v_j + \epsilon_{ij,t}$$ The covariates ($x_{ij,t}$) can vary over time, but the sender effect $a_i$, the receiver effect $b_j$, and the latent positions $u_i$ and $v_j$ are constant: the model assumes a country's position in the latent "sanctioning space" is the same in 1993 as in 2000. ### Dynamic Latent Positions (`dynamic_uv = TRUE`) When you set `dynamic_uv = TRUE`, each actor's latent position evolves over time according to an AR(1) process: $$U_{i,k,t} = \rho_{uv} \, U_{i,k,t-1} + \epsilon_{i,k,t}, \quad \epsilon_{i,k,t} \sim N(0, \sigma_{uv}^2)$$ $$V_{j,k,t} = \rho_{uv} \, V_{j,k,t-1} + \eta_{j,k,t}, \quad \eta_{j,k,t} \sim N(0, \sigma_{uv}^2)$$ Note the **innovation parameterisation**: $\sigma_{uv}^2$ is the variance of the period-to-period innovation, not the marginal/stationary variance of $U_{i,k,t}$ (which is $\sigma_{uv}^2 / (1 - \rho_{uv}^2)$), and `fit$sigma_uv` holds posterior draws of the innovation **standard deviation**. The autoregressive parameter $\rho_{uv}$, estimated from the data, controls persistence: near 1 means positions change slowly, near 0 means they are essentially re-drawn each period. Use this when the underlying community structure is evolving -- alliances shift, trading blocs realign. ### Dynamic Additive Effects (`dynamic_ab = TRUE`) When you set `dynamic_ab = TRUE`, the sender and receiver effects evolve the same way: $$a_{i,t} = \rho_{ab} \, a_{i,t-1} + \epsilon_{i,t}, \quad \epsilon_{i,t} \sim N(0, \sigma_{ab}^2)$$ $$b_{j,t} = \rho_{ab} \, b_{j,t-1} + \eta_{j,t}, \quad \eta_{j,t} \sim N(0, \sigma_{ab}^2)$$ The same innovation parameterisation applies: `fit$sigma_ab` holds draws of the innovation standard deviation, and the implied stationary variance $\sigma_{ab}^2 / (1 - \rho_{ab}^2)$ doubles as the prior on the initial state $a_{i,1}$. Substantively, `dynamic_ab` captures changes in how much an actor participates, not in who they connect with (that is `dynamic_uv`'s territory). The two can be used alone or together; `dynamic_ab` is often the place to start, since changes in overall activity are common and relatively easy to estimate. ### Dynamic Regression Coefficients (`dynamic_beta`) The third option lets the regression coefficients themselves evolve. With `dynamic_beta = TRUE` (or a subset like `dynamic_beta = "dyad"`), each selected coefficient follows an AR(1): $$\beta_{k,t} = \rho_{\beta,b(k)} \, \beta_{k,t-1} + \epsilon_{k,t}, \quad \epsilon_{k,t} \sim N(0, \sigma_{\beta,b(k)}^2 \Lambda_{kk})$$ where $b(k)$ is the block of coefficient $k$ (intercept / dyad / row / col), so coefficients in a block share $\rho$ and $\sigma$, while the design-based scale $\Lambda_{kk}$ (a $g$-prior analogue) keeps each innovation variance commensurate with the data. The whole time path of $\beta$ is drawn jointly via forward-filter / backward-sample (FFBS; Carter & Kohn 1994). One storage difference from `dynamic_uv` / `dynamic_ab`: `fit$sigma_beta` and `fit$rho_beta` hold last-iteration scalars, while the full posterior draws live in `fit$SIGMA_BETA` and `fit$RHO_BETA` -- summarise those (e.g. `colMeans(fit$SIGMA_BETA)`). `dynamic_beta` accepts `FALSE` (default), `TRUE` (all coefficients), block shortcuts (`"intercept"`, `"dyad"`, `"row"`, `"col"`) or coefficient names, integer column indices, or a length-`p` logical mask; see `?lame`. When the intercept (or a nodal coefficient) is dynamic, an automatic sum-to-zero constraint on `a` (and `b`) keeps everything jointly identified. `dynamic_beta` works for every family, unipartite and bipartite, directed and symmetric, and combines with `dynamic_ab` / `dynamic_uv`. #### A worked `dynamic_beta` example Simulate a small longitudinal network where a single dyadic covariate's effect grows linearly over time, and fit with `dynamic_beta = "dyad"`. ```{r dynamic-beta-example, message=FALSE} library(lame) set.seed(2026) n_db <- 15 T_db <- 5 beta_t_true <- seq(-0.5, 1.0, length.out = T_db) # truth: rising effect X_db <- replicate(T_db, matrix(rnorm(n_db * n_db), n_db, n_db), simplify = FALSE) Y_db <- vector("list", T_db) for (t in seq_len(T_db)) { Yt <- beta_t_true[t] * X_db[[t]] + matrix(rnorm(n_db * n_db, 0, 0.4), n_db, n_db) diag(Yt) <- NA rownames(Yt) <- colnames(Yt) <- paste0("a", seq_len(n_db)) Y_db[[t]] <- Yt } names(Y_db) <- paste0("t", seq_len(T_db)) X_db_arr <- lapply(X_db, function(x) array(x, c(n_db, n_db, 1))) fit_db <- lame(Y_db, Xdyad = X_db_arr, family = "normal", R = 0, nscan = 200, burn = 50, odens = 5, dynamic_beta = "dyad", verbose = FALSE) # fit$BETA is a 3-D array [n_stored x p x T] dim(fit_db$BETA) # coef() returns a p x T matrix of posterior means coef(fit_db) # confint() returns one row per coef[t] with 95% credible interval head(confint(fit_db), 8) ``` **Migration warning for `amen` / static-`lame` users.** With `dynamic_beta` active, `fit$BETA` becomes a 3-D `[n_stored, p, T]` array rather than the 2-D `[n_stored, p]` matrix. The amen-era idiom `apply(fit$BETA, 2, mean)` still returns a length-`p` vector -- but it now silently averages across both iterations **and** time, hiding the time variation you fit the model to recover. Use `coef(fit)` instead: it returns a `[p, T]` matrix for dynamic fits and a length-`p` vector for static ones. The per-period posterior mean for `X1_dyad` should track the linear ramp `beta_t_true = c(-0.5, -0.125, 0.25, 0.625, 1.0)`, with the AR(1) prior providing mild shrinkage relative to fitting a separate `ame()` per period. #### Choosing a prior family: `dynamic_beta_kind` `lame()` exposes `dynamic_beta_kind` for selecting the prior family on the time-varying coefficients: > **Mean-reverting** coefficient? Use the default `dynamic_beta_kind = "ar1"`. > **Permanent / unit-root drift** (no mean reversion)? Use `"rw1"`. > **Smooth with curvature** (second-order smoothness)? Use `"rw2"`. > **Smooth with a known length-scale** of variation? Use `"matern32"` (set the scale via `prior = list(matern32_length_scale = ...)`; defaults to `max(2, T/4)`). > If the stationarity warning fires after fitting AR(1) -- the posterior on $\rho_\beta$ concentrated at or above 0.97 -- refit with `"rw1"` (or `"rw2"` if acceleration matters). `"random_walk"` is an alias for `"rw1"`; the `"rw2"` and `"matern32"` samplers use a slower R-level joint Gaussian update rather than the C++ FFBS. Refitting `fit_db` under each kind and comparing the recovered `X1_dyad` path against the truth: ```{r dynamic-beta-kinds} kinds <- c("rw1", "rw2", "matern32") paths <- sapply(kinds, function(k) { fit_k <- lame(Y_db, Xdyad = X_db_arr, family = "normal", R = 0, nscan = 200, burn = 50, odens = 5, dynamic_beta = "dyad", dynamic_beta_kind = k, prior = if (k == "matern32") list(matern32_length_scale = 2) else NULL, verbose = FALSE) coef(fit_k)["X1_dyad", ] }) round(rbind(ar1 = coef(fit_db)["X1_dyad", ], t(paths), truth = beta_t_true), 3) ``` All four kinds recover essentially the same ramp -- the point to take away, not a tie-breaker: on a smooth, gap-free trend the priors are nearly indistinguishable in-sample, and they diverge only when the data match their assumption. Choose with the decision tree above and the post-fit stationarity warning, not by eyeballing this table. When `time_index` is supplied with unequal gaps, the AR(1) / RW1 conditional variance is scaled by the gap (suppressing $\Lambda_{kk}$ for readability): $$q_t = \sigma_\beta^2 \cdot \frac{1 - \rho_\beta^{2\Delta_t}}{1 - \rho_\beta^2}$$ so quarterly-then-annual data gets the correct prior; `"rw2"` and `"matern32"` support unequal gaps through their joint precision-matrix construction. ```{r dynamic-beta-time-index} # observations at t = 1, 2, 4, 8 (a doubling gap structure) fit_gaps <- lame(Y_db, Xdyad = X_db_arr, family = "normal", R = 0, nscan = 200, burn = 50, odens = 5, dynamic_beta = "dyad", time_index = c(1, 2, 4, 8, 16)[seq_len(T_db)], verbose = FALSE) fit_gaps$time_index round(coef(fit_gaps), 3) ``` #### Hierarchical pooling: `dynamic_beta_pool` When two or more coefficient blocks are dynamic, each block has its own AR(1) hyperparameters by default, and short panels ($T \le 5$) can leave them badly under-identified. `dynamic_beta_pool` shrinks block-level hyperparameters toward a shared mean: `"rho"` for common persistence, `"sigma"` for a common innovation scale, `"both"` for full pooling -- e.g. `lame(..., dynamic_beta = c("intercept", "dyad"), dynamic_beta_pool = "both")`. With a single dynamic block the hierarchical layer is skipped entirely and the fit is identical to the default `"none"`. See `?lame` for the diagnostic comparison of pooled vs independent hyperparameter posteriors. #### Per-actor time-varying slopes: `dynamic_beta_per_actor` For questions about whether one actor responds to a covariate differently than another -- and whether that gap evolves -- set `dynamic_beta_per_actor = "row"` (or `"col"`) and pick the covariate via `per_actor_covariate_idx`, e.g. `lame(..., dynamic_beta_per_actor = "row", per_actor_covariate_idx = 1L)`. Deviations are sum-to-zero centred at every period, so the population coefficient stays in `coef(fit)` while the per-actor deviations land in `fit$THETA_ACTOR` (and their posterior means in `fit$theta_actor_mean`, labelled by actor and period). See `?per_actor_slopes`. #### Diagnosing your dynamic fit `summary(fit)` prints the *Dynamic coefficients per period* block and the per-block posterior-mean `rho_beta`, and warns when a block's posterior on $\rho_\beta$ is concentrated at or above 0.97 -- the signal to refit with `dynamic_beta_kind = "rw1"`. `detect_change_point(fit, threshold_bf = 5)` is a heuristic regime-switch diagnostic: it scales each coefficient's largest one-period jump by $\sigma_\beta$ and compares it to a prior-null reference, reporting a bounded tail-ratio score in the `bf` column (not a formal Bayes factor). On the smooth linear ramp above there is no break, so we expect a score of 0: ```{r change-point-demo} detect_change_point(fit_db, threshold_bf = 5) ``` Read `bf` below 3 as weak, 3 to 10 as worth inspecting, and above 10 as a jump the AR(1) prior struggles to generate; `t_hat` is a prompt for visual inspection, not a formal break-time estimate. The diagnostic is deliberately conservative (only stark breaks trip it reliably at short $T$); see `?detect_change_point`. For prior elicitation *before* fitting, `dynamic_beta_prior_summary(T = 10, kind = "ar1")` simulates prior paths so you can check the implied roughness and range for your `T`. Finally, the Gibbs analog of HMC divergences: every fit exposes per-block Metropolis failure counts as `fit$mh_counters`, and a warning fires automatically when any block exceeds a 10% failure rate. Treat anything above 5% as suspicious. #### Temporal-trend posterior-predictive check: `gof_temporal()` Where `detect_change_point()` targets *coefficients*, `gof_temporal(fit, stat = "density", n_rep = 200)` targets the *network statistic itself*: it fits a linear time trend to a chosen statistic (density / mean / reciprocity / transitivity), computes the same slope on `simulate(fit)` replicates, and returns a two-sided posterior-predictive p-value; small `p_pp` (below 0.05) means the observed trend is incompatible with the fitted model. Reach for it when the substantive statistic drifts across the panel; the informative worked example lives in the [overview vignette](lame-overview.html); see `?gof_temporal`. #### Inspecting the priors actually in effect: `prior_summary()` `prior_summary(fit)` prints the hyperparameters actually used, defaults filled in -- for a dynamic fit that includes the per-block AR(1) priors alongside the variance components and the g-prior on `beta`. Run it once per fit to catch typos in a `prior = list(...)` override; see `?prior_summary` and the [cross-sectional vignette](cross_sec_ame.html) for the general tour. #### Multi-chain $\hat R$ via `lame_parallel()` and `rhat_dynamic_beta()` A single chain only ever buys you within-chain split-$\hat R$. For a between-chain diagnostic on the per-period coefficient path, fit several chains with `lame_parallel(..., n_chains = 4, combine_method = "list")` and call `rhat_dynamic_beta()` on the resulting list. Require both `rhat_mvt < 1.01` (the multivariate statistic catches chains that agree per period but disagree on the trajectory) and `rhat_max_univariate < 1.01` before quoting a path; short vignette-length chains typically fail this bar, so use production lengths. `combine_ame_chains()` pools chains and routes through `posterior::summarise_draws()` for per-(coefficient, period) ESS. `ame_parallel()` is the cross-sectional twin; see `?rhat_dynamic_beta`, `?lame_parallel`, and the [cross-sectional vignette](cross_sec_ame.html). ### Forecasting When at least one dynamic component is active, `predict(fit, h = K)` propagates the state-space model forward $K$ periods, sampling $(\rho_\beta, \sigma_\beta)$ from the posterior so the forecast carries hyperparameter uncertainty. The [forecasting vignette](forecasting.html) is the full worked reference (forecast scales, per-draw arrays, credible intervals, counterfactual covariates, exposure offsets, near-unit-root warning). One kind-specific note lives here: under `"rw2"` / `"matern32"` the forward recursion pins $\rho_\beta$ at 1, so forecast variance grows linearly in $h$ as under RW1 -- the smoothing structure shapes the *in-sample* path only. ### Visualising the coefficient path: `autoplot.lame` For a fit with `dynamic_beta` active, `autoplot(fit)` returns a `ggplot` ribbon plot per coefficient (posterior median line, 95% interval, faceted by coefficient). Under `dynamic_beta = "dyad"` the intercept remains static, so pass `coefs =` to restrict the plot to the genuinely time-varying coefficients: ```{r autoplot-example, fig.width = 6, fig.height = 4, fig.alt="Ribbon plot of the posterior path of the dyadic coefficient X1_dyad across five time periods, with the median as a line and the 95 percent credible interval as a shaded band, rising steadily from about minus 0.5 to 1."} library(ggplot2) autoplot(fit_db, coefs = "X1_dyad", probs = c(0.025, 0.5, 0.975)) + labs(title = "Posterior coefficient path", subtitle = "Median line, 95% credible interval ribbon", x = "Time period", y = "Coefficient value") ``` ### LOO / WAIC via `save_log_lik = TRUE` Pass `save_log_lik = TRUE` and `lame` attaches a `[n_stored, n_obs]` pointwise log-likelihood matrix to `fit$log_lik`, which registered S3 methods hand directly to `loo::loo(fit)` and `loo::waic(fit)`: ```{r save-log-lik-eval, eval = FALSE} fit_ll <- lame(Y_db, Xdyad = X_db_arr, family = "normal", R = 0, nscan = 2000, burn = 500, odens = 25, dynamic_beta = "dyad", save_log_lik = TRUE, verbose = FALSE) dim(fit_ll$log_lik) # [n_stored, n_observations] loo_db <- loo::loo(fit_ll) print(loo_db) ``` The code is left unevaluated because a useful PSIS-LOO calculation needs more draws than a vignette should generate. Run it with a converged fit before using `elpd_loo`. If the matrix would not fit in RAM, `save_log_lik = "chunked"` streams it to disk (see the [forecasting vignette](forecasting.html)). On the log-likelihood scale: for `normal`, `binary`, `cbin`, `poisson`, and `ordinal` the stored pointwise log-likelihood is the **exact family-specific Y density** on the response scale, so `elpd_loo` is directly comparable to a Stan / brms `loo()` on the same family and data (modulo satisfactory $\hat R$ / ESS). For the rank likelihood `frn`, the exact marginal needs GHK Monte Carlo (`log_lik_method = "observed_ghk"`, longitudinal path only); the default falls back to the augmented-Z normal density with a one-time warning. Inspect `fit$log_lik_method` for the branch used. ### Prior Specifications The AR(1) coefficients and innovation variances have sensible default priors: - $\rho_{uv} \sim \text{TruncNormal}(0.9, 0.1, 0, 1)$, centered on high persistence - $\rho_{ab} \sim \text{TruncNormal}(0.8, 0.15, 0, 1)$, slightly less persistent - $\rho_{\beta} \sim \text{TruncNormal}(0.8, 0.15, 0, 0.999)$ - $\sigma_{uv}^2, \sigma_{ab}^2, \sigma_{\beta}^2 \sim \text{InverseGamma}(2, 1)$ All can be customized via the `prior` argument (e.g. `prior = list(rho_beta_mean = ..., sigma_beta_scale = ...)`) if you have strong beliefs about the rate of change. ## Genuine Dynamics: What They Look Like We start with an example where dynamics are truly present: 25 actors over 5 periods whose latent positions evolve via an AR(1) with $\rho = 0.85$ and innovation SD $\sigma_{uv} = 0.3$. Initial positions are drawn at SD 1.5 -- above the implied stationary SD of about 0.57 -- seeding strong latent structure that drifts toward stationarity, with a negative intercept giving realistic density (roughly 30%). ```{r, message=FALSE} library(lame) set.seed(6886) n_dyn <- 25 n_per <- 5 R_dyn <- 2 # evolving latent positions with AR(1) dynamics U <- matrix(rnorm(n_dyn * R_dyn, 0, 1.5), n_dyn, R_dyn) V <- matrix(rnorm(n_dyn * R_dyn, 0, 1.5), n_dyn, R_dyn) Y_dyn <- list() for(t in 1:n_per) { if(t > 1) { U <- 0.85 * U + matrix(rnorm(n_dyn * R_dyn, 0, 0.3), n_dyn, R_dyn) V <- 0.85 * V + matrix(rnorm(n_dyn * R_dyn, 0, 0.3), n_dyn, R_dyn) } eta <- -1 + U %*% t(V) Y_t <- matrix(rbinom(n_dyn * n_dyn, 1, pnorm(eta)), n_dyn, n_dyn) diag(Y_t) <- NA rownames(Y_t) <- colnames(Y_t) <- paste0("A", 1:n_dyn) Y_dyn[[t]] <- Y_t } # iterations are small so the vignette builds quickly # use burn >= 1000 and nscan >= 5000 for real analyses. fit_dyn_real <- lame(Y_dyn, R = 2, dynamic_uv = TRUE, dynamic_ab = TRUE, family = "binary", burn = 100, nscan = 1000, odens = 10, verbose = FALSE, plot = FALSE) summary(fit_dyn_real) ``` The trajectory plot shows how actors move through the latent space. With more than eight actors the legend becomes an unreadable rainbow, so pass `highlight =` a few actor names to grey out the rest (colour-blind-safe Okabe--Ito palette for the highlighted set). With genuine temporal structure, the highlighted paths trace short, coherent arcs rather than random jumps: ```{r uv-traj-real, fig.width=7, fig.height=5, fig.alt="Trajectory plot of actors' paths through the 2D latent space across five periods under genuine AR(1) dynamics; four named actors (A1, A5, A10, A15) are coloured with the Okabe-Ito palette while the rest are light grey, and the highlighted paths show coherent, gradual movement."} uv_plot(fit_dyn_real, plot_type = "trajectory", highlight = c("A1", "A5", "A10", "A15")) ``` Always check convergence: dynamic models have more parameters than static ones, so adequate mixing requires longer chains. Under the binary (probit) family the error variance `ve` is fixed at 1 for identification, so exclude it by its display label: ```{r trace-dyn-real, fig.width=7, fig.height=5, fig.alt="MCMC trace and density plots for the sampled parameters of the dynamic-effects fit on truly correlated data, excluding the identification-fixed error variance; well-mixed traces indicate the sampler ran long enough."} trace_plot(fit_dyn_real, exclude = "Error Variance") ``` ## What Happens Without Temporal Structure? To see how the dynamic model behaves when there is no signal, we fit the same model to independent networks -- a null baseline showing what the prior alone produces when the data are uninformative. ```{r} set.seed(6886) n <- 30 n_periods <- 5 # independent binary networks with no temporal structure Y_list <- list() for(t in 1:n_periods) { Y_t <- matrix(rbinom(n*n, 1, 0.2), n, n) diag(Y_t) <- NA rownames(Y_t) <- colnames(Y_t) <- paste0("Actor", 1:n) Y_list[[t]] <- Y_t } # fit with a custom prior on rho_uv to illustrate prior sensitivity prior_custom <- list( rho_uv_mean = 0.95, # expect very slow change in latent positions rho_uv_sd = 0.05 # tight prior ) fit_null <- lame( Y = Y_list, R = 2, dynamic_uv = TRUE, dynamic_ab = TRUE, family = "binary", prior = prior_custom, burn = 100, nscan = 1000, odens = 10, # save per-draw latent trajectories for latent_positions() / # procrustes_align(per_draw = TRUE) below posterior_opts = list(save_UV = TRUE), verbose = FALSE, plot = FALSE ) summary(fit_null) ``` We can compare the estimated persistence parameters across the two settings: ```{r} cat("Independent data rho_uv:", round(mean(fit_null$rho_uv), 3), "\n") cat("Correlated data rho_uv:", round(mean(fit_dyn_real$rho_uv), 3), "\n") ``` Even with a deliberately tight prior centred at 0.95, the two posteriors diverge: the correlated fit lands in the upper AR(1) range (expect small-sample upward bias at $T = 5$, so the posterior can sit above the true 0.85), while the null fit collapses well below the prior centre because the data carry no period-to-period correlation to anchor against. The AR(1) parameter does real work even on a modest panel, but the point estimate of `rho_uv` is no substitute for the trajectory plots (compare the coherent paths above with the erratic paths below). ### Visualizing the Null Case The trajectory plot for independent data looks different. Using the same highlight device as above, the highlighted paths jump from period to period without persistent direction, and the whole cloud clumps near the origin: ```{r uv-traj-null, fig.width=7, fig.height=5, fig.alt="Trajectory plot of latent positions for the null model with four actors highlighted in colour against the remaining 26 in light grey; the highlighted paths jump directionlessly and the whole cloud clumps near the origin, unlike the smooth dispersed paths of the genuine-dynamics fit."} uv_plot(fit_null, plot_type = "trajectory", highlight = c("Actor1", "Actor5", "Actor10", "Actor15")) ``` The additive effects behave the same way: trajectories reflecting noise rather than real shifts in actor activity. ```{r ab-traj-null, fig.width=7, fig.height=5, fig.alt="Trajectory plot of sender additive effects for the null model fit on independent networks; per-actor lines fluctuate without persistent trend, reflecting noise rather than temporal signal."} ab_plot(fit_null, effect = "sender", plot_type = "trajectory") ``` ## Comparing Model Specifications A natural workflow is to fit several specifications and compare their goodness-of-fit. We use the independent data from above, where dynamic effects should add nothing. ```{r} # four small fits to compare specifications. nscan is intentionally # small here so the vignette stays under CRAN's vignette compute # budget; for a real comparison use nscan >= 2000 per fit and run them # under `ame_parallel()` with 4 chains. fit_static <- lame(Y_list, R = 2, family = "binary", burn = 30, nscan = 150, odens = 5, verbose = FALSE, plot = FALSE) fit_uv <- lame(Y_list, R = 2, dynamic_uv = TRUE, family = "binary", burn = 30, nscan = 150, odens = 5, verbose = FALSE, plot = FALSE) fit_ab <- lame(Y_list, R = 2, dynamic_ab = TRUE, family = "binary", burn = 30, nscan = 150, odens = 5, verbose = FALSE, plot = FALSE) fit_full <- lame(Y_list, R = 2, dynamic_uv = TRUE, dynamic_ab = TRUE, family = "binary", burn = 30, nscan = 150, odens = 5, verbose = FALSE, plot = FALSE) ``` For each specification we compute posterior predictive p-values (how often simulated statistics exceed observed): values near 0.5 indicate good fit; values pinned near 0 or 1 flag a statistic the model systematically under- or over-generates. ```{r} # compute p-values for each model and each GOF statistic compute_pvals <- function(fit) { gof <- fit$GOF sapply(names(gof), function(stat) { mat <- gof[[stat]] obs <- mat[, 1] # first column = observed sims <- mat[, -1] # remaining = posterior predictive mean(colMeans(sims) >= mean(obs)) }) } gof_comparison <- rbind( Static = compute_pvals(fit_static), Dynamic_UV = compute_pvals(fit_uv), Dynamic_AB = compute_pvals(fit_ab), Full = compute_pvals(fit_full) ) round(gof_comparison, 3) ``` `sd.rowmean` and `sd.colmean` sit near 1 for *every* specification -- an artifact of fitting a heterogeneity model to homogeneous iid data, not evidence about dynamics -- so compare specifications on the `dyad.dep`, `cycle.dep`, and `trans.dep` columns instead. There, no specification is clearly better, as expected for data generated without temporal structure; with short chains these p-values are noisy, so read the table as "the dynamic terms add nothing", not as a ranking. ## When to Use Dynamic Effects **Use `dynamic_ab`** when actors' overall activity levels change over time: countries cycle through isolationist and interventionist periods, users churn in and out of platforms. **Use `dynamic_uv`** when the underlying community structure is shifting -- a stronger claim: not just that actors are more or less active, but that *who connects with whom* is changing (political realignment, generational turnover). **Use `dynamic_beta`** when a covariate's *effect* on tie formation changes over time while actors and community structure are stable: the question is *how strongly X translates into Y*, not who is connecting (`dynamic_uv`) or who is active (`dynamic_ab`). **Use multiple together** when several types of change happen at once -- the most flexible but most data-hungry specification; with short panels or sparse networks the dynamic parameters may not be well-identified. **Stick with static** when you have few time periods, when the network is genuinely stable, or when you primarily care about *average* covariate effects. ## Dealing with Rotational Indeterminacy One subtlety of dynamic latent space models: the latent space is only identified up to an orthogonal rotation **at each time point** *and* **at each posterior draw**. Even under smooth true dynamics, raw $U_t$ estimates can "jump" between periods, and posterior means can shrink toward zero across differently rotated draws. `procrustes_align()` applies the orthogonal-Procrustes rotation (Gower 1975) to align each period's positions to the previous one: ```{r} # align latent positions across time aligned <- procrustes_align(fit_null) str(aligned$U) # 3D array: actors x dimensions x time ``` By default this aligns the **posterior-mean** trajectory `fit$U` across time only. Rotation indeterminacy is a per-draw property, though, so the stricter standard (Hoff 2005; Sewell & Chen 2015) is to align each posterior draw to a reference before summarising: `procrustes_align(fit, per_draw = TRUE)` does that on the per-draw cube a dynamic fit stores when you opt in via `posterior_opts = list(save_UV = TRUE)` -- as we did for `fit_null`. Without the opt-in it falls back to mean-trajectory alignment with a note; see `?procrustes_align`. `latent_positions(fit, align = TRUE)` returns the aligned positions as a tidy data frame ready for `ggplot2` (custom trajectory plots, overlaying external events). Because this fit saved its per-draw trajectories, the `posterior_sd` column is computed from the draws rather than returned as `NA` -- though those SDs fold in some rotational wobble across draws: ```{r} lp <- latent_positions(fit_null, align = TRUE) head(lp) ``` ## References 1. **Hoff, PD (2021)**. Additive and Multiplicative Effects Network Models. *Statistical Science* 36, 34--50. 2. **Hoff, PD (2005)**. Bilinear mixed-effects models for dyadic data. *Journal of the American Statistical Association* 100(469), 286--295. 3. **Sewell, D. K., & Chen, Y. (2015)**. Latent space models for dynamic networks. *Journal of the American Statistical Association*, 110(512), 1646-1657. 4. **Carter, C. K., & Kohn, R. (1994)**. On Gibbs sampling for state space models. *Biometrika* 81(3), 541--553. The FFBS algorithm behind `dynamic_beta_kind = "ar1"` / `"rw1"`. 5. **Gower, J. C. (1975)**. Generalized Procrustes analysis. *Psychometrika* 40(1), 33--51. The orthogonal-Procrustes rotation underlying `procrustes_align()`.