行列を対象にしたカイ２乗検定

実験で質的データ同士の相関を見る必要があり、カイ２乗検定のプログラムを探していました。ベクトルを対象にしたもの（scipy.stats.chisquare）はすぐに見つかりましたが、一番使うであろう行列を対象にしたカイ２乗検定のプログラムが見つからなかったのでPythonで書いてみました。カイ２乗検定についての説明はハンバーガー統計学にようこそ！が分かりやすいのでそちらを見てください。

test_chisquare.pyの例では、

435	165
265	135

という行列にカイ二乗検定を行ないました。結果、カイ二乗値 (ch2) は4.46、カイ２乗値が4.46よりも大きい確率 (pval) は0.03になりました。つまり有意水準5%で2つの母数の間に有意差があることが言えます。

chisquare.py

# -*- coding: utf-8 -*-
import scipy as sp
from scipy import stats

def chisquare_test(array):
    '''行列に対してカイ二乗検定を行なう'''
    (num_row, num_col) = sp.shape(array)
    df = (num_row-1) * (num_col-1)

    ch2 = get_chisquare(array)
    pval = 1 - stats.chi2.cdf(ch2, df)
    return ch2, pval

def get_chisquare(array):
    '''カイ二乗値を計算する'''
    sum_cols = array.sum(axis=0)
    sum_rows = array.sum(axis=1)
    N        = array.sum()
    
    chi2 = 0.0
    (num_row, num_col) = sp.shape(array)
    for i in xrange(num_row):
        for j in xrange(num_col):
            observation = array[i,j]
            expectation = 1.0 * sum_rows[i] * sum_cols[j] / N
            chi2 += 1.0 * (observation - expectation) ** 2 / expectation
    return chi2

test_chisquare.py

# -*- coding: utf-8 -*-
from nose.tools import *
import chisquare
import scipy as sp

def test_chisquare_test():
    array = sp.array([[435, 165],
                      [265, 135]])
    (ch2, pval) = chisquare.chisquare_test(array)
    assert_almost_equal(ch2, 4.46, 2)
    assert_almost_equal(pval, 0.03, 2)

if __name__ == "__main__":
    test_chisquare_test()