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)
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, burn = 10, nscan = 100, odens = 1,
           family = "normal", verbose = FALSE, chains = 1)
summary(fit)

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 cell is treated as an unobserved tie and handled by data augmentation. This helper sets the affected Y cells to NA: the row of a sender whose covariate is missing, the column of a receiver whose covariate is missing, and each cell with a missing dyadic covariate. ame uses it for unipartite and bipartite networks; lame does the same per period (see get_design_rep).

Usage

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

Arguments

Value

Y with covariate-missing cells set to NA.

Map a coefficient vector onto a design array's slice layout

Description

symmetric fits store one merged nodal coefficient per covariate (the sampler averages the row and column coefficients at storage time), so colnames(fit$BETA) is shorter than dim(fit$X)[3] and a positional contraction of beta against the design pairs coefficients with the wrong slices. this returns a coefficient per design slice, matched by name, with the nodal coefficient <v>_node used for both <v>_row and <v>_col. a vector that already matches the design slice for slice is returned unchanged.

Usage

.beta_on_design(beta, design_names)

Arguments

Value

a numeric vector of length length(design_names).

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

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.

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

Examples

set.seed(1)
n <- 12
X <- array(rnorm(n * n * 2), dim = c(n, n, 2))
beta <- c(0.5, -0.3)
EZ <- Xbeta(X, beta)
dim(EZ)

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", "both"),
  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)
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X,
           nscan = 100, burn = 10, odens = 1, verbose = FALSE, chains = 1)

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

Examples

set.seed(1)
n <- 12
Y <- matrix(rnorm(n * n), n, n); diag(Y) <- NA
fit <- ame_als(Y, R = 1, family = "normal", verbose = FALSE)
ab_plot(fit, effect = "sender")

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 254 x 254 sociomatrix (Y) and a 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: A 254 x 3 matrix of student attributes, one row per actor of Y in the same order, 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; rows "intercept" and then the names of the covariate slices of Xdyad, or "x1", "x2", ... when the slices are unnamed; columns "t1", ..., "tT"), 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. For a dynamic ALS fit it also carries rho_uv, sigma_uv, rho_ab, and sigma_ab, so the MCMC run starts its AR(1) hyperparameters from the ALS solution. The dyadic rho start is the single pooled correlation in fit$VC["rho"], matching ALS simulation. Period-specific rho_path / residual_rho_path values are residual diagnostics, including in legacy objects labelled dynamic_rho = TRUE, and are not used as covariance-path starting values. A MCMC model requesting dynamic_rho = TRUE expands this scalar initial value; explicit MCMC path starts remain available through its start_vals argument. For a penalised dynamic fit (bias_correct = "none" or "jackknife") read these as starting values rather than estimates: the persistences are the values the ALS smoothness penalty used (the prior means unless prior set them), and the innovation scales are implied by the fitted paths, which the penalty shrinks, so they run below the posterior scale. For a dynamic empirical-Bayes fit (the default for normal, binary and Poisson panels) they are estimates under the lame() hyperpriors. For a binary fit the intercept is shifted so that the implied tie rate, mean(pnorm(EZ)), equals the observed one, with everything else held at the ALS values; the shift is stored in attr(, "level_shift"). The ALS intercept is a working-response level: under the default IRLS it sits within a few hundredths of the level the sampler has to reach, and under non_normal_method = "transform" it can sit a probit unit or more above it, the sampler's slowest direction.

Examples

set.seed(1)
n <- 12
Y <- matrix(rnorm(n * n), n, n); diag(Y) <- NA
fit_als <- ame_als(Y, R = 1, family = "normal", verbose = FALSE)
sv <- als_start_vals(fit_als)
# chains = 1 keeps the example short; the default runs two chains and compares them
fit_mcmc <- ame(Y, R = 1, start_vals = sv,
                nscan = 20, burn = 5, odens = 1, verbose = FALSE, chains = 1)

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 = NULL,
  chains = .lame_chains_default(),
  use_sparse_matrices = FALSE,
  method = c("mcmc", "als"),
  bootstrap = 0L,
  bootstrap_type = c("parametric", "block"),
  bootstrap_block_length = 1L,
  bootstrap_seed = NULL,
  multistart = c("none", "cheap", "full"),
  save_log_lik = FALSE,
  ordinal_cutpoints = c("data_induced", "explicit"),
  print,
  ...,
  log_lik_method = c("observed_exact", "observed_approx"),
  model.name = NULL,
  bias_correct = NULL
)

Arguments

Details

This command provides posterior inference for parameters in AME models of cross-sectional relational data, assuming one of six 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 the prior 2\,\mathrm{Beta}(2, 2) - 1 on (-1, 1) (density proportional to 1-\rho^2, the prior lame uses), 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. A rank above floor(n/3) is fitted like any other (on 18-actor panels with a rank-2 truth, ranks 6 to 17 ran without numerical failure; the chain mixed more slowly and the fitted values lost accuracy gradually as the rank grew past the truth), so compare ranks by their fit rather than by a cap. 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 unipartite fits with both resolved actor margins enabled, reporting centers static latent positions and transfers the removed row, column and grand-mean components to the actor effects and intercept. A symmetric shared actor absorbs half the grand constant in each margin when no intercept is fitted. Restricted unipartite models retain sampled coordinates in their stored factor draws. This reporting map preserves each draw's predictor and leaves the sampler state unchanged. 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: 2\,\mathrm{Beta}(2, 2) - 1 prior on (-1, 1), density p(\rho) \propto 1-\rho^2 (the prior of lame), 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 scales Sab0 the same way (one scale for both sides). 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. The data-scaled default sets the prior on the scale of Y; prior_summary() prints it.

To use a fixed identity scale instead, pass prior = list(Sab0 = diag(2), eta0 = 4). Changing this prior can affect other jointly estimated quantities as well as va, vb and cab. Matching these prior settings alone does not establish equivalence to another package's posterior.

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 with data-induced cutpoints. Explicit cutpoints anchored at zero support an intercept.

"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-lognormal model for counts (negative or non-integer observed values are an error): Y \sim \mathrm{Poisson}(e \exp(z)) with z \sim N(\eta, \sigma^2). Given \eta and \sigma^2, a cell's mean count is e \exp(\eta + \sigma^2/2). YPM, fitted() and predict(type = "response") report it at the posterior mean \eta of the linear predictor and the posterior mean \sigma^2 of the residual variance, the form the fast estimators report. Predictive quantiles and event probabilities come from posterior predictive draws (simulate(), predict(type = "distribution")). Zero exposure gives a structural zero count (cross-sectional fits use unit exposure).

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 and symmetric Xdyad slices: the regression term of a symmetric fit is symmetric, so an asymmetric slice is an error (class lame_symmetric_xdyad) that names the slice. Replace it by (X + t(X)) / 2. A value missing in one triangle of a slice is not an asymmetry: it is filled from the other triangle (see Missing covariates). lame applies the same checks to every period, and the ALS routes (method = "als", ame_als) refuse an asymmetric slice the same way. family = "ordinal" with symmetric = TRUE is supported via the unit-variance unordered-dyad samplers for either cutpoint convention, which mirror upper-triangle draws to the lower triangle so Z = t(Z) holds at every sweep. The rank-nomination families ("cbin", "frn") model each actor's own ranked nominations and are directed by construction, so symmetric = TRUE is rejected for them.

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 check any(is.finite(Y - t(Y))) is FALSE, and the call stops with an error instead of mirroring the populated triangle. Symmetrize Y explicitly (for example Y[lower.tri(Y)] <- t(Y)[lower.tri(Y)]) before calling; lame() applies the same check to every period.

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.

Directed normal MCMC models with estimated reciprocity use the prior 2\,\mathrm{Beta}(2, 2) - 1 on the full open interval (-1,1), with density proportional to 1 - \rho^2. If every observed reciprocal response pair is exactly equal, the posterior can be improper at \rho = 1. A conservative structural check counts reciprocal pairs and bounds the antisymmetric actor/regression design dimension. With multiplicative rank R, it also includes a conservative n * R factor bound. The check refuses cases where the pair count exceeds that combined bound by four or more (the posterior is improper from that excess under this prior; under the arc-sine prior of versions before 1.4.0 it was improper from an excess of one); it does not certify propriety of other designs. For an undirected network, use symmetric = TRUE; to fix reciprocity at zero in a directed model, use dcor = FALSE. An exactly symmetric Y fitted with symmetric = FALSE is allowed: the two cells of a dyad are equal, so the residuals of reciprocal cells are equal and the posterior of the dyadic correlation concentrates near one (with dcor = FALSE the two cells count as two observations). A message says so and points at symmetric = TRUE, which fits the undirected model.

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.

Missing covariates

A cell whose dyadic covariate, sender (row) covariate or receiver (column) covariate is missing is an unobserved tie, for unipartite and bipartite networks alike: it is set to NA in Y before the design is built, the sampler augments it like any other missing cell, and the covariate coefficients are estimated from the complete cells. A message gives the number of such cells. A missing value is never modelled as 0, which would bias the coefficient of a covariate that is not centred. In an undirected network (symmetric = TRUE) a dyad has one covariate value: a hole in one triangle of a dyadic covariate is filled from the other triangle (a message gives the number of values), so the slice stays symmetric and the dyad stays observed, and a dyad missing in both triangles is unobserved. A missing node covariate removes the actor's row and column. A covariate that is missing everywhere stops the fit. lame (per period) and the ALS routes (ame_als) apply the same rule and give the same messages.

Collinear covariates

A regression design whose columns are exactly aliased over the observed cells (a covariate that duplicates another, or one that is a linear combination of others, counting the intercept and the sender and receiver covariates) is an error of class lame_design_rank. It names the aliased columns and the combination each equals: the data cannot separate their coefficients, so the posterior along that direction would be the prior's alone. A column is aliased when it is a combination of the columns before it to a relative 1e-8 (as lm() reports them); a design that is only nearly collinear is fitted.

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

Chains and convergence policy

ame() runs chains chains (default 2; options(lame.chains = k) sets another default). Chain k is seeded with seed + 1000 * (k - 1), so chains = 1 is the first chain of a two-chain fit. On Unix the chains run forked, on at most cores workers (all cores but one by default); on Windows they run one after the other. The draws do not depend on the number of workers. The result is the pooled fit of combine_ame_chains, with fit$chains, fit$chain_seeds, fit$rhat (the R-hat of the monitored set), the state each chain stopped in, and the chain_indicator that as_draws() and trace_plot() use. Two chains cost twice the CPU time and about the same wall time on a multicore machine.

At the end of the post-burn run every chain is checked, and the chains are compared with each other. Within a chain, the first and last thirds of the stored draws of the intercept (or first coefficient), va, vb, ve and rho are compared with batch-means standard errors (ten batches per third); a quantity is flagged when the two means differ by more than 3 standard errors and by more than half a posterior standard deviation of the whole series. Between chains, the split rank-normalised R-hat (Vehtari et al., 2021) of the monitored set is compared with 1.05: every regression coefficient, va, vb, cab, ve, rho, and summaries of the additive and latent structure recorded at every stored draw (the rms of a and b, the largest actor effect, the rms of the latent product UV', the norm of G of a bipartite fit and the split of scale between U and V of a directed fit). The dyadic products are not monitored. Quantities a fit reports as not estimated are left out, and a run with fewer than 60 stored draws per chain is not checked. The gate is computed on the draws the returned fit holds (the last block of every chain), and summary() and compute_mcmc_diagnostics report the same statistic of the same draws: the Rhat column of summary() is fit$rhat for every monitored coefficient and variance component.

When a chain drifts or an R-hat exceeds 1.05, every chain is extended by a block of nscan sweeps from its own final state (the seeds advance by the sweeps already run); the flagged block counts as extra burn-in, its draws are discarded, and both checks run again. Up to prior$chain_extensions extensions are made (default 3, so at most four times nscan post-burn sweeps per chain). The returned fit has the same number of stored draws as a fit that needed no extension, and fit$auto_extend records the extensions, the sweeps discarded and the quantities that triggered each. If the last block is still flagged, ame() stops with an error of class lame_chain_unsettled: its fits field holds the chains, fit the pooled fit of the last block, rhat the monitored table and moving the offending quantities, and lame_resume(e$fit, nscan = nscan) continues every chain (lame_resume(e$fits[[k]], nscan = nscan) one). prior = list(chain_extensions = 0) refuses at the first flag; chains = 1 runs one chain with the drift check only. There is no mode that returns a flagged fit. The measured false-extension rate and power are in lame.

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_nrm)
# the default runs two chains and compares them on the regression
# coefficients and variance components
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, burn = 100, nscan = 600, odens = 10,
           verbose = FALSE)
summary(fit)
fit$rhat  # between-chain split R-hat; every value is below 1.05
# the defaults (burn = 500, nscan = 10000, odens = 25) are the length of a final run

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. For normal, binary and Poisson data the default is an empirical-Bayes / variational estimator of the model that ame samples (see bias_correct); its vcov() and confint() are Laplace approximations.

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"),
  bias_correct = NULL,
  bootstrap = 0L,
  bootstrap_type = c("parametric", "block"),
  bootstrap_block_length = 1L,
  bootstrap_seed = NULL,
  verbose = TRUE,
  seed = 6886,
  dcor = !symmetric,
  prior = NULL
)

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 (an exactly rank-deficient design is refused with an error of class "lame_design_rank" that names the aliased columns); an exact solve of the additive effects (a, b) (Gauss-Seidel sweeps started from the direct solution of their normal equations, and an exact eigen solve wherever the sweeps do not reach the normal equations, as on a weakly connected observed-dyad graph; a residual above the tolerance after the exact solve stops the fit with an error of class "lame_als_numerical"); 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 converge in more iterations. A fixed-effects fit that has not reached its criterion at max_iter is refitted once with four times that limit, and refused with an error of class "lame_als_convergence" if it still has not: the message gives the stop (iteration limit, or no step that lowers the objective), the score against its threshold, and the remedy. The last iterate is not returned. An empirical-Bayes fit that does not settle is refused with the same class at once, with the state it stopped in and the remedy.

For fixed-effects Gaussian fits (bias_correct = "none") and fixed-transform fits, convergence additionally requires a centered first-order SSE score no larger than sqrt(tol). A premature plateau or exhausted outer budget receives up to five joint refinements, each bounded by max_iter optimization iterations. These preserve the same objective and observation mask; a requested single sweep remains unfinished. gaussian_diagnostics records the initial and final scores, objective values and refinement counts. Multi-start fits retain these checks for every start in multistart_gaussian_diagnostics; warm refits and the Gaussian bootstrap use the same certificate. A first-order certificate attests a stationary point, a local optimum of the non-convex factor model when R > 0.

With bias_correct = "none" a family = "normal" fit 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; the default for a normal response is the empirical-Bayes estimator of Normal data (the Gaussian empirical-Bayes estimator) below. 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.

