Your First AME Model

Why Network Models?

Imagine you’re studying friendships in a high school. You have data on who nominated whom as a friend, plus information about each student (gender, race, grade). A natural first instinct is to run a logistic regression: does sharing the same gender predict friendship?

The problem is that friendships aren’t independent observations. Some students are more social (they nominate lots of friends), some are more popular (they receive lots of nominations), and friendships tend to be reciprocated: if Alice names Bob, Bob is more likely to name Alice. A standard regression ignores all of this, and your standard errors will be wrong.

The Additive and Multiplicative Effects (AME) model handles these dependencies directly. It gives each actor a sender effect (\(a_i\), how social they are), a receiver effect (\(b_j\), how popular they are), and a position in a latent space (\(u_i\), \(v_j\)) that captures who tends to connect with whom beyond what the covariates explain. The additive effects enter as \(a_i + b_j\); the latent positions enter multiplicatively as the dot product \(u_i'v_j\). Think of it as a regression that takes network structure seriously.

The Data

We’ll analyze a friendship network from the Add Health study, a longitudinal study of adolescents in the United States. The addhealthc3 dataset in lame contains a directed friendship nomination network along with student characteristics.

library(lame)
library(ggplot2)
set.seed(6886)

# load the Add Health friendship network
data(addhealthc3)

# binary friendship: any nomination counts
Y <- (addhealthc3$Y > 0) * 1
X_nodes <- addhealthc3$X

n <- nrow(Y)
cat("Students:", n, "\n")
#> Students: 32
cat("Friendships:", sum(Y, na.rm = TRUE), "\n")
#> Friendships: 127
cat("Network density:", round(mean(Y, na.rm = TRUE), 3), "\n")
#> Network density: 0.128

Notice the na.rm = TRUE calls: network data often has missing entries (the diagonal is NA because self-ties are undefined). The model handles missing values internally via data augmentation, so NAs can stay in the matrix – including dyads built from missing covariates (here race and grade have a handful), which ame() likewise treats as unobserved.

Before modeling: how much do students vary in their number of friends?

out_degree <- rowSums(Y, na.rm = TRUE)
in_degree <- colSums(Y, na.rm = TRUE)

degree_df <- data.frame(
    Degree = c(out_degree, in_degree),
    Type = rep(c("Nominations sent", "Nominations received"), each = n)
)

ggplot(degree_df, aes(x = Degree)) +
    geom_histogram(binwidth = 1) +
    facet_wrap(~Type, ncol = 2) +
    labs(x = "Number of ties", y = "Count") +
    theme_bw() +
    theme(
        panel.border = element_blank(),
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "top",
        plot.title = element_text(hjust = 0),
        plot.subtitle = element_text(hjust = 0),
        plot.caption = element_text(hjust = 0),
        strip.background = element_rect(fill = "black", color = "black"),
        strip.text = element_text(color = "white", hjust = 0)
    )

Some students are much more social than others, and some much more popular – exactly what the sender and receiver random effects will capture.

Building Covariates

The classic question in friendship networks is homophily: do birds of a feather flock together? We’ll test whether students are more likely to be friends with others of the same gender, race, and grade.

We deliberately include both same_grade (binary) and grade_diff (continuous: grades apart) to illustrate a common modeling pitfall, collinearity: they measure nearly the same thing, and we’ll see the consequences in the results below.

# outer() compares every pair of students: 1 if they share the attribute, 0 if not
same_female <- outer(X_nodes[, "female"], X_nodes[, "female"], "==") * 1
same_race <- outer(X_nodes[, "race"], X_nodes[, "race"], "==") * 1
same_grade <- outer(X_nodes[, "grade"], X_nodes[, "grade"], "==") * 1
# absolute grade difference, a continuous covariate
grade_diff <- abs(outer(X_nodes[, "grade"], X_nodes[, "grade"], "-"))

# stack into an n x n x p array
Xdyad <- array(NA, dim = c(n, n, 4))
Xdyad[, , 1] <- same_female
Xdyad[, , 2] <- same_race
Xdyad[, , 3] <- same_grade
Xdyad[, , 4] <- grade_diff
dimnames(Xdyad)[[3]] <- c('same_female', 'same_race', 'same_grade', 'grade_diff')
for (k in 1:4) diag(Xdyad[, , k]) <- NA

# nodal covariates (sender and receiver characteristics)
Xrow <- X_nodes[, c("female", "grade")]
Xcol <- X_nodes[, c("female", "grade")]

Fitting the Model

Now the fun part. We fit a binary probit AME model with:

• Covariates for homophily effects
• Sender and receiver random effects (sociality / popularity)
• Dyadic correlation (reciprocity)
• 2-dimensional latent space (residual clustering)

fit <- ame(
    Y,
    Xdyad = Xdyad,
    Xrow = Xrow,
    Xcol = Xcol,
    R = 2, # 2D latent space
    family = "binary", # probit model for 0/1 data
    rvar = TRUE, # sender random effects
    cvar = TRUE, # receiver random effects
    dcor = TRUE, # dyadic correlation (reciprocity)
    burn = 100, # short burn-in for this example
    nscan = 500, # short run for this example
    odens = 25, # thinning
    verbose = FALSE,
    posterior_opts = posterior_options(
        save_ab = TRUE, save_UV = TRUE,
        thin_ab = 1, thin_UV = 1
    ),
    gof = TRUE
)

What Did We Find?

summary(fit)
#> 
#> ── AME model summary ───────────────────────────────────────────────────────────
#> 
#> ── Call 
#> Y ~ intercept + dyad(same_female, same_race, same_grade, grade_diff) + row(female, grade) + col(female, grade) + a[i] + b[j] + rho*e[ji] + U[i,1:2] %*% V[j,1:2], family = 'binary'
#> 
#> ── Regression coefficients 
#>                  Estimate StdError z_value p_value CI_lower CI_upper  Rhat
#> intercept          -3.547    1.184  -2.996   0.003   -5.765   -1.815 1.035
#> female_row         -0.655     0.38  -1.724   0.085   -1.276    0.023 0.995
#> grade_row           0.033    0.066   0.494   0.621   -0.099    0.134 1.046
#> female_col         -0.567    0.245  -2.311   0.021   -1.115   -0.164 1.022
#> grade_col            0.29    0.091   3.191   0.001    0.093    0.483 1.012
#> same_female_dyad    0.256    0.176   1.456   0.145   -0.092    0.549 0.974
#> same_race_dyad      0.147    0.311   0.475   0.635   -0.413    0.755 1.012
#> same_grade_dyad     0.223    0.249   0.894   0.371   -0.318    0.695 1.037
#> grade_diff_dyad    -0.652    0.103  -6.343       0   -0.825   -0.444 1.107
#>                     ESS    
#> intercept            40  **
#> female_row           40   .
#> grade_row            40    
#> female_col       25.095   *
#> grade_col            40  **
#> same_female_dyad     40    
#> same_race_dyad       40    
#> same_grade_dyad      40    
#> grade_diff_dyad      40 ***
#> ---
#> 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  Rhat    ESS
#> va     0.414    0.168 1.112 40.000
#> cab    0.037    0.113 1.006 40.000
#> vb     0.266    0.187 1.037 40.000
#> rho    0.705    0.128 1.139 13.268
#> ve     1.000    0.000    NA  0.000
#>   va = sender variance, cab = sender-receiver covariance, vb = receiver
#>   variance, rho = dyadic correlation, ve = residual variance

Let’s unpack the key results:

Homophily effects. The grade_diff coefficient is strongly negative and its credible interval excludes zero: students further apart in grade are much less likely to be friends. The same_female coefficient is positive (gender homophily) but its 95% credible interval includes zero: the data do not separate the effect from zero. The same_race coefficient captures race homophily after controlling for grade and gender.

The collinearity pitfall. same_grade and grade_diff measure nearly the same thing – a grade difference of 0 is the same-grade event – so the model cannot cleanly split the grade signal between them. It still identifies the dominant effect: grade_diff is strongly negative and stable from run to run (about -0.65, never near zero). What suffers is the partition: same_grade is left to absorb a small, noisy residual bump on top of the linear grade-distance trend. The lesson is about interpretation, not estimate instability: include one or the other, not both, so the grade effect lands on a single, clean coefficient.

Variance components. The sender variance (va) and receiver variance (vb) quantify how much students differ in sociality and popularity. The dyadic correlation (rho) captures reciprocity: values near 1 mean that if A nominates B, B almost always nominates A back. Reciprocity is strong here.

For binary networks with a probit link, the coefficients are on the latent scale, and the textbook “\(0.4 \times \beta\)” rule for converting them to probability changes is only sharp when the baseline probability is near 0.5 – which this sparse network is not. Rather than juggling correction factors, compute the quantity you actually want with predict(fit, type = "response") and contrast predicted probabilities under counterfactual covariate settings.

Did the Sampler Converge?

The model is estimated via MCMC (Markov chain Monte Carlo), so trace plots are the first check that the sampler explored the posterior thoroughly.

trace_plot(fit, params = "beta", ncol = 3)

The regression coefficients are not the whole story, though. ame() runs two independently seeded chains by default (chains = 2; the seeds are seed and seed + 1000). At the end of the run it checks the intercept, va, vb, ve and rho of each chain for drift, and it compares the two chains on the coefficients, the variance components and summaries of the additive and latent effects (split rank-normalised \(\hat R\) below 1.05; fit$rhat lists them). Chains that are not settled are extended together (up to three times nscan more sweeps) and then refused with an error of class lame_chain_unsettled, so a fit that is returned has passed both checks. The trace of the dyadic correlation rho is the one to look at next:

trace_plot(fit, params = "variance", include = "Dyadic correlation")

What to look for: Traces should fluctuate around a stable mean (long trends mean non-convergence) and densities should be smooth and unimodal. This fit stores 40 post-burn-in draws from 2 chains, enough to illustrate the workflow; the runs are short, and the gate above, not the length of the trace, is what a fit has to pass.

library(posterior)
draws <- posterior::as_draws(fit) # draws_array [iter, chain, var]
# mean, sd, quantiles, rhat and ESS for each parameter
posterior::summarise_draws(draws)
#> # A tibble: 14 × 10
#>    variable          mean  median     sd    mad      q5      q95   rhat ess_bulk
#>    <chr>            <dbl>   <dbl>  <dbl>  <dbl>   <dbl>    <dbl>  <dbl>    <dbl>
#>  1 intercept      -3.55   -3.42   1.18   1.10   -5.15   -1.85     1.03      34.5
#>  2 female_row     -0.655  -0.635  0.380  0.345  -1.22    0.00783  0.995     60.4
#>  3 grade_row       0.0328  0.0390 0.0664 0.0670 -0.0820  0.132    1.05      23.4
#>  4 female_col     -0.567  -0.578  0.245  0.196  -1.07   -0.201    1.02      48.8
#>  5 grade_col       0.290   0.288  0.0908 0.0578  0.115   0.452    1.01      41.3
#>  6 same_female_d…  0.256   0.276  0.176  0.172  -0.0217  0.457    0.974     64.1
#>  7 same_race_dyad  0.147   0.0934 0.311  0.367  -0.301   0.573    1.01      36.2
#>  8 same_grade_dy…  0.223   0.187  0.249  0.159  -0.213   0.615    1.04      38.8
#>  9 grade_diff_dy… -0.652  -0.648  0.103  0.102  -0.813  -0.467    1.11      16.3
#> 10 va              0.414   0.357  0.168  0.138   0.218   0.685    1.11      21.6
#> 11 cab             0.0368  0.0161 0.113  0.121  -0.114   0.176    1.01      40.0
#> 12 vb              0.266   0.215  0.187  0.0616  0.125   0.594    1.04      30.1
#> 13 rho             0.705   0.714  0.128  0.123   0.504   0.897    1.14      12.9
#> 14 ve              1       1      0      0       1       1       NA         NA  
#> # ℹ 1 more variable: ess_tail <dbl>

# four longer chains for a convergence assessment
fit_mc <- ame(
    Y,
    Xdyad = Xdyad,
    Xrow = Xrow,
    Xcol = Xcol,
    family = "binary",
    R = 2,
    burn = 1000,
    nscan = 5000,
    odens = 25,
    chains = 4,
    cores = 4, # forked on Unix; sequential on Windows
    verbose = FALSE
)
posterior::summarise_draws(posterior::as_draws(fit_mc))

str(fit$mh_counters)
#> List of 6
#>  $ beta         : num 0
#>  $ Sab          : num 0
#>  $ s2           : num 0
#>  $ rho          : num 0
#>  $ UV           : num 0
#>  $ rho_collapsed: num 0

prior_summary(fit)
#> 
#> ── Recorded prior settings (ame fit) ──
#> 
#> Regression coefficients have a zero-mean design-based Gaussian g-prior,
#> independent of residual s2, with a proportional numerical ridge and additional
#> intercept variance; see rbeta_ab_fc().
#> `g` (top-level argument) = 992
#> Note: `g` is set at top level on `ame()` -- not inside `prior = list(...)`.
#> `Sab0` = "matrix(c(0.2032, 0, 0, 0.2032), 2, 2)"
#> `eta0` = "8"
#> `Suv0` = "matrix(c(0.2032, 0, 0, 0, 0, 0.2032, 0, 0, 0, 0, 0.2032, 0, 0, 0, 0,
#> 0.2032), 4, 4)"
#> `kappa0` = "13"

# two candidates with independent residuals and saved log-likelihoods
fit_ll <- ame(
    Y,
    Xdyad = Xdyad,
    Xrow = Xrow,
    Xcol = Xcol,
    family = "binary",
    R = 2,
    dcor = FALSE,
    burn = 1000,
    nscan = 5000,
    odens = 25,
    save_log_lik = TRUE,
    verbose = FALSE
)

# alternative model without the dyadic covariates
fit_ll_null <- ame(
    Y,
    Xrow = Xrow,
    Xcol = Xcol,
    family = "binary",
    R = 2,
    dcor = FALSE,
    burn = 1000,
    nscan = 5000,
    odens = 25,
    save_log_lik = TRUE,
    verbose = FALSE
)

loo_dyad <- loo::loo(fit_ll)
loo_null <- loo::loo(fit_ll_null)

loo_dyad # check the Pareto k table first

loo::loo_compare(list(with_dyad = loo_dyad, no_dyad = loo_null))

Does the Model Fit the Data?

Goodness-of-fit (GOF) checks whether the model can reproduce structural features of the observed network – degree heterogeneity, reciprocity, clustering – beyond dyad-level prediction.

gof_plot(fit)

In each panel the grey histogram is the posterior-predictive distribution from networks simulated out of the fitted model, and the dashed black vertical line is the observed value. See ?gof_plot for the statistics = ... aliases and the exact column names stored in fit$GOF.

The observed values for sender and receiver degree heterogeneity and for dyadic dependence (essentially empirical reciprocity) fall inside their posterior-predictive histograms – those three features are explicitly modeled by the random effects, the variance components, and \(\rho\). The two triad-level statistics tell a different story: observed triadic dependence sits toward the upper end of its histogram, and observed transitivity in the upper tail. AME under-predicts triangle closure on this network, the expected signature of conditional dyadic independence: the latent space soaks up some triangle structure through clustering of the \(u_i\), but there is no explicit closure term. If transitivity is your substantive target, fit an ERGM; if it is a diagnostic concern only, increasing R sometimes helps modestly.

You can also compute GOF after the fact using the gof() function, which accepts custom statistics. The built-in trans.dep panel is a correlation-style clustering measure; ERGM users will typically want the raw clustering coefficient (sna::gtrans, the scalar summary that gwesp(decay = 0, fixed = TRUE) targets), which takes a few lines as a custom stat:

# clustering coefficient: transitive triples over two-paths
trans_ratio <- function(Y) {
    Yb <- Y; Yb[is.na(Yb)] <- 0 # missing ties count as 0
    YY <- Yb %*% Yb
    triangles <- sum(YY * Yb) # transitive triples
    two_path <- sum(YY) - sum(diag(YY)) # two-paths (potential triangles)
    c(trans_ratio = triangles / max(two_path, 1))
}

gof_custom <- gof(fit, custom_gof = trans_ratio, nsim = 50, verbose = FALSE)
obs <- gof_custom[1, "trans_ratio"]
sim <- gof_custom[-1, "trans_ratio"]

# posterior predictive p-value (right tail): near 0 means under-predicted
c(
    observed = round(obs, 3),
    sim_mean = round(mean(sim), 3),
    pp_p_right = round(mean(sim >= obs), 3)
)
#>   observed   sim_mean pp_p_right 
#>      0.485      0.401      0.060

Read this the way you would read ergm::gof(): pp_p_right is the fraction of simulated networks at or above the observed value, so a value near 0 means the model under-predicts the statistic. Here it is 0.06 across 50 joint posterior replicates. These replicates retain uncertainty in the actor effects and latent positions as well as the regression coefficients. See ?gof for adding further custom statistics (raw triangle counts, degree assortativity, and so on).

Visualizing the Latent Space

The multiplicative effects (\(u_i'v_j\)) place each student in a 2D latent space: students near each other in sender space (triangles) nominate similar friends, and students near each other in receiver space (circles) are nominated by similar students. Clusters likely share some unobserved characteristic (same social group, same extracurriculars) that drives friendship choices beyond gender, race, and grade.

uv_plot(fit, layout = "biplot", show.edges = FALSE, label.nodes = FALSE)

Who Are the Most Social / Popular Students?

The additive effects decompose individual heterogeneity into sender effects (sociality) and receiver effects (popularity):

p1 <- ab_plot(
    fit,
    effect = "sender",
    sorted = TRUE,
    title = "Sender effects (sociality)"
)
p2 <- ab_plot(
    fit,
    effect = "receiver",
    sorted = TRUE,
    title = "Receiver effects (popularity)"
)

library(patchwork)
p1 / p2

Students with large positive sender effects nominate more friends than expected given their covariates; those with large positive receiver effects receive more nominations.

Making Predictions

The model gives you predicted probabilities for every possible tie, useful for link prediction and for understanding dyad-level fit. For sparse networks, predicting “no tie” for every dyad already gets most predictions right, so raw accuracy flatters the model – judge it against the modal-class baseline.

pred_resp <- predict(fit, type = "response")

# observed and predicted ties as vectors
Y_vec <- as.vector(Y)
pred_vec <- as.vector(pred_resp)
keep <- !is.na(Y_vec)

# classification at a 0.5 threshold
pred_binary <- (pred_vec[keep] > 0.5) * 1
confusion <- table(Actual = Y_vec[keep], Predicted = pred_binary)
knitr::kable(confusion, caption = "Confusion matrix (threshold = 0.5)")

accuracy <- sum(diag(confusion)) / sum(confusion)
baseline <- max(mean(Y_vec[keep]), 1 - mean(Y_vec[keep])) # accuracy of always predicting the modal class
cat(
    "\nAccuracy:",
    round(accuracy, 3),
    " | Baseline (predict modal class):",
    round(baseline, 3),
    " | Lift:",
    round(accuracy - baseline, 3),
    "\n"
)
#> 
#> Accuracy: 0.915  | Baseline (predict modal class): 0.872  | Lift: 0.043

At density ~13% the modal-class baseline is ~0.87, so the relevant question is whether the model’s accuracy is materially above 0.87, not whether it crosses 0.5.

For link prediction, AUC-style metrics evaluate the model’s ability to rank true ties above non-ties across all thresholds. For sparse networks PR-AUC is the more informative summary: a model that always predicts “no tie” already has AUROC near 0.5 but PR-AUC near the network density, so the gap between PR-AUC and density is the metric of substantive interest.

# AUROC with pROC
roc_obj <- pROC::roc(
    response = Y_vec[keep],
    predictor = pred_vec[keep],
    quiet = TRUE
)
auroc <- as.numeric(pROC::auc(roc_obj))

# ROC and precision-recall AUCs with precrec; PR-AUC suits rare ties
ev <- precrec::evalmod(scores = pred_vec[keep], labels = Y_vec[keep])
precrec::auc(ev) # AUROC and PR-AUC side by side
#>   modnames dsids curvetypes      aucs
#> 1       m1     1        ROC 0.9208775
#> 2       m1     1        PRC 0.7479777

cat(
    "AUROC:",
    round(auroc, 3),
    "\n",
    "Baseline AUROC (random ranker):",
    0.5,
    "\n",
    "Baseline PR-AUC (density of Y):",
    round(mean(Y_vec[keep]), 3),
    "\n"
)
#> AUROC: 0.921 
#>  Baseline AUROC (random ranker): 0.5 
#>  Baseline PR-AUC (density of Y): 0.128

Here PR-AUC is 0.748 against a density baseline of 0.128, and AUROC is 0.921 against 0.5 for a random ranker. These are in-sample ranking summaries. The held-out estimates below check predictions on masked dyads; independent-chain convergence and a larger evaluation are needed before making a substantive predictive claim.

# 80/20 stratified mask of off-diagonal dyads
set.seed(6886)
off_diag <- which(row(Y) != col(Y) & !is.na(Y))
pos_idx <- off_diag[Y[off_diag] == 1]
neg_idx <- off_diag[Y[off_diag] == 0]
test_idx <- c(
    sample(pos_idx, round(0.2 * length(pos_idx))),
    sample(neg_idx, round(0.2 * length(neg_idx)))
)

Y_train <- Y
Y_train[test_idx] <- NA # mask test dyads from the likelihood

fit_train <- ame(
    Y_train,
    Xdyad = Xdyad,
    Xrow = Xrow,
    Xcol = Xcol,
    R = 2,
    family = "binary",
    rvar = TRUE,
    cvar = TRUE,
    dcor = TRUE,
    burn = 30,
    nscan = 100,
    odens = 5, # short run for the vignette
    verbose = FALSE,
    gof = FALSE
)
#> ℹ 92 `Y` cells had a missing covariate and are treated as unobserved.
#> ℹ The covariate coefficients are estimated from the complete dyads only.

pred_train <- predict(fit_train, type = "response")

# evaluate_heldout() scores the predictions on a logical test mask
test_mask <- matrix(
    FALSE,
    nrow(Y),
    ncol(Y),
    dimnames = dimnames(Y)
)
test_mask[test_idx] <- TRUE
heldout_metrics <- evaluate_heldout(
    y_obs = Y,
    y_pred = pred_train,
    mask = test_mask,
    family = "binary"
)
heldout_metrics
#>   n_eval     auroc     auprc      brier   logloss
#> 1    198 0.8245087 0.5611221 0.08211174 0.2768877
cat(
    "Test dyads:",
    length(test_idx),
    " | positives:",
    sum(Y[test_idx] == 1),
    "\n"
)
#> Test dyads: 198  | positives: 25

Simulating from the Model

You can generate new networks from the fitted posterior. A posterior predictive check asks: does a network simulated from posterior draws look like the observed one on a feature you care about (here, overall density)? Observed values inside the simulated distribution mean the model reproduces that feature; values in a tail mean it is missing something.

sims <- simulate(fit, nsim = 100)

# simulated and observed densities
sim_densities <- sapply(sims$Y, function(y) mean(y, na.rm = TRUE))
obs_density <- mean(Y, na.rm = TRUE)

ggplot(data.frame(density = sim_densities), aes(x = density)) +
    geom_histogram(bins = 20) +
    geom_vline(
        xintercept = obs_density,
        linewidth = 1,
        linetype = "dashed"
    ) +
    labs(
        title = "Posterior predictive check: network density",
        subtitle = "Dashed line = observed; histogram = simulated from the model",
        x = "Network density (mean tie probability)", y = "Replicate count"
    ) +
    theme_bw() +
    theme(
        panel.border = element_blank(),
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "top",
        plot.title = element_text(hjust = 0),
        plot.subtitle = element_text(hjust = 0),
        plot.caption = element_text(hjust = 0)
    )

Practical Tips

Choosing R (latent dimensions). Start with R = 0 (no latent space), then try R = 1 and R = 2 and compare GOF. There is no cap on R below the network size: ranks above floor(n/3) fit, but the chain mixes more slowly and the fitted values lose accuracy gradually past the structure in the data; in practice R = 2 or R = 3 is the right default.

MCMC settings. The defaults (burn = 500, nscan = 10000, odens = 25) are the length of a final run: a chain that is still moving at the end is extended up to three times and then refused (prior = list(chain_extensions = 0) refuses at the first flag). burn = 500, nscan = 2000, odens = 25 is faster for exploratory work, and the same check applies to it. Check convergence with trace_plot() and between-chain split R-hat.

Variance components. Include rvar = TRUE and cvar = TRUE unless you have a specific reason not to. Include dcor = TRUE for directed networks where reciprocity is plausible.

What’s Next?

• Multiple time periods? See the lame overview for longitudinal models
• Two types of nodes? See the bipartite vignette
• Evolving network structure? See the dynamic effects vignette
• Other data types? See ?ame for the supported family options

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.