ゼロから始めるプログラミングの学習日記

未経験の状態からプログラミング初めます

課題：pythonで（簡単な）仮説検定をしたい

パッケージのインストールはpip3 installで

pythonでの仮説検定も基本的にはRと同じでpackage頼りです。
肝になる部分は「どのpackageをどこで使うの？」ということだと思います。
簡単な仮説検定まで一気に紹介していきましょう。
スクリプトの詳しい内容や検定の説明については、統計の参考書や他のサイトに任せて、
「なんだか分からないけどできるわ」を目指しました。

検定する前の準備

「よし、検定しよう！」と思っても何かしらデータがなければ意味がありません。
ここで、データの整理の仕方を確認していきましょう。

必要そうなpackageはimportしていく

Rと違う点は、まずpythonファイル上でpackageをimportする指示を出すところです。
installさえされていれば、import ほにゃららと書くだけで非常に簡単です。

import pandas as pd  
# dataframeを作るのに必要  
import numpy as np  
# 平均・中央値・分散・偏差などが計算できる
import statistics  
# importする必要があるのかはよく分からなかったけど最頻値を計算するのに使った
import matplotlib.pyplot as plt  
# ヒストグラムなど図表を描く
from scipy import stats  
# ほとんどの仮説検定するのに使った
from statsmodels.stats.multicomp import pairwise_tukeyhsd
# チューキークレイマー検定をするのに使った
import scikit_posthocs as sp 
# スティール・ドゥワス検定をするのに使いたかったがうまくrunできず
お作法がありそうです

例えば、

import pandas

とすれば、pandasがimportされ指示が出せますが、pdとしてimportするのが普通みたいです。
numpy (ナンピーじゃなくてナンパイなんだってさ)の場合もnpとしてimportします。
上述の通りですが、名前を付けてimportする場合は、as ほにゃほにゃとしてimportしておきます。 正直、この辺はよく分かっていないので、周囲の人に聞くか仕様書に合わせて書くのが無難そうですね。
（他の人が見たときに「なんじゃこりゃ」となると非常に迷惑をかけるので）

csvファイルの読み込み

pythonでcsvファイルを読み込むとき、２行目の列だけ読み込みたい場合は、

pd.read_csv('practice.csv',usecols=[1])
# ここでは、practiceという名前のcsvファイルの２行目の列を読み込んでいます。
# 2行目なのになんで[1]なの？と思うかもしれませんが、これはpythonが0始まりでカウントするからです。

ちなみに…毎回自分もこんがらがるのですが「行が横（水平方向）」で「列が縦（垂直方向）」らしいです。
英語では「行がrow」で「列がcolumn」です。
じゃぁ、userowsで指定した行が読み込まれるのか？というとそういう訳ではないようです。

データフレームと次元の平坦化

前述のpractice.csvというファイルを読み込む時にこれにAという名前をつけるとします。

A = pd.read_csv('practice.csv',usecols=[1])
# Aという名前にpractice.csvの2行目の列が格納されました。

このうちの1-30行のデータをA1to30として読み込みたい場合、numpyを使います。
numpyはnpとしてimportしているので、

A1to30 = np.array(A[0:29])
# arrayは配列という意味の英語です。
# くどいようですが、0始まりなので0〜29で1から30行を抜き出していることになります。

次に抜き出したA1to30という配列からA1to30dfというデータフレームを作るとします。
A1to30のDataFrameなのでA1to30dfです。

A1to30df = pd.DataFrame(data=A1to30)  
# Dataframeを作るには、pandasを使います。

説明しやすいかな？と思ってAとかA1to30とか直感的な名前を使ってしまったので少し冗長な説明になってしまいましたね。

配列の次元？

