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