Package 'lame'

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)
fit <- ame(YX_nrm$Y, Xdyad = YX_nrm$X,
           nscan = 100, burn = 10, odens = 1, verbose = FALSE)

# Visualize sender effects
ab_plot(fit, effect = "sender")

# Visualize receiver effects without sorting
ab_plot(fit, effect = "receiver", sorted = FALSE)

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.

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

Description

A valued sociomatrix (Y) and 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: 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(addhealthc9)

Format

list

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

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

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 = 1,
  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,
  ...,
  model.name = NULL
)

Arguments

Details

This command provides posterior inference for parameters in AME models of cross-sectional relational data, assuming one of eight 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 an arc-sine prior on $(-1, 1)$ (density proportional to $(1-\rho^2)^{-1/2}$), 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; R > floor(n/3) is rarely identifiable and will issue a warning. 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 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: arc-sine prior on $(-1, 1)$, density $p(\rho) \propto (1-\rho^2)^{-1/2}$, 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 always uses Sab0 = diag(2). 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$. In Social Relations Model simulations the data-scaled default recovers small true variances with modest bias where a fixed diag(2) prior can nearly double them; the data-scaled prior is mildly more conservative in the opposite regime, when the true variance is large relative to the residual scale.

This default differs from the amen package, whose ame(family = "nrm") leaves Sab0 = diag(2) and eta0 = 4. To reproduce amen's additive-effects posterior, pass prior = list(Sab0 = diag(2), eta0 = 4). The choice affects only va, vb and (weakly) cab; ve and rho are unchanged.

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 in this model.

"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 AME model for count data: $Y \sim \mathrm{Poisson}(\exp(z))$ with $z \sim N(\eta, \sigma^2)$, a lognormal-mixed Poisson. The conditional mean given the latent $z$ is $\exp(z)$; the marginal mean is $\exp(\eta + \sigma^2/2)$, not $\exp(\eta)$. On sparse count matrices (only a few dozen nonzero cells) the multiplicative factors are weakly identified: UVPM can differ across seeds while the fit is otherwise stable. For such data compare seeds, check sum(fit$YPM) against the observed total, and prefer R = 0 unless the factors reproduce. The sampler includes a joint level move on the latent log-rates and the intercept (acceptance rate in fit$mh_counters$pois_shift_accept) so the predicted total settles quickly even when most cells are zero; the dyadic variance ve still mixes slowly on such data, so use long chains.

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. family = "ordinal" with symmetric = TRUE is supported via the dedicated sampler in R/rZ_ord_sym_fc.R, which uses the symmetric-doubled precision and mirrors 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 validator any(is.finite(Y - t(Y))) evaluates FALSE and the call proceeds: the sampler then treats the populated triangle as the symmetric data and mirrors it. This is usually intended, but if the upper / lower triangles were meant to differ, the model is silently fitting half the data. Audit anyNA(Y[upper.tri(Y)]) != anyNA(Y[lower.tri(Y)]) before calling if you are unsure.

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.

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.

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

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_bin)
fit <- ame(YX_bin$Y, Xdyad = YX_bin$X, burn = 10, nscan = 100, odens = 1,
           family = "binary", verbose = FALSE)
summary(fit)
# Note: you should run the Markov chain much longer in practice

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.

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

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

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

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.

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)

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 (requires T > 1); with few slices its intervals are necessarily coarse.

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. Replicates that error or return non-finite values are dropped and counted (n_valid / n_total). Standard errors are replicate column standard deviations; confidence intervals use the percentile method.

Usage

ame_als_bootstrap(
  object,
  R = 200,
  type = c("parametric", "block"),
  block_length = 1,
  seed = NULL,
  verbose = TRUE
)

boot_ame(
  object,
  R = 200,
  type = c("parametric", "block"),
  block_length = 1,
  seed = NULL,
  verbose = TRUE
)

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 is well-determined only when the multiplicative singular values are well separated; with near-equal singular values the per-column U/V standard errors reflect an unstable rotation, and the subspace — or, for symmetric models, the eigenvalues L — is the summary of record. For a symmetric model the eigenvalues L of the multiplicative term are aggregated and reported as the primary multiplicative-uncertainty summary.

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); with only a few time slices it is necessarily coarse, its intervals correspondingly imprecise and somewhat anti-conservative, so parametric is the default. A variance component fit with R > 0 can still carry a small residual parametric-bootstrap bias (the low-rank refit re-absorbs simulated noise); summary() flags any point estimate that falls outside its interval. A binary IRLS fit can carry a finite-sample (incidental-parameters) bias in the point estimator itself, which the bootstrap reproduces rather than removes. 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 n_valid, n_total, 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.

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

Y <- replicate(6, { 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)

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
)

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.

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)

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

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

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

Description

Run AME model with multiple parallel chains

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

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)

Description

Convert array to list.

Usage

array_to_list(arrayObj, actorList, sliceLabel)

Arguments

Value

array in list format

Author(s)

Shahryar Minhas

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

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(..., n_chains = N)) 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 n_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.

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

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.

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)
fit <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  autoplot(fit)
}

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.

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

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.

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

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.

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.

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 : https://polisci.la.psu.edu/people/xuc11/

Description

Combine multiple AME chains

Usage

combine_ame_chains(chain_list, diagnostics = TRUE)

Arguments

Value

A single ame object with combined chains

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

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