配列には次元という考え方があります。
例えば、 AさんとBさんとCさんの走った距離について考えましょう。
「Aさんは3km走った」「Bさんは1.5km」「Cさんは6km走った」という場合、1行（行は横でした）だけで管理できそうですね。
これが１次元データです。
ここにAさん、Bさん、Cさんが走るのにかかった時間をデータとして足してみましょう。
「Aさんは3kmを30分で走った」「Bさんは1.5kmを10分で走った」「Cさんは6kmを60分で走った」とすると、
今度は、2行で管理することになります。これが２次元データです。
さらにAさん、Bさん、Cさんが走りながら飲んだ水の量をデータとして足すとどうなるでしょうか。
「Aさんは3kmを30分で走りながら100ml飲んだ」「Bさんは1.5kmを10分で走りながら50ml飲んだ」「Cさんは6kmを60分で走りながら150ml飲んだ」
となると、今度は３行のデータになります。これが３次元データです。
このようにデータの種類（観測点）が増えていくほど、配列の次元も増えていきます。

多次元配列を平坦化する

配列には次元があり、多次元のままでは比較ができません。そこで、多次元のデータを平坦化して１次元にしてやります。
※この辺りは自分もよく分かっていませんが、仮説検定のいくつかは平坦化したデータでしか実施できないものがありました。
practice.csvというファイルのデータAの中から2列目を抽出して、1から30行目を読み込みA1to30として名前をつけます。

A = pd.read_csv('practice.csv',usecols=[1])
A1to30 = np.array(A[0:29])

これでは、行が1から30ある多次元のデータ（独立したデータ）ということになってしまうので、これを平坦化します。

A1to30fl = A1to30.fatten

これでpracticr.csvから抜き出した1から30行のデータが1行の平坦化されたデータとしてA1to30flとして格納されました。

いよいよ統計処理（記述統計？）

平均 / 中央値 / 最頻値 / 分散 / 不偏分散 / 標準偏差の計算

まずは記述統計の範囲（？）です。簡単な部分から触れていきたいと思います。
ここではAというデータに対して処理を行なっていきます。

#ここでは、データフレームAdfから配列Aをすでに準備しているものとします。

import numpy as np
import statistics 
#numpyとstatisticsをimportしています。

np.mean(A) #配列Aの平均
np.median(A) #配列Aの中央値
statistics.mode(Adf) #データフレームAの最頻値
np.var(A,ddof=0) #配列Aの分散
np.var(A,ddof=1) #配列Aの不偏分散
np.std(A,ddof=1) #配列Aの標準偏差

以上です。
データさえ整っていれば、計算は一瞬で終わります。

ヒストグラム・qqプロット

視覚的に正規性を確認する場合、まずはヒストグラムとqqプロットします（マイルールなのかな？）

#ここでは、データフレームAdfから配列Aをすでに準備していて、Aを平坦化した配列Aflが準備されているものとします。

import matplotlib.pyplot as plt
from scipy import stats
#matplotlib.pyplotをpltとしてscipyからstatsをimportしています。

plt.hist(Afl) #ヒストグラムが描けます。これはデータフレームでもOK。

fig = plt.subplots(figsize=(8,8))
stats.probplot(Afl, dist="norm", plot=plt) #qqプロットが描けます。平坦化されたデータのみ。

plt.show() #作成したヒストグラム・qqプロットを表示する。

やはり、データが整ってさえいれば後は簡単ですね。
plt.show()で指示を出せば一気にヒストグラムとqqプロットの結果を表示してくれます。
ヒストグラムは裾の長さが等しい山形に近いか？qqプロットは直線上か？で正規性が視覚的に分かりそうですね。

複数のデータを比べる（箱ひげ図と散布図）

ここまでは、1つの標本について正規性だとか平均だとかを見てきましたが、多数の標本について比べる場合、
まずは箱ひげ図を書いてみたり、散布図を書いて「どういうデータなのかなー」とイメージを持つようにしています。
余計な憶測をすることになるかもしれないけど…

#ここでは、平坦化された配列Afl、Bfl、Cfl、Dflが準備されているものとします。

import matplotlib.pyplot as plt
#matplotlib.pyplotをpltとしてimport

plt.boxplot([Afl, Bfl,Cfl,Dfl]) #箱ひげ図が描ける。

plt.scatter(Afl, Bfl),plt.title("散布図のタイトル"),plt.xlabel("X軸の名称"),plt.ylabel("Y軸の名称"),plt.grid(True) 
#散布図が描ける。複数の処理は,で区切ってやるといい。

箱ひげ図や散布図を見るとなんとなーくこれから比べたいデータの傾向が見えてきます。

