Merge pull request #1517 from exploratory-io/fix/issue-30488

feat(#30488): support auto/polychoric/mixed correlation in do_cor (16.0.14)

Author: hidekoji

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: exploratory
 Type: Package
 Title: R package for Exploratory
-Version: 16.0.13
+Version: 16.0.14
 Date: 2026-07-24
 Authors@R: c(person("Hideaki", "Hayashi", email = "hideaki@exploratory.io", role = c("aut", "cre")), person("Hide", "Kojima", email = "hide@exploratory.io", role = c("aut")), person("Kan", "Nishida", email = "kan@exploratory.io", role = c("aut")), person("Kei", "Saito", email = "kei@exploratory.io", role = c("aut")), person("Yosuke", "Yasuda", email = "double.y.919.quick@gmail.com", role = c("aut")))
 URL: https://github.com/exploratory-io/exploratory_func
@@ -45,6 +45,7 @@ Imports:
     pwr,
     readr
 Suggests:
+    polycor,
     testthat,
     twitteR,
     sentimentr,

R/stats_wrapper.R

@@ -45,7 +45,7 @@ do_cor <- function(df, ..., skv = NULL, fun.aggregate=mean, fill=0){
 #' @param dimension A column you want to use as a dimension to calculate the correlations.
 #' @param value A column for the values you want to use to calculate the correlations.
 #' @param use Operation type for dealing with missing values. This can be one of "everything", "all.obs", "complete.obs", "na.or.complete", or "pairwise.complete.obs"
-#' @param method Method of calculation. This can be one of "pearson", "kendall", or "spearman".
+#' @param method Method of calculation. This can be one of "auto", "pearson", "kendall", "spearman", "polychoric", or "mixed".
 #' @param fun.aggregate  Set an aggregate function when there are multiple entries for the key column per each category.
 #' @param time_unit Unit of time to aggregate key_col if key_col is Date or POSIXct. NULL doesn't aggregate.
 #' @return correlations between pairs of groups
@@ -175,7 +175,7 @@ do_cor.kv_ <- function(df,
 #' @param df data frame in tidy format
 #' @param ... Arguments to select columns to calculate correlation.
 #' @param use Operation type for dealing with missing values. This can be one of "everything", "all.obs", "complete.obs", "na.or.complete", or "pairwise.complete.obs"
-#' @param method Method of calculation. This can be one of "pearson", "kendall", or "spearman".
+#' @param method Method of calculation. This can be one of "auto", "pearson", "kendall", "spearman", "polychoric", or "mixed".
 #' @return correlations between pairs of columns
 #' @export
 do_cor.cols <- function(df, ..., use = "pairwise.complete.obs", method = "pearson",
@@ -205,10 +205,9 @@ do_cor.cols <- function(df, ..., use = "pairwise.complete.obs", method = "pearso
         stop("More than 1 row are required to calculate correlations.")
       }
     }
-    mat <- dplyr::select(df, !!!select_dots) %>%
-      # Convert logical to numeric explicitly, since implicit conversion by as.matrix does not happen if all the columns are logical.
-      dplyr::mutate(across(where(is.logical), as.numeric)) %>%
-      as.matrix()
+    # Keep the source column types here. Auto and mixed correlation need to distinguish
+    # continuous numeric variables from ordinal (factor/logical) variables.
+    mat <- dplyr::select(df, !!!select_dots) %>% as.data.frame()
 
     ret <- do_cor_internal(mat, use, method, diag, output_cols, na.rm=TRUE)
 
@@ -253,50 +252,130 @@ do_cor.cols <- function(df, ..., use = "pairwise.complete.obs", method = "pearso
 }
 
 
+resolve_correlation_method <- function(mat, method) {
+  if (!identical(method, "auto")) {
+    return(method)
+  }
+
+  is_ordinal <- vapply(mat, function(x) is.factor(x) || is.logical(x), logical(1))
+  is_continuous <- vapply(mat, is.numeric, logical(1))
+
+  if (all(is_continuous)) {
+    return("pearson")
+  }
+  if (all(is_ordinal)) {
+    return("polychoric")
+  }
+  if (all(is_continuous | is_ordinal)) {
+    return("mixed")
+  }
+
+  warning("Could not determine all variable types for automatic correlation; using Pearson.")
+  "pearson"
+}
+
+as_numeric_correlation_matrix <- function(df) {
+  as.matrix(as.data.frame(lapply(df, function(x) {
+    if (is.logical(x)) {
+      as.numeric(x)
+    } else if (is.factor(x)) {
+      as.numeric(x)
+    } else {
+      x
+    }
+  }), check.names = FALSE))
+}
+
+as_hetcor_data_frame <- function(df, force_ordinal = FALSE) {
+  ret <- as.data.frame(lapply(df, function(x) {
+    if (force_ordinal || is.factor(x) || is.logical(x)) {
+      ordered(x)
+    } else {
+      x
+    }
+  }))
+  colnames(ret) <- colnames(df)
+  ret
+}
+
 do_cor_internal <- function(mat, use, method, diag, output_cols, na.rm) {
   # Sort the column name.
   # Now that we sort the variables based on correlation result or input order, this might not be as meaningful,
   # but I'm hoping this might still help align the result between Mac and Windows if there are ties in the mean correlations.
   # We use stringr::str_sort() as opposed to base sort() so that the result is consistent on Windows too.
   sorted_colnames <- stringr::str_sort(colnames(mat))
-  mat <- mat[,sorted_colnames]
+  # `as.data.frame.matrix()` otherwise changes numeric matrix column names such
+  # as "1" to "X1", breaking the key-value correlation output contract.
+  mat <- as.data.frame(mat[, sorted_colnames, drop = FALSE], check.names = FALSE)
+  method <- resolve_correlation_method(mat, method)
+
+  # Create a matrix of P-values for Analytics View case.
+  dim <- length(sorted_colnames)
+
+  if (method %in% c("polychoric", "mixed")) {
+    # Polychoric correlation is for ordinal variables (e.g. survey scales 1-5).
+    # stats::cor()/cor.test() do not support it, so we use polycor::hetcor(), which
+    # returns the full correlation matrix and the matrix of standard errors in one call.
+    # Polychoric coerces every column to ordered so every pair uses the ordinal model.
+    # Mixed preserves continuous numerics and converts factor/logical columns to ordered,
+    # allowing hetcor() to select Pearson, polyserial, or polychoric per pair.
+    if (!requireNamespace("polycor", quietly = TRUE)) {
+      stop(
+        "The 'polycor' package is required for polychoric and mixed correlations. ",
+        "Install it before using this method.",
+        call. = FALSE
+      )
+    }
+    ordered_df <- as_hetcor_data_frame(mat, force_ordinal = identical(method, "polychoric"))
+    # hetcor only supports "complete.obs" and "pairwise.complete.obs", while the other
+    # methods (via cor()/cor.test()) also accept "everything", "all.obs", and
+    # "na.or.complete". Map those to pairwise (the do_cor default) so polychoric does not
+    # error where the other methods would succeed. (hetcor's own default is listwise
+    # "complete.obs", which we do not want.)
+    het_use <- if (identical(use, "complete.obs")) "complete.obs" else "pairwise.complete.obs"
+    het <- polycor::hetcor(ordered_df, ML = TRUE, std.err = TRUE, use = het_use)
+    cor_mat <- het$correlations
+    tvalue_mat <- het$correlations / het$std.errors # z value; off-diagonal NA SE stays NA.
+    pvalue_mat <- 2 * stats::pnorm(-abs(tvalue_mat))
+    diag(tvalue_mat) <- Inf # For i=j case, statistic should be Inf and P value 0.
+    diag(pvalue_mat) <- 0
+  } else {
+    numeric_mat <- as_numeric_correlation_matrix(mat)
+    cor_mat <- cor(numeric_mat, use = use, method = method)
+    pvalue_mat <- matrix(NA, dim, dim)
+    tvalue_mat <- matrix(NA, dim, dim)
+    for (i in 2:dim) {
+      for (j in 1:(i-1)) {
+        tryCatch({
+          cor_test_res <- cor.test(numeric_mat[,i], numeric_mat[,j], method = method)
+          pvalue_mat[i, j] <- cor_test_res$p.value
+          tvalue_mat[i, j] <- cor_test_res$statistic
+        }, error = function(e) {
+          if (e$message == "not enough finite observations") {
+            # This is the error cor.test returns when there is not enough non-NA data.
+            # Rather than stopping, set NA as the result, and we will handle it as a not-significant case on the UI.
+            pvalue_mat[i, j] <- NA
+            tvalue_mat[i, j] <- NA
+          }
+          else {
+            stop(e)
+          }
+        })
+        pvalue_mat[j, i] <- pvalue_mat[i, j]
+        tvalue_mat[j, i] <- tvalue_mat[i, j]
+      }
+    }
+    for (i in 1:dim) { # For i=j case, P value should be always 0 and t statistic should be Inf.
+      pvalue_mat[i, i] <- 0
+      tvalue_mat[i, i] <- Inf
+    }
+  }
 
-  cor_mat <- cor(mat, use = use, method = method)
   ret <- mat_to_df(cor_mat,
                    cnames=output_cols[1:3],
                    diag=diag,
                    na.rm=na.rm,
                    zero.rm=FALSE)
-
-  # Create a matrix of P-values for Analytics View case.
-  dim <- length(sorted_colnames)
-  pvalue_mat <- matrix(NA, dim, dim)
-  tvalue_mat <- matrix(NA, dim, dim)
-  for (i in 2:dim) {
-    for (j in 1:(i-1)) {
-      tryCatch({
-        cor_test_res <- cor.test(mat[,i], mat[,j], method = method)
-        pvalue_mat[i, j] <- cor_test_res$p.value
-        tvalue_mat[i, j] <- cor_test_res$statistic
-      }, error = function(e) {
-        if (e$message == "not enough finite observations") {
-          # This is the error cor.test returns when there is not enough non-NA data.
-          # Rather than stopping, set NA as the result, and we will handle it as a not-significant case on the UI.
-          pvalue_mat[i, j] <- NA
-          tvalue_mat[i, j] <- NA
-        }
-        else {
-          stop(e)
-        }
-      })
-      pvalue_mat[j, i] <- pvalue_mat[i, j]
-      tvalue_mat[j, i] <- tvalue_mat[i, j]
-    }
-  }
-  for (i in 1:dim) { # For i=j case, P value should be always 0 and t statistic should be Inf.
-    pvalue_mat[i, i] <- 0
-    tvalue_mat[i, i] <- Inf
-  }
   colnames(pvalue_mat) <- sorted_colnames
   rownames(pvalue_mat) <- sorted_colnames
   colnames(tvalue_mat) <- sorted_colnames

tests/testthat/test_stats_wrapper.R

@@ -120,6 +120,153 @@ test_that("do_cor should skip group with only one row.", {
   expect_equal(nrow(res %>% filter(z==F)), 0)
 })
 
+test_that("do_cor with polychoric method", {
+  skip_if_not_installed("polycor")
+
+  # Polychoric correlation is for ordinal variables. The latent variables behind x and y
+  # are correlated at 0.7, while z is independent of them.
+  set.seed(123)
+  n <- 200
+  cut5 <- function(z) as.integer(cut(z, breaks = c(-Inf, -0.84, -0.25, 0.25, 0.84, Inf)))
+  z1 <- rnorm(n); z2 <- 0.7 * z1 + sqrt(1 - 0.49) * rnorm(n); z3 <- rnorm(n)
+  df <- data.frame(x = cut5(z1), y = cut5(z2), z = cut5(z3))
+  model_df <- df %>% do_cor(`x`, `y`, `z`, method = "polychoric", distinct = FALSE, diag = TRUE, return_type = "model")
+  res <- model_df %>% tidy_rowwise(model, type = 'cor')
+  expect_equal(nrow(res), 9) # All 9 combinations.
+  xy <- res %>% filter(pair.name.x == "x", pair.name.y == "y")
+  expect_true(xy$correlation > 0.5 && xy$correlation < 0.9) # Recovers the true latent 0.7, not inflated.
+  expect_true(xy$p_value < 0.05) # Statistically significant.
+  expect_true(is.finite(xy$statistic)) # z value is populated.
+  xz <- res %>% filter(pair.name.x == "x", pair.name.y == "z")
+  expect_true(abs(xz$correlation) < 0.3) # z is independent of x: near-zero polychoric correlation.
+  expect_true(xz$p_value > 0.05) # The independent pair is not statistically significant.
+  diag_row <- res %>% filter(pair.name.x == "x", pair.name.y == "x")
+  expect_equal(diag_row$correlation, 1) # Diagonal correlation is 1.
+  expect_equal(diag_row$p_value, 0) # Diagonal P value is 0.
+})
+
+test_that("do_cor automatic method selects the specified correlation family", {
+  ordinal <- ordered(rep(1:5, each = 20))
+  expect_identical(resolve_correlation_method(data.frame(x = 1:100, y = 101:200), "auto"), "pearson")
+  expect_identical(resolve_correlation_method(data.frame(x = ordinal, y = rev(ordinal)), "auto"), "polychoric")
+  expect_identical(resolve_correlation_method(data.frame(x = 1:100, y = ordinal), "auto"), "mixed")
+})
+
+test_that("do_cor supports automatic and mixed correlations with factor inputs", {
+  skip_if_not_installed("polycor")
+
+  set.seed(456)
+  n <- 100
+  ordinal <- ordered(cut(rnorm(n), breaks = c(-Inf, -0.84, -0.25, 0.25, 0.84, Inf)))
+  df <- data.frame(continuous = rnorm(n), ordinal = ordinal)
+
+  auto <- df %>% do_cor(continuous, ordinal, method = "auto", distinct = FALSE, diag = TRUE, return_type = "model") %>% tidy_rowwise(model, type = "cor")
+  mixed <- df %>% do_cor(continuous, ordinal, method = "mixed", distinct = FALSE, diag = TRUE, return_type = "model") %>% tidy_rowwise(model, type = "cor")
+
+  expect_equal(nrow(auto), 4)
+  expect_equal(nrow(mixed), 4)
+  expect_true(is.finite(auto$correlation[auto$pair.name.x != auto$pair.name.y][1]))
+  expect_true(is.finite(mixed$correlation[mixed$pair.name.x != mixed$pair.name.y][1]))
+})
+
+test_that("do_cor with polychoric method handles a constant column without error", {
+  skip_if_not_installed("polycor")
+
+  set.seed(123)
+  n <- 100
+  cut5 <- function(z) as.integer(cut(z, breaks = c(-Inf, -0.84, -0.25, 0.25, 0.84, Inf)))
+  z1 <- rnorm(n); z2 <- 0.7 * z1 + sqrt(1 - 0.49) * rnorm(n)
+  df <- data.frame(x = cut5(z1), y = cut5(z2), w = rep(3L, n)) # w is constant.
+  # hetcor warns for the non-estimable pairs involving the constant column; that is expected.
+  model_df <- suppressWarnings(df %>% do_cor(`x`, `y`, `w`, method = "polychoric", distinct = FALSE, diag = TRUE, return_type = "model"))
+  res <- model_df %>% tidy_rowwise(model, type = 'cor')
+  xy <- res %>% filter(pair.name.x == "x", pair.name.y == "y")
+  expect_equal(nrow(xy), 1) # The estimable pair is still computed.
+  expect_true(xy$correlation > 0.5)
+  # The constant-column pair has NA correlation and is dropped by na.rm in mat_to_df.
+  expect_equal(nrow(res %>% filter(pair.name.x == "x", pair.name.y == "w")), 0)
+})
+
+test_that("do_cor with polychoric method handles complex column names", {
+  skip_if_not_installed("polycor")
+
+  set.seed(123)
+  n <- 100
+  cut5 <- function(z) as.integer(cut(z, breaks = c(-Inf, -0.84, -0.25, 0.25, 0.84, Inf)))
+  z1 <- rnorm(n); z2 <- 0.7 * z1 + sqrt(1 - 0.49) * rnorm(n)
+  sname <- "航空 会社 !\"#$%&'()*+, -./:;<=>?@[]^_'{|}~ 表"
+  df <- data.frame(a = cut5(z1), b = cut5(z2), check.names = FALSE)
+  names(df) <- c(sname, "plain")
+  model_df <- df %>% do_cor(tidyselect::everything(), method = "polychoric", distinct = FALSE, diag = TRUE, return_type = "model")
+  res <- model_df %>% tidy_rowwise(model, type = 'cor')
+  expect_true(sname %in% as.character(res$pair.name.x)) # Complex name survives the round trip.
+  pair <- res %>% filter(as.character(pair.name.x) == sname, pair.name.y == "plain")
+  expect_equal(nrow(pair), 1)
+  expect_true(is.finite(pair$correlation))
+})
+
+test_that("do_cor with polychoric method accepts use values hetcor does not support", {
+  skip_if_not_installed("polycor")
+
+  # hetcor only accepts "complete.obs" and "pairwise.complete.obs", but the public
+  # do_cor API (and the other methods via cor()/cor.test()) also accept "everything",
+  # "all.obs", and "na.or.complete". Those must be mapped, not passed through as an error.
+  set.seed(123)
+  n <- 100
+  cut5 <- function(z) as.integer(cut(z, breaks = c(-Inf, -0.84, -0.25, 0.25, 0.84, Inf)))
+  z1 <- rnorm(n); z2 <- 0.7 * z1 + sqrt(1 - 0.49) * rnorm(n)
+  df <- data.frame(x = cut5(z1), y = cut5(z2))
+  for (u in c("everything", "all.obs", "na.or.complete", "complete.obs")) {
+    model_df <- df %>% do_cor(`x`, `y`, method = "polychoric", use = u, distinct = FALSE, diag = TRUE, return_type = "model")
+    res <- model_df %>% tidy_rowwise(model, type = 'cor')
+    xy <- res %>% filter(pair.name.x == "x", pair.name.y == "y")
+    expect_equal(nrow(xy), 1, info = u) # The pair is computed regardless of the use value.
+    expect_true(xy$correlation > 0.5, info = u)
+  }
+})
+
+test_that("do_cor with polychoric method for grouped (repeat-by) data", {
+  skip_if_not_installed("polycor")
+
+  # Repeat By on Analytics View maps to group_by(). Each group must get its own
+  # polychoric correlation. Group A has a positive relationship, group B a negative one.
+  set.seed(123)
+  n <- 100
+  cut5 <- function(z) as.integer(cut(z, breaks = c(-Inf, -0.84, -0.25, 0.25, 0.84, Inf)))
+  mk <- function(rho) {
+    z1 <- rnorm(n); z2 <- rho * z1 + sqrt(1 - rho^2) * rnorm(n)
+    data.frame(x = cut5(z1), y = cut5(z2))
+  }
+  df <- dplyr::bind_rows(cbind(mk(0.7), grp = "A"), cbind(mk(-0.7), grp = "B"))
+  model_df <- df %>% group_by(grp) %>% do_cor(`x`, `y`, method = "polychoric", distinct = FALSE, diag = TRUE, return_type = "model")
+  res <- model_df %>% tidy_rowwise(model, type = 'cor')
+  expect_setequal(unique(as.character(res$grp)), c("A", "B")) # Both groups produced results.
+  a_xy <- res %>% filter(grp == "A", pair.name.x == "x", pair.name.y == "y")
+  b_xy <- res %>% filter(grp == "B", pair.name.x == "x", pair.name.y == "y")
+  expect_true(a_xy$correlation > 0.4) # Positive correlation in group A.
+  expect_true(b_xy$correlation < -0.4) # Negative correlation in group B, computed independently.
+})
+
+test_that("do_cor with polychoric method handles NA values via pairwise complete obs", {
+  skip_if_not_installed("polycor")
+
+  # Survey data routinely has missing responses (NA). The default use="pairwise.complete.obs"
+  # must drop NAs pairwise rather than error, so every pair is still estimated.
+  set.seed(123)
+  n <- 200
+  cut5 <- function(z) as.integer(cut(z, breaks = c(-Inf, -0.84, -0.25, 0.25, 0.84, Inf)))
+  z1 <- rnorm(n); z2 <- 0.7 * z1 + sqrt(1 - 0.49) * rnorm(n); z3 <- rnorm(n)
+  x <- cut5(z1); y <- cut5(z2); z <- cut5(z3)
+  x[1:20] <- NA; y[15:30] <- NA; z[40:60] <- NA # Missing responses at different rows per column.
+  df <- data.frame(x = x, y = y, z = z)
+  model_df <- df %>% do_cor(`x`, `y`, `z`, method = "polychoric", distinct = FALSE, diag = TRUE, return_type = "model")
+  res <- model_df %>% tidy_rowwise(model, type = 'cor')
+  expect_equal(nrow(res), 9) # Every pair is still estimated from the pairwise-complete rows.
+  xy <- res %>% filter(pair.name.x == "x", pair.name.y == "y")
+  expect_true(is.finite(xy$correlation)) # NAs did not break the estimate.
+  expect_true(xy$correlation > 0.4) # The x-y relationship is still recovered despite the NAs.
+})
+
 test_that("test do_svd.kv with fill", {
   test_df <- data.frame(
     rand=runif(20, min = 0, max=10),