Package {lame}

Longitudinal Additive and Multiplicative Effects Models for Networks

Description

An extension of the 'amen' package specifically designed for longitudinal network analysis. This package provides tools for fitting Additive and Multiplicative Effects (AME) models to longitudinal/replicated relational data with several extensions: (1) handling of changing actor compositions across time periods, allowing for networks with different sets of actors at each time point; (2) significant performance improvements through C++ implementations via Rcpp and RcppArmadillo; (3) specialized functions for temporal network dynamics. The package supports six data types: normal (nrm), binary (bin), ordinal (ord), Poisson count (poisson), censored binary (cbin), and fixed-rank nomination (frn). Based on the AME framework originally developed by Hoff (2009) and Hoff, Fosdick, Volfovsky and Stovel (2013).

Details

Estimators. The package offers two estimation routes:

• ame / lame – the Bayesian MCMC estimators, for calibrated posterior inference (cross-sectional and longitudinal respectively).

• ame_als / lame_als – a fast, MCMC-free point estimator by iterative block coordinate descent, with bootstrap uncertainty via ame_als_bootstrap. Use it for rapid model exploration, rank selection and starting values; use the MCMC estimators for final inference.

Author(s)

Shahryar Minhas, Tosin Salau, Cassy Dorff

Maintainer: Shahryar Minhas minhassh@msu.edu

Examples

data(YX_nrm)
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, burn = 10, nscan = 100, odens = 1,
           family = "normal", verbose = FALSE)
summary(fit)

Univariate Carter-Kohn FFBS for one actor's length-T slope path

Description

Given period-wise sufficient statistics (Ht = sum_j X^2_ij / s2, ht = sum_j X_ij R_ij / s2) and AR(1) hyperparameters, returns one joint draw from the posterior N(m, P) where the prior is AR(1) and the observation is the Gaussian likelihood implied by (H, h).

Usage

.actor_ffbs_path(H, h, rho_actor, sigma_actor2)

Arguments

Value

list with theta (length-T draw) and V (length-T marginal posterior variance per period – needed for the exact centering projection).

Propagate covariate missingness into the response

Description

For an AME fit, a dyad whose covariate value is missing cannot contribute a covariate-coefficient observation. Rather than silently imputing the missing covariate as 0 (which biases the coefficient when the covariate is not mean-centred), the dyad is treated as an unobserved tie and handled by data augmentation. This helper sets the affected Y cells to NA; lame already does this internally, and this brings ame into line.

Usage

.ame_propagate_cov_na(Y, Xrow = NULL, Xcol = NULL, Xdyad = NULL)

Arguments

Value

Y with covariate-missing cells set to NA.

Build the full canonical design for a newdata prediction

Description

Build the full canonical design for a newdata prediction

Usage

.build_full_design(newdata, model_names, fitted_design, n, m)

Arguments

Value

an n x m x length(model_names) array whose slices are in model_names order: dyadic slices from newdata, all other slices held at their fitted values (intercept = 1).

Exact variance-weighted sum-to-zero projection for the per-actor block

Description

Given an unconstrained sample theta_star (n_actors x T) and per-actor per-period posterior variances V (same shape), projects to the sum-to-zero manifold per period using the variance-weighted Lagrangian projection (exact under block-diagonal covariance across actors).

Usage

.exact_center_per_actor(theta_star, V)

The fit's internal (canonical) actor order

Description

lame() sorts actors alphabetically when it ingests a list of Y matrices (via list_to_array), so the stored additive / multiplicative effects and design slices are in sorted actor order, which need not match the order a user later supplies in newdata. This returns the row and column actor names the fit was estimated in, tried from the most authoritative source down, so newdata can be realigned by name before prediction.

Usage

.fit_actor_order(object)

Arguments

Value

a list with rows and cols character vectors (or NULL when the order cannot be recovered).

Realign a newdata covariate array to the fit's internal actor order

Description

Reorders the rows and columns of each n x m x p slice by their dimnames so they line up with the fit's stored (sorted) additive / multiplicative effects. Only acts when newdata carries dimnames whose name set exactly matches the fit's; otherwise it falls back to the supplied positional order unchanged (the historical behaviour, correct when the user already passes actors in the fit's order or supplies no names).

Usage

.reorder_newdata_actors(X, row_order, col_order)

Arguments

Value

X with rows / columns permuted into canonical order, or unchanged.

Rotation-drift diagnostic for the canonical (U_t, G_t, V_t) trio

Description

Reports a scalar ratio comparing the variance of the raw G_t entries (per element across t) to the variance of the canonical G_t entries (after per-period SVD). When the raw / canonical variance ratio is large (>= 5), the apparent G_t time-variation is dominated by rotation drift in U_t, V_t rather than real temporal change. The canonicalisation removes that rotation and the canonical entries should be the user-facing summary.

Usage

.rotation_drift_diagnostic(G_cube_raw, U_cube, V_cube)

Arguments

Value

list with ratio (numeric), var_raw (RA x RB), var_canonical (RA x RB), flag (logical: ratio >= 5)

Per-actor sweep: unconstrained per-actor FFBS + exact projection

Description

One full sweep over actors: each actor's path is sampled from its conditional AR(1) posterior via .actor_ffbs_path(), and the joint sample is then projected to satisfy sum_i theta_i(t) = 0 per period using the variance-weighted Lagrangian.

Usage

.sweep_per_actor_exact(H_mat, h_mat, rho_actor, sigma_actor2)

Arguments

Value

n_actors x T matrix of centered draws.

Name-safe linear predictor from a newdata prediction

Description

Convenience wrapper: builds the full canonical design via .build_full_design and contracts it against a named coefficient vector. Because the design is rebuilt in the model's slice order, the positional contraction is correct even with nodal covariates.

Usage

.xbeta_newdata(beta, newdata, fitted_design, n, m)

Arguments

Value

an n x m matrix of the linear predictor.

International relations in the 90s

Description

A relational dataset recording a variety of nodal and dyadic variables on countries in the 1990s, including information on conflicts, trade and other variables. Except for the conflict variable, the variables are averages across the decade.

Format

A list consisting of a socioarray dyadvars of dyadic variables and matrix nodevars of nodal variables. The dyadic variables include

• total number of conflicts;

• exports (in billions of dollars);

• distance (in thousands of kilometers);

• number of shared IGOs (averages across the years);

• polity interaction.

The nodal variables include

• population (in millions);

• gdp (in billions of dollars);

• polity

Source

Michael Ward.

Linear combinations of submatrices of an array

Description

Computes a matrix of expected values based on an array X of predictors and a vector beta of regression coefficients.

Usage

Xbeta(X, beta)

Arguments

Value

An n by n matrix

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Compute Xbeta product for bipartite networks

Description

Computes sum_k beta_k * X_k for a single time slice

Usage

Xbeta_bip_cpp(X, beta)

Arguments

Value

nA x nB matrix

Column covariates

Description

Column covariates loaded by data(vignette_data). See vignette_data for full details.

Format

A list of matrices of column covariates

See Also

vignette_data

Dyadic covariates

Description

Dyadic covariates loaded by data(vignette_data). See vignette_data for full details.

Format

A list of arrays of dyadic covariates

See Also

vignette_data

Row covariates

Description

Row covariates loaded by data(vignette_data). See vignette_data for full details.

Format

A list of matrices of row covariates

See Also

vignette_data

Relational matrix

Description

Relational matrix loaded by data(vignette_data). See vignette_data for full details.

Format

A list of square binary adjacency matrices

See Also

vignette_data

binary relational data and covariates

Description

a synthetic dataset that includes binary relational data as well as information on eight covariates

Usage

data(YX_bin)

Format

The format is: List of 2 $ Y: num [1:100, 1:100] NA 0 0 0 0 0 0 0 0 1 ... $ X: num [1:100, 1:100, 1:8] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 3 .. ..$ : NULL .. ..$ : NULL .. ..$ : chr [1:8] "intercept" "rgpa" "rsmoke" "cgpa" ...

Examples

data(YX_bin)
gof_stats(YX_bin$Y) 

Synthetic longitudinal binary relational data, list-form (latent-scale)

Description

The list-form sibling of YX_bin_long: same synthetic 4-period 50-actor panel reshaped from arrays to lists. The shipped Y entries hold the latent probit-scale predictor (range roughly -20 to 17), not 0/1 ties. To recover the binary tie indicator the dataset name implies, threshold at zero, e.g. Y_bin <- lapply(YX_bin_list$Y, function(Yt) 1 * (Yt > 0)). If passed unthresholded to lame(..., family = "binary") the fit will warn and silently apply the same Y > 0 threshold.

Format

A list with two elements:

Y

List of 4 numeric matrices, each [50, 50], with NA on the diagonal. Latent probit-scale predictor; threshold at 0 for the binary indicator.

X

List of 4 numeric arrays, each [50, 50, 3] of dyadic covariates.

synthetic longitudinal binary relational data (latent-scale)

Description

a synthetic 50-actor, 4-period dataset used by the vignettes. The shipped Y array holds the latent probit-scale predictor (range roughly -24 to 20), not 0/1 ties. To recover the binary tie indicator the dataset name implies, threshold at zero: Y_bin <- (YX_bin_long$Y > 0) * 1. If passed unthresholded to lame(..., family = "binary") the fit will warn and silently apply the same Y > 0 threshold.

Usage

data(YX_bin_long)

Format

A list with two elements:

Y

Numeric array [50, 50, 4]. Latent probit-scale predictor; threshold at 0 to get the binary tie indicator.

X

Numeric array [50, 50, 3, 4] of dyadic covariates.

Examples

data(YX_bin_long)
# threshold latent z to 0/1 before computing binary GOF stats
Yt <- 1 * (YX_bin_long$Y[, , 1] > 0)
diag(Yt) <- NA
gof_stats(Yt)

Censored binary nomination data and covariates

Description

a synthetic dataset that includes relational data where the number of nominations per row is censored at 10, along with information on eight covariates

Usage

data(YX_cbin)

Format

The format is: List of 2 $ Y: num [1:100, 1:100] NA 0 0 0 1 0 0 0 0 3 ... $ X: num [1:100, 1:100, 1:8] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 3 .. ..$ : NULL .. ..$ : NULL .. ..$ : chr [1:8] "intercept" "rgpa" "rsmoke" "cgpa" ...

Examples

data(YX_cbin)
gof_stats(YX_cbin$Y) 

Fixed rank nomination data and covariates

Description

a synthetic dataset that includes fixed rank nomination data as well as information on eight covariates

Usage

data(YX_frn)

Format

The format is: List of 2 $ Y: num [1:100, 1:100] NA 0 0 0 1 0 0 0 0 3 ... $ X: num [1:100, 1:100, 1:8] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 3 .. ..$ : NULL .. ..$ : NULL .. ..$ : chr [1:8] "intercept" "rgpa" "rsmoke" "cgpa" ...

Examples

data(YX_frn)
gof_stats(YX_frn$Y) 

normal relational data and covariates

Description

a synthetic dataset that includes continuous (normal) relational data as well as information on eight covariates

Usage

data(YX_nrm)

Format

The format is: List of 2 $ Y: num [1:100, 1:100] NA -4.05 -0.181 -3.053 -1.579 ... $ X: num [1:100, 1:100, 1:8] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 3 .. ..$ : NULL .. ..$ : NULL .. ..$ : chr [1:8] "intercept" "rgpa" "rsmoke" "cgpa" ...

Examples

data(YX_nrm)
gof_stats(YX_nrm$Y)

ordinal relational data and covariates

Description

a synthetic dataset that includes ordinal relational data as well as information on seven covariates

Usage

data(YX_ord)

Format

The format is: List of 2 $ Y: num [1:100, 1:100] NA 0 3 0 3 1 0 1 1 0 ... $ X: num [1:100, 1:100, 1:7] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 3 .. ..$ : NULL .. ..$ : NULL .. ..$ : chr [1:7] "rgpa" "rsmoke" "cgpa" "csmoke" ...

Examples

data(YX_ord)
gof_stats(YX_ord$Y)

Visualize sender and receiver random effects

Description

Creates a visualization of the additive sender (row) and receiver (column) random effects from an AME or LAME model. Automatically detects whether effects are static or dynamic and provides appropriate visualization options.

Usage

ab_plot(
  fit,
  effect = c("sender", "receiver"),
  sorted = TRUE,
  labels = NULL,
  title = NULL,
  time_point = NULL,
  plot_type = c("snapshot", "trajectory", "faceted", "ribbon"),
  show_actors = NULL
)

Arguments

Details

The additive effects in AME models represent:

Sender effects (a)

Actor-specific tendencies to form outgoing ties. Positive values indicate actors who send more ties than expected; negative values indicate actors who send fewer ties.

Receiver effects (b)

Actor-specific tendencies to receive incoming ties. Positive values indicate actors who receive more ties than expected; negative values indicate actors who receive fewer ties.

For static effects, the plot displays these effects as a dot plot with vertical lines extending from zero to each effect estimate.

For dynamic effects (when fit contains a_dynamic/b_dynamic), additional options are available to visualize how effects evolve over time.

Value

A ggplot2 object that can be further customized

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples

# Fit an AME model
data(YX_nrm)
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X,
           nscan = 100, burn = 10, odens = 1, verbose = FALSE)

# Visualize sender effects
ab_plot(fit, effect = "sender")

# Visualize receiver effects without sorting
ab_plot(fit, effect = "receiver", sorted = FALSE)

Additive-effects plot for an ame_als fit

Description

Sender (a) and receiver (b) additive-effect point estimates as a sorted lollipop chart. Mirrors ab_plot but without posterior intervals: this is a point estimator. If the fit was produced via ame_als(..., bootstrap = N), bootstrap 95\ error bars.

Usage

ab_plot.ame_als(
  fit,
  effect = c("sender", "receiver", "both"),
  top_n = Inf,
  ...
)

Arguments

Value

A ggplot object.

AddHealth community 3 data

Description

A valued sociomatrix (Y) and matrix of nodal attributes (X) for students in community 3 of the AddHealth study.

• Y: A sociomatrix in which the value of the edge corresponds to an ad-hoc measure of intensity of the relation. Note that students were only allowed to nominate up to 5 male friends and 5 female friends.

• X: Matrix of students attributes, including sex, race (1=white, 2=black, 3=hispanic, 4=asian, 5=mixed/other) and grade.

Originally available at http://moreno.ss.uci.edu/data.html#adhealth (site no longer accessible).

Usage

data(addhealthc3)

Format

list

AddHealth community 9 data

Description

A valued sociomatrix (Y) and matrix of nodal attributes (X) for students in community 9 of the AddHealth study.

• Y: A sociomatrix in which the value of the edge corresponds to an ad-hoc measure of intensity of the relation. Note that students were only allowed to nominate up to 5 male friends and 5 female friends.

