AME model fitting routine for longitudinal relational data

An MCMC routine providing a fit to an additive and multiplicative effects (AME) regression model to longitudinal (time-series) relational data of various types. Supports both unipartite (square) and bipartite (rectangular) network structures. For cross-sectional (single time point) networks, use the ame function.

Usage

lame(
  Y,
  Xdyad = NULL,
  Xrow = NULL,
  Xcol = NULL,
  rvar = TRUE,
  cvar = TRUE,
  dcor = !symmetric,
  nvar = TRUE,
  R = 0,
  R_row = NULL,
  R_col = NULL,
  mode = c("unipartite", "bipartite"),
  dynamic_uv = FALSE,
  dynamic_ab = FALSE,
  dynamic_G = FALSE,
  dynamic_beta = FALSE,
  dynamic_rho = FALSE,
  dynamic_beta_kind = c("ar1", "rw1", "rw2", "matern32"),
  dynamic_uv_kind = c("ar1", "snap", "t"),
  family = "normal",
  intercept = !(family == "ordinal"),
  symmetric = FALSE,
  odmax = NULL,
  prior = list(),
  g = NA,
  seed = 6886,
  nscan = 10000,
  burn = 500,
  odens = 25,
  plot = FALSE,
  verbose = FALSE,
  gof = TRUE,
  start_vals = NULL,
  periodic_save = FALSE,
  out_file = NULL,
  save_interval = 0.25,
  model.name = NULL,
  save_log_lik = FALSE,
  posterior_opts = NULL,
  log_lik_path = NULL,
  log_lik_chunk_size = 10000L,
  keep_snap_draws = c("none", "summary", "draws", "chunked"),
  freeze_call = FALSE,
  dynamic_beta_pool = c("none", "rho", "sigma", "both"),
  dynamic_beta_per_actor = NULL,
  per_actor_covariate_idx = 1L,
  per_actor_identifiability = c("center", "exact_center", "drop_population"),
  keep_per_actor = c("auto", "draws", "summary", "none"),
  time_index = NULL,
  period_exposure = NULL,
  max_seconds = Inf,
  checkpoint_path = NULL,
  checkpoint_every = 100L,
  log_lik_method = c("observed_exact", "observed_ghk", "augmented"),
  ordinal_cutpoints = c("data_induced", "explicit"),
  method = c("mcmc", "als"),
  als_stability = c("none", "quick", "validation"),
  als_max_iter = 1000L,
  als_tol = NULL,
  bootstrap = 0L,
  bootstrap_type = c("parametric", "block"),
  bootstrap_block_length = 1L,
  bootstrap_seed = NULL,
  resume_from = NULL,
  print
)

Arguments

Y

a T length list of relational matrices, or a 3D array of dimensions [n_row, n_col, T], where T is the number of time periods. Named lists may have changing actor composition; actors are aligned by row and column names. A longitudinal netify object is also accepted and converted with netify::to_lame(lame = TRUE). When family or mode is omitted, the value inferred by netify is used. See family below for data types.

Xdyad

a T length list of n x n x pd arrays of covariates

Xrow

a T length list of n x pr matrices of nodal row covariates

Xcol

a T length list of n x pc matrices of nodal column covariates

rvar

logical: fit row random effects (asymmetric case)?

cvar

logical: fit column random effects (asymmetric case)?

dcor

logical: fit a dyadic correlation (asymmetric case)?

nvar

logical: fit nodal random effects (symmetric case)?

R

integer: dimension of the multiplicative effects (can be zero)

R_row

integer: for bipartite networks, dimension of row node multiplicative effects (defaults to R)

R_col

integer: for bipartite networks, dimension of column node multiplicative effects (defaults to R)

mode

character: either "unipartite" (default) for square networks or "bipartite" for rectangular networks

dynamic_uv

logical: fit dynamic multiplicative effects (latent factors). The transition model is selected by dynamic_uv_kind. fit$U / fit$V store only the posterior-mean trajectory cubes; draw-level trajectory uncertainty requires refitting with posterior_opts = list(save_UV_draws = TRUE), which attaches the per-iteration cubes as fit$U_draws / fit$V_draws. Default FALSE.

dynamic_ab

logical: fit dynamic additive effects (sender/receiver effects) that evolve over time using AR(1) processes. When TRUE, the row effects (a) and column effects (b) become time-varying, following $a_{i,t} = \rho_{ab} a_{i,t-1} + \epsilon_{i,t}$. This captures temporal heterogeneity in actors' baseline propensities to send and receive ties, allowing for smooth changes in activity levels and popularity over time. For example, an actor's tendency to form outgoing ties might gradually increase or decrease across observation periods. The AR(1) specification ensures temporal smoothness while allowing for actor-specific evolution patterns. Implementation uses conjugate updates where possible and C++ for computational efficiency. Default FALSE.

