2x2の分割表での独立性検定の比較

黒木玄

2017-09-26

$2\times 2$ の分割表における次のような形の独立性の条件を満たしている確率分布を考える:

$$\displaystyle \begin{bmatrix} p_{11} & p_{12} \\ p_{21} & p_{22} \\ \end{bmatrix} = \begin{bmatrix} pq & p(1-q) \\ (1-p)q & (1-p)(1-q) \\ \end{bmatrix} $$

このとき $p=p_{11}+p_{12}$, $q=p_{11}+p_{21}$ を周辺確率 (marginal probabilities) と呼ぶことにする。

周辺確率を動かしながら、$\chi^2$ 検定、G検定、Fisherの正確確率検定の3種類について、p値が $\alpha$ 以下になる確率をプロットしてみる。

p値が $\alpha$ 以下になる確率は $\alpha$ に近い方がよい。 $\alpha$ に近い場合には白色でプロットし、$\alpha$ より小さい場合には青色をどんどん濃くして行き、$\alpha$ より大きい場合には赤色をどんどん濃くして行くことにする。白っぽく見える方がよい。青色はその検定が conservative であることを意味している。

サンプルを生成する確率分布として、

  • 4×Poisson分布
  • 多項分布
  • 2×二項分布
  • 超幾何分布

を用いる。$\alpha=0.05$ (有意水準5%)の場合を扱うことにする。

結論

  • $\chi^2$ 検定が最も安定して白っぽくなりやすい。

  • G検定は濃い赤が出易い。

  • Fisherの正確確率検定は $n$ が小さいとき濃い青がたくさん出る。

追記: 最後の方に独立性の条件を満たさない場合に関するシミュレーション結果も追加した。

独立性の条件を満たさない場合にはもっと意味のある情報が得られるシミュレーションを考えたい。

In [1]:
using PyPlot
using Distributions
In [2]:
# xlog(x,y) = if x == 0 then 0 else x*log(x/y)
#
# Warning:
#
#   Don't use x*log(x/y) instead of xlog(x,y).
#   Because log(0) = -Inf and 0*Inf = NaN.
#
function xlog(x::Float64, y::Float64)
    if x == zero(x)
        return zero(x)
    elseif y == zero(y)
        return Inf
    else
        return x*log(x/y)
    end
end
xlog(x,y) = xlog(Float64(x),Float64(y))

x, y = 0, 1
@show x, y
@show x*log(x/y)
@show xlog(x,y);
(x, y) = (0, 1)
x * log(x / y) = NaN
xlog(x, y) = 0.0
In [3]:
prodprob(p::Float64, q::Float64) = [p*q, p*(1-q), (1-p)*q, (1-p)*(1-q)]

function expected(a::AbstractArray{T,1}) where T
    n = sum(a)
    return [(a[1]+a[2])*(a[1]+a[3])/n, (a[1]+a[2])*(a[2]+a[4])/n, 
            (a[3]+a[4])*(a[1]+a[3])/n, (a[3]+a[4])*(a[2]+a[4])/n]
end

function chisqtest(a::AbstractArray{Int64,1})
    mu = expected(a)
    chisq = sum((a .- mu).^2 ./mu)
    pval = ccdf(Chisq(1), chisq)
    return pval
end

function gtest(a::AbstractArray{Int64,1})
    mu = expected(a)
    g = 2*sum(xlog.(a,mu)) # Don't use a.*log(a./mu)
    pval = ccdf(Chisq(1), g)
    return pval
end

function pvaluehg(d::Hypergeometric, k::Int64)
    c = params(d)
    amax = min(c[1],c[3])
    p0 = pdf(d, k)
    p1 = 0.0
    pval = 0.0
    for j in 0:amax
        p1 = pdf(d, j)
        pval += ifelse(p1  p0, p1, 0.0)
    end
    return min(pval, 1.0)
end

function fishertest(a::AbstractArray{Int64,1})
    d = Hypergeometric(a[1]+a[2], a[3]+a[4], a[1]+a[3])
    return pvaluehg(d, a[1])
end
Out[3]:
fishertest (generic function with 1 method)
In [4]:
ecdf(pval::AbstractArray{Float64,1}, x::Float64) = count(p -> p  x, pval)/length(pval)
ecdf(pval, x) = ecdf(pval, Float64(x))