仮説検定の実施

さて、複数のデータを比べることにしたいので有意差が気になります。
ここで扱うのは、自分が普段使うような仮説検定だけになります。
そもそも自分のための備忘録か…

正規性の検定（シャピロウィルク検定とコルモゴロフ・スミルノフ検定）

データを比べる前に正規性の検定から。
一般にサンプルサイズが小さい時には、シャピロウィルク検定、1000を越える時にはコルモブロフ・スミルノフ検定を用いるようです。

#ここでは、  配列Aとそれを平坦化した配列Aflが準備されているものとします。

from scipy import stats 
#scipyからstatsをimportします。

stats.shapiro(A) #シャピロウィルク検定

stats.kstest(Afl, "norm") #コルモブロフ・スミルノフ検定

statsを使えば正規性の検定は簡単にできますね。p値を確認して正規性がありそうか確認できそうです。

相関に関する検定

相関係数に関して確認しておきたいこと
相関係数の基準 0～0.3未満：ほぼ無関係
0.3～0.5未満：非常に弱い相関
0.5～0.7未満：相関がある
0.7～0.9未満：強い相関
0.9以上：非常に強い相関
決定係数は、相関係数の２乗値で１に近いほど回帰直線の説明力が高い。 これを踏まえて、相関を見ていきます。

#ここでは、平坦化された配列AflとBflが準備されているものとします。  

import numpy as np 
from scipy.stats import spearmanr 
#numpyとscipy.statからspearmanrをimportします。  

np.corrcoef(Afl, Bfl) #ピアソンの積率相関係数  

spearmanr(Afl, Bfl) #スピアマンの順位相関分析

ピアソンは正規性が仮定できる場合に、スピアマンは正規性が仮定できない場合に用いるようです。

等分散性の検定

多群の等分散性の検定のバートレットとルビーンしかないけど、統計的には問題ない（らしい）。

#ここでは、平坦化された配列Afl、Bfl、Cfl、Dflが準備されているものとします。  

from scipy import stats 
#scipyからstatsをimportします。

stats.bartlett(Afl,Bfl,Cfl,Dfl) #バートレット検定

stats.levene(Afl, Bfl,Cfl,Dfl) #ルービン検定

バートレットは正規性が仮定できる時、ルービンは正規性が仮定できない時の検定。
pythonはF検定をカバーしていないようなので、これらで検定します。

正規性・等分散性が仮定できる場合の仮説検定（パラメトリック）

やっと本題というか、なんというか…
ここで一旦区切ろうとも思ったのですが、スクリプト見ているとそこまで長い説明にはならなそうなので、一気に書いていきます。

t検定

まずは、２群の比較です。t検定には対応がある場合とない場合がありますね。
対応があるというのは、
例えば、同じ人の血糖値が薬を飲む前と後で差があるか？というのを調べるような時です。
対応がないというのは、
例えば、２人の血糖値に差があるか？というのを調べるような時です。

#ここでは、平坦化された配列AflとBflが準備されているものとします。

from scipy import stats
#scipyからstatsをimportします。  

stats.ttest_rel(Afl,Bfl) #対応あり

stats.ttest_ind(Afl,Bfl) #対応なし

これだけです。p値で仮説を確かめましょう。

チューキー・クレイマー検定

Rと同じですが、データにそれぞれラベルをつけていく必要があります。
Rと同じとか言ってもあれか…

#ここでは、平坦化された配列Afl、Bfl、CflとDflがそれぞれサンプルサイズ30で準備されているものとします。

import numpy as np
from statsmodels.stats.multicomp import pairwise_tukeyhsd
#numpyとstatsmodels.stats.multicompからpairwise_tukeyhsdをimportします。

name = np.array(["A","B","C","D"]) #データの名前のセットをnameとして作る。
name_rep = np.repeat(name,[len(Afl),len(Bfl),len(Cfl),len(Dfl)]) 
#作ったセットをそれぞれの配列分の長さ分繰り返す配列をname_repとして作る。