dynamic_G

logical (bipartite only). When TRUE the bipartite interaction matrix $G_t$ in $U_t G_t V_t'$ varies by period and is returned as fit$G_cube, an $R_\text{row} \times R_\text{col} \times T$ array. MCMC fits also return posterior summaries for this cube; ALS fits return the penalized point path using the same storage names. The marginal bilinear predictor is identified under the package's canonicalization conventions, while individual entries of $G_t$ can move under equivalent rotations and scalings of $U_t$ and $V_t$. Default FALSE.

dynamic_beta

logical, character, integer, or logical-vector flag selecting which regression coefficients evolve over time via independent AR(1) processes. Default FALSE keeps every coefficient static (the historical behaviour). Accepted forms:

FALSE / NULL: no coefficient is dynamic.
TRUE: every coefficient is dynamic.
character vector of block shortcuts ("intercept", "dyad", "row", "col") or specific coefficient names from colnames(BETA) – those become dynamic.
integer vector: 1-based column indices into colnames(BETA).
logical vector of length p: per-coefficient mask.

Each dynamic-coefficient block (one per distinct intercept / dyad / row / col label that has at least one dynamic coefficient) gets its own AR(1) parameters $\rho_\beta$ (truncated-Normal prior, default mean 0.8, bounded between 0 and 0.999) and $\sigma_\beta^2$ (inverse-Gamma prior, defaults shape 2 / scale 1). The joint posterior over the time path $\beta_t$ is drawn by Forward-Filter / Backward-Sample (FFBS) inside the MCMC loop, conditional on the current static-block beta and on (a, b, U, V). When the intercept or a nodal coefficient is dynamic, a sum-to-zero contrast basis is applied to the additive-effect sampler to keep $a_i + intercept_t$ identified. Requires at least 2 time periods. dynamic_beta composes with multiplicative latent factors (R > 0) and additive random effects for every family, across unipartite/bipartite and directed/symmetric panels and in combination with dynamic_ab / dynamic_uv. Note: with a single dynamic coefficient block and small $T$, the AR(1) parameter $\rho_\beta$ sees only $T - 1$ transitions and is strongly informed by its default truncated-Normal prior (mean 0.8, sd 0.15) rather than the data; if $\rho_\beta$ itself is of interest, check prior sensitivity (see dynamic_beta_prior_summary). With method = "als", selected intercept, dyadic, row-node, and column-node coefficients are supported on normal, binary, and Poisson panels, including named panels where actors enter or exit. Dynamic ALS uses period-specific Xrow or Xcol values for selected dynamic node coefficients; static node coefficients still use per-actor means. Dynamic intercepts are unavailable when family = "ordinal". Default FALSE.

dynamic_rho

logical. For directed unipartite normal models, allow the residual dyadic reciprocity parameter $\rho_t$ to vary by period instead of using one pooled dyadic correlation across the whole panel. MCMC fits store per-draw paths in fit$RHO and the posterior mean path in fit$rho_path. ALS dynamic fits report a residual rho_path that can be passed to MCMC with als_start_vals. Default FALSE.

dynamic_beta_kind

character: state-space prior on the dynamic coefficient block(s). "ar1" (default) gives mean-reverting AR(1) with a truncated-Normal prior on $\rho_\beta$. "rw1" gives a random walk ($\rho_\beta$ pinned at 1, no truncation), appropriate when you expect permanent drift with no mean reversion – e.g. trade-gravity coefficients in a permanently changing world economy. The alias "random_walk" is also accepted for "rw1". When the stationarity warning in summary(fit) fires under the default AR(1), refit with dynamic_beta_kind = "rw1" for cleaner inference. Decision tree: mean-reverting? AR(1). Permanent / unit-root drift? RW1. Smooth with curvature? "rw2" (second-order random walk). Smooth with a known length-scale? "matern32" (Matern 3/2). Note: "rw2" and "matern32" use an R-level joint Gaussian sampler and run substantially slower per iteration than the C++ FFBS used for "ar1" / "rw1" (roughly 2-4x on n = 100, T = 10).

dynamic_uv_kind

Character string selecting the transition model for dynamic multiplicative latent positions when dynamic_uv = TRUE. "ar1" uses Gaussian AR(1) drift, "snap" uses a mixture of AR(1) drift and discontinuous reset transitions, and "t" uses heavy-tailed Student-t innovations represented by local transition scales. On the MCMC path, "snap" and "t" are currently supported for unipartite directed and symmetric models only. With method = "als", "snap" is also supported for normal bipartite panels when it is the only dynamic block and $G$ is static. With method = "als", "t" is supported for bipartite normal, binary, and Poisson panels. ALS Student-t dynamic UV fits attach the final local transition-weight matrices as fit$lambda_u and fit$lambda_v.

family