• X: Matrix of students attributes, including sex, race (1=white, 2=black, 3=hispanic, 4=asian, 5=mixed/other) and grade.

Originally available at http://moreno.ss.uci.edu/data.html#adhealth (site no longer accessible).

Usage

data(addhealthc9)

Format

list

Penalised ALS time-varying coefficient estimate

Description

Computes a fast point estimate of a time-varying coefficient vector \beta_t for a longitudinal network model by solving the first-difference-penalised least-squares problem. No MCMC, no multiplicative effects, no random effects – purely a regression point estimate with a smoothing penalty on the coefficient path. Useful for rapid exploration and for generating starting values for a full lame(dynamic_beta = ...) fit.

Usage

als_dynamic_beta(Y, Xdyad, lambda = 0, intercept = TRUE)

Arguments

Details

\min_{\beta_{1:T}} \sum_t \|y_t - X_t \beta_t\|^2 + \lambda \sum_{t=2}^{T} \|\beta_t - \beta_{t-1}\|^2

At \lambda = 0 this is T independent per-period OLS fits. As \lambda \to \infty the solution converges to a single pooled \beta (constant over time).

Value

A list with beta (a p \times T matrix of estimates), lambda, residual_ss (sum of squared residuals across periods), intercept, and call. Class "als_dynamic_beta".

Relation to other entry points

als_dynamic_beta is a regression-only smoother; the returned object carries time-varying \beta_t alone, no additive (a, b) or multiplicative (U, V) AME components. It is therefore not a special case of lame_als (which always estimates the full AME decomposition) and not a special case of lame(dynamic_beta = TRUE) (which is Bayesian with an AR(1) / RW1 / RW2 / Matern 3/2 state-space prior). Use it when you want a fast, point-only, deterministic smoother for the coefficient path; use lame(dynamic_beta = ...) when you need full posterior uncertainty and additive / multiplicative effects.

See Also

lame for the full Bayesian dynamic-\beta fit (AR(1) / RW1 / RW2 / Matern 3/2); lame_als for the longitudinal AME point estimator (static \beta, full a, b, U, V).

Examples

set.seed(1)
n <- 10; Tn <- 4
X <- replicate(Tn, array(rnorm(n*n*2), c(n, n, 2)), simplify = FALSE)
beta_true <- rbind(seq(-1, 1, length.out = Tn), seq(0.5, -0.5, length.out = Tn))
Y <- vector("list", Tn)
for (t in seq_len(Tn)) {
  Yt <- X[[t]][, , 1] * beta_true[1, t] + X[[t]][, , 2] * beta_true[2, t] +
        matrix(rnorm(n*n, 0, 0.2), n, n)
  diag(Yt) <- NA
  Y[[t]] <- Yt
}
fit_als <- als_dynamic_beta(Y, X, lambda = 0)     # per-period LS
fit_smooth <- als_dynamic_beta(Y, X, lambda = 10) # smoother path

Convert an ALS fit to MCMC starting values

Description

Builds a start_vals list for ame or lame from an ALS point estimate. This is useful when the ALS fit has already found a good latent-space solution and the MCMC chain should start near that solution rather than from diffuse random values.

Usage

als_start_vals(fit, jitter = 0, seed = NULL)

Arguments

Value

A list that can be passed to start_vals.

AME model fitting routine

Description

An MCMC routine providing a fit to an additive and multiplicative effects (AME) regression model to cross-sectional relational data of various types. This function supports both unipartite (square) and bipartite (rectangular) networks. For longitudinal networks, use the lame function. The additive and multiplicative effects framework is due to Hoff (2005, 2021).

Usage

ame(
  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"),
  family = "normal",
  intercept = !(family == "ordinal"),
  symmetric = FALSE,
  odmax = rep(max(apply(Y > 0, 1, sum, na.rm = TRUE)), nrow(Y)),
  prior = list(),
  g = NA,
  seed = 6886,
  nscan = 10000,
  burn = 500,
  odens = 25,
  verbose = TRUE,
  gof = TRUE,
  custom_gof = NULL,
  plot = FALSE,
  start_vals = NULL,
  periodic_save = FALSE,
  out_file = NULL,
  save_interval = 0.25,
  posterior_opts = NULL,
  n_chains = 1,
  cores = 1,
  use_sparse_matrices = FALSE,
  method = c("mcmc", "als"),
  bootstrap = 0L,
  bootstrap_type = c("parametric", "block"),
  bootstrap_block_length = 1L,
  bootstrap_seed = NULL,
  save_log_lik = FALSE,
  ordinal_cutpoints = c("data_induced", "explicit"),
  print,
  ...,
  model.name = NULL
)

Arguments

Details

This command provides posterior inference for parameters in AME models of cross-sectional relational data, assuming one of eight possible data types/models. The function supports both unipartite networks (square adjacency matrices) and bipartite networks (rectangular adjacency matrices with distinct row and column node sets) for single time point analysis.

Model. For a dyad (i, j) the AME linear predictor is

\eta_{ij} = \beta_0 + x_{ij}'\beta + a_i + b_j + u_i' D v_j,

where \beta are regression coefficients on dyadic / nodal covariates, a_i is a row (sender) random effect, b_j is a column (receiver) random effect, and u_i' D v_j is the multiplicative latent-factor term (rank R). The observation model is Y_{ij} \sim F(\eta_{ij}, \theta) with F specified by family (Gaussian for "normal", probit for "binary" / "cbin", etc.). For unipartite (symmetric = FALSE) the residual error has dyad-level correlation \rho between (i, j) and (j, i); for bipartite, dyad correlation is fixed at 0.

Priors (in brief). \beta has a Zellner-style g-prior (g); (a_i, b_i) are jointly Normal with covariance \Sigma_{ab} (Inverse-Wishart prior Sab0 / eta0); u_i, v_j are independent Normal with covariance \Sigma_{uv} (Inverse-Wishart prior with scale kappa0 * Suv0 and kappa0 degrees of freedom); the dyad-correlation \rho has an arc-sine prior on (-1, 1) (density proportional to (1-\rho^2)^{-1/2}), updated with Metropolis steps. See prior_summary(fit) for the priors actually used.

Choosing R. The multiplicative rank R controls the dimensionality of latent homophily / heterogeneity not explained by covariates and additive effects. R = 0 fits an additive-only social-relations model; R = 1 or 2 is typical for small / medium networks; R > floor(n/3) is rarely identifiable and will issue a warning. Latent factors capture unobserved structure (clusters, hub patterns, transitive triangles the covariates miss) and are accessed at fit$U, fit$V.

Identifiability. The latent factor term u_i' D v_j is invariant to rotation and reflection of U, V; the package canonicalises with an SVD so successive draws are interpretable. For visual stability across posterior summaries see procrustes_align and latent_positions.

Theoretical Foundation:

The AME model decomposes network structure into several components:

y_{ij} = \beta'x_{ij} + a_i + b_j + u_i'v_j + \epsilon_{ij}

where:

• \beta'x_{ij}: Fixed effects of dyadic/nodal covariates

• a_i: Additive sender (row) effect for node i

• b_j: Additive receiver (column) effect for node j

• u_i'v_j: Multiplicative interaction between latent factors

• \epsilon_{ij}: Dyadic error term (may be correlated)

This specification generalizes the social relations model (Warner et al. 1979) and latent space models (Hoff et al. 2002) within a unified framework.

Prior Distributions:

The model uses conjugate and semi-conjugate priors where possible:

• Regression coefficients: \beta \sim N(0, g\sigma^2(X'X)^{-1}) (g-prior)

• Additive effects: (a_i, b_i)' \sim N(0, \Sigma_{ab}) jointly

• Covariance: \Sigma_{ab} \sim IW(\eta_0, \eta_0 S_{ab0}) (inverse-Wishart)

• Multiplicative effects: Hierarchical shrinkage via \eta_0

• Dyadic correlation: arc-sine prior on (-1, 1), density p(\rho) \propto (1-\rho^2)^{-1/2}, with Metropolis updates

The inverse-Wishart prior on \Sigma_{ab} allows learning correlation between sender and receiver effects, capturing reciprocity patterns.

Note on the additive-effects variance scale: for a unipartite fit Sab0 defaults to a data-scaled matrix rather than a fixed diag(2). For the "normal" and "poisson" families the default is Sab0 = diag(2) * vscale, where vscale is an empirical-Bayes moment estimate of the sender/receiver variance: the mean of the variances of the centred row means and column means of Y (of log1p(Y) for "poisson"). The "binary", "ordinal", "cbin" and "frn" families apply the same idea to a probit-moment residual, but only when start_vals is not supplied and the design has at least one column; otherwise they fall back to Sab0 = diag(2) with eta0 = 4. A bipartite fit always uses Sab0 = diag(2). Call prior_summary() on a fitted object to see the prior actually used.

The data scaling matters because the inverse-Wishart prior contributes pseudo-data on the scale of eta0 * Sab0 no matter what units Y is in, so a fixed Sab0 = diag(2) pulls the additive-effects variances va/vb upward whenever the true sender/receiver variance is well below 1, most visibly at small n. In Social Relations Model simulations the data-scaled default recovers small true variances with modest bias where a fixed diag(2) prior can nearly double them; the data-scaled prior is mildly more conservative in the opposite regime, when the true variance is large relative to the residual scale.

This default differs from the amen package, whose ame(family = "nrm") leaves Sab0 = diag(2) and eta0 = 4. To reproduce amen's additive-effects posterior, pass prior = list(Sab0 = diag(2), eta0 = 4). The choice affects only va, vb and (weakly) cab; ve and rho are unchanged.

Multiplicative Effects (Latent Factors):

When R > 0, the model includes R-dimensional latent factors:

• Asymmetric case: u_i, v_j \in \mathbb{R}^R with u_i'v_j interaction

• Symmetric case: u_i = v_i with eigendecomposition ULU'

• Captures homophily, transitivity, and community structure

• R chosen via model selection or set to 2-3 for visualization

Estimation Algorithm:

The model uses a Gibbs sampler with the following updates:

1. Sample latent Z given parameters (data augmentation for non-normal families)

2. Update regression coefficients \beta via g-prior conjugate update

3. Update additive effects (a,b) jointly with \beta

4. Update covariance \Sigma_{ab} from inverse-Wishart

5. Update multiplicative effects U,V via Gibbs or Metropolis-Hastings

6. Update dyadic correlation \rho via Metropolis-Hastings

7. Update variance \sigma^2 (for continuous families)

Standard Model Types:

The following data types/models are available:

"normal": A normal AME model (identity link: E[Y] = \eta).

"binary": A binary probit AME model (probit link: P(Y=1) = \Phi(\eta)).

"ordinal": An ordinal probit AME model (cumulative probit link). An intercept is not identifiable in this model.

"cbin": An AME model for censored binary data (probit link with censoring). 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).

Value

Posterior Samples (full MCMC chains):

Posterior Means (averaged over chain):

Metadata:

Optional Posterior Samples (if requested via posterior_options):

Note on the latent-scale matrices: The posterior-mean multiplicative product is stored on the fit (UVPM, or ULUPM for symmetric fits); EZ (the expected latent network) is not stored, to save memory. Accessors:

• reconstruct_EZ(fit) - Returns linear predictor (link scale, not response scale)

• reconstruct_UVPM(fit) - Returns the stored posterior-mean multiplicative product (UVPM / ULUPM) when present, otherwise U\

Generating posterior distributions: Use simulate_posterior(fit, component="UV") to generate posterior samples for components where only means are stored, or use posterior_options() during model fitting to save full posterior samples.

Supported family x mode combinations

Every family is supported under both modes. The bipartite Z-samplers live in R/rZ_bipartite.R and dispatch per family; see also the inline comment in R/lame.R (the rectangular samplers live in R/rZ_bipartite.R).

Symmetric (symmetric = TRUE) fits require a symmetric Y. family = "ordinal" with symmetric = TRUE is supported via the dedicated sampler in R/rZ_ord_sym_fc.R, which uses the symmetric-doubled precision and mirrors upper-triangle draws to the lower triangle so Z = t(Z) holds at every sweep.

Symmetric input with one triangle missing. When symmetric = TRUE and one triangle of Y is fully NA (the user stored only the lower or upper triangle), the symmetry validator any(is.finite(Y - t(Y))) evaluates FALSE and the call proceeds: the sampler then treats the populated triangle as the symmetric data and mirrors it. This is usually intended, but if the upper / lower triangles were meant to differ, the model is silently fitting half the data. Audit anyNA(Y[upper.tri(Y)]) != anyNA(Y[lower.tri(Y)]) before calling if you are unsure.

Data preparation

ame() accepts a matrix directly. For long-format edgelists, bipartite data, and covariates, use netify to build the network object and pass that object as Y. ame() will call netify::to_lame() internally. If you already have an igraph or network object and only need a plain adjacency matrix, as_lame_y is still available as a small convenience helper.

For an undirected/symmetric network, pass symmetric = TRUE; for a rectangular two-mode network (students x courses, donors x candidates), pass mode = "bipartite" or build the netify object with mode = "bipartite". Xrow and Xcol accept either a numeric matrix or a data.frame (coerced internally); Xdyad must be a 3-D array n x n x p of numeric covariates with no Inf/NaN.

When to use AME vs ERGM

AME and ERGM are complementary tools for binary network analysis, not direct substitutes. ERGM is a class of exponential-family models built around explicit network statistics (counts of edges, mutual ties, triangles, geometrically-weighted shared partners, ...). You write the statistics you think matter, ERGM gives you their coefficients. ERGM excels when you have a substantive theory about which configurations drive tie formation.

AME models latent homophily / heterogeneity directly via sender, receiver, and multiplicative latent-factor effects. You don't enumerate triadic terms; the multiplicative-effects rank R captures higher-order structure (clustering, transitivity, hub patterns) implicitly. AME excels when (a) you have dyadic / nodal covariates whose effects you want to interpret cleanly without ERGM degeneracy, (b) higher-order structure is "nuisance" that you want to absorb but not parameterise, or (c) you need a posterior distribution over predictions for forecasting or imputation.

Practical guidance: if your research question is "do nodes that share attribute X tend to form triangles together?", reach for ERGM's gwesp. If your research question is "controlling for unobserved sender / receiver heterogeneity and latent clustering, what is the effect of dyadic covariate X?", reach for AME. The R = 0 additive-only case is the social relations model (Warner, Kenny, Stoto 1979); R >= 1 adds latent space.

ERGM to AME translation

For users coming from statnet::ergm, the rough analogues are:

For family = "binary" the link is probit, so the intercept is on the probit scale; do not compare it to an ERGM edges estimate by simple arithmetic.

Migration from amen

Both amen and lame export ame(); loading both packages fires a startup warning telling you to call lame::ame(...) or amen::ame(...) explicitly.

For default cross-sectional calls, lame::ame() follows the amen::ame() interface: the fit$BETA slot is a 2-D [n_stored, p] matrix in both packages, so scripts that call colMeans(fit$BETA) or apply(fit$BETA, 2, mean) continue to work unchanged. lame additionally accepts family = "binary" (which amen 1.4.5 no longer accepts; amen requires "bin"). The print argument is deprecated in favour of verbose; calls that pass print = ... still work but warn.

The cross-sectional path has no dynamic_beta option (an AR(1) prior on a single-period coefficient is unidentified), so the BETA 3-D shape that lame() can produce never arises from ame(). See lame for the longitudinal path and the silent-aggregation hazard with 2-D apply(fit$BETA, 2, mean) scripts under dynamic_beta = TRUE.

Notes on priors

• The regression-coefficient prior is a Zellner g-prior (beta ~ N(0, g * sigma^2 * (XtX)^-1) where XtX is the design cross-product). g is a top-level argument of ame(), not an entry in prior = list(...) (a common slip).

• There is no per-coefficient prior knob. If you need student_t, horseshoe, or to centre a slope away from zero, AME does not currently expose it; the g-prior structure is the only knob.

• Unknown names in prior = list(...) are warned about (a typo like Sab = ... instead of Sab0 would otherwise be silently dropped).

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

lame for longitudinal models, gof for post-hoc goodness-of-fit computation, gof_plot for visualizing GOF results, latent_positions for extracting latent positions as a tidy data frame, procrustes_align for Procrustes alignment of latent positions, summary.ame for model summaries, coef.ame for coefficient extraction

Examples

data(YX_bin)
fit <- ame(YX_bin$Y, Xdyad = YX_bin$X, burn = 10, nscan = 100, odens = 1,
           family = "binary", verbose = FALSE)
summary(fit)
# Note: you should run the Markov chain much longer in practice

 

Fast (MCMC-free) AME estimation for a cross-sectional network

Description

Fits an additive and multiplicative effects (AME) model to a single cross-sectional network by iterative block coordinate descent, producing a fast point estimate with no MCMC and no credible intervals.

The estimation algorithm adapts the iterative block coordinate descent estimator of the Social Influence Regression (SIR) model of Hoff & Minhas (the algorithm implemented in sir::sir_alsfit() and nicknamed "ALS" in that package) to the AME model. It is a port and adaptation, not original lame methodology.

Usage

ame_als(
  Y,
  Xdyad = NULL,
  Xrow = NULL,
  Xcol = NULL,
  R = 0,
  family = "normal",
  mode = c("unipartite", "bipartite"),
  symmetric = FALSE,
  max_iter = 200,
  tol = 1e-06,
  lowrank_method = c("mm", "als", "hybrid"),
  non_normal_method = c("irls", "transform"),
  link = c("probit", "logit"),
  linear_solver = c("eigen", "qr", "auto"),
  multistart = c("none", "cheap", "full"),
  bootstrap = 0L,
  bootstrap_type = c("parametric", "block"),
  bootstrap_block_length = 1L,
  bootstrap_seed = NULL,
  verbose = TRUE,
  seed = 6886
)

Arguments

Details

The AME decomposition

z_{ij} = \mu + \beta' x_{ij} + a_i + b_j + u_i' v_j + \epsilon_{ij}

is conditionally linear: it is a linear regression in (\mu, \beta, a, b) for fixed multiplicative factors (U, V), and the optimal rank-R (U, V) for fixed (\mu, \beta, a, b) is the truncated SVD of the residual matrix. The estimator therefore cycles, until the residual sum of squares stabilises, through three blocks — each a monotone (objective-non-increasing) update: a joint least-squares solve for the intercept \mu and regression coefficients \beta (with a ginv pseudoinverse fallback for rank-deficient designs); Gauss-Seidel sweeps to convergence for the additive effects (a, b); and a weighted low-rank update for the multiplicative factors (U, V) via SVD (eigen-decomposition when symmetric = TRUE). The weighting makes the multiplicative step correct for unbalanced longitudinal panels (dyads observed at unequal numbers of time points); severely unbalanced panels may need more iterations to converge, so raise max_iter if the fit reports non-convergence.

For family = "normal" this is a least-squares (Gaussian maximum- likelihood) fit — the exact global solution when R = 0, and a local optimum of the non-convex low-rank objective when R > 0. For "binary" and "poisson", the estimator uses an IRLS working-response path by default, so coefficients are on the requested link scale (binary probit/logit, Poisson log). The fixed-transform path is still available as a faster exploratory score, but its coefficients are not calibrated effect sizes. For R > 0 the multiplicative factor block of the IRLS families solves a penalized (MAP) sub-problem: each factor row carries the ridge implied by the AME prior [u_i, v_i] \sim N(0, \Sigma_{uv}), with the per-column prior scale estimated by an EM step (the analog of the MCMC \Sigma_{uv} draw). This is what makes the estimate well-defined: the unpenalized rank-R binary MLE does not exist under quasi-separation, and fitting it directly inflates every coefficient (a runaway that grows with R, with the latent variance diverging). With the penalty the latent variance does not run away and the slope inflation is bounded and stable in R; a residual finite-sample (incidental-parameters) upward bias remains, on the order of +10-15\ posterior mean). This residual is a property of the estimator class, not a code defect: because a nonlinear-likelihood point estimator profiles the 2 R n latent factors at their joint optimum rather than integrating them, it cannot be consistent for the regression coefficient (Neyman-Scott). The bias is dominated by rank over-specification – it is largest when R is set higher than the data support (the factors then absorb noise) and is only ~5\ ame() / lame() path integrates the factors out and is verified calibrated by simulation-based calibration, so binary R > 0 ALS is best used as a fast exploratory / rank-selection score, with ame(..., method = "mcmc") for final coefficient inference (bootstrap = N gives ALS-scale intervals but reproduces, rather than removes, this bias). Rank and ordinal families ("ordinal", "cbin", "frn") are not supported by ALS and raise an informative error; use ame or lame for those likelihoods.