Empirical-Bayes estimator (the default of binary, Poisson and normal fits). A binary or Poisson IRLS fit, or a normal fit, with bias_correct = NULL (or "eb") treats the sender and receiver effects as the MCMC ame / lame do, as draws from a Gaussian prior (a_i, b_i) \sim N(0, \Sigma_{ab}) (an undirected fit has one effect per actor and a scalar variance, a bipartite fit independent sender and receiver variances), instead of as free fixed effects, and estimates the model by empirical Bayes in two stages. The first is the posterior mode with a Laplace correction, in three parts. (1) The IRLS working problem carries the penalty \sum_i \alpha_i' \Sigma_{ab}^{-1} \alpha_i, where \alpha_i is the actor's effects less the part its node covariates explain (the node-covariate coefficients are fixed effects of the prior mean, as in the MCMC model); the factor block is the penalised ridge above, whose scales are updated by EM from the posterior second moments of the factors, an undirected fit's too. At every reweighting step the intercept, the coefficients and all the effects are solved jointly and exactly from dyad-level sums, in one dense system with no centring step. (2) \Sigma_{ab} is updated by an EM step from the effects and their Laplace covariance, under the inverse-Wishart prior that lame() puts on it for binary data, so that the penalty is the MCMC prior; without that prior the update drives the variance to zero on sparse networks and with a latent term. (3) The coefficients solve the Laplace-marginal estimating equations instead of the penalised-fit ones: the score of every cell loses half the derivative of the log determinant of the effects' information (the integral over the effects). The penalised fit alone had a slope too small by up to 4.3 percent with 24 percent ties and 6.8 percent with 7 percent ties (30 networks a cell). The latent factors get the same cell correction, so one objective, the penalised deviance plus that log determinant, is descended. This stage is what fit$irls_diagnostics certifies and what the marginal likelihood fit$eb$log_marginal (see Choosing the rank) is evaluated at.

Variational posterior (the values an empirical-Bayes fit reports). The mode is the starting point of the second stage, not its result. The mode of a latent factor that the data barely identify lies at the edge of a broad posterior and is too large (it fits noise), which inflates the coefficient of a covariate that follows the product and, one rank above the generating one, the slope. The reported values (a, b, U, V, UVPM, the coefficients and the variance components) are those of a Gaussian variational posterior of the sender and receiver effects and of the factors, independent across actors and between the sender and the receiver side of the factors; the coefficients and the variance components remain empirical-Bayes point estimates. Under this posterior the predictor of a cell has a mean, the fitted predictor, and a variance, that of the effects plus that of the product, \hat u_i' \Omega_j \hat u_i + \hat v_j' \Sigma_i \hat v_j + \mathrm{tr}(\Sigma_i \Omega_j) (\Sigma_i and \Omega_j are the posterior covariances of a sender's and a receiver's factor row), and the score and the information of the likelihood are expected over that Gaussian (nine-point Gauss-Hermite quadrature) instead of taken at its mean. Every equation of the fit uses the expected score and information, and the system of a factor row adds the second moment of the other side (mean product plus covariance). The means and these covariances are the stationary point of one objective, the expected log likelihood less the Kullback-Leibler divergence to the prior, at the scale of the prior of the factors; the step of the iteration is halved when the objective would fall, which leaves the fixed point unchanged. A fit without a latent rank is variational too (the effects only). fit$eb$estimator says which estimator produced the values ("variational", or "laplace" for the mode alone), fit$eb$variational records the iteration (convergence, stop reason, number of iterations, the largest step a full update would still take, the lower bound and its parts, and the scale of the prior of the factors) and fit$eb$factor_posterior the covariances. A fit whose mode or whose iteration stops short of its fixed point is not returned: the call stops with an error of class "lame_als_convergence" that reports the state and the remedy (a smaller rank, a larger max_iter, or ame / lame); the values of an iteration that has not settled are never reported as the estimate, nor replaced by the mode.

