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

```{r, include = FALSE}
# the models run when the tarball is built; R CMD check re-knits only the text
run_vignette <- !nzchar(Sys.getenv("_R_CHECK_PACKAGE_NAME_"))
if (!run_vignette) knitr::knit_hooks$set(evaluate.inline = function(code, envir) "")
knitr::opts_chunk$set(
	fig.retina = 1,
	collapse = TRUE,
	comment = "#>",
	fig.width = 6,
	fig.height = 4,
	dpi = 72,
	out.width = "100%",
	fig.align = "center",
	eval = run_vignette
)
```

## 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 **L**ongitudinal **A**dditive and **M**ultiplicative **E**ffects 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:

```{r netify-bridge, eval = FALSE}
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](https://netify-dev.github.io/lame/articles/lame-overview.html) 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`.

```{r quick-start, message=FALSE}
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)
```

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:

```{r s3-methods}
# regression coefficients (posterior means)
coef(fit)

# 95% credible intervals
confint(fit)

# one row per coefficient; works with or without broom loaded
tidy(fit)

# one-row model summary; elpd_loo is NA unless a loo result is cached on the fit
glance(fit)

# 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
dim(Y_hat[[1]]) # 25 x 25
cat(
	"Predicted probability range:",
	round(range(unlist(Y_hat), na.rm = TRUE), 3),
	"\n"
)

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

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](forecasting.html)).

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

```{r trace-plot, fig.height=5, fig.alt="MCMC trace and density plots for the regression coefficients (intercept and similarity_dyad); well-mixed traces and unimodal densities indicate adequate convergence."}
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.

```{r gof-plot, fig.width=8, fig.height=5, fig.alt="Longitudinal goodness-of-fit panels: each facet plots one network statistic (sender/receiver degree heterogeneity, dyadic dependence, triadic dependence, transitivity) across time. The observed series is a solid black line with points; the posterior-predictive median is a dashed grey line; the 95 percent credible interval is a grey band."}
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:

```{r uv-plot, fig.width=7, fig.height=6, fig.alt="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."}
uv_plot(fit)
```

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

```{r ab-plot, fig.height=4, fig.alt="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."}
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](https://netify-dev.github.io/lame/articles/lame-overview.html)
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()`:

```{r cross-sectional, message=FALSE}
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)
```

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

```{r dynamic, message=FALSE}
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)
```

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](https://netify-dev.github.io/lame/articles/dynamic_effects.html) 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.

```{r bipartite, message=FALSE}
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)
```

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 `r round(mean(fit_bip$BETA[, "similarity_dyad"]), 2)` and a 95% interval of [`r paste(round(confint(fit_bip)["similarity_dyad", ], 3), collapse = ", ")`], compared with a simulation truth of 0.6. This short example demonstrates extraction of posterior summaries. The [bipartite vignette](https://netify-dev.github.io/lame/articles/bipartite.html) 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:

- **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](https://netify-dev.github.io/lame/articles/dynamic_effects.html).
- **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](https://netify-dev.github.io/lame/articles/dynamic_effects.html).
- **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](cross_sec_ame.html) for a worked example.

The one code pattern worth spelling out is checkpoint / resume:

```{r checkpoint-demo, eval = FALSE}
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)
}
```

## Where to Go Next

| I want to... | Read this |
|--------------|-----------|
| Understand the full workflow with real data | [lame overview](https://netify-dev.github.io/lame/articles/lame-overview.html) |
| Learn about cross-sectional models in depth | [Your first AME model](cross_sec_ame.html) |
| Model two-mode networks | [Bipartite networks](https://netify-dev.github.io/lame/articles/bipartite.html) |
| Let the network structure evolve over time | [Dynamic effects](https://netify-dev.github.io/lame/articles/dynamic_effects.html) |

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