pythonのpandasによる簡単な統計処理:第3回 F検定,t検定その他

カテゴリー名: [pandasによる簡単な統計処理

2017/11/13

 単純な集計処理だけやっていてもつまらないので、
今回は少し有意性の検証を試みます。

 男女間で身長と体重に有意差があるかをみます。

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

 素材データは、「第1回」 「第2回」と同じ
pt_source.xls です。

 「ID、性別、身長、体重」の4列からなる 400人分のデータです。

    


《このページの目次》


    

はじめに

 pythonをめぐる環境ですが、今回は numpy, scipy.stats もimportします。

    

 python 2.7 において、環境整備でつまずいたので少しメモしておきます。

 WinPython 2.7.10 をインストールした後で、
pandas や numpy をアップグレードしたためかどうか、
t検定のための ttestなどが使えなくなってしまいました(エラー発生)。

 そこで、次の作業をしたら、エラーなく使えるようになりました。

 numexpr のアップグレードは必須ではありませんが、
これをやらないと pythonスクリプトを実行したときに警告メッセージが出ます。

 なお、python 3.5 では上記のような作業をしなくても
トラブルにはなりませんでした。

目次に戻る


1. 有意差を検定する前の準備

 男女間で身長と体重に有意差があるかどうか検証します。

 いわゆる t検定を行いますが、
前提として、まずは男女のグループ分け、
そして、正規分布に即しているかどうかもチェックします。

    

(1) 男女のグループ分け

 まず、データを男女2グループに分類します。

 男女それぞれで人数が違うので、分類した結果を
一つの DataFrame に収納するのはやめます。

 辞書(hhash)に収納することにします。

 また、欠損値を取り除いて収納します。

 以下では、変数dtfに Excelの読み取り結果が既に入っているものとします。
dtf は DataFrame です。

    

◇ groupbyを利用する方法

 groupby() で行うのであれば次のとおり。

dgb = dtf.groupby(u'性別')  # 男女別に分類
dct = dict()  # 空の辞書を用意
for g_name, g_dtf in dgb:
    ser = g_dtf[u'身長']
    dct[g_name] = ser.dropna().reset_index(drop=True)

 上を実行すると、dct[u'女性'] で女性の身長データ(192人分)を参照できます。

 ser.dropna().reset_index(drop=True)
dropna() が欠損値を取り除くためのもの、
reset_index(drop=True) は、行ラベルを付け直すためのものです。

 reset_indexを適用しないと、男女混在の dtf の行ラベルが引き継がれます。
つまり、0, 1, 2, 3, …… というきれいな通し番号にはなりません。

    

◇ 条件式を使う方法

 DataFrame の [……] の中に条件式を書く方法もあります。

male = dtf[dtf[u'性別'] == u'男性']

 上のようにすると、変数 male に男性だけのデータが代入されます。

 male は、「ID、性別、身長、体重」の4列からなる DataFrame です。

 このことを利用した分類スクリプトは次のとおり。

dct = dict()  # 空の辞書を用意
for key in dtf[u'性別'].dropna().unique():
    ser = dtf[dtf[u'性別'] == key][u'身長']
    dct[key] = ser.dropna().reset_index(drop=True)

 dtf[u'性別'].dropna().unique() というのは、今回のケースでは
[u'女性', u'男性'] と実質的に同じです。

目次に戻る


(2) データが正規分布に即しているかどうかの検証

 2グループの平均値の差を検証する t検定は、
標本が正規分布に即していることを前提にしているといわれます。

 標本サイズが 30以上(各グループ 15以上)のときは、
必ずしも正規分布を前提にしなくてもいいという話もあるようですが……

 それはともかく、ここでは正規分布に即しているかどうかを
shapiro() で確認します。

 このメソッドが返す p値(有意確率の値)が 0.05 より小さければ
「正規分布である」という仮説を棄却してもいいようです。
つまり、「正規分布でない」といって差し支えないことになります。

 男性の身長のデータが変数 ser に入っているとき、

wval, pval = stats.shapiro(ser)

 上のようにすると、pval に p値がセットされています。

    

 これ以降に具体的なスクリプトを掲げます。

 「身長、体重」×「男性、女性」の計4つ series を用意して、
それぞれを検証します。

 4つを別々の変数に格納したのでは扱いにくいので、
辞書(hash)を入れ子にして用います。

 dct[u'身長'][u'男性'] に男性の身長を代入するかたちです。

 結果は、Excelファイルとして書き出します。

 先に結果を掲げると次のとおり。

  標本サイズ p.value 正規性
男性・身長 198 0.755964816 Yes
女性・身長 192 0.617443383 Yes
男性・体重 201 0.028328262 No
女性・体重 190 0.002008719 No

 身長は男女とも正規分布に即しているようですが、
体重の方は、正規分布とは認めにくいようです。

 標本サイズは、今回の例では要するに「人数」のことです。

    

 以下、pd01.py 全体を示します。

 1# pd01.py (coding: cp932)
 2import os, sys
 3import pandas as pd
 4import scipy.stats as stats
 5import xlwt
 6
 7from platform import python_version
 8if int(python_version()[0]) < 3:  # python ver 2 の場合
 9    reload(sys)
10    sys.setdefaultencoding('cp932')  # デフォルト文字コードを変更
11
12xls_file = "pt_source.xls"
13dtf = pd.read_excel(xls_file)
14dgb = dtf.groupby(u'性別')  # 男女別に分類
15
16dct = {u'身長': dict(), u'体重': dict()}
17for key in dct.keys():  # 「身長、体重」それぞれを処理
18    for g_name, g_dtf in dgb:  # 「男性、女性」に切り分け
19        ser = g_dtf[key]
20        dct[key][g_name] = ser.dropna().reset_index(drop=True)
21
22ilist = list()  # 行ラベルをこれに記録
23dtf2 = pd.DataFrame()  # 空の DataFrame
24for key1 in dct.keys():  # 身長、体重をたどる
25    for key2 in dct[key1].keys():  # 男性、女性をたどる
26        idx = u"%s・%s" % (key2, key1)  # 行ラベル
27        ilist.append(idx)
28        size = dct[key1][key2].count()  # 標本サイズ
29        wval, pval = stats.shapiro(dct[key1][key2])  # 正規性の検定
30        if pval < 0.05:
31            reg = 'No'
32        else:
33            reg = 'Yes'
34        ser = pd.Series([size, pval, reg])
35        dtf2 = dtf2.append(ser, ignore_index=True)
36dtf2.columns = [u'標本サイズ', 'p.value', u'正規性']
37dtf2.index = ilist
38dtf2.to_excel("pd01.xls", u'shapiroによる正規性の検証')

 python には sprintf() がないようですが、
u"%s・%s" % (key2, key1) のような書き方で同じことができるようです。

目次に戻る


2. 等分散かどうかの検証(F検定)

 t検定に入る前に、男女の身長の分散(ばらつき)が等しいかどうか確認します。

 正規性が認められる「身長」の方を確認します。

 等分散か否かで t検定のやり方が違ってくるためです。

(1) pythonにおけるF検定

 分散が等しいかどうかなどを検証するのは F検定です。

 pythonでは F検定を次ぎの手順で行うようです。

 男性の身長データが変数 male、女性のデータは female に入っているとします。

male_var = np.var(male, ddof=1)  # 男性の不偏分散
female_var = np.var(female, ddof=1)  # 女性の不偏分散
male_df = len(male) - 1  # 男性の自由度
female_df = len(female) - 1  # 女性の自由度
f = male_var / female_var  # F比の値
pval1 = stats.f.cdf(f, male_df, female_df)  # 片側検定のp値・その1
pval2 = stats.f.sf(f, male_df, female_df)  # 片側検定のp値・その2
pval3 = min(pval1, pval2) * 2  # 両側検定のp値

 今回、F比の値は 0.95 という 1 に近い値です(おそらく等分散)。

 また、両側検定の pval3 が 0.7239 という高い値なので、
男女が等分散とみていいとおもいます。

 ちなみに、片側検定の2番目のp値については
pval2 = 1.0 - pval1 という式で求めることもできます。

目次に戻る


(2) F検定における両側検定と片側検定の確認

 今回は、「男性の分散(38.63) < 女性の分散(40.64)」の関係です。

 なので、f = male_var / female_var の値が f < 1.0 になります。

 この場合、f.cdf() で求めた pval1 が 0.362
 f.sf() で求めた pval2 が 0.638 となって
pval1 < pval2 の関係になります。

 一方、F比を算出するときに分子と分母を入れ替えて f > 1.0 にすれば
pval1, pval2 の値が入れ替わって、
pval1 > pval2 の関係になります。

 今回の3種類のp値の意味を逆説的にいうと次のとおり。

 今回は、結果的にどの対立仮説も採択できません。

 つまり帰無仮説の方を採択し、等分散と見てよいことになります。

 片側検定の pval1, pval2 のどちらか一方だけを参照して
「等分散か否か」を判断するのは根拠が弱いような気がしますが、
python関係のWebを見ると、一方だけ(cdfだけ)を参照しているケースが
結構あるような気がします。

 「小さいか否か」 「大きいか否か」ではなく
「等しいか否か」をチェックするなら両側検定だと思いますが、
そうでもないんでしょうか?

 片側検定の p_value が分かれば、両側検定の値も分かるはず、
ということかもしれませんが……

目次に戻る


(3) F検定を関数化

 F検定を関数化してみたので掲載します。

var_test(a, b, alternative="two-sided")

 検定の結果は、構造体 FtestRes にセットする形にしました。

 統計解析ソフトRの関数をまねっこして、
alternativeオプションを指定できるようにしました。

 両側検定のときは "two-sided" (デフォルト)
片側検定なら "less" または "greater" です。
対立仮説を指示するものです。

    

◇ スクリプトを実行したときの出力

 まず、実行結果がどうなるかを掲げます。

 男女の身長の分散について print() で出力。

f = 0.950586713521
p_value = 0.723947458808
alt = 対立仮説: ほんとのF比が 1.0 と等しくない.
var = [38.632291698713011, 40.640470931500872]
df = [197, 191]

 対立仮説の「ほんとのF比」は、母集団のF比というニュアンスのつもり。

 左に表示されている f, p_value, alt, var, df は、構造体のメンバー名です。

    

◇ 構造体や関数の定義

 以下、zipファイルにはいっている pd02.py の抜粋です。

    (前略)
# 構造体の定義(F検定の結果をこれに記録)
from collections import namedtuple
FtestRes = namedtuple("FtestRes", "f p_value alt var df")

# F検定: alternative="two-sided", "less", "greater"
def var_test(a, b, alternative="two-sided"):
    alt_arg = alternative
    aa = np.array(a)
    bb = np.array(b)
    aa_var = np.var(aa, ddof=1)  # aの不偏分散
    bb_var = np.var(bb, ddof=1)  # bの不偏分散
    aa_df = len(aa) - 1  # aの自由度
    bb_df = len(bb) - 1  # bの自由度
    f = aa_var / bb_var  # F比の値
    pval1 = stats.f.cdf(f, aa_df, bb_df)  # 片側検定のp値・その1
    pval2 = stats.f.sf(f, aa_df, bb_df)  # 片側検定のp値・その2
    pval3 = min(pval1, pval2) * 2  # 両側検定のp値
    if alt_arg == 'less':
        alt = u"対立仮説: ほんとのF比が 1.0 より小さい."
        pval = pval1
    elif alt_arg == 'greater':
        alt = u"対立仮説: ほんとのF比が 1.0 より大きい."
        pval = pval2
    elif alt_arg == 'two-sided':
        alt = u"対立仮説: ほんとのF比が 1.0 と等しくない."
        pval = pval3
    else:
        alt = "3 p_values(less, greater, two-sided)"
        pval = [pval1, pval2, pval3]
    fr = FtestRes(f, pval, alt, [aa_var, bb_var], [aa_df, bb_df])
    return fr

    (中略)
fr = var_test(male, female, alternative="two-sided")
for key in fr._fields:
    print("%s = " % key + str(fr._asdict()[key]))

    

 関数 var_test() のオプション alternative は、
デフォルトが "two-sided" なので、それでよければ省略可能。

 alternativeオプションを two-sided, less, greater 以外にした場合、
たとえば all にしたとすると、三つのp_value を得ることができます。

 fr.p_value が、三つの数値が格納されたリストになります。

    

 構造体 fr については、fr.f, fr.p_value のようにして
個々の検定結果を参照できます。

 上のスクリプトでは、すべての要素を出力するようにしています。

 python2 では dct = vars(fr) のようにすると
実質的に構造体を辞書に変換できるようですが、
python3 でそれをやるとエラーになるので
_fields_asdict() をつかいました。

    

 なお、95%信頼区間の下限・上限を求めたいときは
変数frに構造体が代入された後で下のようにします。

i_l, i_u = stats.f.interval(0.95, fr.df[0], fr.df[1])
ci_val = (fr.var[0] / (fr.var[1] * i_u), fr.var[0] / (fr.var[1] * i_l))

 ci_val[0] が下限値、
 ci_val[1] が上限値です。

目次に戻る


3. 平均値の差の検証(t検定)

 身長は、正規分布に即しているようなので t検定にかけてみます。

 男女間に有意差が認められるかどうかの検証です。

 身長について男女が等分散なので、その前提で検定します。

(1) t検定(両側検定)

 scipy.stats には ttest_ind() というメソッドがあります。

 今回の男女比較では、2グループの人数が不揃いなど
対応関係にあるとはいえないので、この ttest_ind をつかいます。

 (対応関係があるときは ttest_rel を使うようです。)

 ttest_ind は、統計量のt値、有意確率のp値の二つを返します。

 男性の身長データが変数 male に、女性のデータが female に入っている場合、

t, p = stats.ttest_ind(male, female, equal_var=True)

 上のように用います。

 equal_var は、等分散なら True、不等分散なら False にします。

 男性の身長の平均値が女性のそれより大きな値なので、
male, female の順番で引数を与えると t値がプラスになり(10.378)、
female, male にすると t値がマイナスになります(-10.378)。

 このことは、p値の解釈の際にかかわってきます。

    

 ttest_ind が返す p値は、両側検定の値です。

 p値が 0.05 より小さいなら、男女の身長の平均値は等しくないと解釈できます。

 今回は p = 1.98017665846e-22 という非常に小さい値なので
「男性が女性と異なる」と判断していいことになります。

目次に戻る


(2) pythonにおけるt検定の片側検定

 「小さいか否か」 「大きいか否か」を検定したいなら片側検定ですが、
それを見るには他の p値が必要です。

 ttest_ind に片側検定を行うオプションがなさそうなので
便法として次ぎのようにします。

t, p = stats.ttest_ind(male, female, equal_var=True)
pval3 = p
pval2 = p / 2.0
pval1 = 1.0 - pval2
if t < 0.0:
    w = pval2
    pval2 = pval1
    pval1 = w

 上を実行して得られる pval3 は、両側検定の p値です。

 pval1 は、「第1引数の値が第2の値より小さい」(対立仮説)に関する値、
 pval2 は、「第1引数の値が第2の値より大きい」(対立仮説)に関する値です。

 pval1 < 0.05 なら、「小さい」の対立仮説を採択できます。

 pval2 < 0.05 なら、「大きい」の対立仮説を採択できます。

 今回のケースでいうと、male, female の順番にした場合、

pval1 = 1.0
pval2 = 9.90088329228e-23
pval3 = 1.98017665846e-22

 上のような結果になるので、次の判断が可能です。

 「男性が女性より小さい」:採択できない。
 「男性が女性より大きい」:採択できる。
 「男性が女性と異なる」:採択できる。

 pval1 は、厳密には 1.0 でないはずですが、
通常の計算式で処理すると 1.0 になります。

 高い精度を保ったまま値を計算する方法があるんでしょうか?

目次に戻る


(3) t検定を関数化

 t検定を関数化してみたので掲載します。

t_test(a, b, equal_var=True, alternative="two-sided")

 検定の結果は、構造体 TtestRes にセットする形にしました。

 equal_var は、等分散なら True、そうでなければ False を指定。

 F検定の関数と同様、alternativeオプションを指定できます。

    

◇ スクリプトを実行したときの出力

 まず、実行結果がどうなるかを掲げます。

 男女の身長の平均値について print() で出力。

t = 10.3775809653
p_value = 1.98017665846e-22
alt = 対立仮説: ほんとの平均値の差が 0.0 と等しくない.
mean = [166.71565656565656, 160.09947916666667]
df = 388

 左に表示されている t, p_value, alt, mean, df は、構造体のメンバー名です。

    

◇ 構造体や関数の定義

 以下、zipファイルにはいっている pd03.py の抜粋です。

    (前略)
# 構造体の定義(t検定の結果をこれに記録)
from collections import namedtuple
TtestRes = namedtuple("TtestRes", "t p_value alt mean df")

# t検定: alternative="two-sided", "less", "greater"
def t_test(a, b, equal_var=True, alternative="two-sided"):
    eqv = equal_var
    alt_arg = alternative
    aa = np.array(a)
    bb = np.array(b)
    aa_mean = np.mean(aa)  # aの平均値
    bb_mean = np.mean(bb)  # bの平均値
    aa_df = len(aa) - 1  # aの自由度
    bb_df = len(bb) - 1  # bの自由度
    t, p = stats.ttest_ind(aa, bb, equal_var=eqv)
    pval3 = p
    pval2 = p / 2.0
    pval1 = 1.0 - pval2
    if t < 0.0:
        w = pval2
        pval2 = pval1
        pval1 = w
    if alt_arg == 'less':
        alt = u"対立仮説: ほんとの平均値の差が 0.0 より小さい."
        pval = pval1
    elif alt_arg == 'greater':
        alt = u"対立仮説: ほんとの平均値の差が 0.0 より大きい."
        pval = pval2
    elif alt_arg == 'two-sided':
        alt = u"対立仮説: ほんとの平均値の差が 0.0 と等しくない."
        pval = pval3
    else:
        alt = "3 p_values(less, greater, two-sided)"
        pval = [pval1, pval2, pval3]
    tr = TtestRes(t, pval, alt, [aa_mean, bb_mean], aa_df+bb_df)
    return tr

    (中略)
tr = t_test(male, female, equal_var=True, alternative="two-sided")
for key in tr._fields:
    print("%s = " % key + str(tr._asdict()[key]))

 equal_var の値は、そのまま ttest_ind に引き渡されます。

 equal_var, alternative は、デフォルト値でよければ省略可能です。

 t検定したときの信頼区間を求めるのにチャレンジしてみましたが、
どうも自信ないのでパスします。

 stats.t.ppf() を使ってごちゃごちゃやると
ごく近い値にはなりますが、いまいちなので……

目次に戻る


4. 正規分布を前提にしない2群の差の検証(U検定)

 体重は、正規分布に即していないようなので、
マン・ホイットニーのU検定を用いてみます。

 男女間の体重に有意差が認められるかどうかの検証です。

 今回のように標本サイズがそれなりにあり、
極端な値に平均値が左右されている訳でもない場合、
t検定(不等分散)を使っていいとおもいますが、
参考までマン・ホイットニーのU検定を試します。

    

◇ マン・ホイットニーのU検定を行うメソッド

 scipy.stats には mannwhitneyu() があるので、それを使います。

 男女それぞれのデータが変数 male, female に代入されている場合

uval, pval = stats.mannwhitneyu(male, female, alternative="two-sided")

 上のように用います。

 結果は、uval に統計量Uの値、
pval には p_value が入ります。

    

 このメソッドでは alternativeオプションを指定できます。

 less, greater, two-sided の三つのどれかを指定しますが、
デフォルトは None になっていて、結果が三つとは違ってきます。

 None の p_value は、two-sided の p_value の半分の値で、
統計量Uは、三つのどれとも違います。

 None がデフォルトになっているのは、旧バージョンとの互換性のためでは?

 実際には三つのうちのどれかを指定するのがいいとおもいます。

    

◇ スクリプトの実行結果

 スクリプト pd04.py の実行結果は下のとおり(print の出力)。

 alternative を less, greater, two-sided の三つとも試しています。

対立仮説: less
U=30267.0, p_value=1.0

対立仮説: greater
U=30267.0, p_value=7.46205078273e-24
「男性の体重が女性の体重より大きい」を採択可能.

対立仮説: two-sided
U=30267.0, p_value=1.49241015655e-23
「男性の体重が女性の体重と等しくない」を採択可能.

    

◇ スクリプト pd04.py の全体像

 以下、体重の検証を行う pd04.py です。

 1# pd04.py (coding: cp932)
 2import os, sys
 3import pandas as pd
 4import numpy as np
 5import scipy.stats as stats
 6
 7from platform import python_version
 8if int(python_version()[0]) < 3:  # python ver 2 の場合
 9    reload(sys)
10    sys.setdefaultencoding('cp932')  # デフォルト文字コードを変更
11
12xls_file = "pt_source.xls"
13dtf = pd.read_excel(xls_file)
14dgb = dtf.groupby(u'性別')  # 男女別に分類
15
16dct = {u'身長': dict(), u'体重': dict()}
17for key in dct.keys():  # 「身長、体重」それぞれを処理
18    for g_name, g_dtf in dgb:  # 「男性、女性」に切り分け
19        ser = g_dtf[key]
20        dct[key][g_name] = ser.dropna().reset_index(drop=True)
21
22male = dct[u'体重'][u'男性']
23female = dct[u'体重'][u'女性']
24alt_dct = {'less': u'より小さい', 'greater': u'より大きい',
25  'two-sided': u'と等しくない'}
26for alt in ['less', 'greater', 'two-sided']:
27    uval, pval = stats.mannwhitneyu(male, female, alternative=alt)
28    print(u'対立仮説: ' + alt)
29    print("U=%s, p_value=%s" % (str(uval), str(pval)))
30    if pval < 0.05:
31        print(u"「男性の体重が女性の体重%s」を採択可能." % alt_dct[alt])
32    print('')

    

◇ 検定の名称について

 マン・ホイットニーのU検定は、
ウィルコクソンの順位和検定といわれるものと
実質的に同じということのようです。

 これらは、対応関係のない2群の違いを検証します。

 一方、対応関係がある場合に用いるのは「ウィルコクソンの符号順位検定」です。

 この「ウィルコクソンの符号順位検定」をpythonで行うときは
stats.wilcoxon() を使います。

 ウィルコクソンという名前がよく出てくるので少々混乱しますが、
python の wilcoxon() は、「符号順位検定」です。

    

 今回は これで終了です。

 身長と体重の相関のチェックとか、
カテゴリ化したクロス集計表へのカイ2乗検定の適用
といったのも試したかったのですが、別の機会にします。

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