Row/column (node) covariates. A node covariate broadcasts to a per-actor constant, which is collinear with the additive sender/receiver effect, so its coefficient is not identified by the objective alone. It is identified here by an explicit constraint: the additive effects are taken orthogonal to the node covariates, and beta_row/beta_col are the corresponding between-actor regression coefficients (the additive effects then carry only the residual heterogeneity). This is the standard estimand under the assumption that the additive effects are uncorrelated with the node covariates; if that assumption is doubtful the coefficient absorbs the covariate-correlated part of the additive heterogeneity. Time-varying node covariates are summarised by their per-actor mean; if a node covariate varies within actor over time the discarded within-actor variation triggers a warning.

Identifiability of the additive and multiplicative terms. For R > 0 the additive term a_i + b_j and the multiplicative term u_i' v_j are not separately identified by the objective alone: a broadcast (row- or column-constant) component can sit in either, since a pure sender effect a 1' is itself rank one. The estimator imposes the standard AME gauge — the multiplicative term is double-centered (zero row and column means), so all broadcast structure is carried by a, b — which makes the reported a, b, U, V unique given the fitted values and independent of the optimisation path. For a unipartite network the double-centering means are taken over the full matrix, which includes the structurally unobserved self-tie diagonal that the model fills in by its low-rank completion; the additive/multiplicative split therefore carries an O(1/n) dependence on that completion. On a disconnected observed-dyad graph the additive effects additionally have a per-component level shift, which is pinned to a precision-weighted minimum-norm gauge so that a, b and the variance components remain reproducible. A dyadic covariate that is itself (near) low-rank can still be partially aliased with the multiplicative term, so keep R modest relative to the covariate structure. This residual aliasing is intrinsic to the AME model — the MCMC estimator resolves it only through its priors.

Choosing R. There is no automatic order-selection criterion (the working-response objective has no likelihood, so AIC/BIC do not apply). Fit a few values — e.g. lapply(0:4, function(r) ame_als(Y, R = r, ...)) — and inspect deviance (the residual sum of squares): it falls steeply while real multiplicative signal is being captured and then flattens, so the “elbow” of that curve is a reasonable choice. With R > 0 the objective is non-convex; use multistart to guard against local optima.

