Fast (MCMC-free) AME Estimation

When to use the fast estimator

ame() and lame() fit AME models by Bayesian MCMC, which provides posterior samples and can take minutes. For rapid model exploration, trying several ranks, screening covariates, or getting starting values, the package also ships a fast, MCMC-free point estimator:

• ame_als(): cross-sectional networks
• lame_als(): longitudinal (replicated) networks
• ame_als_bootstrap(): bootstrap standard errors and intervals

It fits the same model

z_ij = mu + beta’x_ij + a_i + b_j + u_i’v_j + e_ij

by iterative block coordinate descent rather than Gibbs sampling. For normal, binary and Poisson outcomes, static and dynamic, the default is an empirical-Bayes / variational estimator of the model that ame() and lame() sample (see below), which also returns Laplace standard errors. It is intended for fast exploration and initialization; runtime depends on the model, convergence tolerance and any bias-correction refits. Use it when a point estimate is enough; use ame() / lame() when the target is a posterior summary.

set.seed(1)
n <- 100
Y <- matrix(rnorm(n * n), n, n)
diag(Y) <- NA
rownames(Y) <- colnames(Y) <- paste0("a", sprintf("%03d", 1:n))

t_als <- system.time(
    fit_als <- ame_als(Y, R = 2, family = "normal", verbose = FALSE)
)
stopifnot(fit_als$converged)
t_mcmc <- system.time(
    fit_mcmc <- ame(
        Y,
        R = 2,
        family = "normal",
        burn = 500,
        nscan = 2000,
        odens = 1,
        seed = 1,
        verbose = FALSE,
        plot = FALSE
    )
)
packageVersion("lame")
data.frame(
    method = c("ALS", "MCMC: 2500 iterations"),
    elapsed_seconds = c(t_als["elapsed"], t_mcmc["elapsed"])
)

A first fit

set.seed(1)
n <- 40
a <- rnorm(n, 0, 0.5); b <- rnorm(n, 0, 0.5)
u <- rnorm(n, 0, 0.8); v <- rnorm(n, 0, 0.8) # a rank-one latent term
Xd <- matrix(rnorm(n * n), n, n)
Y <- 0.5 + 0.8 * Xd + outer(a, b, "+") + tcrossprod(u, v) + matrix(rnorm(n * n), n, n)
diag(Y) <- NA

fit <- ame_als(Y, Xdyad = Xd, R = 1, family = "normal", verbose = FALSE)
fit
#> AME fit via iterative block coordinate descent (fast, MCMC-free)
#> Family: normal | Mode: unipartite | R = 1
#> Network: 40 x 40, 1 time slice
#> Converged: TRUE in 8 iterations
#> Estimator: empirical Bayes (exact Gaussian cell with the variance of the noise
#> by EM, Gaussian prior on the actor effects, covariance by EM, variational
#> posterior of the effects and factors, coefficients and scales by empirical
#> Bayes); no bootstrap or jackknife correction is applied
#> 
#> Coefficients:
#> intercept   Xd_dyad 
#>    0.6046    0.7948 
#> 
#> Variance components:
#>     va    cab     vb    rho     ve 
#> 0.2074 0.0291 0.2027 0.0180 1.0954 
#> 
#> va/vb: variances of the sender/receiver effects | cab: their covariance | rho:
#> dyadic residual reciprocity | ve: residual variance.
#> va, cab and vb are the fitted covariance of the actor effects (the model's
#> random-effect variances). The effects a and b are not re-centred; the intercept
#> is the population intercept.
#> 
#> Laplace standard errors (penalised Fisher information, empirical-Bayes
#> actor-effect prior):
#> intercept   Xd_dyad 
#>    0.1118    0.0264

coef(), fitted(), residuals(), predict() and plot() all work as for any model object. The printed variance components (va, vb, cab, rho, ve) are point estimates. Their construction accounts for the fitted mean model and, where supported, finite-sample bias. They do not carry posterior uncertainty; raw actor-variance summaries are retained separately.

round(cbind(estimate = coef(fit), se = fit$eb$se), 3)
#>           estimate    se
#> intercept    0.605 0.112
#> Xd_dyad      0.795 0.026
round(fit$VC, 3)
#>    va   cab    vb   rho    ve 
#> 0.207 0.029 0.203 0.018 1.095

Choosing the rank R