The scale of the prior of the factors. The rows of the factors have a Gaussian prior whose covariance the data inform weakly where there are few events per loading. The fit estimates it by the variational update of the inverse-Wishart prior ame() and lame() put on it (kappa0 degrees of freedom and scale matrix kappa0 * Suv0, the data-scaled default of ame), the inverse of the posterior mean precision, (\Psi + SE) / (\mathrm{df} + n), with SE the second moments of the stacked sender and receiver rows. What second moments to use decides the result. The mean-field covariance of a row counts the uncertainty of the other side twice (the second moment of the other side in its system is the mean product plus the covariance) and is below the covariance of the exact posterior; the update evaluated there drives the scale to the floor of its hyperprior wherever the data are weak, and the product and the coefficients that follow it with it. Against an exact Gibbs sampler of the same model run at the same scale (fixed hyperparameters, 8 networks of 32 to 100 actors), the second moments of the rows were 4 to 19 percent below the sampler's at the mean-field covariances (on average 9.5). The default evaluates the update at the cavity covariances, \Sigma_i = (\Lambda + \sum_j \nu_{ij} \hat v_j \hat v_j')^{-1}, where \nu_{ij} is the cavity weight of the cell, the curvature of the log likelihood integrated over its Gaussian predictor (the information the cell adds about the row once the rest of the predictor is integrated out), and these were from 9 percent below to 5 percent above the sampler's (on average 3 below); they are the covariances fit$eb$factor_posterior reports (the covariances of the fixed point of the updates, from which the lower bound and a refit start, are kept beside them as Sig_mf and Om_mf). The covariance of the stacked rows of a directed network is isotropic, one scale and one correlation between the sender and receiver position of an actor (the sampler's prior is a full 2R \times 2R matrix: the scale carries one dimension of its uncertainty, see below); a bipartite fit has the priors of its sampler, the receiver factor at unit scale and, for the sender columns, the scales of P = U G V' with standard normal G: the sender rows seen by the likelihood have covariance G'G, a Wishart matrix with R degrees of freedom, and the scales of the columns are its ordered eigenvalues, which repel each other (the density of the eigenvalues carries the product of their differences), so a column the data say nothing about falls to the small eigenvalues of that matrix rather than staying at the scale of the prior; and an undirected fit the prior of its sampler on each column (U L U' with L_r \sim N(0, 1) and column precision \psi_r (gamma, shape 1, rate 1/2), so that a column of the factor seen by the likelihood has variance |L_r| / \psi_r). The scale is recorded in fit$eb$variational$scale (the rule, the scale, the correlation of a directed fit, and the standard deviation of its logarithm). options(lame.eb_scale = "plugin") selects the update at the mean-field covariances with a full 2R \times 2R covariance and the closed-form rescaling of the factors that it needs, the fit of the earlier release. The option, and options(lame.eb_scale_uncertainty) below, also govern the static factors of a normal, binary or Poisson dynamic fit (lame(method = "als")), not the persistence of its moving positions.

Dyadic correlation. A directed binary probit fit with the empirical-Bayes estimator and dcor = TRUE (the default for directed networks) takes the two cells of a dyad as one observation: the latent errors of y_{ij} and y_{ji} are bivariate normal with correlation \rho, and the likelihood of the pair is the bivariate normal probability \Phi_2(s_1 \eta_{ij}, s_2 \eta_{ji}; s_1 s_2 \rho), s = 2y - 1. The variational posterior is taken in that model. The expected log likelihood of a pair, its score and its information (a 2 \times 2 matrix: the cells of a dyad are not independent given the predictors) are computed under the Gaussian posterior of the two predictors, which are correlated through the effects of the two actors, by a Gauss-Hermite rule of five nodes a side. The equations of the coefficients and of the actor effects carry the cross information of the pairs, the rows of the factors see the exact score of the pair, the cavity weights that set the scale of the factors are those of the pair probability integrated over the predictors (a bivariate normal probability in closed form), and \rho is estimated with the other hyperparameters: at every iteration it moves to the posterior mean of the quadratic approximation of the expected log likelihood in \rho under the prior of the samplers, 2\,\mathrm{Beta}(2, 2) - 1 (density proportional to 1 - \rho^2), the one ame and lame share and the fast fit takes from the same source. The value is fit$VC["rho"], its posterior standard deviation and the method are in fit$rho_estimation, and the record of the pair model (the correlation, its score and curvature, the number of pairs, the quadrature) is fit$eb$pair. Against the sampler (80 simulated networks of 60 and 100 actors, one and four periods, ranks 0 and 2, with correlations of 0, 0.3, 0.6 and 0.8) the correlation is within 0.02 of the sampler's posterior mean on average and within one posterior standard deviation (0.04 to 0.1) in 79 of the 80 networks, where the post hoc correlation of the independent-cell fit is 0.02 to 0.08 too high; the slope of a covariate is within 0.2 sampler standard deviations of the sampler's on average, and its standard error 1.01 to 1.02 times the sampler's standard deviation, with or without the pair model. The factors of an independent-cell fit absorb reciprocity, which the pair model attributes to \rho. The cost is the pair quadrature, in proportion to the number of dyads and periods: 1 to 2.5 times the independent-cell fit for one period (60 to 150 actors, ranks 0 and 2) and 2 to 4.5 times for four.

Counts (Poisson). A Poisson fit with the default bias_correct (or "eb") is the same estimator for the model that ame and lame sample: a count is y_{ij} \sim \mathrm{Poisson}(e^{z_{ij}}) with z_{ij} = \eta_{ij} + e_{ij} and e_{ij} \sim N(0, v_e), where \eta is the predictor of the sections above (coefficients, sender and receiver effects, latent product). The fixed-effects Poisson fit (bias_correct = "none", the default of the earlier releases) has no finite estimate for the effect of an actor without a count on a side, which runs to minus infinity (-17.5 after 71 reweighting steps on a network of 30 actors with one empty row and one empty column), and it fits the marginal mean with a post hoc overdispersion. The empirical-Bayes fit gives such an actor a finite effect, shrunk towards the prior mean as in the sampler. The effects have the Gaussian prior of the binary fit, with the scale and the degrees of freedom that ame() and lame() give counts: the variance across actors of the centred log row and column rates (a half-count pseudo-observation keeps an empty actor finite) and \mathrm{round}(4 + 3 n / 100) degrees of freedom, n the number of actors; the factors have the priors of the binary fit. The variance of the noise v_e is a hyperparameter, estimated as the covariance of the effects is, by type-II maximum likelihood: it is the value that maximises the Laplace marginal likelihood of the model (fit$eb$log_marginal, the mode of the effects and the factors solved at that value) plus the log density of the prior that ame() and lame() give it (the inverse-gamma of one degree of freedom on the scale one). It is fit$VC["ve"] and fit$eb$ve, and the objective it was found on is fit$eb$ve_search$curve. The stationarity condition of the maximum is the EM step that ame() and lame() draw v_e from, (1 + \sum E[e^2]) / (N - 1) (the sum running over the observed cells, the dyads once for an undirected network, of the second moment of the noise given the count), with the second moment taken, to the order of the Laplace approximation, under the posterior of the predictor of the cell given all the parameters; the step at the mode alone, or under the mean-field posterior of the predictor, which is shrunk where the data inform the latent product weakly, is too large (a count network with fewer than three nonzero counts per loading: 1.3 to 4.6 posterior standard deviations above the sampler's, where the maximum of the marginal likelihood is within 1). A count network without noise beyond the Poisson has a positive estimate, as the sampler's posterior does: the inverse-gamma prior holds v_e off zero (on four plain Poisson panels of 22 actors and five periods the estimate was 0.034 to 0.046, within 0.0015 of the sampler's posterior mean; 0.040 on average over 60 panels). The maximum is found by continuation on \log v_e: the mode is solved at 0.05, where every column of the latent product is alive, and then at larger values until the objective falls, then at the vertex of the parabola through the three points about its maximum, each warm from the nearest solved mode (typically 8 to 10 modes, 5 to 30 reweighting steps each); a mode solved cold at a large v_e loses a weakly supported column and keeps it lost at every smaller value, on a lower branch of the marginal likelihood. The variational stage holds v_e at the maximum, and a column the mode lost is restarted from the score of the cells before it iterates. The cell has everything the binary cell has. Its noise is integrated out, so that the log likelihood, the score and the information of a count at the predictor \eta are those of the Poisson-lognormal distribution (the information is positive: the log likelihood is concave in \eta); they are the cumulants of the posterior of the noise given the count, computed by an adaptive Gauss-Hermite rule centred at its mode (14 to 64 nodes, by the size of v_e), with the derivative of the information that the Laplace correction needs. The variational stage takes the expectation of the cell over the Gaussian posterior of the predictor with the rule of the binary fit (nine nodes where the predictor is uncertain, fewer where its variance is below 0.3), and the cavity weight of a cell is the information of the Poisson-lognormal cell whose noise variance is v_e plus the variance of the predictor. The linear predictor EZ is on the latent log scale of the model (there is no shift by half the noise variance), the mean of a count is e^{\eta + v_e / 2}, which is what fit$fitted, fit$YPM, fitted() and predict() hold (the plug-in mean, as the samplers report it; the posterior predictive mean of a count is not finite in this model, and plug_in = FALSE returns the value the integration over the predictor gives, e^{m + (v + v_e) / 2}). A network without a count above zero has no finite log rate, and a dyadic covariate of one sign that is non-zero only where every count is zero has no finite coefficient: both are refused (class "lame_als_separation"). Rank evidence (fit$eb$log_marginal) is the Laplace value of the Poisson-lognormal model.

Dyadic correlation of counts. A directed count network with dcor = TRUE (the default for directed networks, as in ame and lame) takes the two counts of a dyad as one observation at the fitted state: the noises (e_{ij}, e_{ji}) are bivariate normal with variance v_e and correlation \rho, and the pair has the bivariate Poisson-lognormal probability. The effects, the factors and v_e are fitted with independent cells (they agree with the sampler's, see below); the pair stage then gives (1) \rho, the posterior mean, under the prior the samplers share, of the profile of the pair probability at the fitted predictors, reported as fit$VC["rho"] with its standard deviation in fit$rho_estimation, and (2) the information of the pairs: the pair probability integrated over the Gaussian posterior of the two predictors is the pair probability with the noise covariance plus their covariance, and its curvature in the two means (two diagonal weights and a cross weight, negative for a positive correlation) enters the covariance of the coefficients, so that the standard error of a covariate that is the same in the two cells of a dyad is not understated (fit$eb$pair records the correlation, its precision, the number of pairs and the cross weights). The pair stage costs 13 to 22 evaluations of the pair probability of every dyad; dcor = FALSE leaves it out and keeps the post hoc moment estimate of the correlation of the latent errors (fit$rho_working holds it in either case). Undirected and bipartite networks have no pair of cells.

Evidence for counts. Against the sampler (29 simulated count networks of 30 to 50 actors, one and four periods; directed, bipartite and undirected, ranks 0 to 3, mean counts 0.3 to 31, noise variances 0.05 to 2, empty actors, node covariates and reciprocity 0 to 0.8; ame() / lame(), 30,000 scans after the burn-in) the slopes are within 0.35 posterior standard deviations of the posterior means (median 0.05; the fixed-effects fit 2.8 and 0.6) and the intercepts within 0.38, the standard errors are 0.99 to 1.09 times the posterior standard deviations of the slopes (mean 1.04) and 0.97 to 1.02 times those of the intercept of a cross-section, va and vb are within 0.3 posterior standard deviations and ve within 0.5 (median 0.07; the fixed-effects fit 6.9 and 1.7), the actor effects differ from the posterior means by a median 0.029 of their standard deviation (the fixed-effects fit 0.17, and 2.3 with empty actors), and the correlation of a directed network is within 0.027 of the sampler's (median 0.010 and at most 1.2 of its posterior standard deviations). Without the pair stage the standard error of a slope was 0.77 to 0.87 of the sampler's for dense counts with a covariate that is the same in both cells of a dyad and a reciprocity of 0.6 to 0.8. On 1,500 simulated networks (25 designs of 60: 30 to 60 actors, mean counts 0.1 to 20, noise variances 0.1 to 1.5, ranks 0 to 3, directed, bipartite and undirected, one and four periods, reciprocal designs with a symmetric covariate, and ranks above the generating one) the fit converged in every network (the fixed-effects fit in 94.9 percent, 67 percent at worst), the slope had a bias of 0.6 percent on average (4.7 at most, with 30 actors, mean count 0.1 and noise variance 1.5; the fixed-effects fit -1.8, and -12.1 at worst, at rank 3) and its 95 percent interval covered the truth in 95.5 percent of the networks (88 to 100 by design); its sampling standard deviation was 0.76 times the fixed-effects fit's (median) and the root mean squared error of the actor effects 0.79 times (0.18 at the sparsest). The central 50 and 90 percent predictive intervals of 17,150 held-out cells (10 percent of the cells of 70 networks of seven designs) covered 51 to 54 and 90 to 91 percent of the observed counts, with no departure of the randomised probability integral transform from uniformity except at rank 3 with 3.0 nonzero counts per loading (Kolmogorov-Smirnov p = 0.007). The lowest rank within one unit of the largest Laplace value was the generating rank in 142 of 144 networks of 50 actors, and the rank of the largest value in 143 (each rank is compared at the noise variance that maximises its own marginal likelihood).

Counts with little information about the latent product. With fewer than about three nonzero counts per factor loading (fit$eb$variational$information["per_loading"]) the variational posterior shrinks the product. Below 2.1 units of the information-weighted signal per loading (fit$eb$variational$information["signal_per_loading"]) the fit withholds the product: fit$U, fit$V and fit$UVPM are NA, and print() states the statistic and the threshold (see Choosing the rank of a normal fit). The noise variance does not take up what the shrunk product leaves: on 22 networks of 40 and 50 actors at ranks 2 and 3 (1.1 to 6.6 nonzero counts per loading) ve is within 0.41 posterior standard deviations of the sampler's (mean 0.18; the step of the variational iteration itself was above it by 1.3 to 4.6, 2.6 on average, below three nonzero counts per loading), on 60 networks of 50 actors at rank 3 and 2.8 nonzero counts per loading it averaged 0.412 for a truth of 0.40 (0.406 at 4.2), and at a rank above the generating rank it is within 0.6 posterior standard deviations of the sampler's (32 networks; 4 percent below it on average). The variational posterior shrinks a product the data determine a little more than the sampler does: its norm is 0.95 to 0.98 of the norm of the sampler's posterior mean product (lame(), whose prior of the factors is the fit's) from 2.6 nonzero counts per loading and 0.80 to 0.90 up to 2.1 (against ame(), whose prior of the factors is tighter: 0.89 to 0.98 and 0.65 to 0.85); it is as close to the generating product as the sampler's posterior mean is (root mean squared difference 0.99 times the sampler's on average, 0.92 to 1.06 on the 22 networks), the mean of the predictor of a cell differs from the sampler's by 0.15 of the sampler's posterior standard deviation (root mean square over the cells, averaged over the networks; at most 0.23), and its variance is 0.92 of the sampler's (0.74 to 1.03). With about one nonzero count per loading or fewer the iteration still settles: on 54 networks of 40 actors at ranks 1 to 3 with 0.3 to 3.6 nonzero counts per loading every fit converges, none collapses a column of its product, and one, at rank 3 with 0.6 nonzero counts per loading, retires a column (see A rank the data do not support); ve averages 0.39, 0.35 and 0.39 at ranks 1, 2 and 3 for a truth of 0.40. The product correlates with the truth at 0.83 on average at rank 3 and 2.8 nonzero counts per loading (the fixed-effects fit 0.71).

Normal data (the Gaussian empirical-Bayes estimator). A normal fit with the default bias_correct (or "eb") is the same estimator for the model that ame and lame sample, y_{ij} = \mu + \beta' x_{ij} + a_i + b_j + u_i' v_j + e_{ij} with e_{ij} \sim N(0, s^2), and it is the estimator of the binary and count fits with the one cell in which everything is exact. The cell of a normal response has the working response y itself, the weight 1 / s^2 and the score (y - \eta) / s^2; its information does not depend on the predictor, so the Laplace correction of the coefficients vanishes and the Laplace approximation of the coefficients and the effects is not an approximation: at R = 0 the fit is the posterior mean (the best linear unbiased prediction) of the linear Gaussian model at the fitted hyperparameters, and vcov() is its generalised least squares covariance (checked against an independent computation to 1e-6 in the tests). The effects have the Gaussian prior (a_i, b_i) \sim N(0, \Sigma_{ab}) and the inverse-Wishart prior on \Sigma_{ab} that ame() and lame() give a normal response (the scale is the mean variance of the row and the column means of the centred response, \mathrm{round}(4 + 3 n / 100) degrees of freedom), the coefficients are flat, and the node covariates are fixed effects of the prior mean of the effects. The variance of the noise s^2 is a hyperparameter with the sampler's inverse-gamma prior of one degree of freedom on the scale prior$s20, by default the variance of the observed responses. s^2 and \Sigma_{ab} are the fixed points of the EM images that the samplers draw from (the posterior mean of the inverse-gamma and inverse-Wishart full conditionals), formed from the posterior second moments, so that they are net of the degrees of freedom the effects and the factors use: s^2 = (s_{20} + \sum E[e^2]) / (N - 1) with E[e^2] = (y - m)^2 + v (m and v the mean and variance of the predictor of the cell under the posterior). The factors of a latent rank are the variational posterior of the binary fit, with its cavity rule for the scale of their prior and the priors of the samplers (a directed network the shared prior of ame() and lame(), whose scale is that of the effects). The fitted values are the posterior means of the predictor; the response is Gaussian, so they are also the predictive means.

Dyadic correlation of normal data. A directed normal network with dcor = TRUE (the default) takes the two cells of a dyad as one bivariate normal observation with covariance s^2 [[1, \rho], [\rho, 1]]. Everything is in closed form: the information of a pair is the constant matrix [[1, -\rho], [-\rho, 1]] / (s^2 (1 - \rho^2)), whose cross weight enters the equations of the coefficients and of the effects (the estimates are generalised least squares), the rows of the factors see the exact score of the pair, the cavity weights of the factor scale are those of the pair density integrated over the predictors, and \rho takes its posterior mean under the prior of the samplers (density proportional to 1 - \rho^2), computed on a fine grid from the sufficient statistics of the pairs and not by a quadratic approximation. The variance of the noise is estimated from the whitened pairs. The value is fit$rho (fit$VC["rho"]), its posterior standard deviation and method are in fit$rho_estimation, and the record of the pair model in fit$eb$pair. dcor = FALSE fits independent cells and reports the correlation of the residuals of the pairs as before; undirected and bipartite networks have no pairs.

Standard errors of normal fits. vcov() is the Laplace covariance of the coefficients jointly with the effects and the factors, built from the information of the Gaussian cell: the weight 1 / s^2 and the cross weight of the pairs, not the cavity weight that sets the scale of the factors. The information of a Gaussian cell is the same whatever the predictor, so integrating the variance of the predictor into the weight counts it twice: with the cavity weights the standard errors of the dyadic coefficients were 1.05 to 1.13 times the sampler's posterior standard deviations at ranks 1 to 3, with the information of the cell 0.97 to 1.02. Node-covariate and intercept coefficients are covered too, and the part of the uncertainty that comes from the scale of the prior of the factors is added as for the other families. Beyond the designs above the intervals keep their level: on 1,000 simulated normal panels (200 a cell: 40 actors with 8 and 12 periods at ranks 0 and 2, and 120 actors with 4 periods at rank 2; a dyadic covariate, a row and a column covariate, dyadic correlation 0.3) the 95 percent Laplace intervals covered the dyadic coefficient in 0.943 of the panels (cells 0.915 to 0.965), the row and column covariates in 0.959 and 0.956 (0.935 to 0.975) and the intercept in 0.950 (0.935 to 0.955); the Monte Carlo standard error of a pooled figure is 0.007 and of a cell 0.015 to 0.020.

Choosing the rank of a normal fit. fit$eb$log_marginal is the Laplace approximation to the log marginal likelihood of the independent-cell model at the posterior mode, solved from the variational state with the variance of the noise held at its fitted value and to a scaled score of 0.01 (the log marginal likelihood moves by at most m s^2 / 2 from the mode, 0.05 for the 1,000 effects and factor entries of 170 actors at rank 2); it is exact at R = 0, and it is the value of Choosing the rank. The mode is the one fit$irls_diagnostics certifies. fit$eb$variational$information measures the signal of the fitted product in noise variances per loading (for a binary or count fit, signal_per_loading: the product weighted by the information of each cell, which is the same quantity with the fisher weight of the cell in place of 1 / ve). Below 2.1, in every family, the data do not determine a product of that rank: the sampler's posterior mean product is shrunk towards zero and the variational mean is not, so its size is not the posterior's. Such a fit withholds the product: fit$U, fit$V (and fit$L), fit$UVPM and the positions of latent_positions are NA, uv_plot stops with an error of class "lame_product_withheld", and fit$eb$product_withheld records the reason, the statistic, the threshold and the log marginal likelihood of every rank from 0 to the fitted one with the rank it chooses, which print() and summary() state. fitted(), predict() and residuals() are those of the fit (the predictor of the variational posterior), and so are the coefficients, their standard errors and the variance components. The parametric bootstrap draws its replicates from the fit and reports no interval for the positions. The dynamic estimator of a normal panel (lame(method = "als") with dynamic_ab, dynamic_beta or dynamic_uv) records the same statistic in fit$eb$information and withholds in the same way, with the same threshold for static factors and 1.3 per loading of a period for moving positions (the products of every period: fit$UVPM, fit$Omat, fit$Oarr and the positions with their standard deviations); it has no rank evidence, so its record has none.

Scale of a normal response. The iteration works in the units of the response: its tolerance is relative to the standard deviation of the response and the floors of the variances are relative to its variance. The priors of the factors are the samplers' and, like them, written in the unit of the response (R/normal_scale.R): the rows of a directed network have the typical size of the additive effects (the square root of the mean variance of the row and column means of the centred response) as the scale of their prior, the interaction matrix of a bipartite network and the weights of an undirected one have the variance of the response as the variance of theirs. So the fit of c * y is the fit of y in other units: on simulated networks of every kind at ranks 0 and 2 the coefficients, the variance components, the effects and the product of c * y are c (or c^2) times those of y to 1e-13 for c of 1e-3 and 1e3, an undirected panel without factors to 5e-5 (its stop rule), and the standard errors, which carry two refits of the scale of the factors, to 1e-3. A response without variation, or with a standard deviation outside 1e-6 to 1e8, is refused with the remedy (classes "lame_als_data" and "lame_als_scale").

The earlier normal fit. bias_correct = "none" returns the fixed-effects fit of the earlier releases unchanged: the least-squares fit by block coordinate descent (the global solution at R = 0, a local optimum of the non-convex low-rank objective otherwise), with the conditional sandwich covariance of the intercept and the dyadic coefficients, which holds the additive and multiplicative effects fixed.

Validation of normal fits. Against the samplers' posterior means (55 simulated networks of 30 to 50 actors a side, directed, undirected and bipartite, ranks 0 to 3, 38 against ame and 17 against lame fits, 12 of one period and 5 of four periods, of 20,000 scans after 3,000 of burn-in), the 245 coefficients lie within 0.28 posterior standard deviations (mean 0.04), the Laplace standard errors are 0.91 to 1.05 times the posterior standard deviations (mean 1.00; the posterior standard deviation of the intercept of a lame() fit is that of the level of the centred effects, 2.4 to 5.7 times smaller than the standard error of the population intercept, and is left out of this comparison), the variance of the noise is within 0.95 posterior standard deviations, the variances of the sender and receiver effects within 0.16 and the dyadic correlation within 0.45. On 103 simulated networks of 40 actors a side (the 55 above with 48 more of weak and moderate signal; 89 with a latent rank, 1 to 3, directed, undirected and bipartite; against ame and lame fits of 20,000 scans) the norm of the double-centred product is 0.993 to 1.146 of the sampler's posterior mean product where the fit has at least 2.1 noise variances of signal per loading (39 networks; correlation of the products 0.984 to 1.000). Below (50 networks) the undirected fits are 1.37 to 6.9 times the sampler's (0.31 to 1.73 per loading), the directed and bipartite ones 1.08 to 5.0 below 0.7 per loading and 1.01 to 1.25 from 0.7: the sampler's posterior mean is shrunk towards zero and the mean field is not. Every network from 1.99 per loading has a norm within 0.9 to 1.15 of the sampler's, and the largest that does not is at 1.73 (pooled with the binary and count networks below, the threshold is 2.1). Evaluating the update of the scale of the prior at the posterior mean of the inverse-Wishart (options(lame.eb_scale_stat = "mean")) or at the mean-field second moments does not close the gap (median 2.25 and 1.74 times the sampler's below one per loading, against 2.14). The coefficients of the fits below 2.1 lie within 0.30 posterior standard deviations of the sampler's, with standard errors 0.93 to 1.05 times its posterior standard deviations. On 118 binary and count networks of 40 actors a side with a latent product (directed, undirected and bipartite, ranks 1 and 2, against ame fits of 20,000 scans; 71 binary networks of 8 and 25 percent ties and 47 Poisson-lognormal ones) the training events per loading do not measure the information about the product: a bipartite binary fit with 5.6 events per loading had 4.4 to 5.7 times the norm of the sampler's product, and a count fit with 13 nonzero counts per loading 7.5 times it. The information-weighted signal does: below 2.1 (79 networks) the binary products are up to 6.2 times the sampler's (median 1.47; undirected up to 6.2, bipartite 5.7, directed 2.9) and the count products up to 7.5 (median 1.87), and these are withheld; from 2.1 the norm is 0.83 to 1.00 times the sampler's for a binary fit (19 networks: the variational posterior shrinks the product somewhat more than the sampler does) and 0.94 to 1.12 for counts (20 networks). The largest network outside the rule is at 2.03 per loading and the next at 2.17. In repeated sampling (21 designs of 100 to 200 networks, 4,100 networks and 19,700 coefficient fits: every kind, ranks 0 to 3, 20 to 80 actors, one and four periods, 20 percent missing cells, a rank above the generating one and a weak product) the 95 percent Laplace intervals cover 95.2 percent (dyadic coefficients 95.1, node covariates 95.9, intercept 94.3; designs from 94.1 to 96.7). The rank with the largest fit$eb$log_marginal is the generating rank in all of the networks (252 networks of 40 actors, generating ranks 0 to 2) with at least 8 noise variances of signal per loading, in 95 percent from 4 to 8, 89 percent from 2 to 4, 50 percent from 1 to 2, and in none below 1, where it picks a smaller rank; it picked a larger rank than the generating one in one network (a rank-0 network). Of 300 random designs (10 to 60 actors, one and three periods, up to 30 percent missing cells, fitted ranks 0 to 3 against generating ranks 0 to 2), all 300 fit: the two that the iteration did not settle at rank 3 (undirected panels of 14 and 10 actors over generating ranks of 0 and 1, a limit cycle of the scales of the columns) retire two columns (the rank-1 fit) and converge at rank 3 under the priors in the unit of the response.

Speed. The quadratures of the empirical-Bayes fit (the expected working problem of a cell and of a pair of cells, the cavity weights, the bivariate normal probability of the pair likelihood, the log-domain integral for its small values) and the sums over the panel that its equations need run in compiled code; the reference implementations are kept in R and agree with the kernels to rounding (typically 1e-12, and 1e-9 at worst on realistic predictors). The kernels split their cells and pairs over OpenMP threads where the platform has them: options(lame.als_threads = k) sets the number of threads (default 1), and the estimates do not depend on it. Do not combine threads with forked workers such as parallel::mclapply(). The dense solves of the effect step (a Cholesky inverse of the stacked sender and receiver effects, 600 by 600 at 300 actors) are the cost of a large panel and run in the BLAS of the R installation, whose speed sets that cost.

Standard errors. vcov() and confint() return the Laplace approximation to the covariance of the coefficients: the inverse penalised information of all the parameters (coefficients, effects and factors) at the cavity weights (for a normal fit the information of the Gaussian cell, Standard errors of normal fits), with the second moments of the factors in the curvature of the other side, plus the part of the uncertainty that comes from the scale of the prior of the factors. The logarithm of the scale has, to the second order, a normal posterior whose precision is the information of the data about it (the derivative of the log marginal likelihood in the log scale, (\mathrm{tr}(\Sigma^{-1} SE) - n p) / 2 at the cavity second moments, from two refits at the scale times e^{\pm 0.15}) plus that of the hyperprior; the coefficients move with it by J = d\theta / d \log s, and the covariance gains J J' / c. This is one dimension of the uncertainty of the scale (the sampler's covariance of the rows is a full matrix, with 10 free parameters at R = 2); options(lame.eb_scale_uncertainty = FALSE) leaves it out.

