pythonのprastクラスによる集計（３）　カイ２乗検定

　今回はクロス集計の次ということで、カイ２乗検定を取り上げます。

　処理するcsvファイルは第１回・第２回と同じ data01.csv です。

　当Webページで紹介するスクリプトや素材データ一式は、
prast03.zip という圧縮ファイルに同梱しておきます。

《このページの目次》

• はじめに
• 1.　scipyライブラリの chi2_contingency() によるカイ２乗検定
• 2.　Prastクラスの chisq_test() によるカイ２乗検定
• 3.　スクリプトの例
• 4.　コクランの規則とwarning
• 5.　ごく少数の行または列を対象外にする
  • (1)　別名割当てを変更する
  • (2)　欠損値として設定する
• 6.　クラスに属さない chisq_test()

はじめに

　カイ２乗検定は、２つの事象の間に有意な関連性があるかどうかを検証するものです。

　検定の素材としてクロス集計表を与えます。

　性別×意見のクロス集計表を与えた場合、
この２つの事象が独立である(互いに影響しあう訳ではない)
それとも関連性がある(影響しあう)のどちらなのかを見ます。

　要するに、男女で意見が違うかどうかをみます。

　「カイ」というのはギリシャ文字のχ(小文字)のことですが、
ここではカナで記します。

　χは、アルファベットでいうと x(エックス)に該当するようです。

　カイ２乗検定の結果、２つの事象に関連性があるということになった場合、
表の中のどのセルが有意に大きな値が(あるいは小さな値か)を見るのに
残差分析を行います。

　そこで、Prastクラスの chisq_test() では、残差関連の情報も返すようにしました。

目次に戻る

1.　scipyライブラリの chi2_contingency() によるカイ２乗検定

　pythonのscipyライブラリには、カイ２乗検定を行うための
chi2_contingency() が用意されています。

　戻り値として統計量、有意確率、自由度、期待度数を得ることができます。

　ですが、残差に関する表が得られません。

　scipyに用意されているメソッドの呼び出し方は次のとおり。

import numpy as np
import scipy.stats as stats
chi2, p, df, expected = stats.chi2_contingency(observed)

　引数として与えるobserved(観測度数)は、numpy の array(２次元配列)です。

　戻り値は次のとおり。

• chi2: 統計量(カイ２乗値)
• p: 有意確率。これが0.05より小さいと「２つの事象に関連性あり」と見る。
• df: 自由度
• expected: 期待度数。numpyのarray(２次元配列)。
  　「仮に関連性がないとしたら観測度数がこの値になるはず」という理論上の値。

　観測度数と期待度数のズレをみて、「２つの事象が独立である」という仮説を
棄却できるほどのズレなのか、それほどでもないのかを検証するのがカイ２乗検定。

　棄却できる(関連性がある)と判断される場合は残差分析を行います。

　observed, expected の２つを材料として
残差情報を次のように求めることができます。

res = observed - expected  # 残差
residuals = res / np.sqrt(expected)  # 標準化残差
tsum = observed.sum()*1.0  # 観測度数の計
relm = 1.0 - np.sum(observed, axis=1) / tsum  # 行要素
celm = 1.0 - np.sum(observed, axis=0) / tsum  # 列要素
mouter = np.outer(relm, celm)  # 行要素×列要素：行列の外積
stdres = residuals / np.sqrt(mouter)  # 調整済み残差

　residuals(標準化残差), stdres(調整済み残差)は
どちらも numpyのarray(２次元配列)です。

　カイ２乗検定で２つの事象が独立ではない(何らかの関連性がある)となった場合
最後に求めた調整済み残差が分析の重要な手がかりになります。

　なお、２つの事象の連関の度合いを示すクラメール係数は次のように算出。

cramer = np.sqrt(chi2 / (tsum*(min(observed.shape)-1)))

　クラメール係数は 0〜1 の値をとり、連関が強いほど大きな値になります。

目次に戻る

2.　Prastクラスの chisq_test() によるカイ２乗検定

　Prastクラスに chisq_test() を設けました。

　カイ２乗検定を行い、残差関係の情報とクラメール係数も込みで返します。

Prast.chisq_test(gkey1, gkey2)
dod = psx.chisq_test("gender", "opinion")

　引数には２つの列ラベルを指定します(２つとも必須)。