A fit by empirical Bayes (the default for normal, binary and Poisson data) reports fit$eb$log_marginal, the Laplace approximation, at the posterior mode, to the log marginal likelihood of the data with the effects and the factors integrated out. It does not fall with every added dimension as the deviance does, so the rank with the largest value is the one the data support and the difference between two ranks is an approximate log Bayes factor:

ev <- sapply(
    0:3,
    function(r)
        ame_als(
            Y,
            Xdyad = Xd,
            R = r,
            family = "normal",
            verbose = FALSE
        )$eb$log_marginal
)
data.frame(R = 0:3, log_marginal = round(ev - ev[1], 1))
#>   R log_marginal
#> 1 0          0.0
#> 2 1         69.8
#> 3 2         54.6
#> 4 3         39.7

The fixed-effects fits (bias_correct = "none") have no such value. Fit a few ranks and look for an elbow in the deviance / SSE; if the curve declines smoothly with no clear elbow, keep R small and confirm the chosen rank with an ame() fit.

For binary data, in a simulation of 1,040 networks the rank with the largest value was the generating rank in every network with at least 10 events per factor loading and in 95 percent of those with 3 to 10; with fewer events the data do not determine the rank. Take the lowest rank whose value is within one unit of the largest: it named the generating rank in 96, 54 and 60 percent of the networks with 3 to 10, 1 to 3 and fewer than 1 events per loading, against 95, 49 and 35 percent for the largest value. For normal data the measure of the information is the signal of the product in noise variances per loading (fit$eb$variational$information). In a simulation of 252 networks of 40 actors the rank with the largest value was the generating rank in every network with at least 8 noise variances of signal per loading, in 95 percent of those with 4 to 8 and in 89 percent with 2 to 4. Below 2.1 the data do not determine the product at that rank: the sampler’s posterior mean product is shrunk towards zero and the variational mean is not (1.01 to 6.9 times its norm on 50 simulated networks, against 0.99 to 1.15 from 2.1). Such a fit withholds the product: fit$U, fit$V, fit$UVPM and the positions of latent_positions() are NA, fit$eb$product_withheld records the statistic and the log marginal likelihood of every rank up to the fitted one, and the print names the rank they choose. The coefficients, their standard errors, the variance components, fitted(), predict() and residuals() are reported as usual.

Binary and count fits follow the same rule. Their statistic, signal_per_loading, weighs the squared product by the information of each cell. The number of events per loading does not measure it: a binary fit with 5.6 ties per loading can still be all noise. On 118 simulated networks the products below 2.1 were up to 6.2 (binary) and 7.5 (counts) times the sampler’s norm, and they are withheld. From 2.1 a count product is 0.94 to 1.12 times the sampler’s. A binary product is 0.83 to 1.00 times it, because the variational posterior shrinks a binary product somewhat more than the sampler does. The slopes and their standard errors are unaffected.

The dynamic estimator of a normal panel (lame(method = "als") with dynamic effects, coefficients or positions) withholds in the same way. Static factors use the threshold of 2.1. Moving positions use 1.3 per loading of a period: below it their products were up to 11 times lame()’s on simulated panels.

Non-normal families

For binary and poisson data the default non_normal_method = "irls" runs an iteratively reweighted least squares loop, returning coefficients on the likelihood link scale (probit for binary, log for poisson). The alternative non_normal_method = "transform" is faster but its coefficients are on an uncalibrated rank scale, good for the sign and ranking of effects but not their magnitude.

A fixed-effects binary fit estimates every sender and receiver effect and every latent position from its own actor’s cells, which stretches the latent scale outward in finite samples (the incidental-parameters problem) and leaves an actor without a tie with no finite effect. The default bias_correct = NULL therefore fits a binary model by empirical Bayes: the actor effects get the Gaussian prior that ame() and lame() place on them, its covariance is estimated from the data, and the coefficients solve the Laplace-marginal estimating equations. A penalised reweighted fit gives the posterior mode, and a variational stage then replaces the mode of the effects and of the latent factors by the mean of a Gaussian variational posterior, with the variance of the latent product carried through the likelihood: the mode of a factor that the data barely identify is too large, and that inflates the slope of a covariate that follows it. The scale of the prior of the factors, which the data inform weakly where there are few events per loading, is estimated from the cavity covariances of the rows, which approximate the covariances of their exact posterior, and not from the smaller mean-field ones that drive it to zero; the standard errors carry its uncertainty and fit$eb$variational$scale records it. There is no refitting. Every actor keeps a finite effect (an actor without a tie is shrunk, not dropped), vcov() and confint() return Laplace standard errors for every coefficient, node-covariate coefficients included, and va, cab and vb are the estimated covariance of the random effects. The effects a and b are not re-centred, and the intercept is the population intercept. fit$eb$estimator says which estimator produced the values. A latent rank that the data do not support is not refused: the fit retires the weakest latent column, repeats at the next rank down until its iteration settles, and says so (fit$R is the rank of the fit, fit$R_requested the rank asked for and fit$eb$retirement the record); an iteration that does not settle with no column left is refused. See ?ame_als for the simulation results. The estimator covers directed, bipartite and undirected networks, with or without a latent rank. fitted() and predict() of a binary fit return predictive probabilities: the link integrated over the posterior of the linear predictor instead of applied to its mean, which keeps the probability of a rare held-out event away from zero; plug_in = TRUE returns the plug-in values. A count fit returns the plug-in mean count, as the samplers do (the posterior predictive mean of a count is not finite in this model); plug_in = FALSE returns the mean that integrates the predictor.