function randPoisson(n::Int64, p::AbstractArray{Float64}, N::Int64)
    return [
        rand(Poisson(n*p[1]),N)';
        rand(Poisson(n*p[2]),N)';
        rand(Poisson(n*p[3]),N)';
        rand(Poisson(n*p[4]),N)';
    ]
end

function randMultinomial(n::Int64, p::AbstractArray{Float64}, N::Int64)
    return rand(Multinomial(n,p), N)
end
    
function randBinomial(n::Int64, p::AbstractArray{Float64}, N::Int64)
    m = Int64(round(n*(p[1]+p[2])))
    q1 = p[1]/(p[1]+p[2])
    q2 = p[3]/(p[3]+p[4])
    return [rand(Multinomial(m, [q1, 1.0-q1]), N); rand(Multinomial(n-m, [q2, 1.0-q2]), N)]
end

function randHypergeometric(n::Int64, p::AbstractArray{Float64}, N::Int64)
    m1 = Int64(round(n*(p[1]+p[2])))
    m2 = Int64(round(n*(p[3]+p[4])))
    n1 = Int64(round(n*(p[1]+p[3])))
    a = rand(Hypergeometric(m1, m2, n1), N)'
    b = m1 .- a
    c = n1 .- a
    d = m2 .- c
    return [a; b; c; d]
end

function pvaluesby(sampler::T, n::Int64;
        N = 10^4, alpha = 0.05) where T<:Function
    px = collect(0.05:0.05:0.50)
    py = px
    np = length(px)
    prob_chisq  = Array{Float64,2}(np,np)
    prob_g      = Array{Float64,2}(np,np)
    prob_fisher = Array{Float64,2}(np,np)
    a = Array{Float64,2}(4,N)
    
    for i in 1:np
        for j in 1:np
            if i > j
                prob_chisq[i,j]  = prob_chisq[j,i]
                prob_g[i,j]      = prob_g[j,i]
                prob_fisher[i,j] = prob_fisher[j,i]
            else
                a = sampler(n, prodprob(px[i],py[j]), N)
                prob_chisq[i,j]  = ecdf([chisqtest(a[:,i])  for i in 1:size(a,2)], alpha)
                prob_g[i,j]      = ecdf([gtest(a[:,i])      for i in 1:size(a,2)], alpha)
                prob_fisher[i,j] = ecdf([fishertest(a[:,i]) for i in 1:size(a,2)], alpha)
            end
        end
    end
    return alpha, px, prob_chisq, prob_g, prob_fisher
end

function plotcomparisontest(sampler, n; N = 10^4, alpha = 0.05)
    alpha, px, prob_chisq, prob_g, prob_fisher = pvaluesby(sampler, n; N = N, alpha = alpha)
    py = px
    np = length(px)
    ps = [0.0;px]
    vmin = 0.0
    vmax = 2*alpha
    cmap = "RdBu_r"

    figure(figsize=(8,6.4))

    ax1 = subplot2grid((16,20), (0,0), rowspan=7, colspan=8)
    pcolormesh(ps, ps, prob_chisq, cmap=cmap, vmin=vmin, vmax=vmax)
    colorbar(label="probability of p-value \$\\leqq\$ $alpha")
    title("\$\\chi^2\$-test")
    xlabel("marginal probability", fontsize=8)
    ylabel("marginal probability", fontsize=8)

    ax2 = subplot2grid((16,20), (0,12), rowspan=7, colspan=8)
    pcolormesh(ps, ps, prob_g, cmap=cmap, vmin=vmin, vmax=vmax)
    colorbar(label="probability of p-value \$\\leqq\$ $alpha")
    title("G-test")
    xlabel("marginal probability", fontsize=8)
    ylabel("marginal probability", fontsize=8)

    ax3 = subplot2grid((16,20), (9,0), rowspan=7, colspan=8)
    pcolormesh(ps, ps, prob_fisher, cmap=cmap, vmin=vmin, vmax=vmax)
    colorbar(label="probability of p-value \$\\leqq\$ $alpha")
    title("Fisher's exact test")
    xlabel("marginal probability", fontsize=8)
    ylabel("marginal probability", fontsize=8)

    suptitle("sampler = $(typeof(sampler)),  n = $n,  \$\\alpha\$ = $alpha")
end
Out[4]:
plotcomparisontest (generic function with 1 method)

4×Poisson分布でサンプルを生成する場合

