--- title: "Bipartite Network Analysis" author: "Cassy Dorff, Shahryar Minhas, and Tosin Salau" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Bipartite Network Analysis} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, dpi = 96, out.width = "100%", fig.align = "center" ) ``` ```{r setup} library(lame) library(ggplot2) set.seed(6886) ``` ## 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](cross_sec_ame.html)): $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: | Family | Bipartite supported? | Notes | |-----------|----------------------|-------| | `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()` ### Coming from long-format data? If your data lives in a long-format table (`patient_id`, `drug`, `prescribed`), build a bipartite `netify` object and pass it directly to `ame()`: ```{r tibble_bridge, eval = FALSE} rx_net <- netify::netify( rx, actor1 = "patient_id", actor2 = "drug", weight = "prescribed", mode = "bipartite", symmetric = FALSE, missing_to_zero = FALSE ) fit <- ame(rx_net, family = "binary", mode = "bipartite", R = 2) ``` Row and column names carry through to `fit$APM` (one number per patient) and `fit$BPM` (one number per drug) downstream. Use `missing_to_zero = TRUE` only when absent patient-drug rows are observed non-prescriptions rather than unobserved dyads. ### 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. ```{r simulate_cross_sectional} # 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") cat("Enrollment rate:", round(mean(Y_bipartite), 2), "\n") ``` ### 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`. ```{r fit_cross_sectional, message=FALSE, warning=FALSE} # 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) ``` 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. ```{r visualize_cross_sectional, fig.height=5, fig.alt="Bipartite latent-space biplot placing students and courses in the same 2D coordinate system, with proximity indicating a higher predicted enrollment probability."} 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: ```{r interaction_heatmap, fig.height=6, fig.alt="Heatmap of the fitted multiplicative structure U G V' across student-course pairs using a diverging colour-blind-safe ColorBrewer RdBu palette with the scale capped at the 99th percentile of absolute fitted values: red tiles mark positive latent affinity (predicted enrollment above baseline), blue tiles mark negative latent affinity (predicted below baseline), and near-white tiles mark a posterior estimate near zero."} 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. ```{r simulate_longitudinal} 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") cat("Dimensions per period:", dim(Y_list[[1]]), "\n") cat("Average edge weight:", round(mean(sapply(Y_list, mean)), 2), "\n") cat("True latent AR(1) coefficient:", true_rho, "\n") ``` ### 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. ```{r fit_longitudinal_static, message=FALSE} # 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) ``` ### 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. ```{r fit_longitudinal_dynamic, message=FALSE} # 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) ``` The dynamic model reports estimated AR(1) persistence parameters for the latent positions (`rho_uv`) and additive effects (`rho_ab`). ```{r dynamic_rho_recovery} 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))) ``` `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. ```{r fit_dynamic_G, message=FALSE} 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 ) # 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") } ``` 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. ```{r visualize_longitudinal, fig.height=5, fig.alt="Trajectory plot tracing each actor's latent position from the first to the last period; paths 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. ```{r visualize_longitudinal_highlight, fig.height=5, fig.alt="Trajectory plot with three highlighted actors (User1, User5, Item3) coloured and named in the legend while all other trajectories are greyed out for context."} 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. ```{r trace_dynamic, fig.height=6, fig.alt="MCMC trace plots and posterior densities for the sampled dynamic bipartite model parameters (intercept and variance components); well-mixed caterpillar traces indicate adequate sampling."} 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](cross_sec_ame.html). ## 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. ```{r model_comparison, fig.height=4, fig.alt="Side-by-side goodness-of-fit panels comparing posterior predictive coverage of bipartite sender and receiver degree heterogeneity for the static versus dynamic model."} # 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). ```{r model_comparison_boxplot, fig.height=4, fig.alt="Boxplots of posterior predictive bipartite statistics (row-mean SD and column-mean SD) faceted by static and dynamic models, with horizontal reference lines at observed values."} 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. ```{r dimension_selection, eval=FALSE} # 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?" ## Extracting Positions for Custom Analysis If you need the latent positions in a tidy format, `latent_positions()` returns a data frame with one row per actor-dimension-time combination. The `posterior_sd` column is populated here because `fit_cross` stored thinned U/V draws via `posterior_opts = posterior_options(save_UV = TRUE)`; without that option the column is `NA` and a refit hint is printed. ```{r} lp <- latent_positions(fit_cross) head(lp) # filter to just the row nodes (students) lp_students <- lp[lp$type == "U", ] cat("Student positions:", nrow(lp_students), "rows\n") ``` This is especially useful for bipartite networks where U (row nodes) and V (column nodes) have different numbers of actors.