character: one of "normal","binary","ordinal","cbin","frn","poisson" - see the details below

intercept

logical: fit model with an intercept?

symmetric

logical: is the sociomatrix symmetric?

odmax

a scalar integer or vector of length n giving the maximum number of nominations that each node may make - used for "frn" and "cbin" families

prior

a list containing hyperparameters for the prior distributions. Available options and their defaults:

Sab0: Prior scale matrix for the additive-effects covariance. A 2x2 matrix where Sab0\[1,1\] is the prior variance for row effects, Sab0\[2,2\] is the prior variance for column effects, and off-diagonals control correlation between row and column effects. For the "normal" and "poisson" families, and for unipartite "binary", this defaults to diag(2) scaled to the observed sender/receiver heterogeneity in Y rather than to a fixed diag(2); see ame for the rationale. Pass diag(2) explicitly for a fixed unit-scale prior.
eta0: Prior degrees of freedom for the additive-effects covariance $\Sigma_{ab}$ (default: round(4 \* vdfmlt) for "binary" fits, where vdfmlt is a probit-moment variance multiplier estimated from Y, and round(4 + 3 \* n/100) otherwise, where n is the number of actors). Higher values impose stronger shrinkage of the row/column effects toward the prior scale. The multiplicative-effects prior is fixed in the longitudinal sampler; Suv0 and kappa0 apply to ame only.
etaab: Prior degrees of freedom for covariance of additive effects (default: 4 + 3 \* n/100). Controls shrinkage of row/column random effects.
ab_min_observed: Minimum observed sender/receiver cells used before a static additive effect is left unshrunk (default: max(10, T)). Actors with less support are pulled toward the prior mean for that role.
rho_uv_mean: For dynamic_uv=TRUE: Prior mean for UV AR(1) parameter (default: 0.9). Values close to 1 indicate high temporal persistence.
rho_uv_sd: For dynamic_uv=TRUE: Prior SD for UV AR(1) parameter (default: 0.1). Controls uncertainty about temporal dependence.
sigma_uv_shape: For dynamic_uv=TRUE: Shape parameter for inverse-gamma prior on UV innovation variance (default: 2).
sigma_uv_scale: For dynamic_uv=TRUE: Scale parameter for inverse-gamma prior on UV innovation variance (default: 1).
uv_max_abs: For dynamic_uv=TRUE: loose bound on the raw latent coordinate scale after the product-preserving normalization step (default: 50). This guards the unidentified coordinate scale without changing the fitted $U_t V_t'$ bilinear term.
rho_ab_mean: For dynamic_ab=TRUE: Prior mean for additive effects AR(1) parameter (default: 0.8). Controls temporal smoothness of sender/receiver effects.
rho_ab_sd: For dynamic_ab=TRUE: Prior SD for additive effects AR(1) parameter (default: 0.15).
sigma_ab_shape: For dynamic_ab=TRUE: Shape parameter for inverse-gamma prior on additive effects innovation variance (default: 2).
sigma_ab_scale: For dynamic_ab=TRUE: Scale parameter for inverse-gamma prior on additive effects innovation variance (default: 1).
rho_beta_mean: For dynamic_beta: Prior mean for the per-block AR(1) parameter on time-varying regression coefficients (default: 0.8). Closer to 1 = smoother evolution.
rho_beta_sd: For dynamic_beta: Prior SD for the per-block AR(1) parameter (default: 0.15).
rho_beta_lower: For dynamic_beta: Lower truncation bound on the AR(1) parameter (default: 0). Pass a negative value to allow negative autoregression.
rho_beta_upper: For dynamic_beta: Upper truncation bound (default: 0.999). Closer to 1 admits near-unit-root behaviour.
sigma_beta_shape: For dynamic_beta: Shape parameter for the inverse-Gamma prior on the per-block innovation variance (default: 2).
sigma_beta_scale: For dynamic_beta: Scale parameter for the inverse-Gamma prior (default: 1).
sigma_beta_init: For dynamic_beta: Initial value of the per-block innovation standard deviation (default: 0.25). Affects mixing, not the stationary distribution.
beta0_mean: For dynamic_beta: Mean of the Gaussian prior on the dynamic coefficients at t = 0 (default: 0).
beta0_var: For dynamic_beta: Variance of the Gaussian prior on the dynamic coefficients at t = 0 (default: 10). Weakly informative.

Common usage: prior = list(Sab0 = diag(c(1, 1)), eta0 = 10) for stronger shrinkage, or prior = list(rho_uv_mean = 0.95) for higher temporal persistence, or prior = list(rho_beta_mean = 0.95, sigma_beta_scale = 0.1) for very smooth time-varying coefficients with tight innovations.

g

optional scalar or vector for g-prior on regression coefficients. Default is p^2 where p is the number of regression parameters. The g-prior controls the variance of regression coefficients: larger values allow for larger coefficient values. Can be a vector of length p for parameter-specific control.

