Assoc graph weights

NAMESPACE

@@ -136,6 +136,7 @@ importFrom(utils,data)
 importFrom(utils,menu)
 importFrom(utils,type.convert)
 importFrom(vcd,assoc)
+importFrom(vcd,assocstats)
 importFrom(vcd,is.structable)
 importFrom(vcd,labeling_border)
 importFrom(vcd,mosaic)

R/assoc_graph.R

@@ -15,16 +15,24 @@
 #' @param result Type of result to return: \code{"igraph"} (default) returns an
 #'   \code{\link[igraph:igraph-package]{igraph}} object; \code{"matrix"} returns the
 #'   adjacency matrix; \code{"edge_list"} returns a two-column character matrix of edges.
+#' @param measure Type of association measure for edge weights (only for \code{glm} method):
+#'   \code{"none"} (default) produces an unweighted graph;
+#'   \code{"chisq"} computes partial chi-squared statistics (deviance change when each
+#'   edge is removed from the model);
+#'   \code{"cramer"} computes Cramer's V from the marginal two-way table for each edge.
 #' @param \dots Additional arguments (currently unused).
 #'
 #' @return Depending on \code{result}:
 #'   \itemize{
 #'     \item \code{"igraph"}: An \code{igraph} undirected graph object of class
 #'           \code{c("assoc_graph", "igraph")}, with vertex names corresponding to
-#'           the variable names.
-#'     \item \code{"matrix"}: A symmetric adjacency matrix (0/1) with variable names as
-#'           row and column names.
-#'     \item \code{"edge_list"}: A two-column character matrix, each row an edge.
+#'           the variable names. When \code{measure != "none"}, edge weights are
+#'           stored as \code{E(g)$weight} and the measure name as \code{g$measure}.
+#'     \item \code{"matrix"}: A symmetric adjacency matrix with variable names as
+#'           row and column names. Contains 0/1 when unweighted, or association
+#'           strength values when \code{measure} is specified.
+#'     \item \code{"edge_list"}: When unweighted, a two-column character matrix (from, to).
+#'           When \code{measure} is specified, a data frame with columns from, to, and weight.
 #'   }
 #'
 #' @details
@@ -82,13 +90,33 @@
 #' mod.cond <- glm(Freq ~ (cigarette + alcohol + marijuana)^2 +
 #'                         (alcohol + sex + race)^2,
 #'                 data = DaytonSurvey, family = poisson)
+#'
+#' # define groups for the model
+#' gps <- list(c("cigarette", "marijuana"),
+#'             "alcohol",
+#'             c("sex", "race"))
+#'
 #' assoc_graph(mod.cond)
 #' plot(assoc_graph(mod.cond),
