From af09d22b39b483a5976c6180d294930f5ccb3412 Mon Sep 17 00:00:00 2001 From: Alex Gallinat Date: Sat, 3 Oct 2026 11:52:28 -0400 Subject: [PATCH] Fix silently wrong participation coefficients in calculate_roles() What looked like a coverage gap -- nodes missing from the roles table -- was also a correctness bug affecting the nodes that remained. With communities = NULL, calculate_roles() delegated to find_modules() at its default min_size = 3, which discards small modules. A node whose module was discarded is absent from module_table, so membership[node] is NA. In the participation-coefficient loop, table() silently drops those NAs from a neighbour's module tally, but the denominator was still the node's full degree. The fractions therefore summed to less than 1 and P = 1 - sum(...) came out too high for any *retained* node with a discarded neighbour. Measured on a 150-node sample_gnp with cluster.method = "walktrap": 29 nodes absent, and of the 121 reported, 18 had a wrong P and 7 a wrong role (e.g. P 0.75 instead of 0.50, R3 instead of R2). Two fixes: - find_modules() is now called with min_size = 1. Roles are a per-node measure, so filtering modules by size here was never meaningful, and it is what introduced the NAs. Every vertex now receives a role. - The denominator is now the number of neighbours actually tallied rather than the full degree. These agree whenever every neighbour has a module, and the former stays well defined when some do not, so a partial membership yields a defined P instead of an inflated one. Partial coverage can still arise legitimately: spinglass cannot run on a disconnected graph, so find_modules() restricts it to the largest connected component. That case now warns with a count and example vertices instead of returning a short table silently. Two further latent defects on the user-supplied membership path, both found while testing the above: - A membership vector containing NAs errored with "Invalid vertex names": NA was treated as a module and induced_subgraph() was handed NA node names. NA modules are now skipped, and those vertices get role = NA. - An *unnamed* membership vector of the correct length was accepted, then had index names ("1", "2", ...) invented for it, so every subsequent vertex lookup failed unless the graph's own names happened to match. Unnamed vectors are now treated as positional, in vertex order, which is what igraph::membership() returns. Tests go from 371 to 386, including a regression test that compares the find_modules() path against the equivalent explicit membership and asserts identical P and roles. Co-Authored-By: Claude Opus 5 (1M context) --- R/calculate_roles.R | 72 ++++++++++++++++++---- man/calculate_roles.Rd | 15 ++++- tests/testthat/test-calculate_roles.R | 86 ++++++++++++++++++++++++--- 3 files changed, 151 insertions(+), 22 deletions(-) diff --git a/R/calculate_roles.R b/R/calculate_roles.R index ff58376..2640a1e 100644 --- a/R/calculate_roles.R +++ b/R/calculate_roles.R @@ -21,10 +21,23 @@ #' \item{`plot`}{A `ggplot2` object, or `NULL` when `plot = FALSE`. The element is #' always present, so the return shape does not depend on the arguments.} #' \item{`roles_definitions`}{A data frame describing the seven role types and their conditions.} -#' \item{`result`}{A data frame with node-level information: node name, module, z-score, participation coefficient, and assigned role.} +#' \item{`result`}{A data frame with node-level information: node name, module, +#' z-score, participation coefficient, and assigned role. It has one row per +#' graph vertex. The exception is `cluster.method = "spinglass"`, which can only +#' be run on the largest connected component of a disconnected graph; vertices +#' outside it have no module and are absent, which raises a warning. `z`, `p` and +#' `role` are `NA` for any vertex whose module, or whose neighbours' modules, are +#' unknown.} #' } #' #' @details +#' When `communities` is `NULL`, community detection is delegated to +#' [find_modules()] with `min_size = 1`, so no module is discarded for being small +#' and every vertex receives a role. This matters for correctness as well as +#' coverage: the participation coefficient of a node is computed from the module +#' memberships of its neighbours, so dropping a neighbour's module silently distorts +#' the coefficient of the node that remains. +#' #' The node roles are defined as: #' #' \tabular{ll}{ @@ -69,7 +82,15 @@ calculate_roles <- function(graph, # Extract membership vector if (is.null(communities)) { - modules <- find_modules(graph, method = cluster.method, plot = FALSE, return_subgraphs = FALSE) + # min_size = 1 is deliberate. find_modules() defaults to min_size = 3, which + # discards small modules -- and a node whose module was discarded is absent + # from module_table, so its membership is NA. That silently corrupted the + # participation coefficient of *retained* neighbours, because the NA was + # dropped from the neighbour tally while the full degree was still used as the + # denominator. Roles are a per-node measure, so there is no reason to filter + # modules by size here at all. + modules <- find_modules(graph, method = cluster.method, min_size = 1, + plot = FALSE, return_subgraphs = FALSE) membership <- stats::setNames(modules$module_table$module, modules$module_table$node) } else if ("communities" %in% class(communities)) { @@ -102,12 +123,27 @@ calculate_roles <- function(graph, )) - vnames <- names(membership) - if (is.null(vnames)) { - vnames <- as.character(seq_len(length(membership))) + # An unnamed membership vector is positional, in vertex order, which is what + # igraph::membership() returns. Naming it from the graph is required: inventing + # index names ("1", "2", ...) instead meant every later vertex lookup failed with + # "Invalid vertex names" whenever the graph's own names were anything else. + if (is.null(names(membership))) { + names(membership) <- igraph::V(graph)$name } + vnames <- names(membership) - degrees <- igraph::degree(graph, mode = "all")[vnames] + # Warn rather than silently return a partial table. With communities = NULL this + # should now cover every vertex; the exception is cluster.method = "spinglass", + # which find_modules() can only run on the largest connected component. + missing_nodes <- setdiff(igraph::V(graph)$name, vnames) + if (length(missing_nodes) > 0) { + warning(sprintf( + "%d of %d vertices have no module assignment and are absent from the result: %s%s", + length(missing_nodes), igraph::vcount(graph), + paste(utils::head(missing_nodes, 5), collapse = ", "), + if (length(missing_nodes) > 5) ", ..." else "" + )) + } # Initialize roles dataframe roles_df <- tibble::tibble( @@ -118,11 +154,15 @@ calculate_roles <- function(graph, role = NA_character_ ) - # Compute within-module z-score - for (mod in unique(roles_df$module)) { - mod_nodes <- roles_df$node[roles_df$module == mod] + # Compute within-module z-score. NA modules are skipped rather than treated as a + # module of their own: a membership vector containing NAs would otherwise select + # NA node names and make induced_subgraph() fail with "Invalid vertex names". + # Those nodes keep z = NA and so are classified with role = NA. + for (mod in unique(roles_df$module[!is.na(roles_df$module)])) { + idx <- which(roles_df$module == mod) + mod_nodes <- roles_df$node[idx] if (length(mod_nodes) <= 1) { - roles_df$z[roles_df$node %in% mod_nodes] <- 0 + roles_df$z[idx] <- 0 next } @@ -145,11 +185,17 @@ calculate_roles <- function(graph, } neighbor_names <- igraph::V(graph)$name[nbrs] - k_i <- degrees[node] neighbor_modules <- membership[neighbor_names] k_i_m <- table(neighbor_modules) - sum_frac_sq <- sum((k_i_m / k_i)^2) - roles_df$p[i] <- 1 - sum_frac_sq + + # The denominator must be the number of neighbours actually tallied, not the + # node's full degree: table() drops neighbours with NA membership, so using + # the full degree makes the fractions sum to less than 1 and inflates P. + # These agree whenever every neighbour has a module, and stay well defined + # when some do not. + k_i <- sum(k_i_m) + + roles_df$p[i] <- if (k_i == 0) NA_real_ else 1 - sum((k_i_m / k_i)^2) } # Classify roles based on z and p diff --git a/man/calculate_roles.Rd b/man/calculate_roles.Rd index c6e6740..3ba6b1e 100644 --- a/man/calculate_roles.Rd +++ b/man/calculate_roles.Rd @@ -38,7 +38,13 @@ A list with three elements: \item{\code{plot}}{A \code{ggplot2} object, or \code{NULL} when \code{plot = FALSE}. The element is always present, so the return shape does not depend on the arguments.} \item{\code{roles_definitions}}{A data frame describing the seven role types and their conditions.} -\item{\code{result}}{A data frame with node-level information: node name, module, z-score, participation coefficient, and assigned role.} +\item{\code{result}}{A data frame with node-level information: node name, module, +z-score, participation coefficient, and assigned role. It has one row per +graph vertex. The exception is \code{cluster.method = "spinglass"}, which can only +be run on the largest connected component of a disconnected graph; vertices +outside it have no module and are absent, which raises a warning. \code{z}, \code{p} and +\code{role} are \code{NA} for any vertex whose module, or whose neighbours' modules, are +unknown.} } } \description{ @@ -50,6 +56,13 @@ Nodes are assigned to one of seven role categories (R1–R7) based on their loca If no community structure is provided, modules are automatically detected using the specified clustering method. The function can optionally produce a 2D role plot (z vs. P) highlighting the canonical role regions. +When \code{communities} is \code{NULL}, community detection is delegated to +\code{\link[=find_modules]{find_modules()}} with \code{min_size = 1}, so no module is discarded for being small +and every vertex receives a role. This matters for correctness as well as +coverage: the participation coefficient of a node is computed from the module +memberships of its neighbours, so dropping a neighbour's module silently distorts +the coefficient of the node that remains. + The node roles are defined as: \tabular{ll}{ diff --git a/tests/testthat/test-calculate_roles.R b/tests/testthat/test-calculate_roles.R index b864f8e..e17bd7f 100644 --- a/tests/testthat/test-calculate_roles.R +++ b/tests/testthat/test-calculate_roles.R @@ -92,13 +92,83 @@ test_that("a malformed communities argument is rejected", { "clustering object, a membership vector or NULL") }) -test_that("falling back to find_modules() can drop nodes from small modules", { - # Documented quirk rather than endorsed behaviour: with communities = NULL the - # function delegates to find_modules(), whose default min_size = 3 discards - # small modules, so those nodes never reach the roles table. - res <- calculate_roles(g, cluster.method = "louvain", plot = FALSE) +test_that("every vertex receives a role", { + # find_modules() is called with min_size = 1 so that no module is discarded. + # Previously it ran at the default min_size = 3 and nodes in small modules never + # reached the roles table. + for (m in c("louvain", "walktrap", "infomap")) { + set.seed(7) + res <- calculate_roles(g, cluster.method = m, plot = FALSE) + expect_setequal(res$result$node, igraph::V(g)$name) + expect_false(any(duplicated(res$result$node))) + } +}) + +test_that("the participation coefficient is unaffected by module-size filtering", { + # Regression guard for a silent numerical bug. When a module was discarded, its + # nodes' membership became NA; table() dropped those NAs from a neighbour's tally + # while the node's full degree was still used as the denominator, so P came out + # too high for *retained* nodes. On a 150-node graph with walktrap this gave 18 + # wrong coefficients and 7 wrong roles. + set.seed(3) + gw <- igraph::sample_gnp(150, 0.02, directed = FALSE) + igraph::V(gw)$name <- paste0("b", seq_len(150)) + + set.seed(7) + via_find_modules <- calculate_roles(gw, cluster.method = "walktrap", plot = FALSE) + set.seed(7) + via_membership <- calculate_roles( + gw, communities = igraph::membership(igraph::cluster_walktrap(gw)), plot = FALSE + ) + + expect_equal(nrow(via_find_modules$result), igraph::vcount(gw)) + + a <- via_find_modules$result[order(via_find_modules$result$node), ] + b <- via_membership$result[order(via_membership$result$node), ] + expect_equal(a$p, b$p) + expect_equal(a$role, b$role) +}) + +test_that("the participation coefficient stays within bounds for pendant nodes", { + # A degree-1 node has all its edges inside one module, so P must be exactly 0. + star <- igraph::make_star(10, mode = "undirected") + igraph::V(star)$name <- paste0("s", seq_len(10)) + res <- calculate_roles(star, communities = rep(1L, 10), plot = FALSE) - expect_lte(nrow(res$result), igraph::vcount(g)) - expect_true(all(res$result$node %in% igraph::V(g)$name)) - expect_false(any(duplicated(res$result$node))) + leaves <- res$result[res$result$node != "s1", ] + expect_true(all(leaves$p == 0)) + expect_true(all(res$result$p >= 0 & res$result$p <= 1)) +}) + +test_that("vertices with no module assignment are reported, not dropped silently", { + # spinglass cannot run on a disconnected graph, so find_modules() restricts it to + # the largest connected component. The resulting partial coverage now warns. + set.seed(42) + gd <- igraph::sample_gnp(80, 0.02, directed = FALSE) + igraph::V(gd)$name <- paste0("s", seq_len(80)) + expect_gt(igraph::components(gd)$no, 1) + + expect_warning( + res <- calculate_roles(gd, cluster.method = "spinglass", plot = FALSE), + "vertices have no module assignment" + ) + expect_lt(nrow(res$result), igraph::vcount(gd)) +}) + +test_that("a membership vector containing NAs is handled rather than erroring", { + # Previously this failed with "Invalid vertex names": the NA module was treated as + # a module, and induced_subgraph() was handed NA node names. + set.seed(1) + gn <- igraph::sample_pa(40, directed = FALSE) + igraph::V(gn)$name <- paste0("n", seq_len(40)) + m <- igraph::membership(igraph::cluster_louvain(gn)) + m[1:5] <- NA + + res <- suppressWarnings(calculate_roles(gn, communities = m, plot = FALSE)) + + expect_equal(nrow(res$result), 40L) + expect_equal(sum(is.na(res$result$module)), 5L) + # Undefined inputs propagate to an undefined role rather than a wrong one. + expect_true(all(is.na(res$result$role[is.na(res$result$module)]))) + expect_true(all(res$result$p >= 0 & res$result$p <= 1, na.rm = TRUE)) })