value = np.array([Afl,Bfl,Cfl,Dfl]) #Afl、Bfl、CflとDflを一つの配列にvalueとしてまとめる。
#valueはこのままでは4次元配列。
valuefl = value.reshape(1,120)
#valueを1行120列（４配列がそれぞれサンプルサイズ30なので120）にvalueflとして平坦化します。

pairwise_tukeyhsd(valuefl[0],name_rep).summary() #チューキー・クレイマー検定

チューキー・クレイマーをするには少しテクニック？というかなんか小手先で対応する必要がありました。
Rで多重比較するときもラベルつけたりベクトルをいじったりするので、それと同じだと思います。
あとは、summary()をつけないと結果が表示されないので注意です。

正規性・等分散性が仮定できない場合の仮説検定（ノンパラメトリック）

いつも感動するのは、正規性も等分散性も仮定できないのに標本から母数が推測できる仕組みを作った人がいるんだな、ということです。
なんか、雑多なデータというか規則性がないように見えるデータにも法則を見出すことができるのは素直にすごいと思います。

welchのt検定

正規性は仮定できるけど等分散性は仮定できないような場合にするt検定です。

#ここでは、平坦化された配列AflとBflが準備されているものとします。

from scipy import stats
#scipyからstatsをimportします。  

stats.ttest_ind(Afl,Bfl, equal_var=False) #welchのt検定

内容はすごく簡単で、t検定にequal_varがfalseですよ、つまり「等分散じゃないですよ」と指示するだけです。

ノンパラメトリックだけど対応あるなし２群比較

正規性も等分散性も仮定できないけれど、2群の比較をしたい場合に
対応がある場合、ウィルコクソンの符号順位検定
対応がない場合、マンホイットニーのU検定をします。

#ここでは、平坦化された配列AflとBflが準備されているものとします。

from scipy import stats
#scipyからstatsをimportします。

stats.wilcoxon(Afl, Bfl, alternative="two-sided").pvalue #ウィルコクソンの符号順位検定

stats.mannwhitneyu(Afl, Bfl, alternative="two-sided").pvalue #マンホイットニーのU検定

ここまでで２群の比較はほぼ説明し終わりました。

ノンパラメトリックの多重比較

と、言いたいところですが、
チューキー・クレイマーのノンパラ版的なスティール・ドゥワスの検定はpackageがあるようですが実行できず。
多群と１つの対照群を比較するダネットの検定は探してもどうにも出てきませんでした。残念。
今回は、ここまでです。少し悔しいですが、RならサクッとできるのでやっぱりRの方が好き！という感じです。

まとめ

データの準備からノンパラメトリックな２群の比較までサラッと小手先で描いてみました。
Rをいじっている時やExcelでデータをなんやかんやしている時から感じていることですが、
データをうまく準備できさえすれば、あとはpackage任せです。
これでいいとは思えませんが、自分には０からスクリプトを書く技量もありませんし、とりあえずはこれでp値を見ていこうかな？と思います。
そもそもR studioが使えさえすれば見向きもしなかったpythonなので、一つ良い勉強になりました。
ノンパラメトリックな多重比較はRでやりたいし、そもそも検定ならRのほうが直感的で使いやすい気がしました。
また勉強し直したら全く逆のことを思っている可能性もありますが…

ちょっと調べるとpythonは予測することに長けていて、機械学習だとかそういうところで真価を発揮するようです。
Rは逆にこういった仮説検定とかすでにあるデータの説明が得意だとか。
今回仮説検定をpythonでできるようにしてみましたが、なんだかどのサイトもさっぱりしているというか、Rほど体系的にまとまっていないような印象を受けたのは、そのためかもしれません。

いつもなら「このサイトを引用したよ」とか「このサイトが分かりやすいよ」という引用元を示すのですが、
あまりにも多くのサイトを参考にしながらまとめたので、もう何がなんだか、どのサイトを見たんだか見ていないんだか分からないので、
今回は出来ませんでした…
もし、剽窃だと感じて「これはここにあるじゃないか！」と思う人がいたら是非コメントしてください。
ソース元を書く方が役に立つことの方が多いので。
あー疲れた！今回はおしまい。
次は、VSCodeでRの統計処理をしようかなと思ってます。
macが古すぎて、Rstudioが使えないのでね…