Julia: HypothesisTests.ChisqTest() の問題点

Julia の HypothesisTests.ChisqTest() は色々問題がある。

今までもいくつか指摘したところであるが，今回は独立性の検定（χ二乗検定）に引数として 2 つのベクトルを与える場合について書く。

ChisqTest() の定義で，引数に 2 変数のベクトルを与える場合は以下の 2 つの関数が定義されている。

https://github.com/JuliaStats/HypothesisTests.jl/blob/413a901612d5ca9551fe3e28a9a0646dd4b3c98d/src/power_divergence.jl#L349-L370

function ChisqTest(x::AbstractVector{T}, y::AbstractVector{T}, levels::Levels{T}) where T<:Integer
    d = counts(x, y, levels)
    PowerDivergenceTest(d, lambda=1.0)
end

function ChisqTest(x::AbstractVector{T}, y::AbstractVector{T}, k::T) where T<:Integer
    d = counts(x, y, k)
    PowerDivergenceTest(d, lambda=1.0)
end

第 1 の問題は，引数で与えることができるのは整数ベクトルに限られるということである。

文字列ベクトルの場合も多いと思うが，一旦整数ベクトルに直してから ChisqTest() を呼ぶとか，freqtable() で集計した結果を ChisqTest() に渡すしかない。

しかし，もっと問題なのは levels とか k という意味不明の引数である。

Julia の常ではあるが，オンラインヘルプはおろか，ネット上にも説明がないし，使用例もない。

内部を見ると，集計表を作るのは counts() であるが，この関数のオンラインヘルプも，そもそも二重集計の例については触れられていない。

試行錯誤でやってみると，x, y の値のうち，第 3 引数で示した数値以下のものについて集計するようだ。
第 3 引数は，x, y の両方に適用されるので，カテゴリー数の違う変数の場合は　ChisqTest() は適用できなくなってしまう。

まず，x, y ともに 4 つのカテゴリー（1,2,3,4）を持つ場合にやってみよう。

using HypothesisTests
using FreqTables

using Random; Random.seed!(123)
x = rand(1:4, 100)
y = rand(1:4, 100)
result = ChisqTest(x, y, 4)
println("chisq. = $(result.stat), df = $(result.df), p.value = $(pvalue(result))")

chisq. = 7.673743524028154, df = 9, p.value = 0.5673295593772689

ここで，自前の関数を定義しておく。

2 つの整数ベクトルを引数に取る chisq_test() と，結果の出力をする summary() を定義しておく。

function summary(result)
    println("chisq. = $(round(result.stat, digits=5)), df = $(result.df), p.value = $(round(pvalue(result), digits=5))")
end

function chisq_test(x::AbstractVector{T}, y::AbstractVector{T}) where T<:Integer
    d = freqtable(x, y)
    result = PowerDivergenceTest(d, lambda=1.0)
    summary(result)
end

chisq_test() の結果は，ChisqTest() と同じになる。

chisq_test(x, y)

chisq. = 7.67374, df = 9, p.value = 0.56733

次に，第 1 引数はカテゴリー数が 3 個，第 2 引数はカテゴリー数が 4 個になる整数ベクトルを与える場合を検討する。

まず，正しい結果を，ChisqTest() に freqtables() で得られた集計表を与えることで求める。

using Random; Random.seed!(456)
x = rand(1:3, 100)
y = rand(1:4, 100)
result = ChisqTest(freqtable(x, y))
println("chisq. = $(result.stat), df = $(result.df), p.value = $(pvalue(result))")

chisq. = 5.209486804755349, df = 6, p.value = 0.5172394227637924

次に，ChisqTest() に 2 つのベクトルを与える場合を見てみる。 第 3 引数に 3, 4 のいずれを与えても（もちろん，その他の値を与えても）正しい答えは得られない。

自由度（df）の数値を見ると，(第 3 引数の数値 - 1)^2 なので，PowerDivergenceTest() を呼出すときの集計表 d = counts(x, y, levels) は levels × levels の正方行列になっていることがわかる。当然，それは不適切である。

result = ChisqTest(x, y, 3)
println("chisq. = $(result.stat), df = $(result.df), p.value = $(pvalue(result))")

chisq. = 3.4795508838987104, df = 4, p.value = 0.480994482228378

result = ChisqTest(x, y, 4)
println("chisq. = $(result.stat), df = $(result.df), p.value = $(pvalue(result))")

chisq. = NaN, df = 9, p.value = NaN

chisq_test() は正しい結果を返す。

関数定義は以下の通り。

function chisq_test(x::AbstractVector{T}, y::AbstractVector{T}) where T<:Integer
    d = freqtable(x, y)
    result = PowerDivergenceTest(d, lambda=1.0)
    summary(result)
end

chisq_test(x, y)

chisq. = 5.20949, df = 6, p.value = 0.51724

次に，2 個の文字列ベクトルを与える場合を見てみる。

ChisqTses() は，文字列ベクトルを扱えない。事前に freqtable() で集計した結果で呼ぶ必要がある。

x = ["Yes", "Yes", "No", "No", "No", "Yes", "No", "Yes", "Yes", "No",
     "No", "Yes", "Yes", "No", "Yes", "Yes", "No", "No", "No", "No"]
y = ["Hi", "Lo", "Med", "Med", "Lo", "Hi", "Lo", "Hi", "Lo", "Lo",
     "Med", "Med", "Lo", "Med", "Med", "Lo", "Lo", "Hi", "Lo", "Med"]
tbl = freqtable(x, y)
tbl

2×3 Named Matrix{Int64}
Dim1 ╲ Dim2 │  Hi   Lo  Med
────────────┼──────────────
No          │   1    5    5
Yes         │   3    4    2

result = ChisqTest(tbl)
println("chisq. = $(result.stat), df = $(result.df), p.value = $(pvalue(result))")

chisq. = 2.2190155523488855, df = 2, p.value = 0.32972121777778757

chisq_test(x, y) は問題無く動く。

関数定義は以下の通り。

function chisq_test(x::AbstractVector{T}, y::AbstractVector{T}) where T<:String
    d = freqtable(x, y)
    result = PowerDivergenceTest(d, lambda=1.0)
    summary(result)
end

chisq_test(x, y)

chisq. = 2.21902, df = 2, p.value = 0.32972

どちらかが整数ベクトル，もう一方が文字列ベクトルの場合も，以下の関数定義でカバーできる。

function chisq_test(x::AbstractVector{T}, y::AbstractVector{U}) where {T<:String,U<:Integer}
    d = freqtable(x, y)
    result = PowerDivergenceTest(d, lambda=1.0)
    summary(result)
end

function chisq_test(x::AbstractVector{T}, theta0::Vector{U} = ones(length(x))/length(x)) where {T<:Integer,U<:AbstractFloat}
    result = PowerDivergenceTest(reshape(x, length(x), 1), lambda=1.0, theta0=theta0)
    summary(result)
end

using Random; Random.seed!(789)
x = rand(["foo", "bar", "baz"], 100)
y = rand(1:4, 100)

result = ChisqTest(freqtable(x, y))
println("chisq. = $(result.stat), df = $(result.df), p.value = $(pvalue(result))")

chisq. = 8.798699647157212, df = 6, p.value = 0.18521956943031795

chisq_test(x, y)

chisq. = 8.7987, df = 6, p.value = 0.18522

chisq_test(y, x)

chisq. = 8.7987, df = 6, p.value = 0.18522

2 つのベクトルを与える chisq_test() を 4 通り定義したのだが，1 つの定義だけで済むはず。