　chisq_test("gender", "opinion") とすれば、
性別と意見のクロス集計表(合計欄がないもの)を内部で生成し、
それをカイ２乗検定にかけます。

　戻り値のdod(OrderedDict)のkey, valueは下のとおり。

• name: ２つの列ラベル→　性別/意見
• chi2: カイ２乗値→　2.98163830557
• df: 自由度→　2
• p_value: 有意確率→　0.225188116494
• observed: 観測度数、DataFrame　素材としたクロス集計表
• expected: 期待度数、DataFrame
• residuals: 標準化残差、DataFrame
• stdres: 調整済み残差、DataFrame
• cramer: クラメール係数→　0.0871026171678

　今回、p_valueが0.05より大きいので有意性が認められません。

　つまり、性別と意見の間に有意な関連性があるとはいえないと判断します。

　この場合、残差分析は行わず、調整済み残差は参照しません。

　戻り値に含まれる DataFrame の小数点以下の桁数ですが、
cround (Prastクラスのメンバー変数)の値で調整することはできません。

　桁数を変更したいときは、それぞれに pythonのround()を適用します。

dod = psx.chisq_test("gender", "opinion")
for key in ['expected', 'residuals', 'stdres']:
    dod[key] = dod[key].round(3)

　上のようにすると、３つの DataFrame の小数点以下の桁数が３桁以内になります。

目次に戻る

3.　スクリプトの例

　カイ２乗検定のサンプルスクリプトを掲げます。

　地域×意見を検定にかけて、その結果を検証します。

　以下、スクリプトです。

 # chisq01.py (coding: cp932)
 import pandas as pd
 from prast import Prast
 
 dtf = pd.read_csv("data01.csv")
 psx = Prast(dtf, "data01_c.txt", "data01_i.txt", "cp932")
 dod = psx.chisq_test("area", "opinion")
 for key in ['expected', 'residuals', 'stdres']:
     dod[key] = dod[key].round(3)  # 小数点以下の桁数を調整
for key in dod.keys():
    print(key)
    print(dod[key])
    print('')

　上のスクリプトを実行すると、p_value が 9.85593038005e-07 と算出されます。

　これは 0.05 より小さいので、地域と意見の間に何らかの関連性があると見ます。

　そこで、調整済み残差を確認すると下のとおり。

　有意水準を５％とした場合、各セルの値 v が
v＜-1.96 なら該当のセルが有意に小さい
v＞1.96 なら有意に大きいということになります。

　有意水準を１％とした場合は 1.96 を 2.58 に置き換えます。

　上の調整済み残差をみると、
「海岸」では賛成が有意に多く、反対が有意に少ない。
「川沿い」では賛成が有意に少なく、反対が有意に多い。
ということが分かります。

　その他のセルでは有意性が認められません。

目次に戻る

4.　コクランの規則とwarning

　カイ２乗検定を行う場合、期待度数の表をみて
５未満のセルの個数が全体の20％未満でないと
適切な検定結果が得られない、という考え方があります。

　これをコクランの規則というようです。これには異論もあるようですが……

　data01.csv では欠損値がごく少数なので
欠損値を組み入れたクロス集計表を chisq_test() にかけると
warningが標準出力に出力されます。

　性別×意見の場合、次のようなメッセージになります。

警告: '期待度数が5未満のセルの割合' >= 20% | 50.0%(6/12)

　セルの総数 12個のうち、６個が期待度数５未満であり、
その割合は 50％となり、20％以上になっているよ、という警告です。

　chisq_test() の戻り値 dod でも dod['warning'] が設けられ
その値として上記の警告メッセージを得ることができます。

　dod['warning'] がなければコクランの規則に抵触しないということになります。

目次に戻る

5.　ごく少数の行または列を対象外にする

　コクランの規則に抵触しないようにするためには
セルの値がごく少数の行または列を対象外とすることになるでしょうか。

　今回「地域×意見」では欠損値を含めなければ抵触しませんが、
あえて「C:海岸」と「R:川沿い」を対象外にしてみます。

　方法は２つです。

• 別名割当てから C, R を外す。
• C, R を欠損値にしてしまう。

(1)　別名割当てを変更する

　Prastクラスのメソッドで集計する場合、
別名が定義されているなら、定義にないものは対象外になります。