seed

random seed for the MCMC sampler (default 6886). The sampler is seeded internally with this value, so results are reproducible by default and an external set.seed() call has no effect on the chain – pass a different seed here to vary the draws (e.g. when running multiple chains). The caller's .Random.seed is restored on exit, so fitting never perturbs your RNG stream.

nscan

number of iterations of the Markov chain (beyond burn-in)

burn

burn in for the Markov chain

odens

output density for the Markov chain

plot

logical: plot results while running?

verbose

logical: print progress while running? Default FALSE.

gof

logical: calculate goodness of fit statistics?

start_vals

list of parameter starting values for the MCMC chain. Use als_start_vals to start MCMC from an ALS fit.

periodic_save

logical: indicating whether to periodically save MCMC results

out_file

character vector indicating name and path in which file should be stored if periodic_save is selected. For example, on an Apple OS out_file="~/Desktop/ameFit.rda".

save_interval

quantile interval indicating when to save during the post-burn-in period.

model.name

optional string for model selection output

save_log_lik

one of FALSE (default), TRUE, or "chunked". When TRUE, stores the per-iteration pointwise log-likelihood matrix on fit$log_lik (an [n_stored, n_obs] double matrix). When "chunked", streams the log-lik values to per-column-chunk binary files under log_lik_path so the in-memory cost during MCMC is just one chunk's width; the chunks are recovered later via read_log_lik(fit). Required for loo::loo(fit) and loo::waic(fit).

posterior_opts

optional list of posterior draw-storage options, usually built with posterior_options. Recognised names (unknown names trigger a warning): save_UV, save_UV_draws, save_ab, thin_UV, thin_ab. save_UV = TRUE (the ame spelling) and save_UV_draws = TRUE are aliases on this path. With dynamic_uv = TRUE they store the per-iteration dynamic latent-position cubes at every stored (odens) iteration on fit$U_draws / fit$V_draws, each an [actor, dim, period, draw] array. On a static fit with a positive multiplicative rank they store the per-iteration latent-position matrices on fit$U_samples / fit$V_samples ([actor, dim, draw] arrays, as ame does). Bipartite fits with positive latent ranks additionally store the interaction-matrix draws on fit$G_samples ([R_row, R_col, draw]), except under dynamic_G = TRUE whose per-period cube already gets posterior summaries. save_ab is not available on the longitudinal path and is ignored with a warning; thin_UV / thin_ab are accepted for posterior_options compatibility but not applied here (odens already thins the stored chain). Requesting latent draw storage with a zero multiplicative rank warns and stores nothing. Default NULL.

log_lik_path

directory to write log-lik chunks to when save_log_lik = "chunked". Default NULL creates a session-scoped subdirectory under tempdir().

log_lik_chunk_size

column-width of each on-disk chunk when save_log_lik = "chunked". Larger chunks mean fewer files (and one readBin per chunk on read) but more memory during MCMC. Default 10000L.

keep_snap_draws

controls storage of post-burn-in MCMC snap-shift indicators when dynamic_uv = TRUE and dynamic_uv_kind = "snap". "none" keeps only the posterior mean snap probabilities. "summary" stores cell-level count, mean, and variance summaries. "draws" stores retained binary indicator arrays as fit$snap_draws and, for directed fits, fit$snap_draws_v, with dimensions draw x actor x time and first period set to NA. "chunked" currently stores the same arrays in memory and labels the storage mode for downstream tooling.

freeze_call

logical: if TRUE, store a snapshot of the evaluated Y, Xdyad, Xrow, Xcol on fit$data_snapshot so that a later update(fit, ...) refits against the same data even if the caller has mutated those objects in their workspace. Memory cost equals the size of the data; default FALSE.

dynamic_beta_pool

one of "none" (default), "rho", "sigma", or "both". When two or more coefficient blocks are dynamic, non-"none" values update shared hyperparameters across blocks for the selected persistence and / or innovation-scale component. With a single dynamic block the value is recorded but there is no other block to pool against.

dynamic_beta_per_actor

optional, one of NULL (default), "row", or "col". When non-NULL, fits actor-specific time-varying deviations for one dyadic covariate on the selected margin. The default per_actor_identifiability = "center" keeps the population coefficient in coef(fit) and enforces a per-period sum-to-zero constraint on the actor deviations. Draws are stored in fit$THETA_ACTOR when retained; summary mode stores fit$theta_actor_mean and fit$theta_actor_sd.

per_actor_covariate_idx

positive integer; index into the dyadic covariate cube to slope on for the per-actor extension. Default 1L.

per_actor_identifiability

one of "center" (default) or "drop_population". "center" preserves the population coefficient and constrains per-actor deviations to sum to zero per period via pairwise contrast FFBS. "drop_population" is reserved and currently falls back to "center".

