---
title: "lame Overview"
author: "Cassy Dorff, Shahryar Minhas, and Tosin Salau"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{lame Overview}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
	collapse = TRUE,
	comment = "#>",
	fig.width = 6,
	fig.height = 4,
	dpi = 96,
	out.width = "100%",
	fig.align = "center",
	dev = "png"
)
```

## Package Overview

The `lame` package provides tools for fitting **L**ongitudinal **A**dditive and **M**ultiplicative **E**ffects models to network data observed over time. If you study relationships between actors (countries trading with each other, legislators co-sponsoring bills, students forming friendships across semesters) and you have repeated observations of those relationships, `lame` is designed for you.

The modeling approach builds on the Additive and Multiplicative Effects (AME) framework developed by Peter Hoff, whose [`amen`](https://pdhoff.github.io/amen/) package provides the foundational cross-sectional implementation. `lame` extends that framework with:

- **Longitudinal data**: panel networks observed across multiple time periods, potentially with actors entering and leaving the sample.
- **Dynamic effects**: latent positions and baseline activity levels can evolve via AR(1) processes (`dynamic_uv` and `dynamic_ab`).
- **C++ acceleration**: core sampling routines in Rcpp/RcppArmadillo.
- **Bipartite networks**: two-mode networks (e.g., countries and treaties) with separate latent spaces for row and column nodes.
- **ggplot2-based diagnostics**: `trace_plot`, `gof_plot`, `ab_plot`, and `uv_plot` return standard ggplot2 figures you can theme and extend.
- **Standard S3 methods**: `coef()`, `confint()`, `fitted()`, `residuals()`, `predict()`, and `simulate()` work as you would expect from any R model object.

The package works with [`netify`](https://netify-dev.github.io/netify/) for data preparation.

### Migrating 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`.

## Application: Dutch College Friendships

To see how `lame` works in practice, we analyze the **Dutch college friendship network** (van de Bunt, van Duijn, and Snijders 1999), a panel of 32 students who arrived as strangers and rated each other's friendship at seven time points. (We drop the cold-start first wave, when the students barely knew each other and almost no ties exist; six waves remain.) The substantive question is whether students sort by gender, smoking status, and academic program, and whether residual clustering remains once those covariates are absorbed.

The dataset includes:

- **Y**: directed friendship ratings on a -1 to 4 scale (-1 a negative/troubled relationship, 0 no or uncertain tie, 1-4 increasing friendship). We binarize any positive rating to a tie (`(Y > 0) * 1`), so the rare -1 ratings fold in with 0 as "no tie".
- **X**: three node-level attributes: `male`, `smoker`, `program`.
- We construct three dyadic **homophily indicators** from the node attributes: `same_male`, `same_smoker`, `same_program`.

```{r setup}
library(lame)
library(ggplot2)
set.seed(6886)

data("dutchcollege")

n <- nrow(dutchcollege$Y)
T_all <- dim(dutchcollege$Y)[3]
actor_names <- sprintf("S%02d", seq_len(n))

# binarize and drop the cold-start wave (t = 1 has almost no ties)
Y <- lapply(2:T_all, function(t) {
	Yt <- (dutchcollege$Y[, , t] > 0) * 1
	diag(Yt) <- NA
	rownames(Yt) <- colnames(Yt) <- actor_names
	Yt
})
names(Y) <- paste0("t", 2:T_all)

# nodal covariates passed identically as sender and receiver
X_node <- dutchcollege$X
rownames(X_node) <- actor_names
Xrow <- lapply(seq_along(Y), function(t) X_node)
Xcol <- Xrow

# dyadic homophily indicators
same_male    <- outer(X_node[, "male"],    X_node[, "male"],    "==") * 1
same_smoker  <- outer(X_node[, "smoker"],  X_node[, "smoker"],  "==") * 1
same_program <- outer(X_node[, "program"], X_node[, "program"], "==") * 1
Xdyad_one <- array(0, dim = c(n, n, 3),
                   dimnames = list(actor_names, actor_names,
                                   c("same_male", "same_smoker", "same_program")))
Xdyad_one[, , 1] <- same_male
Xdyad_one[, , 2] <- same_smoker
Xdyad_one[, , 3] <- same_program
Xdyad <- lapply(seq_along(Y), function(t) Xdyad_one)
```

