パッケージのインストールはpip3 installで
pythonでの仮説検定も基本的にはRと同じでpackage頼りです。
肝になる部分は「どのpackageをどこで使うの?」ということだと思います。
簡単な仮説検定まで一気に紹介していきましょう。
スクリプトの詳しい内容や検定の説明については、統計の参考書や他のサイトに任せて、
「なんだか分からないけどできるわ」を目指しました。
検定する前の準備
「よし、検定しよう!」と思っても何かしらデータがなければ意味がありません。
ここで、データの整理の仕方を確認していきましょう。
必要そうなpackageはimportしていく
Rと違う点は、まずpythonファイル上でpackageをimportする指示を出すところです。
installさえされていれば、import ほにゃららと書くだけで非常に簡単です。
import pandas as pd
import numpy as np
import statistics
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import scikit_posthocs as sp
お作法がありそうです
例えば、
import pandas
とすれば、pandasがimportされ指示が出せますが、pdとしてimportするのが普通みたいです。
numpy (ナンピーじゃなくてナンパイなんだってさ)の場合もnpとしてimportします。
上述の通りですが、名前を付けてimportする場合は、as ほにゃほにゃとしてimportしておきます。
正直、この辺はよく分かっていないので、周囲の人に聞くか仕様書に合わせて書くのが無難そうですね。
(他の人が見たときに「なんじゃこりゃ」となると非常に迷惑をかけるので)
csvファイルの読み込み
pythonでcsvファイルを読み込むとき、2行目の列だけ読み込みたい場合は、
pd.read_csv('practice.csv',usecols=[1])
ちなみに…毎回自分もこんがらがるのですが「行が横(水平方向)」で「列が縦(垂直方向)」らしいです。
英語では「行がrow」で「列がcolumn」です。
じゃぁ、userowsで指定した行が読み込まれるのか?というとそういう訳ではないようです。
データフレームと次元の平坦化
前述のpractice.csvというファイルを読み込む時にこれにAという名前をつけるとします。
A = pd.read_csv('practice.csv',usecols=[1])
このうちの1-30行のデータをA1to30として読み込みたい場合、numpyを使います。
numpyはnpとしてimportしているので、
A1to30 = np.array(A[0:29])
次に抜き出したA1to30という配列からA1to30dfというデータフレームを作るとします。
A1to30のDataFrameなのでA1to30dfです。
A1to30df = pd.DataFrame(data=A1to30)
説明しやすいかな?と思ってAとかA1to30とか直感的な名前を使ってしまったので少し冗長な説明になってしまいましたね。
配列の次元?
配列には次元という考え方があります。
例えば、 AさんとBさんとCさんの走った距離について考えましょう。
「Aさんは3km走った」「Bさんは1.5km」「Cさんは6km走った」という場合、1行(行は横でした)だけで管理できそうですね。
これが1次元データです。
ここにAさん、Bさん、Cさんが走るのにかかった時間をデータとして足してみましょう。
「Aさんは3kmを30分で走った」「Bさんは1.5kmを10分で走った」「Cさんは6kmを60分で走った」とすると、
今度は、2行で管理することになります。これが2次元データです。
さらにAさん、Bさん、Cさんが走りながら飲んだ水の量をデータとして足すとどうなるでしょうか。
「Aさんは3kmを30分で走りながら100ml飲んだ」「Bさんは1.5kmを10分で走りながら50ml飲んだ」「Cさんは6kmを60分で走りながら150ml飲んだ」
となると、今度は3行のデータになります。これが3次元データです。
このようにデータの種類(観測点)が増えていくほど、配列の次元も増えていきます。
多次元配列を平坦化する
配列には次元があり、多次元のままでは比較ができません。そこで、多次元のデータを平坦化して1次元にしてやります。
※この辺りは自分もよく分かっていませんが、仮説検定のいくつかは平坦化したデータでしか実施できないものがありました。
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というデータに対して処理を行なっていきます。
import numpy as np
import statistics
np.mean(A)
np.median(A)
statistics.mode(Adf)
np.var(A,ddof=0)
np.var(A,ddof=1)
np.std(A,ddof=1)
以上です。
データさえ整っていれば、計算は一瞬で終わります。
ヒストグラム・qqプロット
視覚的に正規性を確認する場合、まずはヒストグラムとqqプロットします(マイルールなのかな?)
import matplotlib.pyplot as plt
from scipy import stats
plt.hist(Afl)
fig = plt.subplots(figsize=(8,8))
stats.probplot(Afl, dist="norm", plot=plt)
plt.show()
やはり、データが整ってさえいれば後は簡単ですね。
plt.show()で指示を出せば一気にヒストグラムとqqプロットの結果を表示してくれます。
ヒストグラムは裾の長さが等しい山形に近いか?qqプロットは直線上か?で正規性が視覚的に分かりそうですね。
複数のデータを比べる(箱ひげ図と散布図)
ここまでは、1つの標本について正規性だとか平均だとかを見てきましたが、多数の標本について比べる場合、
まずは箱ひげ図を書いてみたり、散布図を書いて「どういうデータなのかなー」とイメージを持つようにしています。
余計な憶測をすることになるかもしれないけど…
import matplotlib.pyplot as plt
plt.boxplot([Afl, Bfl,Cfl,Dfl])
plt.scatter(Afl, Bfl),plt.title("散布図のタイトル"),plt.xlabel("X軸の名称"),plt.ylabel("Y軸の名称"),plt.grid(True)
箱ひげ図や散布図を見るとなんとなーくこれから比べたいデータの傾向が見えてきます。
仮説検定の実施
さて、複数のデータを比べることにしたいので有意差が気になります。
ここで扱うのは、自分が普段使うような仮説検定だけになります。
そもそも自分のための備忘録か…
正規性の検定(シャピロウィルク検定とコルモゴロフ・スミルノフ検定)
データを比べる前に正規性の検定から。
一般にサンプルサイズが小さい時には、シャピロウィルク検定、1000を越える時にはコルモブロフ・スミルノフ検定を用いるようです。
from scipy import stats
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以上:非常に強い相関
決定係数は、相関係数の2乗値で1に近いほど回帰直線の説明力が高い。
これを踏まえて、相関を見ていきます。
import numpy as np
from scipy.stats import spearmanr
np.corrcoef(Afl, Bfl)
spearmanr(Afl, Bfl)
ピアソンは正規性が仮定できる場合に、スピアマンは正規性が仮定できない場合に用いるようです。
等分散性の検定
多群の等分散性の検定のバートレットとルビーンしかないけど、統計的には問題ない(らしい)。
from scipy import stats
stats.bartlett(Afl,Bfl,Cfl,Dfl)
stats.levene(Afl, Bfl,Cfl,Dfl)
バートレットは正規性が仮定できる時、ルービンは正規性が仮定できない時の検定。
pythonはF検定をカバーしていないようなので、これらで検定します。
正規性・等分散性が仮定できる場合の仮説検定(パラメトリック)
やっと本題というか、なんというか…
ここで一旦区切ろうとも思ったのですが、スクリプト見ているとそこまで長い説明にはならなそうなので、一気に書いていきます。
t検定
まずは、2群の比較です。t検定には対応がある場合とない場合がありますね。
対応があるというのは、
例えば、同じ人の血糖値が薬を飲む前と後で差があるか?というのを調べるような時です。
対応がないというのは、
例えば、2人の血糖値に差があるか?というのを調べるような時です。
from scipy import stats
stats.ttest_rel(Afl,Bfl)
stats.ttest_ind(Afl,Bfl)
これだけです。p値で仮説を確かめましょう。
チューキー・クレイマー検定
Rと同じですが、データにそれぞれラベルをつけていく必要があります。
Rと同じとか言ってもあれか…
import numpy as np
from statsmodels.stats.multicomp import pairwise_tukeyhsd
name = np.array(["A","B","C","D"])
name_rep = np.repeat(name,[len(Afl),len(Bfl),len(Cfl),len(Dfl)])
value = np.array([Afl,Bfl,Cfl,Dfl])
valuefl = value.reshape(1,120)
pairwise_tukeyhsd(valuefl[0],name_rep).summary()
チューキー・クレイマーをするには少しテクニック?というかなんか小手先で対応する必要がありました。
Rで多重比較するときもラベルつけたりベクトルをいじったりするので、それと同じだと思います。
あとは、summary()をつけないと結果が表示されないので注意です。
正規性・等分散性が仮定できない場合の仮説検定(ノンパラメトリック)
いつも感動するのは、正規性も等分散性も仮定できないのに標本から母数が推測できる仕組みを作った人がいるんだな、ということです。
なんか、雑多なデータというか規則性がないように見えるデータにも法則を見出すことができるのは素直にすごいと思います。
welchのt検定
正規性は仮定できるけど等分散性は仮定できないような場合にするt検定です。
from scipy import stats
stats.ttest_ind(Afl,Bfl, equal_var=False)
内容はすごく簡単で、t検定にequal_varがfalseですよ、つまり「等分散じゃないですよ」と指示するだけです。
ノンパラメトリックだけど対応あるなし2群比較
正規性も等分散性も仮定できないけれど、2群の比較をしたい場合に
対応がある場合、ウィルコクソンの符号順位検定
対応がない場合、マンホイットニーのU検定をします。
from scipy import stats
stats.wilcoxon(Afl, Bfl, alternative="two-sided").pvalue
stats.mannwhitneyu(Afl, Bfl, alternative="two-sided").pvalue
ここまでで2群の比較はほぼ説明し終わりました。
ノンパラメトリックの多重比較
と、言いたいところですが、
チューキー・クレイマーのノンパラ版的なスティール・ドゥワスの検定はpackageがあるようですが実行できず。
多群と1つの対照群を比較するダネットの検定は探してもどうにも出てきませんでした。残念。
今回は、ここまでです。少し悔しいですが、RならサクッとできるのでやっぱりRの方が好き!という感じです。
まとめ
データの準備からノンパラメトリックな2群の比較までサラッと小手先で描いてみました。
Rをいじっている時やExcelでデータをなんやかんやしている時から感じていることですが、
データをうまく準備できさえすれば、あとはpackage任せです。
これでいいとは思えませんが、自分には0からスクリプトを書く技量もありませんし、とりあえずはこれでp値を見ていこうかな?と思います。
そもそもR studioが使えさえすれば見向きもしなかったpythonなので、一つ良い勉強になりました。
ノンパラメトリックな多重比較はRでやりたいし、そもそも検定ならRのほうが直感的で使いやすい気がしました。
また勉強し直したら全く逆のことを思っている可能性もありますが…
ちょっと調べるとpythonは予測することに長けていて、機械学習だとかそういうところで真価を発揮するようです。
Rは逆にこういった仮説検定とかすでにあるデータの説明が得意だとか。
今回仮説検定をpythonでできるようにしてみましたが、なんだかどのサイトもさっぱりしているというか、Rほど体系的にまとまっていないような印象を受けたのは、そのためかもしれません。
いつもなら「このサイトを引用したよ」とか「このサイトが分かりやすいよ」という引用元を示すのですが、
あまりにも多くのサイトを参考にしながらまとめたので、もう何がなんだか、どのサイトを見たんだか見ていないんだか分からないので、
今回は出来ませんでした…
もし、剽窃だと感じて「これはここにあるじゃないか!」と思う人がいたら是非コメントしてください。
ソース元を書く方が役に立つことの方が多いので。
あー疲れた!今回はおしまい。
次は、VSCodeでRの統計処理をしようかなと思ってます。
macが古すぎて、Rstudioが使えないのでね…