Pythonで分散分析（対応なし・二元配置）

# <font color = "Green">Python</font>で<font color = "deeppink">分散分析</font>（対応なし・二元配置）

# 目的
<font color = "BLUE">SPSS</font>で実施していた対応なしの二元配置分散分析を無料の<font color = "Green">Python</font>で実行すること．



## 参考にしたwebページの紹介
本記事の内容のほとんどは<a href="https://www.marsja.se/three-ways-to-carry-out-2-way-anova-with-python/" target="_blank">Erik Marsja「Three ways to do a two-way ANOVA with Python」</a>の記事を参考にさせていただきました． 
また，二要因分散分析の復習は<a href="https://bit.ly/2EN45V4" target="_blank">統計検定を理解せずに使っている人のためにⅢ</a>のページを参考にしました．
途中で行う多重比較についてはこの<a href="https://bit.ly/2ShbAGT" target="_blank">Pythonで統計学を学ぶ(6)</a>というページから関数を使わせてもらっています．   
## 以下に書いてあること
Pythonで対応なしの二元配置分散分析を行った道筋を，サンプルデータを使いながらまとめてみました．

### サンプルデータについて
サンプルデータは600人の若年男性を対象に100m走のタイムを計測した設定です．このとき，<font color = "RED"><b>「シューズの違い」及び「走ったトラックの違い」の2つの要因が100m走のタイムに及ぼす影響</b></font>を，pythonで二元配置分散分析をしながら統計的にみていきます．

要因は
<b>シューズ条件：A社製シューズ[A]，B社製シューズ[B]，N社製シューズ[N]の3水準
トラック条件：グラウンド[Ground]，タータン[Tartan]の2水準</b>（陸上トラックに使用されているコルクのような見た目の材質は「タータン」と呼ぶんだと陸上経験者の友人から教わりました）
の2つであり，条件数は2要因×3水準=6条件です．各条件につき100名が走っている設定です．

|  対象  	| 若年男性600名                                           	| 
|:------:	|---------------------------------------------------------	|
|  要因  	| シューズ（3水準），トラック（2水準）                    	|
|  指標  	| 100m走ったタイム[s]                                     	|
| 条件数 	| （シューズ3水準）×（トラック2水準）＝6条件（600人集めて各条件に付き100名ずつ走ってもらった設定）	|


### ライブラリのインポート
```python
	# 数値計算に使うライブラリ
	import pandas as pd
	import numpy as np
	import scipy as sp
	from scipy import stats as st

	# グラフを描画するライブラリ
	from matplotlib import pyplot as plt
	import seaborn as sns
	sns.set()
	
	# 統計モデルを推定するライブラリ
	import statsmodels.formula.api as smf
	import statsmodels.api as sm
	import statsmodels.stats.anova as anova # 分散分析やるライブラリ
	from statsmodels.stats.multicomp import pairwise_tukeyhsd # Tukeyの多重比較やるライブラリ

	from statsmodels.graphics.factorplots import interaction_plot
	from statsmodels.formula.api import ols
	from statsmodels.stats.anova import anova_lm

	# 表示桁数の指定
	%precision 3

	# グラフをjupyter notebook内に表示させるための指定
	%matplotlib inline
```
### サンプルデータの読み込み  

<a href="https://github.com/TaigaU121/Sampledata.git" target="_blank">サンプルデータ</a>も上げておきます．


```python
	file = "Runtime.csv"
	data = pd.read_csv(file)
```

```python
	data.head(10) #データの中身を確認（先頭10行のみ）
```

```:結果
Unnamed: 0	ShoesMaker	Track	Time
0	0	A	Ground	15.624345
1	1	A	Ground	13.388244
2	2	A	Ground	13.471828
3	3	A	Ground	12.927031
4	4	A	Ground	14.865408
5	5	A	Ground	11.698461
6	6	A	Ground	15.744812
7	7	A	Ground	13.238793
8	8	A	Ground	14.319039
9	9	A	Ground	13.750630
```

### グラフの表示
```python
	fig = interaction_plot(data.ShoesMaker, data.Track, data.Time, # (横軸にする要因（水準多い方），系列にする要因（水準少ない方），評価指標)
             colors=['red','blue'], markers=['D','^'], ms=10)
```