A rank the data do not support. A sparse binary or count network at a rank above what it supports (or a network with no signal at all) leaves the iteration at that rank with no stationary point the data determine: a column the data barely inform keeps a live scale, its rows are fitted at the scale of the prior and their second moments are the prior's, so that the posterior mode drifts (the largest score of the stationarity certificate rises with every accepted step) or the variational iteration stops improving at a residual of a few thousandths, or enters a limit cycle in the scales of the columns. Such a fit is not refused. The five-node rule of the pair likelihood of a directed binary network integrates a dyad to less than the gain of a step near the fixed point where the rows are weakly informed, which is a limit of the quadrature as often as of the data: a stall of the pair rule is first refined, to seven and then to nine nodes a side, from the state it reached (fit$eb$pair$refined). An iteration that is still not settled retires its weakest column, the one with the smallest norm of the product \|u_k\| \|v_k\| (the smallest |l_k| of an undirected fit), and the fit of the next rank down starts from the other columns; this repeats while the iteration does not settle. The fit is the one of the rank it reached: fit$R is that rank, fit$R_requested the rank asked for, fit$eb$retirement the record, and the evidence, the standard errors and the factors are those of the effective rank; the fit says so with a message of class "lame_als_retirement" when it is made, and print() and summary() say what was retired. A fit that still does not settle with no column left is refused with an error of class "lame_als_convergence" that names the rank asked for, and so is one that ran out of a max_iter below the default of 200 (a bound the caller chose: the iteration may only need more steps, so it is not reduced in rank). Three sparse binary networks of 50 actors at ranks 2 and 3, YX_bin at rank 4 and a network without a product at rank 2 settle at the rank asked for with the five-node rule, and a count network of 40 actors at rank 3 settles by retiring a column (rank 2). The data determine none of the six products (0.00 to 0.39 units of signal per loading), and every one is withheld.

Predictions. The response of a normal fit is its linear predictor plus Gaussian noise, so fitted() and predict() return the posterior mean of the predictor (fit$fitted), which is the mean of the response; the rest of this paragraph is about binary and count fits. fitted() and predict() of a binary empirical-Bayes fit return predictive probabilities: the link integrated over the posterior of the linear predictor rather than applied at its mean. The predictor of a cell is uncertain, so the plug-in probability of a held-out event can be orders of magnitude too small. A Poisson fit returns the plug-in mean count e^{\eta + v_e / 2} by default, as the samplers do: the posterior predictive mean of a count is not finite in this model (the inverse-gamma tail of v_e, and the inverse-Wishart tail of the scale of the factors for a cell the data do not inform), so the value of the integrated mean e^{m + (v + v_e) / 2} would be a property of the quadrature and of the prior; it stays available as plug_in = FALSE, and the predictive law of a held-out count (its quantiles and event probabilities, which are finite) is what the held-out scores below use. On the 170-state panel of Chapter 11 (10 periods, 15 percent of the cells held out, 39 and 45 events on the two masks) the log loss of the held-out cells at R = 1 to 3 was 0.00427, 0.00440 and 0.00450 (first mask) and 0.00506 at each rank (second mask, where every fit at R = 1 to 3 retires its columns and returns the rank-0 fit) under the plug-in values, and 0.003930, 0.003934 and 0.003972, and 0.004568, under the predictive ones (sampler, one chain, from its draws: 0.003902, 0.003886 and 0.003897, and 0.004480, 0.004448 and 0.004449); the smallest probability of a held-out event was 7e-7 to 3e-6 (first mask) and 7e-9 under the plug-in values and 4.5e-5 to 5.9e-5 and 5e-6 under the predictive ones. At R = 0 the predictive log loss, 0.004054 and 0.004568 on the two masks, is the sampler's (0.004049 and 0.004585). The sampler's own predictor, replaced by a Gaussian with its mean and variance, loses 0.003963, 0.003942 and 0.003947 on the first mask against 0.003902, 0.003886 and 0.003897 from the draws, so at ranks 1 and 2 the gap of the fast fit to the sampler is in the shape of the sampler's posterior of the predictor, not in its moments (at rank 3 the fast fit is 0.00003 above the Gaussian replacement). The predictor is N(m, v): m is the fitted predictor, the posterior mean; v is the variance g' S g of the cell's additive predictor (the coefficients and the two effects, with g its gradient and S the inverse of their penalised information at the expected weights) plus the variance of the product. The probability is \Phi(m / \sqrt{1 + v}) for a probit link and \mathrm{plogis}(m / \sqrt{1 + \pi v / 8}) for a logit link. plug_in = TRUE returns the link applied to the fitted predictor (what fit$fitted holds and simulate() uses). Forming the covariance costs one dense inverse, of the size of the one behind vcov(), per call.

Choosing the rank. fit$eb$log_marginal is the Laplace approximation to the log marginal likelihood of the data with the effects and factors integrated, at the posterior mode that the variational stage starts from (for a normal fit the mode at the fitted variance of the noise, Choosing the rank of a normal fit), with that mode's coefficients and scales: the deviance and the penalties of the effects and the factors, with the log determinant of the information of all of them, on the scale of the likelihood (an undirected network counts each dyad once; the hyperpriors of the scales are left out). The deviance falls with the rank; this does not, so the difference between two ranks fitted to the same data is an approximate log Bayes factor. The sign and label multiplicity of the factors is not counted, which favours the lower rank. On 1,040 simulated binary networks (directed, bipartite and undirected; generating rank 0 to 3, fitted ranks 0 to 3; 60 and 100 actors, one and four periods, 7 to 25 percent ties and two sparser designs) 17 of the 4,160 fits were refused as not converged (class "lame_als_convergence"), all in the two sparsest designs without a latent term; of the 1,033 networks fitted at all four ranks, the rank with the largest value was the generating rank in all of the 309 with at least 10 training events per factor loading, in 96 percent of the 386 with 3 to 10, in 59 percent of the 234 with 1 to 3 (35 percent were too low) and in 38 percent of the 104 with fewer than 1 (27 percent were too high, 36 percent too low); the lowest rank within one unit of the largest value was right in 100, 96, 62 and 43 percent. With little information per loading the data do not determine the rank. Take the lowest rank whose value is within one unit of the largest: it named the generating rank more often than the largest value below 3 events per loading (62 against 59 and 43 against 38 percent) and as often from 3 to 10 (96 percent). The variational lower bound, fit$eb$variational$elbo, is not the value to use for this. A Gaussian posterior cannot hold the broad posterior of positions that the data barely identify, and the bound's slack grows with the rank: on the 1,033 networks its largest value was the generating rank in 35, 56, 96, 100 and 100 percent (fewer than 1, 1 to 3, 3 to 10, 10 to 30 and more than 30 events per loading), never too high but too low for 96 percent of the networks with a latent term and fewer than 1 event per loading, where the Laplace value found the rank in 38 percent (and was too high for 64 percent of the networks without a latent term and fewer than 1 event per loading). On the 170-state panel of Chapter 11 the Laplace value was -1137.8, -1119.3, -1117.2 and -1121.1 at R = 0 to 3 on the first mask, so ranks 1 to 3 are preferred to 0 (by 17 to 21 units, most at rank 2), as the sampler's held-out log loss prefers them, while the lower bound falls with the rank (-1144.7, -1152.9, -1170.6 and -1195.9). On the second mask every fit at R = 1 to 3 retires its columns and returns the rank-0 fit (Laplace value -1117.3, lower bound -1125.1), so neither value can prefer a latent term there, although the sampler's held-out log loss does. The value is the evidence of the independent-cell model at the posterior mode: a fit with a dyadic correlation (Dyadic correlation) reports the same log_marginal as the fit with dcor = FALSE. On 15 simulated networks of 60 actors with a latent correlation of 0.6 (generating ranks 0 to 2) the rule above picked the generating rank in 10 and a higher one in none.

What picks the rank where the Laplace value cannot is the held-out predictive log score: mask a share of the observed cells (15 percent), refit at each rank, and score the masked cells with fitted() (the integrated predictive; a variational fit's predictive is not shrunk to the plug-in probability, and its surplus columns are not shrunk away, so the score falls above the generating rank). On the 1,033 networks it found the generating rank in 51, 72, 85, 91 and 92 percent (the Laplace value in 38, 59, 96, 100 and 100), and where there are fewer than 3 events per loading it was right for 64 percent of the networks without a latent term (the Laplace value 72, with the one-unit rule 87; the lower bound 100) and 66 percent of those with one (Laplace 45, the lower bound 32): no other criterion that needs no refit was right for both. The rule that uses the held-out score below 3 events per loading (fit$eb$variational$information["per_loading"] of the fit at the rank the one-unit rule names, rank 1 where it names 0) and the Laplace value with the one-unit rule above found the generating rank in 85 percent of the 1,033 networks (the Laplace value alone 83 and 84, the lower bound 82, the held-out score alone 81). On the 170-state panel the held-out log loss at R = 0 to 3 was 0.004054, 0.003930, 0.003934 and 0.003972 on the first mask and 0.004568 at every rank on the second, whose fits at R = 1 to 3 retire to rank 0 (the sampler's, one chain: 0.004049, 0.003902, 0.003886 and 0.003897, and 0.004585, 0.004480, 0.004448 and 0.004449; the sampler's loss is lowest at rank 2 on both masks, with rank 3 within 0.00001 of it).

The variational stage needs a median of 46 iterations (87 with fewer than 1 event per loading), and 9 in 10 fits finish within 130. The 170-state panel of Chapter 11 takes 29 seconds at R = 0 and 1.5 to 3 minutes at R = 1 to 3 (one start, one thread; 8 minutes at R = 2 with multistart = "full"). The joint solve is dense, so the cost grows with the cube of the number of actors. An actor without a tie (or without a non-tie) on a side keeps a finite, shrunken effect and stays in the fit, so there is no fit$not_estimated. The effects a and b are not re-centred: the intercept is the population intercept (the level of a dyad between two actors with zero effects, at the means of the node covariates) and the effects have mean zero a priori; without node covariates the posterior-mean effects of the variational fit sum to zero, because the intercept is stationary. va, cab and vb are the fitted covariance of the effects, the EM estimate of \Sigma_{ab}, not descriptive summaries of the estimated effects.

Accuracy. Against the truth, on the 1,040 simulated networks above (dyadic slope 0.7; sender and receiver sd 0.6, correlation 0.3; variance 0.6 a latent dimension), fitted at the generating rank (above zero), the relative bias of the slope was, by training events per loading (fewer than 1, 1 to 3, 3 to 10, 10 to 30, more than 30), -3.9, -2.4, -0.8, -0.1 and 0.4 percent, with root mean square error 10.4, 7.1, 4.8, 2.8 and 1.8 percent (the posterior mode: -0.7, 0.8, 1.1, 0.7 and 0.8; 9.8, 7.2, 4.9, 2.9 and 1.9). One rank above the generating one the bias was -0.5, 1.4, 2.2, 1.3 and 1.0 percent (mode 0.5, 5.1, 5.2, 2.6 and 1.8), and with no latent term in the data and a latent rank fitted, with fewer than 3, 3 to 10, 10 to 30 and more than 30 events per actor, it was 17.5, 5.7, 4.0 and 1.1 percent (mode 12.4, 8.5, 6.7 and 2.1). At the generating rank the fit withheld the product in 85, 70 and 16 percent of the fits with fewer than 1, 1 to 3 and 3 to 10 events per loading and in none above; the products it kept correlated 0.57, 0.78, 0.87, 0.93 and 0.96 with the generating one and had norm 43, 54, 57, 58 and 52 against the truth's 81, 72, 67, 63 and 54 (the mode, every fit: correlation 0.48, 0.68, 0.85, 0.93 and 0.96, norm 49, 60, 61, 62 and 54 against 69, 62, 63, 63 and 54), and the held-out log loss was no higher than the mode's in any band. Against the sampler's posterior mean, on 101 fits with fewer than 1 training event per loading and 96 with 1 to 3 (directed, bipartite and undirected networks of 60 and 100 actors, ranks 1 to 3, one lame() run of two chains each), the slope of the fit lay on average 0.09 (95 percent interval -0.001 to 0.18) and 0.10 (0.05 to 0.16) sampler standard deviations from it (relative to the slope 1.1 and 1.2 percent), and the standard error was 1.05 and 1.03 times the sampler's standard deviation. The fit withheld the product in 90 and 88 of these fits; of the 135 directed fits, 132 withheld it, two with fewer than 1 event per loading retired every column and one kept it. The 9 and 8 products it kept (bipartite and undirected, and the one directed) had 0.74 to 1.08 and 0.92 to 1.07 times the norm of the sampler's posterior mean product. For undirected networks the slope lay 0.40 and 0.34 standard deviations above the sampler's. Where the sampler finds no structure the variational mean keeps a sign pattern from the mode: in 4 of 24 networks of 60 actors its norm is 2.3 to 7.8 times the sampler's posterior-mean norm, which is the small remainder (0.11 to 0.37) of a cloud of mutually near-orthogonal, sign-symmetric draws of about the size of the fitted product (mean cosine between draws 0.004 to 0.11); a unimodal Gaussian posterior cannot report that mean, and the fit withholds these four products (0.30 to 1.18 units of signal per loading). The slope on those four (|z| \le 0.83), its standard error (1.08 to 1.10 times the sampler's standard deviation) and the held-out log loss (0.011 below to 0.002 above the sampler's) are within the usual bands. For bipartite networks of 60 actors fitted at R = 3 (5 networks, 0.65 to 1.22 events per loading) the slope lies within 1.52 sampler standard deviations of the sampler's (|z| \le 0.41 on four) and the fit withholds the product, whose variational mean has 0.76 to 1.21 times the norm of the sampler's posterior mean product. The scales of the sender columns are the eigenvalues of a Wishart matrix, as in the sampler, so that a column the data do not inform falls to the small eigenvalues: on these 5 networks the variance of the predictor of the held-out cells is 0.57 to 1.02 times the sampler's, and the held-out log loss 0.003 below to 0.014 above its. In cells of 40 networks (directed 30 to 100 actors, bipartite and undirected 60 actors, one and four periods, ranks 0 to 2, 14 cells) the slope bias was -2.2 to 1.2 percent (mode -0.5 to 3.4), the 95 percent interval from vcov() covered the slope in 80 to 100 percent of the networks (mean 95) and the intercept in 88 to 98 (mean 92), and the standard error was 0.78 to 1.43 times the standard deviation of the slope. On the 170-state panel the coefficients of the full-system fit (R = 2) lay within 0.31 standard deviations of the sampler's posterior mean, with standard errors 0.86 to 1.23 times the sampler's; the fit withholds the product (0.26 units of signal per loading, below 2.1), whose variational mean has norm 14.4 against the 15.9 of the sampler's posterior mean product (correlation 0.93). On the collinear-covariate design (40 actors, 8 periods, R = 2, truth 0.5, 24 networks) the covariate that follows the latent product was 28.0 percent too large (sampler 26.2, mode 13.5, no latent rank 93): the prior shrinks the product and the covariate absorbs what it shrinks. The independent-cell fit (dcor = FALSE) has no dyadic correlation. On a strongly reciprocal panel (the school network of Chapter 9, 32 actors, R = 2; ame() with the book's specification, two chains of 50,000 scans, puts rho at 0.67, posterior standard deviation 0.16) its coefficients lie within 0.72 sampler standard deviations of the sampler's, with standard errors 1.07 to 1.14 times the sampler's, and its post hoc correlation is 0.98, 2.0 posterior standard deviations above the sampler's. The pair model (Dyadic correlation), the default for a directed network, estimates rho at 0.74 (0.43 posterior standard deviations from the sampler's), and its coefficients lie within 0.74 sampler standard deviations, with standard errors 0.76 to 0.97 times the sampler's. Both fits withhold the product of this panel. Choose the rank as described above.