The friendship network's average density rises across the panel as students settle in, though the trajectory is not monotone (density bounces between roughly 0.40 and 0.63 across the six waves):

```{r}
sapply(Y, function(y) round(mean(y, na.rm = TRUE), 2))
```

### Fitting the Model

We fit a binary probit AME model with sender random effects (`rvar`), receiver random effects (`cvar`), dyadic correlation (`dcor`), and a one-dimensional latent space (`R = 1`). The random effects capture sender heterogeneity (some students simply nominate more friends) and receiver heterogeneity (some students are nominated more often), while the latent space picks up residual clustering: students who befriend similar people for reasons the observed covariates do not capture. (Why `R = 1` and not 2? See the Latent Space section below.)

```{r, message=FALSE, warning=FALSE}
fit <- lame(
	Y = Y,
	Xdyad = Xdyad,           # dyadic homophily indicators
	Xrow = Xrow,             # sender covariates
	Xcol = Xcol,             # receiver covariates
	family = "binary",       # binary probit model
	rvar = TRUE,             # sender random effects
	cvar = TRUE,             # receiver random effects
	dcor = TRUE,             # dyadic correlation (reciprocity)
	R = 1,                   # 1-D multiplicative latent space
	symmetric = FALSE,       # friendships are directed
	burn = 30,               # short burn-in for this worked example
	nscan = 150,             # compact post-burn-in run for the vignette
	odens = 10,              # thinning -> 15 stored draws
	posterior_opts = list(save_UV = TRUE),  # keep U/V draws: gives
	                                        # latent_positions() real
	                                        # posterior SDs below
	save_log_lik = TRUE,     # pointwise log-lik: feeds loo() at the end
	verbose = FALSE,
	plot = FALSE
)
```

### Interpreting the Results

```{r}
summary(fit)
```

A lot lands at once, so let's walk through it.

**Intercept.** The probit intercept is the latent linear predictor when *every covariate is zero*. The covariates here are uncentered (`program` is coded 2-4 and is never zero), so that corner sits well outside the observed data; the familiar `qnorm(density)` rule of thumb instead describes the predictor at the *average* covariate profile. The two reconcile once you add back the mean covariate contribution:

```{r intercept-decomp}
bhat        <- colMeans(fit$BETA)
slice_means <- apply(fit$X[[1]], 3, mean, na.rm = TRUE)   # mean of each design column
centroid_lp <- sum(bhat * slice_means[names(bhat)])       # predictor at the average profile
c(intercept     = unname(bhat["intercept"]),
  centroid_lp   = centroid_lp,
  qnorm_density = qnorm(mean(unlist(Y), na.rm = TRUE)))
```

The intercept itself is strongly negative, but the predictor at the average covariate profile (`centroid_lp`) sits close to `qnorm(density)` -- the uncentered covariates account for nearly the entire gap. The intercept absorbs the mean covariate contribution; it is not a baseline probability. Center the covariates if you want it to read as one, and read predicted probabilities off `predict(fit, type = "response")`.

**Homophily.** The three dyadic coefficients (`same_male`, `same_smoker`, `same_program`) measure whether students of the same type are more likely to be friends, holding the sender, receiver, and latent positions fixed. Positive values mean homophily; credible intervals that exclude zero indicate that the direction of the association is well-identified. Program homophily tends to be the strongest of the three, consistent with students bonding around shared coursework.

**Nodal covariates.** `male_row` / `male_col` measure whether male students are unusually active senders or popular receivers, and similarly for `smoker` and `program`. With only 32 students, the credible intervals on individual nodal effects are wide; treat the sign as suggestive, the magnitude with caution.

**Variance components.** `va` (sender variance) and `vb` (receiver variance) quantify how much students differ in sociability and popularity beyond what the covariates explain. The dyadic correlation `rho` captures reciprocity; in friendship networks `rho > 0` is the rule rather than the exception.

### Checking Convergence

The model is estimated via MCMC, so convergence checks come first. `trace_plot` shows the sampled values over iterations (top panels) and the corresponding posterior densities (bottom panels) for each parameter.

