Stabilize multi-way cluster SE diagnostics

Co-authored-by: kaiemjoy <16113030+kaiemjoy@users.noreply.github.com>

Author: Copilot

NEWS.md

@@ -29,9 +29,10 @@
   their snapshots, differ between macOS, Windows, and Linux). Simulated
   values change slightly as a result of this fix. (#447)
 * Cluster-robust standard errors now treat multiple `cluster_var` columns as
-  multi-way clustering instead of collapsing them to a single interaction, and
-  they are no longer allowed to be smaller than the corresponding model-based
-  standard errors. (#543)
+  multi-way clustering instead of collapsing them to a single interaction.
+  One-way subset terms now use the CR1 small-sample correction by default,
+  the floor to model-based variance is optional, and debug output now exposes
+  the multi-way variance decomposition used in the final standard error. (#543)
 * Corrected default axis labels in `strat_ests_barplot()` (`xlab`) and
   `strat_ests_scatterplot()` (`ylab`) to say "seroincidence" rather than
   "seroconversion"/"incidence".

R/compute_cluster_robust_var.R

@@ -16,10 +16,14 @@
 .compute_cluster_robust_var <- function(
     fit,
     cluster_var,
-    stratum_var = NULL) {
+    stratum_var = NULL,
+    small_sample = c("none", "CR1"),
+    floor_to_standard = FALSE,
+    debug_cluster = FALSE) {
+  small_sample <- match.arg(small_sample)
   pop_data_list <- attr(fit, "pop_data")
   pop_data_combined <- do.call(rbind, pop_data_list)
-  standard_var_log_lambda <- 1 / fit$hessian |> as.numeric()
+  standard_var_log_lambda <- as.numeric(1 / fit$hessian)[1]
 
   cluster_var_combinations <- unlist(
     lapply(seq_along(cluster_var), function(n_vars) {
@@ -29,7 +33,7 @@
   )
 
   n_vars_per_subset <- vapply(cluster_var_combinations, length, integer(1))
-  robust_var_log_lambda <- 0
+  decomp_rows <- vector("list", length(cluster_var_combinations))
 
   for (i in seq_along(cluster_var_combinations)) {
     cluster_vars_subset <- cluster_var_combinations[[i]]
@@ -46,20 +50,73 @@
     subset_var_log_lambda <- .compute_cluster_var_oneway(
       fit = fit,
       cluster_ids = cluster_ids,
-      pop_data_combined = pop_data_combined
+      pop_data_combined = pop_data_combined,
+      small_sample = small_sample
+    )
+    subset_sign <- (-1)^(n_vars_per_subset[i] + 1)
+    decomp_rows[[i]] <- tibble::tibble(
+      subset = paste(cluster_vars_subset, collapse = " + "),
+      order = n_vars_per_subset[i],
+      sign = subset_sign,
+      subset_variance = subset_var_log_lambda,
+      signed_term = subset_sign * subset_var_log_lambda
     )
-    robust_var_log_lambda <- robust_var_log_lambda +
-      (-1)^(n_vars_per_subset[i] + 1) * subset_var_log_lambda
   }
 
-  # Use a conservative floor so cluster-robust variance does not fall below the
-  # model-based variance when clustering adds no effective dependence (for
-  # example, singleton clusters) or finite-sample noise would otherwise shrink
-  # the standard error.
-  robust_var_log_lambda <- max(
-    standard_var_log_lambda,
-    robust_var_log_lambda
-  )
+  decomp_terms <- dplyr::bind_rows(decomp_rows)
+  robust_raw <- sum(decomp_terms$signed_term)
+  floor_applied <- isTRUE(floor_to_standard) &&
+    robust_raw < standard_var_log_lambda
+  robust_final <- if (isTRUE(floor_to_standard)) {
+    max(standard_var_log_lambda, robust_raw)
+  } else {
+    robust_raw
+  }
+
+  if (debug_cluster) {
+    cli::cli_inform("Cluster-robust variance decomposition:")
+
+    for (i in seq_len(nrow(decomp_terms))) {
+      term_label <- if (
+        length(cluster_var) == 2 &&
+          decomp_terms$order[i] == 2
+      ) {
+        "V_intersection"
+      } else {
+        paste0(
+          "V_",
+          gsub(" \\+ ", "_", decomp_terms$subset[i])
+        )
+      }
+      term_message <- glue::glue(
+        "{term_label} ({decomp_terms$subset[i]}) = ",
+        "{signif(decomp_terms$subset_variance[i], 6)}"
+      )
+
+      cli::cli_inform(term_message)
+    }
 
-  return(robust_var_log_lambda)
+    cli::cli_inform("V_raw = {signif(robust_raw, 6)}")
+    cli::cli_inform("V_final = {signif(robust_final, 6)}")
+
+    if (isTRUE(floor_to_standard)) {
+      floor_message <- glue::glue(
+        "Variance floor relative to standard variance ",
+        "{signif(standard_var_log_lambda, 6)}: ",
+        "{if (floor_applied) 'applied' else 'not applied'}"
+      )
+      cli::cli_inform(floor_message)
+    }
+  }
+
+  structure(
+    robust_final,
+    cluster_decomp = list(
+      standard_var = standard_var_log_lambda,
+      robust_raw = robust_raw,
+      robust_final = robust_final,
+      terms = decomp_terms,
+      floor_applied = floor_applied
+    )
+  )
 }

R/compute_cluster_var_oneway.R

@@ -10,7 +10,9 @@
 .compute_cluster_var_oneway <- function(
     fit,
     cluster_ids,
-    pop_data_combined) {
+    pop_data_combined,
+    small_sample = c("none", "CR1")) {
+  small_sample <- match.arg(small_sample)
   pop_data_list <- attr(fit, "pop_data")
   sr_params_list <- attr(fit, "sr_params")
   noise_params_list <- attr(fit, "noise_params")
@@ -19,7 +21,8 @@
   epsilon <- 1e-6
 
   unique_clusters <- unique(cluster_ids)
-  cluster_scores <- numeric(length(unique_clusters))
+  n_clusters <- length(unique_clusters)
+  cluster_scores <- numeric(n_clusters)
 
   for (i in seq_along(unique_clusters)) {
     cluster_id <- unique_clusters[i]
@@ -57,7 +60,17 @@
   }
 
   score_variance <- sum(cluster_scores^2)
-  hessian <- fit$hessian
+  hessian <- as.numeric(fit$hessian)[1]
 
-  score_variance / (hessian^2)
+  if (!is.finite(hessian) || hessian <= 0) {
+    return(NA_real_)
+  }
+
+  var_log_lambda <- score_variance / (hessian^2)
+
+  if (small_sample == "CR1" && n_clusters > 1) {
+    var_log_lambda <- var_log_lambda * (n_clusters / (n_clusters - 1))
+  }
+
+  var_log_lambda
 }

R/summary.seroincidence.R

@@ -8,6 +8,13 @@
 #' @param object a [list()] outputted by [stats::nlm()] or [est_seroincidence()]
 #' @param coverage desired confidence interval coverage probability
 #' @param verbose whether to produce verbose messaging
+#' @param small_sample small-sample correction for cluster-robust variance.
+#'   Use `"CR1"` to apply the one-way CR1 factor to each subset term, or
+#'   `"none"` to leave the one-way terms unadjusted.
+#' @param floor_to_standard whether to floor cluster-robust variance at the
+#'   model-based variance.
+#' @param debug_cluster whether to print the cluster-robust variance
+#'   decomposition when clustering is used.
 #' @param ... unused
 #'
 #' @return a [tibble::tibble()] containing the following:
@@ -71,7 +78,11 @@ summary.seroincidence <- function(
     object,
     coverage = .95,
     verbose = TRUE,
+    small_sample = c("none", "CR1"),
+    floor_to_standard = FALSE,
+    debug_cluster = FALSE,
     ...) {
+  small_sample <- match.arg(small_sample)
   start <- object |> attr("lambda_start")
   antigen_isos <- object |> attr("antigen_isos")
   cluster_var <- object |> attr("cluster_var")
@@ -106,13 +117,17 @@ summary.seroincidence <- function(
     var_log_lambda <- .compute_cluster_robust_var(
       fit = object,
       cluster_var = cluster_var,
-      stratum_var = stratum_var
+      stratum_var = stratum_var,
+      small_sample = small_sample,
+      floor_to_standard = floor_to_standard,
+      debug_cluster = debug_cluster
     )
   } else {
     # Standard variance from Hessian
     var_log_lambda <- 1 / object$hessian |> as.vector()
   }
 
+  cluster_decomp <- attr(var_log_lambda, "cluster_decomp")
   # Ensure var_log_lambda is a scalar
   var_log_lambda <- as.numeric(var_log_lambda)[1]
   se_log_lambda <- sqrt(var_log_lambda)
@@ -140,5 +155,9 @@ summary.seroincidence <- function(
     "summary.seroincidence" |>
     union(class(to_return))
 
+  if (!is.null(cluster_decomp)) {
+    attr(to_return, "cluster_decomp") <- cluster_decomp
+  }
+
   return(to_return)
 }

R/summary.seroincidence.by.R

@@ -17,6 +17,7 @@
 #' (see help for [optim()] for details).
 #' Default = `FALSE`.
 #' @param confidence_level desired confidence interval coverage probability
+#' @inheritParams summary.seroincidence
 #' @return
 #' A `summary.seroincidence.by` object, which is a [tibble::tibble],
 #' with the following columns:
@@ -69,7 +70,11 @@ summary.seroincidence.by <- function(
     show_deviance = TRUE,
     show_convergence = TRUE,
     verbose = FALSE,
+    small_sample = c("none", "CR1"),
+    floor_to_standard = FALSE,
+    debug_cluster = FALSE,
     ...) {
+  small_sample <- match.arg(small_sample)
   alpha <- 1 - confidence_level
   quantiles <- c(alpha / 2, 1 - alpha / 2)
 
@@ -87,7 +92,10 @@ summary.seroincidence.by <- function(
     lapply(
       FUN = summary.seroincidence,
       coverage = confidence_level,
-      verbose = verbose
+      verbose = verbose,
+      small_sample = small_sample,
+      floor_to_standard = floor_to_standard,
+      debug_cluster = debug_cluster
     ) |>
     bind_rows(.id = "Stratum")