--- title: "Fast (MCMC-free) AME Estimation" author: "Cassy Dorff, Shahryar Minhas, and Tosin Salau" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 fig_width: 7 fig_height: 5 vignette: > %\VignetteIndexEntry{Fast (MCMC-free) AME Estimation} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(lame) ``` ## When to use the fast estimator `ame()` and `lame()` fit AME models by Bayesian MCMC, which gives calibrated posterior inference but can take minutes. For rapid model exploration, trying several ranks, screening covariates, or getting starting values, the package also ships a fast, MCMC-free point estimator: - `ame_als()`: cross-sectional networks - `lame_als()`: longitudinal (replicated) networks - `ame_als_bootstrap()`: bootstrap standard errors and intervals It fits the same model z_ij = mu + beta'x_ij + a_i + b_j + u_i'v_j + e_ij by **iterative block coordinate descent** rather than Gibbs sampling. It is typically tens to hundreds of times faster than `ame()`, and the point estimates are usually close. Use it when a point estimate is enough; use `ame()` / `lame()` when the target is a posterior summary. ### Wall-clock benchmark The numbers below are indicative for normal-family AME with `R = 2`; absolute timings depend heavily on CPU, BLAS, and whether you are on a high-clock workstation or a virtualised / WSL host. On the dev box (WSL2, single-threaded BLAS) the same benchmark ran roughly 2-4× slower than the workstation reference column, so treat these as order-of- magnitude. The point estimator scales O(n² · iter) per BCD pass and is dominated by the linear system at each block; the MCMC scales similarly per iteration but pays the constant of running thousands of iterations. The "ratio" column is what you actually feel: the ALS fit completes before MCMC finishes its burn-in for any reasonable n. | n | `ame()` MCMC (4000 iter) | `ame_als()` | Speed ratio | |-----|--------------------------|-------------|-------------| | 50 | 5-25 s | 0.1-0.5 s | 25-250× | | 100 | 15-50 s | 0.3-1 s | 30-200× | | 200 | 60-120 s | 2-3 s | 30-60× | | 500 | 7-15 min | 15-30 s | 25-40× | These ranges are wide on purpose: under a single-threaded BLAS on a virtualised host, the small-n speed ratios we have measured are several times larger than the "≈ 30×" figure in older docs, because the MCMC's per-iteration overhead is fixed and the ALS solve is nearly free at n ≤ 100. Treat the ratio as "ALS finishes before MCMC's burn-in" rather than a number to plan against. For binary/poisson the **point-estimate** speedup ratio is much smaller (often only 3-10× on small networks) because the IRLS path runs 3-5 reweighting iterations per fit, and because the MCMC binary sampler is itself cheaper per iteration than the normal sampler. Bipartite and unipartite ALS run at broadly similar cost for comparable dimensions; which is faster depends on the exact shape and hardware, so do not read much into small differences. **The bootstrap can dominate.** A `bootstrap = 100` parametric bootstrap runs the full ALS fit 101 times. For a binary `ame_als(R = 1, bootstrap = 100)` at n = 50, we measured ≈ 135 s -- slower than the equivalent MCMC fit. If you only need the point estimate, skip the bootstrap (`bootstrap = 0`, the default) and use the sandwich `vcov()`/`confint()` on the fit. If you need interval estimates at scale, compare the bootstrap cost with a regular `ame()` MCMC run; for some binary networks the MCMC fit is cheaper -- and for binary `R > 0` it is also the statistically preferred choice, since the bootstrap inherits the point estimator's finite-sample upward bias (see Non-normal families below). To benchmark on your own hardware: ```{r benchmark_self, eval = FALSE} set.seed(1) n <- 100 Y <- matrix(rnorm(n*n), n, n); diag(Y) <- NA rownames(Y) <- colnames(Y) <- paste0("a", sprintf("%03d", 1:n)) t_als <- system.time(ame_als(Y, R = 2, family = "normal", verbose = FALSE)) t_mcmc <- system.time(ame(Y, R = 2, family = "normal", burn = 500, nscan = 4000, odens = 25, verbose = FALSE, plot = FALSE)) cat("ALS: ", round(t_als["elapsed"], 2), "s\n", "MCMC: ", round(t_mcmc["elapsed"], 2), "s\n", "ratio:", round(t_mcmc["elapsed"] / t_als["elapsed"], 1), "x\n") ``` ### What ALS can do vs. what only MCMC does ALS is a fast point estimator. It covers most static AME workflows and several dynamic workflows, while posterior-specific features still require the Bayesian MCMC fits. The current coverage is: | Capability | `ame()` / `lame()` (MCMC) | Fast point estimators | |-------------------------------------------|---------------------------|-----------------------| | Families: normal, binary | yes | yes | | Family: poisson | yes | yes | | Families: ordinal / cbin / frn | yes | no (error) | | Unipartite / bipartite mode | yes | yes | | `symmetric = TRUE` | yes | yes | | `R = 0..min(n)-1` | yes | yes | | Xrow / Xcol / Xdyad covariates | yes | yes | | Reproducibility via `seed=` | yes | yes (deterministic ALS) | | `coef`, `vcov`, `fitted`, `residuals` | yes | yes | | `confint` | posterior quantile | sandwich Wald OR bootstrap | | `predict(newdata=)` | yes | yes | | `summary`, `print`, `latent_positions` | yes | yes | | `nobs` | yes | yes | | `simulate` | yes | yes | | `gof_plot` | posterior predictive | bootstrap-based GOF check | | `ab_plot`, `uv_plot` | yes | yes | | `prior_summary` | yes | yes (reports "no priors") | | `dynamic_ab` / selected `dynamic_beta` | yes | normal, binary, and Poisson panels via dynamic ALS; named changing-composition panels are aligned | | Dynamic node-covariate coefficient paths | yes | period-specific node values when selected by `dynamic_beta`; static node effects use actor means | | `dynamic_uv` (smooth AR(1) / t) | yes | AR(1) and t for directed, symmetric, and bipartite panels | | `dynamic_G` | yes | bipartite normal, binary, and Poisson panels via dynamic ALS; named changing-composition panels are aligned | | `dynamic_uv_kind = "snap"` | yes | `lame_snap_als()` for normal unipartite and bipartite panels | | Custom priors (`prior = list(...)`, `g`) | yes | no Bayesian priors; dynamic ALS reads selected `prior$rho_*_mean` and `prior$lambda_*_als` values as smoothing controls | | Multi-chain (`n_chains`, via `ame_parallel`) | yes | no (use `bootstrap = N`) | | `posterior_opts` (save U/V/a/b samples) | yes | no (point estimator) | | `trace_plot` (Rhat / ESS) | yes | no (deterministic) | | `custom_gof`, `periodic_save` | yes | no | When you call the unified front door with MCMC-only arguments, the ALS dispatcher warns and lists exactly which ones it ignored -- for the arguments each entry point actually accepts. Two notes on the edges: - `nscan` / `burn` / `odens` / `g` are accepted-then-ignored with a warning by both `ame(method = "als")` and `lame(method = "als")`. - The static ALS dispatcher warns and ignores MCMC tuning arguments such as `nscan`, `burn`, `odens`, prior settings not used by the fast path, and `g`. Dynamic ALS has no posterior priors, but it does use selected entries in `prior` as deterministic smoothing controls (`prior$rho_ab_mean`, `prior$rho_beta_mean`, `prior$rho_uv_mean`, `prior$rho_G_mean`, and `prior$lambda_*_als`). Dynamic requests are stricter: supported `dynamic_ab`, selected `dynamic_beta`, AR(1) or t `dynamic_uv`, and bipartite `dynamic_G` requests route to dynamic ALS for normal, binary, and Poisson panels. Named changing-composition panels are aligned to the union actor set and actor-entry gaps break the smoothing penalties. Student-t dynamic UV fits attach final local transition-weight matrices as `fit$lambda_u` and `fit$lambda_v`. `lame(method = "als", dynamic_uv = TRUE, dynamic_uv_kind = "snap")` routes to `lame_snap_als()` for supported normal unipartite and bipartite panels; unsupported dynamic stacks raise a specific error rather than being silently ignored. Use `als_max_iter`, `als_tol`, and `als_stability` for longer dynamic ALS runs and start-sensitivity checks. - Passing longitudinal-only arguments such as `dynamic_uv`, `dynamic_ab`, or `dynamic_beta` to `ame()` (cross-sectional) raises the usual "longitudinal-only argument" error before ALS runs, because they are not valid `ame()` arguments at all. - `n_chains` is an `ame()` / `ame_parallel()` argument; it is not a `lame()` formal, so `lame(method = "als", n_chains = ...)` raises an "unused argument" error rather than a warn-and-ignore. Use `bootstrap = N` for ALS uncertainty regardless of entry point. ## A first fit ```{r} set.seed(1) n <- 40 a <- rnorm(n, 0, 0.5); b <- rnorm(n, 0, 0.5) Xd <- matrix(rnorm(n * n), n, n) Y <- 0.5 + 0.8 * Xd + outer(a, b, "+") + matrix(rnorm(n * n), n, n) diag(Y) <- NA fit <- ame_als(Y, Xdyad = Xd, R = 1, family = "normal", verbose = FALSE) fit ``` `coef()`, `fitted()`, `residuals()`, `predict()` and `plot()` all work as for any model object. The printed variance components (`va`, `vb`, `cab`, `rho`, `ve`) are descriptive summaries of the fitted effects, not random-effect variance components. ## Choosing the rank R There is no information criterion for `R`. Fit a few values and look for an elbow in the deviance / SSE: ```{r} devs <- sapply(0:3, function(r) ame_als(Y, Xdyad = Xd, R = r, family = "normal", verbose = FALSE)$deviance) data.frame(R = 0:3, deviance = round(devs, 1)) ``` If the curve declines smoothly with no clear elbow, keep `R` small and confirm the chosen rank with an `ame()` fit. ## Non-normal families For `binary` and `poisson` data the default `non_normal_method = "irls"` runs an iteratively reweighted least squares loop, returning coefficients on the calibrated link scale (probit for binary, log for poisson). The alternative `non_normal_method = "transform"` is faster but its coefficients are on an uncalibrated rank scale, good for the sign and ranking of effects but not their magnitude. One caveat for binary fits with `R > 0`: the factor block is a penalized (MAP-style) point estimate, and a residual finite-sample (incidental-parameters) upward bias in the coefficients remains -- roughly +10-15% at n = 50, shrinking with n, and largest when `R` is set higher than the data support. The bootstrap reproduces rather than removes this bias, and `ame_als()` prints a reminder when `verbose = TRUE`. The binary examples in this vignette use `R = 0`, which does not trigger this caveat; for final binary coefficient inference at `R > 0`, use the MCMC `ame()`. ```{r} set.seed(2) Yb <- 1 * (Y > median(Y, na.rm = TRUE)); diag(Yb) <- NA fb <- ame_als(Yb, Xdyad = Xd, R = 0, family = "binary", verbose = FALSE) coef(fb) ``` ## Uncertainty The point estimator carries no uncertainty. Use the bootstrap: ```{r} bt <- ame_als_bootstrap(fb, R = 100, type = "parametric", seed = 1, verbose = FALSE) summary(bt) ``` The `Boot SE` column is the raw standard deviation across replicates, whereas the interval reported by `confint()` is percentile-based by default. A small number of unstable binary IRLS refits can inflate the raw standard deviation without moving the percentile limits much, so use the percentile interval rather than constructing a Wald interval from `Boot SE`. One row of the variance-components table is constant by construction: for the binary family the latent residual variance `ve` is fixed at 1 by probit identification, so its bootstrap SE is exactly zero and its interval is the point [1, 1] -- that row is expected, not a failed bootstrap. See `?ame_als` for the `ve` vs `ve_working` distinction. `type = "parametric"` (the default) simulates fresh outcomes from the fitted model and refits; `type = "block"` resamples time slices and is available for longitudinal fits with several time points. `confint()` and `vcov()` work on the bootstrap object. A fast analytic alternative for the **regression coefficients only** is the conditional sandwich covariance, `vcov()` / `confint()` applied directly to the fit. It is anti-conservative (it holds the additive and multiplicative effects fixed); the bootstrap is the recommended tool for full inference. ### ALS bootstrap vs MCMC posterior on the same data A useful sanity check on the fast estimator is to fit the same data with both engines and overlay the coefficient intervals. With a calibrated link (`non_normal_method = "irls"` for binary, the default) the ALS bootstrap should give a coefficient story qualitatively similar to the MCMC posterior on cases the model fits well; large differences mean the posterior fit is the better reference for that dataset. ```{r als-vs-mcmc, fig.width = 7, fig.height = 3.5, fig.alt="Side-by-side coefficient plot comparing the ALS bootstrap and MCMC posterior on the same binary network. Each coefficient (intercept and the single dyadic covariate) has two horizontal intervals: one for the ALS bootstrap, one for the MCMC posterior. A dashed vertical line at zero marks the no-effect reference; intervals on the same side of zero and with overlapping range indicate the two estimators agree."} # wrap Xd into a 3-D array with an explicit slice name so the ALS and MCMC # paths produce identical coefficient names (the auto-naming default differs # between the two engines) Xd_arr <- array(Xd, dim = c(nrow(Xd), ncol(Xd), 1), dimnames = list(NULL, NULL, "Xd")) # refit the binary ALS fit with the bootstrap attached so tidy() returns # bootstrap-based intervals fb_b <- ame_als(Yb, Xdyad = Xd_arr, R = 0, family = "binary", verbose = FALSE, bootstrap = 100, bootstrap_seed = 1) # fit the same data with MCMC using a short chain for a quick comparison fit_mcmc <- ame(Yb, Xdyad = Xd_arr, family = "binary", R = 0, burn = 30, nscan = 100, odens = 5, verbose = FALSE, plot = FALSE, gof = FALSE, seed = 1) # collect the common coefficient interval columns without relying on # vignette-time s3 dispatch through broom coef_interval_table <- function(fit) { est <- coef(fit) ci <- confint(fit) terms <- intersect(names(est), rownames(ci)) data.frame( term = terms, estimate = unname(est[terms]), conf.low = unname(ci[terms, 1]), conf.high = unname(ci[terms, 2]), row.names = NULL ) } als_tdy <- coef_interval_table(fb_b) mcmc_tdy <- coef_interval_table(fit_mcmc) df <- rbind( cbind(estimator = "ALS bootstrap", als_tdy[, c("term", "estimate", "conf.low", "conf.high")]), cbind(estimator = "MCMC posterior", mcmc_tdy[, c("term", "estimate", "conf.low", "conf.high")]) ) library(ggplot2) ggplot(df, aes(x = estimate, y = term, colour = estimator, shape = estimator)) + geom_vline(xintercept = 0, linetype = "dashed", colour = "grey60") + geom_pointrange(aes(xmin = conf.low, xmax = conf.high), position = position_dodge(width = 0.45), size = 0.5) + scale_colour_manual(values = c("ALS bootstrap" = "#0072B2", "MCMC posterior" = "#D55E00")) + scale_shape_manual(values = c("ALS bootstrap" = 16, "MCMC posterior" = 17)) + labs(x = "coefficient (probit scale)", y = NULL, colour = NULL, shape = NULL) + theme_bw() + theme(panel.border = element_blank(), axis.ticks = element_blank(), legend.position = "top") ``` Two intervals on the same side of zero with substantial overlap means ALS is a fine point estimate for that coefficient; intervals that disagree on sign or location point to the MCMC fit for that dataset. ## Longitudinal data `lame_als()` is the longitudinal counterpart. Its effects are static (pooled across time), as in a non-dynamic `lame()` fit. ```{r} set.seed(3) Yl <- replicate(5, { # intercept 0.4 + sender/receiver structure + noise (no dyadic # covariate is passed below, so only the intercept is recoverable) m <- 0.4 + outer(a, b, "+") + matrix(rnorm(n * n), n, n) diag(m) <- NA m }, simplify = FALSE) lf <- lame_als(Yl, R = 1, family = "normal", verbose = FALSE) coef(lf) ``` ## Summary | Task | Function | |------|----------| | Fast cross-sectional point estimate | `ame_als()` | | Fast longitudinal point estimate | `lame_als()` | | Bootstrap standard errors / intervals | `ame_als_bootstrap()` | | Calibrated posterior inference | `ame()` / `lame()` | The fast estimator returns point estimates and bootstrap or sandwich uncertainty. The MCMC estimator returns posterior summaries.