Repository: JuliaStats/HypothesisTests.jl
Pull request #146 (2019-01-26 07:25:35)

Add partial correlation test
NOTE: This will have merge conflicts with #145.

Imported from https://github.com/ararslan/MultivariateTests.jl. Requires StatsBase 0.27.0, which introduced `partialcor`, also taken from MultivariateTests. Test values validated against SAS when originally written.

mschauer: I tried to compare with https://github.com/mschauer/CausalInference.jl/blob/f354ca8ce9052293023f7ccefda299bf946dd977/src/skeleton.jl#L105, but could not match up turns, does the atan belong there somewhere or do I have it wrong?
ararslan: > but could not match up turns

I'm not familiar with that terminology. Can you elaborate?
mschauer: Apparently I meant to write: match up the *terms* in both expressions
ararslan: Oh, haha, that makes more sense. I don't know about the use of trig functions there; my original implementation was based on the book and notes from a course on multivariate statistics. The online version of the notes that mentions this test is here: https://newonlinecourses.science.psu.edu/stat505/node/47/. They just formulate the test statistic basically like a regular correlation.
simonster: `(log1p(test.r) - log1p(-test.r)) / 2` is the same as `atanh(test.r)` if that helps.
ararslan: If you can think of a test case in which such catastrophic cancelation could occur, I'd be happy to add it.
mschauer: If you consider writing a test I would rather replace by the difference by `atanh`, because such tests have already been written for `atanh`. I checked
```
f(r) = (log1p(r) - log1p(-r)) / 2
julia> using Test

julia> @test f(-1+eps()) == atanh(-1+eps())
Test Passed
```
Again I don't think anything goes wrong, I just wanted to report I did due diligence.
nalimilan: @ararslan Do you want to finish this? Looks like it's quite close.
ararslan: I think I was blocked on coming up with an appropriate test as Moritz mentioned.
mschauer: I think this would do
ararslan: You mean the PR as-is?
ararslan: Rebased.
ararslan: Looks like there are some numerical differences >0.001.
mschauer: I understand this know, the test and the CI are unrelated... the test is exact using a different test statistic which is (only under the null hypothesis rho=0) t distributed and the CI is approximate using fishers z as statistic which is approximately normal under a range of partial correlation values (this is the story for correlation estimation and it seems that this translates to partial correlations.)
ararslan: That's vaguely what I had thought but I didn't remember well enough to say definitively.
mschauer: I'll merge. One question: this is the right thing to do if you want to do a correlation tests. Should we define an alias emphasising  testing for (non-partial) Pearson correlation coefficients?
ararslan: Hm, we could define it instead as `CorrelationTest` where the two-argument method performs a normal correlation test and the three-argument method does the partial correlation test.
ararslan: I'm not sure when I'll have time to make further changes to this PR, unfortunately.
ararslan: Rebased and renamed to `CorrelationTest`. The two-argument method that performs the "regular" correlation test is not implemented here.
mschauer: @aguang We will merge this without, but then it should be fairly straight forward to implement the two argument version on top of it. @ararslan Thank you for bearing with me, the princess is not in another castle this time.

Patch:
diff --git a/Project.toml b/Project.toml
index 9ca7863..94b9699 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,6 +11,9 @@ Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
+[compat]
+StatsBase = ">=0.27.0"
+
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/REQUIRE b/REQUIRE
index 00891c5..2a4362d 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1,6 +1,6 @@
 julia 0.7
 Distributions 0.16.3
 Roots
-StatsBase 0.25.0
+StatsBase 0.27.0
 Rmath 0.5.0
 Combinatorics 0.7.0
diff --git a/docs/src/multivariate.md b/docs/src/multivariate.md
index 4b587a5..c136da7 100644
--- a/docs/src/multivariate.md
+++ b/docs/src/multivariate.md
@@ -16,3 +16,9 @@ This is equivalent to Box's ``M``-test for two groups.
 ```@docs
 BartlettTest
 ```