In [5]:
@time plotcomparisontest(randPoisson, 25, N=10000, alpha=0.05)
 
 5.861474 seconds (15.16 M allocations: 706.183 MiB, 1.42% gc time)
Out[5]:
PyObject <matplotlib.text.Text object at 0x000000003281B5F8>
In [6]:
@time plotcomparisontest(randPoisson, 50, N=10000, alpha=0.05)
 
 5.360726 seconds (13.12 M allocations: 599.222 MiB, 0.84% gc time)
Out[6]:
PyObject <matplotlib.text.Text object at 0x0000000032C3A5F8>
In [7]:
@time plotcomparisontest(randPoisson, 100, N=10000, alpha=0.05)
 
 9.175158 seconds (13.12 M allocations: 599.223 MiB, 0.34% gc time)
Out[7]:
PyObject <matplotlib.text.Text object at 0x00000000332EE7F0>
In [8]:
@time plotcomparisontest(randPoisson, 200, N=10000, alpha=0.05)
 
16.607463 seconds (13.12 M allocations: 599.223 MiB, 0.20% gc time)
Out[8]:
PyObject <matplotlib.text.Text object at 0x00000000337B84E0>

多項分布でサンプルを生成する場合

In [9]:
@time plotcomparisontest(randMultinomial, 25, N=10000, alpha=0.05)
 
 3.602437 seconds (13.85 M allocations: 617.143 MiB, 1.24% gc time)
Out[9]:
PyObject <matplotlib.text.Text object at 0x0000000033C20E80>
In [10]:
@time plotcomparisontest(randMultinomial, 50, N=10000, alpha=0.05)
 
 5.340411 seconds (13.67 M allocations: 607.617 MiB, 0.76% gc time)
Out[10]:
PyObject <matplotlib.text.Text object at 0x00000000341464A8>
In [11]:
@time plotcomparisontest(randMultinomial, 100, N=10000, alpha=0.05)
 
 9.183502 seconds (13.67 M allocations: 607.614 MiB, 0.43% gc time)
Out[11]:
PyObject <matplotlib.text.Text object at 0x00000000348B85F8>
In [12]:
@time plotcomparisontest(randMultinomial, 200, N=10000, alpha=0.05)
 
16.333810 seconds (13.67 M allocations: 607.614 MiB, 0.24% gc time)
Out[12]:
PyObject <matplotlib.text.Text object at 0x00000000318B2668>

2×二項分布でサンプルを生成する場合

In [13]:
@time plotcomparisontest(randBinomial, 25, N=10000, alpha=0.05)
 
 3.455279 seconds (14.48 M allocations: 655.758 MiB, 1.54% gc time)
Out[13]:
PyObject <matplotlib.text.Text object at 0x0000000031D944E0>
In [14]:
@time plotcomparisontest(randBinomial, 50, N=10000, alpha=0.05)
 
 5.334195 seconds (14.22 M allocations: 649.629 MiB, 0.89% gc time)
Out[14]:
PyObject <matplotlib.text.Text object at 0x00000000356A3A58>
In [15]:
@time plotcomparisontest(randBinomial, 100, N=10000, alpha=0.05)
 
 9.082301 seconds (14.22 M allocations: 649.629 MiB, 0.54% gc time)
Out[15]:
PyObject <matplotlib.text.Text object at 0x000000003ABBF4A8>
In [16]:
@time plotcomparisontest(randBinomial, 200, N=10000, alpha=0.05)
 
16.318503 seconds (14.22 M allocations: 649.629 MiB, 0.31% gc time)
Out[16]:
PyObject <matplotlib.text.Text object at 0x000000003AC6BC18>

超幾何分布でサンプルを生成する場合

In [17]:
@time plotcomparisontest(randHypergeometric, 25, N=10000, alpha=0.05)
 
 3.441200 seconds (13.19 M allocations: 615.505 MiB, 1.23% gc time)
Out[17]:
PyObject <matplotlib.text.Text object at 0x000000003C671828>
In [18]:
@time plotcomparisontest(randHypergeometric, 50, N=10000, alpha=0.05)
 
 5.225723 seconds (13.12 M allocations: 611.826 MiB, 0.78% gc time)
Out[18]:
PyObject <matplotlib.text.Text object at 0x000000003CC134A8>
In [19]:
@time plotcomparisontest(randHypergeometric, 100, N=10000, alpha=0.05)