　別名割当てファイル(data01_i.txt)を書き換えてもいいのですが、
一時的に外すだけなら put_indexes() で定義し直す方がいいとおもいます。

　具体的には下のようなスクリプトです。

 # chisq02.py (coding: cp932)
 import pandas as pd
 from prast import Prast
 
 dtf = pd.read_csv("data01.csv")
 psx = Prast(dtf, "data01_c.txt", "data01_i.txt", "cp932")
 iod_org = psx.get_indexes()  # 念のためオリジナルのindexesを保持
 psx.put_indexes('''
 area
H,丘陵
L,湖畔
U,都市部
''')

dod = psx.chisq_test("area", "opinion")
for key in dod.keys():
    print(key)
    print(dod[key])
    print('')
psx.put_indexes(iod_org)  # indexesを元に戻す

(2)　欠損値として設定する

　もう１つの方法として、対象外にしたいものを欠損値にする場合は
recona() を使います。

（前略）
psx = Prast(dtf, "data01_c.txt", "data01_i.txt", "cp932")
for x in ['C', 'R']:
    psx.recona("area", x)
dod = psx.chisq_test("area", "opinion")
（後略）

　この方法の場合、欠損値を元に戻すのが少しやっかいです。

　以前は C だったものと R だったものが一律に NaN になってしまうからです。

　元に戻したいときは、欠損値に置き換える前の area の列を記録しておいて
検定処理終了後にその列を元に戻すということになるでしょうか。

目次に戻る

6.　クラスに属さない chisq_test()

　Prastクラスに属さない chisq_test() も定義してあります。

import prast as ps
dod = ps.chisq_test(dtf)

　上のように引数(DataFrame)を１つだけ指定します。

　既存のクロス集計表があって、それをカイ２乗検定にかけたいときは
こちらの chisq_test() の方を使います。

　引数のクロス集計表は、pandasのDataFrameを与えます。
numpyのarrayではありません。

　行ラベルも列ラベルもない、数字だけからなる行列を検定する例を掲げます。

 # chisq03.py (coding: cp932)
 import os
 import pandas as pd
 import prast as ps
 
 csv_data = '''\
 17,1,5
 4,6,14
 7,13,9
8,9,2
'''
tfile = "tmp.csv"  # 一時的に設けるファイルの名前
with open(tfile, "w") as f:
    f.write(csv_data)
dtf = pd.read_csv(tfile, header=None)
os.remove(tfile)  # 一時的ファイルを削除
dod = ps.chisq_test(dtf)
for key in dod.keys():
    print(key)
    print(dod[key])
    print('')

　行列をcsvファイルとして書き出し、
それを read_csv() で読み込んで DataFrameを生成します。

　csvデータには列ラベルがない(ヘッダがない)ので
pd.read_csv(tfile, header=None) のようにheaderオプションを指定。

　行ラベルと列ラベルが付いている行列を検定するケースも載せておきます。

　read_csv() のオプションとして index_col=0 を指定します。

 # chisq04.py (coding: cp932)
 import os
 import pandas as pd
 import prast as ps
 
 csv_data = '''\
 ,Ａ,Ｂ,Ｃ
 Ｗ,17,1,5
 Ｘ,4,6,14
Ｙ,7,13,9
Ｚ,8,9,2
'''
tfile = "tmp.csv"  # 一時的に設けるファイルの名前
with open(tfile, "w") as f:
    f.write(csv_data)
dtf = pd.read_csv(tfile, index_col=0, encoding="cp932")
os.remove(tfile)  # 一時的ファイルを削除
dtf.index = list(dtf.index)
dod = ps.chisq_test(dtf)
for key in dod.keys():
    print(key)
    print(dod[key])
    print('')

　行ラベルと列ラベルはアルファベットですが、半角ではなく全角文字です。

　そのため read_csv() で encoding オプションを指定しています。

行ラベルがあるcsvを読み込むときは index_colオプションを付けて
pd.read_csv(tfile, index_col=0, encoding="cp932") とします。

〜　以上　〜

Copyright (C) T. Yoshiizumi, 2018 All rights reserved.

• 第１回　度数と割合の表示
• 第２回　クロス集計
• 第４回　分布のチェック
• 第５回　t検定
• 第６回　複数回答の集計・前編
• 第７回　複数回答の集計・後編
• python関連記事
• トップページへ