Uncertainty is obtained separately, by the bootstrap; see ame_als_bootstrap. (The SIR paper's own primary standard errors are Hessian-based, classical and sandwich/robust. For AME a Hessian-based variance is awkward on two counts: without an explicit gauge fix the rotational invariance of the multiplicative factors leaves the joint Hessian rank-deficient, and even with a gauge fixed the Gaussian working-response approximation used for the non-normal families leaves a Hessian-based variance miscalibrated. The bootstrap side-steps both, so it is preferred here.)

Value

An object of class "ame_als": a list with the point estimates mu, beta, a, b, U, V (L for symmetric models), the linear predictor EZ, response-scale fitted values, working-scale residuals, convergence information, and the variance-component vector VC with five descriptive entries:

va, vb

empirical variances of the sender and receiver additive effects.

cab

covariance of the sender and receiver effects (NA for symmetric or bipartite models).

rho

dyadic residual reciprocity — the correlation between the residuals of (i,j) and (j,i) — not the sender/receiver correlation cor(a, b).

ve

residual variance, on the model's own scale and degrees-of-freedom-corrected (\mathrm{SSE}/(n_{obs} - \mathrm{df})) so it matches the MCMC posterior-mean s2. For a binary IRLS fit the probit/logit model fixes the latent error variance at 1 by identification, so ve is reported as 1; the working-scale GLM dispersion (about 1 under correct specification) is kept separately as ve_working.

These are descriptive summaries of the point estimates, not random-effect variance components. See ame_als_bootstrap for uncertainty and vcov.ame_als for a fast analytic covariance of the regression coefficients.

Coverage relative to ame / lame

The ALS estimator is a fast, frequentist point estimator. It covers most static AME workflows, and the top-level lame(..., method = "als") dispatcher covers several dynamic workflows, but posterior-specific features still require the MCMC estimator:

• Families. ALS supports normal, binary, and poisson. For ordinal, cbin, and frn, fall back to ame / lame.

• Dynamic effects. lame_als itself fits a static model pooled across time slices. The top-level dispatcher lame(..., method = "als") routes supported dynamic requests to a dynamic point estimator for normal, binary, and poisson panels, including named panels where actors enter or exit: dynamic_ab, selected intercept/dyadic/node dynamic_beta, and AR(1) or Student-t dynamic_uv for 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. bipartite dynamic_g is available on the dynamic als path for normal, binary, and poisson panels. changing actor composition requires row and column names on every slice so actors can be aligned; smoothing penalties are broken across actor-entry gaps. rank/censored dynamic families remain on the mcmc path. static ame_als() / lame_als() fits still use a single latent rank R; the dynamic bipartite ALS dispatcher honours separate R_row and R_col values.

• Priors. ALS has no priors. prior = list(...) and g = ... are MCMC-only; the dispatcher ame(..., method = "als") warns and ignores them.

• Posterior quantities. No $BETA / $VC posterior draws on the point fit. Use bootstrap = N for a sampling distribution analogue; the draws are bootstrap, not Bayesian.

• Multi-chain. Not applicable; n_chains is dropped.

• Convergence diagnostics. No Rhat / ESS / trace_plot.

Features that work the same on ALS fits: coef, vcov (sandwich on regression block), confint (auto-routes bootstrap intervals when present, sandwich Wald otherwise), predict, fitted, residuals, summary, print, nobs, simulate.ame_als, gof_plot.ame_als, ab_plot.ame_als, uv_plot, latent_positions.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

References

Minhas, S. and Hoff, P. D. (2025). Decomposing Network Dynamics: Social Influence Regression. Political Analysis. The iterative block coordinate descent estimator adapted here originates with that work (implemented in sir::sir_alsfit()).

See Also

lame_als for longitudinal networks, ame_als_bootstrap for bootstrap uncertainty, ame for the full MCMC estimator.

Examples

Y <- matrix(rnorm(400), 20, 20); diag(Y) <- NA
fit <- ame_als(Y, R = 1, family = "normal", verbose = FALSE)
coef(fit)

Bootstrap uncertainty for the fast AME estimator

Description

Computes bootstrap standard errors and percentile confidence intervals for a fast AME fit produced by ame_als or lame_als. Two strategies are available:

parametric

(default) Simulates fresh outcomes from the fitted model and refits. For normal and IRLS (binary, poisson) fits the simulation is on the calibrated response scale; for a non-normal transform fit – whose estimator is a Gaussian fit to a fixed transformed response – it is on the Gaussian working scale. The fitted residual variance is df-corrected (mean-model degrees of freedom) at fit time, so simulating with it keeps the replicates centred on the point estimate. Available for both cross-sectional and longitudinal fits.

block

Resamples time slices with replacement. With block_length = 1 slices are resampled independently; block_length > 1 draws contiguous blocks (a moving-block bootstrap), which preserves short-run temporal dependence. Longitudinal fits only (requires T > 1); with few slices its intervals are necessarily coarse.

Every replicate is refit by ame_als_refit, warm-started from the original point estimate so that replicates do not drift to different local optima. Replicates that error or return non-finite values are dropped and counted (n_valid / n_total). Standard errors are replicate column standard deviations; confidence intervals use the percentile method.

Usage

ame_als_bootstrap(
  object,
  R = 200,
  type = c("parametric", "block"),
  block_length = 1,
  seed = NULL,
  verbose = TRUE
)

boot_ame(
  object,
  R = 200,
  type = c("parametric", "block"),
  block_length = 1,
  seed = NULL,
  verbose = TRUE
)

Arguments

Details

Choice of inference. The Social Influence Regression paper of Hoff & Minhas (2025) derives its primary standard errors from the observed Hessian (classical -H^{-1} and the sandwich/robust estimator H^{-1} S H^{-1}). Two features of the AME model make a Hessian-based variance awkward here. First, without an explicit gauge fix the rank-R multiplicative term is identified only up to a full R\times R rotation/reflection (U \to U R, V \to V R^{-\top}), leaving the joint Hessian rank-deficient. Second, even with a gauge fixed, the non-normal families are fit on a Gaussian working response, so a Hessian computed from that working objective is not calibrated to the family likelihood. The bootstrap side-steps both issues and is the recommended uncertainty tool here, mirroring sir::boot_sir().

Alignment-sensitive quantities. The regression coefficients beta, the additive effects a, b and the variance components are rotation-invariant and are aggregated directly. The multiplicative factors U, V are not: each replicate is Procrustes-aligned to the original fit before its standard errors are computed. The object stores both the raw and aligned replicate factors so the effect of alignment can be inspected. The alignment is well-determined only when the multiplicative singular values are well separated; with near-equal singular values the per-column U/V standard errors reflect an unstable rotation, and the subspace — or, for symmetric models, the eigenvalues L — is the summary of record. For a symmetric model the eigenvalues L of the multiplicative term are aggregated and reported as the primary multiplicative-uncertainty summary.

The parametric bootstrap is a full parametric bootstrap of the AME model: each replicate draws fresh additive random effects a, b from their fitted dispersion (va, vb, cab) and fresh residuals, holding mu, beta and the multiplicative term fixed. Regenerating the additive effects (rather than holding the fitted a, b fixed) is what gives the intercept and node-covariate standard errors their actor-level sampling-variability component.

Assumptions. The block bootstrap treats the time slices as exchangeable replicates of the static-effects model — appropriate for that model, but not for strongly trended or serially dependent series (use the dynamic lame there); with only a few time slices it is necessarily coarse, its intervals correspondingly imprecise and somewhat anti-conservative, so parametric is the default. A variance component fit with R > 0 can still carry a small residual parametric-bootstrap bias (the low-rank refit re-absorbs simulated noise); summary() flags any point estimate that falls outside its interval. A binary IRLS fit can carry a finite-sample (incidental-parameters) bias in the point estimator itself, which the bootstrap reproduces rather than removes. The MCMC ame / lame path gives posterior summaries when that is the target.

Value

An object of class "boot_ame" with components including coefs (replicate intercept + regression coefficients), se, ci_lo, ci_hi, point_est, param_names; vc_* for the variance components; a_coefs, b_coefs, se_a, se_b; U_aligned, V_aligned, U_raw, V_raw, se_U, se_V (when R > 0); and n_valid, n_total, type, family.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

References

Minhas, S. and Hoff, P. D. (2025). Decomposing Network Dynamics: Social Influence Regression. Political Analysis. The block and parametric bootstrap follow the inference scheme of sir::boot_sir() developed for the SIR estimator.

See Also

ame_als and lame_als, both of which accept bootstrap = N, bootstrap_type, bootstrap_block_length, bootstrap_seed to do this work in a single call when you are about to fit the model anyway. ame_als_bootstrap() is the post-hoc path (bootstrap an already-fit object without refitting from scratch). ame_als_refit for warm-start refits; confint.boot_ame for the bootstrap CI extractor.

Examples

Y <- replicate(6, { m <- matrix(rnorm(225), 15, 15); diag(m) <- NA; m },
               simplify = FALSE)
fit <- lame_als(Y, R = 1, family = "normal", verbose = FALSE)
# post-hoc bootstrap of an existing fit:
bt <- ame_als_bootstrap(fit, R = 50, type = "block",
                            seed = 1, verbose = FALSE)
# equivalent one-shot call:
# fit_b <- lame_als(Y, R = 1, family = "normal", verbose = FALSE,
#                   bootstrap = 50, bootstrap_type = "block",
#                   bootstrap_seed = 1)
print(bt)

Refit a fast AME model with a warm start

Description

Refits an AME model by iterative block coordinate descent, initialised (“warm-started”) from an existing ame_als fit. This is the workhorse of ame_als_bootstrap: every bootstrap replicate is refit from the original point estimate rather than from a cold random start, which prevents replicates from converging to different local optima or rotations and is essential for meaningful bootstrap standard errors.

Usage

ame_als_refit(
  object,
  Y_new = NULL,
  X_new = NULL,
  Z_new = NULL,
  max_iter = 30,
  tol = 1e-05,
  verbose = FALSE
)

Arguments

Details

The estimation algorithm is the iterative block coordinate descent estimator of the Social Influence Regression model of Hoff & Minhas (sir::sir_alsfit()), adapted to the AME model, with each bootstrap replicate warm-started from the original point estimate.

Value

An object of class "ame_als"; see ame_als.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

References

Minhas, S. and Hoff, P. D. (2025). Decomposing Network Dynamics: Social Influence Regression. Political Analysis. The iterative block coordinate descent estimator refit here originates with that work (implemented in sir::sir_alsfit()).

See Also

ame_als, ame_als_bootstrap.

Examples

Y <- matrix(rnorm(400), 20, 20); diag(Y) <- NA
fit <- ame_als(Y, R = 1, family = "normal", verbose = FALSE)
refit <- ame_als_refit(fit, verbose = FALSE)
coef(refit)

Display memory usage information for AME models

Description

Shows estimated memory usage for networks of given size. Memory optimization is automatic, so this is informational only.

Usage

ame_memory_settings(n_nodes, R = 2)

Arguments

Details

Memory levers available in the package:

• Run compact_ame() on a fitted model to drop empty slots and, for genuinely sparse posterior means, use sparse storage via use_sparse_matrices = TRUE

• Increase odens in ame()/lame() to store fewer posterior draws

• Pass posterior_opts = list(thin_UV = ..., thin_ab = ...) to thin the stored latent-factor and additive-effect draws

Value

Invisibly returns memory estimates

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Calculate memory usage of AME model components

Description

Calculate memory usage of AME model components

Usage

ame_memory_usage(fit = NULL, detailed = TRUE)

Arguments

Value

Invisibly returns data.frame with memory usage

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

AME model fitting options

Description

Configure options for fitting AME models. Memory efficiency is now handled automatically based on network size.

Usage

ame_options(
  parallel_chains = 1,
  verbose = TRUE,
  odens = 25,
  use_sparse_matrices = FALSE
)

Arguments

Details

Memory optimization features:

• Redundant matrices (EZ, UVPM) are never stored - they can be reconstructed if needed

• use_sparse_matrices = TRUE: Converts large matrices to sparse format

• Posterior samples are always thinned appropriately

When to use sparse matrices:

• Your network has < 10\

• Memory usage is a critical concern

• You're willing to trade computational speed for memory efficiency

Note: For dense networks (most edges observed), sparse matrices will be slower and may use more memory than dense storage.

Value

List of options to pass to ame()

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples

# Configure options
opts <- ame_options(verbose = TRUE, odens = 25)
opts

Run AME model with multiple parallel chains

Description

Run AME model with multiple parallel chains

Usage

ame_parallel(
  Y,
  n_chains = 4,
  cores = n_chains,
  combine_method = c("pool", "list"),
  fitter = c("auto", "ame", "lame"),
  ...
)

Arguments

Value

If combine_method = "pool": A single ame object with pooled chains If combine_method = "list": A list of ame objects, one per chain

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

lame_multi for the unrelated multi-panel wrapper that fits K distinct networks with shared regression coefficients. ame_parallel / lame_parallel run K MCMC chains of the same model (for R-hat / ESS diagnostics); lame_multi runs one chain across K panels with pooled beta.

Examples

# Run 2 chains sequentially
data(YX_nrm)
fit_parallel <- ame_parallel(YX_nrm$Y, Xdyad = YX_nrm$X,
                             n_chains = 2, cores = 1,
                             nscan = 100, burn = 10, odens = 1,
                             verbose = FALSE)

Convert array to list.

Description

Convert array to list.

Usage

array_to_list(arrayObj, actorList, sliceLabel)

Arguments

Value

array in list format

Author(s)

Shahryar Minhas

Generic dispatcher for posterior::as_draws on lame fits

Description

Lightweight S3 generic so calls of the form as_draws(fit) dispatch through R's S3 system even when the posterior package is not loaded. When posterior is loaded, its generic of the same name is resolved first by R's namespace search; this fallback only fires for bare-namespace use.

Usage

as_draws(x, ...)

Arguments

Value

An object dispatched by the relevant method (typically a posterior::draws_array for an ame / lame fit).

Convert an AME / LAME fit to a posterior draws object

Description

Reshapes the BETA + VC posterior draws stored on an ame or lame fit into a posterior::draws_array (or draws_df-equivalent) suitable for downstream packages that consume that representation. Multi-chain fits (ame(..., n_chains = N)) are reshaped with chain as a separate dimension; single-chain fits are returned as a single-chain draws_array.

Usage

## S3 method for class 'ame'
as_draws(x, include = c("beta", "vc"), ...)

## S3 method for class 'lame'
as_draws(x, include = c("beta", "vc"), ...)

## S3 method for class 'ame_als'
as_draws(x, ...)

Arguments

Details

The reshape uses fit$chain_indicator (set when n_chains > 1) to keep chain identity separate so within-chain Rhat / ESS estimates and chain-coloured traceplots resolve correctly.

Value

A posterior::draws_array (when posterior is available) or otherwise a plain 3-D array with named dimnames.

Convert a graph object to a lame-ready adjacency matrix

Description

Convenience wrapper that takes an igraph or network object and returns the kind of matrix that ame / lame expect: a numeric adjacency matrix with NA on the diagonal (self-ties not modelled) and (if available) actor names preserved on the row/column dimnames.

Usage

as_lame_y(x, na_diag = TRUE)

Arguments

Details

Plain matrices and data.frames are accepted too; data.frames are coerced to a numeric matrix and only the diagonal is rewritten to NA.

Value

A numeric matrix suitable to pass as Y to ame / lame.

Examples

# plain matrices are returned with the diagonal set to NA
m <- matrix(rbinom(25, 1, 0.4), 5, 5,
            dimnames = list(letters[1:5], letters[1:5]))
as_lame_y(m)

if (requireNamespace("igraph", quietly = TRUE)) {
  g <- igraph::sample_gnp(8, 0.3)
  Y <- as_lame_y(g)
  dim(Y)
}

autoplot method for ALS fits

Description

Coefficient point-and-interval plot for an ame_als or lame_als fit. The body builds the same frame autoplot.lame uses, sourcing the point estimate from coef(fit) and the interval from confint(fit) (sandwich or bootstrap, depending on which is available). which = "beta" is the only mode supported; which = "uv" / "ab" return an informative error.

Usage

## S3 method for class 'ame_als'
autoplot(object, which = c("beta", "uv", "ab"), conf.level = 0.95, ...)

## S3 method for class 'lame_als'
autoplot(object, which = c("beta", "uv", "ab"), conf.level = 0.95, ...)

Arguments

Value

A ggplot object.

Ribbon plot of time-varying coefficients (or coefplot for static fits)

Description

For a lame fit with dynamic_beta on, returns a faceted ggplot of the posterior mean coefficient path per period with a 95 percent credible-interval ribbon. For a static fit (no dynamic_beta), falls back to a tidy()-driven horizontal coefplot with posterior-mean point estimate and credible-interval bars so that autoplot(fit) returns a ggplot regardless of fit type.

Usage

## S3 method for class 'lame'
autoplot(
  object,
  which = c("beta", "ab", "uv"),
  probs = c(0.025, 0.5, 0.975),
  coefs = NULL,
  ...
)

## S3 method for class 'ame'
autoplot(
  object,
  which = c("beta", "ab", "uv"),
  probs = c(0.025, 0.5, 0.975),
  coefs = NULL,
  ...
)

Arguments

Value

A ggplot2 object that can be further customised.

Examples

data(YX_bin_list)
fit <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  autoplot(fit)
}

Bipartite network helper functions

Description

Bipartite network helper functions

Usage

init_bipartite_startvals(
  Y,
  family,
  nA,
  nB,
  RA,
  RB,
  Tn,
  Xlist = NULL,
  odmax = NULL
)

Value

List of starting values for bipartite MCMC

fitted/residuals are not defined for a bootstrap object

Description

A boot_ame stores replicate-level distributional summaries, not per-cell predictions or residuals. These methods exist only to refuse the call clearly (instead of falling through to stats::fitted.default / stats::residuals.default, which would silently return NULL). Apply fitted / residuals to the underlying ame_als fit instead.

Usage

## S3 method for class 'boot_ame'
fitted(object, ...)

## S3 method for class 'boot_ame'
residuals(object, ...)

Arguments

Value

Never returns; raises an error.

Validate input data format for lame function

Description

Internal validation function that checks the format and consistency of input data for longitudinal AME models. Ensures that network data and covariates are properly formatted as lists with consistent dimensions across time periods.

Usage

check_format(Y, Xdyad=NULL, Xrow=NULL, Xcol=NULL)

Arguments

Details

Validates input data for longitudinal network analysis:

• Verifying Y is a non-empty list of matrices

• Checking that all covariates (if provided) are lists of appropriate length

• Validating data types for all elements

• Warning about dimension inconsistencies across time periods

The function uses informative error messages via the cli package to help users identify and correct data formatting issues.

Value

Invisible TRUE if all checks pass. Throws an error with an informative message if any validation fails.

Note

This is an internal function primarily used by lame() but exported for advanced users who want to validate their data before model fitting.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Extract beta path from a penalised-ALS object

Description

Extract beta path from a penalised-ALS object

Usage

## S3 method for class 'als_dynamic_beta'
coef(object, ...)

Arguments

Value

The p \times T matrix of estimates.

Extract model coefficients from AME model

Description

Returns posterior means of regression coefficients from a fitted AME or LAME model.

Usage

## S3 method for class 'ame'
coef(object, ...)

## S3 method for class 'lame'
coef(object, ...)

Arguments

Details

For a static fit (the default, and any model with dynamic_beta = FALSE), coefficients are returned as a named numeric vector computed as colMeans(fit$BETA).

For a dynamic fit (lame(..., dynamic_beta = ...) where some coefficient is time-varying), fit$BETA is a 3-dimensional array [n_stored, p, T] and coef.lame returns a [p, T] matrix of per-period posterior means. Rownames are the coefficient names; colnames are the period labels (from names(Y) or t1, t2, ...). Static coefficients in a dynamic fit are constant across the columns.

For binary models, these are on the probit (latent) scale. Use predict.ame with type = "response" to get predicted probabilities.

What coef() does not return. The multiplicative latent positions U, V are not part of the coefficient vector; they live on fit$U and fit$V (or as 3-D arrays [n, R, T] when dynamic_uv is on). The additive sender / receiver effects a, b are on fit$APM and fit$BPM. For a tidy frame of latent positions use latent_positions; for sender / receiver lollipops use ab_plot.

Value

Named numeric vector (static fit) or p x T matrix (dynamic_beta fit) of posterior mean coefficients.

See Also

vcov.ame for the posterior covariance matrix, confint.ame for credible intervals, summary.ame for a full summary table

Extract coefficients from a fast AME fit

Description

Returns the point-estimated regression coefficients (intercept first) of an ame_als fit. There is no posterior here; for uncertainty use ame_als_bootstrap.

Usage

## S3 method for class 'ame_als'
coef(object, ...)

Arguments

Value

A named numeric vector of coefficients.

Point estimates from a fast AME bootstrap

Description

Returns the intercept and regression coefficient point estimates carried by a boot_ame object (the estimates the bootstrap quantifies). For their bootstrap standard errors and intervals see confint.boot_ame and vcov.boot_ame.

Usage

## S3 method for class 'boot_ame'
coef(object, ...)

Arguments

Value

A named numeric vector of coefficient point estimates.

Cold War data

Description

Positive and negative relations between countries during the cold war

Format

A list including the following dyadic and nodal variables:

• cc: a socioarray of ordinal levels of military cooperation (positive) and conflict (negative), every 5 years;

• distance: between-country distance (in thousands of kilometers);

• gdp: country gdp in dollars every 5 years;

• polity: country polity every 5 years.

Source

