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

データ分析
Sponsored Link

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

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

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

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

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

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

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

プログラミングを勉強したい!! プログラミングに興味がある!!
そんなあなたへのおすすめ記事
無料体験あり、キャッシュバック有りプログラミングスクールおすすめはこちら

 

対応のない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値を計算するモジュールはないのでプログラムを書くわけですが、ちょっとまだ書けないので以下の記事を参考にしました。

F検定(概要とpython実装)
F検定について対象とする独立な2群の分散が等しい(等分散)であるかを調べる検定。(下記図参照)コード※ scipyモジュールには、 F検定を直接実行するメソッドがないので一つ一つ計算する。用いるデータは「対応のないデータ」と仮定し、有意水準

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

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

A_var = 0.126                    # Aの不偏標準分散
B_var = 0.113                    # Bの不偏標準分散
A_df = len(c_saltA)-1           # A の自由度
B_df = len(c_saltB)-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組の平均身長よりも低いとは言えない」が採用されます。

プログラミングを勉強したい!! プログラミングに興味がある!!
そんなあなたへのおすすめ記事
無料体験あり、キャッシュバック有りプログラミングスクールおすすめはこちら

 

関連記事:

検定の手順とPythonを使った2項分布の検定
検定、つまり統計学的検定は統計の基本として知って知っておいたほうがいい知識です。 じゃ検定は何をするのかというと、例えば「製品AとBの間に優劣があるのかどうか」とか「テストの成績がクラスごとに差があるのかどうか」など、何かと何かを比べ...
統計的検定:Python によるカイ2乗検定の実装
カイ二乗検定について、例題を解きながら学習することにします。 通常のt検定とカイ2乗検定の違いは カイ二乗検定はカテゴリカルデータで独立性を検定 t検定は、連続データで平均の差を検定 になります。 といわれ...

 

Sponsored Link

コメント

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