Without the empirical-Bayes prior (bias_correct = "none" and the two explicit corrections below) a binary fit is the maximum-likelihood fixed-effects probit (or logit) fit, penalised at R > 0: each reweighting step caps the working-response increment at 30 and gives a capped cell the weight at which it still contributes its score, and the loop runs to its fixed point (on 32 simulated networks of 60 and 120 actors the coefficient and deviance matched glm() to 5e-6 and 3e-8). In a fixed-effects fit an actor with a single outcome value on a side (a sender without a tie, a receiver with only ties, and so on) has its effect there at the boundary of the likelihood and its cells on that side carry no information about the other parameters: those cells are left out of the fit (and of fit$Y), the effect is NA in a / b with its linear predictor and fitted values, the actors are listed in fit$not_estimated with a message, and va / vb describe the estimated effects. On simulated 250-actor networks with 1 to 50 such senders the coefficient, va and vb stayed where they were without them.

Split-sample bias correction. A binary IRLS fit is a fixed-effects probit (or logit) fit: every additive effect and every latent position is estimated from the cells of its own actor, and maximising over all of them stretches the latent scale outward in finite samples (the incidental-parameters problem), while the factor penalty pulls the product term in. On simulated 30 x 20 bipartite panels with five periods and R = 2 (20 seeds), the uncorrected dyadic coefficient was 3.3\ 15\ too flat; at 60 x 40 the stretch was 2 to 5\ panels 2 to 6\ split-sample jackknife wherever the jackknife holds (see the size boundary below; the empirical-Bayes default needs no such correction). The observed cells are split into two interleaved halves, the model is refit on each half, warm-started from the full fit with the same settings, and every identified quantity is reported as 2 \times full minus the mean of the two halves: a half carries twice the stretch of the full panel, so the combination cancels it to first order. A half that has not converged within max_iter is continued from its last iterate with five times that limit; if it still has not converged the call is an error of class "lame_jk_nonconvergence" (the class of the dynamic count jackknife), since the correction never combines a last iterate. The split is the checkerboard of the cells with i + j + t even and odd (every dyad keeps the periods of one parity in each half, every actor about half of its cells); when it leaves some actor without a tie or without a non-tie in a half, a random split of the dyad-periods (drawn from a fixed internal seed; the caller's random number stream is left unchanged) and then the split into odd and even periods are tried, the latter with weights for its unequal halves. Corrected are the intercept and regression coefficients, the additive effects a and b, the double-centred product surface (U, V, L are the rank-R factorisation of the corrected surface) and va, vb, cab; EZ, the fitted values and the residuals are rebuilt from them; ve and the IRLS weights are those of the full fit, and rho is re-estimated from the paired binary likelihood at the corrected marginal probabilities. On the panels above the corrected coefficient ratio was 0.990 (Monte Carlo standard error 0.010) at 30 x 20 and 1.002 at 60 x 40, the effect spreads within 1.5\ 2\ product slope 0.97 to 1.04. The correction costs two half-data refits, 1.6 to 3.4 times the time of the uncorrected fit, and widens the spread across samples slightly (standard deviation of the coefficient 0.043 to 0.047 at 30 x 20). fit$raw keeps the uncorrected coefficients, a, b, VC, product surface O and deviance (fit$deviance is that of the corrected EZ); fit$bias_correct records whether the correction was applied, the split used and the splits tried, the half fits' convergence flags and iteration counts, and the shrink factors of the effects and the product. bias_correct = "none" returns the uncorrected fit and bias_correct = "jackknife" requests the correction for another family.

A single cross-section at R = 0 (ame_als()) is split by the checkerboard of its cells (i + j even and odd, each actor keeping half of its cells in each half). Each half's observed graph then falls into two sender-receiver components whose relative level is not identified, so every cell's half value is taken from the half that observes that cell's block. On simulated directed networks (20 to 40 seeds) the corrected slope, va and vb were 0.990, 0.982 and 0.954 of the truth at 75 actors with 36\ 100 and 0.999, 0.979 and 0.996 at 150 actors with 24\ 1.017 to 1.029 and up to 1.09). At R > 0 a single cross-section is not corrected: the checkerboard over-corrected wherever the uncorrected estimate needs it (100 actors, R = 2: slope 1.024 to 0.971, vb 1.090 to 0.950; 50 actors: va 1.209 to 0.791), and so did a random split. Such a cross-section is corrected by the parametric bootstrap instead (see Parametric-bootstrap bias correction).

Parametric-bootstrap bias correction. bias_correct = "bootstrap" simulates the bias instead of splitting the sample, which reaches the panels the split-sample correction does not (a single cross-section with R > 0 or below 25 rarer-outcome cells per actor, a panel of fewer than five periods at R > 0, an actor without both outcomes in a half). Networks are drawn from the fitted model with fresh additive effects from the fitted va, vb, cab and the intercept, coefficients, node-covariate parts, product surface and reciprocity held at their fitted values, and are fitted by the same uncorrected estimator, warm-started from the fit. The bias of each coefficient, of va, vb, cab and of each cell of the product surface is the mean of its replicate values less the value it was generated with (the variance and covariance of the drawn effects for va, vb, cab); the additive effects are divided by the mean slope of the replicate effects on the drawn ones. The bias grows with the spread of the effects, so it is larger at the stretched estimate than at the truth, and one subtraction over-corrects small networks (30 actors, R = 0, 150 networks: slope -2.4\ subtracting the bias at the fitted effects themselves gave 0.85); a second round of replicates at the once-corrected parameters, drawn with the same random numbers, gives the slope s of the bias, and the estimate is corrected by b_0 / (1 + s), the secant step toward the parameters whose simulated estimate equals the fitted one. Each round has 30 replicates from a fixed internal seed (the fit is reproducible and the caller's random number stream is left unchanged), so the correction costs about 60 refits. On simulated directed probit networks (mean -1, about 27\ ties, dyadic slope 0.7, sender and receiver sd 0.6, factor variance 0.6 at R = 1; 100 to 150 networks a cell) the relative bias of the slope went from 12.2\ at 30 and R = 1, 6.8\ 0.2\ 1.9\ panels of four periods at R = 1, and on 40 x 30 bipartite cross-sections (60 networks) 9.3\ 2.4\ va and vb went from 1.10 to 1.51 of the truth to 0.94 to 1.08 (Monte Carlo standard error 0.02 to 0.05), the intercept bias from -0.03 to -0.11 to within 0.02, and the slope of the product surface on the generating one from 0.71 to 0.79 (30 actors), 0.86 to 0.93 (50) and 0.94 to 0.98 (100); the rest of its flattening is the factor penalty. The 95\ interval estimate plus or minus 1.96 replicate standard deviations of the first round (fit$bias_correct$replicate_sd) covered the slope in 93 to 99\ bootstrap-corrected fits among 100-actor cross-sections at R = 0), against 72 to 92\ uncorrected estimate and 49 to 82\ (vcov()). The correction took 4 to 8 seconds at R = 0 and 40 seconds (30 actors) to 9 minutes (200 actors) at R = 1 with twelve such fits running at once, against 0.2 to 9 seconds uncorrected. fit$bias_correct records method = "bootstrap", the replicates of each round that refit and converged, the first-round biases, the secant factors and the corrections applied; fit$raw keeps the uncorrected values. A refit of the corrected fit (ame_als_refit, and every replicate of ame_als_bootstrap) removes these stored corrections rather than simulating again. Fewer than 80\ and converging is an error of class "lame_als_size" (1 of 100 networks at 30 actors and R = 1), and so is an undirected fit with R > 0, whose fitted product is flattened beyond what the replicates reproduce: on undirected 50-actor networks with a rank-1 term of variance 1 (60 networks) the product slope was 0.68 and the correction moved the dyadic slope from 1.0\ variance 0.36 (50 networks) the factor penalty removed the product in every fit and the slope went from 6.5\ Undirected fits at R = 0 are corrected (50 actors: slope 6.6\ -0.8\

Size boundary. A fixed-effects binary fit is supported where the bias of what it reports is within a tolerance, 2\ va and vb (the internal constant .ae_bias_tol; the floors below follow from it). The size that governs the bias is the number of observed cells per actor (over all periods, on the side with fewer: senders or receivers, rows or columns) times the share of the rarer outcome (ties in a sparse network), called rarer-outcome cells per actor below, and the bias also grows with the spread of the additive effects, measured on the fit as the larger of the sd of the estimated sender and receiver effects over the latent noise sd (1 for probit, \pi / \sqrt{3} for logit). On simulated networks (single cross-sections of 60 to 500 actors at 7 to 40\ four periods, directed, bipartite and undirected, probit and logit, additive-effect sd 0.5 to 2, 3 to 10 seeds a cell) the uncorrected bias fell as c / size, with c growing with the rank and the spread. The split-sample correction is used where it can be computed: panels with two or more periods at R = 0 and five or more at R > 0, and single cross-sections at R = 0, in each case with every actor keeping a tie and a non-tie in both halves and, at R = 0, from 25 rarer-outcome cells per actor (below that it over-corrected va and vb to 0.75 and 0.87 of the truth at 14, and to 0.86 and 0.89 at 18). The empirical-Bayes default does not use this rule, which governs the fixed-effects fits of the families and settings the estimator does not cover: a fit the split-sample correction does not reach is reported uncorrected when its size is at or above the floor for its rank and spread, and refused below it (the parametric bootstrap does not reach an undirected fit with a multiplicative term); a fit whose spread is above 2.25, or whose rank is above 8, is refused, corrected or not, since neither the uncorrected nor the corrected estimate was measured there. The floors (rarer-outcome cells per actor; in brackets those at 3\ tolerance up to one Monte Carlo standard error:

"Refused" marks a band where no measured size (up to about 175 cells per actor) met the tolerance; ranks above 8 are not measured and refused. For a directed cross-section with 36\ are about 204, 218 and 326 actors at R = 0, 1-2 and 3-8, and half that with two periods. A fit above the floor records in fit$bias_correct method = "none" with the size, the floor, the spread and a basis giving the measured bias band there, and print() / summary() say that no correction is needed; there is no warning. Above the measured spread or rank, binary ALS stops with an error of class "lame_als_size" that names the spread or rank and the measured band; the MCMC ame() / lame() integrate the effects and positions out and are the calibrated route there. The parametric bootstrap also stops with that class when fewer than 80\ replicates refit and converge (networks simulated from the fit too sparse for the fixed-effects estimator). bias_correct = "none" returns the uncorrected fit at any size, recorded as requested (fit$bias_correct$requested); bias_correct = "jackknife" requires the correction and is an error (class "lame_jk_split") where it cannot be computed. A bootstrap of a binary fit that carries no correction is an error (class "lame_als_bootstrap_bias"), whatever its size: each replicate is refit at the biased estimate and reproduces the bias, which falls more slowly than the standard error (ame_als_bootstrap). The default empirical-Bayes fit and the fits corrected by the jackknife or the parametric bootstrap have one.

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.

Missing covariates. As in ame and lame, a cell whose dyadic covariate, sender covariate or receiver covariate is missing in its period is unobserved: it is set to NA in Y before the fit, a message gives the number of cells, and the coefficients are estimated from the complete cells. A missing value is not modelled as 0, which biases the coefficient of a covariate that is not centred. A missing value on the diagonal of a square unipartite covariate is not counted, since the diagonal is unobserved anyway. In an undirected fit (symmetric = TRUE) a dyad has one covariate value: a hole in one triangle of a dyadic covariate is filled from the other triangle, and a dyad missing in both is unobserved. The per-actor mean of a node covariate runs over the periods in which it is observed. An actor missing a node covariate in every period has no observed cell on that side: the empirical-Bayes fit leaves its effect there to the prior (zero when the actor has no cell on either side), and a fixed-effects fit reports it as not estimated (NA in a / b, listed in fit$not_estimated as "no observed cell"). A covariate that is missing everywhere stops the fit, as it does in ame and lame.

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. An empirical-Bayes fit (the default for normal, binary and Poisson responses) has a rank criterion, fit$eb$log_marginal (see Choosing the rank). For the fixed-effects fits 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 for the normal family and the transform route, the binomial or Poisson deviance of the fitted values for an IRLS fit, in each case of the reported EZ; the optimiser's final objective, which adds the factor penalty of an IRLS fit, is objective): 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.

A completely separating unpenalized linear predictor has no finite estimate; the IRLS route refuses it with class lame_als_separation. Use ame() or lame() with proper priors for such data.

Uncertainty. An empirical-Bayes fit (the default for normal, binary and Poisson responses) reports Laplace standard errors for every coefficient (vcov(), confint(), fit$eb$se). For a fixed-effects fit (bias_correct = "none") 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. When R > 0 the fitted multiplicative surface is exposed as UVPM (U L U'} for symmetric fits, which also carry it as \code{ULUPM}; \code{U V' otherwise), matching the slot names on MCMC fits. 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. R is the rank of the fit; where the data did not support the rank asked for and columns were retired, R_requested (present only then) is the rank that was asked for (see A rank the data do not support). For IRLS, irls_iterations counts accepted outer reweighting steps; iterations counts the final inner solve (the reweighting steps, for an empirical-Bayes fit). irls_diagnostics records the score and variance checks, their threshold, the stopping reason, and each accepted step's improvement under a common factor penalty. These checks concern the uncorrected fit; a split-sample correction also records each half's checks in bias_correct$irls_diagnostics. A changing empirical factor variance means successive penalized objective values need not decrease; an empirical-Bayes step is also accepted when it lowers the residual of the estimating equations that the checks evaluate. A fixed-effects fit that does not pass the checks within max_iter, nor within four times it, is an error of class "lame_als_convergence": no fit is returned. (An empirical-Bayes fit that does not pass them warns with that class and has converged = FALSE.) For an empirical-Bayes fit (bias_correct = "eb") eb is a list: sab, the fitted covariance of the actor effects (a 2 \times 2 matrix; an undirected fit has one variance and repeats it on the diagonal), Sab0 and eta0, the scale and degrees of freedom of its inverse-Wishart prior, gamma, the coefficients of the node covariates in the prior mean of the effects, sab_change, the relative change of sab in the last update, sig2_u and sig2_v, the fitted scales of the factor columns (one value per column and side; the isotropic scale of a directed fit repeated, the unit receiver scale of a bipartite fit, one scale per column of an undirected fit), mu_int, a_raw and b_raw, the intercept and the effects as the posterior is built from them (the effects are not re-centred), log_marginal, the Laplace log marginal likelihood at the posterior mode (see Choosing the rank in Details), vcov and vcov_cond, the Laplace covariance of the coefficients (intercept, dyadic and node-covariate) jointly with the effects and the latent factors and conditional on the factors, se, the square roots of the diagonal of vcov, se_method and vcov_error, the reason when the covariance could not be computed (vcov() then falls back to the sandwich), estimator, which estimator produced the reported values ("variational", the default, or "laplace", the posterior mode alone), variational, the record of the variational iteration (converged, termination, iterations, residual, the largest step a full update would still take, and gain, the objective's gain at the last step; step, the step length; tolerance; elbo, the evidence lower bound at the plug-in scales, with ell, kl_actors and kl_factors; information, the training events, the loadings and the events per loading of a binary or count fit with a latent rank, with signal, the double-centred fitted product squared and weighted by the fisher weight of each observed cell, signal_per_loading and threshold, 2.1, below which the product is withheld (a normal fit: signal, the squared norm of the double-centred fitted product summed over the periods, loadings, per_loading, the signal in noise variances per loading, and threshold); scale, the rule that set the scale of the prior of the factors ("cavity" or "plugin"), the scale s (the mean over the columns of a bipartite or undirected fit), the correlation rho of the sender and receiver position of an actor (directed) and sd_log, the standard deviation of the log scale that vcov() carries; start, mode_log_marginal), factor_posterior, the covariances of the posterior of the factor rows (arrays Sig and Om, or Sig and the factor F with its signs sg for an undirected fit; the cavity covariances, with the covariances of the fixed point of the updates beside them as Sig_mf and Om_mf) and the prior of the fit, suv, the fitted 2R \times 2R covariance of the factor rows of a directed network, laplace, the parts of log_marginal, and, for a directed binary probit fit with a dyadic correlation, pair: the correlation rho, its posterior standard deviation sd, the score and curvature of the expected log likelihood in it, the number of dyads with both cells observed n_pairs, the nodes a side of the quadrature quad (zero for a normal pair, which has none), the rules the quadrature was refined to after a stall refined (empty where it was not) and the expected cross weights of the pairs; retirement, present only when columns were retired: a list with the requested_rank, the effective_rank, whether the fit of the effective rank settled, and steps, one record per retired column (the rank it was retired at, its column, its strength and share of the strongest column, and the stage and termination of the iteration that did not settle); for a normal fit, ve, the variance of the noise (also VC["ve"] and fit$s2) and ve_prior, the degrees of freedom and the scale of its prior; and, for a Poisson fit, ve, the variance of the lognormal noise of the count (also VC["ve"], fit$s2 and fit$pois_latent_var), and ve_search, the search that found it: ve, objective (the Laplace log marginal likelihood plus the log prior density of ve at the maximum), curve (a data frame of the values tried: ve, objective, whether the mode there solved its equations and its iterations) and evaluations. A fit with multistart = "cheap" or "full" (and R > 0) also records, for each start, multistart_converged, multistart_criterion and multistart_sse (the value of that criterion), the kept start multistart_selected and, for IRLS fits, multistart_deviance and multistart_distance (the distance of the fitted predictor from the kept start's, in deviance units); and for every multistart fit multistart_coefficients, a starts by coefficients matrix of what each start reports (intercept, dyadic and node-covariate coefficients, with the names and on the scale of coef(fit)) whose row for the kept start is coef(fit) (the uncorrected fit$raw$coefficients of a bias-corrected fit, since the correction is made on the kept start only), and multistart_product_distance, the Frobenius norm of the difference between the UVPM each start reports (before any bias correction) and the kept start's, over the norm of the kept start's: zero for the kept start. A kept UVPM of no size, when the factors of the starts shrink to nothing, is measured against sqrt(.Machine$double.eps), so that starts that agree there are at distance zero. The surfaces themselves are not stored. The variance-component vector has five descriptive entries:

va, vb

variances of the sender and receiver additive effects. The empirical variance of the estimated effects carries their estimation error as well, so the first-order sampling variance of the additive block — ve times the centred trace of the inverse two-way weighted normal equations — is subtracted, floored at zero, and va/vb estimate the variance of the effects themselves, the quantity ame / lame report. The uncorrected empirical variances are kept as vc_raw and the correction as vc_noise; it is computed conditional on the fitted multiplicative term. For an empirical-Bayes fit they are the fitted variances of the effects themselves (fit$eb$sab).

cab

covariance of the estimated sender and receiver effects, uncorrected (NA for symmetric or bipartite models); for an empirical-Bayes fit, their fitted covariance.

rho

dyadic reciprocity, distinct from the sender/receiver correlation cor(a, b). For a directed binary probit fit with the empirical-Bayes estimator and dcor = TRUE it is the correlation of the fitted model (Dyadic correlation), and rho_estimation holds its posterior standard deviation and method = "variational pair likelihood". For other directed binary IRLS fits, this is the Gaussian copula correlation estimated from the paired binary likelihood at the reported marginal probabilities. With a probit link it is the latent error correlation. The working-residual correlation is retained as rho_working; rho_estimation records the number of reciprocal pairs, likelihood, and boundary or unidentified state. Without informative reciprocal pairs rho is NA and conditional network simulation is refused. Poisson IRLS fits estimate latent correlation from the weighted count covariance identity \mathrm{cov}(Y_{ij}/\mu_{ij}-1,Y_{ji}/\mu_{ji}-1) = e^{\rho v_e}-1. Independent Poisson noise contributes to the marginal variances but not this covariance. (This is the estimate of a fit without the pair stage, bias_correct = "none" or dcor = FALSE; a directed empirical-Bayes count fit estimates the correlation from the pairs, see Dyadic correlation of counts.) With zero latent variance, correlation is unidentified and independent count simulation remains valid. Normal and transformed-response fits use reciprocal residual correlation.

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. For a poisson IRLS fit ve is the variance of the lognormal latent layer of the model ame / lame fit (y \sim \mathrm{Pois}(e^{z}), z \sim N(\eta, v_e)); an empirical-Bayes fit estimates it by maximising the Laplace marginal likelihood (see Counts), and the fixed-effects fit (bias_correct = "none") from the count overdispersion \mathrm{var}[y] = \mathrm{E}[y] + (e^{v_e} - 1)\mathrm{E}[y]^2 and 0 for equidispersed counts; the working-response mean square is again ve_working.

Poisson level scale. The fixed-effects IRLS fit (bias_correct = "none") solves the marginal log-mean \mathrm{E}[y] = e^{\eta}, which is the latent predictor plus v_e/2. The reported intercept and EZ subtract that Jensen term, so they are on the same latent scale as ame / lame; the shift is stored as pois_latent_var and added back by fitted, predict(type = "response") and simulate, which therefore do not change. The dyadic slopes are unaffected: the term is constant across dyads. The empirical-Bayes fit works on the latent scale throughout and needs no shift.

These are summaries of the point estimates, not draws from a posterior. 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 estimator for normal, binary, and Poisson panels (an empirical-Bayes fit by default; see bias_correct in lame), including named panels where actors enter or exit: dynamic_ab, selected intercept/dyadic/node dynamic_beta, AR(1) dynamic_uv for directed, symmetric, and bipartite panels, and Student-t dynamic_uv for directed and symmetric panels (bipartite t requests are refused). dynamic_uv_kind = "snap" has no ALS route and is refused. 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. A time-varying bipartite interaction matrix (dynamic_G = TRUE) is refused on every route: use dynamic_uv = TRUE and read fit$UVPM. A dynamic normal fit whose latent product the data do not determine withholds it (fit$eb$product_withheld; see Choosing the rank of a normal fit in ame_als). Changing actor composition requires row and column names on every slice so actors can be aligned; smoothing penalties are broken across actor-entry gaps (the empirical-Bayes fit treats an absent actor-period as a free AR state). 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. ame_als(), lame_als() and lame(..., method = "als") read one entry of prior, s20, the scale of the prior of the noise variance of a normal or Poisson fit, and the dynamic empirical-Bayes fits read the hyperprior entries of lame in their EM. Every other entry of the samplers' prior list and g = ... are MCMC-only; the dispatcher ame(..., method = "als") warns and ignores prior and g (pass prior = list(s20 = ) to ame_als()). The empirical-Bayes fits apply the default additive-effect prior of ame / lame with its covariance estimated.

• 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: an ALS fit is one point fit, chains has no effect, and n_chains is ignored with a warning.

• Convergence diagnostics. No Rhat / ESS / trace_plot.

Features that work the same on ALS fits: coef, vcov (Laplace covariance of an empirical-Bayes fit, sandwich on the regression block of a fixed-effects fit), confint (auto-routes bootstrap intervals when present, Wald intervals 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, and only when the panel supplies at least ceiling(T / block_length) >= 10 resampling blocks: below that the ends of the interval are the extreme replicates (a plain percentile interval covered the dyadic coefficient in 0.83 of the panels at five blocks), so the call warns and uses parametric instead. The returned object records the scheme actually used in type and the one asked for in requested_type. The interval end points are read at expanded percentile levels, and the scheme returns no standard error or interval for the parameters that vary over actors and not over time; see Block bootstrap below.

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. An empirical-Bayes fit (fit$bias_correct$method == "eb", the default for normal, binary and Poisson fits; see ame_als) has every replicate refit with the same estimator, re-estimating the covariance of the actor effects from the fitted one, and draws the replicate effects uncentred, since the fitted intercept is the population intercept; the intervals then describe the empirical-Bayes estimator, and agree with the Laplace standard errors of vcov.ame_als where those are calibrated. A fit carrying the split-sample jackknife (fit$bias_correct$method == "jackknife") has every replicate corrected the same way, so the intervals describe the corrected estimator; a replicate whose split cannot be refit is dropped and counted like any failed replicate. A fit carrying the parametric-bootstrap correction (fit$bias_correct$method == "bootstrap") has the bias it estimated removed from every replicate rather than re-estimated, so a replicate costs one refit; the bias moves slowly with the parameters, and the replicate spread is that of the corrected estimator to first order. A binary replicate in which an actor has a single outcome value on a side records that actor's effect there as -Inf (no tie) or Inf (only ties), its maximum-likelihood value, so the actor's interval is open on that side and its standard error is Inf; an actor the original fit leaves out is NA throughout. A dynamic ALS fit is refit dynamically, on the parametric route only: a block resample reorders the time axis the temporal penalties are defined along. Only a penalised normal-family dynamic fit (bias_correct = "none") has a bootstrap, and its time-varying factors are not summarised. Every replicate redraws the additive paths (an actor level with the fitted level variance plus an AR(1) deviation), so the standard errors carry the variability of the actors: the intercept's interval is for the population intercept (it covered 0.90 and 0.82 of 100 panels of 20 actors and 8 periods and 40 of 30 and 12, against 0.08 and 0.07 when the fitted paths were held fixed, which covers the realized intercept only) and the dyadic coefficient's covered 0.91 and 0.95. The bootstrap of lame (bootstrap = N) draws its replicates the same way. The default dynamic empirical-Bayes fit of every family refuses a bootstrap, with an error of class "lame_dynamic_eb_unsupported", and carries standard errors in vcov() and confint() instead (Laplace ones for a binary or Poisson fit, the exact Gaussian marginal ones for a normal fit). A penalised binary or Poisson dynamic fit refuses it with an error of class "lame_als_bootstrap_unsupported". These refusals come before the replicate count is checked. Standard errors are replicate column standard deviations; confidence intervals use the percentile method.

Replicates. At least 10 valid replicates are needed (R below 10 is an error of class "lame_als_bootstrap_replicates"; bootstrap = N in ame_als and lame_als needs 30, and a 95 percent percentile interval is best read from 200 or more). A replicate whose refit has not converged is retried on the same data with four times the budget and, failing that, dropped: it never enters an interval. A replicate that errors, returns non-finite values or cannot carry the split-sample correction is dropped too, and every dropped replicate is counted (n_failed, with the reason in fail_reason). A dropped replicate is usually an extreme one, so the nominal level of a percentile interval can move by at most the share dropped; more than 5 percent of the R replicates dropped is an error of class "lame_als_bootstrap_convergence", raised at the replicate that crosses the cap, with no interval from the rest.

One optimum per replicate. The replicates of a fit with a multiplicative term (R > 0) are refit from the point estimate and also from the start the original fit used (the cold spectral start and, for a fit with multistart, its random starts); the converged refit with the lowest objective is kept, so a replicate cannot sit at a worse local optimum than the original procedure finds on its data. The start each replicate kept is in boot_refit_start, and optimum_check says what was compared. Without a multiplicative term the optimum is unique and a replicate is refit once. (An empirical-Bayes refit starts from the variational state of the fit and is not compared with other starts.) The point estimate of every coefficient must lie inside the central 95 percent of its own replicates: if it does not, the replicates are distributed around another optimum, or the estimator is biased by more than two standard errors, and the call is an error of class "lame_als_bootstrap_optimum" that carries the bootstrap as its boot field. A variance component whose estimate lies outside its replicate interval is the bias of its estimator in a regime a percentile interval cannot describe (the variance of the lognormal layer of a count network without extra-Poisson noise is the usual case): its standard error and interval are NA, it is named in vc_uncentred with the estimate and the replicate range, and the rest of the bootstrap stands.

Block bootstrap. A statistic that varies over B resampling blocks has a percentile interval that falls short of its level: in the simulations below a 95 percent interval covered the dyadic coefficient in 0.90 of the panels over ten to thirty blocks. The end points of a block interval are therefore read at the expanded percentile levels \Phi(\sqrt{B / (B - 1)}\, t_{B - 1, \alpha / 2}) (Hesterberg, 2015): 0.9 and 99.1 percent for a 95 percent interval at ten blocks, 1.6 and 98.4 percent at twenty. n_blocks records B, and confint(), tidy() and the interval columns of print() and summary() use these levels. The number of replicates matters more here than for a parametric bootstrap, since the end points sit in the tails of the replicates: use 200 or more.

A block resample redraws the time slices and holds the actors fixed, so it carries the variability of the data over time and none of the variability over actors. The parameters whose sampling variability is that over actors have no standard error and no interval under type = "block": they are NA, named in withheld with the reason, and print() and summary() say so. type = "parametric" gives them. They are the node-covariate coefficients, va, vb and cab, and, for an empirical-Bayes fit, the intercept (the population intercept) and the additive effects and factors, which are shrunk toward their prior mean: the replicates vary around the shrunk estimates, not around the effects. The residual variance ve has none either: a refit to a resample, which repeats some slices and omits others, absorbs more of the noise than the original fit, so its replicates sit below the estimate. The intercept of a fixed-effects fit is the mean of its linear predictor over the observed actors (the additive effects sum to zero), and its interval is for that quantity and not for a population intercept.

Coverage of the 95 percent interval was measured on directed normal-family panels of 15 to 40 actors and 10 to 30 periods with a time-varying dyadic covariate, static node covariates, random actor effects (standard deviation 0.6) and reciprocal errors (correlation 0.3), at block lengths 1 to 3, with 100 to 200 replicates a bootstrap and 80 to 150 panels a cell (1720 bootstraps of simulated panels in all, 530 of them with 200 replicates). Coverage of the 95 percent interval, over all panels with ten or more blocks (plain percentile), at the levels above with 200 replicates (expanded), and for the parametric scheme on the same panels:

The first five rows keep their intervals (for a fit by the fixed-effects estimator, bias_correct = "none"). The parametric column is the interval of the parametric scheme with its own corrections. The intercept of a fixed-effects fit is the mean of its linear predictor, so a replicate's intercept carries the projection of the redrawn actor effects on the node covariates, which the realized intercept does not: the replicate intercept is taken relative to it (the interval was 2.6 times the sampling standard deviation, coverage 1.00, before; now 0.92 of 500 panels, standard error 0.94 of the sampling standard deviation), and the node-covariate coefficients keep that variation (0.90). The variances va and vb are biased and their replicates reproduce the bias (a percentile interval covered 0.82 to 0.89 of the fixed-effects panels and 0.65 to 0.75 of the empirical-Bayes ones), so their parametric intervals are bias-corrected: the basic interval on the log scale for a fixed-effects fit (0.92 to 0.98 at 15, 20 and 40 actors) and the basic interval for an empirical-Bayes fit, whose estimate is biased upward by an amount that does not scale with it (0.91 to 0.97 at 20 and 30 actors). On binary empirical-Bayes panels of 20 actors and 12 periods (80 panels) the block intervals covered the dyadic coefficient in 0.91 of the panels and rho in 0.94, and the intercept in 0.19, the node-covariate coefficients in 0.25 and 0.26, va, vb and cab in 0.43, 0.28 and 0.33 and the additive effects in 0.67 and 0.63 (0.93 for the parametric scheme). The same holds for the empirical-Bayes normal fit, the default for the normal family (20 actors, 12 to 20 periods, 160 panels, 200 replicates): the block intervals covered the dyadic coefficient in 0.96 to 1.00 of the panels and rho in 0.93 to 0.98, and the intercept in 0.03 to 0.13, the node-covariate coefficients in 0.05 to 0.10, va, vb and cab in 0.02 to 0.18 and the additive effects in 0.27 to 0.39; ve covered 0.91 to 0.95 and has no interval. The fixed-effects rows above are the fits made with bias_correct = "none".

Fixed-effects fits. A binary fit by the fixed-effects estimator (bias_correct = "none") has no bootstrap: an error of class "lame_als_bootstrap_bias", raised before anything is fitted by ame_als and lame_als. The dyadic coefficient of such a fit carries the incidental-parameters bias, every replicate is refit at the biased estimate and reproduces it, and the percentile interval sits about twice the bias from the truth. A 95 percent interval covered the dyadic coefficient in 0.54, 0.72, 0.76 and 0.84 of 100 simulated panels of 30 actors and one period, 60 and two, 100 and two, and 30 and four; the bias relative to the standard error falls slowly with the number of actors (0.85, 0.71, 0.57, 0.40) and the fit passes the size rule of ame_als at 100 actors and two periods. The default empirical-Bayes estimator, or the split-sample or parametric-bootstrap correction, is the fit to bootstrap. The 95 percent interval of a corrected fit covered the dyadic coefficient in 0.940 (jackknife) and 0.940 (bootstrap correction) of 150 simulated panels of 30 actors and four periods, and in 0.945 and 0.953 of 145 and 150 panels of 60 actors and two (200 replicates; Monte Carlo standard error 0.019). Over those four cells the intervals for the dyadic coefficient, va and vb covered 0.944, 0.940 and 0.943 of 595 panels (standard error 0.010), and the standard error of the dyadic coefficient was 0.99 to 1.04 times its sampling standard deviation; 5 of the 600 panels were refused by the refit rule or the size rule. A replicate of a jackknife-corrected fit repeats the split-sample correction, and a replicate of a bootstrap-corrected fit has the stored bias removed. Scaling that stored bias by each replicate's own estimate changed the coverage of 60 panels in each of two cells by at most 0.02 and is not used; repeating the bootstrap correction in every replicate costs about 60 refits a replicate. A fixed-effects Poisson fit has no incidental-parameters bias in its dyadic coefficient (0.93 to 0.94) or its variance components, but its intercept does: it has no interval under either scheme (0.61, 0.85 and 0.73 of panels of 30 actors and one period, 60 and two, and 100 and two).

Usage

ame_als_bootstrap(
  object,
  R = 200,
  type = c("parametric", "block"),
  block_length = 1,
  seed = NULL,
  verbose = TRUE,
  max_iter = 200L,
  tol = 1e-06
)

boot_ame(
  object,
  R = 200,
  type = c("parametric", "block"),
  block_length = 1,
  seed = NULL,
  verbose = TRUE,
  max_iter = 200L,
  tol = 1e-06
)

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 identifies a column only when its singular value is separated from its neighbours': two adjacent singular values (eigenvalues, for a symmetric model) count as separated when their gap in the point estimate is at least 5 bootstrap standard deviations of that gap (sv_separation$threshold; the rotation the alignment applies to the pair is then under 4 degrees at the median and under 10 at the 90th percentile, against 9 and 25 when the two are equal). The U and V columns that are not separated from a neighbour have NA for the column-wise standard errors (se_U, se_V, and their unaligned _raw versions) and NA intervals in confint(). What is reported for every fit with a multiplicative term, separated or not, is the uncertainty of the subspaces the factors span, which a rotation of the columns does not change: the principal angles, in degrees, between each replicate's row-factor (and column-factor) space and the point estimate's (subspace). For a symmetric model the eigenvalues L of the multiplicative term are aggregated and reported as well.

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). The block bootstrap needs at least ten resampling blocks, so parametric (the default) is substituted, with a warning, on a panel that supplies fewer than ten blocks. The finite-sample (incidental-parameters) stretch of a binary IRLS fit is removed from the point estimate and from every replicate by the empirical-Bayes estimator (bias_correct, on by default for binary, Poisson and normal fits) or, when requested, the split-sample jackknife or the parametric-bootstrap correction; a binary fit that carries none of them has no bootstrap (Fixed-effects fits above). 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 the factors are static; NA for the columns whose singular values are not separated), sv_separation (the singular values of the point estimate, the gaps and their z-scores in bootstrap standard deviations, the threshold, and which columns are resolved) and subspace (the principal angles, in degrees, of the replicate factor spaces to the point estimate's, and their largest angle by side); n_valid, n_total, n_failed, fail_reason (per replicate: "", "nonconverged", "split", "error" or "nonfinite"), max_fail_share, boot_refit_start (the start each replicate kept), optimum_check, n_other_start (replicates that kept a start other than the warm one), vc_uncentred (the variance components with no interval because their replicates are not centred on the estimate, when there are any), withheld (a data frame of the parameters with no standard error or interval, their kind — "coefficient" or "vc" or "effects" — and the reason; NULL when every parameter has one), n_blocks (the number of resampling blocks of a block bootstrap, NA otherwise), estimator_eb (whether the fit is empirical Bayes, which sets the interval of va and vb), 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.

Hesterberg, T. C. (2015). What teachers should know about the bootstrap: resampling in the undergraduate statistics curriculum. The American Statistician 69(4), 371–386.

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

# twelve periods supply the ten resampling blocks the block scheme needs
Y <- replicate(12, { 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,
  bias_correct = NULL
)

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.

A fit from one of the dynamic ALS routes (lame(..., method = "als") with dynamic_uv, dynamic_ab or dynamic_beta) stores trajectories rather than a single value per parameter. Such a fit is refit through the dynamic fitter that produced it, seeded from the stored trajectories, and the result is again a dynamic fit of the same class. Z_new is not available on that route, which runs its own working-response loop.

A static fit carrying the split-sample jackknife (see ame_als) is refit with it: the refit on the new data is followed by its own two half-sample refits, warm-started from it with the same max_iter and tol, and the corrected values are returned, with the refit's uncorrected values in $raw.

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

Examples

ame_memory_settings(n_nodes = 12, R = 1)

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

Examples

set.seed(1)
n <- 12
Y <- matrix(rnorm(n * n), n, n); diag(Y) <- NA
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- ame(Y, R = 1, nscan = 20, burn = 5, odens = 1, verbose = FALSE, chains = 1)
usage <- ame_memory_usage(fit, detailed = FALSE)

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

ame_parallel(Y, n_chains = k, ...) is ame(Y, chains = k, ...): the same chains (seeds seed + 1000 * (k - 1)), the same drift check, between-chain gate and extensions, and the same pooled fit. It is kept for the older argument names; combine_method = "list" returns the fits of the chains instead of the pooled fit.

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. A call whose chains are not settled after the extensions stops with an error of class lame_chain_unsettled (see ame).

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

Examples

n <- 12
arr <- array(rnorm(n * n * 3), dim = c(n, n, 3))
rownames(arr) <- colnames(arr) <- paste0("a", 1:n)
actors <- rep(list(paste0("a", 1:n)), 3)
Y_list <- array_to_list(arr, actors, sliceLabel = paste0("t", 1:3))
length(Y_list)

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

Examples

set.seed(1)
n <- 12
Y <- matrix(rnorm(n * n), n, n); diag(Y) <- NA
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- ame(Y, R = 1, nscan = 20, burn = 5, odens = 1, verbose = FALSE, chains = 1)
draws <- as_draws(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(..., chains = N), two chains by default) 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 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)
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- lame(lapply(YX_bin_list$Y, function(y) 1 * (y > 0)), YX_bin_list$X, family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE, chains = 1)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  autoplot(fit)
}

Bipartite network helper functions

Description

Initialize bipartite starting values

Usage

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

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

Examples

n <- 12
Y <- lapply(1:3, function(t) matrix(rnorm(n * n), n, n))
check_format(Y)

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, Department of Political Science, Pennsylvania State University.

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. Conditional sampler health counts are summed when every chain supplies them. The original conditional counts and joint-move diagnostics are retained in sampler_diagnostics_by_chain; move diagnostics from a single chain are not presented as pooled values. If a legacy chain lacks the numerical/bound decomposition, the combined sampler_health is NULL; its available per-chain records remain accessible.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples

set.seed(1)
n <- 12
Y <- matrix(rnorm(n * n), n, n); diag(Y) <- NA
# each fit is one chain (chains = 1); combine_ame_chains() pools them
fit1 <- ame(Y, R = 1, nscan = 20, burn = 5, odens = 1, seed = 1, verbose = FALSE, chains = 1)
fit2 <- ame(Y, R = 1, nscan = 20, burn = 5, odens = 1, seed = 2, verbose = FALSE, chains = 1)
pooled <- combine_ame_chains(list(fit1, fit2))

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

Examples

set.seed(1)
n <- 12
Y <- matrix(rnorm(n * n), n, n); diag(Y) <- NA
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- ame(Y, R = 1, nscan = 20, burn = 5, odens = 1, verbose = FALSE, chains = 1)
fit_small <- compact_ame(fit)

Compute MCMC convergence diagnostics for multiple chains

Description

Computes 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

Details

rhat is the split rank-normalised R-hat of Vehtari et al. (2021), the larger of the bulk and tail values, as posterior::rhat() reports it. It is the statistic of the between-chain gate of ame and lame, computed on the first min(n_samples) draws of every chain, so for a two-chain fit it equals fit$rhat for every monitored coefficient and variance component, and it is the Rhat column of summary(). It is NA for a quantity that is constant in every chain and for chains with fewer than 8 draws. rhat_classic is the Gelman-Rubin statistic of the whole chains, without splitting or rank normalisation. The printed verdict compares rhat with the gate's threshold, 1.05.

For lame_robust fits, the posterior package is required. Diagnostics include static coefficients, time intercepts, actor paths and levels, estimated hyperparameters, every off-diagonal latent product, and persistent pair precisions. Unused rank-zero factor hyperparameters and fixed Gaussian precisions are excluded. Chains must have identical inputs and priors, different seeds, and equal lengths of at least four retained draws. Robust pass requires both rhat and rhat_classic below 1.05 and finite positive ESS, bulk ESS and tail ESS for every included parameter. The individual R-hat and effective sample size values are reported for every included parameter. The existing ame/lame diagnostic scope remains coefficients and variance components.

Value

A list with rhat, rhat_classic, ess, param_names (over BETA and VC), n_chains and n_samples. For robust fits, additionally returns ess_bulk, ess_tail, a full diagnostics data frame, and pass.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples

set.seed(1)
n <- 12
Y <- matrix(rnorm(n * n), n, n); diag(Y) <- NA
# each fit is one chain (chains = 1); compute_mcmc_diagnostics() compares them
fit1 <- ame(Y, R = 1, nscan = 20, burn = 5, odens = 1, seed = 1, verbose = FALSE, chains = 1)
fit2 <- ame(Y, R = 1, nscan = 20, burn = 5, odens = 1, seed = 2, verbose = FALSE, chains = 1)
diag_out <- compute_mcmc_diagnostics(list(fit1, fit2))
diag_out$rhat

Comtrade data

Description

Ten years (1996 to 2005) of bilateral trade among 30 countries in six SITC Rev. 1 commodity classes aggregated at the first level (AG1).

Format

A four-way array of dimension 30 x 30 x 6 x 10 (exporter, importer, commodity class, year) with dimnames giving the country names, the classes (chemicals; crude materials, inedible, except fuels; food and live animals; machinery and transport equipment; manufactured goods classified chiefly by material; miscellaneous manufactured articles) and the years. Each value is the change from the previous year in the natural logarithm of exports measured in 2000 US dollars. For example, the 1996 slice records changes from 1995 to 1996. Values are growth measures, not log levels or raw trade volumes; negative values indicate decreases. The 1,800 missing cells are the within-country diagonals. The first index is the reporting (exporting) country.

Source

Hoff, P. D. (2011), Section 4 of Separable covariance arrays via the Tucker product, with applications to multivariate relational data, doi:10.1214/11-BA606. Original processed array: https://www2.stat.duke.edu/~pdh10/Code/hoff_2011_ba/. Underlying trade and inflation data: 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); for an empirical-Bayes fit (the default for normal, binary and Poisson responses), for every coefficient from the Laplace covariance.

Usage

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

Arguments

Details

An empirical-Bayes fit (the default for normal, binary and Poisson fits) has Wald intervals from its Laplace covariance, which propagates the uncertainty of the actor effects and the latent factors and covers the node-covariate coefficients too.

The sandwich intervals of the other fits 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. The rows of the U and V columns whose singular values are not separated from a neighbour's, of a variance component with no interval (vc_uncentred), and of a parameter the bootstrap returns no interval for (withheld), are NA; see ame_als_bootstrap.

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 selected dynamic coefficient, compares the largest absolute first difference divided by its fitted innovation standard deviation \sigma_\beta \sqrt{\Lambda_{jj}} with a simulated reference. Each reference path has the same number of periods as the fitted path, starts from the first state of the matched posterior draw, and follows that draw's persistence and innovation scale: it is the largest jump a path that begins where the fitted path begins makes by chance. The design-based innovation scale and the averaging of symmetric node coefficients are retained. This is a posterior-conditioned transition reference, not an unconditional prior predictive check or a marginal-likelihood Bayes factor.

Usage

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

Arguments

Details

The bf column is the tail ratio \Pr(M^{post} > m^*) / 0.05, where m^* is the reference's 95th percentile. Its range is zero to twenty, in steps of 20 over the number of stored draws. Values above threshold_bf flag unusually abrupt posterior paths for inspection; they do not establish a change point. The reference supports "ar1" and "rw1" on equally spaced periods. RW2, Matern, irregular periods and actor-specific coefficient laws are refused rather than assigned an AR(1) reference. Saved named training designs, the dynamic mask, and aligned transition draws are required.

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

Measured properties

The table gives the percentage of simulated panels with bf > threshold_bf at the default threshold, from 200 fits of each path type and panel shape. Every panel has one standard-normal dyadic covariate whose coefficient follows the named path, is fitted with dynamic_beta = "dyad" (dynamic_beta_kind = "rw1" for random-walk paths, "ar1" otherwise) and keeps 40 stored draws (12 for the binary panels, which have 50 actors; the others have 15). The outcome is normal, binary or Poisson. Smooth paths are a constant, a linear ramp and two stationary AR(1) paths (persistence 0.7, innovation standard deviation 0.15 and 0.4); in-model paths are random walks and AR(1) paths drawn from the prior of the fit; a step is a jump of 0.6 or 1.2 (from 0.2 to 0.8 or to 1.4) at the middle of the panel. Normal outcomes have residual standard deviation 0.4; binary outcomes are probit with intercept -0.3; Poisson outcomes have log-rate offset -0.5 and no AR(1) paths drawn from the prior. Worst path type is the highest rate among the smooth and in-model path types of the row.

Pooled over panel shapes, smooth paths were flagged in 1.2 percent of the fits and in-model paths in 2.3 percent. At a threshold of 9 the pooled in-model rate is 5.0 percent, with no margin, and at 10 it is 2.3 percent, so 10 is the default. The tail ratio judges a jump against the innovation scale the fit estimates from the same path, so a panel of four or five periods rarely shows a single jump as unusual (about a fifth of the panels with a step of 1.2 are flagged), while panels of eight or more periods flag it in most cases.

Examples

data(YX_bin_list)
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- lame(lapply(YX_bin_list$Y, function(y) 1 * (y > 0)), YX_bin_list$X, family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE, chains = 1)
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", fit = NULL)

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

For the normal family mean_logdens is the mean over the scored cells of \log \phi(y; \hat y, v_e), the normal log density at the prediction with the fit's residual variance v_e: the posterior mean of VC[, "ve"] for an MCMC fit, the fitted VC["ve"] for a fast fit. It is NA when fit is not supplied.

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
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- ame(Y_train, R = 0, family = "binary",
           burn = 15, nscan = 60, odens = 5, verbose = FALSE, plot = FALSE, chains = 1)
y_pred <- predict(fit, type = "response")
evaluate_heldout(Y, y_pred, mask, family = "binary")

Extract fitted values from AME model

Description

Returns the stored response-scale summary (YPM). For binary models, these are predicted probabilities between 0 and 1. For normal models, these are predicted continuous values. For MCMC Poisson models, these are mean counts, \exp(\eta + \sigma^2 / 2), at the posterior mean \eta of the linear predictor and the posterior mean \sigma^2 of the residual variance, the form the fast estimators report. Predictive quantiles and event probabilities come from predict(type = "distribution").

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, plug_in = NULL, ...)

Arguments

Details

For a binary empirical-Bayes fit (the default) the fitted values are predictive probabilities: the expectation of the response over the posterior of the linear predictor, \Phi(m / \sqrt{1 + v}) for a probit link, with m the posterior mean and v the posterior variance of the predictor (see predict.ame_als). For a Poisson fit they are the plug-in mean counts, e^{\eta + v_e / 2} at the fitted predictor \eta (the posterior mean) and the variance v_e of the noise of the count, as the samplers report them and as the dynamic empirical-Bayes fit does: the posterior predictive mean of a count is not finite in this model (the inverse-gamma tail of v_e, and the inverse-Wishart tail of the scale of the factors for a cell the data do not inform), so its value would be a property of the numerical integration and of the prior; the plug-in mean is finite everywhere. plug_in = FALSE returns the predictive mean of a count that the integration over the predictor gives, e^{m + (v + v_e) / 2}, and plug_in = TRUE the plug-in values of a binary fit (the response mean at the fitted predictor); the plug-in values are what fit$fitted and fit$YPM hold and simulate() uses. residuals() follows fitted(). The mean of a normal response is its predictor, so a normal fit's fitted values are the posterior means of the predictor, with or without plug_in.

A dynamic fit (lame(method = "als") with dynamic_ab, dynamic_beta or dynamic_uv) follows the same rules. The empirical-Bayes binary fit returns the predictive probability, with v the variance of the cell's predictor: the additive effects given the coefficients (the prior variance for an actor with no observed cell), the product of the positions, and the coefficients' part at the design asked for; plug_in = TRUE returns \Phi(m), which fit$fitted keeps. A fit stored without the variance record is an error of class "lame_predictive_missing". The empirical-Bayes count fit returns the plug-in mean e^{\eta + v_e / 2}, the normal fit the predictor, and a penalised fit (bias_correct = "none") the response mean at its predictor (\Phi(\eta), e^{\eta}). The dynamic route has no exposure and no forecasts (predict(h > 0) is an error).

Value

A matrix or list of matrices of fitted values.

Extract fitted values from LAME model

Description

Returns the stored response-scale summaries. For MCMC Poisson models these are mean counts, e_t \exp(\eta + \sigma^2 / 2), at the posterior mean \eta of the linear predictor and the posterior mean \sigma^2 of the residual variance, with e_t the exposure of the period (period_exposure; one without); the form the fast estimators report. Predictive quantiles and event probabilities come from simulate() and from predict(fit, h = 1, interval = "prediction").

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,
  newdata = NULL,
  newdesign = NULL,
  newexposure = 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 averages the Gaussian CDF over matched predictor and residual-variance draws. For discrete families ("binary", "poisson", "ordinal") the randomised PIT of Czado, Gneiting & Held (2009) is used. Poisson forecasts integrate the Gaussian observation effect and future exposure. Ordinal scoring requires saved explicit posterior cutpoints; data-induced joint rank simulation does not define a fixed future threshold law. For rank families ("cbin", "frn") PIT is not implemented; the function returns NULL with an informational note. Robust fits from lame_robust() use their genuine future latent states, matched persistent pair precisions and Gaussian or transient Student residual CDF. Fitted temporary observation precisions are never reused for future outcomes.

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 (NA; a dependent-network null reference is unavailable), cover_95 (fraction of PIT values in [0.025, 0.975]), matched draw_index, predictive cdf_lower/cdf_upper, and scored actor/period cells. Network cells need not be independent. The KS statistic, histogram and coverage are descriptive. No global model-adequacy p-value is reported without a calibrated joint network reference. ks_inference records this restriction.

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 <- lapply(YX_bin_list$Y, function(y) 1 * (y > 0))[1:3]
Y_test <- lapply(YX_bin_list$Y, function(y) 1 * (y > 0))[4]
X_train <- YX_bin_list$X[1:3]
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- lame(Y_train, Xdyad = X_train,
            family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE, chains = 1)
pit <- forecast_pit(fit, y_future = Y_test)
pit$ks_stat
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.

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,
symmetric = FALSE)