A directed binary probit fit also models the dyadic correlation, as ame() and lame() do with dcor = TRUE (the default for directed networks): the two cells of a dyad are one bivariate probit observation, fit$VC["rho"] is the correlation of that model (its posterior standard deviation is in fit$rho_estimation), and the latent factors are not asked to absorb reciprocity. dcor = FALSE fits independent cells, as the earlier release did. The pair likelihood costs about one to two and a half times the independent-cell fit for one period and two to four and a half times for four periods.

A Poisson fit by the default bias_correct = NULL is the same estimator for the model that ame() and lame() sample for counts: a count is Poisson(exp(z)) with z the linear predictor plus a normal noise of variance ve, which the fit estimates, as it does the covariance of the effects, by maximising the marginal likelihood (fit$VC["ve"]). The cell integrates the noise out, so the likelihood, the score and the information of a count are those of the Poisson-lognormal distribution, and the fixed-effects Poisson fit’s failure with an actor that has no count (its effect runs to minus infinity) does not arise: the effect is finite and shrunk, as in the sampler. The linear predictor EZ is on the latent log scale, the mean of a count is exp(EZ + ve / 2), which fitted() and predict() return as the samplers do (plug_in = FALSE gives the mean that integrates the uncertainty of the predictor). A directed count fit also takes the two counts of a dyad as one observation at its fitted state (dcor = TRUE, the default): fit$VC["rho"] is the correlation of the pair model, and the information of the pairs enters vcov(), which matters for a covariate that is the same in both cells of a dyad. With a latent rank and fewer than about three nonzero counts per loading (fit$eb$variational$information, which print() reports) the data barely inform the product: the fit shrinks it more than the sampler does (its norm is 0.95 to 0.98 of that of the lame() posterior mean product from 2.6 nonzero counts per loading and 0.80 to 0.90 up to 2.1), the slopes and their intervals are not affected, and ve is within 0.41 posterior standard deviations of the sampler’s on 22 networks of 1.1 to 6.6 nonzero counts per loading. A binary or count fit whose iteration does not settle at the rank asked for retires its weakest latent column and fits the next rank down (fit$R is the rank fitted, fit$R_requested the rank asked for; ?ame_als gives the numbers). bias_correct = "none" returns the fixed-effects fit of the earlier releases. The next example has an actor that sends no count and one that receives none:

set.seed(4)
nc <- 30
xc <- matrix(rnorm(nc * nc), nc, nc)
ac <- rnorm(nc, 0, 0.7); bc <- rnorm(nc, 0, 0.7)
Yc <- matrix(
    rpois(
        nc * nc,
        exp(
            -0.3 + 0.5 * xc + outer(ac, bc, "+") +
                rnorm(nc * nc, 0, 0.6)
        )
    ),
    nc,
    nc
)
# one actor sends no counts and another receives none
Yc[2, ] <- 0; Yc[, 5] <- 0
diag(Yc) <- NA
fc <- ame_als(Yc, Xdyad = xc, R = 0, family = "poisson", verbose = FALSE)
round(c(coef(fc), ve = fc$VC[["ve"]]), 3)
#> intercept   xc_dyad        ve 
#>    -0.647     0.434     0.384
round(c(sender = fc$a[2], receiver = fc$b[5]), 2)
#>   sender.row2 receiver.row5 
#>         -2.33         -2.41
fx <- tryCatch(
    ame_als(
        Yc,
        Xdyad = xc,
        R = 0,
        family = "poisson",
        bias_correct = "none",
        verbose = FALSE
    ),
    lame_als_convergence = function(e) e
)
class(fx)[1]
#> [1] "lame_als_convergence"
fx$state$termination
#> [1] "no improving conditional step"

