Repository: JuliaStats/HypothesisTests.jl
Pull request #334 (2025-04-01 11:37:35)

Replace Rmath with StatsFuns
Closes #331. Closes #330.

devmotion: @ArnoStrouwen maybe you could review this PR? You're surely the person most familiar with `StatsFuns.wilcoxcdf` and `StatsFuns.signrankcdf` 🙂 
devmotion: I added more tests in https://github.com/JuliaStats/HypothesisTests.jl/pull/335 but I haven't been able to validate them with R.
andreasnoack: @devmotion would you be able to resolve the conflicts here? Then it should be good to go, no?

Patch:
diff --git a/Project.toml b/Project.toml
index cfa5ac6..0872dac 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,22 +8,22 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Rmath = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 
 [compat]
 Combinatorics = "1"
 DelimitedFiles = "<0.0.1, 1"
 Distributions = "0.25.90"
-Rmath = "0.5, 0.6, 0.7, 0.8"
 Roots = "0.7, 0.8, 1.0, 2"
 StableRNGs = "1"
 Statistics = "1"
 StatsAPI = "1.6"
 StatsBase = "0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34"
+StatsFuns = "1.4"
 Test = "<0.0.1, 1"
 julia = "1.10"
 
diff --git a/src/mann_whitney.jl b/src/mann_whitney.jl
index d1be2b3..9899d98 100644
--- a/src/mann_whitney.jl
+++ b/src/mann_whitney.jl
@@ -93,8 +93,8 @@ from the same population as `y` is greater than an observation drawn from the sa
 population as `x` against the alternative hypothesis that these probabilities are not
 equal.
 
-When there are no tied ranks, the exact p-value is computed using the `pwilcox` function
-from the `Rmath` package. In the presence of tied ranks, a p-value is computed by exhaustive
+When there are no tied ranks, the exact p-value is computed using the `wilcoxcdf` and `wilcoxccdf`
+functions from the `StatsFuns` package. In the presence of tied ranks, a p-value is computed by exhaustive
 enumeration of permutations, which can be very slow for even moderately sized data sets.
 
 Implements: [`pvalue`](@ref)
@@ -152,15 +152,15 @@ function StatsAPI.pvalue(x::ExactMannWhitneyUTest; tail=:both)
     check_tail(tail)
 
     if x.tie_adjustment == 0
-        # Compute exact p-value using method from Rmath, which is fast but
+        # Compute exact p-value using method from StatsFuns, which is fast but
         # cannot account for ties
         if tail == :both
-            p = pwilcox(min(x.U, x.nx * x.ny - x.U), x.nx, x.ny, true)
+            p = wilcoxcdf(x.nx, x.ny, min(x.U, x.nx * x.ny - x.U))
             min(2 * p, 1.0)
         elseif tail == :left
-            pwilcox(x.U, x.nx, x.ny, true)
+            wilcoxcdf(x.nx, x.ny, x.U)
         else # tail == :right
-            pwilcox(x.U - 1, x.nx, x.ny, false)
+            wilcoxccdf(x.nx, x.ny, x.U - 1)
         end
     else
         # Compute exact p-value by enumerating possible ranks in the tied data
diff --git a/src/wilcoxon.jl b/src/wilcoxon.jl
index f3e5f79..e122dc7 100644
--- a/src/wilcoxon.jl
+++ b/src/wilcoxon.jl
@@ -85,8 +85,8 @@ Perform a Wilcoxon exact signed rank U test of the null hypothesis that the dist
 `x` (or the difference `x - y` if `y` is provided) has zero median against the alternative
 hypothesis that the median is non-zero.
 
-When there are no tied ranks, the exact p-value is computed using the `psignrank` function
-from the `Rmath` package. In the presence of tied ranks, a p-value is computed by exhaustive
+When there are no tied ranks, the exact p-value is computed using the `signrankcdf` and `signrankccdf`
+functions from the `StatsFuns` package. In the presence of tied ranks, a p-value is computed by exhaustive
 enumeration of permutations, which can be very slow for even moderately sized data sets.
 
 Implements: [`pvalue`](@ref), [`confint`](@ref)
@@ -139,17 +139,17 @@ function StatsAPI.pvalue(x::ExactSignedRankTest; tail=:both)
     if n == 0
         1.0
     elseif x.tie_adjustment == 0
-        # Compute exact p-value using method from Rmath, which is fast but cannot account for ties
+        # Compute exact p-value using method from StatsFuns, which is fast but cannot account for ties
         if tail == :both
             if x.W <= n * (n + 1)/4
-                2 * psignrank(x.W, n, true)
+                2 * signrankcdf(n, x.W)
             else
-                2 * psignrank(x.W - 1, n, false)
+                2 * signrankccdf(n, x.W - 1)
             end
         elseif tail == :left
-            psignrank(x.W, n, true)
+            signrankcdf(n, x.W)
         else
-            psignrank(x.W - 1, n, false)
+            signrankccdf(n, x.W - 1)
         end
     else
         # Compute exact p-value by enumerating all possible ranks in the tied data
@@ -253,7 +253,7 @@ function calculate_ci(x::AbstractVector, level::Real=0.95; tail=:both)
     n = length(x)
     m = div(n * (n + 1), 2)
     k_range = 1:div(m, 2)
-    l = [1 - 2 * psignrank(i, n, true) for i in k_range]
+    l = [1 - 2 * signrankcdf(n, i) for i in k_range]
     k = argmin(abs.(l .- c))
     vals = Float64[]
     enumerated = enumerate(x)


Test patch:
diff --git a/src/HypothesisTests.jl b/src/HypothesisTests.jl
index cc33ed8..feefd49 100644
--- a/src/HypothesisTests.jl
+++ b/src/HypothesisTests.jl
@@ -27,7 +27,7 @@ module HypothesisTests
 using Statistics, Random, LinearAlgebra
 using Distributions, Roots, StatsBase
 using Combinatorics: combinations, permutations
-using Rmath: pwilcox, psignrank
+using StatsFuns: wilcoxcdf, wilcoxccdf, signrankcdf, signrankccdf
 using Printf: @printf
 
 import StatsAPI