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.
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
)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, vbempirical variances of the sender and receiver additive effects.
cabcovariance of the sender and receiver effects
(NA for symmetric or bipartite models).
rhodyadic residual reciprocity --- the correlation between
the residuals of \((i,j)\) and \((j,i)\) --- not the sender/receiver
correlation cor(a, b).
veresidual 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.
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.
an n_row x n_col matrix or n_row x n_col x pd
array of dyadic covariates, or NULL.
an n_row x pr matrix of row/sender covariates, or NULL.
an n_col x pc matrix of column/receiver covariates, or NULL.
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.
one of "normal", "binary", or
"poisson". The rank and censoring families are MCMC-only.
"unipartite" (square) or "bipartite" (rectangular).
logical; fit a symmetric (undirected) model. Unipartite only.
maximum number of block coordinate descent iterations (default 200).
convergence tolerance on the relative change in residual sum of
squares (default 1e-6).
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.
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 for non_normal_method = "irls" with a binary
family: "probit" (default; matches ame/lame)
or "logit". poisson always uses the log link;
ignored otherwise.
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.
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.
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.
character: "parametric" (default) or
"block" -- the bootstrap scheme to use when bootstrap > 0.
integer: block length for the block
bootstrap; only used when bootstrap_type = "block".
optional integer seed for the bootstrap (the point
fit uses seed).
logical; print progress (default TRUE).
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.
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.
Cassy Dorff, Shahryar Minhas, Tosin Salau
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.)
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()).
lame_als for longitudinal networks,
ame_als_bootstrap for bootstrap uncertainty,
ame for the full MCMC estimator.
Y <- matrix(rnorm(400), 20, 20); diag(Y) <- NA
fit <- ame_als(Y, R = 1, family = "normal", verbose = FALSE)
coef(fit)
Run the code above in your browser using DataLab