Xun Cao : https://polisci.la.psu.edu/people/xuc11/

Combine multiple AME chains

Description

Combine multiple AME chains

Usage

combine_ame_chains(chain_list, diagnostics = TRUE)

Arguments

Value

A single ame object with combined chains

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Optimize AME model output for memory efficiency

Description

Optimizes AME model storage by converting matrices to sparse format.

Usage

compact_ame(fit, use_sparse_matrices = FALSE)

Arguments

Details

When use_sparse_matrices = TRUE, the posterior-mean matrices YPM and EZ are converted to sparse storage via the Matrix package, but only when a matrix's nonzero density is below 0.5 – converting a dense posterior-mean matrix would grow memory, so denser matrices are left as ordinary dense matrices. For lame fits, where YPM/EZ are per-period lists, the check and conversion are applied per element. The additive-effect summaries APM/BPM are named numeric vectors and are never converted. Names and dimnames are preserved.

Value

Memory-optimized AME model object

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Compute X'X and X'y for bipartite covariate regression

Description

Replaces the O(T x p^2 x n^2) nested R loop for bipartite XtX and Xty computation with a single C++ call.

Usage

compute_XtX_Xty_bip_cpp(Xlist, resid, p)

Arguments

Value

List with XtX (p x p) and Xty (p vector)

Compute GOF statistics for bipartite networks

Description

Compute GOF statistics for bipartite networks

Usage

compute_gof_bipartite(Y_obs, Y_sim, family)

Value

List of observed and simulated GOF statistics

Compute MCMC convergence diagnostics for multiple chains

Description

Computes Gelman-Rubin R-hat and effective sample size (ESS) across two or more independently-seeded chains, for both the regression coefficients (BETA) and the variance components (VC). ESS is computed on the pooled draws, so genuinely non-converged chains report a low ESS rather than an inflated one.

Usage

compute_mcmc_diagnostics(chain_list)

Arguments

Value

A list with rhat, ess, param_names (over BETA and VC), n_chains and n_samples.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Comtrade data

Description

Eleven years of import and export data between 229 countries. The data use the SITC Rev. 1 commodity classification, aggregated at the first level (AG1).

Format

A list consisting of a socioarray Trade and a vector dollars2010 of inflation rates. The socioarray gives yearly trade volume (exports and imports) in dollars for 10 different commodity classes for eleven years between 229 countries. This gives a five-way array. The first index is the reporting country, so Trade[i,j,t,k,1] is what i reports for exports to j, but in general this is not the same as Trade[j,i,t,k,2], what j reports as importing from i.

Source

https://comtrade.un.org/, https://www.measuringworth.com/

Bayesian credible intervals for AME model parameters

Description

Returns posterior equal-tailed quantile-based credible intervals (not highest-posterior-density, HPD). Built directly from quantile(object$BETA, c(alpha/2, 1-alpha/2)) and quantile(object$VC, ...). These are Bayesian credible intervals, not frequentist confidence intervals.

Usage

## S3 method for class 'ame'
confint(object, parm = NULL, level = 0.95, ...)

## S3 method for class 'lame'
confint(object, parm = NULL, level = 0.95, ...)

Arguments

Value

Matrix with one row per parameter and two columns (e.g. "2.5%", "97.5%").

Note on interval type

These are equal-tailed quantile intervals, not HPD. For an HPD interval use e.g. coda::HPDinterval on the columns of object$BETA and object$VC directly.

Confidence intervals for a fast AME fit

Description

Wald confidence intervals for the intercept and dyadic-covariate coefficients of an ame_als fit, built from the conditional sandwich covariance (vcov.ame_als).

Usage

## S3 method for class 'ame_als'
confint(object, parm = NULL, level = 0.95, ...)

Arguments

Details

These intervals are a fast convenience. They are conditional (the additive and multiplicative effects are held fixed) and therefore anti-conservative, and they cover only the regression coefficients the sandwich covariance is defined for – not the node-covariate, additive or multiplicative parameters. For fully-propagated intervals on all parameters, use ame_als_bootstrap and confint.boot_ame.

Value

A matrix with one row per coefficient and lower/upper bound columns.

See Also

ame_als_bootstrap for bootstrap intervals on all parameters.

Confidence intervals from a fast AME bootstrap

Description

Bootstrap confidence intervals for the intercept and regression coefficients of a fast AME fit.

Usage

## S3 method for class 'boot_ame'
confint(
  object,
  parm = NULL,
  level = 0.95,
  ci_type = c("percentile", "basic"),
  which = c("all", "beta", "vc", "a", "b", "U", "V"),
  ...
)

Arguments

Value

A matrix with one row per parameter and lower/upper bound columns.

Assemble the dyadic design socioarray for an AME model

Description

Stacks the row, column and dyadic covariates into a single n \times n \times p array in which slice k holds the dyad-level values of the kth predictor. A sender attribute is broadcast down the rows, a receiver attribute is broadcast across the columns, and a dyadic covariate is copied verbatim. An intercept slice of ones is prepended unless it is switched off or would duplicate an already-constant predictor.

Usage

design_array(Xrow=NULL,Xcol=NULL,Xdyad=NULL,intercept=TRUE,n,warn=TRUE)

Arguments

Value

an n x n x (pr+pc+pd+intercept) 3-way array

Author(s)

lame authors

Computes the design socioarray of covariate values

Description

Computes the design socioarray of covariate values for an AME fit

Usage

design_array_listwisedel(Xrow=NULL,Xcol=NULL,Xdyad=NULL,intercept=TRUE,n)

Arguments

Value

an n x n x (pr+pc+pd+intercept) 3-way array

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Detect potential change points in a dynamic_beta posterior path

Description

For each dynamic coefficient in a dynamic_beta fit, compares the posterior distribution of the maximum scaled first-difference M = \max_t |\beta_t - \beta_{t-1}| / \sigma_\beta to a Monte-Carlo approximation of the same statistic under the AR(1) prior. The returned bf column is the tail-ratio score \Pr(M^{post} > m^*) / 0.05, where m^* is the 95\ quantile of M^{prior}. It is kept under the historical column name for compatibility, but it is not a marginal-likelihood Bayes factor. Use it to surface posterior temporal jumps that the AR(1) prior cannot comfortably accommodate.

Usage

detect_change_point(
  fit,
  coefs = NULL,
  n_prior_sims = 2000L,
  threshold_bf = 10,
  seed = NULL
)

Arguments

Value

A data frame with one row per dynamic coefficient: coef, bf (the heuristic tail-ratio score), m_post_mean (posterior mean of M), m_prior_q95 (95\ M), t_hat (period of the largest scaled jump), warn (logical, bf > threshold_bf).

Examples

data(YX_bin_list)
fit <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE)
detect_change_point(fit)

Dutch college data

Description

Longitudinal relational measurements and nodal characteristics of Dutch college students, described in van de Bunt, van Duijn, and Snijders (1999). The time interval between the first four measurements was three weeks, whereas the interval between the last three was six weeks.

Format

A list consisting of a socioarray Y and a matrix X of static nodal attributes. The relational measurements range from -1 to 4, indicating the following:

• -1 a troubled or negative relationship

• 0 don't know

• 1 neutral relationship

• 2 friendly

• 3 friendship

• 4 best friends

Source

Originally available at http://moreno.ss.uci.edu/data.html#vdb (site no longer accessible).

Summarise the implied prior on a time-varying coefficient path

Description

Draws ndraws sample paths of length n_periods from the dynamic_beta prior (given a kind, AR(1) hyperparameters, and an inverse-gamma on the innovation variance) and reports the implied distribution of common summary statistics: maximum absolute first difference, roughness (sum of squared first differences), path range, and correlation with time. Use this before fitting to confirm that your prior is loose / tight enough.

Usage

dynamic_beta_prior_summary(
  n_periods = 10L,
  ndraws = 5000L,
  kind = c("ar1", "rw1", "rw2", "matern32"),
  rho_mean = 0.8,
  rho_sd = 0.15,
  rho_lower = 0,
  rho_upper = 0.999,
  sigma_shape = 2,
  sigma_scale = 1,
  matern32_length_scale = 2,
  threshold = 1,
  seed = NULL
)

Arguments

Details

For kind = "ar1", rho_mean/rho_sd are translated to a Beta prior on the standardised (\rho - \rho_{lower}) / (\rho_{upper} - \rho_{lower}). For kind = "rw1", \rho = 1 is fixed and only \sigma_\beta^2 is sampled. For kind = "rw2", the second-difference variance is sampled and the first two values use a diffuse N(0, 10) initial prior.

Value

A list with components:

summary

Data frame with quantile rows for max-first-diff, roughness, path range, trend correlation.

prob_max_diff_gt_threshold

Empirical probability that the maximum absolute first difference exceeds threshold.

rho

Length ndraws vector of sampled rho values (NA for kinds with fixed rho).

sigma

Length ndraws vector of sampled sigma values.

paths

ndraws x n_periods matrix of sample paths.

Examples

# Default prior: AR(1) with rho_mean = 0.8, rho_sd = 0.15
s <- dynamic_beta_prior_summary(n_periods = 10, ndraws = 2000, seed = 1)
s$summary
# what's the probability that consecutive beta_t differ by more than 1?
s$prob_max_diff_gt_threshold

# Tighter prior on innovation variance
s_tight <- dynamic_beta_prior_summary(sigma_scale = 0.1, seed = 1)
s_tight$summary

Edgelist to sociomatrix

Description

Construct a sociomatrix from an edgelist. Each row of el gives the two endpoints of an edge and, optionally, a weight in a third column (unweighted edges default to 1).

Usage

el2sm(el,directed=TRUE,nadiag=all(el[,1]!=el[,2]))

Arguments

Value

a sociomatrix

Author(s)

lame authors

Examples

Y<-matrix(rpois(10*10,.5),10,10) ; diag(Y)<-NA
E<-sm2el(Y)
el2sm(E) - Y

Held-out predictive evaluation for an ame / lame fit

Description

Computes family-appropriate held-out predictive scores given a fit and a logical mask of cells to score. The function does not split the data for you: it expects you to have either (a) refit on a training subset and now want to score the held-out cells of the same matrix, or (b) have predicted probabilities you want scored. Both AUROC and PR-AUC are reported when applicable, alongside Brier and mean log density.

Usage

evaluate_heldout(y_obs, y_pred, mask, family = "binary")

Arguments

Details

Workflow. The typical pattern is: mask a random sample of dyads to NA in Y, refit (lame() handles NA internally via data augmentation), call predict(fit, type = "response"), then pass that prediction alongside the original Y and the held-out mask to this function. See the examples.

Dependencies. AUROC / PR-AUC use precrec when available; if not installed, only Brier and mean log-density are computed and a one-line note is emitted.

Value

A one-row data frame with columns appropriate to the family: n_eval, plus auroc + auprc + brier + logloss (binary / cbin), or rmse + mae + mean_logdens (normal), or mean_logdens + rmse (poisson).

Examples

set.seed(1)
n <- 25; Y <- matrix(rbinom(n*n, 1, 0.3), n, n); diag(Y) <- NA
rownames(Y) <- colnames(Y) <- paste0("a", sprintf("%02d", 1:n))
# mask 20% of dyads
mask <- matrix(FALSE, n, n)
obs_idx <- which(!is.na(Y))
set.seed(1)
mask[sample(obs_idx, floor(0.2 * length(obs_idx)))] <- TRUE
Y_train <- Y; Y_train[mask] <- NA
fit <- ame(Y_train, R = 0, family = "binary",
           burn = 15, nscan = 60, odens = 5, verbose = FALSE, plot = FALSE)
y_pred <- predict(fit, type = "response")
evaluate_heldout(Y, y_pred, mask, family = "binary")

Forward-filter / backward-sample for vec(G_t) under AR(1) state prior

Description

One Carter-Kohn sweep over the per-period observations of the form vec(E_t) = H_t g_t + eps_t, with state transition g_t = rho * g_{t-1} + eta_t.

Usage

ffbs_vecG(E_cube, U_cube, V_cube, s2, rho_G, sigma_G2)

Arguments

Value

list with G_cube (RA x RB x T) and vecG_path (p x T draws of the vectorised state).

Extract fitted values from AME model

Description

Returns the posterior mean of the network on the response scale (YPM). For binary models, these are predicted probabilities between 0 and 1. For normal models, these are predicted continuous values. For Poisson models, these are predicted counts.

Usage

## S3 method for class 'ame'
fitted(object, ...)

Arguments

Value

An n x n matrix (unipartite) or nA x nB matrix (bipartite) of fitted values on the response scale. Diagonal entries are NA for unipartite networks.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

predict.ame for predictions with type control, residuals.ame for residuals

Extract fitted values from a fast AME fit

Description

Returns response-scale fitted values: a single matrix for a cross-sectional (ame_als) fit, or a list of matrices for a longitudinal (lame_als) fit.

Usage

## S3 method for class 'ame_als'
fitted(object, ...)

Arguments

Value

A matrix or list of matrices of fitted values.

Extract fitted values from LAME model

Description

Extract fitted values from LAME model

Usage

## S3 method for class 'lame'
fitted(object, ...)

Arguments

Value

List of fitted value matrices (one per time point)

Probability-integral-transform calibration check for h-step forecasts

Description

For a fit with at least one dynamic component, evaluates how well the h-step posterior-predictive distribution covers the actually observed outcomes at the corresponding period(s). A well-calibrated forecast produces approximately Uniform(0, 1) PIT values.

Usage

forecast_pit(fit, y_future, h = NULL, n_draws = NULL)

Arguments

Details

Workflow: fit on periods 1:(T - h), forecast h steps ahead, score against the held-out T - h + 1 : T observations. Returns a numeric vector of PIT values (one per observed dyad in the held-out periods), plus a one-row summary with the Kolmogorov-Smirnov statistic against Uniform(0, 1) and a fraction-of-observed-dyads coverage diagnostic.

Families. For family = "normal", the PIT is computed analytically from the per-draw forecast variance. For discrete families ("binary", "poisson", "ordinal") the randomised PIT of Czado, Gneiting & Held (2009) is used. For rank families ("cbin", "frn") PIT is not implemented; the function returns NULL with an informational note.

Plotting. The companion plot() method on the returned object renders a histogram with a uniform reference line. Use ggplot2::ggplot(pit$pit) + geom_histogram() for custom plots.

Value

A list of class forecast_pit with components pit (numeric vector), ks_stat (KS statistic against Uniform), ks_p (KS p-value), cover_95 (fraction of PIT values in [0.025, 0.975]).

References

Czado, C., Gneiting, T., & Held, L. (2009). Predictive model assessment for count data. Biometrics, 65(4), 1254-1261.

Examples

data(YX_bin_list)
Y_train <- YX_bin_list$Y[1:3]
Y_test <- YX_bin_list$Y[4]
X_train <- YX_bin_list$X[1:3]
fit <- lame(Y_train, Xdyad = X_train,
            family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE)
