群間の平均値比較

5月から職種がデータアナリストになります。
検証する仮説として、平均値等しい、違うのでは系がよくありそうだなと感じたため、Pythonで実装練習をしてみました。

対応なしデータ(iris)にて、以下①②③手順で行いました。
説明変数は、「sepal length, sepal width, petal length, petal width」の4つあります。今回平均値比較に使用したのは、petal widthです。
petal widthにしたのは、ヒストグラムみて、ああ、帰無仮説を棄却できるなあと思ったからです。
①正規性の確認
②等分散性の確認
③平均値比較

今回使ったirisの散布図行列/target変数ごとのヒストグラムは以下です。

## データロード
from sklearn.datasets import load_iris
iris = load_iris()

import pandas as pd
df = pd.DataFrame(iris.data, columns=iris.feature_names)
df['target'] = iris.target
df.loc[df['target'] == 0, 'target'] = "setosa"
df.loc[df['target'] == 1, 'target'] = "versicolor"
df.loc[df['target'] == 2, 'target'] = "virginica"

import seaborn as sns
sns.pairplot(df, hue = "target")

①正規性の確認
(i)Q-Qプロットを使って平均値比較　
※直線上に点が乗ってるほど、正規分布に当てはまりよい、と解釈します。
上から、setosa, versicolor, virsinica のプロットです。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

nsample = 288
np.random.seed(1)


#左寄りである点にデータ集中してそう。
stats.probplot(df_setosa_PW, plot=plt)
plt.show()

stats.probplot(df_versicolor_PW, plot=plt)
plt.show()


#stats.probplot(x2, plot=plt)

#setosaよりは、正規性たかそう
stats.probplot(df_virginica_PW, plot=plt)
plt.show()

見た感じ、virsinica以外の2つは、正規分布ぽくなさそうですね。
例えばsetosaのpetal widthは、0.2付近にめっちゃあつまってます。

(ii) シャピロウィルク検定
これもQQプロット同様、正規性を満たすかどうかを判断するものです。
QQプロットはヒストグラムのように、視覚的に判断するのに使いますが、
こちらは、p-valueが出力されるため、0.05未満→正規分布じゃなさそう。0.05以上→正規分布といってもよさそう。のように、使えます。

from scipy import stats
import numpy as np
from matplotlib import pyplot as plt


#シャピロ-ウィルク検定実施
result_setosa_PW = stats.shapiro(df_setosa_PW)
result_versicolor_PW = stats.shapiro(df_versicolor_PW)
result_virginica_PW = stats.shapiro(df_virginica_PW)

print('p value(setosa):', result_setosa_PW.pvalue)
print('p value(versicolor):', result_versicolor_PW.pvalue)
print('p value(virginica):', result_virginica_PW.pvalue)

以下、出力です。virginicaのみが、p value 0.05以上です。つまり、virginicaのみ正規分布といえる。残り2つは正規分布といえない、という、QQプロットと同じ結論にいたります。

p value(setosa): 8.65842082475865e-07
p value(versicolor): 0.027278218418359756
p value(virginica): 0.08695744723081589

②等分散性の確認
次は、等分散性の確認をします。setosaとversicolorの2つで行います。
これは、2群間の分散が等しいかどうかを検証するためのものです。
この結果により、平均値比較を行うt検定の種類が変わります。


本来これは、正規性が確認できたうえで使うものです。
よって、本当は①より得られた結論で使うべきではありませんが、実装練習のため、無理やりつかってみました。

#不偏分散を計算

var_setosa_PW = np.var(df_setosa_PW)
var_versicolor_PW = np.var(df_versicolor_PW)
#var_virginica_PW = np.var(df_virginica_PW)

#F分布の自由度
free_setosa_PW = np.size(df_setosa_PW) - 1
free_versicolor_PW = np.size(df_versicolor_PW) - 1
#free_virginica_PW = np.size(df_setosa_PW) - 1

s1 = np.max([var_setosa_PW ,var_versicolor_PW])
s2 = np.min([var_setosa_PW ,var_versicolor_PW])

#F値:フィッシャーの分散比
F = s1 / s2

#F検定
## F(α/2)[m-1, n-1]
pval_left = stats.f.cdf(F, free_setosa_PW, free_versicolor_PW)
pval_right = stats.f.sf(F, free_setosa_PW, free_versicolor_PW)

##両側検定のP値
pval_both = min(F_left, F_right) * 2  # 両側検定のp値

print("p value:",pval_both)


if pval_both < 0.05:
   print("等分散でない")

else:
   print("等分散である")

p値が0.05未満となり、等分散でないという結論にいたりました。

p value: 2.1287944125767378e-05
等分散でない

③平均値比較

・比較したい群数(target変数のラベル種類の数)が2か、3以上か
・①②の結果で異なります。

今回は、①にて「正規性があるとはいえない」の結果が得られたことと、
setosa, versicolor, virsinicaの3群比較をしたいことから、

クラスカルウォリス検定を実施します。

result = stats.kruskal(df_setosa_PW, df_versicolor_PW, df_versicolor_PW)
print("result:",result)

if result.pvalue < 0.05:
   print("H0棄却:平均は異なる")
else: print("H0採択:平均は同じ")

結論、平均値は異なる　※3群中、どこかの平均値が異なる
の結論に至りました。

result: KruskalResult(statistic=101.06153041974676, pvalue=1.134403356921334e-22)
H0棄却:平均は異なる

↑のクラスカルウォリスでは、どの群間に差があったかは不明です。
多重比較を使うと、どの群に差があるかを知ることができます。

まだ触れたことがないため、お勉強いたします。。

ここまで読んでくださり、まことにありがとうございました！！(^▽^)/