結果
![二元配置分散分析.png](https://qiita-image-store.s3.amazonaws.com/0/319114/846b8b3a-c3ff-fcae-3897-c4fef7f6f56b.png)

### グラフの解釈
<b>平均時間が短い＝速く走れている</b>
ということですので，横軸のシューズごとにみると，N社のシューズが他のシューズと比べて速く走れるシューズのようです．
また系列ごと（赤い線と青い線）にみると，Taratanで走ったほうが速く走れそうです．  

加えて，N社製のシューズかつTaratanで走ると，なんだかますます速く走れそうであること（相乗効果がありそうなこと）が，グラフの傾き具合から言えるかもしれません．換言すれば，シューズ×トラックの交互作用もあるかもしれません．
そこで，二元配置の分散分析を実施します．

### 二元配置分散分析の実施
```python
	formula = 'Time ~ (Track) + (ShoesMaker) + (Track):(ShoesMaker)'
	model = ols(formula, data).fit()
	aov_table = anova_lm(model, typ=2)
	aov_table.columns = ["平方和","自由度","F値","p値"] #列名を日本語に差し替え
	print(aov_table) #被験者間効果の検定結果を出力
```
```:結果

	                         平方和  自由度      　 F値  　　      p値
	Track                345.591633    1.0  344.753535  4.987607e-61
	ShoesMaker           211.505336    2.0  105.496206  6.228932e-40
	ShoesMaker:Track      16.505274    2.0    8.232623  2.973539e-04
	Residual             595.444018  594.0         NaN           NaN
```
シューズやトラックで主効果がみられました．さらにシューズ条件とトラック条件の間で交互作用がみられました（ <i>p</i> < 0.05）．
<b>交互作用がみられた場合，全群の多重比較を実施</b>していくのが王道なようです．
交互作用が無かった場合は，各要因ごとの水準間の平均値で（他方の要因は一旦置いといて）群間比較すれば良いそうです．
補足：<a href="https://bit.ly/2VErKwb" target="_blank">anova_lmについて</a>


### 交互作用あり→全群の多重比較（今回はシューズA時のトラックの違いのみ）

多重比較については<a href="https://qiita.com/TaigaU121/items/12c480f51a026ca9f333" target="_blank">対応なしの一元配置分散分析</a>のときも少しまとめましたので，今回はシューズAのときのトラック条件間の比較のみメモしておきます．

```python	
MakerATrackG = data[(data.ShoesMaker=="A")&(data.Track=="Ground")] # dataからShoesMakerがAかつTrackがGroundのものを抽出
MakerATrackT = data[(data.ShoesMaker=="A")&(data.Track=="Tartan")] # dataからShoesMakerがAかつTrackがTartanのものを抽出
MakerATrackG.head(10) # 中身の確認
```
```:結果
Unnamed: 0	ShoesMaker	Track	Time
0	0	A	Ground	15.624345
1	1	A	Ground	13.388244
2	2	A	Ground	13.471828
3	3	A	Ground	12.927031
4	4	A	Ground	14.865408
5	5	A	Ground	11.698461
6	6	A	Ground	15.744812
7	7	A	Ground	13.238793
8	8	A	Ground	14.319039
9	9	A	Ground	13.750630
```

### Tukeyの多重比較を実施

##### 多重比較するための関数を用意
以下の関数を使います( <a href="https://bit.ly/2ShbAGT" target="_blank">Pythonで統計学を学ぶ(6)</a> 内の「 In [42]　# 参考: 関数にしてみた」より)

```python
	# 関数(グループ数の羅列，グループ名)
	def tukey_hsd( ind, *args ): # 第1引数:名称のリスト（index）, 第2引数以降: データ(*args: 複数の引数をタプルとして受け取る)
  	  from statsmodels.stats.multicomp import pairwise_tukeyhsd
	
  	  data_arr = np.hstack( args ) # 配列の結合をhstackでしているらしい．水平方向（新しく列を増やしている方向）に結合．
  	  ind_arr = np.array([])
   	 for x in range(len(args)):
  	    ind_arr = np.append(ind_arr, np.repeat(ind[x], len(args[x]))) # ind_arrの配列の最後にnp.repeat（)を加える
  	  # np.repeat(A,N)は配列A内の各要素をN回繰り返す
    	
 	   print(pairwise_tukeyhsd(data_arr,ind_arr))
```
### A社製シューズを履いた時のトラック間の違いを比較
```python
tukey_hsd(list('GT') , MakerATrackG['Time'],MakerATrackT['Time'])
```
```:結果
Multiple Comparison of Means - Tukey HSD,FWER=0.05
=============================================
group1 group2 meandiff  lower   upper  reject
---------------------------------------------
  G      T    -1.0809  -1.3553 -0.8066  True 
---------------------------------------------
```
以上より，シューズAのときのGurand-Tartan間に有意差がありそうです．
同じことをシューズBのときのTrackごと，シューズNのときのTrackごとに実施し，要因を変えてGroundのときのシューズごと，TartanのときのTrackごとに多重比較を実施してやれば全群の多重比較が網羅できると思います．



#結論
Pythonで対応なしの二元配置の分散分析ができました！


##### 終わりに
ここまでやっといてなんですが，<a href="https://docs.scipy.org/doc/scipy/reference/stats.html" target="_blank">scipy公式（多分）</a>は

>For many more stat related functions install the software R and the interface package rpy.
統計にまつわるその他の関数については，Rやrpyというパッケージ（Rとやり取りできるらしい）をインストールしましょう

とのことでしたので，分散分析関連もRの力を借りた方が良いのかもしれません（無念）．