pit <- forecast_pit(fit, y_future = Y_test)
pit$ks_p
plot(pit)

formula() is not defined for an ame() / lame() fit

Description

ame() and lame() use named arguments (Xrow, Xcol, Xdyad, symmetric, etc.) rather than a formula interface. The summary() output renders a pseudo-formula for diagnostic display, but it is not a valid R formula and cannot be passed back to a fitting function.

Usage

## S3 method for class 'ame'
formula(x, ...)

## S3 method for class 'lame'
formula(x, ...)

Arguments

Value

Never returns; raises an error.

Compute EZ when beta is time-varying

Description

Rebuild the (n x n x T) or (nA x nB x T) EZ cube using a per-period beta vector. Mirrors get_EZ_cpp / get_EZ_bip_cpp but accepts a (T x p) beta matrix where each row is the beta for that period.

Usage

get_EZ_dynamic_beta_cpp(
  Xlist,
  beta_full_path,
  a_mat,
  b_mat,
  U_cube,
  V_cube,
  G,
  bipartite,
  symmetric
)

Arguments

Value

n_a x n_b x T cube of EZ values.

Create design array for replicate data

Description

Create design array for replicate data

Usage

get_design_rep(Y, Xdyad, Xrow, Xcol, actorSet, intercept, n, N, pr, pc, pd)

Arguments

Value

returns list of design array values necessary for ame_repL

Author(s)

Shahryar Minhas

Get fitted object from MCMC results

Description

Get fitted object from MCMC results

Usage

get_fit_object(
  APS,
  BPS,
  UVPS,
  YPS,
  BETA,
  VC,
  GOF,
  Xlist,
  actorByYr,
  colActorByYr = NULL,
  start_vals,
  symmetric,
  tryErrorChecks,
  model.name = NULL,
  U = NULL,
  V = NULL,
  dynamic_uv = FALSE,
  dynamic_ab = FALSE,
  bip = FALSE,
  rho_ab = NULL,
  rho_uv = NULL,
  family = NULL,
  odmax = NULL,
  nA = NULL,
  nB = NULL,
  n_time = NULL,
  Y_obs = NULL,
  G = NULL,
  dynamic_beta = FALSE,
  beta_dynamic_mask = NULL,
  beta_dynamic_groups = NULL,
  rho_beta = NULL,
  sigma_beta = NULL,
  RHO_BETA = NULL,
  SIGMA_BETA = NULL,
  dynamic_rho = FALSE,
  RHO = NULL,
  rho_path = NULL
)

Arguments

Value

Fitted AME object

Author(s)

Shahryar Minhas

Get fitted object from MCMC results

Description

Get fitted object from MCMC results

Usage

get_start_vals(start_vals, Y, family, xP, rvar, cvar, R, odmax = NULL)

Arguments

Value

List of starting values for MCMC

Author(s)

Shahryar Minhas

S3 generic for glance

Description

Light-weight fallback so glance(fit) dispatches through S3 even when broom or generics is not loaded.

Usage

glance(x, ...)

Arguments

Value

A one-row data frame.

Glance method for fitted ame / lame objects

Description

One-row data frame summarising model-level statistics, in the broom idiom. Used by modelsummary::modelsummary() and similar tabling tools to populate the lower goodness-of-fit panel of a regression table.

Usage

## S3 method for class 'ame'
glance(x, ...)

## S3 method for class 'lame'
glance(x, ...)

Arguments

Value

A one-row data frame with columns:

• nobs – number of observed dyads (NA cells excluded).

• n_actors – number of distinct actors; for bipartite fits this is the row-actor count (a 12 x 9 fit reports 12).

• n_row_actors, n_col_actors – row and column actor counts for bipartite fits (NA for unipartite).

• n_periods – number of time periods (1 for ame).

• n_stored – number of stored MCMC draws.

• family – outcome family, such as "normal" or "binary".

• mode – "unipartite" or "bipartite".

• R – latent-space dimension (or max of R_row, R_col for bipartite).

• dynamic_uv, dynamic_ab, dynamic_beta – logicals; whether each component is time-varying.

• elpd_loo – leave-one-out expected log predictive density. Populated only when a loo object has been cached on the fit: fit with save_log_lik = TRUE, then attach it via fit$loo <- loo(fit). NA otherwise (calling loo(fit) alone does not modify the fit object).

Examples

data(YX_bin_list)
fit <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
            nscan = 100, burn = 20, odens = 5, verbose = FALSE)
glance(fit)

Glance method for fitted ame_als / lame_als objects

Description

One-row data frame summarising an ALS fit, compatible with broom / modelsummary. Reports observation count, actor count, latent dimension, family, mode, ALS convergence flag, and iteration count.

Usage

## S3 method for class 'ame_als'
glance(x, ...)

## S3 method for class 'lame_als'
glance(x, ...)

Arguments

Value

One-row data frame with columns nobs, n_actors, n_periods, family, mode, R, converged, iterations, se_source.

Examples

data(YX_bin_list)
Y1 <- 1 * (YX_bin_list$Y[[1]] > 0); diag(Y1) <- NA
fit <- ame_als(Y = Y1, Xdyad = YX_bin_list$X[[1]],
               family = "binary", R = 1, verbose = FALSE)
glance(fit)

Compute GOF statistics from saved posterior samples

Description

Computes goodness-of-fit statistics after model estimation by generating posterior predictive networks from the saved MCMC samples. This is useful when the model was fitted with gof = FALSE to speed up MCMC sampling, or when you want to evaluate custom GOF statistics without re-running the model.

Usage

gof(fit, Y = NULL, custom_gof = NULL, nsim = 100, verbose = TRUE)

Arguments

Details

This function requires that the model was estimated with posterior sampling of the parameters needed to generate posterior predictive datasets. Specifically, it needs:

• BETA: regression coefficients

• VC: variance components

• For models with random effects: samples of a, b

• For models with latent factors: U_samples, V_samples

To enable posterior sampling during model estimation, use: posterior_opts = posterior_options(save_UV = TRUE, save_ab = TRUE)

Computing GOF post-hoc has several advantages:

• Faster MCMC sampling (no GOF overhead)

• Can experiment with different GOF statistics without re-running model

• Can control number of posterior predictive simulations independently

Value

A matrix of GOF statistics with the same format as if gof=TRUE was used during model estimation. First row contains observed statistics, subsequent rows contain posterior predictive statistics.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

ame, lame, gof_plot

Examples

# Run model without GOF during fitting
data(YX_nrm)
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, R = 2, gof = FALSE,
           nscan = 100, burn = 10, odens = 1, verbose = FALSE)

# Compute GOF post-hoc
gof_result <- gof(fit)

Visualize goodness-of-fit statistics for AME and LAME models

Description

Plots observed network statistics against their posterior predictive distributions to assess model fit.

Usage

gof_plot(
  fit,
  type = c("auto", "static", "longitudinal"),
  statistics = NULL,
  credible.level = 0.95,
  ncol = 2,
  point.size = 2,
  line.size = 1,
  title = NULL,
  ...
)

Arguments

Details

Overview:

Goodness-of-fit (GOF) assessment is crucial for network models because standard residual diagnostics are often inadequate for capturing network dependencies. This function implements posterior predictive checking by:

1. Computing key network statistics from the observed data

2. Generating multiple networks from the model's posterior predictive distribution

3. Computing the same statistics on simulated networks

4. Visualizing the comparison to identify model inadequacies

Network Statistics Evaluated:

For unipartite (square) networks:

sd.row (Out-degree heterogeneity)

Standard deviation of row means. High values indicate substantial variation in how active nodes are as senders/initiators. If the model underestimates this, it may be missing important sender effects or covariates.

sd.col (In-degree heterogeneity)

Standard deviation of column means. High values indicate substantial variation in node popularity as receivers. Underestimation suggests missing receiver effects or popularity-related covariates.

dyad.dep (Reciprocity/Mutuality)

Correlation between Y[i,j] and Y[j,i]. Positive values indicate reciprocity (mutual ties are more likely). The AME model captures this through the dyadic correlation parameter rho. Poor fit here points to adjusting the dcor parameter.

triad.dep (Transitivity/Clustering)

Measures tendency for triadic closure (friend of a friend is a friend). Calculated as correlation between Y[i,j] and sum(Y[i,k]*Y[k,j])/sqrt(n-2). AME captures this through multiplicative effects (U,V). Poor fit suggests increasing the latent dimension R.

For bipartite (rectangular) networks:

sd.row (Type A activity variation)

Standard deviation of row means for Type A nodes. Indicates heterogeneity in how actively Type A nodes connect to Type B nodes.

sd.col (Type B popularity variation)

Standard deviation of column means for Type B nodes. Indicates heterogeneity in how popular Type B nodes are with Type A nodes.

four.cycles (Bipartite clustering)

Count of 4-cycles (rectangular paths A1-B1-A2-B2-A1). High values indicate that pairs of Type A nodes tend to connect to the same Type B nodes. Captured through bipartite multiplicative effects with appropriate R_row, R_col.

Interpretation Guide:

Histogram plots (static models):

• Dashed orange vertical line (Okabe-Ito #D55E00): observed statistic value (dual-encoded on colour and linetype)

• Grey histogram: distribution from posterior predictive simulations

• Good fit: orange line falls within the bulk of the histogram

• Poor fit: orange line in the tail or outside the distribution

Common model inadequacies and solutions:

Observed sd.row/sd.col too high

Model underestimates degree heterogeneity. Solutions: Add row/column covariates (Xrow, Xcol), enable random effects (rvar=TRUE, cvar=TRUE), or increase their variance.

Observed dyad.dep outside distribution

Reciprocity not captured well. Solutions: For positive reciprocity, ensure dcor=TRUE. For negative, consider transformation or different family.

Observed triad.dep too high

Clustering/transitivity underestimated. Solutions: Increase latent dimension R, add network covariates that capture homophily, or consider including community structure covariates.

Multiple statistics showing poor fit

Fundamental model misspecification. Solutions: Change family, add missing covariates, or consider different model class.

Mathematical Details:

The posterior predictive p-value for statistic s is:

p = P(s(Y_{rep}) \geq s(Y_{obs}) | Y_{obs})

where Y_rep is drawn from the posterior predictive distribution. Values near 0 or 1 indicate poor fit. The function visualizes the full distribution rather than just p-values for richer diagnostics.

Customization:

The function accepts custom GOF statistics through the model's custom_gof argument. These are automatically included in the plot. Custom statistics should capture network features important for your specific application.

Computational Notes:

GOF computation involves simulating multiple networks, which can be computationally intensive. The model uses the networks simulated during MCMC if gof=TRUE was specified. For post-hoc GOF, use the gof() function which generates new simulations.

Good model fit is indicated when:

• Observed statistics fall within the central 95% of posterior predictive distributions

• No systematic patterns of over/underestimation across statistics

• Custom statistics (if provided) also show adequate fit

• For longitudinal models, observed trajectories track within credible bands

Value

A ggplot2 object that can be further customized

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples

# Fit an AME model
data(YX_nrm)
fit_ame <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, gof = TRUE,
               nscan = 100, burn = 10, odens = 1, verbose = FALSE)

# Basic GOF plot
gof_plot(fit_ame)

# Plot only degree-related statistics
gof_plot(fit_ame, statistics = c("sd.row", "sd.col"))

Goodness-of-fit check for an ame_als fit

Description

Bootstrap-style analogue of the MCMC gof_plot: draws nsim simulated networks from the fitted ALS model, computes the standard network statistics (sd.rowmean, sd.colmean, dyad.dep, cycle.dep, trans.dep for unipartite; sd.rowmean, sd.colmean, four.cycles for bipartite) on each replicate, and overlays the observed value on a histogram of replicate values.

Usage

gof_plot.ame_als(fit, nsim = 100, seed = NULL, ...)

Arguments

Details

Uses simulate.ame_als for the replicates and is therefore conditional on the fitted point estimates: the noise is resampled, but mu, beta, a, b, U, V are held fixed. The MCMC gof_plot.ame also averages over posterior movement in those parameters.

Value

A ggplot (or patchwork) object.

Goodness of fit statistics

Description

Calculates goodness of fit statistics for relational data matrices, evaluating second-order (dyadic) and third-order (triadic) dependence patterns. Can handle both unipartite and bipartite networks.

Usage

gof_stats(Y, mode = NULL, custom_gof = NULL)

Arguments

Details

The function computes network statistics that capture different aspects of network structure beyond simple density. These statistics are particularly useful for:

• Model checking: comparing observed statistics to those from simulated networks

• Model selection: choosing between models that better capture network dependencies

• Descriptive analysis: summarizing key structural features of the network

For bipartite networks with square dimensions (nA = nB), you must explicitly specify mode="bipartite" to ensure correct statistics are calculated.

Missing values in Y are handled by pairwise deletion for correlations and are excluded from matrix products in triadic calculations.

Value

A named numeric vector containing goodness-of-fit statistics. For unipartite networks:

sd.rowmean

Standard deviation of row means. Measures the heterogeneity in out-degree centrality (sender effects).

sd.colmean

Standard deviation of column means. Measures the heterogeneity in in-degree centrality (receiver effects).

dyad.dep

Dyadic dependence/reciprocity correlation.

cycle.dep

Cyclic/transitive triadic dependence.

trans.dep

Transitive triadic dependence.

For bipartite networks:

sd.rowmean

Standard deviation of row means (sender heterogeneity).

sd.colmean

Standard deviation of column means (receiver heterogeneity).

four.cycles

Count of four-cycles in the bipartite network.

If custom_gof is provided, additional statistics will be included with their user-specified names.

Naming note. The columns above are the internal names used by gof_stats() and stored on fit$GOF. gof_plot exposes a shorter user-facing alias for some of them; the mapping is:

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples

data(YX_nrm) 

# Auto-detect unipartite
gof_stats(YX_nrm$Y) 

# Bipartite (rectangular) networks are auto-detected; mode can also
# be given explicitly (required when a bipartite network is square)
Y_bip <- matrix(rnorm(120), 10, 12)
gof_stats(Y_bip, mode = "bipartite")

# Custom GOF function
my_stat <- function(Y) { c(my_measure = sum(Y > 0, na.rm = TRUE)) }
gof_stats(YX_nrm$Y, custom_gof = my_stat)

Goodness of fit statistics for bipartite networks

Description

Calculates goodness of fit statistics specifically designed for bipartite networks, evaluating degree heterogeneity and higher-order dependencies.

Usage

gof_stats_bipartite(Y, warn_square = TRUE)

Arguments

Details

For bipartite networks, reciprocity and triadic closure are not meaningful concepts since edges only exist between the two node sets. Instead, this function focuses on:

• Degree heterogeneity in both node sets

• Four-cycles as the simplest higher-order dependence pattern

Value

