SymPy でやった。簡単だった。
またまた HypothesisTests.ChisqTest() の件であるが，Julia は Yates の連続性の補正は眼中にないらしい。いろいろなページにも一切出てこない。
で，以下のように拡張した。
using HypothesisTests
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::AbstractMatrix{T}; correct::Bool=false) where T<:Integer
nrows, ncols = size(x)
if correct && nrows == 2 && ncols == 2
n = sum(x)
a, c, b, d = float.(x)
stat = n*max(0, abs(a*d - b*c) - 0.5n)^2 / ((a+b)*(c+d)*(a+c)*(b+d))
result = PowerDivergenceTest(1, ones(4), stat, 1, x, n, ones(4), x, x, x)
else
result = PowerDivergenceTest(x, lambda=1.0)
end
summary(result)
end
julia> x = [10 3; 4 12]
2×2 Matrix{Int64}:
10 3
4 12
julia> result = ChisqTest(x);
julia> println("chisq. = $(result.stat), df = $(result.df), p.value = $(pvalue(result))")
chisq. = 7.743956043956044, df = 1, p.value = 0.005389260671694261
julia> chisq_test(x)
chisq. = 7.74396, df = 1, p.value = 0.00539
julia> chisq_test(x, correct=true)
chisq. = 5.80415, df = 1, p.value = 0.01599
Julia の HypothesisTests.ChisqTest() は色々問題がある。
今までもいくつか指摘したところであるが，今回は独立性の検定（χ二乗検定）に引数として 2 つのベクトルを与える場合について書く。
ChisqTest() の定義で，引数に 2 変数のベクトルを与える場合は以下の 2 つの関数が定義されている。
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 つの定義だけで済むはず。
それが起こったのは，たぶん 2022 年 1 月末のある日。
以前にはなんのエラーもなく実行できていたプログラムでエラーが生じる。FreqTables.freqtable() で作成したクロス集計表の列や行を抽出できなくなっていた。
再現すると以下のようである。
まず以下のようなデータを読む。
io = IOBuffer("""
ID,Weight,Height,Age,BloodType,Gender
125,44.6,157.6,17,A,Male
321,57.3,158.8,26,B,Male
437,54.3,162.2,22,AB,Female
426,45.4,160.8,31,O,Male
243,48.5,174.5,34,B,Female
""")
CSV.read() の第 1 引数は IOBuffer() の戻り値をとることができる。
using CSV, DataFrames
df = CSV.read(io, DataFrame)
もう気づいてしまった人もいるかもしれないが，先に進む。
df.Gender と df.BloodType のクロス集計表を作る。
using FreqTables
tbl = freqtable(df.Gender, df.BloodType)
2×4 Named Matrix{Int64}
Dim1 ╲ Dim2 │ "A" "AB" "B" "O"
────────────┼───────────────────────
"Female" │ 0 1 1 0
"Male" │ 1 0 1 1
集計表 tbl の 1 行目，2 列目を取り出してみる。
以下はいずれもエラーになる。
MethodError: no method matching indices(::OrderedCollections.OrderedDict{String7, Int64}, ::String)
tbl["Female", :]
tbl[:, "AB"]
以下はなんの問題もない。
tbl[1, :]
4-element Named Vector{Int64}
Dim2 │
──────┼──
"A" │ 0
"AB" │ 1
"B" │ 1
"O" │ 0
tbl[:, 2]
2-element Named Vector{Int64}
Dim1 │
─────────┼──
"Female" │ 1
"Male" │ 0
あれこれやって，半日が過ぎた（大げさ）。
幸い，以前実行した結果が残っているので，ダメ元で今の出力結果と比較した。
気がついた。
BloodType と Gender の型が String3 と String7 になっている。
今までは両方とも String だった。
BloodType と Gender の型を String に変換しよう。
df[!, :BloodType] = string.(df.BloodType)
df[!, :Gender] = string.(df.Gender)
df
これでもう一度集計表を作り，行と列を取り出してみよう。
tbl = freqtable(df.Gender, df.BloodType)
2×4 Named Matrix{Int64}
Dim1 ╲ Dim2 │ A AB B O
────────────┼───────────────
Female │ 0 1 1 0
Male │ 1 0 1 1
tbl["Female", :]
4-element Named Vector{Int64}
Dim2 │
──────┼──
A │ 0
AB │ 1
B │ 1
O │ 0
tbl[1, :]
4-element Named Vector{Int64}
Dim2 │
──────┼──
A │ 0
AB │ 1
B │ 1
O │ 0
tbl[:, "AB"]
2-element Named Vector{Int64}
Dim1 │
───────┼──
Female │ 1
Male │ 0
tbl[:, 2]
2-element Named Vector{Int64}
Dim1 │
───────┼──
Female │ 1
Male │ 0
前と同じようにちゃんと動いている。
そもそも String3 と String7 とはなんぞやと，ぐぐってみると以下のようなことがわかった。
・ 固定長文字列(Fixed width strings) ということで，String1, String3, String7, ..., String(2^8-1) がある。
・ かなり前からあるが，CSV.jl 0.9.1 ではすでに取り込まれていた。2022/02/03 では CSV v0.10.2 だ。
・ パフォーマンスが優れている。
・ まだサポートしていないパッケージがある。まさに freqtable() がそうだった。
そのうちに freqtable() も固定長文字列に対応するのだろうか。
PVアクセスランキング にほんブログ村