keep_per_actor

one of "auto" (default; full draws when memory cost < 250 MB else streaming summary), "draws" (full [n_iter, n_actors, T] cube), "summary" (streaming posterior mean + sd only), or "none" (only hyperparameter chains).

time_index

optional numeric vector of length $T$ giving observation times. Default NULL treats periods as equally spaced. Strictly-increasing values are required when supplied. Unequal gaps scale the AR(1) / RW1 transition variance and are also used by the RW2 and Matérn 3/2 dynamic-beta precision builders.

period_exposure

optional non-negative numeric vector of length $T$ giving period-level exposure offsets. Wired for Poisson: when supplied with any value != 1, the Poisson observation likelihood becomes $Y_{ij,t} \sim \text{Poisson}(e_t \cdot \exp(Z_{ij,t}))$, while the latent $Z$ retains its existing semantic of "unexposed log-rate". For other families, non-trivial period_exposure is an error (rescale Y or use a covariate offset). NULL (default) or all-ones uses the unscaled likelihood path.

max_seconds

optional positive scalar; if the MCMC wall-clock time exceeds this many seconds, the chain terminates cleanly and fit$terminated_early is set to TRUE.

checkpoint_path

optional file path. When set, the chain periodically writes a snapshot of BETA, VC, the RNG state, and the original call to this file. Use lame_resume or pass the same path back as resume_from = path to continue.

checkpoint_every

positive integer; iterations between checkpoint writes. Default 100L.

log_lik_method

one of "observed_exact" (default), "observed_ghk", or "augmented". Selects the argument for selecting which pointwise log-lik is stored on fit$log_lik. "observed_exact" uses the closed-form marginal log-likelihood for normal, binary, cbin, poisson, and ordinal fits. "observed_ghk" uses a randomized-Halton GHK Monte Carlo marginal for rank-family fits ("frn") and the exact observed likelihood when no GHK step is needed. "augmented" uses the augmented-data Gaussian-on-$Z$ contribution. Note the GHK cost grows exponentially with the per-row dimension $D$ (the number of observed dyads in a sender row): the per-row Monte Carlo budget escalates to min(2048, max(64, 8 * 2^D)) draws across a bank of 8 randomized Halton sequences, once per stored draw, so "observed_ghk" on a unipartite "frn" fit with more than roughly 12 actors (where $D = n - 1$) is expensive; rows with $D > 15$ fall back to a fast pairwise-normal approximation. Power users can fix the per-row budget exactly via prior$ghk_n_mc (a single positive integer), which both raises the 2048 ceiling and lowers the 8 * 2^D accuracy floor; lower budgets increase the downward Jensen bias of the log-lik estimate.

ordinal_cutpoints

character: cutpoint convention for family = "ordinal". "data_induced" (default) uses the data-induced cutpoints; "explicit" samples explicit cutpoints via a Cowles (1996) Metropolis-Hastings update. Ignored for other families.

method

character: "mcmc" (default, the Bayesian MCMC fit) or "als" (the fast, MCMC-free point-estimation path). With no dynamic arguments, method = "als" uses the pooled-static lame_als estimator. For normal, binary, and Poisson panels, method = "als" also supports a dynamic penalized point-estimation path for dynamic_ab, selected dynamic_beta, smooth AR(1) dynamic_uv in directed, symmetric, and bipartite panels, bipartite dynamic_G, and Student-t dynamic_uv in directed, symmetric, and bipartite panels. The snap-only dynamic_uv = TRUE, dynamic_uv_kind = "snap" case routes to lame_snap_als for supported normal unipartite and bipartite panels. Node-covariate coefficients use the same orthogonal additive-effect decomposition as lame_als; dynamic node coefficients use period-specific node values when selected by dynamic_beta, while static node coefficients use per-actor means. Named changing-composition panels are aligned to the union actor set and actor-entry gaps break the dynamic smoothing penalties. Posterior draws, rank/censored families, and bipartite snap-shift fits with node covariates, dynamic coefficients, or dynamic G_t are not handled by the snap ALS shortcut.

als_stability

character (only used when method = "als" and a dynamic ALS path is selected): "none" (default), "quick", or "validation". The non-"none" presets rerun the dynamic ALS fit from jittered starts and attach fit$stability. Smooth dynamic ALS fits report fitted-surface and coefficient-path differences across starts; snap ALS fits report snap-score, snap-class, top-period, and fitted-surface differences.

als_max_iter

positive integer (only used when method = "als"): maximum block-coordinate iterations for the ALS fit. Default 1000.

als_tol

optional positive scalar (only used when method = "als"): convergence tolerance for the ALS objective and fitted values. NULL keeps each ALS estimator's built-in default.

bootstrap