+
+## Partial correlation test
+
+```@docs
+CorrelationTest
+```
diff --git a/src/correlation.jl b/src/correlation.jl
new file mode 100644
index 0000000..d40b1b6
--- /dev/null
+++ b/src/correlation.jl
@@ -0,0 +1,78 @@
+# Correlation
+# TODO: Implement `CorrelationTest` for two arguments
+
+using Statistics: clampcor
+
+export CorrelationTest
+
+"""
+    CorrelationTest(x, y, Z)
+
+Perform a t-test for the hypothesis that ``\\text{Cor}(x,y|Z=z) = 0``, i.e. the partial
+correlation of vectors `x` and `y` given the matrix `Z` is zero.
+
+Implements `pvalue` for the t-test.
+Implements `confint` using an approximate confidence interval based on Fisher's
+``z``-transform.
+
+See also `partialcor` from StatsBase.
+
+# External resources
+
+* [Partial correlation on Wikipedia](https://en.wikipedia.org/wiki/Partial_correlation#As_conditional_independence_test)
+  (for the construction of the confidence interval)
+* [Section testing using Student's t-distribution from Pearson correlation coefficient on Wikipedia](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#Testing_using_Student's_t-distribution)
+* [K.J. Levy and S.C. Narula (1978): Testing Hypotheses concerning Partial Correlations: Some Methods and Discussion. International Statistical Review 46(2).](https://www.jstor.org/stable/1402814)
+"""
+struct CorrelationTest{T<:Real} <: HypothesisTest
+    r::T
+    n::Int
+    k::Int
+    t::T
+
+    # Error checking is done in `partialcor`
+    function CorrelationTest(x::AbstractVector, y::AbstractVector, Z::AbstractMatrix)
+        r = partialcor(x, y, Z)
+        n, k = size(Z)
+        t = r * sqrt((n - 2 - k) / (1 - r^2))
+        return new{typeof(r)}(r, n, k, t)
+    end
+
+    function CorrelationTest(x::AbstractVector, y::AbstractVector, z::AbstractVector)
+        r = partialcor(x, y, z)
+        n = length(z)
+        t = r * sqrt((n - 3) / (1 - r^2))
+        return new{typeof(r)}(r, n, 1, t)
+    end
+end
+
+testname(p::CorrelationTest) =
+    string("Test for nonzero ", p.k == 0 ? "partial " : "", " correlation")
+
+function population_param_of_interest(p::CorrelationTest)
+    param = p.k == 0 ? "Partial correlation" : "Correlation"
+    (param, zero(p.r), p.r)
+end
+
+StatsBase.nobs(p::CorrelationTest) = p.n
+StatsBase.dof(p::CorrelationTest) = p.n - 2 - p.k
+
+function StatsBase.confint(test::CorrelationTest{T}, alpha::Float64=0.05) where T
+    dof(test) > 1 || return (-one(T), one(T))  # Otherwise we can get NaNs
+    q = quantile(Normal(), 1 - alpha / 2)
+    fisher = atanh(test.r)
+    bound = q / sqrt(dof(test) - 1)
+    elo = clampcor(tanh(fisher - bound))
+    ehi = clampcor(tanh(fisher + bound))
+    return (elo, ehi)
+end
+
+default_tail(::CorrelationTest) = :both
+pvalue(test::CorrelationTest; tail=:both) = pvalue(TDist(dof(test)), test.t, tail=tail)
+
+function show_params(io::IO, test::CorrelationTest, indent="")
+    println(io, indent, "number of observations:          ", nobs(test))
+    println(io, indent, "number of conditional variables: ", test.k)
+    println(io, indent, "t-statistic:                     ", test.t)
+    println(io, indent, "degrees of freedom:              ", dof(test))
+end


Test patch:
diff --git a/src/HypothesisTests.jl b/src/HypothesisTests.jl
index 12b45a1..fdfcab7 100644
--- a/src/HypothesisTests.jl
+++ b/src/HypothesisTests.jl
@@ -189,4 +189,5 @@ include("hotelling.jl")
 include("bartlett.jl")
 include("wald_wolfowitz.jl")
 include("f.jl")
+include("correlation.jl")
 end
diff --git a/test/correlation.jl b/test/correlation.jl
new file mode 100644
index 0000000..4c27ed2
--- /dev/null
+++ b/test/correlation.jl
@@ -0,0 +1,31 @@
+using HypothesisTests
+using Test
+using DelimitedFiles
+using StatsBase
+
+@testset "Partial correlation" begin
+    # Columns are information, similarities, arithmetic, picture completion
+    wechsler = readdlm(joinpath(@__DIR__, "data", "wechsler.txt"))[:,2:end]
+    w = CorrelationTest(wechsler[:,1], wechsler[:,2], wechsler[:,3:4])
+    let out = sprint(show, w)
+        @test occursin("reject h_0", out) && !occursin("fail to", out)
+    end
+    let ci = confint(w)
+        @test first(ci) ≈ 0.4963917 atol=1e-6
+        @test last(ci) ≈ 0.8447292 atol=1e-6
+    end
+    @test nobs(w) == 37
+    @test dof(w) == 33
+    @test pvalue(w) < 0.00001
+
+    X = [ 2 1 0
+          4 2 0
+         15 3 1
+         20 4 1]
+    x = CorrelationTest(view(X,:,1), view(X,:,2), view(X,:,3))
+    @test occursin("fail to reject", sprint(show, x))
+    @test confint(x) == (-1.0, 1.0)
+    @test nobs(x) == 4
+    @test dof(x) == 1
+    @test pvalue(x) ≈ 0.25776212 atol=1e-6
+end
diff --git a/test/data/wechsler.txt b/test/data/wechsler.txt
new file mode 100644
index 0000000..913211f
--- /dev/null
+++ b/test/data/wechsler.txt
@@ -0,0 +1,37 @@
+ 1  7  5  9  8
+ 2  8  8  5  6
+ 3 16 18 11  9
+ 4  8  3  7  9
+ 5  6  3 13  9
+ 6 11  8 10 10
+ 7 12  7  9  8
+ 8  8 11  9  3
+ 9 14 12 11  4
+10 13 13 13  6
+11 13  9  9  9
+12 13 10 15  7
+13 14 11 12  8
+14 15 11 11 10
+15 13 10 15  9
+16 10  5  8  6
+17 10  3  7  7
+18 17 13 13  7
+19 10  6 10  7
+20 10 10 15  8
+21 14  7 11  5
+22 16 11 12 11
+23 10  7 14  6
+24 10 10  9  6
+25 10  7 10 10
+26  7  6  5  9
+27 15 12 10  6
+28 17 15 15  8
+29 16 13 16  9
+30 13 10 17  8
+31 13 10 17 10
+32 19 12 16 10
+33 19 15 17 11
+34 13 10  7  8
+35 15 11 12  8
+36 16  9 11 11
+37 14 13 14  9
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index ade795d..629a098 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -27,6 +27,7 @@ include("hotelling.jl")
 include("kolmogorov_smirnov.jl")
 include("kruskal_wallis.jl")
 include("mann_whitney.jl")
+include("correlation.jl")
 include("permutation.jl")
 include("power_divergence.jl")
 include("show.jl")
@@ -34,4 +35,4 @@ include("t.jl")
 include("wald_wolfowitz.jl")
 include("wilcoxon.jl")
 include("z.jl")
-include("f.jl")
\ No newline at end of file
+include("f.jl")