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

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

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

Usage

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

Arguments

Y: a square (unipartite) or rectangular (bipartite) relational matrix. A cross-sectional netify object is also accepted and converted with netify::to_lame(); when family or mode is omitted, the value inferred by netify is used.
Xdyad: an n_row x n_col matrix or n_row x n_col x pd array of dyadic covariates, or NULL.
Xrow: an n_row x pr matrix of row/sender covariates, or NULL.
Xcol: an n_col x pc matrix of column/receiver covariates, or NULL.
R: integer dimension of the multiplicative effects (default 0). The covariate coefficients are conditional on this choice. The multiplicative term $u_i'v_j$ is a flexible high-variance regressor that can correlate with the dyadic covariates, so the estimated beta can shift – and occasionally change sign – as R increases. Comparing against an R = 0 fit is a useful check on whether the covariate story is being driven by the latent rank.
family: one of "normal", "binary", or "poisson". The rank and censoring families are MCMC-only.
mode: "unipartite" (square) or "bipartite" (rectangular).
symmetric: logical; fit a symmetric (undirected) model. Unipartite only.
max_iter: maximum number of block coordinate descent iterations (default 200).
tol: convergence tolerance on the relative change in residual sum of squares (default 1e-6).
lowrank_method: inner solver for the multiplicative (low-rank) block: "mm" (default) weighted majorise-minimise; "als" alternating least squares; "hybrid" runs both and keeps the lower-objective result. "als"/"hybrid" converge faster than "mm" on strongly unbalanced longitudinal panels and are available for directed and bipartite models (symmetric fits always use "mm"). All three minimise the same objective, so the point estimate is unchanged for balanced data.
non_normal_method: for the non-normal ALS families, "irls" (default for binary/poisson) runs iteratively reweighted least squares, giving a fast approximate GLM AME fit with coefficients on the requested link scale (Poisson log, binary logit/probit). "transform" fits one fixed Gaussian working response (log(y+1) for Poisson, rank-normal scores for binary); its coefficients are on an uncalibrated working scale and are mainly useful for direction/ranking checks. A directed R > 0 IRLS fit uses the hybrid low-rank solver internally because the IRLS weights are unbalanced. Uncertainty for either path comes from the bootstrap or sandwich covariance.
link: link for non_normal_method = "irls" with a binary family: "probit" (default; matches ame/lame) or "logit". poisson always uses the log link; ignored otherwise.
linear_solver: solver for the regression block: "eigen" (default) eigendecomposes the normal equations; "qr" uses a QR factorisation of the observed design, which is more stable for ill-conditioned covariates; "auto" picks "qr" when the design is ill-conditioned but full rank. All give the same answer for well-conditioned designs.
multistart: for R > 0 (a non-convex objective), "none" (default) fits from a single deterministic start; "cheap" (4 starts) and "full" (8 starts) also try random low-rank starts and keep the lowest-SSE fit, warning when the starts reach materially different optima. Reproducible given seed; the global RNG stream is left unchanged.
bootstrap: integer: if > 0, additionally run bootstrap replicates of the parametric or block bootstrap (via ame_als_bootstrap) after the point fit and attach the result as fit$bootstrap. The downstream accessors (confint.ame_als, summary, print) then surface bootstrap intervals instead of the anti-conservative sandwich Wald intervals. Default 0 (no bootstrap) – bootstrap is expensive, so it is not imposed on a user who just wants a quick fit.
bootstrap_type: character: "parametric" (default) or "block" – the bootstrap scheme to use when bootstrap > 0.
bootstrap_block_length: integer: block length for the block bootstrap; only used when bootstrap_type = "block".
bootstrap_seed: optional integer seed for the bootstrap (the point fit uses seed).
verbose: logical; print progress (default TRUE).
seed: random seed (default 6886). The block coordinate descent is deterministic, so the point estimate is reproducible regardless; the argument is retained for API consistency with ame.

Value

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

va, vb: empirical variances of the sender and receiver additive effects.
cab: covariance of the sender and receiver effects (NA for symmetric or bipartite models).
rho: dyadic residual reciprocity — the correlation between the residuals of $(i,j)$ and $(j,i)$ — not the sender/receiver correlation cor(a, b).
ve: residual variance, on the model's own scale and degrees-of-freedom-corrected ($\mathrm{SSE}/(n_{obs} - \mathrm{df})$) so it matches the MCMC posterior-mean s2. For a binary IRLS fit the probit/logit model fixes the latent error variance at 1 by identification, so ve is reported as 1; the working-scale GLM dispersion (about 1 under correct specification) is kept separately as ve_working.

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

Details

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

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

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

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

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

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

Coverage relative to `ame` / `lame`

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

Families. ALS supports normal, binary, and poisson. For ordinal, cbin, and frn, fall back to ame / lame.
Dynamic effects. lame_als itself fits a static model pooled across time slices. The top-level dispatcher lame(..., method = "als") routes supported dynamic requests to a dynamic point estimator for normal, binary, and poisson panels, including named panels where actors enter or exit: dynamic_ab, selected intercept/dyadic/node dynamic_beta, and AR(1) or Student-t dynamic_uv for directed, symmetric, and bipartite panels. The snap-only dynamic_uv = TRUE, dynamic_uv_kind = "snap" case routes to lame_snap_als for supported normal unipartite and bipartite panels. Node-covariate coefficients use the same orthogonal additive-effect decomposition as lame_als; dynamic node coefficients use period-specific node values when selected by dynamic_beta, while static node coefficients use per-actor means. bipartite dynamic_g is available on the dynamic als path for normal, binary, and poisson panels. changing actor composition requires row and column names on every slice so actors can be aligned; smoothing penalties are broken across actor-entry gaps. rank/censored dynamic families remain on the mcmc path. static ame_als() / lame_als() fits still use a single latent rank R; the dynamic bipartite ALS dispatcher honours separate R_row and R_col values.
Priors. ALS has no priors. prior = list(...) and g = ... are MCMC-only; the dispatcher ame(..., method = "als") warns and ignores them.
Posterior quantities. No $BETA / $VC posterior draws on the point fit. Use bootstrap = N for a sampling distribution analogue; the draws are bootstrap, not Bayesian.
Multi-chain. Not applicable; n_chains is dropped.
Convergence diagnostics. No Rhat / ESS / trace_plot.

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

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

Author

Cassy Dorff, Shahryar Minhas, Tosin Salau

Examples

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