The empirical-Bayes fit keeps the actors with no count finite: their effects are shrunk toward the rest of the network, as in the sampler. The fixed-effects fit has no finite solution for these data, because the effect of an actor without a count runs to minus infinity. The call stops with an error of class lame_als_convergence that gives the stop, the score against its threshold and the remedy (the default estimator, or a smaller R). A fixed-effects fit is returned only when it has converged, never at its last iterate.

bias_correct = "none" returns the fixed-effects fit, and "jackknife" and "bootstrap" ask for the two corrections of that fit. The split-sample jackknife refits the model on two interleaved halves of the cells and reports each estimate as twice the full value minus the mean of the two halves; it needs two or more periods at R = 0 and five or more at R > 0, or a single cross-section at R = 0 with at least 25 rarer-outcome cells per actor. The parametric-bootstrap correction simulates the bias from the fitted model (about 60 refits) and reaches the networks the jackknife does not. Both keep the uncorrected values in fit$raw. The bias of the fixed-effects fit falls as one over the number of observed cells per actor times the share of the rarer outcome; ?ame_als gives the size from which it is within 2 percent on the coefficients and 5 percent on va / vb.

set.seed(2)
nb <- 80
ab <- rnorm(nb, 0, 0.5); bb <- rnorm(nb, 0, 0.5)
Xb <- matrix(rnorm(nb * nb), nb, nb)
Yb <- 1 * (0.8 * Xb + outer(ab, bb, "+") + matrix(rnorm(nb * nb), nb, nb) > 0)
diag(Yb) <- NA
fb <- ame_als(Yb, Xdyad = Xb, R = 0, family = "binary", verbose = FALSE)
coef(fb)
#>  intercept    Xb_dyad 
#> 0.07406687 0.79542159
fb$bias_correct$method
#> [1] "eb"
sqrt(diag(vcov(fb)))
#>  intercept    Xb_dyad 
#> 0.08727896 0.02305176

set.seed(6)
nd <- 30; Td <- 6
nmd <- paste0("a", sprintf("%02d", 1:nd))
ad <- rnorm(nd, 0, 0.6); bd <- rnorm(nd, 0, 0.6)
Yd <- lapply(seq_len(Td), function(t) {
    eta <- -0.6 + outer(ad, bd, "+") + matrix(rnorm(nd * nd, 0, 0.3), nd)
    m <- 1 * (eta + matrix(rnorm(nd * nd), nd) > 0)
    diag(m) <- NA
    dimnames(m) <- list(nmd, nmd)
    m
})
fd <- lame(
    Yd,
    family = "binary",
    method = "als",
    dynamic_ab = TRUE,
    verbose = FALSE
)
fd$bias_correct$method
#> [1] "eb"
round(c(rho_ab = fd$rho_ab, sigma_ab = fd$sigma_ab), 2)
#>   rho_ab sigma_ab 
#>     0.70     0.14
fd0 <- lame(
    Yd,
    family = "binary",
    method = "als",
    dynamic_ab = TRUE,
    bias_correct = "none",
    verbose = FALSE
)
round(c(rho_ab = fd0$rho_ab, sigma_ab = fd0$sigma_ab), 2)
#>   rho_ab sigma_ab 
#>     0.80     0.19

set.seed(7)
nn <- 30; Tn <- 6
nmn <- paste0("a", sprintf("%02d", 1:nn))
an <- rnorm(nn, 0, 0.6); bn <- rnorm(nn, 0, 0.6)
xn <- lapply(
    seq_len(Tn),
    function(t)
        array(rnorm(nn * nn), c(nn, nn, 1), dimnames = list(nmn, nmn, "x"))
)
Yn <- lapply(seq_len(Tn), function(t) {
    m <- 0.3 + 0.7 * xn[[t]][, , 1] + outer(an, bn, "+") +
        matrix(rnorm(nn * nn, 0, 0.8), nn)
    diag(m) <- NA
    dimnames(m) <- list(nmn, nmn)
    m
})
fn <- lame(
    Yn,
    Xdyad = xn,
    family = "normal",
    method = "als",
    dynamic_ab = TRUE,
    verbose = FALSE
)
fn$bias_correct$method
#> [1] "eb"
round(cbind(estimate = coef(fn), se = fn$eb$se), 3)
#>           estimate    se
#> intercept    0.568 0.011
#> x_dyad       0.697 0.011
round(c(fn$VC[c("ve", "rho")], rho_ab = fn$rho_ab, sigma_ab = fn$sigma_ab), 3)
#>       ve      rho   rho_ab sigma_ab 
#>    0.644   -0.004    0.655    0.120
fn0 <- lame(
    Yn,
    Xdyad = xn,
    family = "normal",
    method = "als",
    dynamic_ab = TRUE,
    bias_correct = "none",
    verbose = FALSE
)
round(c(rho_ab = fn0$rho_ab, sigma_ab = fn0$sigma_ab), 2)
#>   rho_ab sigma_ab 
#>     0.80     0.16