```{r overview-trace, fig.width=8, fig.height=8, dpi=100, dev="png", fig.alt="MCMC trace plots (top) and marginal posterior density plots (bottom) for each regression coefficient of the Dutch college fit; well-mixed fuzzy-caterpillar traces around a stable mean and smooth unimodal densities are the signatures of a converged chain."}
trace_plot(fit, params = "beta")
```

Trace plots should bounce around a stable mean without long-term trends or sticky regions, and density plots should be smooth and unimodal. Numerically, common MCMC checks are **split-$\hat R$ < 1.01** for monitored parameters and **bulk / tail ESS $\ge$ 400 per chain**.

With `r nrow(fit$BETA)` stored samples, this run cannot meet those thresholds: the number of stored draws caps the ESS, and several $\hat R$ values in the `summarise_draws` table below land well above 1.01. Read the table as a demonstration of how to check these diagnostics, not as evidence of convergence -- for a real analysis, lengthen the chain (`burn` in the hundreds, `nscan` in the thousands) until every monitored parameter clears both thresholds. Which parameter mixes worst is seed-dependent in a short chain, so treat this fit as a compact worked example rather than a final diagnostic run.

### Bayesian-ecosystem diagnostics: `posterior::as_draws()`

`lame` registers `as_draws()` methods so any fit drops straight into the Stan-era diagnostic ecosystem (`posterior`, `bayesplot`, `tidybayes`): `posterior::summarise_draws()` gives per-parameter posterior summaries, **split-$\hat R$**, and **bulk / tail ESS** in one call.

```{r as-draws-demo, eval = requireNamespace("posterior", quietly = TRUE)}
library(posterior)
draws <- posterior::as_draws(fit)              # draws_array [iter, chain, var]
posterior::summarise_draws(draws)              # rhat, ess_bulk, ess_tail per param
```

A single-chain fit like this one has `chain = 1`, so $\hat R$ reduces to within-chain split-$\hat R$; for $\hat R$ across independently initialised chains, see the multi-chain section in the dynamic-effects vignette.

### Goodness of Fit

The `gof_plot` function compares observed network statistics to their posterior predictive distributions: can the model reproduce emergent structural features like degree heterogeneity, transitivity, and reciprocity, not just the dyad-level relationships?

```{r overview-gof, fig.width=8, fig.height=6, dpi=100, dev="png", fig.alt="Longitudinal goodness-of-fit panels for the Dutch college fit: each facet plots one network statistic (Sender Degree Heterogeneity, Receiver Degree Heterogeneity, Dyadic Dependence, Triadic Dependence, Transitivity) across the six time periods. The observed series is a solid orange line with filled orange points; the posterior-predictive median is a dashed dark line; the 95 percent credible interval is a grey ribbon. The three elements differ on both colour and linetype, so the comparison survives grayscale printing and colour-blind viewing. Observed series falling outside the band flag structural features the model is missing."}
gof_plot(fit)
```

For a longitudinal model, `gof_plot` shows these statistics across time, encoding three series per panel with colour and linetype redundantly: the **observed series is a solid orange line with filled orange points**, the **posterior-predictive median is a dashed dark line**, and the **95% credible interval is a grey ribbon** (width set by `credible.level`). Observed values inside the ribbon are features the model reproduces; values outside flag structure it is missing -- itself a substantive finding, not a failure to report.

Reading the figure: the fit comes closest on the cyclic-triad structure (Triadic Dependence sits inside the band at four of the six periods) but under-predicts sender degree heterogeneity (observed above the band at every period) and transitivity (above the band at five of six periods) -- some students nominate many more friends than the additive $a_i$ plus latent $u_i$ can absorb, and triangle closure is systematically too low. Receiver heterogeneity drifts in and out of the band across the panel, while Dyadic Dependence (reciprocity) swings from one side of it to the other without settling inside.

This is the standard signature of a latent-space model on a friendship network. AME assumes conditional dyadic independence, so triangle closure must be soaked up by clustering of the latent positions; ERGMs model it directly via terms like `gwesp`. The two classes are complementary: ERGM for triangle-closure questions, AME for stable, non-degenerate estimates with explicit actor heterogeneity and a coherent posterior to forecast or simulate from.

#### Posterior-predictive temporal-trend test: `gof_temporal()`