Description

Computes Gelman-Rubin 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

Value

A list with rhat, ess, param_names (over BETA and VC), n_chains and n_samples.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

Description

Replaces the O(T x p^2 x n^2) nested R loop for bipartite XtX and Xty computation with a single C++ call.

Usage

compute_XtX_Xty_bip_cpp(Xlist, resid, p)

Arguments

Value

List with XtX (p x p) and Xty (p vector)

Description

Eleven years of import and export data between 229 countries. The data use the SITC Rev. 1 commodity classification, aggregated at the first level (AG1).

Format

A list consisting of a socioarray Trade and a vector dollars2010 of inflation rates. The socioarray gives yearly trade volume (exports and imports) in dollars for 10 different commodity classes for eleven years between 229 countries. This gives a five-way array. The first index is the reporting country, so Trade[i,j,t,k,1] is what i reports for exports to j, but in general this is not the same as Trade[j,i,t,k,2], what j reports as importing from i.

Source

https://comtrade.un.org/, https://www.measuringworth.com/

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.

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

Usage

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

Arguments

Details

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

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.

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

Description

For each dynamic coefficient in a dynamic_beta fit, compares the posterior distribution of the maximum scaled first-difference $M = \max_t |\beta_t - \beta_{t-1}| / \sigma_\beta$ to a Monte-Carlo approximation of the same statistic under the AR(1) prior. The returned bf column is the tail-ratio score $\Pr(M^{post} > m^*) / 0.05$, where $m^*$ is the 95\ quantile of $M^{prior}$. It is kept under the historical column name for compatibility, but it is not a marginal-likelihood Bayes factor. Use it to surface posterior temporal jumps that the AR(1) prior cannot comfortably accommodate.

Usage

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

Arguments

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

Examples

data(YX_bin_list)
fit <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE)
detect_change_point(fit)

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

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

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

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

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

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
fit <- ame(Y_train, R = 0, family = "binary",
           burn = 15, nscan = 60, odens = 5, verbose = FALSE, plot = FALSE)
y_pred <- predict(fit, type = "response")
evaluate_heldout(Y, y_pred, mask, family = "binary")

Description

Returns the posterior mean of the network on the response scale (YPM). For binary models, these are predicted probabilities between 0 and 1. For normal models, these are predicted continuous values. For Poisson models, these are predicted counts.

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

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

Arguments

Value

A matrix or list of matrices of fitted values.

Description

Extract fitted values from LAME model

Usage

## S3 method for class 'lame'
fitted(object, ...)

Arguments

Value

List of fitted value matrices (one per time point)

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)

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 is computed analytically from the per-draw forecast variance. For discrete families ("binary", "poisson", "ordinal") the randomised PIT of Czado, Gneiting & Held (2009) is used. For rank families ("cbin", "frn") PIT is not implemented; the function returns NULL with an informational note.

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 (KS p-value), cover_95 (fraction of PIT values in [0.025, 0.975]).

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 <- YX_bin_list$Y[1:3]
Y_test <- YX_bin_list$Y[4]
X_train <- YX_bin_list$X[1:3]
fit <- lame(Y_train, Xdyad = X_train,
            family = "binary", R = 0,
            dynamic_beta = "dyad",
            nscan = 60, burn = 15, odens = 5, verbose = FALSE)
pit <- forecast_pit(fit, y_future = Y_test)
pit$ks_p
plot(pit)

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.

Description

Create design array for replicate data

Usage

get_design_rep(Y, Xdyad, Xrow, Xcol, actorSet, intercept, n, N, pr, pc, pd)

Arguments

Value

returns list of design array values necessary for ame_repL

Author(s)

Shahryar Minhas

Description

Rebuild the (n x n x T) or (nA x nB x T) EZ cube using a per-period beta vector. Mirrors get_EZ_cpp / get_EZ_bip_cpp but accepts a (T x p) beta matrix where each row is the beta for that period.

Usage

get_EZ_dynamic_beta_cpp(
  Xlist,
  beta_full_path,
  a_mat,
  b_mat,
  U_cube,
  V_cube,
  G,
  bipartite,
  symmetric
)

Arguments

Value

n_a x n_b x T cube of EZ values.

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
)

Arguments

Value

Fitted AME object

Author(s)

Shahryar Minhas

Description

Get fitted object from MCMC results

Usage

get_start_vals(
  start_vals,
  Y,
  family,
  xP,
  rvar,
  cvar,
  R,
  odmax = NULL,
  bip = FALSE
)

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.

Value

List of starting values for MCMC

Author(s)

Shahryar Minhas

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.

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)
fit <- lame(YX_bin_list$Y, YX_bin_list$X, family = "binary", R = 0,
            nscan = 100, burn = 20, odens = 5, verbose = FALSE)
glance(fit)

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_bin_list)
Y1 <- 1 * (YX_bin_list$Y[[1]] > 0); diag(Y1) <- NA
fit <- ame_als(Y = Y1, Xdyad = YX_bin_list$X[[1]],
               family = "binary", R = 1, verbose = FALSE)
glance(fit)

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

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)

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 matrix of GOF statistics with the same format as if gof=TRUE was used during model estimation. First row contains observed statistics, subsequent rows contain posterior predictive statistics.

Author(s)

Cassy Dorff, Shahryar Minhas, Tosin Salau

See Also

ame, lame, gof_plot