-#'      groups = list(c("cigarette", "marijuana"),
-#'                    "alcohol",
-#'                    c("sex", "race")),
-#'      main = "{Mar,Cig} indep of {Race,Sex} | Alc",
-#'      layout = igraph::layout_nicely)
+#'      groups = gps,
+#'      layout = igraph::layout_nicely,
+#'      main = "{R,S} indep {M,C} | A")
+#'
+#' # Weighted graph: partial chi-squared
+#' g <- assoc_graph(mod.cond, measure = "chisq")
+#' g
+#' plot(g, edge.label = TRUE,
+#'      groups = gps,
+#'      layout = igraph::layout_nicely,
+#'      main = "Partial chi-squared weights")
+#'
+#' # Cramer's V (marginal)
+#' g2 <- assoc_graph(mod.cond, measure = "cramer")
+#' g2
+#' plot(g2, edge.label = TRUE,
+#'      groups = gps,
+#'      layout = igraph::layout_nicely,
+#'      main = "Cramer's V weights")
 #'
 assoc_graph <- function(x, ...) {
   UseMethod("assoc_graph")
@@ -112,11 +140,60 @@ assoc_graph.loglm <- function(x, result = c("igraph", "matrix", "edge_list"), ..
 }
 
 #' @rdname assoc_graph
+#' @importFrom vcd assocstats
 #' @export
-assoc_graph.glm <- function(x, result = c("igraph", "matrix", "edge_list"), ...) {
+assoc_graph.glm <- function(x, result = c("igraph", "matrix", "edge_list"),
+                            measure = c("none", "chisq", "cramer"), ...) {
   result <- match.arg(result)
+  measure <- match.arg(measure)
   margins <- .glm_to_margins(x)
-  .margins_to_assoc_graph(margins, result = result)
+
+  if (measure == "none") {
+    return(.margins_to_assoc_graph(margins, result = result))
+  }
+
+  # Compute weights — need the igraph object regardless of result type
+
+  g <- .margins_to_assoc_graph(margins, result = "igraph")
+
+  weights <- switch(measure,
+    chisq  = .compute_chisq_weights(g, x),
+    cramer = .compute_cramer_weights(g, x)
+  )
+
+  g$measure <- measure
+
+  if (igraph::ecount(g) > 0) {
+    igraph::E(g)$weight <- weights
+  }
+
+  if (result == "igraph") {
+    return(g)
+  }
+
+  if (result == "matrix") {
+    all_vars <- igraph::V(g)$name
+    nv <- length(all_vars)
+    adj <- matrix(0, nv, nv, dimnames = list(all_vars, all_vars))
+    if (igraph::ecount(g) > 0) {
+      el <- igraph::as_edgelist(g)
+      for (i in seq_len(nrow(el))) {
+        adj[el[i, 1], el[i, 2]] <- weights[i]
+        adj[el[i, 2], el[i, 1]] <- weights[i]
+      }
+    }
+    return(adj)
+  }
+
+  if (result == "edge_list") {
+    if (igraph::ecount(g) == 0) {
+      return(matrix(character(0), ncol = 2, dimnames = list(NULL, c("from", "to"))))
+    }
+    el <- igraph::as_edgelist(g)
+    out <- data.frame(from = el[, 1], to = el[, 2], weight = weights,
+                      stringsAsFactors = FALSE)
+    return(out)
+  }
 }
 
 
@@ -175,6 +252,61 @@ assoc_graph.glm <- function(x, result = c("igraph", "matrix", "edge_list"), ...)
 }
 
 
+# --- Helper: margins -> glm formula with hierarchy (using * notation) ---
+
+.margins_to_glm_formula <- function(margins, response) {
+  rhs_terms <- vapply(margins, function(m) {
+    if (length(m) == 1) m else paste(m, collapse = " * ")
+  }, character(1))
+  stats::as.formula(paste(response, "~", paste(rhs_terms, collapse = " + ")))
+}
+
+
+# --- Helper: partial chi-squared weights (deviance change per edge) ---
+
+.compute_chisq_weights <- function(g, model) {
+  el <- igraph::as_edgelist(g)
+  if (nrow(el) == 0) return(numeric(0))
+
+  dev_full <- stats::deviance(model)
+  resp <- as.character(stats::formula(model)[[2]])
+
+  weights <- numeric(nrow(el))
+  for (i in seq_len(nrow(el))) {
+    # Remove this edge
+    eid <- igraph::get.edge.ids(g, c(el[i, 1], el[i, 2]))
+    g_reduced <- igraph::delete_edges(g, eid)
+
+    # Convert reduced graph -> margins -> glm formula
+    margins_reduced <- .graph_to_margins(g_reduced)
+    new_formula <- .margins_to_glm_formula(margins_reduced, resp)
+
+    # Refit and compute deviance change
+    mod_reduced <- stats::update(model, formula. = new_formula)
+    weights[i] <- stats::deviance(mod_reduced) - dev_full
+  }
+  weights
+}
+
+
+# --- Helper: Cramer's V weights (marginal 2-way tables) ---
+
+.compute_cramer_weights <- function(g, model) {
+  el <- igraph::as_edgelist(g)
+  if (nrow(el) == 0) return(numeric(0))
+
+  mf <- stats::model.frame(model)
+  freq <- stats::model.response(mf)
+
+  weights <- numeric(nrow(el))
+  for (i in seq_len(nrow(el))) {
+    tab <- stats::xtabs(freq ~ mf[[el[i, 1]]] + mf[[el[i, 2]]])
+    weights[i] <- vcd::assocstats(tab)$cramer
+  }
+  weights
+}
+
+
 # --- Helper: extract generating class (margins) from a glm formula ---
 
 .glm_to_margins <- function(object) {
@@ -225,12 +357,22 @@ print.assoc_graph <- function(x, ...) {
 
   if (ne > 0) {
     el <- igraph::as_edgelist(x)
-    edge_strings <- paste0(el[, 1], " -- ", el[, 2])
+    w <- igraph::E(x)$weight
+    if (!is.null(w)) {
+      edge_strings <- paste0(el[, 1], " -- ", el[, 2], " (", round(w, 2), ")")
+    } else {
+      edge_strings <- paste0(el[, 1], " -- ", el[, 2])
+    }
     cat("Edges:", paste(edge_strings, collapse = ", "), "\n")
   } else {
     cat("Edges: (none -- mutual independence)\n")
   }
 
+  # Show measure if present
+  if (!is.null(x$measure) && x$measure != "none") {
+    cat("Measure:", x$measure, "\n")
+  }
+
   # Show bracket notation
   margins <- .graph_to_margins(x)
   cat("Model:", loglin2string(margins), "\n")

dev/LRstats-test.R

@@ -3,54 +3,6 @@
 
 library(vcdExtra)
 
-# --- Modified methods with label argument ---
-
-LRstats.glmlist <- function(object, ..., saturated = NULL, sortby = NULL,
-                            label = c("name", "formula"), abbrev = FALSE) {
-  label <- match.arg(label)
-  ns <- sapply(object, function(x) length(x$residuals))
-  if (any(ns != ns[1L]))
-    stop("models were not all fitted to the same size of dataset")
-  nmodels <- length(object)
-  if (nmodels == 1)
-    return(LRstats.default(object[[1L]], saturated = saturated))
-
-  rval <- lapply(object, LRstats.default, saturated = saturated)
-  rval <- do.call(rbind, rval)
-
-  if (label == "formula") {
-    rownames(rval) <- get_models(object, abbrev = abbrev)
-  }
-
-  if (!is.null(sortby)) {
-    rval <- rval[order(rval[, sortby], decreasing = TRUE), ]
-  }
-  rval
-}
-
-LRstats.loglmlist <- function(object, ..., saturated = NULL, sortby = NULL,
-                              label = c("name", "formula"), abbrev = FALSE) {
-  label <- match.arg(label)
-  ns <- sapply(object, function(x) length(x$residuals))
-  if (any(ns != ns[1L]))
-    stop("models were not all fitted to the same size of dataset")
-  nmodels <- length(object)
-  if (nmodels == 1)
-    return(LRstats.default(object[[1L]], saturated = saturated))
-
-  rval <- lapply(object, LRstats.default, saturated = saturated)
-  rval <- do.call(rbind, rval)
-
-  if (label == "formula") {
-    rownames(rval) <- get_models(object, abbrev = abbrev)
-  }
-
-  if (!is.null(sortby)) {
-    rval <- rval[order(rval[, sortby], decreasing = TRUE), ]
-  }
-  rval
-}
-
 # --- Examples ---
 
 cat("\n=== loglmlist: Sequential joint independence (Titanic) ===\n")