The figure read above is qualitative. `gof_temporal()` gives a scalar check on the temporal *trend* in a chosen statistic: it fits an OLS slope to the observed per-period series, recomputes that slope in posterior-predictive replicate panels drawn via `simulate(fit)`, and returns a two-sided posterior-predictive p-value -- near 0 when the observed slope is unreachable under the model, well away from 0 when it is central in the predictive distribution.

```{r overview-gof-temporal}
# density rises across the panel on net (the network gets denser as
# students settle in, though the trajectory is not monotone), so it is the
# natural temporal-trend target.
gt_density <- gof_temporal(fit, stat = "density", n_rep = 100, seed = 6886)
gt_density   # prints stat, observed slope, n_rep, and the p_pp
```

A `p_pp` of 0 says the static fit cannot reproduce the empirical net densification slope (about +0.025 per period in our run) -- the expected diagnostic for a model whose parameters are constant across periods, and consistent with the per-period mass shift in `gof_plot()`. `dynamic_beta` or time-varying covariates are the remedy; see the [Dynamic Effects](dynamic_effects.html) vignette.

### Latent Space

The multiplicative effects capture association patterns beyond what the
covariates and additive effects explain. We fit `R = 1`, so each student's
sender position ($u_i$) and receiver position ($v_i$) is a single number --
and a one-dimensional space should be *plotted* as one dimension (a 2-D
scatter of a 1-D space piles every point onto a line and looks degenerate
even when the positions are informative). Why not just fit `R = 2`?
Extra dimensions are not free: they can absorb signal the dyadic
covariates would otherwise pick up, and a one-dimensional space is far
easier to display and interpret. Parsimony wins here; to let the data
adjudicate a rank choice, compare fits with `loo::loo_compare()` out
of sample, as done for `R = 1` versus `R = 0` at the end of this
vignette.

The honest display for a 1-D space is a ranked dot plot: every student's
sender and receiver position with a 95% interval, ordered by sender
position. The intervals come from the saved U/V draws
(`posterior_opts = list(save_UV = TRUE)` in the fit above).

```{r overview-uv, fig.width=7, fig.height=6.5, dpi=100, dev="png", fig.alt="Ranked dot plot of the one-dimensional latent positions for all 32 students: for each student, an orange point marks the posterior-mean sender position and a blue point the receiver position, each with a horizontal 95 percent interval; students are ordered by sender position, and a dashed vertical line marks zero."}
lp1 <- subset(latent_positions(fit), dimension == 1)
ordv <- with(subset(lp1, type == "U"), actor[order(value)])
lp1$actor <- factor(lp1$actor, levels = ordv)

ggplot(lp1, aes(x = value, y = actor, color = type)) +
	geom_vline(xintercept = 0, linetype = 2, color = "grey60") +
	geom_errorbar(aes(xmin = value - 2 * posterior_sd,
	                  xmax = value + 2 * posterior_sd),
	              width = 0, alpha = 0.5, orientation = "y") +
	geom_point(size = 1.8) +
	scale_color_manual(values = c(U = "#D55E00", V = "#0072B2"),
	                   labels = c(U = "sender (u)", V = "receiver (v)")) +
	labs(x = "Latent position (dimension 1)", y = NULL, color = NULL) +
	theme_bw() +
	theme(panel.border = element_blank(), panel.grid.major.y = element_blank(),
	      axis.ticks = element_blank(), legend.position = "top")
```

Three things to read off this plot. First, the dimension is genuinely
informative: the extremes (around $+1.7$ at the top versus $-1.1$ at the
bottom) are separated by roughly twice the width of a
typical 95% interval, so the ordering is signal, not noise. Second, sender
and receiver positions are strongly related (their correlation is about
0.55): students
near the top both seek out and attract the same friendship cluster. Third,
the exceptions are the interesting actors -- students whose blue and orange
points diverge are attractive to a cluster they do not themselves reach
toward (or vice versa), a pattern no covariate in the model encodes.
These static positions are pooled across all six waves; if you suspect
students *moved* through this space as the year progressed, that is what
`dynamic_uv = TRUE` estimates (see the
[dynamic effects vignette](dynamic_effects.html)).


## Extracting Latent Positions

If you need the estimated latent positions in a tidy format (for custom plots, merging with external data, or exporting to other tools), the `latent_positions()` function returns a data frame:

