python による統計的検定の実装:t検定、F検定、ウェルチの検定

数学的理論と計算

検定を行う際によく出てくるのが「対応がある」「対応がない」という用語です。

自分も初めは何のことかわからなかったので、まず簡単に説明します。

まず「対応がある」というのは、同じ集団に違った処置を施し、その結果を検定する場合などに使われます。

例えば、ある新薬の効果を確かめる場合、同じ患者の集団に新薬Aと新薬Bを投与した結果に違いを検定する場合です。

逆に「対応がない」場合は、まったく別の集団の特性を検定する場合。

例えばA校の数学のテスト成績とB校のテスト成績の統計量を検定する場合などですね。

というわけで、さっそく対応のない2集団の統計量を検定してみましょう。

Sponsored Link

対応のない2集団の統計量検定

例によって問題は演習で身につける統計学入門の123ページの例題を使っていきましょう。

2種類のチーズAとBについて各種メーカーの製品をそれぞれ36個と34個集めその食塩濃度を測った時標本軽金と不偏標本分散はA:2.04と0.126。Bで1.75と0.113。

AとBの食塩濃度の平均に有意差はあるか有意水準5%で検定。

まず仮説を以下のように立てます。

  • 帰無仮説 チーズAの平均食塩濃度 = チーズBの平均食塩濃度
  • 対立仮説 チーズAの平均食塩濃度 ≠ チーズBの平均食塩濃度

さっそく計算していきましょう。

データは、書籍のデータリンクからダウンロードしたものを読み込みます。

# 必要なライブラリーのインポート
import pandas as pd
import numpy as np

# エクセルファイルのインポート(書籍のデータリンクからダウンロードしたデータ)
c_salt = pd.read_excel('Ex6-1_z_test.xlsx', header = 0)
c_salt.head()  # Aの最初の5行を表示

実行すると、以下のようなDataFrameができます。

次に読み込んだデータに不備がないかを確認していきます。

まずデータの数を確認してみます。

len(c_salt)    # len()関数でデータの長さを確認する。

36         # データが36行あることがわかる。

データ数はAとBで36個あるようです。

データは36あるようですが、データがしっかりそろっているのか分かりませんので、欠損値を確認します。

c_salt.isna().any()    # 欠損値があるかを確認する。
No.    False       # 欠損値無し
A      False       # 欠損値無し
B       True       # 欠損値あり
dtype: bool

結果を見ると、どうやらBに欠損値があるようですので、欠損値の数を確認してみましょう。

c_salt.isna().sum()   # 欠損値の個数を確認する。

No.    0
A      0
B      2        # 欠損値2個
dtype: int64

Bに欠損値が2個あることがわかります。

この欠損値がある行をそのまま削除することもできますが、値を0で埋めておくことにします。

c_salt=c_salt.fillna(0)  # 欠損値を0で埋める。
c_salt.isna().sum()    # 再び欠損値の個数を確認する

No.    0
A      0
B      0
# 欠損値無し

ここから分析しやすいように、DataFrame の整形を行います。

参考サイト:

【Python】欠損値(NaN)を含む行列を抽出・削除・変換・置換する方法|dropnaメソッドの使い方総まとめ
Pythonでcsvやexcelファイルを読み込むときに、空白セルがあると欠損値(NaN)として読み込まれます。 そのような、欠損値(NaN)が含まれる行や列を探したり、削除したり、変換する方法を実例を用いて解説しています。 使用する元デー

DataFrameの整形

ここから分析しやすいようにデータを整形していきます。

まず列Aの値を抜き出します。

c_saltA = c_salt.loc[:, "A"]       # Aの列を抜き出す
c_saltA.head()              # Aの最初の5行を表示

0    2.7
1    2.0
2    1.8
3    1.5
4    2.0
Name: A, dtype: float64

Bも同じように抜き出します。

c_saltB = c_salt.loc[:, "B"]          # Bの列を抜き出す
c_saltB.head()                        # Bの最初の5行を表示

0    2.4
1    1.3
2    2.2
3    1.7
4    1.9
Name: B, dtype: float64

参考サイト:

ここで実際の検定を行うわけですが、検定方法にもいくつかの種類があるわけですね。

とりあえず以下の三つ。

  • Z検定(母分散既知)
  • t検定(母分散未)
  • ウェルチの検定(母分散に違いがある)

カイ二乗検定は質的データ(カテゴリカルデータ)に対しての検定なので、別の記事で解説

今回のケースの場合、不偏標本分散は分かっていますが母分散は分からないので、Z検定以外の検定を行います。

t検定かウェルチ検定かを決める場合、母分散に違いがあるかを確認する必要がありますで、母分散に差があるのかどうか確認するために、F検定を行います。

F検定

サンプルAとBの不偏標本分散は分かっているので、それをもとにF検定を行います。

F値を計算するモジュールはないのでプログラムを書くわけですが、ちょっとまだ書けないので以下の記事を参考にしました。

datadriven-rnd.com

この場合の帰無仮説は「AとBの分散に差はない」であり、対立仮説は「AとBの分散には差がある」になります。

#有意水準を0.05、帰無仮説を「2群間の分散に差がないこと」とする。
A_var = 0.126
B_var = 0.113
A_df = len(c_saltA_n)-1 # A の自由度
B_df = len(c_saltB_n)-1 # B の自由度

#Fを求める場合は、数値が大きい方を分子にする。
bigger_var = np.max([A_var, B_var])
smaller_var = np.min([A_var, B_var])
f = bigger_var/smaller_var
one_sided_pval1 = stats.f.cdf(f, A_df, B_df)
one_sided_pval2 = stats.f.sf(f, A_df, B_df)
two_sided_pval = min(one_sided_pval1, one_sided_pval2)*2

