Getting Started with lame
Cassy Dorff, Shahryar Minhas, and Tosin Salau
2026-11-03
What Is lame?
Networks of trade, conflict, alliance, friendship, and sanction share
a common inferential difficulty: the tie between two actors is rarely
independent of the ties around it. A more central actor sends more,
attracts more, and shifts the incentives of every actor sharing a
neighbour. The lame package fits
Longitudinal Additive and
Multiplicative Effects models to this
kind of data. It estimates covariate effects on the directed tie scale
while accounting for sender heterogeneity, receiver heterogeneity,
reciprocity, and the residual relational structure that covariates alone
cannot capture.
The model decomposes each tie into three pieces. Every actor has a
tendency to send ties (\(a_i\)) and to
receive ties (\(b_j\)), reflecting how
active and how popular they are. Each actor also occupies a latent
position (\(u_i\) as a sender, \(v_j\) as a receiver) so that actors with
similar positions tend to form similar tie patterns. Covariates shift
the baseline probability of ties on top of those actor-level pieces.
When observations span multiple periods, any of these components can
evolve over time via an AR(1) process – short for “autoregressive of
order 1”, meaning each period’s value is a slightly noisy copy of the
previous period’s value (so positions drift smoothly rather than jumping
around).
Coming from a tidy edgelist? Use netify
If your data lives in a long-format edgelist, build a
netify object and pass it directly to lame()
or ame(). netify handles the matrix
construction, actor ordering, bipartite row/column sets, changing actor
composition, and dyadic or nodal covariates:
library(netify)
library(lame)
netlet <- netify::netify(
edges_long,
actor1 = "from",
actor2 = "to",
time = "year",
weight = "tie",
dyad_vars = "distance",
output_format = "longit_list",
missing_to_zero = FALSE
)
fit <- lame(netlet, family = "binary", R = 2, verbose = FALSE)
Use missing_to_zero = TRUE only when unlisted dyads are
real zeros; if an unlisted dyad means the relation was not observed,
keep it FALSE so those cells enter lame as
missing values. The underlying shape contract is still simple –
lame() consumes a named list of matrices plus optional
per-period covariate arrays – and the overview
article walks through an executed netify round-trip on real
data.
Coming from amen?
Both amen and lame export an
ame() function; if both packages are attached, a startup
message reminds you that unqualified ame(...) calls
dispatch to whichever was attached last, so write
lame::ame(...) to be explicit. Old amen
scripts largely run unchanged – the most common porting fix is
family = "nrm" becoming family = "normal" (the
short forms are accepted as aliases, with a message), and
print = is deprecated in favour of verbose =.
The full mapping of argument and output differences lives in the
Migration from amen sections of ?ame and
?lame.
A 5-Minute Example
Let’s fit a model to a small longitudinal binary network with one
dyadic covariate (a bilateral similarity score that drives tie
formation). The truth is intercept = -1.0 and the covariate
effect is beta = 0.7.
library(lame)
set.seed(6886)
# three periods of a 25-actor directed network driven by a dyadic covariate
n <- 25; T_periods <- 3
true_intercept <- -1.0
true_beta <- 0.7
# zero-padded names keep the input order under the alphabetical actor sort;
# matching dimnames on Y and every covariate array keep them aligned
actor_names <- sprintf("N%02d", seq_len(n))
X_list <- lapply(seq_len(T_periods), function(t) {
# [actor x actor x covariate] array; the slice name labels the coefficient
x <- matrix(rnorm(n * n), n, n)
array(
x,
dim = c(n, n, 1),
dimnames = list(actor_names, actor_names, "similarity")
)
})
Y_list <- lapply(seq_len(T_periods), function(t) {
# probit ties: tie probability pnorm(intercept + beta * x)
eta <- true_intercept + true_beta * X_list[[t]][, , 1]
Y <- matrix(rbinom(n * n, 1, pnorm(eta)), n, n)
diag(Y) <- NA # self-ties are undefined in a unipartite network
rownames(Y) <- colnames(Y) <- actor_names
Y
})
# longitudinal AME model with the similarity covariate
fit <- lame(
Y = Y_list,
Xdyad = X_list, # one similarity matrix per period
R = 2, # 2D latent space
family = "binary", # probit for 0/1 networks
burn = 20, # short burn-in for this example
nscan = 100, # short run for this example
odens = 5, # thinning
verbose = FALSE, # no progress output
plot = FALSE # no live diagnostic plots during sampling
)
summary(fit)
#>
#> ── Longitudinal AME model summary ──────────────────────────────────────────────
#>
#> ── Call
#> Y ~ dyad(similarity) + a[i] + b[j] + rho*e[ji] + U[i,1:2] %*% V[j,1:2], family = 'binary'
#>
#> ── Model
#> • Time periods: 3
#> • Family: binary
#> • Mode: unipartite
#> ! This is a STATIC fit pooled across 3 time periods: U, V, a, b are
#> time-invariant, and per-period predictions vary only through per-period
#> covariates.
#> ℹ For time-varying effects, refit with `dynamic_uv = TRUE` and/or `dynamic_ab =
#> TRUE`.
#>
#> ── Regression coefficients
#> Estimate StdError z_value p_value CI_lower CI_upper
#> intercept -1.068 0.041 -26.032 0 -1.146 -1.003 ***
#> similarity_dyad 0.732 0.049 15.062 0 0.645 0.825 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ℹ Stars are a visual hint from the posterior mean / SD only; for inference use
#> the credible intervals.
#>
#> ── Variance components
#> Estimate StdError
#> va 0.065 0.021
#> cab -0.005 0.014
#> vb 0.061 0.019
#> rho -0.113 0.077
#> ve 1.000 0.000
#> va = sender variance, cab = sender-receiver covariance, vb = receiver
#> variance, rho = dyadic correlation, ve = residual variance
The intercept lands near the true -1.0 and
similarity_dyad lands near the true 0.7, the
basic sanity check that the sampler is recovering known parameters. The
_dyad suffix on the coefficient name flags that the
covariate is a pair-level quantity; sender-only and receiver-only
covariates carry _row and _col suffixes
respectively. The remaining variance components describe residual
actor-level and dyad-level structure that the covariate alone does not
explain:
va: variance of the sender random effects, capturing
how much actors differ in their overall tendency to send ties.vb: variance of the receiver random effects, capturing
how much they differ in their tendency to receive.cab: covariance between sender and receiver effects
within an actor, indicating whether prolific senders are also frequent
receivers.rho: within-dyad residual correlation, capturing
reciprocity beyond what the additive effects and covariates
explain.
A note on MCMC settings. The three key parameters
are burn (iterations discarded as burn-in),
nscan (post-burn-in iterations kept for inference), and
odens (thinning: keep every odens-th sample to
reduce autocorrelation). With the settings above, we store
nscan / odens = 100 / 5 = 20 posterior samples (kept small
for vignette build time). For a final run, aim for at least 1000 stored
samples with adequate effective sample sizes
(e.g. burn = 1000, nscan = 25000, odens = 25).
What Can You Do with a Fitted Model?
lame objects work with all the standard R methods you’d
expect:
# regression coefficients (posterior means)
coef(fit)
#> intercept similarity_dyad
#> -1.0683169 0.7322389
# 95% credible intervals
confint(fit)
#> 2.5% 97.5%
#> intercept -1.1461086 -1.0032381
#> similarity_dyad 0.6447467 0.8253964
# one row per coefficient; works with or without broom loaded
tidy(fit)
#> # A tibble: 2 × 7
#> term estimate std.error statistic p.value conf.low conf.high
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 intercept -1.07 0.0410 -26.0 0 -1.15 -1.00
#> 2 similarity_dyad 0.732 0.0486 15.1 0 0.645 0.825
# one-row model summary; elpd_loo is NA unless a loo result is cached on the fit
glance(fit)
#> # A tibble: 1 × 13
#> nobs n_actors n_row_actors n_col_actors n_periods n_stored family mode R
#> <int> <int> <int> <int> <int> <int> <chr> <chr> <int>
#> 1 1800 25 NA NA 3 40 binary unip… 2
#> # ℹ 4 more variables: dynamic_uv <lgl>, dynamic_ab <lgl>, dynamic_beta <lgl>,
#> # elpd_loo <dbl>
# predicted tie probabilities: a list with one n x n matrix per period
Y_hat <- predict(fit, type = "response")
length(Y_hat) # 3, one matrix per period
#> [1] 3
dim(Y_hat[[1]]) # 25 x 25
#> [1] 25 25
cat(
"Predicted probability range:",
round(range(unlist(Y_hat), na.rm = TRUE), 3),
"\n"
)
#> Predicted probability range: 0 0.975
# residuals, observed minus predicted, in the same shape
resid_list <- residuals(fit)
cat("Residual SD:", round(sd(unlist(resid_list), na.rm = TRUE), 3), "\n")
#> Residual SD: 0.358
The two regression coefficients (intercept and
similarity_dyad) should land close to the simulation truth.
Predicted probabilities span a wide range because dyads with high
similarity have an intercept-plus-effect linear predictor much larger
than dyads with low similarity. Use confint(fit) for the
95% credible intervals on the coefficients.
type = "response" gives predictions back on the natural
scale of the outcome (probabilities for family = "binary"
and predicted values for "normal"). For MCMC
"poisson" fits it is the mean count of a cell,
exp(EZ + ve / 2) times the exposure of the period, at the
posterior means of the linear predictor and of the residual variance.
Predictive count quantiles and event probabilities come from
simulate(). type = "link" gives the underlying
linear predictor (the probit-scale value before the \(\Phi(\cdot)\) transform for binary), which
is what you want if you plan to compute marginal effects by hand. For
h-step-ahead forecasts on a longitudinal fit, use
predict(fit, h = K) (see the forecasting vignette).
Checking Your Model
Two diagnostics matter for any AME fit. The first is whether the MCMC
sampler explored the posterior adequately. The trace plots should mix
freely around a stable mean, without trends or stuck regions, and the
marginal densities should be unimodal.
trace_plot(fit, params = "beta")
For this fit the two regression parameters, the intercept and
similarity_dyad, mix around the simulation truths of \(-1.0\) and \(0.7\).
The second diagnostic is whether the model reproduces the structural
features of the observed network. Posterior-predictive goodness-of-fit
plots simulate networks from the fitted model and compare their
structural statistics to those of the observed data. The observed series
is a solid black line with points; the posterior-predictive median is a
dashed grey line and the 95% credible interval is the grey band.
Panel titles use human-readable names (Sender degree heterogeneity,
Transitivity, and so on) rather than the underlying column codes such as
sd.rowmean or trans.dep; the internal codes
are what fit$GOF returns directly. Because we generated
this data purely from eta = -1 + 0.7 * X (no planted degree
heterogeneity, reciprocity, or clustering), the model reproduces almost
every structural statistic: across the five statistics and three
periods, fourteen of the fifteen observed values sit inside their 95%
posterior-predictive bands. The lone marginal excursion (transitivity at
one period) sits just past the band edge on a statistic whose values are
all near zero, i.e. Monte-Carlo noise rather than a structural misfit.
On real friendship or trade data the picture is different: at least one
panel, typically transitivity or sender degree heterogeneity, drifts
clearly outside the band. That misfit is informative: it tells the
analyst which higher-order features AME absorbs through its actor and
latent structure and which remain unexplained. It is a substantive
finding about the data, not grounds to discard the fit.
Visualizing Network Structure
The latent space is where you look for actors playing similar roles
in the network:
And the additive-effects plot is where you look for unusually active
or popular actors:
ab_plot(fit, effect = "sender")
A note on this example: the simulation planted no latent structure
and no sender/receiver heterogeneity (ties depend only on the similarity
covariate), and the near-zero va and vb in the
summary above confirm the model found essentially none. The apparent
clusters in the latent-space plot and the longer stems here are
therefore sampling noise: these two figures show the displays and carry
no substantive finding. On real data these are the plots where the
substantive story lives; see the overview
article for a friendship network in which the latent space captures
genuine structure.
Cross-Sectional Models
If you have a single network (not a time series), use
ame() instead of lame():
fit_cs <- ame(
Y = Y_list[[1]], # just one time period
Xdyad = X_list[[1]], # n x n x 1 similarity array
R = 2,
family = "binary",
burn = 20,
nscan = 100,
odens = 5,
verbose = FALSE
)
coef(fit_cs)
#> intercept similarity_dyad
#> -1.2326924 0.7697791
The output and methods are the same. The difference is that
lame() pools information across time periods, giving you
more precise estimates when the structure is stable.
Dynamic Effects
When you believe the network structure is changing over time
(alliances shifting, friendships evolving), you can let the latent
positions and additive effects drift through time. The drift is modelled
as an AR(1) process: each period’s value is a noisy
copy of the previous period’s value, with a persistence parameter \(\rho \in (-1, 1)\) controlling how strongly
past predicts present (\(\rho\) near 1
= slow drift, \(\rho\) near 0 = each
period almost independent).
fit_dyn <- lame(
Y = Y_list,
Xdyad = X_list,
R = 2,
dynamic_ab = TRUE, # time-varying sociality/popularity
dynamic_uv = TRUE, # time-varying latent positions
family = "binary",
burn = 20,
nscan = 100,
odens = 5,
verbose = FALSE,
plot = FALSE
)
summary(fit_dyn)
#>
#> ── Longitudinal AME model summary ──────────────────────────────────────────────
#>
#> ── Call
#> Y ~ dyad(similarity) + a[i] + b[j] + rho*e[ji] + U[i,1:2] %*% V[j,1:2], family = 'binary'
#>
#> ── Model
#> • Time periods: 3
#> • Family: binary
#> • Mode: unipartite
#> • Dynamic latent positions: enabled (kind = ar1, rho_uv = 0.723)
#> rho_uv: posterior sd 0.055 is 0.71 of the prior sd 0.077: the data inform it
#> only weakly.
#> • Dynamic additive effects: enabled (rho_ab = 0.622)
#> rho_ab: posterior sd 0.111 is 0.88 of the prior sd 0.126: the data barely
#> inform it.
#>
#> ── Regression coefficients
#> Estimate StdError z_value p_value CI_lower CI_upper
#> intercept -1.193 0.05 -23.781 0 -1.271 -1.103 ***
#> similarity_dyad 0.805 0.051 15.938 0 0.728 0.888 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ℹ Stars are a visual hint from the posterior mean / SD only; for inference use
#> the credible intervals.
#>
#> ── Variance components
#> Estimate StdError
#> va 0.069 0.017
#> cab -0.001 0.016
#> vb 0.064 0.021
#> rho -0.126 0.130
#> ve 1.000 0.000
#> va = sender variance, cab = sender-receiver covariance, vb = receiver
#> variance, rho = dyadic correlation, ve = residual variance
#> Under `dynamic_ab`, va and vb are the variances of the actor levels; rho_ab
#> and sigma_ab describe the per-period deviations around them.
The model estimates how persistent the latent positions are (\(\rho_{uv}\)) and how persistent the
additive effects are (\(\rho_{ab}\)).
Values near 1 indicate slow evolution; values near 0 indicate
near-independent positions across periods. The simulated data here use a
single latent structure across all three periods, varying only the
dyadic similarity covariate, so the persistence estimates reflect a mix
of the prior and whatever stable-position signal the chain extracts from
three time points. The dynamic
effects vignette provides a longer-panel example in which temporal
evolution is genuinely present in the simulation and the persistence
parameters are informative.
Bipartite Networks
For two-mode networks (students and courses, countries and treaties),
pass a rectangular matrix and set mode = "bipartite". Below
we simulate a 15 row × 10 column network where a dyadic similarity score
again drives the ties.
set.seed(42) # reproducible simulation
nA <- 15; nB <- 10
row_names <- sprintf("R%02d", seq_len(nA))
col_names <- sprintf("C%02d", seq_len(nB))
X_bip <- lapply(1:3, function(t) {
x <- matrix(rnorm(nA * nB), nA, nB)
array(
x,
dim = c(nA, nB, 1),
dimnames = list(row_names, col_names, "similarity")
)
})
Y_bip <- lapply(1:3, function(t) {
eta <- -0.8 + 0.6 * X_bip[[t]][, , 1]
Y <- matrix(rbinom(nA * nB, 1, pnorm(eta)), nA, nB)
rownames(Y) <- row_names; colnames(Y) <- col_names
Y
})
fit_bip <- lame(
Y = Y_bip,
Xdyad = X_bip,
mode = "bipartite",
R = 2,
family = "binary",
burn = 20,
nscan = 100,
odens = 5,
verbose = FALSE,
plot = FALSE
)
summary(fit_bip)
#>
#> ── Longitudinal AME model summary ──────────────────────────────────────────────
#>
#> ── Call
#> Y ~ dyad(similarity) + a[i] + b[j] + U[i,1:2] %*% G %*% V[j,1:2]', family = 'binary'
#>
#> ── Model
#> • Time periods: 3
#> • Family: binary
#> • Mode: bipartite
#> ! This is a STATIC fit pooled across 3 time periods: U, V, a, b are
#> time-invariant, and per-period predictions vary only through per-period
#> covariates.
#> ℹ For time-varying effects, refit with `dynamic_uv = TRUE` and/or `dynamic_ab =
#> TRUE`.
#>
#> ── Regression coefficients
#> Estimate StdError z_value p_value CI_lower CI_upper
#> intercept -0.946 0.084 -11.277 0 -1.102 -0.763 ***
#> similarity_dyad 0.644 0.101 6.373 0 0.494 0.835 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ℹ Stars are a visual hint from the posterior mean / SD only; for inference use
#> the credible intervals.
#>
#> ── Variance components
#> Estimate StdError
#> va 0.064 0.026
#> cab 0.000 0.000
#> vb 0.076 0.030
#> rho 0.000 0.000
#> ve 1.000 0.000
#> va = sender variance, cab = sender-receiver covariance, vb = receiver
#> variance, rho = dyadic correlation, ve = residual variance
#> ℹ Bipartite model: rho and cab are fixed to 0.
Notice that the bipartite summary reports cab = 0.000
and rho = 0.000. These two parameters describe
relationships between the sender and receiver of a tie, which
only makes sense when the same set of actors can play both roles. In a
bipartite network the rows and columns are different kinds of entities
(rows = students, columns = courses; rows = donors, columns =
candidates), so a row actor is never also a column actor. There is no
within-actor sender / receiver covariance and no reciprocity from \(i \to j\) to \(j
\to i\) to estimate, and the package fixes both at zero. The
similarity_dyad coefficient has a posterior mean of 0.64
and a 95% interval of [0.494, 0.835], compared with a simulation truth
of 0.6. This short example demonstrates extraction of posterior
summaries. The bipartite
vignette walks through the full workflow.
Supported Data Types
"normal" |
Continuous |
Trade volumes, survey ratings |
"binary" |
0/1 |
Friendships, alliances, sanctions |
"ordinal" |
Ordered categories |
Conflict intensity (none/threat/action) |
"poisson" |
Counts |
Number of co-sponsored bills |
"cbin" |
Censored binary |
Friendships with nomination limits |
"frn" |
Fixed rank nomination |
“Name your top 5 friends” |
Power-User Features
A few features that show up in advanced workflows but rarely in a
first model. Each gets a one-line pointer here rather than a full
demo:
- Choosing the rank. There is no cap on
R below the network size: ranks above
floor(n/3) fit without numerical failure, but the chain
mixes more slowly and the fitted values lose accuracy gradually once the
rank passes the structure in the data. Compare
R = 0, 1, 2, 3 by their fit; R = 2 or
R = 3 is the usual default.
- Checkpoint / resume for long runs. Pass
max_seconds and checkpoint_path when a run may
hit a wall-clock limit, then continue with resume_from; see
?lame_resume for the resume-cycle semantics.
- K-panel joint posterior.
lame_multi()
fits K parallel networks and pools the per-panel beta posteriors into
one precision-weighted shared posterior; see
?lame_multi.
- Memory-conscious
loo().
save_log_lik = "chunked" streams the log-score matrix to
disk instead of RAM. loo(fit) reads it when the saved
observed likelihood-factor contract passes; marginal and latent
diagnostic scores remain available through read_log_lik().
See the dynamic
effects vignette.
- Multi-chain diagnostics.
lame() runs
two chains by default (chains = 2) and stores their \(\hat R\) in fit$rhat;
rhat_dynamic_beta() gives a multivariate between-chain
\(\hat R\) for coefficient paths from a
list of fits run under different seeds; see the dynamic
effects vignette.
- Held-out predictive scoring. Mask dyads to
NA, refit (the sampler imputes them), and score the masked
cells with evaluate_heldout(). AUROC / PR-AUC require
precrec (or AUROC alone via pROC); both are
optional and the helper degrades gracefully when neither is installed.
See the cross-sectional vignette for a
worked example.
The one code pattern worth spelling out is checkpoint / resume:
ck <- tempfile(fileext = ".rds")
fit1 <- lame(
Y = Y_list,
Xdyad = X_list,
R = 2,
family = "binary",
nscan = 5000,
burn = 200,
odens = 25,
chains = 1, # a checkpoint is one chain
max_seconds = 30, # stop after 30 s wall clock
checkpoint_path = ck, # write state here
verbose = FALSE
)
if (isTRUE(fit1$terminated_early)) {
# 2,000 more sweeps, keeping every 25th
fit2 <- lame(resume_from = ck, nscan = 2000)
}
References
Hoff, PD (2021). Additive and Multiplicative Effects Network Models.
Statistical Science, 36, 34–50.
Minhas, S., Dorff, C., Gallop, M. B., Foster, M., Liu, H., Tellez,
J., & Ward, M. D. (2022). Taking dyads seriously. Political
Science Research and Methods, 10(4), 703–721.