err <- tryCatch(
    lame(
        Yn,
        Xdyad = xn,
        family = "normal",
        method = "als",
        dynamic_ab = TRUE,
        als_max_iter = 2,
        verbose = FALSE
    ),
    lame_als_convergence = function(e) e
)
class(err)[1]
#> [1] "lame_als_convergence"
class(err$fit)[1]
#> [1] "lame_dynamic_als"

Uncertainty

A fit by the default empirical-Bayes estimator (normal, binary or Poisson) reports Laplace standard errors for its coefficients (vcov(), confint(), summary()): the inverse penalised information of all the parameters, factors included, at the cavity weights of the variational posterior (for a normal response, the information of the Gaussian cell), plus the part of the uncertainty that comes from the scale of the prior of the factors. A dynamic empirical-Bayes fit has Laplace standard errors too (for a normal response the exact Gaussian marginal ones given the hyperparameters) and refuses the bootstrap. For the variance components, the effects and the factors of a static fit, and for every other fit, use the bootstrap:

bt <- ame_als_bootstrap(
    fb,
    R = 30,
    type = "parametric",
    seed = 1,
    verbose = FALSE
)
summary(bt)
#> 
#> ── Bootstrap results: fast AME estimator ───────────────────────────────────────
#> Type: parametric | Family: binary | Mode: unipartite
#> Replicates: 30 valid of 30 (100%); every one converged
#> Optimum check of each replicate: none: an empirical-Bayes refit starts from the
#> variational state of the fit
#> Network: 80 x 80, 1 time slice, R = 0
#> ────────────────────────────────────────────────────────────────────────────────
#> Regression coefficients
#>           Estimate Boot SE 2.5%    97.5%  sig
#> intercept 0.0741   0.0796  -0.0471 0.2126    
#> Xb_dyad   0.7954   0.026   0.7544  0.838  *  
#> * = 95% percentile CI excludes zero
#> ────────────────────────────────────────────────────────────────────────────────
#> Variance components
#>     Estimate Boot SE    2.5%  97.5%
#> va    0.3230  0.0478  0.2328 0.3959
#> cab  -0.0310  0.0407 -0.1310 0.0246
#> vb    0.3203  0.0623  0.2141 0.4184
#> rho  -0.0263  0.0417 -0.0966 0.0575
#> ve    1.0000  0.0000  1.0000 1.0000
#> ────────────────────────────────────────────────────────────────────────────────
#> Additive effects (rotation-invariant, aggregated directly)
#> Sender a: mean |SE| = 0.5266
#> Receiver b: mean |SE| = 0.5301

The Boot SE column is the standard deviation across replicates, whereas confint() reports percentile intervals by default. These summarize different features of the bootstrap distribution and can differ substantially when it is skewed. Every replicate behind these numbers has converged: a refit that has not converged within the iteration budget is retried at four times the budget and otherwise dropped, never kept in an interval, and a bootstrap that drops more than 5 percent of its replicates stops with an error of class lame_als_bootstrap_convergence and returns no interval from the rest (bt$n_failed and bt$fail_reason count the dropped ones). The original fit must have converged too. ame_als_bootstrap() needs at least 10 valid replicates, bootstrap = N in a fit needs N >= 30, and 200 or more give stable intervals. The penalised dynamic normal fit (bias_correct = "none") has the same bootstrap and the same rules. A coefficient outside the central 95 percent of its own replicates is an error of class lame_als_bootstrap_optimum.

One row of the variance-components table is constant by construction: for the binary family the latent residual variance ve is fixed at 1 by probit identification, so its bootstrap SE is exactly zero and its interval is the point [1, 1] – that row is expected, not a failed bootstrap. See ?ame_als for the ve vs ve_working distinction.