print('F: ', round(f, 3))
print('p-value: ', round(two_sided_pval, 3))
F:    1.115
p-value:  0.749

p>0.05なので、帰無仮説の「AとBの分散に差はない」が採用されます。

分散に差がないことがわかったので、さっそくt検定を行ってみましょう。

ここはt検定のモジュールがあるので、シンプルにそれを使って計算します。

statsmodelsを使ったstudentのt検定

分散に差がないことが分かったので、モジュールstatsmodelsでstudent のt検定を行います。

alternative で両側「two-sided」片側「smaller またはlarger 」を指定できます。

usevar では「pooled」でsutudentのt検定(分散に差がない)「unequal」でウェルチのt検定(分散に差がある)を指定できる。

# 必要なモジュールのインポート
from statsmodels.stats.weightstats import ttest_ind

# 対応無し両側t検定を行う場合
ttest_ind(c_saltA, c_saltB, alternative='two-sided', usevar='pooled')

(3.700279473352106, 0.0004252714975510123, 70.0)

検定の結果は t値が3.700、p値が 0.000425になります。

最後の70は自由度ということなのですが、なぜ70なのかよくわかりません。

この場合、p<0.05なので「AとBの塩分濃度に差がないとは言えない」という帰無仮説が採用されます。

参考サイト:

Pythonのt検定ではstatsmodelsを使う以外に、scipyのstatsを使っても求めることができます。

実際にやってみます。

import scipy.stats as stats 

# 対応無しstudentのt検定を行う場合
stats.ttest_ind(c_saltTA_n, c_saltTB_n)

Ttest_indResult(statistic=3.700279473352104, pvalue=0.00042527149755101505)

ということで同じ結果を得ることができました。

statsmodelsを使ったウェルチの検定

次に2群間で分散に差がある場合の検定「ウェルチの検定」をやってみます。

例によって問題は演習で身につける統計学入門の130p、の例題12から。

ある男子高校で3年A組とB組の生徒からそれぞれ10人と12人を任意に選び、慎重を測定。A組の平均身長はB組の平均身長よりも低いか有意水準5% で検定。

まずはデータをダウンロード。

hight = pd.read_excel('Ex6-5 t test.xlsx', header = 3)
hight.head()

データの確認を行います。

len(hight['A'])             # コラムAのデータの長さを確認。
12

len(hight['B'])             # コラムBのデータの長さを確認。
12

hight.isna().sum()            # 欠損値があるか確認。

No.    0
A      2
B      0
dtype: int64

                     # コラムAに欠損値が2個ある。
hight_comp=hight.fillna(hight.mean()) # 欠損値をコラムAの平均で埋める。
hight_comp.isna().sum()

No.    0
A      0
B      0
dtype: int64

次にDataFrameから、列ごとのデータを抜き出します。

まずはコラムAのデータ。

# コラムAのデータを抜き出す。
hight_compA = hight_comp.loc[:, "A"] 
hight_compA.head()
0    167.0
1    171.0
2    160.0
3    172.0
4    165.0
Name: A, dtype: float64

次にコラムBのデータを抜き出します。

# コラムBのデータを抜き出す。

hight_compB = hight_comp.loc[:, "B"]
hight_compB.head()
0    177
1    175
2    165
3    190
4    171
Name: B, dtype: int64

F検定による分散の差の確認

では、F値を確認して2群間の分散に差があるかどうかを見てみましょう。

#有意水準を0.05、帰無仮説を「2群間の分散に差がないこと」とする。

A_var = np.var(hight_compA, ddof=1)       # Aの分散を代入
B_var = np.var(hight_compB, ddof=1)      # Bの分散を代入
A_df = len(hight_compA)-1            # A の自由度
B_df = len(hight_compB)-1            # B の自由度

#Fを求める場合は、数値が大きい方を分子にする。
bigger_var = np.max([A_var, B_var])
smaller_var = np.min([A_var, B_var])
f = bigger_var/smaller_var
one_sided_pval1 = stats.f.cdf(f, A_df, B_df)
one_sided_pval2 = stats.f.sf(f, A_df, B_df)
two_sided_pval = min(one_sided_pval1, one_sided_pval2)*2

print('F: ', round(f, 3))
print('p-value: ', round(two_sided_pval, 3))

F:    5.704
p-value:  0.008

結果、F値5.704、p値 0.008。

つまりp<0.05なので、帰無仮説は棄却され「2群間の分散に差がないとは言えない」が採用されます。

次にウェルチのt検定で、平均身長に差があるかを見てみましょう。

ウェルチの片側検定

大小関係を検定するので、片側検定を行うことにします。

参考サイト:

片側検定と両側検定の違いや使い分けは?有意水準や棄却域はどう設定? | いちばんやさしい、医療統計
統計学的検定を勉強していくと、「片側検定」と「両側検定」という用語に出会います。 あなたはこの「片側検定と両側検定の違い」を説明できますか? 違いを一言でいうと、「興味のある方向が1つだけかそうじゃないか」ということです。 …おそらく、これ

この場合使用するモジュールは、片側検定対応のstatsmodelsになり、片側右辺を使って検定を行います。

つまり「alternative=’smaller’」さらに、分散に差があるので「usevar=’unequal’」になります。

# 必要なモジュールのインポート

from statsmodels.stats.weightstats import ttest_ind
ttest_ind(hight_compA, hight_compB, alternative='smaller', usevar='unequal')

(-1.0042201353885785, 0.16573595640046906, 14.74210201305821)

結果はt値が-1.004、p値が0.165になります。

つまりp値>0.05なので、帰無仮説「A組の平均身長はB組の平均身長よりも低い」が棄却され、

「A組の平均身長はB組の平均身長よりも低いとは言えない」が採用されます。

Sponsored Link

コメント

タイトルとURLをコピーしました