分散分析
被験者間１要因（独立なサンプル）、被験者内１要因（繰り返し／対応のあるデータ）、２要因分散分析のPythonスクリプト例を作成した。スクリプトファイルは、１要因の場合のものを圧縮ファイルANOVA.zipにまとめた。Jupyter Notebook用のものは圧縮ファイルANOVA_jn.zipにまとめた。それぞれダウンロード・解凍して自由に利用されたい。２要因のものは、２要因のところに用意している。
被験者間１要因デザインのデータ例として、表１のものを考える。条件として水準１、水準２、水準３の３水準がせっていされ、各水準内の条件で独立に３個のデータが収集されている。
|
表１ 被験者間１要因分散分析データ例 |
|||
|
条件 |
データ |
||
|
水準０ |
５ |
７ |
６ |
|
水準１ |
９ |
８ |
１０ |
|
水準２ |
６ |
３ |
３ |
表１のデータの分散分析のためのスクリプトを、＜岡本安晴「データを使いこなすための統計入門」２２世紀アート＞のpp. 99-107の説明に従って、以下のスクリプトを作成した。
import numpy as np
import scipy.stats as ss
Y = [[5, 7, 6],
[9, 8, 10],
[6, 3, 3]]
a = len(Y)
N = len(Y[0])
print("a = %d"%a)
print("N = %d"%N)
Y = np.array(Y)
Yd = Y.mean(axis = 1)
print('Yd = ', Yd)
Ydd = Y.mean()
print('Ydd = ', Ydd)
SSA = ((Yd - Ydd) ** 2).sum() * N
print('SSA = ', SSA)
SSS_A = 0.0
for j in range(a):
SSS_A += ((Y[j] - Yd[j]) **
2).sum()
print('SSS_A = ', SSS_A)
SST = SSA + SSS_A
dfA = a - 1
dfS_A = a * (N - 1)
MSA = SSA / dfA
MSS_A = SSS_A / dfS_A
F = MSA / MSS_A
p = ss.f.cdf(F, dfA, dfS_A)
print()
print("%10s"%"Source",
end = " ")
print("%12s"%"SS",
end = " ")
print("%5s"%"df",
end = " ")
print("%12s"%"MS",
end = " ")
print("%12s"%"F")
print("%10s"%"A",
"%12.3f"%SSA, "%5d"%dfA, "%12.3f"%MSA,
"%12.3f"%F,
" p = %-10.5f"%(1.0 - p))
print("%10s"%"S/A",
"%12.3f"%SSS_A, "%5d"%dfS_A, "%12.3f"%MSS_A)
print("%10s"%"Total",
"%12.3f"%SST, "%5d"%(a * N - 1))
表１のデータをリストに
Y = [[5, 7, 6],
[9, 8, 10],
[6, 3, 3]]
と設定している。各行が各水準のデータに対応していることに注意。リストYに他のデータを設定すれば、そのデータの分析ができる。スクリプトは圧縮ファイルanovaOneFacInd.zipをダウンロード・解凍して、自由に利用されたい。
上のプログラムを実行すると以下の結果（端末画面からコピー・貼り付け）を得る。
====== RESTART:
D:\yasuharu\www\LaCoocan\Python\IntroStat\ANOVA\anova_files\anovaOneFacInd.py
======
a = 3
N = 3
Yd = [6. 9. 4.]
Ydd = 6.333333333333333
SSA = 38.0
SSS_A = 10.0
Source
SS df
MS
F
A
38.000 2
19.000
11.400 p =
0.00904
S/A
10.000 6
1.667
Total
48.000 8
>>>
被験者内１要因デザインの分散分析のデータ例として表２のものを考える。条件として水準０、水準１、水準２の３水準があり、被験者は各水準全ての条件においてデータが採られる。
|
表２ 被験者内１要因分散分析データ例 |
|||
|
|
条件 |
||
|
被験者 |
水準０ |
水準１ |
水準２ |
|
被験者０ |
１ |
３ |
２ |
|
被験者１ |
５ |
９ |
７ |
|
被験者２ |
６ |
１２ |
９ |
表２のデータを分散分析するためのスクリプトを、＜岡本安晴「データを使いこなすための統計入門」２２世紀アート＞のpp. 107-115の説明に従って、次のように作成した。
import numpy as np
import scipy.stats as ss
Y = [[1, 3, 2],
[5, 9, 7],
[6, 12, 9]]
N = len(Y)
a = len(Y[0])
print("a = %d N = %d"%(a, N))
Y = np.array(Y)
Yid = Y.mean(axis = 1)
print(Yid)
Ydj = Y.mean(axis = 0)
print(Ydj)
Ydd = Y.mean()
print(Ydd)
SS_A = ((Ydj - Ydd) ** 2).sum() * N
SS_S = ((Yid - Ydd) ** 2).sum() * a
SS_SA = ((Y - np.array([Yid]).T - Ydj + Ydd) **
2).sum()
SS_T = SS_A + SS_S + SS_SA
df_A = a - 1
df_S = N - 1
df_SA = (a - 1) * (N - 1)
MS_A = SS_A / df_A
MS_S = SS_S / df_S
MS_SA = SS_SA / df_SA
F = MS_A / MS_SA
p = ss.f.cdf(F, df_A, df_SA)
print()
print("%10s"%"Source",
"%12s"%"SS", "%5s"%"df",
"%12s"%"MS", "%12s"%"F")
print("%10s"%"S",
"%12.3f"%SS_S, "%5d"%df_S, "%12.3f"%MS_S)
print("%10s"%"A",
"%12.3f"%SS_A, "%5d"%df_A, "%12.3f"%MS_A,
"%12.3f"%F,
" p = %-10.5f"%(1.0 - p))
print("%10s"%"SA",
"%12.3f"%SS_SA, "%5d"%df_SA, "%12.3f"%MS_SA)
表２のデータをリストに
Y = [[1, 3, 2],
[5, 9, 7],
[6, 12, 9]]
と設定している。他のデータをYに設定すれば、そのデータの分析が行える。スクリプトファイルは圧縮ファイルanovaOneFacRel.zipとして用意したが、ダウンロード解凍して自由に利用されたい。
上のスクリプトを実行すると、以下の出力を得る（端末画面のコピー・貼り付け）。
====== RESTART: D:/yasuharu/www/LaCoocan/Python/IntroStat/ANOVA/anova_files/anovaOneFacRel.py
======
a = 3 N = 3
[2. 7. 9.]
[4. 8. 6.]
6.0
Source
SS df
MS
F
S
78.000 2 39.000
A
24.000 2
12.000 12.000 p = 0.02041
SA
4.000 4
1.000
>>>
2要因分散分析
２要因分散分析のスクリプト例は別のウェブサイトに用意した。