type = "parametric" (the default) simulates fresh outcomes from the fitted model and refits. type = "block" resamples time slices and holds the actors fixed; it is available for longitudinal fits with at least ten resampling blocks, and its intervals are read at expanded percentile levels that restore the coverage a plain percentile interval loses over few blocks (?ame_als_bootstrap). Because the actors are fixed, a block bootstrap returns no standard error or interval for the parameters that vary over actors: the node-covariate coefficients, va, vb, cab and ve (and, for an empirical-Bayes fit, the intercept and the shrunk additive effects). They are NA, bt$withheld and the printed summary say why, and type = "parametric" gives them. A binary fit by bias_correct = "none" has no bootstrap, because its dyadic coefficient is biased and every replicate reproduces the bias; bootstrap the default empirical-Bayes fit, or one corrected with bias_correct = "jackknife" or "bootstrap", instead. A fixed-effects Poisson fit gets no interval for its intercept (NA), for the same reason. confint() and vcov() work on the bootstrap object.

An empirical-Bayes fit (the default for normal, binary and Poisson data) returns its Laplace covariance from vcov() / confint() for every coefficient, node covariates included. For a fixed-effects fit (bias_correct = "none") the analytic alternative for the regression coefficients only is the conditional sandwich covariance, which holds the additive and multiplicative effects fixed; the bootstrap is the tool for inference from such a fit.

# a named covariate slice gives both engines the same coefficient names
Xd_arr <- array(
    Xb,
    dim = c(nrow(Xb), ncol(Xb), 1),
    dimnames = list(NULL, NULL, "Xd")
)

# binary ALS fit with 30 bootstrap replicates
fb_b <- ame_als(
    Yb,
    Xdyad = Xd_arr,
    R = 0,
    family = "binary",
    verbose = FALSE,
    bootstrap = 30,
    bootstrap_seed = 1
)

# MCMC fit on the same data with a short chain
fit_mcmc <- ame(
    Yb,
    Xdyad = Xd_arr,
    family = "binary",
    R = 0,
    burn = 30,
    nscan = 100,
    odens = 5,
    verbose = FALSE,
    plot = FALSE,
    gof = FALSE,
    seed = 1
)

# coefficient estimates and intervals as a data frame
coef_interval_table <- function(fit) {
    est <- coef(fit)
    ci <- confint(fit)
    terms <- intersect(names(est), rownames(ci))
    data.frame(
        term = terms,
        estimate = unname(est[terms]),
        conf.low = unname(ci[terms, 1]),
        conf.high = unname(ci[terms, 2]),
        row.names = NULL
    )
}

als_tdy <- coef_interval_table(fb_b)
mcmc_tdy <- coef_interval_table(fit_mcmc)
df <- rbind(
    cbind(
        estimator = "ALS bootstrap",
        als_tdy[, c("term", "estimate", "conf.low", "conf.high")]
    ),
    cbind(
        estimator = "MCMC posterior",
        mcmc_tdy[, c("term", "estimate", "conf.low", "conf.high")]
    )
)

library(ggplot2)
ggplot(
    df,
    aes(
        x = estimate,
        y = term,
        colour = estimator,
        shape = estimator
    )
) +
    geom_vline(xintercept = 0, linetype = "dashed", colour = "grey60") +
    geom_pointrange(
        aes(xmin = conf.low, xmax = conf.high),
        position = position_dodge(width = 0.45),
        size = 0.5
    ) +
    scale_colour_manual(
        values = c(
            "ALS bootstrap" = "#377EB8",
            "MCMC posterior" = "#E41A1C"
        )
    ) +
    scale_shape_manual(
        values = c(
            "ALS bootstrap" = 16,
            "MCMC posterior" = 17
        )
    ) +
    labs(
        x = "Coefficient (probit scale)", y = NULL,
        colour = NULL, shape = NULL
    ) +
    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)
    )

Longitudinal data

lame_als() is the longitudinal counterpart. Its effects are static (pooled across time), as in a non-dynamic lame() fit.

set.seed(3)
Yl <- replicate(5, {
    # intercept 0.4 plus sender and receiver effects and noise
    m <- 0.4 + outer(a, b, "+") + matrix(rnorm(n * n), n, n)
    diag(m) <- NA
    m
}, simplify = FALSE)
lf <- lame_als(Yl, R = 1, family = "normal", verbose = FALSE)
coef(lf)
#> intercept 
#> 0.4974971

Summary

The fast estimator returns point estimates with Laplace (the default for normal, binary and Poisson data), bootstrap or sandwich uncertainty. The MCMC estimator returns posterior summaries.