カテゴリー名： ［pandasによる簡単な統計処理］
単純な集計処理だけやっていてもつまらないので、
今回は少し有意性の検証を試みます。
男女間で身長と体重に有意差があるかをみます。
当Webページで紹介するスクリプトや素材データ一式は、
pandas03.zip という圧縮ファイルに同梱しておきます。
素材データは、「第１回」 「第２回」と同じ
pt_source.xls です。
「ＩＤ、性別、身長、体重」の４列からなる 400人分のデータです。
pythonをめぐる環境ですが、今回は numpy, scipy.stats もimportします。
python 2.7 において、環境整備でつまずいたので少しメモしておきます。
WinPython 2.7.10 をインストールした後で、
pandas や numpy をアップグレードしたためかどうか、
t検定のための ttestなどが使えなくなってしまいました(エラー発生)。
そこで、次の作業をしたら、エラーなく使えるようになりました。
python -m pip install "numpy-1.13.3+mkl-cp27-cp27m-win32.whl" --upgradepython -m pip install numexpr --upgrade numexpr のアップグレードは必須ではありませんが、
これをやらないと pythonスクリプトを実行したときに警告メッセージが出ます。
なお、python 3.5 では上記のような作業をしなくても
トラブルにはなりませんでした。
男女間で身長と体重に有意差があるかどうか検証します。
いわゆる t検定を行いますが、
前提として、まずは男女のグループ分け、
そして、正規分布に即しているかどうかもチェックします。
まず、データを男女２グループに分類します。
男女それぞれで人数が違うので、分類した結果を
一つの DataFrame に収納するのはやめます。
辞書(hhash)に収納することにします。
また、欠損値を取り除いて収納します。
以下では、変数dtfに Excelの読み取り結果が既に入っているものとします。
dtf は DataFrame です。
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 は、「ＩＤ、性別、身長、体重」の４列からなる 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'男性'] と実質的に同じです。
２グループの平均値の差を検証する t検定は、
標本が正規分布に即していることを前提にしているといわれます。
標本サイズが 30以上(各グループ 15以上)のときは、
必ずしも正規分布を前提にしなくてもいいという話もあるようですが……
それはともかく、ここでは正規分布に即しているかどうかを
shapiro() で確認します。
このメソッドが返す p値(有意確率の値)が 0.05 より小さければ
「正規分布である」という仮説を棄却してもいいようです。
つまり、「正規分布でない」といって差し支えないことになります。
男性の身長のデータが変数 ser に入っているとき、
wval, pval = stats.shapiro(ser)
上のようにすると、pval に p値がセットされています。
これ以降に具体的なスクリプトを掲げます。
「身長、体重」×「男性、女性」の計４つ series を用意して、
それぞれを検証します。
４つを別々の変数に格納したのでは扱いにくいので、
辞書(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) のような書き方で同じことができるようです。
t検定に入る前に、男女の身長の分散(ばらつき)が等しいかどうか確認します。
正規性が認められる「身長」の方を確認します。
等分散か否かで t検定のやり方が違ってくるためです。
分散が等しいかどうかなどを検証するのは 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値・その１
pval2 = stats.f.sf(f, male_df, female_df) # 片側検定のp値・その２
pval3 = min(pval1, pval2) * 2 # 両側検定のp値
今回、F比の値は 0.95 という 1 に近い値です(おそらく等分散)。
また、両側検定の pval3 が 0.7239 という高い値なので、
男女が等分散とみていいとおもいます。
ちなみに、片側検定の２番目のp値については
pval2 = 1.0 - pval1 という式で求めることもできます。
今回は、「男性の分散(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 の関係になります。
今回の３種類のp値の意味を逆説的にいうと次のとおり。
pval1(実際は 0.362)<0.05 なら「男性が女性より小さい」(対立仮説)と見てよい。pval2(実際は 0.638)<0.05 なら「男性が女性より大きい(対立仮説)と見てよい。pval3(実際は 0.724)<0.05 なら「男性が女性と異なる」(対立仮説)と見てよい。今回は、結果的にどの対立仮説も採択できません。
つまり帰無仮説の方を採択し、等分散と見てよいことになります。
片側検定の pval1, pval2 のどちらか一方だけを参照して
「等分散か否か」を判断するのは根拠が弱いような気がしますが、
python関係のWebを見ると、一方だけ(cdfだけ)を参照しているケースが
結構あるような気がします。
「小さいか否か」 「大きいか否か」ではなく
「等しいか否か」をチェックするなら両側検定だと思いますが、
そうでもないんでしょうか？
片側検定の p_value が分かれば、両側検定の値も分かるはず、
ということかもしれませんが……
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値・その１
pval2 = stats.f.sf(f, aa_df, bb_df) # 片側検定のp値・その２
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] が上限値です。
身長は、正規分布に即しているようなので t検定にかけてみます。
男女間に有意差が認められるかどうかの検証です。
身長について男女が等分散なので、その前提で検定します。
scipy.stats には ttest_ind() というメソッドがあります。
今回の男女比較では、２グループの人数が不揃いなど
対応関係にあるとはいえないので、この 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 という非常に小さい値なので
「男性が女性と異なる」と判断していいことになります。
「小さいか否か」 「大きいか否か」を検定したいなら片側検定ですが、
それを見るには他の 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 は、「第１引数の値が第２の値より小さい」(対立仮説)に関する値、
pval2 は、「第１引数の値が第２の値より大きい」(対立仮説)に関する値です。
pval1 < 0.05 なら、「小さい」の対立仮説を採択できます。
pval2 < 0.05 なら、「大きい」の対立仮説を採択できます。
今回のケースでいうと、male, female の順番にした場合、
pval1 = 1.0
pval2 = 9.90088329228e-23
pval3 = 1.98017665846e-22
上のような結果になるので、次の判断が可能です。
「男性が女性より小さい」：採択できない。
「男性が女性より大きい」：採択できる。
「男性が女性と異なる」：採択できる。
pval1 は、厳密には 1.0 でないはずですが、
通常の計算式で処理すると 1.0 になります。
高い精度を保ったまま値を計算する方法があるんでしょうか？
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() を使ってごちゃごちゃやると
ごく近い値にはなりますが、いまいちなので……
体重は、正規分布に即していないようなので、
マン・ホイットニーのU検定を用いてみます。
男女間の体重に有意差が認められるかどうかの検証です。
今回のように標本サイズがそれなりにあり、
極端な値に平均値が左右されている訳でもない場合、
t検定(不等分散)を使っていいとおもいますが、
参考までマン・ホイットニーの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 です。
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検定は、
ウィルコクソンの順位和検定といわれるものと
実質的に同じということのようです。
これらは、対応関係のない２群の違いを検証します。
一方、対応関係がある場合に用いるのは「ウィルコクソンの符号順位検定」です。
この「ウィルコクソンの符号順位検定」をpythonで行うときは
stats.wilcoxon() を使います。
ウィルコクソンという名前がよく出てくるので少々混乱しますが、
python の wilcoxon() は、「符号順位検定」です。
今回は これで終了です。
身長と体重の相関のチェックとか、
カテゴリ化したクロス集計表へのカイ２乗検定の適用
といったのも試したかったのですが、別の機会にします。
Copyright (C) T. Yoshiizumi, 2017 All rights reserved.