integer (only used when method = "als"): number of bootstrap replicates. 0 (default) skips the bootstrap; N > 0 runs N replicates and attaches the result so that confint returns bootstrap intervals.

bootstrap_type

character (only used when method = "als"): "parametric" (default) or "block".

bootstrap_block_length

integer: block length for the block bootstrap.

bootstrap_seed

optional integer seed for the bootstrap.

resume_from

optional path to a checkpoint file produced by a previous lame(..., checkpoint_path = path) call. When non-NULL, lame() short-circuits all other input parsing and delegates to lame_resume with the user-supplied overrides forwarded. Pass nscan = K to request K additional stored draws on the continuation. Note: checkpoint_path and max_seconds cannot currently be overridden on a resume (they are stripped before re-evaluating the saved call); open both at the original lame() call if you need iterative checkpointing or a time budget on a chain of resumes. Pass verbose = FALSE on resume (the sampling progress bar is gated on burn != 0 and the continuation forces burn = 0). Default NULL.

print

Deprecated. Use verbose instead.

Value

BETA

posterior samples of regression coefficients. A 2-dimensional matrix [n_stored, p] when all coefficients are static. When dynamic_beta flags any coefficient as time-varying, BETA is a 3-dimensional array [n_stored, p, T] whose third dimension is the time period. Static coefficients are still present with their values replicated across the third dimension. coef(fit) collapses this to a [p, T] posterior-mean matrix.

Migration from amen. Under amen::ame() the BETA slot was always 2-D. Scripts that compute apply(fit$BETA, 2, mean) on an amen fit will silently aggregate across periods when run against a lame() fit with dynamic_beta active (the second margin is the coefficient index in both shapes; the new third margin is time). Use length(dim(fit$BETA)) to detect the shape, or call coef(fit) which returns a [p, T] matrix in either case.

VC

posterior samples of the variance parameters

APM

posterior mean of additive row effects a

BPM

posterior mean of additive column effects b

U

posterior mean of multiplicative row effects u. For dynamic_uv=TRUE, this is a 3D array (n x R x T)

V

posterior mean of multiplicative column effects v (asymmetric case). For dynamic_uv=TRUE, this is a 3D array (n x R x T)

UVPM

posterior mean of UV

ULUPM

posterior mean of ULU (symmetric case)

L

posterior mean of L (symmetric case)

EZ

estimate of expectation of Z matrix. For mode = "bipartite", EZ is a list of per-period matrices whose row/column dimnames are currently NULL. Use names(fit$APM) / names(fit$BPM) (or dimnames(fit$YPM[[t]])) to recover the row/column actor ordering, which is the lexicographic sort of the input actor names (so "r10" sorts before "r2" unless you zero-pad).

YPM

posterior mean of Y (for imputing missing values)

GOF

observed (first row) and posterior predictive (remaining rows) values of four goodness-of-fit statistics. See gof for post-hoc computation and gof_plot for visualization.

start_vals

Final parameter values from MCMC, can be used as the input for a future model run.

model.name

Name of the model (if provided)

Details

This command provides posterior inference for parameters in AME models of longitudinal relational data, assuming one of eight possible data types/models. The model supports both unipartite networks (square adjacency matrices) and bipartite networks (rectangular adjacency matrices with distinct row and column node sets) across multiple time points.

Dynamic Effects Implementation:

The dynamic_uv and dynamic_ab parameters enable time-varying latent representations through autoregressive processes. These extensions are particularly useful for understanding how network structure evolves over time.

Dynamic Multiplicative Effects (dynamic_uv=TRUE): The latent factors U and V evolve according to the transition selected by dynamic_uv_kind. The default "ar1" model uses Gaussian AR(1) drift: $$U_{i,k,t} = \rho_{uv} U_{i,k,t-1} + \epsilon_{i,k,t}$$ where $\epsilon_{i,k,t} \sim N(0, \sigma_{uv}^2)$, i indexes actors, k indexes latent dimensions, and t indexes time. The parameter $\rho_{uv}$ controls temporal persistence (values near 1 indicate slow evolution). This captures time-varying homophily, latent community structure, and transitivity dynamics. The "snap" transition uses a mixture of AR(1) drift and discontinuous reset transitions, while "t" uses heavy-tailed Student-t innovations represented by local transition scales. On the MCMC path, snap and t transitions are supported for unipartite directed and symmetric models. With method = "als", Student-t dynamic UV is also supported for bipartite normal, binary, and Poisson panels, and snap ALS is supported for normal unipartite and bipartite panels, including named panels where actors enter or exit.

Key references:

Sewell & Chen (2015): Introduced dynamic latent space models with actor-specific evolution rates
Durante & Dunson (2014): Nonparametric Bayesian approach allowing flexible evolution of network structure
Hoff (2011): Hierarchical multilinear models providing theoretical foundation for temporal dependencies

