Change co-occurrence convergence analysis to analysis of binary Markov chain.

Update summary.bgm and print.summary.bgm accordingly.

Author: MaartenMarsman

R/bgms-methods.R

@@ -70,8 +70,9 @@ summary.bgms <- function(object, ...) {
     }
 
     if (!is.null(object$posterior_summary_pairwise_allocations)) {
-      out$pairwise_allocations <- object$posterior_summary_pairwise_allocations
-      out$allocations_mean <- object$posterior_mean_allocations
+      out$allocations <- object$posterior_summary_pairwise_allocations
+      out$mean_allocations <- object$posterior_mean_allocations
+      out$mode_allocations <- object$posterior_mode_allocations
       out$num_blocks <- object$posterior_num_blocks
     }
 
@@ -118,23 +119,24 @@ print.summary.bgms <- function(x, digits = 3, ...) {
     ind <- head(x$indicator, 6)
     ind[] <- lapply(ind, function(col) ifelse(is.na(col), "", round(col, digits)))
     print(ind)
-    if (nrow(x$indicator) > 6)
-      cat("... (use `summary(fit)$indicator` to see full output)\n")
+    if (nrow(x$indicator) > 6) cat("... (use `summary(fit)$indicator` to see full output)\n")
     cat("Note: NA values are suppressed in the print table. They occur when an indicator\n")
     cat("was constant (all 0 or all 1) across all iterations, so sd/mcse/n_eff/Rhat\n")
     cat("are undefined; `summary(fit)$indicator` still contains the NA values.\n\n")
   }
 
-  if (!is.null(x$pairwise_allocations)) {
+  if (!is.null(x$allocations)) {
     cat("Pairwise node co-clustering proportion:\n")
-    print(round(head(x$pairwise_allocations, 6), digits = digits))
-    if (nrow(x$pairwise_allocations) > 6) cat("... (use `summary(fit)$allocations` to see full output)\n")
+    print(round(head(x$allocations, 6), digits = digits))
+    if (nrow(x$allocations) > 6) cat("... (use `summary(fit)$allocations` to see full output)\n")
     cat("\n")
   }
 
-  if (!is.null(x$allocations_mean)) {
-    cat("Mean posterior node allocation vector :\n")
-    print(round(head(x$allocations_mean, 6), digits = digits))
+  if (!is.null(x$mean_allocations)) {
+    cat("Mean posterior node allocation vector:\n")
+    print(round(head(x$mean_allocations, 6), digits = digits))
+    cat("Mode posterior node allocation vector:\n")
+    print(round(head(x$mode_allocations, 6), digits = digits))
     cat("\n")
   }
 
@@ -185,7 +187,6 @@ coef.bgms <- function(object, ...) {
 }
 
 
-
 .warning_issued <- FALSE
 warning_once <- function(msg) {
   if (!.warning_issued) {

R/mcmc_summary.R

@@ -249,8 +249,8 @@ summarize_alloc_pairs = function(allocations, node_names = NULL) {
   Pairs  = t(combn(seq_len(no_variables), 2))
   nparam = nrow(Pairs)
 
-  result = matrix(NA, nparam, 5)
-  colnames(result) = c("mean", "sd", "mcse", "n_eff", "Rhat")
+  result = matrix(NA, nparam, 9)
+  colnames(result) = c("mean", "sd", "mcse", "n0->0", "n0->1", "n1->0", "n1->1", "n_eff", "Rhat")
 
   # helper to construct a "time-series"
   get_draws_pair = function(i, j) {
@@ -265,16 +265,32 @@ summarize_alloc_pairs = function(allocations, node_names = NULL) {
   for (p in seq_len(nparam)) {
     i = Pairs[p, 1]; j = Pairs[p, 2]
     draws = get_draws_pair(i, j)
-    vec   = as.vector(draws)
-    phat  = mean(vec)
-    sdev  = sd(vec)
 
-    est   = compute_rhat_ess(draws)
-    n_eff = as.numeric(est$ess)
-    Rhat  = as.numeric(est$rhat)
+    vec = as.vector(draws)
+    T = length(vec)
+    g_next = vec[-1]
+    g_curr = vec[-T]
+
+    p_hat = mean(vec)
+    sd = sqrt(p_hat * (1 - p_hat))
+    n00 = sum(g_curr == 0 & g_next == 0)
+    n01 = sum(g_curr == 0 & g_next == 1)
+    n10 = sum(g_curr == 1 & g_next == 0)
+    n11 = sum(g_curr == 1 & g_next == 1)
+
+    if (n01 + n10 == 0) {
+      n_eff = mcse = R = NA_real_
+    } else {
+      a = n01 / (n00 + n01)
+      b = n10 / (n10 + n11)
+      tau_int = (2 - (a + b)) / (a + b)
+      n_eff = T / tau_int
+      mcse = sd / sqrt(n_eff)
+      est = compute_rhat_ess(draws)
+      R = est$rhat
+    }
 
-    mcse  = if (is.finite(n_eff) && n_eff > 0) sdev / sqrt(n_eff) else NA
-    result[p, ] = c(phat, sdev, mcse, n_eff, Rhat)
+    result[p, ] = c(p_hat, sd, mcse, n00, n01, n10, n11, n_eff, R)
   }
   if (is.null(node_names)) {
     rn = paste0(Pairs[,1], "-", Pairs[,2])

R/output_utils.R

@@ -136,7 +136,7 @@ prepare_output_bgm = function(
         allocations = lapply(out, `[[`, "allocations"),
         node_names = data_columnnames
       )
-      results$posterior_coclustering_matrix = sbm_convergence$co_occur_matrix
+      results$posterior_mean_coclustering_matrix = sbm_convergence$co_occur_matrix
       # calculate the estimated clustering and block probabilities
       sbm_summary = posterior_summary_SBM(allocations = lapply(out, `[[`, "allocations"),
                                           arguments = arguments)  # check if only arguments would work