```{r}
lp <- latent_positions(fit)
head(lp)
```

Each row gives one actor's position on one latent dimension at one time point, with a `type` column distinguishing the sender position (`"U"`) from the receiver position (`"V"`). The `posterior_sd` column is populated because the fit saved its U/V draws (`posterior_opts = list(save_UV = TRUE)`); without that option it is `NA` with a one-time message telling you how to refit. A static fit returns a single slice (labelled `time = 1`) that applies to every period; dynamic fits (`dynamic_uv = TRUE`) return each time point separately, and `procrustes_align()` removes the arbitrary per-period rotations so that trajectories across time are interpretable (see the [dynamic effects vignette](dynamic_effects.html)).

## Tidy-data round-trip via netify

`lame` uses the `netify` package for tidy network data: build a netify
object from an edgelist and pass it directly to `lame()` or `ame()`. On the
way out, `tidy()`, `autoplot()`, and `prediction_draws_long()` give you
broom-/ggplot-/marginaleffects-ready output.

### Entry: tidy edgelist → `lame()`

```{r tidyverse-entry, eval = requireNamespace("netify", quietly = TRUE), message = FALSE}
library(netify)

# a tidy edgelist is often how network data arrives. Here we melt the
# Dutch-college wave list back into one (from, to, year, tie) row per
# ordered dyad, carrying a dyadic covariate along for the ride.
edges_long <- do.call(rbind, lapply(seq_along(Y), function(t) {
	Yt  <- Y[[t]]
	idx <- which(!is.na(Yt), arr.ind = TRUE)
	data.frame(
		from = rownames(Yt)[idx[, 1]],
		to   = colnames(Yt)[idx[, 2]],
		year = names(Y)[t],
		tie  = as.integer(Yt[idx]),
		same_program = Xdyad[[t]][, , "same_program"][idx]
	)
}))
head(edges_long)

netlet <- netify::netify(
	edges_long,
	actor1 = "from", actor2 = "to", time = "year",
	weight = "tie", dyad_vars = "same_program",
	symmetric = FALSE, mode = "unipartite",
	missing_to_zero = FALSE,
	output_format = "longit_array"
)

# a netify object drops straight into lame(); reduced iterations here so the
# vignette builds quickly (use burn >= 500, nscan >= 4000 for real analyses).
fit_netlet <- lame(
	netlet,
	family = "binary", R = 1,
	burn = 50, nscan = 200, odens = 5,
	verbose = FALSE
)

# to inspect what lame() will see, to_lame() returns the pieces directly
pieces <- netify::to_lame(netlet, lame = TRUE)
names(pieces)
```

`fit_netlet` is a live fit taken straight off the netify object -- no
reshaping, no hand-built lists -- and `names(pieces)` shows the inputs
`lame()` actually consumes (`Y`, `Xdyad`, `Xrow`, `Xcol`, plus the
`mode` / `family` metadata netify carries). The payoff is the alignment:
`netify()` rebuilt, for every period, a tie matrix and a covariate slice
keyed to the same actor order -- exactly the `rownames`-and-`array()`
bookkeeping done by hand in the setup chunk at the top of this vignette.
When your data arrives as an edgelist, as panel network data usually
does, the edgelist is the only object you touch.

### Exit: `tidy()`, `autoplot()`, `prediction_draws_long()`

The rest of the round-trip runs against the `fit` object built from the
Dutch-college data above:

```{r tidy-fit}
# tidy(fit) -> term / estimate / std.error / statistic / p.value /
# conf.low / conf.high; one row per coefficient (or per
# coefficient x period for dynamic_beta fits). The generic ships with
# lame so this runs without broom installed; when broom is on the path
# `broom::tidy(fit)` dispatches to the same method via
# generics::tidy registration.
tidy(fit)
```

### Coefficient tables: `glance()` + `tidy()` -> `modelsummary` / `gt` / `kableExtra`

Because `tidy()` and `glance()` are registered against their `generics`
counterparts, any broom-aware table package consumes a `lame` fit
directly -- most compactly `modelsummary::modelsummary(fit)`:

