Introduction
Many real-world networks are bipartite: two distinct
types of nodes, with ties forming only between types, never within them.
Students enroll in courses, countries sign treaties, legislators join
committees. The adjacency matrix is rectangular rather than square, and
the row and column nodes play fundamentally different roles.
Both ame() (cross-sectional) and lame()
(longitudinal) support bipartite networks. The key difference from the
unipartite case is that the model uses separate latent spaces for the
two node types, connected by an interaction matrix \(G\).
The Bipartite AME Model
In a bipartite network, the probability of a tie between row node
\(i\) and column node \(j\) is modeled as:
\[y_{ij} = \beta'x_{ij} + a_i + b_j +
u_i'Gv_j + \epsilon_{ij}\]
The additive effects work as in the unipartite case (see the cross-sectional vignette): \(a_i\) captures how “active” row node \(i\) is, \(b_j\) how “popular” column node \(j\) is. The multiplicative term \(u_i'Gv_j\) replaces the unipartite
\(u_i'v_j\): row and column nodes
live in potentially different-dimensional latent spaces (\(R_\text{row}\) and \(R_\text{col}\)), and \(G\) is a \(R_\text{row} \times R_\text{col}\)
interaction matrix.
How \(G\) is handled differs between
the cross-sectional and longitudinal entry points:
- In
ame() (cross-sectional bipartite),
\(G\) is initialised to a
rank-truncated scaled identity and held fixed during
MCMC. The rotation and scaling of the interaction are absorbed into
\(U\) and \(V\), so \(G\) itself is not estimated; what is
identified is the overall multiplicative term \(U G V'\).
- In
lame() (longitudinal bipartite),
the default MCMC estimator samples \(G\) each iteration from its conditional
Gaussian. With dynamic_G = FALSE a single \(G\) is shared across periods; with
dynamic_G = TRUE a separate \(G_t\) is drawn per period (see the
dynamic_G section below, which also covers the fast ALS
path).
In either case, the row-side positions (\(U\), shape \(n_\text{row} \times R_\text{row}\)) and
column-side positions (\(V\), shape
\(n_\text{col} \times R_\text{col}\))
can be chosen at different ranks.
Two terms are absent in the bipartite setting: the dyadic correlation
\(\rho\) (a row actor is never also a
column actor, so there is no reciprocal tie to correlate), and the joint
covariance of \(a\) and \(b\), which instead have independent
variances since they describe different types of actors.
Supported families in bipartite mode
Every family that works for a square network also works for a
rectangular one. The rectangular Z-samplers live in
R/rZ_bipartite.R and treat each cell as independent given
the linear predictor (there is no reciprocal-cell coupling in a
bipartite graph). The supported set:
normal |
yes |
the unrestricted-Gaussian baseline. |
binary |
yes |
probit link on tie / no-tie. |
ordinal |
yes |
rectangular ordinal-probit sampler (rZ_ord_bip_fc); use
data-induced cutpoints (default) or pass
ordinal_cutpoints = "explicit" for the Cowles MH variant on
the unipartite path. |
cbin |
yes |
constrained-binary with odmax interpreted as the
per-row outdegree cap. |
frn |
yes |
fixed-rank-nominations; each row’s top-odmax ranks are
observed. |
poisson |
yes |
rectangular MH on Z; supports the period_exposure
offset on the longitudinal path. |
The same table applies to lame() in bipartite mode.
Time-varying interaction matrices (dynamic_G = TRUE) are
covered in their own section below, on both the MCMC and ALS paths.
Cross-Sectional Analysis with ame()
Simulating Data
We simulate a bipartite network of 30 students and 20 courses, each
with a 2-dimensional latent position; the interaction matrix \(G\) determines how those dimensions combine
to predict enrollment.
# simulate a student-course enrollment network
n_students <- 30
n_courses <- 20
# true latent positions (unobserved in practice)
U_true <- matrix(rnorm(n_students * 2, 0, 0.8), n_students, 2)
V_true <- matrix(rnorm(n_courses * 2, 0, 0.8), n_courses, 2)
# interaction matrix: dimension 1 has positive affinity,
# dimension 2 has negative (students high on dim 2 avoid courses high on dim 2)
G_true <- matrix(c(1, 0.5, 0.5, -1), 2, 2)
# generate enrollment probabilities and binary outcomes
# negative intercept keeps enrollment rate realistic (~25-30%)
eta <- -0.8 + U_true %*% G_true %*% t(V_true)
prob <- pnorm(eta)
Y_bipartite <- matrix(rbinom(n_students * n_courses, 1, prob),
n_students, n_courses)
rownames(Y_bipartite) <- paste0("Student", 1:n_students)
colnames(Y_bipartite) <- paste0("Course", 1:n_courses)
cat("Network dimensions:", dim(Y_bipartite), "\n")
#> Network dimensions: 30 20
cat("Enrollment rate:", round(mean(Y_bipartite), 2), "\n")
#> Enrollment rate: 0.29
Fitting the Model
Fitting a bipartite model requires setting
mode = "bipartite" and specifying the latent dimensions for
each node type separately via R_row and
R_col.
# burn and nscan are small so the vignette builds quickly
# for real analyses, use burn >= 1000 and nscan >= 5000.
fit_cross <- ame(
Y = Y_bipartite,
mode = "bipartite",
R_row = 2, # latent dimensions for students
R_col = 2, # latent dimensions for courses
family = "binary",
burn = 100,
nscan = 500,
odens = 5,
verbose = FALSE,
# save thinned U/V draws so latent_positions() can report posterior SDs
posterior_opts = posterior_options(save_UV = TRUE)
)
summary(fit_cross)
#>
#> === AME Model Summary ===
#>
#> Call:
#> [1] "Y ~ a[i] + b[j] + U[i,1:2] %*% G %*% V[j,1:2]', family = 'binary'"
#>
#> Regression coefficients:
#> ------------------------
#> Estimate StdError z_value p_value CI_lower CI_upper
#> intercept -0.813 0.209 -3.881 0 -1.228 -0.451 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Note: stars are a visual hint from posterior mean / SD only; for inference use the credible intervals.
#>
#> Variance components:
#> -------------------
#> Estimate StdError
#> va 0.301 0.085
#> cab 0.000 0.000
#> vb 0.406 0.123
#> ve 1.000 0.000
#> rho 0.000 0.000
#> (va = sender, cab = sender-receiver covariance, vb = receiver,
#> rho = dyadic correlation, ve = residual variance)
#> Note: bipartite model (rho fixed to 0, cab fixed to 0)
The summary shows the intercept (we included no covariates) and the
variance components. The fitted object also stores student latent
positions (U), course latent positions (V),
the (fixed) interaction matrix (G), and additive effects
(APM for student activity, BPM for course
popularity).
Visualizing the Latent Space
uv_plot with layout = "biplot" places both
node types in the same latent space: students and courses that appear
close together have a higher predicted probability of a tie.
uv_plot(fit_cross, layout = "biplot") +
ggtitle("Bipartite Latent Space: Students and Courses")
Because ame() holds \(G\) fixed (its entries carry no information
– rotation and scaling are absorbed into \(U\) and \(V\)), the estimated quantity to inspect is
the overall multiplicative term \(U G
V'\), which we visualise directly:
G_mult <- fit_cross$U %*% fit_cross$G %*% t(fit_cross$V)
G_df <- data.frame(
Student = rownames(G_mult)[row(G_mult)],
Course = colnames(G_mult)[col(G_mult)],
Fitted = as.vector(G_mult)
)
# order axes numerically (Student1, Student2, ...) instead of the default
# alphabetical factor order (Student1, Student10, Student11, ...)
G_df$Student <- factor(G_df$Student, levels = paste0("Student", n_students:1))
G_df$Course <- factor(G_df$Course, levels = paste0("Course", 1:n_courses))
# rdbu reversed (direction = -1): positive => red, negative => blue,
# zero => white. By default scale_fill_distiller() centres white at the
# midpoint of the DATA range, not at zero -- so we pass symmetric limits
# c(-mx, mx) to force the white midpoint onto zero, which is the right
# choice for a signed effect. We cap the fill at the 99th percentile of
# |UGV'| (clamping the few more extreme cells to the cap) so a single
# outlying cell cannot wash the rest of the map to near-white. RdBu is
# colour-blind-safe per ColorBrewer's protanopia / deuteranopia checks.
mx <- quantile(abs(G_df$Fitted), 0.99, na.rm = TRUE)
G_df$Fitted_capped <- pmax(pmin(G_df$Fitted, mx), -mx)
ggplot(G_df, aes(x = Course, y = Student, fill = Fitted_capped)) +
geom_tile() +
scale_fill_distiller(palette = "RdBu", direction = -1,
limits = c(-mx, mx)) +
labs(title = "Fitted Multiplicative Structure (UGV')",
subtitle = "Red = positive, blue = negative; colour scale capped at the 99th percentile",
x = "Course", y = "Student",
fill = "U G V'") +
theme_bw() +
theme(panel.border = element_blank(),
axis.ticks = element_blank(),
legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.y = element_text(size = 6))
Longitudinal Analysis with lame()
When bipartite networks are observed over multiple time periods
(enrollments across semesters, country–treaty memberships across
decades), lame() pools information across time while
optionally allowing the latent structure to evolve.
Simulating Longitudinal Data
We simulate a panel of user–item networks over 5 periods whose latent
positions drift via an AR(1) process with \(\rho = 0.9\). We use the
normal family throughout this section because the bipartite
dynamic sampler mixes well for continuous outcomes at vignette-scale
chains; family = "binary" runs on the same machinery but
wants burn >= 1000 and nscan >= 5000.
One parameterisation note: the chunk uses the
stationary-variance form (\(U_t = \rho\,U_{t-1} + \sqrt{1 -
\rho^2}\,\epsilon\) with \(\epsilon
\sim N(0, I)\), marginal variance held at 1), while the sampler
uses the innovation-variance form (\(U_t = \rho\,U_{t-1} + \epsilon\) with \(\epsilon \sim N(0, \sigma_{uv}^2)\)); the
two coincide when \(\sigma_{uv}^2 = 1 -
\rho^2\), so don’t mix them. The data-generating \(G\) is constant across periods, matching
dynamic_G = FALSE; the time-varying alternative appears
further down.
n_periods <- 5
n_users <- 20
n_items <- 15
true_rho <- 0.9
G_long <- matrix(c(1, 0.3, 0.3, -0.8), 2, 2)
U_t <- vector("list", n_periods)
V_t <- vector("list", n_periods)
U_t[[1]] <- matrix(rnorm(n_users * 2), n_users, 2)
V_t[[1]] <- matrix(rnorm(n_items * 2), n_items, 2)
for(t in 2:n_periods) {
U_t[[t]] <- true_rho * U_t[[t-1]] +
sqrt(1 - true_rho^2) * matrix(rnorm(n_users * 2), n_users, 2)
V_t[[t]] <- true_rho * V_t[[t-1]] +
sqrt(1 - true_rho^2) * matrix(rnorm(n_items * 2), n_items, 2)
}
Y_list <- list()
for(t in 1:n_periods) {
eta_t <- U_t[[t]] %*% G_long %*% t(V_t[[t]])
# continuous-valued interactions (normal family): a numerically stable
# outcome for demonstrating the dynamic bipartite sampler end to end.
Y_list[[t]] <- eta_t +
matrix(rnorm(n_users * n_items, 0, 0.5), n_users, n_items)
rownames(Y_list[[t]]) <- paste0("User", 1:n_users)
colnames(Y_list[[t]]) <- paste0("Item", 1:n_items)
}
names(Y_list) <- paste0("T", 1:n_periods)
cat("Time periods:", length(Y_list), "\n")
#> Time periods: 5
cat("Dimensions per period:", dim(Y_list[[1]]), "\n")
#> Dimensions per period: 20 15
cat("Average edge weight:", round(mean(sapply(Y_list, mean)), 2), "\n")
#> Average edge weight: -0.01
cat("True latent AR(1) coefficient:", true_rho, "\n")
#> True latent AR(1) coefficient: 0.9
Static Model
The simplest longitudinal model holds latent positions and additive
effects constant over time, pooling all periods into a single, more
precise estimate per actor.
# iterations are small so the vignette builds quickly; use burn >= 1000
# and nscan >= 5000 for real analyses.
fit_static <- lame(
Y = Y_list,
mode = "bipartite",
R_row = 2,
R_col = 2,
family = "normal",
dynamic_uv = FALSE,
dynamic_ab = FALSE,
burn = 100,
nscan = 500,
odens = 5,
verbose = FALSE,
plot = FALSE
)
summary(fit_static)
#>
#> === Longitudinal AME Model Summary ===
#>
#> Call:
#> [1] "Y ~ a[i] + b[j] + U[i,1:2] %*% G %*% V[j,1:2]', family = 'normal'"
#>
#> Time periods: 5
#> Family: normal
#> Mode: bipartite
#>
#> Note: STATIC fit pooled across 5 time periods --
#> U, V, a, b are time-invariant; per-period predictions vary
#> only through per-period covariates. For time-varying effects,
#> refit with dynamic_uv = TRUE and/or dynamic_ab = TRUE.
#>
#> Regression coefficients:
#> ------------------------
#> Estimate StdError z_value p_value CI_lower CI_upper
#> intercept 0.007 0.024 0.277 0.782 -0.036 0.057
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Note: stars are a visual hint from posterior mean / SD only; for inference use the credible intervals.
#>
#> Variance components:
#> -------------------
#> Estimate StdError
#> va 0.018 0.007
#> cab 0.000 0.000
#> vb 0.019 0.008
#> rho 0.000 0.000
#> ve 0.897 0.032
#> (va = sender, cab = sender-receiver covariance, vb = receiver,
#> rho = dyadic correlation, ve = residual variance)
#> Note: bipartite model (rho fixed to 0, cab fixed to 0)
Dynamic Model
When the latent structure changes over time (tastes shift, items gain
or lose popularity), the dynamic model lets latent positions and
additive effects drift via AR(1) processes: \(\rho\) close to 1 means slow, gradual
change; close to 0, near-independent structure from period to
period.
# same reduced iterations as above for vignette speed.
fit_dynamic <- lame(
Y = Y_list,
mode = "bipartite",
R_row = 2,
R_col = 2,
family = "normal",
dynamic_uv = TRUE,
dynamic_ab = TRUE,
burn = 100,
nscan = 500,
odens = 5,
verbose = FALSE,
plot = FALSE
)
summary(fit_dynamic)
#>
#> === Longitudinal AME Model Summary ===
#>
#> Call:
#> [1] "Y ~ a[i] + b[j] + U[i,1:2] %*% G %*% V[j,1:2]', family = 'normal'"
#>
#> Time periods: 5
#> Family: normal
#> Mode: bipartite
#> Dynamic latent positions: enabled (rho_uv = 0.97 )
#> Dynamic additive effects: enabled (rho_ab = 0.389 )
#>
#> Regression coefficients:
#> ------------------------
#> Estimate StdError z_value p_value CI_lower CI_upper
#> intercept -0.001 0.017 -0.035 0.972 -0.037 0.029
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Note: stars are a visual hint from posterior mean / SD only; for inference use the credible intervals.
#>
#> Variance components:
#> -------------------
#> Estimate StdError
#> va 0.022 0.005
#> cab 0.000 0.000
#> vb 0.017 0.004
#> rho 0.000 0.000
#> ve 0.287 0.021
#> (va = sender, cab = sender-receiver covariance, vb = receiver,
#> rho = dyadic correlation, ve = residual variance)
#> Note: bipartite model (rho fixed to 0, cab fixed to 0)
The dynamic model reports estimated AR(1) persistence parameters for
the latent positions (rho_uv) and additive effects
(rho_ab).
cat(sprintf("True rho_uv = %.2f | posterior mean = %.3f | 95%% CI = [%.3f, %.3f]\n",
true_rho,
mean(fit_dynamic$rho_uv),
quantile(fit_dynamic$rho_uv, 0.025),
quantile(fit_dynamic$rho_uv, 0.975)))
#> True rho_uv = 0.90 | posterior mean = 0.970 | 95% CI = [0.949, 0.988]
rho_uv answers whether the latent structure
persists: a posterior near 1 justifies pooling across time and
forecasting from the trajectory; a posterior near 0 says positions
refresh each period and a dynamic model adds little over separate
per-period fits. Here the posterior concentrates firmly at the high end,
correctly identifying the strong persistence we simulated. One
calibration note specific to bipartite panels: the likelihood identifies
only the product \(U G V'\), so the
innovation variance that pins down the exact value of \(\rho_{uv}\) is weakly determined and the
estimate is pushed toward the top of its range. Use it to answer
“persistent or not” – which it does reliably – rather than to
distinguish, say, 0.90 from 0.95, and corroborate it with the trajectory
plot and the rotation-drift diagnostic introduced below.
rho_ab, the persistence of the additive sender/receiver
effects, comes back well below rho_uv here – expected, not
a warning sign. Our simulation put no persistent sender or receiver
activity into the data (the outcome is \(U_t G
V_t'\) plus noise), so va and vb
are small and rho_ab is estimated from what is essentially
noise. When your data do carry real period-to-period activity levels,
rho_ab is read exactly as rho_uv; in general,
interpret each persistence parameter alongside its own variance
component rather than comparing the two to each other.
Time-varying Interaction Matrix (dynamic_G)
When the way the row and column latent spaces interact is itself
shifting over time, not just the actor positions, set
dynamic_G = TRUE. By default, lame() runs a
Carter-Kohn FFBS on vec(G_t) under an AR(1) state-space
prior (persistence and innovation variance sampled jointly) and attaches
three outputs:
fit$G_cube – the final draw’s per-period \(G_t\) (\(R_\text{row} \times R_\text{col} \times
T\)).fit$G_cube_post_mean – the posterior-mean \(G_t\) averaged across all stored draws,
same shape; this is the recommended summary.fit$G_cube_post_sd – the elementwise posterior SD
across draws, same shape.
The hyperparameter chains live at fit$RHO_G and
fit$SIGMA_G2. A rotation-drift diagnostic
(fit$G_rotation_drift) compares the variance of the raw
\(G_t\) entries to their
SVD-canonicalised counterparts; when the ratio is large the apparent
time variation is dominated by latent-rotation drift rather than real
change in the interaction matrix.
With method = "als", dynamic_G = TRUE
estimates the same per-period surface as a penalized point path for
normal, binary, and Poisson panels, aligning named changing-composition
panels to the union actor sets. The object still stores
G_cube, G_cube_post_mean (the point path), and
G_cube_post_sd (zero – ALS has no posterior draws);
smoothing settings are reported through rho_G,
RHO_G, SIGMA_G2, and
lambda_G_als. Use MCMC when you need posterior draws or the
rank/censored families.
fit_dynG <- lame(
Y = Y_list,
mode = "bipartite",
R_row = 2,
R_col = 2,
family = "normal",
dynamic_uv = TRUE,
dynamic_ab = TRUE,
dynamic_G = TRUE,
burn = 100,
nscan = 500,
odens = 5,
verbose = FALSE,
plot = FALSE
)
#> Warning: `dynamic_G` rotation-drift diagnostic flagged the fit (ratio = 6.63).
#> ℹ Apparent G_t variation is dominated by U/V rotation drift, not real temporal
#> change.
#> ℹ Use `fit$G_cube_post_mean` (canonical reporting) rather than per-draw
#> `fit$G_cube`.
# fit$G_cube_post_mean has dimensions R_row x R_col x n_periods
if (!is.null(fit_dynG$G_cube_post_mean)) {
cat("G_cube_post_mean dimensions:", dim(fit_dynG$G_cube_post_mean), "\n")
cat("Frobenius norm of posterior-mean G_t per period:\n")
print(round(apply(fit_dynG$G_cube_post_mean, 3,
function(g) sqrt(sum(g^2))), 3))
cat("Rotation-drift ratio (>= 5 flags latent-rotation domination):",
round(fit_dynG$G_rotation_drift$ratio, 2), "\n")
} else {
cat("G_cube was not attached (no bipartite + RA>0 + RB>0 case).\n")
}
#> G_cube_post_mean dimensions: 2 2 5
#> Frobenius norm of posterior-mean G_t per period:
#> [1] 6.920 7.005 5.607 4.697 4.341
#> Rotation-drift ratio (>= 5 flags latent-rotation domination): 6.63
The Frobenius norms do drift – 6.920, 7.005, 5.607, 4.697, 4.341,
down roughly a third from first period to last – which, read alone,
would look like a genuinely time-varying interaction. They are not: the
rotation-drift ratio is 6.63, above the ~5 threshold, and
lame() printed the warning shown above. Only the product
\(U_t G_t V_t'\) is identified, so
rotation and scale can slosh between the latent positions and \(G_t\) without changing the fit; we
simulated \(G\) as constant, and that
sloshing is all the drifting norms pick up. The rule: when the ratio
exceeds ~5, do not treat norm drift as evidence of temporal change;
report G_cube_post_mean rather than the per-draw
G_cube, and keep dynamic_G = TRUE only when
the norms drift and the ratio stays below ~5 and the
static-\(G\) fit shows worse GOF
coverage.
Visualizing Temporal Evolution
With dynamic effects, uv_plot can show trajectories:
each line traces one actor’s latent position from the first to the last
period. With 35 actors a per-actor legend and labels would swamp the
panel, so we suppress both and read the plot for the overall pattern:
positions drift gradually and stay clustered rather than scattering.
uv_plot(fit_dynamic, plot_type = "trajectory", label.nodes = FALSE) +
guides(color = "none") +
ggtitle("Dynamic latent trajectories")
To inspect specific actors, pass their names to
highlight =: the named actors are coloured (and identified
in the legend) while every other trajectory is greyed out.
uv_plot(fit_dynamic, plot_type = "trajectory",
highlight = c("User1", "User5", "Item3"),
label.nodes = FALSE) +
ggtitle("Highlighted actor trajectories")
Convergence Diagnostics
Dynamic models add parameters (\(\rho\), innovation variances) that need
adequate samples to be well-estimated. We exclude the dyadic correlation
and sender-receiver covariance from the trace panels because they are
structural constants in bipartite mode (fixed at zero; for a binary
bipartite fit the probit error variance is likewise fixed at 1), so the
diagnostics show only parameters that are actually sampled.
trace_plot(fit_dynamic,
exclude = c("Dyadic Correlation", "Sender-Receiver Covariance"))
Look for traces that mix well (no long flat stretches or slow drifts)
and densities that are smooth and unimodal. The short chains used here
don’t guarantee convergence; run longer chains in practice. For
numerical diagnostics – posterior::summarise_draws() for
split-\(\hat R\) and ESS,
ame_parallel(..., n_chains = 4, combine_method = "pool")
for between-chain \(\hat R\),
prior_summary(), and loo::loo() after
refitting with save_log_lik = TRUE – see the worked
end-to-end stack in the cross-sectional
vignette.
Comparing Static and Dynamic Models
Is the dynamic model worth the added complexity? The goodness-of-fit
plots provide one comparison: if the dynamic model’s
posterior-predictive distributions better cover the observed statistics,
the extra flexibility is paying off. In each panel the observed
statistic at each period is a solid Okabe-Ito orange
(#D55E00) line with filled points, the
posterior-predictive median a dashed dark grey line,
and the grey ribbon the 95% posterior-predictive credible
interval – dual-encoded on colour and linetype so the
comparison stays legible in greyscale and for colour-blind readers. We
restrict the panels to the two degree-heterogeneity statistics: “Sender
Degree Heterogeneity” (row-mean SD) and “Receiver Degree Heterogeneity”
(column-mean SD). The third bipartite GOF statistic, the four-cycle
census, binarizes at \(|y| > 0\), so
on a continuous (normal-family) network every cell is an edge and the
count is the constant \(\binom{n_A}{2}\binom{n_B}{2} = 19950\) for
the data and every posterior draw; it is informative only for families
with structural zeros, such as the binary network in the cross-sectional
section above.
# gof plots for each model; four-cycles is omitted because it is
# constant on this dense continuous network (see text above)
gof_plot(fit_static, statistics = c("sd.row", "sd.col"))
gof_plot(fit_dynamic, statistics = c("sd.row", "sd.col"))
We can also compare the two informative statistics directly: the SD
of row means (how much users vary in activity) and the SD of column
means (how much items vary in popularity).
gof_static <- fit_static$GOF
gof_dynamic <- fit_dynamic$GOF
# extract posterior predictive samples (exclude column 1, which is observed)
# each element is a matrix [n_time x n_mcmc], so colMeans averages across time
extract_ppc <- function(gof, stat) {
mat <- gof[[stat]]
colMeans(mat[, -1, drop = FALSE])
}
# observed values (column 1, averaged across time periods)
obs_vals <- c(
mean(gof_static$sd.rowmean[, 1]),
mean(gof_static$sd.colmean[, 1])
)
gof_df <- data.frame(
value = c(
extract_ppc(gof_static, "sd.rowmean"),
extract_ppc(gof_dynamic, "sd.rowmean"),
extract_ppc(gof_static, "sd.colmean"),
extract_ppc(gof_dynamic, "sd.colmean")
),
model = rep(rep(c("Static", "Dynamic"),
each = ncol(gof_static$sd.rowmean) - 1), 2),
statistic = rep(c("Row Mean SD", "Column Mean SD"),
each = 2 * (ncol(gof_static$sd.rowmean) - 1))
)
obs_df <- data.frame(
statistic = c("Row Mean SD", "Column Mean SD"),
observed = obs_vals
)
ggplot(gof_df, aes(x = model, y = value)) +
geom_boxplot() +
geom_hline(data = obs_df, aes(yintercept = observed),
linetype = 2) +
facet_wrap(~statistic, scales = "free_y") +
labs(title = "GOF Comparison: Static vs Dynamic",
subtitle = "Dashed Line = Observed Value",
x = "Model", y = "Simulated Statistic") +
theme_bw() +
theme(
panel.border = element_blank(),
axis.ticks = element_blank(),
legend.position = "top",
strip.background = element_rect(fill = "black", color = "black"),
strip.text = element_text(color = "white", hjust = 0)
)
The boxplot alone makes the two models look closer than they are. In
gof_plot(fit_static) the posterior-predictive median is
essentially flat across periods, missing the observed line as it swings
(observed sender heterogeneity runs 0.35, 0.49, 0.18, 0.36, 0.37); in
gof_plot(fit_dynamic) the median tracks that zigzag
closely. The residual variance agrees: ve falls from 0.897
in the static fit to 0.287 in the dynamic fit, close to the 0.25 noise
floor we simulated (noise sd 0.5). The boxplot looks more even because
colMeans averages each statistic across the five periods
within each draw, washing out exactly the period-to-period tracking that
is the dynamic model’s advantage. Even so, the dynamic boxes sit closer
to the dashed observed line on both statistics (on Column Mean SD,
roughly 0.19 versus 0.22 against an observed value of about 0.14), and
neither box covers it – both models over-predict time-averaged degree
heterogeneity, the dynamic one less so. Time-averaged degree statistics
hide temporal fit, so pair the boxplot with ve, the
per-period gof_plot(), and loo::loo_compare()
before concluding a dynamic fit is not earning its complexity.
Choosing Latent Dimensions
The latent dimensions (R_row, R_col)
control how rich the multiplicative structure is; a practical way to
choose them is to fit at several ranks and compare GOF statistics.
# include R=0 as a no-latent-space baseline
dims_to_test <- list(
c(0, 0),
c(1, 1),
c(2, 2)
)
# compute observed GOF statistics from the data
obs_gof <- gof_stats(Y_bipartite, mode = "bipartite")
gof_results <- list()
for(i in seq_along(dims_to_test)) {
fit_temp <- ame(
Y = Y_bipartite,
mode = "bipartite",
R_row = dims_to_test[[i]][1],
R_col = dims_to_test[[i]][2],
family = "binary",
burn = 200,
nscan = 2500,
odens = 5,
verbose = FALSE
)
# posterior predictive tail probability: proportion of simulated
# statistics at or above the observed value (one-sided)
gof_results[[i]] <- c(
R_row = dims_to_test[[i]][1],
R_col = dims_to_test[[i]][2],
pval_rowmean = mean(fit_temp$GOF[, "sd.rowmean"] >= obs_gof["sd.rowmean"]),
pval_colmean = mean(fit_temp$GOF[, "sd.colmean"] >= obs_gof["sd.colmean"]),
pval_fourcycles = mean(fit_temp$GOF[, "four.cycles"] >= obs_gof["four.cycles"])
)
}
do.call(rbind, gof_results)
These are one-sided posterior predictive tail probabilities, not
formal test p-values: values near 0.5 mean the observed statistic lies
in the bulk of the predictive distribution, while values near 0 or 1
flag systematic over- or under-prediction. For ame(),
fit$GOF stores the observed statistic in row 1, so the
calculation includes that one point. The relationship to R
is not monotonic, so choose the smallest rank at which the important
statistics fall away from the endpoints. If degree statistics remain at
an endpoint across ranks while four-cycle statistics are well
reproduced, adding dimensions is not addressing the source of the
mismatch.
Practical Guidance
When to use bipartite models. Use bipartite mode
whenever ties only form between two fundamentally different types of
nodes. A rectangular adjacency matrix almost certainly needs it, and
even a square one (say, 20 students and 20 courses) does when rows and
columns represent different kinds of entities.
Interpreting the output. The additive effects
(APM, BPM) are the most directly
interpretable: which row nodes are unusually active and which column
nodes unusually popular, after controlling for covariates and latent
structure. The latent positions (U, V) and
interaction matrix (G) capture residual association,
answering “beyond the covariates and overall activity levels, which row
nodes tend to connect to which column nodes?”