Arguments

Value

returns list of design array values necessary for ame_repL, plus intercept: the intercept flag actually used (FALSE when the covariates span a constant across the observed panel). A cell whose covariate is missing is set to NA in the returned Y and its design value to 0.

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,
  L_dyn = NULL,
  period_exposure = 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,
  bip = FALSE,
  unit_var = 1
)

Arguments

Details

In addition to the slots listed below, any of the dynamic hyperparameter scalars rho_uv, sigma_uv, rho_ab, and sigma_ab present in start_vals (for example from als_start_vals) are passed through so a dynamic MCMC run can warm-start them, as are the adapted step sizes level_move_sd and dev_move_sd of the binary scale move, so that a run continued with lame_resume keeps the step it had tuned rather than re-adapting from the default.

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.

Examples

set.seed(1)
n <- 12
Y <- matrix(rnorm(n * n), n, n); diag(Y) <- NA
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- ame(Y, R = 1, nscan = 20, burn = 5, odens = 1, verbose = FALSE, chains = 1)
glance(fit)

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)
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- lame(lapply(YX_bin_list$Y, function(y) 1 * (y > 0)), YX_bin_list$X, family = "binary", R = 0,
            nscan = 100, burn = 20, odens = 5, verbose = FALSE, chains = 1)
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_nrm)
fit <- ame_als(Y = YX_nrm$Y, Xdyad = YX_nrm$X[, , 2:3],
               family = "normal", 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

Gaussian bipartite checks use weighted four.cycle.dep; a nonzero count would be fixed by the missing-cell pattern. Each replicate is evaluated on the same observed cells as its corresponding data slice, aligning actor labels before applying the mask. Named data periods are matched to the fit; unnamed periods must be in fitted order. All custom functions receive individual matrices for both observed and replicated data, and their values are averaged over periods with equal weight.

For the ordinal family the posterior predictive draws are quantile-matched to the observed category frequencies (the cutpoints are induced from the data), so the simulated marginal distribution always equals the observed one: the GOF statistics and YPM cannot show a marginal misfit for ordinal outcomes and only assess the dependence structure (degree spread, dyadic and triadic dependence).

This function requires joint posterior replication and refuses fits that lack aligned saved draws; it never substitutes conditional mean trajectories. 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, thin_UV = 1, thin_ab = 1). At rank zero, omit save_UV.

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 (nsim + 1) x n_stat matrix of GOF statistics: the first row holds the observed statistics, the remaining rows the posterior predictive ones. For an ame fit this is the format of fit$GOF; a lame fit stores its own goodness of fit as a per-statistic list of period-by-draw matrices instead, so the two are not interchangeable. Attributes replication_method and draw_index record joint replication and the sampled MCMC rows.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