```{r modelsummary-fit, eval = requireNamespace("modelsummary", quietly = TRUE) && requireNamespace("broom", quietly = TRUE)}
# glance() returns the one-row model summary modelsummary uses for its
# lower panel: nobs, n_actors, n_periods, n_stored, family, mode, R,
# dynamic_uv / dynamic_ab / dynamic_beta, elpd_loo (NA unless save_log_lik=TRUE).
broom::glance(fit)

# pass the fit straight to modelsummary -- tidy() drives the upper
# (coefficients) panel, glance() drives the lower (GOF) panel.
modelsummary::modelsummary(
	list("Dutch college" = fit),
	statistic = "conf.int",
	gof_map = c("nobs", "n_actors", "n_periods", "n_stored",
	            "family", "R", "dynamic_uv", "dynamic_ab")
)
```

Multiple fits compose into a side-by-side table via a named list, one
column per fit. Below we compare the full specification against a
no-latent-space variant (`R = 0`):

```{r fit-nouv, message = FALSE, warning = FALSE}
fit_noUV <- lame(
	Y = Y, Xdyad = Xdyad, Xrow = Xrow, Xcol = Xcol,
	family = "binary", rvar = TRUE, cvar = TRUE, dcor = TRUE,
	R = 0,                          # no latent space
	symmetric = FALSE,
	burn = 30, nscan = 150, odens = 10,
	save_log_lik = TRUE,            # reused by loo_compare() below
	verbose = FALSE, plot = FALSE
)
```

```{r modelsummary-multi, eval = requireNamespace("modelsummary", quietly = TRUE) && requireNamespace("broom", quietly = TRUE), message = FALSE, warning = FALSE}
modelsummary::modelsummary(
	list("AME (R = 1)" = fit, "Additive only (R = 0)" = fit_noUV),
	statistic = "conf.int",
	gof_map = c("nobs", "n_actors", "n_periods", "n_stored",
	            "family", "R", "dynamic_uv", "dynamic_ab")
)
```

Reading down the two columns shows how a latent space can redistribute
signal: when residual clustering is real, homophily coefficients move
toward zero in the AME column because structure that was loaded onto
`same_program` (etc.) is instead absorbed by the latent positions. On
this small panel the movement is modest, with widely overlapping
credible intervals -- the direction is what to look for, not the exact
shift.

```{r autoplot-fit, fig.width = 6, fig.height = 4, fig.alt="Horizontal coefficient plot returned by autoplot.lame: each row is a regression coefficient, with a point at the posterior mean and a horizontal segment marking the 95 percent credible interval; a dashed reference line at zero indicates where an effect is indistinguishable from null."}
# autoplot.lame returns a horizontal coefplot for static fits and a
# ribbon-per-period for dynamic_beta fits, so the same call works for
# both. ggplot layers compose on top.
autoplot(fit) +
	ggtitle("Dutch college friendship: posterior coefficients")
```

The coefplot is the visual form of the `tidy(fit)` table above: the
three homophily terms sit clear of the dashed zero line --
`same_program` largest (0.57), `same_male` next (0.45), `same_smoker`
smallest (0.16) -- while the nodal terms are mixed, with
`smoker_col` and `program_col` straddling the line. The
same generic dispatches on `which = "uv"` (latent positions -- the
ranked dot plot in the Latent Space section above) and `which = "ab"`
(sender / receiver random effects):

```{r autoplot-ab, fig.width = 6, fig.height = 4, fig.alt="Lollipop plot of sender random effects: one stem per student, length equal to the posterior-mean sender effect a_i; long positive stems mark students who nominate far more friends than their covariates predict."}
# which = "ab" delegates to ab_plot() (sender side by default; use
# ab_plot(fit, effect = "receiver") for the column side). Each lollipop
# is a student's sender random effect a_i: the longest positive stems
# are students who nominate far more friends than their covariates
# predict -- the heterogeneity that va summarises as one number.
autoplot(fit, which = "ab")
```

```{r pdl, message = FALSE}
# prediction_draws_long() returns a long-format data frame with
# .chain / .iteration / .draw / period / period_label / i / j /
# actor_i / actor_j / .value -- ready for tidybayes / marginaleffects.
pdl <- prediction_draws_long(fit, type = "response", n_draws = 50)
# the diagonal (i == j) is a self-tie, undefined in a friendship network, but
# the linear predictor is still numerically defined there -- drop it so the
# frame covers only real (off-diagonal) pairs.
pdl <- pdl[pdl$i != pdl$j, ]
head(pdl, 3)
dim(pdl)
```