Dynamic Additive Effects (dynamic_ab=TRUE): The sender (a) and receiver (b) effects evolve as: $$a_{i,t} = \rho_{ab} a_{i,t-1} + \epsilon_{i,t}$$ $$b_{i,t} = \rho_{ab} b_{i,t-1} + \eta_{i,t}$$ where $\epsilon_{i,t}, \eta_{i,t} \sim N(0, \sigma_{ab}^2)$. This models time-varying individual activity levels (outdegree) and popularity (indegree).

Dynamic additive effects are useful when sender activity and receiver popularity change over time rather than staying fixed across the whole panel. A dynamic_ab fit returns the posterior-mean paths on fit$a_dynamic / fit$b_dynamic (actor x period) and their per-period posterior standard deviations on fit$a_dynamic_sd / fit$b_dynamic_sd; ab_plot(fit, plot_type = "ribbon") uses the latter to draw a credible band around each actor's path.

Prior Specification for Dynamic Parameters:

$\rho_{uv}, \rho_{ab} \sim TruncNormal(mean, sd, 0, 1)$: Ensures stationarity of AR(1) process
$\sigma_{uv}^2, \sigma_{ab}^2 \sim InverseGamma(shape, scale)$: Controls innovation variance
Default priors ($\rho_{uv}$ mean=0.9, $\rho_{ab}$ mean=0.8) favor smooth evolution
Adjust rho_*_mean closer to 1 for slower evolution, closer to 0 for more rapid changes

Computational Considerations:

Dynamic effects increase computation by roughly 30 to 50 percent per iteration
Memory usage scales as O(nRT) for dynamic_uv, O(n*T) for dynamic_ab
C++ implementation is substantially faster than a pure R loop
Convergence diagnostics: Monitor rho and sigma parameters carefully
Effective sample sizes typically lower due to temporal correlation
Recommend burn >= 1000 and nscan >= 20000 for dynamic models

Model Selection Guidelines: Use both dynamic_uv and dynamic_ab when:

Networks show clear temporal trends in density or clustering
Individual node behavior changes systematically over time
Community structure evolves (merging, splitting, drift)

Use only dynamic_uv when:

Latent structure/communities change but individual effects are stable
Focus is on evolving homophily or clustering patterns
Network shows structural reconfiguration over time

Use only dynamic_ab when:

Individual heterogeneity varies but overall structure is stable
Actors' activity/popularity changes over observation period
Focus is on individual-level temporal dynamics

Bipartite Network Models:

When mode="bipartite", the model handles rectangular adjacency matrices Y with dimensions n_A x n_B, where n_A and n_B represent the number of row and column nodes respectively.

Static Bipartite Case: The model uses separate latent factor matrices:

U: n_A x R_row matrix of row node latent positions
V: n_B x R_col matrix of column node latent positions
G: R_row x R_col interaction matrix mapping between latent spaces
Multiplicative term: U G V' captures bipartite community structure

Dynamic Bipartite Case: When dynamic_uv=TRUE for bipartite networks: $$U_{i,k,t} = \rho_{uv} U_{i,k,t-1} + \epsilon_{i,k,t}$$ $$V_{j,k,t} = \rho_{uv} V_{j,k,t-1} + \eta_{j,k,t}$$ where i indexes row nodes, j indexes column nodes, k indexes latent dimensions.

When dynamic_G = TRUE the bipartite interaction matrix $G_t$ varies by period (a separate $R_\text{row} \times R_\text{col}$ matrix at every time slice) and is returned as fit$G_cube. The MCMC estimator samples $G_t$ with a Carter-Kohn/FFBS update; method = "als" estimates a penalized point path for normal, binary, and Poisson bipartite panels, including named changing-composition panels. The marginal $U_t G_t V_t'$ linear predictor is the identified object; individual $G_t$ entries can shift under equivalent rotations and scalings of $U_t$ and $V_t$. dynamic_G = TRUE is bipartite only.

Key Differences from Unipartite Models:

No dyadic correlation (rho): Bipartite edges are inherently directed
Separate dimensions: R_row and R_col can differ for row/column spaces
Rectangular structure: Network density patterns differ from square matrices
Community interpretation: Captures affiliation patterns between node types

Standard AME Model Types:

The following describes the six standard data types/models available:

"normal": A normal AME model.

"binary": A binary probit AME model.

"ordinal": An ordinal probit AME model. An intercept is not identifiable in this model.

"cbin": An AME model for censored binary data. The value of 'odmax' specifies the maximum number of links each row may have.

"frn": An AME model for fixed rank nomination networks. A higher value of the rank indicates a stronger relationship. The value of 'odmax' specifies the maximum number of links each row may have.