ame, lame, gof_plot

Examples

# Run model without GOF during fitting
data(YX_nrm)
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, R = 2, gof = FALSE,
           nscan = 100, burn = 10, odens = 1, verbose = FALSE,
           posterior_opts = posterior_options(save_ab = TRUE, save_UV = TRUE,
                                              thin_ab = 1, thin_UV = 1), chains = 1)

# 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. Used for families with structural zeros, including binary outcomes.

four.cycle.dep (Weighted bipartite dependence)

For Gaussian outcomes, averages products of centered, standardized cell values around fully observed rectangles. This measures weighted dependence; a nonzero-edge count would be fixed by the observed-cell pattern.

Interpretation Guide:

Histogram plots (static models):

• Dashed black vertical line: observed statistic value (dual-encoded on colour and linetype)

• Grey histogram: distribution from posterior predictive simulations

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

• Poor fit: dashed 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)
# chains = 1 keeps the example short; the default runs two chains and compares them
fit_ame <- ame(YX_nrm$Y, Xdyad = YX_nrm$X, gof = TRUE,
               nscan = 100, burn = 10, odens = 1, verbose = FALSE, chains = 1)

# 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, replacing the count with weighted four.cycle.dep for Gaussian outcomes) on each replicate, and compares the observed value with the replicates.