Because the `.chain` / `.iteration` / `.draw` / `.value` columns follow
the `tidybayes` / `marginaleffects` convention, this frame composes
directly with `ggdist`, `tidybayes`, and `marginaleffects` -- no
reshaping step needed.

## Model comparison: `loo::loo_compare()` as a one-liner

Does the latent space earn its keep over additive effects alone? The
standard out-of-sample yardstick is the expected log pointwise
predictive density (`elpd_loo`) from `loo::loo()`, which requires
`save_log_lik = TRUE` at fit time; the mechanics (and the Pareto-$k$
diagnostics to read before trusting any `elpd` number) are covered in
[Your First AME Model](cross_sec_ame.html). Both fits above already
carry the pointwise log-likelihood, so `loo::loo_compare()` ranks them
in one call:

```{r loo-compare-demo, eval = requireNamespace("loo", quietly = TRUE), message = FALSE, warning = FALSE}
# both fits already carry save_log_lik = TRUE, so this reuses them -- no
# refitting. Pass a NAMED list so the rows are labelled by model rather
# than model1/model2.
cmp <- loo::loo_compare(list(no_latent = loo::loo(fit_noUV),
                             latent_R1 = loo::loo(fit)))
cmp
```

```{r loo-compare-read, echo = FALSE, eval = requireNamespace("loo", quietly = TRUE)}
best  <- rownames(cmp)[1]
worse <- rownames(cmp)[2]
dval  <- abs(round(cmp[2, "elpd_diff"], 1))
sval  <- round(cmp[2, "se_diff"], 1)
```

**How to read the output.** The top row is the preferred model
(`elpd_diff = 0`); each subsequent row reports the elpd *difference*
relative to it, with a standard error. Here `` `r best` `` is
preferred, beating `` `r worse` `` by about `` `r dval` `` elpd units
(SE `` `r sval` ``); a ratio of `|elpd_diff|/se_diff` near or above 2
is a reasonably confident preference, so the latent space is doing
real predictive work on this network. When `|elpd_diff|` is instead
close to its SE, the simpler model is the defensible choice.

This result sits comfortably beside the earlier `modelsummary` table,
where the `R = 1` and `R = 0` columns told the same homophily story
(the nodal terms shift more, with wide overlapping intervals). The two
are not in tension: the latent space leaves the homophily coefficients
intact yet buys roughly `` `r dval` `` elpd units out of sample,
because the predictive gain is carried by the estimated latent
positions themselves -- the $u_i' v_j$ term in the linear predictor.
Residual dyadic structure shows up in held-out predictive density long
before it shows up in a coefficient table.

## Reproducibility: `sessionInfo()` and `renv.lock`

An analysis a reviewer (or your future self) can reproduce records the
seed and the package versions it was fit under. The pattern:

```{r repro-seed, eval = FALSE}
fit <- lame(..., seed = 6886)          # the sampler seed lives in lame()
saveRDS(list(fit = fit, session = sessionInfo()),
        file = "lame_fit_replication.rds")
```

For long-term reproducibility, commit an `renv.lock` so reviewers can
recreate the toolchain with `renv::restore()`. This vignette's own
session record:

```{r repro-sessioninfo}
sessionInfo()
```

## What's Next?

This vignette covered the core workflow: fitting, convergence, GOF, and visualization. For more specialized topics:

- **Single networks (no time series)?** See [Your First AME Model](cross_sec_ame.html) for a detailed cross-sectional walkthrough
- **Two types of nodes?** See [Bipartite Networks](bipartite.html)
- **Evolving network structure?** See [Dynamic Effects](dynamic_effects.html)
- **Just want the quick version?** See [Getting Started](lame.html)

**References**:

Hoff, PD (2021) Additive and Multiplicative Effects Network Models. Statistical Science 36, 34--50.

van de Bunt, G. G., van Duijn, M. A. J., & Snijders, T. A. B. (1999). Friendship networks through time: An actor-oriented dynamic statistical network model. *Computational & Mathematical Organization Theory*, 5(2), 167--192.

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. doi:10.1017/psrm.2021.56