"poisson": An overdispersed Poisson AME model for count data: $Y \sim \mathrm{Poisson}(\exp(z))$ with $z \sim N(\eta, \sigma^2)$, a lognormal-mixed Poisson. The conditional mean given the latent $z$ is $\exp(z)$; the marginal mean is $\exp(\eta + \sigma^2/2)$, not $\exp(\eta)$.

Author

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples


data(YX_bin_list)
fit<-lame(YX_bin_list$Y,YX_bin_list$X,burn=5,nscan=5,odens=1,family="binary")
#> Warning: `family` = "binary" but `Y` contains values other than 0/1.
#> ℹ `Y` will be thresholded to `1 * (Y > 0)`; if you meant counts, use "poisson",
#>   or "ordinal"/"normal" as appropriate.
# you should run the Markov chain much longer than this

# \donttest{
## Time-varying regression coefficients (dynamic_beta).
## Make every dyadic coefficient evolve as an AR(1):
fit_dyn <- lame(YX_bin_list$Y, YX_bin_list$X,
                family = "binary", R = 0,
                nscan = 60, burn = 15, odens = 5,
                dynamic_beta = "dyad")
#> Warning: `family` = "binary" but `Y` contains values other than 0/1.
#> ℹ `Y` will be thresholded to `1 * (Y > 0)`; if you meant counts, use "poisson",
#>   or "ordinal"/"normal" as appropriate.
dim(fit_dyn$BETA)        # [n_stored, p, T] -- 3-D when dynamic
#> [1] 12  4  4
coef(fit_dyn)            # [p, T] posterior-mean coefficient paths
#>                   t1         t2         t3         t4
#> intercept 0.07826748 0.07826748 0.07826748 0.07826748
#> X1_dyad   0.49111731 0.44428714 0.44571817 0.45857251
#> X2_dyad   0.57504492 0.53712751 0.54922098 0.54269183
#> X3_dyad   0.70705884 0.66310812 0.71517498 0.67996897
confint(fit_dyn)         # per-period 95% credible intervals
#>                     2.5%     97.5%
#> intercept[t1] 0.04120203 0.1162258
#> X1_dyad[t1]   0.39659728 0.5574730
#> X2_dyad[t1]   0.46336940 0.6543630
#> X3_dyad[t1]   0.57432997 0.8100819
#> intercept[t2] 0.04120203 0.1162258
#> X1_dyad[t2]   0.33614094 0.5199627
#> X2_dyad[t2]   0.42941619 0.6286343
#> X3_dyad[t2]   0.54589670 0.7464563
#> intercept[t3] 0.04120203 0.1162258
#> X1_dyad[t3]   0.36008943 0.5079309
#> X2_dyad[t3]   0.45680753 0.5965079
#> X3_dyad[t3]   0.60084123 0.7833216
#> intercept[t4] 0.04120203 0.1162258
#> X1_dyad[t4]   0.35522172 0.5317025
#> X2_dyad[t4]   0.43112989 0.6066965
#> X3_dyad[t4]   0.55522092 0.7628115
summary(fit_dyn)         # prints a "Dynamic coefficients per period" block
#> 
#> === Longitudinal AME Model Summary ===
#> 
#> Call:
#> [1] "Y ~ dyad(X1, X2, X3) + a[i] + b[j] + rho*e[ji], family = 'binary'"
#> 
#> Time periods: 4 
#> Family: binary 
#> Mode: unipartite 
#> Dynamic regression coefficients: enabled, kind = 'ar1' (rho_beta = dyad : 0.9)
#> 
#> Regression coefficients:
#> ------------------------
#>           Estimate StdError z_value p_value CI_lower CI_upper    
#> intercept    0.078    0.024   3.314   0.001    0.041    0.116 ***
#> X1_dyad       0.46    0.052   8.766       0    0.362    0.523 ***
#> X2_dyad      0.551    0.056   9.817       0    0.445    0.621 ***
#> X3_dyad      0.691    0.069  10.036       0    0.569    0.771 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Note: stars are a visual hint from posterior mean / SD only; for inference use the credible intervals.
#> 
#> Dynamic coefficients per period:
#> -------------------------------
#>            Mean   Min   Max Drift Drift_pct Dynamic
#> intercept 0.078 0.078 0.078 0.000     0.000       N
#> X1_dyad   0.460 0.444 0.491 0.047    10.182       Y
#> X2_dyad   0.551 0.537 0.575 0.038     6.881       Y
#> X3_dyad   0.691 0.663 0.715 0.052     7.531       Y
#> 
#> Variance components:
#> -------------------
#>     Estimate StdError
#> va     0.184    0.056
#> cab    0.091    0.038
#> vb     0.164    0.074
#> rho    0.161    0.051
#> ve     1.000    0.000
#>   (va = sender, cab = sender-receiver covariance, vb = receiver,
#>    rho = dyadic correlation, ve = residual variance)
# }