Merge c01968d into 7459957

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

docs/make.jl

@@ -9,7 +9,8 @@ makedocs(
         "methods.md",
         "parametric.md",
         "nonparametric.md",
-        "time_series.md"
+        "time_series.md",
+        "multivariate.md",
     ]
 )
 

docs/src/multivariate.md

@@ -0,0 +1,7 @@
+# Multivariate tests
+
+## Partial correlation test
+
+```@docs
+PartialCorTest
+```

src/HypothesisTests.jl

@@ -177,4 +177,5 @@ include("augmented_dickey_fuller.jl")
 include("jarque_bera.jl")
 include("durbin_watson.jl")
 include("permutation.jl")
+include("partialcor.jl")
 end

src/partialcor.jl

@@ -0,0 +1,61 @@
+# Partial correlation
+
+export PartialCorTest
+
+"""
+    PartialCorTest(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` and `confint`. See also [`partialcor`](@ref).
+"""
+struct PartialCorTest <: HypothesisTest
+    r::Real
+    n::Int
+    k::Int
+    t::Real
+
+    # Error checking is done in `partialcor`
+    function PartialCorTest(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(r, n, k, t)
+    end
+
+    function PartialCorTest(x::AbstractVector, y::AbstractVector, z::AbstractVector)
+        r = partialcor(x, y, z)
+        n = length(z)
+        t = r * sqrt((n - 3) / (1 - r^2))
+        return new(r, n, 1, t)
+    end
+end
+
+testname(::PartialCorTest) = "Test for partial correlation"
+population_param_of_interest(p::PartialCorTest) = ("Partial correlation", zero(p.r), p.r)
+
+StatsBase.nobs(p::PartialCorTest) = p.n
+StatsBase.dof(p::PartialCorTest) = p.n - 2 - p.k
+
+function StatsBase.confint(test::PartialCorTest, alpha::Float64=0.05)
+    dof(test) > 1 || return (-one(test.r), one(test.r)) # Otherwise we can get NaNs
+    q = quantile(Normal(), 1 - alpha / 2)
+    fisher = (log1p(test.r) - log1p(-test.r)) / 2
+    bound = q / sqrt(test.n - 3 - test.k)
+    lo = 2 * (fisher - bound)
+    hi = 2 * (fisher + bound)
+    elo = Statistics.clampcor(expm1(lo) / (exp(lo) + 1))
+    ehi = Statistics.clampcor(expm1(hi) / (exp(hi) + 1))
+    return (elo, ehi)
+end
+
+default_tail(::PartialCorTest) = :both
+pvalue(test::PartialCorTest; tail=:both) = pvalue(TDist(dof(test)), test.t, tail=tail)
+
+function show_params(io::IO, test::PartialCorTest, 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/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

test/partialcor.jl

@@ -0,0 +1,32 @@
+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 = PartialCorTest(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 = PartialCorTest(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
+

test/runtests.jl

@@ -31,3 +31,4 @@ include("show.jl")
 include("t.jl")
 include("wilcoxon.jl")
 include("z.jl")
+include("partialcor.jl")