streamline kendall cor functions

R/kendall_correlation.R

@@ -31,14 +31,14 @@
 #' @export
 kendall_cor <- function(x, y) {
   arr <- cbind(x, y)
-  storage.mode(arr) <- "double"
+  if (storage.mode(arr) != "double") { storage.mode(arr) <- "double" }
   arr <- arr[complete.cases(arr), ]
 
-  kw <- kendall_warnings(arr)
+  n <- ncol(arr)
 
-  if (isFALSE(kw)) {
-    return(NA_real_)
-  }
+  kw <- kendall_warnings(arr, n)
+
+  if (isFALSE(kw)) { return(NA) }
 
   kendall_cor_(arr)
 }
@@ -70,29 +70,22 @@ kendall_cor <- function(x, y) {
 #' kendall_cor_test(x, y)
 #' 
 #' @export
-kendall_cor_test <- function(x,y,
+kendall_cor_test <- function(x, y,
   alternative = c("two.sided", "greater", "less")) {
   alternative <- match.arg(alternative)
 
   arr <- cbind(x, y)
-  storage.mode(arr) <- "double"
+  if (storage.mode(arr) != "double") { storage.mode(arr) <- "double" }
   arr <- arr[complete.cases(arr), ]
 
-  kw <- kendall_warnings(arr)
+  n <- nrow(arr)
+  m <- ncol(arr)
 
-  if (isFALSE(kw)) {
-    return(NA)
-  }
+  kw <- kendall_warnings(arr, m)
 
-  r <- kendall_cor_(arr)
-  n <- nrow(arr)
+  if (isFALSE(kw)) { return(NA) }
 
-  if (n < 2) {
-    stop("not enough finite observations")
-  }
-
-  # r = correlation
-  # n = number of observations
+  r <- kendall_cor_(arr)
 
   if (n < 50) {
     q <- round((r + 1) * n * (n - 1) / 4)
@@ -145,21 +138,21 @@ kendall_cor_test <- function(x,y,
   )
 }
 
-kendall_warnings <- function(arr) {
-  if (ncol(arr) != 2) {
+kendall_warnings <- function(arr, n) {
+  if (n != 2) {
     stop("x and y must be uni-dimensional vectors")
   }
 
-  if (nrow(arr) < 2) {
-    stop("x and y must have at least 2 observations")
+  if (n < 2) {
+    stop("x and y must have at least 2 non-null observations")
   }
 
-  if (sd(arr[, 1]) == 0) {
+  if (var(arr[, 1]) == 0) {
     warning("x has zero variance")
     return(FALSE)
   }
 
-  if (sd(arr[, 2]) == 0) {
+  if (var(arr[, 2]) == 0) {
     warning("y has zero variance")
     return(FALSE)
   }

src/05_kendall_correlation.cpp

@@ -3,6 +3,9 @@
 // Data)
 // Filzmoser, Fritz, and Kalcher 2023 (pcaPP package)
 
+// note: the len < 2 conditions are commented out because the R function checks
+// for this condition before calling the C++ functions
+
 #include <algorithm>
 #include <cmath>
 #include <cpp11.hpp>
@@ -13,9 +16,9 @@
 using namespace cpp11;
 
 uint64_t insertion_sort_(double *arr, size_t len) {
-  if (len < 2) {
-    return 0;
-  }
+  // if (len < 2) {
+  //   return 0;
+  // }
 
   size_t maxJ = len - 1, i;
   uint64_t swapCount = 0;
@@ -70,9 +73,9 @@ static uint64_t merge_(double *from, double *to, size_t middle, size_t len) {
 }
 
 uint64_t merge_sort_(double *x, double *buf, size_t len) {
-  if (len < 2) {
-    return 0;
-  }
+  // if (len < 2) {
+  //   return 0;
+  // }
 
   if (len < 10) {
     return insertion_sort_(x, len);

tests/testthat/test-kendall.R

@@ -6,7 +6,7 @@ test_that("kendall", {
   expect_equal(kendall_cor(x, x), 1)
 
   x <- rep(1, 3)
-  expect_equal(kendall_cor(x, x), cor(x, x, method = "kendall"))
+  expect_warning(kendall_cor(x, x), "zero variance")
 
   x <- c(1, 0, 2)
   y <- c(5, 3, 4)
@@ -34,6 +34,16 @@ test_that("kendall", {
   k2 <- cor.test(x, y, method = "kendall", alternative = "two.sided")
   expect_equal(k1$statistic, unname(k2$estimate))
   expect_equal(k1$p_value, k2$p.value)
+
+  k1 <- kendall_cor_test(x, y, alternative = "greater")
+  k2 <- cor.test(x, y, method = "kendall", alternative = "greater")
+  expect_equal(k1$statistic, unname(k2$estimate))
+  expect_equal(k1$p_value, k2$p.value)
+
+  k1 <- kendall_cor_test(x, y, alternative = "less")
+  k2 <- cor.test(x, y, method = "kendall", alternative = "less")
+  expect_equal(k1$statistic, unname(k2$estimate))
+  expect_equal(k1$p_value, k2$p.value)
 
   x <- rnorm(1e3)
   y <- rpois(1e3, 2)
@@ -48,7 +58,6 @@ test_that("kendall", {
   }
   t_kendall <- median(t_kendall)
 
-
   t_cor <- c()
   for (i in 1:100) {
     t1 <- Sys.time()