A named numeric vector containing bipartite-specific goodness-of-fit statistics:

sd.rowmean

Standard deviation of row means. Measures the heterogeneity in out-degree from set A nodes (sender effects). Higher values indicate more variation in how active A nodes are.

sd.colmean

Standard deviation of column means. Measures the heterogeneity in in-degree to set B nodes (receiver effects). Higher values indicate more variation in how popular B nodes are.

four.cycles

Count of four-cycles (also called 4-paths or squares) in the bipartite network. A four-cycle occurs when two nodes from set A (e.g., i and k) both connect to the same two nodes in set B (e.g., j and l), forming a closed path: i->j->k->l->i. This measures the tendency for pairs of A-nodes to share multiple common B-node connections, capturing a form of clustering specific to bipartite networks. High four-cycle counts indicate that connections are not random but show patterns of shared preferences or co-occurrence. For example, in a user-item network, many four-cycles suggest that users who like one item tend to also like other items that co-occur with it.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples

# Create a random bipartite network
Y <- matrix(rnorm(10*12), 10, 12)

# Calculate GOF statistics
gof_stats_bipartite(Y)

Goodness of fit statistics for unipartite networks

Description

Calculates goodness of fit statistics specifically for unipartite (square) networks, evaluating second-order (dyadic) and third-order (triadic) dependence patterns.

Usage

gof_stats_unipartite(Y)

Arguments

Details

This function computes network statistics that capture different aspects of network structure beyond simple density. These statistics are particularly useful for evaluating how well a model captures the observed network patterns.

The dyadic dependence statistic captures reciprocity - the tendency for relationships to be mutual. The triadic statistics capture different forms of triadic closure that are common in social networks.

Missing values in Y are handled by pairwise deletion for correlations and are excluded from matrix products in triadic calculations.

Value

A named numeric vector containing five goodness-of-fit statistics:

sd.rowmean

Standard deviation of row means. Measures the heterogeneity in out-degree centrality (sender effects).

sd.colmean

Standard deviation of column means. Measures the heterogeneity in in-degree centrality (receiver effects).

dyad.dep

Dyadic dependence/reciprocity correlation. Measures the correlation between Y[i,j] and Y[j,i], capturing reciprocity patterns.

cycle.dep

Cyclic triadic dependence. Measures the tendency for directed cycles (i->j->k->i) in the network.

trans.dep

Transitive triadic dependence. Measures the tendency for transitivity (if i->j and j->k, then i->k) in the network.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples

# Create a random unipartite network
Y <- matrix(rnorm(100), 10, 10)
diag(Y) <- NA
gof_stats_unipartite(Y)

Posterior-predictive temporal-trend test

Description

For a fitted lame object, computes a network statistic (density, reciprocity, or transitivity) at each observed period, fits a least-squares linear trend on period index, and compares the observed slope to slopes from posterior-predictive replicates. The two-sided value is p_pp = 2 * min(p_up, 1 - p_up), which runs from 0 (observed slope in the extreme tail) to 1 (observed slope dead-centre). A static fit on truly trending data yields p_pp near 0; a dynamic fit that captures the trend yields p_pp near 1.

Usage

gof_temporal(
  fit,
  stat = c("auto", "density", "mean", "reciprocity", "transitivity"),
  n_rep = 500,
  seed = NULL
)

Arguments

Value

A list with

• stat (chosen statistic name),

• slope_obs (observed slope of stat on period index),

• slope_rep (length-n_rep vector of replicate slopes),

• p_pp (two-sided posterior-predictive p-value),

• stat_obs_by_t, stat_rep_by_t (per-period statistics)

Examples

data(YX_bin_list)
fit <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE)
gof_temporal(fit, stat = "density", n_rep = 50)

Initialize dynamic additive effects with AR(1) structure

Description

Initialize dynamic additive effects with AR(1) structure

Usage

init_dynamic_ab_cpp(n, Tn, rho_ab, sigma_ab, mean_a = 0, mean_b = 0)

Arguments

Value

List with initialized a and b matrices

Initialize dynamic latent positions with AR(1) structure

Description

Initialize dynamic latent positions with AR(1) structure

Usage

init_dynamic_positions(n, R, Tn, rho_uv, sigma_uv)

Arguments

Value

3D array of latent positions (n x R x Tn)

AME model fitting routine for longitudinal relational data

Description

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

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

Value

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

ame for cross-sectional models, lame_als for the fast MCMC-free point estimator, lame_snap_als for the approximate dynamic snap-shift point estimator, als_dynamic_beta for a regression-only penalised smoother on the time-varying coefficient path, lame_resume for resuming saved checkpoints, gof for post-hoc goodness-of-fit computation, gof_plot for visualizing GOF results, latent_positions for extracting latent positions as a tidy data frame, procrustes_align for Procrustes alignment of latent positions, summary.lame for model summaries, coef.lame for coefficient extraction. Long-format edgelists and covariates can be prepared with netify::netify() and passed directly to lame(), or converted explicitly with netify::to_lame(lame = TRUE).

Examples

data(YX_bin_list)
fit<-lame(YX_bin_list$Y,YX_bin_list$X,burn=5,nscan=5,odens=1,family="binary")
# you should run the Markov chain much longer than this


## 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")
dim(fit_dyn$BETA)        # [n_stored, p, T] -- 3-D when dynamic
coef(fit_dyn)            # [p, T] posterior-mean coefficient paths
confint(fit_dyn)         # per-period 95% credible intervals
summary(fit_dyn)         # prints a "Dynamic coefficients per period" block

Fast (MCMC-free) AME estimation for a longitudinal network

Description

Fits an additive and multiplicative effects (AME) model to a longitudinal (replicated) network by iterative block coordinate descent, producing a fast point estimate with no MCMC and no credible intervals. This is the longitudinal counterpart of ame_als; the estimated effects (mu, beta, a, b, U, V) are static (pooled across time), as in a non-dynamic lame fit.

The estimation algorithm adapts the iterative block coordinate descent estimator of the Social Influence Regression model of Hoff & Minhas (sir::sir_alsfit()) to the AME model; it is a port and adaptation, not original lame methodology. See ame_als for the algorithm details.

Usage

lame_als(
  Y,
  Xdyad = NULL,
  Xrow = NULL,
  Xcol = NULL,
  R = 0,
  family = "normal",
  mode = c("unipartite", "bipartite"),
  symmetric = FALSE,
  max_iter = 200,
  tol = 1e-06,
  lowrank_method = c("mm", "als", "hybrid"),
  non_normal_method = c("irls", "transform"),
  link = c("probit", "logit"),
  linear_solver = c("eigen", "qr", "auto"),
  multistart = c("none", "cheap", "full"),
  bootstrap = 0L,
  bootstrap_type = c("parametric", "block"),
  bootstrap_block_length = 1L,
  bootstrap_seed = NULL,
  verbose = TRUE,
  seed = 6886
)

Arguments

Value

An object of class "ame_als"; see ame_als.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

References

Minhas, S. and Hoff, P. D. (2025). Decomposing Network Dynamics: Social Influence Regression. Political Analysis. The iterative block coordinate descent estimator adapted here originates with that work (implemented in sir::sir_alsfit()).

See Also

ame_als, ame_als_bootstrap, lame for the full MCMC estimator, lame_snap_als for the approximate dynamic snap-shift point estimator, als_dynamic_beta for a regression-only smoother that estimates only a time-varying \beta_t (no a, b, U, V; not a special case of this function).

Examples

Y <- replicate(4, { m <- matrix(rnorm(400), 20, 20); diag(m) <- NA; m },
               simplify = FALSE)
fit <- lame_als(Y, R = 1, family = "normal", verbose = FALSE)
coef(fit)

Multi-panel lame() with shared coefficients

Description

Fits K independent lame() models (one per panel) and pools the per-panel beta posteriors into a precision-weighted shared posterior. Returns a list with the per-panel fits, the pooled beta posterior, and the panel-specific deviations.

Usage

lame_multi(Y_list, Xdyad_list, ...)

Arguments

Details

This is an R-level wrapper: it fits each panel with its own MCMC and pools the results afterwards. The pooling is exact when the panels are conditionally independent given beta, which is the standard assumption.

Value

A list with

• fits: list of K per-panel lame fits.

• beta_shared: pooled posterior mean of beta (per-period when dynamic).

• beta_deviations: list of K panel-specific deviations from beta_shared.

• K: number of panels.

Class "lame_multi".

See Also

lame_parallel for the unrelated multi-chain wrapper that runs K MCMC chains of the same model (used for R-hat / ESS diagnostics and pooled effective sample size). lame_multi is for K distinct panels with shared regression coefficients; lame_parallel is for K chains of one model.

Examples

data(YX_bin_list)
fit_multi <- lame_multi(
  Y_list = list(YX_bin_list$Y, YX_bin_list$Y),
  Xdyad_list = list(YX_bin_list$X, YX_bin_list$X),
  family = "binary", R = 0,
  nscan = 100, burn = 25, odens = 5, verbose = FALSE)
dim(fit_multi$beta_shared)

Run LAME (longitudinal AME) with multiple parallel chains

Description

Thin wrapper around ame_parallel that forces fitter = "lame". Use this for longitudinal data; in particular, multi-chain dynamic_beta / dynamic_uv / dynamic_ab fits are reached through here. The combined fit's $BETA is a 3-D array when dynamic_beta is on, and posterior::rhat(as_draws(fit)) gives R-hat across chains for every per-period coefficient.

Usage

lame_parallel(
  Y,
  n_chains = 4,
  cores = n_chains,
  combine_method = c("pool", "list"),
  ...
)

Arguments

Value

Same as ame_parallel: a combined lame fit (when combine_method = "pool") or a list of fits.

Examples

data(YX_bin_list)
fit_pll <- lame_parallel(YX_bin_list$Y, Xdyad = YX_bin_list$X,
                         family = "binary", n_chains = 2, cores = 1,
                         nscan = 50, burn = 10, odens = 5,
                         dynamic_beta = "dyad", verbose = FALSE)
dim(fit_pll$BETA)        # [iter * n_chains, p, T]
fit_pll$chain_indicator  # length iter * n_chains
if (requireNamespace("posterior", quietly = TRUE)) {
  posterior::summarise_draws(posterior::as_draws(fit_pll))
}

Resume a lame() MCMC run from a checkpoint

Description

For a fit started with lame(..., checkpoint_path = "X.rds") that terminated early (either via max_seconds or an external interruption), continues the chain from the most recent checkpoint. The implementation is pragmatic: it loads the saved RNG state and re-invokes lame() with the original call arguments. The result is a fresh fit that picks up where the previous one left off in the random-number stream; the underlying MCMC counter restarts at 1.

Usage

lame_resume(path, nscan_more = NULL, ..., .envir = NULL)

Arguments

Details

Equivalent consolidated entry point. lame(resume_from = path, ...) short-circuits to this function with the user-supplied overrides forwarded; pass nscan = K on the resume call to request K additional stored draws. Both call shapes are supported. The lame_resume(path, ...) form is safer for nested calls: the consolidated form re-evaluates the saved lame() call in parent.frame(), which is the lame() frame whose required formal arguments (Y, Xdyad, etc.) were not supplied on the resume call. Prefer lame_resume(path, ...) when calling from inside other functions or from non-global scopes.

Use nscan_more = K to override the nscan value to K for the continuation. Other arguments can be overridden by passing them to ....

Value

A fitted lame object.

Examples

ck <- tempfile(fileext = ".rds")
data(YX_bin_list)
# short run with very-aggressive max_seconds to force early termination
fit1 <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
             nscan = 5000, burn = 50, odens = 5,
             checkpoint_path = ck, checkpoint_every = 50L,
             max_seconds = 0.5, verbose = FALSE)
if (isTRUE(fit1$terminated_early)) {
  fit2 <- lame_resume(ck, nscan_more = 200)
  dim(fit2$BETA)
}

Fast approximate dynamic snap-shift AME estimator

Description

Fits a fast point-estimator approximation to the dynamic snap-shift latent-factor model for longitudinal normal-valued networks. It targets the same drift-versus-reset transition estimand as lame with dynamic_uv = TRUE and dynamic_uv_kind = "snap", but returns ALS snap scores rather than MCMC draws.

Usage

lame_snap_als(
  Y,
  Xdyad = NULL,
  Xrow = NULL,
  Xcol = NULL,
  R = 2L,
  R_row = NULL,
  R_col = NULL,
  family = "normal",
  mode = c("unipartite", "bipartite"),
  symmetric = FALSE,
  max_iter = 200L,
  tol = 1e-06,
  snap_kappa = 2,
  snap_pi_prior = c(a = 1, b = 9),
  snap_update = c("soft", "hard", "annealed"),
  snap_damping = 0.7,
  estimate_rho_uv = TRUE,
  estimate_sigma_uv = TRUE,
  rho_uv = NULL,
  sigma_uv = NULL,
  hyper_update = c("robust", "em"),
  drift_quantile = 0.05,
  drift_min_transitions = 50,
  rho_prior_mean = 0.98,
  rho_prior_weight = 25,
  align = c("sequential", "global", "none"),
  threshold = 0.5,
  min_sigma = 1e-04,
  sigma_floor_fraction = 0.75,
  ridge = 1e-08,
  snap_stability_tol = 0.05,
  snap_convergence = c("quantile", "max", "classification"),
  snap_delta_quantile = 0.95,
  snap_class_change_tol = 0.005,
  unstable_top_n = 10L,
  stability = c("none", "quick", "validation"),
  verbose = TRUE,
  seed = 6886
)

Arguments

Details

