こんにちは。からあげ博士(@phd_karaage)です。今回は分散分析を取り扱おうと思います。それはなぜか。単純にTAで分散分析に関する授業のレポートの採点が回ってきたからですね。模範解答を作って、その上で採点をしていく訳ですが、せっかくなら分散分析の手法的な部分について記事にまとめておこうと思ったからです。
データをいろいろ見ていくと、平均や標準偏差は分かった、じゃあこの母集団の平均値に差があるの？という問題に直面します。ぱっと見で平均値の棒グラフを見てみて差があるかどうかなんて、あまり当てにならない直感でしかない訳です。逆に消費者にそれだけの知識がないだろうと踏んだ、最低で最悪なグラフがこの世には蔓延る訳ですが。
という訳でirisデータセットを使って、Rとpythonを使った分散分析について紹介しようと思います。これまでにirisデータセットを使った話はこちら。
目次
分散分析(ANOVA)とは？
分散分析とはなにか。それは得られたデータの変動が、処理に由来するのか、誤差に由来するのか評価を行う手法です。irisデータセットでも、品種ごとに平均値が微妙に異なりましたね。この平均値の差異が品種によって生じたものか、それとも誤差（環境効果など）によって生じたものか評価してあげる必要があります。最も単純な分散分析は単純に母平均の差を比較するt検定（スチューデントのt検定）と言われています。
分散分析の流れとしては、F値を求めてそれに対して検定するという流れになります。
実際にpythonで分散分析を実行してみた
from sklearn.datasets import load_iris import pandas as pd iris_dataset = load_iris() data_array = iris_dataset["data"] line_num = iris_dataset["target"] line = iris_dataset["target_names"] feature = iris_dataset["feature_names"] df_data = pd.DataFrame(data_array, columns = feature) df_line = pd.DataFrame(line_num) df = pd.concat([df_data, df_line], axis = 1) df = df.rename(columns = {0: "Target"}) df_line_0 = df[df.Target == 0] df_line_1 = df[df.Target == 1] df_line_2 = df[df.Target == 2]
それでは下準備として、品種ごとにデータを分けましょう。データの読み込み部分については、これまで紹介してきた方法と同じです。品種が同じものだけ各データフレームから抽出し、それぞれdf_line_0とdf_line_2として格納しておきます。
今回はとりあえずがく片の長さ（sepal length (cm))について分散分析を行ってみましょう。
from scipy.stats import f_oneway sep_len_0 = df_line_0["sepal length (cm)"] sep_len_1 = df_line_1["sepal length (cm)"] sep_len_2 = df_line_2["sepal length (cm)"] result = f_oneway(sep_len_0, sep_len_1, sep_len_2) print(result)
今回は一元配置分散分析を行うので、scipyから stat.f_onewayをインポートしておきます。f_onewayは返り値にF値とp値を取るので、そのままresultに格納しておいて最後にprintで表示しています。
結果はこんな感じ。F値が119.264で、p値が1.67×10の-31乗ということでp値は1%有意水準を遥かに下回っていますね。ということで帰無仮説が棄却できます。この一元配置分散分析における帰無仮説は「すべてのデータが同じ母平均を持つ集団から得られた」という仮説なので、帰無仮説が棄却されるということは、「それぞれのデータは異なる母平均を持つ集団から得られた」と言える訳です。このデータはそれぞれ異なる品種から取られている訳ですから、なるほど。という結果ですね。
ちなみに分散分析は基本的に3群以上のデータについて検定するときに用います。このf_onewayは2群だけを入れても動作しますが、結局それはt検定の結果と変わらないということになります。
ただ、これだと分散分析の結果は見れても、分散分析表を見せてはくれません。pythonでもstatmodelsを使うと求められるようですが、ちょっと面倒なのでここでは扱いません。というか使ったことがない。
統計使うならRのほうが便利だよね？
絶妙にかゆいところに手が届くようで届かないpythonの統計解析ですが、Rを使うと案外楽勝だったりします。自分もなるべくpythonを使いたい派の人間ですが、統計解析諸々が必要になるとRを使います。RとRStudioは便利。
という訳でRでもデータを用意しておきましょう。同じirisデータセットがRの中にもベース関数として入っており簡単に取り出すことができます。何もライブラリを読み込まなくてもいいというのはいいのか悪いのか、なんとも言い難いところではありますが。
data_iris <- iris summary(aov(lm(data_iris$Sepal.Length ~ data_iris$Species)))
なにも考えずとも、一行でデータフレームを抽出できます。最高かよ。そして次の1行で完結しています。やはり統計解析言語と言われるだけあって楽勝です。
返ってきた結果はこんな感じ。一撃で分散分析表を求めてくれるあたり賢い。F値はpythonで求めたものと同じ。p値はpythonで求めたものとは異なりますが、十分に小さいと出ています。Rではpythonがp値を正確に表示する一方で、十分に小さい値の場合、2.2e-16としか返してくれなくなります。この挙動は果たしていいものか。
さてこのコードを読んで、ん？と思った人もいるかもしれませんね。ここではlmを使った回帰を行っています。回帰を行い、それについて分散分析をaov行ってその結果をsummaryで表示していると。やっていることはそんな感じです。
今回線形回帰lmを使って線形モデルにあてはめて、その結果について分散分析を行っています。単純にaovを使おうとすると、Speciesが数値データではなく、文字列なので怒られてしまうんですね。だからとりあえずモデルにあてはめようということをしている訳です。
線形回帰と分散分析はある種共通のことをしている訳なので、こんなこともできるよということですね。
統計解析はRのほうが楽だと思う
研究室ではもっぱらpythonを使うことが多く、統計解析を行うときに関してだけRを使ってますが、Rを使いこなせるようになることで様々な統計的手法を身に着けることができるのではないかなという気がしています。線形回帰も、分散分析も、主成分分析も、ここでは取り扱っていない統計的手法の多くはRだとするっと書いて実行できるものが多いです。
割と使い心地がいいなあと思うpythonに逃げがちではありますが、こうしてRを使ってみると、Rってやっぱいいなと思うこともあります。今回ちょっとそれを再確認しましたね。データサイエンスをやるなら、pythonもいいけどRも身につけましょう。なんてことを思ったのでした。