Usage

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

Arguments

Details

A cross-sectional fit gives one histogram per statistic with the observed value overlaid. A longitudinal fit gives, as gof_plot does for an MCMC lame fit, the statistics of every time slice: for each statistic, the observed value in each period against the median and the central 95 percent of the simulated values in that period.

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.

Examples

set.seed(1)
n <- 12
Y <- matrix(rnorm(n * n), n, n); diag(Y) <- NA
fit <- ame_als(Y, R = 1, family = "normal", verbose = FALSE)
gof_plot(fit, nsim = 20, seed = 1)

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, family = 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", family = "normal")

# 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, family = NULL)

Arguments

Details

For Gaussian outcomes, four.cycle.dep is the mean product of four globally centered, sample-standardized values around each fully observed rectangle. It measures dependence among edge magnitudes. It returns NA when no complete rectangle is observed or the observed values have zero variance. Missing cells are not zeros. A nonzero-edge count is unsuitable for Gaussian predictive checks: continuous draws are nonzero almost surely, so fixing the observation mask fixes their four-cycle count regardless of model fit.

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.cycle.dep

For Gaussian outcomes, weighted dependence over complete observed rectangles, as defined below. Replaces the count.

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, family = "normal")

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 min(1, 2 * min(p_up, p_down)), where both empirical tails include ties. The value runs from 0 (observed slope in the extreme tail) to 1 (observed slope central in the predictive distribution).

Usage

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

Arguments

Details

Observed and replicated statistics use the same actor labels and observed cells in every period. Transitivity counts only fully observed triples. Joint posterior draws are required; missing storage is refused. Fit with posterior_options(save_ab = TRUE, save_UV = TRUE, thin_ab = 1, thin_UV = 1), omitting latent storage at rank zero. A constant observed series is still tested when the predictive slopes vary.

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

• n_rep (number of replicate networks drawn),

• degenerate (TRUE when the replicate slopes have no spread or a statistic is undefined, in which case p_pp is NA),

• reason and n_undefined describe an unestimated check,

• replication_method and draw_index identify joint draws

Examples

data(YX_bin_list)
# chains = 1 keeps the example short; the default runs two chains and compares them
fit <- lame(lapply(YX_bin_list$Y, function(y) 1 * (y > 0)), YX_bin_list$X, family = "binary", R = 0,
            dynamic_beta = "dyad",
            posterior_opts = posterior_options(save_ab = TRUE, thin_ab = 1),
            nscan = 60, burn = 15, odens = 5, verbose = FALSE, chains = 1)
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

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,
  als_lowrank_method = NULL,
  als_linear_solver = NULL,
  als_multistart = NULL,
  multistart = c("none", "cheap", "full"),
  bias_correct = NULL,
  bootstrap = 0L,
  bootstrap_type = c("parametric", "block"),
  bootstrap_block_length = 1L,
  bootstrap_seed = NULL,
  chains = .lame_chains_default(),
  cores = NULL,
  resume_from = NULL,
  print,
  sampler = c("gibbs", "tempered")
)

Arguments