The raw snap-ALS intercept mu is not separately identified from the uncentered multiplicative term: a constant level can be absorbed into \mathrm{mean}(U_t V_t') (bipartite: \mathrm{mean}(U_t G V_t')), so fit$mu on its own is arbitrary. After convergence the estimator therefore computes a post-hoc identified level: fit$mu_t stores the per-period level \mu + \mathrm{mean}(U_t V_t') (named by period) and fit$mu_identified is the mean of fit$mu_t. Interpret the network level through mu_identified; the latent trajectories and the covariate coefficients are unaffected by this shift and remain reliable outputs of this estimator.

On longer panels, the default align = "sequential" keeps the period-to-period movement in view. align = "global" can help on very short panels, but on longer panels it can pull each period back toward the pooled static fit and make early jumps look too strong. If you have a fixed drift scale, pass rho_uv and sigma_uv with estimate_rho_uv = FALSE and estimate_sigma_uv = FALSE.

The convergence output separates the usual score summary from the worst moving actor-period. By default, convergence uses the 95th percentile of score changes and the share of class changes. final_max_snap_delta and unstable_transitions still show the largest local moves. Read snap_prob as an ALS snap score. It is useful for rankings and heuristic classifications; it is not a Bayesian posterior probability.

For unipartite fits, the tracked objective (fit$objective_trace, fit$convergence$final_objective) is the exact -2\times variational free energy of the snap/drift mixture model, including the Gaussian log-normalizers, the responsibility entropy, and the beta prior on the snap rate. The exact block updates descend it monotonically; small increases can come from damped snap scores or the robust hyperparameter updates. Because it carries log-variance normalizers it can be negative and is not a deviance; fit$param_trace$sse tracks raw data fit directly.

In bipartite mode, the fitted multiplicative term is U_t G V_t' with one static interaction matrix G. The returned snap_prob matrix contains row-actor snap scores and snap_prob_v contains column-actor snap scores. Bipartite snap ALS does not estimate node covariates, dynamic coefficients, or a dynamic G_t; those combinations need a separate model.

Value

An object of class "lame_snap_als" with dynamic latent positions, snap scores, fitted values, residuals, and convergence diagnostics.

Examples

set.seed(1)
Y_bip <- lapply(seq_len(3), function(t) {
  m <- matrix(rnorm(6 * 5), 6, 5)
  rownames(m) <- paste0("r", seq_len(6))
  colnames(m) <- paste0("c", seq_len(5))
  m
})
fit_bip <- lame_snap_als(Y_bip, R = 1, mode = "bipartite",
                         max_iter = 3, verbose = FALSE)
dim(fit_bip$snap_prob_v)

Extract latent positions as a tidy data frame

Description

Extracts multiplicative latent factor positions (U and V) from a fitted ame, lame or ame_als model and returns them as a tidy data frame suitable for plotting and analysis. Optionally applies Procrustes alignment for dynamic models and includes posterior standard deviations when posterior samples are available.

Usage

latent_positions(object, ...)

## S3 method for class 'ame'
latent_positions(object, align = FALSE, ...)

## S3 method for class 'lame'
latent_positions(object, align = TRUE, ...)

## S3 method for class 'ame_als'
latent_positions(object, align = FALSE, ...)

Arguments

Value

A data frame with columns:

actor

Character. Actor name (from rownames of U or V).

dimension

Integer. Latent dimension index (1 to R).

time

Character. Time period label. Dynamic fits use the time labels from the input; static (cross-sectional) fits return "1" for every row so downstream filtering by time behaves the same in both cases.

value

Numeric. The posterior mean latent position.

posterior_sd

Numeric. Posterior standard deviation of the latent position, or NA if posterior samples are not available. To enable, fit the model with posterior_opts = posterior_options(save_UV = TRUE).

type

Character. "U" for sender/row positions, "V" for receiver/column positions. Symmetric models have only "U".

Returns a zero-row data frame with correct column names if R = 0.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

procrustes_align for standalone Procrustes alignment, uv_plot for visualizing latent positions, posterior_options for enabling posterior sampling of U/V

Examples

data(YX_nrm)
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, R = 2,
           burn = 5, nscan = 5, odens = 1, verbose = FALSE)
lp <- latent_positions(fit)
head(lp)

Lazega's law firm data

Description

Several nodal and dyadic variables measured on 71 attorneys in a law firm.

Format

A list consisting of a socioarray Y and a nodal attribute matrix X.

The dyadic variables in Y include three binary networks: advice, friendship and co-worker status.

The categorical nodal attributes in X are coded as follows:

• status (1=partner, 2=associate)

• office (1=Boston, 2=Hartford, 3=Providence)

• practice (1=litigation, 2=corporate)

• law school (1=Harvard or Yale, 2=UConn, 3=other)

seniority and age are given in years, and female is a binary indicator.

Source

Originally available at http://moreno.ss.uci.edu/data.html#lazega (site no longer accessible).

Exact rolling-origin leave-future-out cross-validation

Description

For a fitted lame object with T periods, refits the model on the first t - 1 periods (for each t in periods) and computes the expected log predictive density (elpd) of period t under the refit. Returns the summed elpd across all leave-out periods, along with per-period contributions.

Usage

lfo(fit, periods = NULL, refit = TRUE, ...)

Arguments

Value

A list with elpd_lfo (total summed elpd), pointwise (per-dyad log-density at each leave-out period; a list of numeric vectors, one per period – unlist(pointwise) gives a flat vector suitable for loo::loo_compare()-style stacking), p_lfo (effective number of parameters), per_period (data frame with period, elpd, n_obs), and periods (the periods evaluated).

Examples

data(YX_bin_list)
fit <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE)
lfo_res <- lfo(fit, periods = 4L, refit = TRUE,
               nscan = 50, burn = 10, odens = 5, verbose = FALSE)
print(lfo_res)

Convert list to array

Description

Convert list to array

Usage

list_to_array(actors, Y, Xdyad, Xrow, Xcol)

Arguments

Value

transforms Y, Xdyad, Xrow, and Xcol to arrays

Author(s)

Shahryar Minhas

Convert bipartite list data to array format

Description

Convert bipartite list data to array format

Usage

list_to_array_bipartite(
  rowActorSet,
  colActorSet,
  Y_list,
  Xdyad = NULL,
  Xrow = NULL,
  Xcol = NULL
)

Value

List with components Y, Xdyad, Xrow, Xcol in array format

Log-likelihood is not directly exposed for ame() / lame() fits

Description

ame() and lame() produce a posterior sample, not a maximum- likelihood point. A pointwise log-likelihood is computable from the posterior draws but is not stored on the fit object, so logLik() (and the AIC / BIC generics that dispatch through it) error out informatively rather than return a misleading number.

Usage

## S3 method for class 'ame'
logLik(object, ...)

## S3 method for class 'lame'
logLik(object, ...)

Arguments

Details

For Bayesian model comparison use posterior-predictive checks via gof / gof_plot, or compute WAIC / LOO yourself from the per-draw log-likelihoods (e.g. via the loo package on the BETA / VC chains).

Value

Never returns; raises an error.

Log-likelihood is not defined for a fast AME fit

Description

The fast estimator minimises a (working-response) least-squares objective, not a family likelihood, so it has no log-likelihood and AIC / BIC are undefined. Calling logLik() raises an informative error rather than returning a misleading number.

Usage

## S3 method for class 'ame_als'
logLik(object, ...)

Arguments

Value

Never returns; raises an error.

Numerically stable log(Phi(hi) - Phi(lo))

Description

Computes log of the standard-normal CDF difference. Naive log(pnorm(hi) - pnorm(lo)) returns -Inf whenever both arguments are above ~7 or below ~-7 (catastrophic cancellation). Uses pnorm(., log.p = TRUE) and a stable log-subtraction.

Usage

log_phi_diff(hi, lo)

Arguments

Value

numeric, log(Phi(hi) - Phi(lo))

Generic dispatcher for loo / waic on ame / lame fits

Description

Lightweight S3 generics so calls of the form loo(fit) dispatch through R's S3 system even when the loo package is not loaded. When loo is loaded, its generic resolves first; this fallback only fires for bare-namespace use.

Usage

loo(x, ...)

waic(x, ...)

Arguments

Value

A loo / waic object.

Approximate leave-one-out cross-validation for AME / LAME fits

Description

S3 method for loo that uses the per-iteration pointwise log-likelihood stored on the fit object (fit$log_lik) when the model was fit with save_log_lik = TRUE. Returns the standard loo object with Pareto-k diagnostics.

Usage

## S3 method for class 'ame'
loo(x, ...)

## S3 method for class 'lame'
loo(x, ...)

## S3 method for class 'ame_als'
loo(x, ...)

Arguments

Details

What log_lik measures. For family in {normal, binary, cbin, poisson, 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 loo() output from Stan / brms fit to the same family. For the rank likelihood frn the exact marginal needs GHK Monte Carlo (Halton sequence); on the longitudinal lame() path you can opt in with log_lik_method = "observed_ghk", on the cross-sectional ame() path the fallback is the augmented-Z normal approximation (with a one-time warning). Inspect fit$log_lik_method on any fit to see which branch was used.

Chunked log-lik portability. When fit with save_log_lik = "chunked", the on-disk chunk files default to tempdir(), which is cleared at the end of the R session. If you intend to saveRDS() the fit and reload it in a fresh session, supply an explicit persistent log_lik_path (e.g. "./loglik_chunks") so the chunks survive the round trip.

Value

A loo object.

Examples

data(YX_nrm)
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, R = 0,
           nscan = 60, burn = 15, odens = 5,
           save_log_lik = TRUE, verbose = FALSE)
if (requireNamespace("loo", quietly = TRUE)) {
  loo_res <- loo::loo(fit)
  print(loo_res)
}

Symmetric square root of a matrix

Description

Computes the symmetric square root of a positive definite matrix

Usage

mhalf(M)

Arguments

Value

a matrix H such that H^2 equals M

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Number of observed dyads in an AME / LAME fit

Description

Number of finite (non-missing) cells in Y, summed over time slices for longitudinal fits. Useful as a denominator for sample-size reporting and as the basis for AIC/BIC if you also have a log-likelihood (which ame() / lame() do not directly expose; see logLik.ame).

Usage

## S3 method for class 'ame'
nobs(object, ...)

## S3 method for class 'lame'
nobs(object, ...)

Arguments

Value

Integer: number of observed dyads in Y.

Number of observed dyads in an ame_als fit

Description

Mirrors nobs.ame so the MCMC and ALS fits expose the same accessor.

Usage

## S3 method for class 'ame_als'
nobs(object, ...)

Arguments

Value

Integer: number of observed dyads.

ERGM-style covariate helpers for ame() / lame()

Description

nodematch(x) returns an n x n matrix with 1 where x[i] == x[j] and 0 otherwise – the AME analogue of ERGM's nodematch term, suitable for Xdyad as a single homophily covariate. absdiff(x) returns the absolute difference |x[i] - x[j]|, the AME analogue of ERGM's absdiff. nodefactor(x) returns the dyadic matrix of x[i] + x[j] for a numeric x (or, for a factor / character x, a named list with one matrix per level giving the dyadic co-membership count (x[i] == level) + (x[j] == level), with entries 0, 1, or 2, the ERGM nodefactor semantics).

Usage

nodematch(x, na_diag = TRUE)

absdiff(x, na_diag = TRUE)

nodefactor(x, na_diag = TRUE)

Arguments

Details

All helpers return a matrix or list of matrices that is ready to wrap into an n x n x p Xdyad array via simplify2array or array().

Value

An n x n numeric matrix (or, for a factor x passed to nodefactor, a named list of such matrices).

Examples

n <- 12
grp <- sample(letters[1:3], n, replace = TRUE)
age <- rnorm(n, 40, 10)
# build a single same-group homophily covariate
Xdyad <- array(nodematch(grp), dim = c(n, n, 1),
               dimnames = list(NULL, NULL, "same_group"))
# combine homophily + age difference
Xdyad <- array(c(nodematch(grp), absdiff(age)),
               dim = c(n, n, 2),
               dimnames = list(NULL, NULL, c("same_group", "age_diff")))

Post-MCMC per-actor time-varying slopes

Description

Computes a smoothed per-actor time-varying slope coefficient on a slice of the combined dyadic design array (fit$Xlist). For kind = "row", each row-actor i gets a length-T slope path \beta_{i,t} on the across-column mean of design slice covariate_idx, fit by ridge-penalised least squares on the residual after the main MCMC linear predictor; kind = "col" uses the across-row mean for each column-actor.

Usage

per_actor_slopes(fit, kind = c("row", "col"), covariate_idx = 1L, lambda = 1)

Arguments

Value

A list with slopes (an n_{actors} \times T matrix), kind, covariate_idx, lambda, and label. Class "per_actor_slopes".

Examples

data(YX_bin_list)
fit <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
            nscan = 100, burn = 25, odens = 5, verbose = FALSE)
# post-MCMC per-row-actor slopes on the first dyadic covariate
pas <- per_actor_slopes(fit, kind = "row", lambda = 1)
dim(pas$slopes)

Simple diagnostic plot for AME model fit

Description

Provides a quick visual summary of an AME model fit, focusing on key convergence diagnostics. For more detailed diagnostics, use the specific plotting functions: trace_plot(), gof_plot(), ab_plot(), and uv_plot().

Usage

## S3 method for class 'ame'
plot(x, ...)

Arguments

Details

By default, this function simply calls trace_plot() to show MCMC trace plots and posterior distributions for key parameters. This provides a quick check of model convergence and mixing.

For more detailed visualizations, use the specialized functions:

trace_plot()

MCMC diagnostics and posterior distributions

gof_plot()

Goodness-of-fit assessment

ab_plot()

Additive sender/receiver effects

uv_plot()

Multiplicative latent factors

Value

A ggplot2 object from trace_plot() (invisibly)

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

ame, trace_plot, gof_plot, ab_plot, uv_plot

Examples

# Fit an AME model
data(YX_nrm)
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, R = 2, gof = TRUE,
           nscan = 100, burn = 10, odens = 1, verbose = FALSE)

# Quick diagnostic plot (shows trace plots)
plot(fit)

Plot the convergence of a fast AME fit

Description

Plots the block coordinate descent deviance/SSE history – a convergence diagnostic for an ame_als or lame_als fit. For substantive plots see uv_plot, ab_plot and gof_plot.

Usage

## S3 method for class 'ame_als'
plot(x, ...)

Arguments

Value

x, invisibly.

Plot diagnostics for a LAME model fit

Description

Creates diagnostic plots for a LAME (Longitudinal Additive and Multiplicative Effects) model, including MCMC diagnostics, parameter evolution over time, and longitudinal goodness-of-fit checks. This is the default plot method for LAME objects.

Usage

## S3 method for class 'lame'
plot(
  x,
  which = c(1, 2, 3, 4),
  time.points = NULL,
  ask = FALSE,
  pages = c("single", "multiple"),
  ...
)

Arguments

Details

The function produces a multi-panel plot containing:

MCMC trace plots

Shows mixing and convergence of key parameters

Posterior distributions

Density plots of regression coefficients and variance components

Longitudinal GOF

Time series of observed network statistics with posterior predictive intervals

Effects over time

Evolution of additive effects across time periods (if applicable)

Network snapshots

Visualization of network at selected time points

The plot adapts to the longitudinal structure:

• Shows temporal trends in network statistics

• Highlights composition changes if actors enter/exit

• Displays credible intervals for time-varying statistics

Value

NULL (invisibly). Plots are displayed as side effects.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

lame, trace_plot, gof_plot, ab_plot, uv_plot

Examples

# Create simple longitudinal network data
set.seed(6886)
n <- 10
nms <- paste0("n", 1:n)
Y_list <- list(
  matrix(rnorm(n * n), n, n, dimnames = list(nms, nms)),
  matrix(rnorm(n * n), n, n, dimnames = list(nms, nms))
)
diag(Y_list[[1]]) <- diag(Y_list[[2]]) <- NA
fit <- lame(Y_list, family = "normal",
            nscan = 50, burn = 10, odens = 1, verbose = FALSE, plot = FALSE)

# default combined plot
plot(fit)