Getting Started with lame

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 vignette walks through an executed netify round-trip on real data.

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)

# simulate 3 time periods of a 25-node directed network with a
# dyadic similarity covariate driving the ties.
n <- 25; T_periods <- 3
true_intercept <- -1.0
true_beta      <- 0.7

# zero-padded names ("N01" ... "N25") sort the same way alphabetically as
# positionally, so `lame()`'s internal alphabetic actor sort leaves the
# row order unchanged. with non-padded names ("N1", "N10", "N2", ...) the
# stored output comes back in sorted order, not your input order -- the fit
# is still correct as long as Y and the X arrays share the same names.
# explicit dimnames on Y and every covariate array are the safe way to
# guarantee that alignment.
actor_names <- sprintf("N%02d", seq_len(n))

X_list <- lapply(seq_len(T_periods), function(t) {
    # 3-D array [actor x actor x covariate]; the third-dim name becomes
    # the coefficient label downstream
    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) {
    # build the linear predictor eta = intercept + beta * X (the same
    # arithmetic a logistic / probit regression would do). then push it
    # through pnorm() -- the standard-normal CDF -- to turn it into a
    # tie probability in [0, 1]. that CDF link is what "probit" means,
    # and it is the link `family = "binary"` uses inside `lame()`. finally,
    # draw a 0/1 tie from a Bernoulli with that probability.
    #
    # the `[, , 1]` peels the first (and only) covariate slice off the
    # 3-D array, leaving an n x n matrix. third-dim index = covariate.
    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
})

# fit a longitudinal AME model with the dyadic 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,              # compact burn-in for this example
    nscan = 100,            # compact post-burn-in run for the vignette
    odens = 5,              # thinning
    verbose = FALSE,        # suppress the progress bar / iteration log
    plot = FALSE            # don't pop up live MCMC diagnostic plots during sampling
                            # `lame()` draws live diagnostics when plot = TRUE;
                            #  `ame()` accepts plot = for signature parity but ignores it
                            #  -- the single-period sampler has no live plotting
)

summary(fit)
#> 
#> === Longitudinal AME Model Summary ===
#> 
#> Call:
#> [1] "Y ~ dyad(similarity) + a[i] + b[j] + rho*e[ji] + U[i,1:2] %*% V[j,1:2], family = 'binary'"
#> 
#> Time periods: 3 
#> Family: binary 
#> Mode: unipartite 
#> 
#> Note: STATIC fit pooled across 3 time periods --
#>   U, V, a, b are time-invariant; 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.112    0.034 -32.307       0   -1.169    -1.06 ***
#> similarity_dyad    0.742    0.051   14.51       0    0.646    0.813 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Note: stars are a visual hint from posterior mean / SD only; for inference use the credible intervals.
#> 
#> Variance components:
#> -------------------
#>     Estimate StdError
#> va     0.055    0.012
#> cab    0.000    0.018
#> vb     0.065    0.034
#> rho   -0.052    0.099
#> ve     1.000    0.000
#>   (va = sender, cab = sender-receiver covariance, vb = receiver,
#>    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.1116591       0.7415237

# 95% credible intervals
confint(fit)
#>                      2.5%      97.5%
#> intercept       -1.168612 -1.0596470
#> similarity_dyad  0.646386  0.8126144

# broom-style one-row-per-coefficient frame; ships with lame so it works
# without broom installed and dispatches through broom::tidy(fit) when
# broom is loaded. glance(fit) gives the one-row model summary that
# modelsummary uses for its lower panel. See the overview vignette for
# the full modelsummary / tidybayes / autoplot round-trip.
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.11     0.0344     -32.3       0   -1.17     -1.06 
#> 2 similarity_dyad    0.742    0.0511      14.5       0    0.646     0.813

# n_row_actors / n_col_actors are bipartite-only and are NA on a unipartite
# fit like this one. elpd_loo is NA unless the model was fit with
# save_log_lik = TRUE and the loo result cached via fit$loo <- loo(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       20 binary unip…     2
#> # ℹ 4 more variables: dynamic_uv <lgl>, dynamic_ab <lgl>, dynamic_beta <lgl>,
#> #   elpd_loo <dbl>

# predicted probabilities for every dyad at every time point.
# for a `lame()` fit (panel data), `predict()` returns a *list of length T*,
# where T is the number of time periods. Each element is an n x n matrix of
# posterior-mean predicted tie probabilities (between 0 and 1, since
# family = "binary"). For a single-period `ame()` fit, predict() returns
# the n x n matrix directly, not wrapped in a list.
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 1

# residuals: same list-of-matrices shape, observed minus predicted
resid_list <- residuals(fit)
cat("Residual SD:", round(sd(unlist(resid_list), na.rm = TRUE), 3), "\n")
#> Residual SD: 0.353

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", expected counts for "poisson", expected values for "normal"). 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")

MCMC trace and density plots for the regression coefficients (intercept and similarity_dyad); well-mixed traces and unimodal densities indicate adequate convergence.

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 dual-encoded as an Okabe-Ito orange (#D55E00) solid line with points; the posterior-predictive median is grey-dashed and the 95% credible interval is the grey ribbon. The colour-plus-linetype encoding stays legible in greyscale and for colour-blind readers.

gof_plot(fit)

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:

uv_plot(fit)

Circular-layout latent-space plot: each actor's sender position is a triangle and receiver position a circle on concentric rings; actors placed near each other share similar tie patterns.

And the additive-effects plot is where you look for unusually active or popular actors:

ab_plot(fit, effect = "sender")

Lollipop plot of sender effects, sorted by posterior mean; each actor is a stem from zero to its point estimate, with positive values marking unusually active senders. In this null example all stems are short and reflect sampling noise, not planted heterogeneity.

A caveat for this particular 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 – treat these two figures as a tour of the displays, not a substantive finding. On real data these are the plots where the substantive story lives; see the overview vignette 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]],     # 3-D array [n, n, 1] with "similarity" slice name
    R = 2,
    family = "binary",
    burn = 20,
    nscan = 100,
    odens = 5,
    verbose = FALSE
)

coef(fit_cs)
#>       intercept similarity_dyad 
#>       -1.231792        0.776326

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:
#> [1] "Y ~ dyad(similarity) + a[i] + b[j] + rho*e[ji] + U[i,1:2] %*% V[j,1:2], family = 'binary'"
#> 
#> Time periods: 3 
#> Family: binary 
#> Mode: unipartite 
#> Dynamic latent positions: enabled (rho_uv = 0.886 )
#> Dynamic additive effects: enabled (rho_ab = 0.698 )
#> 
#> Regression coefficients:
#> ------------------------
#>                 Estimate StdError z_value p_value CI_lower CI_upper    
#> intercept         -1.106    0.056 -19.673       0   -1.177       -1 ***
#> similarity_dyad    0.759    0.051  14.935       0    0.659    0.843 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Note: stars are a visual hint from posterior mean / SD only; for inference use the credible intervals.
#> 
#> Variance components:
#> -------------------
#>     Estimate StdError
#> va     0.071    0.027
#> cab   -0.003    0.014
#> vb     0.072    0.025
#> rho   -0.089    0.105
#> ve     1.000    0.000
#>   (va = sender, cab = sender-receiver covariance, vb = receiver,
#>    rho = dyadic correlation, ve = residual variance)

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)          # seed the data simulation so the recovery is reproducible
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:
#> [1] "Y ~ dyad(similarity) + a[i] + b[j] + U[i,1:2] %*% G %*% V[j,1:2]', family = 'binary'"
#> 
#> Time periods: 3 
#> Family: binary 
#> Mode: bipartite 
#> 
#> Note: STATIC fit pooled across 3 time periods --
#>   U, V, a, b are time-invariant; 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.911    0.073 -12.555       0    -1.05   -0.796 ***
#> similarity_dyad    0.604     0.06   10.01       0    0.503    0.709 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Note: stars are a visual hint from posterior mean / SD only; for inference use the credible intervals.
#> 
#> Variance components:
#> -------------------
#>     Estimate StdError
#> va     0.394    0.124
#> cab    0.000    0.000
#> vb     0.519    0.191
#> rho    0.000    0.000
#> ve     1.000    0.000
#>   (va = sender, cab = sender-receiver covariance, vb = receiver,
#>    rho = dyadic correlation, ve = residual variance)
#>   Note: bipartite model (rho fixed to 0, cab 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 should land near the simulation truth of 0.6 (with this seed the point estimate is about 0.6, a little above the truth, and the 95% interval comfortably covers 0.6). The bipartite vignette walks through the full workflow.

Supported Data Types

Family	Data	Example
`"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:

The R > n/3 warning. Both ame() and lame() warn when the latent-space rank exceeds floor(n/3) (or floor(min(nA, nB)/3) in bipartite mode); past that, the multiplicative effects absorb structure that belongs to the additive effects. Advisory, not an error – R = 2 or R = 3 is the right 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-likelihood matrix to disk instead of RAM, and loo(fit) reads it back transparently; see the dynamic effects vignette.
Multi-chain diagnostics. lame_parallel(..., n_chains = 4) plus rhat_dynamic_beta() give a between-chain $\hat R$ ; 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,
    max_seconds = 30,          # stop after 30 s wall clock
    checkpoint_path = ck,      # write state here
    verbose = FALSE
)
if (isTRUE(fit1$terminated_early)) {
    # `nscan` on the resume call = additional stored draws
    fit2 <- lame(resume_from = ck, nscan = 2000)
}

Where to Go Next

I want to…	Read this
Understand the full workflow with real data	lame overview
Learn about cross-sectional models in depth	Your first AME model
Model two-mode networks	Bipartite networks
Let the network structure evolve over time	Dynamic effects

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.

Cassy Dorff, Shahryar Minhas, and Tosin Salau

2026-07-25