被験者間１要因・被験者内１要因分散分析
Pythonスクリプトファイルとサンプルデータファイルを圧縮ファイルanova_one_ind_one_rel.zipにまとめた。Jupyter Notebook用は圧縮ファイルanova_one_ind_one_rel_jn.zipとした。それぞれダウンロード・解凍して自由に利用されたい。
被験者間要因（非繰り返し要因）Ａ、および被験者内要因（繰り返し要因）Ｂのモデルを次のようにおく。
ここで、主効果Ａ
主効果Ｂ
交互作用ＡＢ
ij番目の被験者固有の効果
ＡjとＢkとの交互作用による効果を調整した、ij番目の被験者と要因Ｂの水準kとの交互作用
である。
データ値
の平方和を次のように分解する。
ここで
ただし、
である。
これらの平方和の期待値は以下のようになる。
ここで、
である。
以上より、分散分析表は表１のように作成される。
表１ 分散分析表
|
変動要因 |
平方和 |
自由度 |
平均平方和 |
Ｆ比 |
|
被験者間 |
|
|
|
|
|
Ａ |
|
|
|
|
|
残差（Ａ） |
|
|
|
|
|
被験者内 |
|
|
|
|
|
Ｂ |
|
|
|
|
|
Ａ×Ｂ |
|
|
|
|
|
残差（Ｂ） |
|
|
|
|
|
全体 |
|
|
|
|
上のモデルに基づいて、リスト１に示すスクリプトを用意した。このスクリプトのファイルとサンプルデータファイルはまとめて圧縮ファイルfiles.zipとした。
リスト１ 被験者間１要因被験者内１要因分散分析Pythonスクリプト例
import csv
import numpy as np
fin_name = input('Data file
name(*.csv) = ')
with open(fin_name, 'r') as f:
data = [v for v in
csv.reader(f)]
fout_name = input('Output file
name(*.txt) = ')
fout = open(fout_name, 'w')
fout.write('Data file name =
{}\n'.format(fin_name))
print('Data file...')
for v in data:
for w in v:
print('{0: >7s}'.format(w), end = '')
print()
ck_cnt = data[0][1:]
cnts = np.bincount(ck_cnt)
a = len(cnts) - 1
ck_cnt = data[1][1:]
cnts = np.bincount(ck_cnt)
b = len(cnts) - 1
N = len(data) - 2
print('a = ', a)
fout.write('\na = {}\n'.format(a))
print('b = ', b)
fout.write('b = {}\n'.format(b))
print('N = ', N)
fout.write('N = {}\n'.format(N))
Y = np.zeros((N, a, b))
for i in range(N):
for j in range(a):
for
k in range(b):
Y[i][j][k] = float(data[i+2][1 + k + j*b])
print('\nYijk =')
fout.write('\nYijk =\n')
for i in range(N):
for j in range(a):
for
k in range(b):
print('{0:>7.1f}'.format(Y[i][j][k]), end = '')
fout.write('{0:>7.1f}'.format(Y[i][j][k]))
print()
fout.write('\n')
Yppp = Y.mean()
print('\nY... =
{0:<10.2f}'.format(Yppp))
fout.write('\nY... =
{0:<10.2f}\n'.format(Yppp))
Yijp = np.zeros((N,a))
for i in range(N):
for j in range(a):
Yijp[i][j]
= Y[i,j,:].mean()
Ypjk = np.zeros((a, b))
for j in range(a):
for k in range(b):
Ypjk[j][k] = Y[:,j,k].mean()
print('\nY.jk =')
fout.write('\nY.jk =\n')
for j in range(a):
for k in range(b):
print('{0:>10.2f}'.format(Ypjk[j][k]), end = '')
fout.write('{0:>10.2f}'.format(Ypjk[j][k]))
print()
fout.write('\n')
Ypjp = np.zeros(a)
for j in range(a):
Ypjp[j] = Y[:,j,:].mean()
print('Y.j. = ')
fout.write('\nY.j. =\n')
for j in range(a):
print('{0:>10.2f}'.format(Ypjp[j]),
end = '')
fout.write('{0:>10.2f}'.format(Ypjp[j]))
print()
fout.write('\n')
Yppk = np.zeros(b)
for k in range(b):
Yppk[k] = Y[:,:,k].mean()
print('Y..k =')
fout.write('\nY..k =\n')
for k in range(b):
print('{0:>10.2f}'.format(Yppk[k]),
end = '')
fout.write('{0:>10.2f}'.format(Yppk[k]))
print()
fout.write('\n')
SS_btwnS = 0.0
for i in range(N):
for j in range(a):
SS_btwnS += (Yijp[i][j] - Yppp)**2
SS_btwnS *= b
SS_A = b * N * np.sum((Ypjp -
Yppp)**2)
SS_SA = 0.0
for i in range(N):
for j in range(a):
SS_SA += (Yijp[i][j] - Ypjp[j])**2
SS_SA *= b
SS_wthnS = 0.0
for i in range(N):
for j in range(a):
SS_wthnS += np.sum((Y[i,j,:] - Yijp[i][j])**2)
SS_B = a * N * np.sum((Yppk -
Yppp)**2)
SS_AB = 0.0
for j in range(a):
for k in range(b):
SS_AB += (Ypjk[j][k] - Ypjp[j] - Yppk[k] + Yppp)**2
SS_AB *= N
SS_T = np.sum((Y - Yppp)**2)
SS_SB_A = 0.0
for i in range(N):
for j in range(a):
for
k in range(b):
SS_SB_A += (Y[i][j][k] - Yijp[i][j] - Ypjk[j][k] + Ypjp[j])**2
df_btwnS = a * N - 1
df_A = a - 1
df_SA = a * (N - 1)
df_wthnS = a * N * (b - 1)
df_B = b - 1
df_AB = (a - 1) * (b - 1)
df_SB_A = a * (N - 1) * (b - 1)
df_T = a * N * b - 1
MS_A = SS_A / df_A
MS_SA = SS_SA / df_SA
MS_B = SS_B / df_B
MS_AB = SS_AB / df_AB
MS_SB_A = SS_SB_A / df_SB_A
F_A = MS_A / MS_SA
F_B = MS_B / MS_SB_A
F_AB = MS_AB / MS_SB_A
import scipy.stats as ss
p_a = 1.0 - ss.f.cdf(F_A, df_A, df_SA)
p_b = 1.0 - ss.f.cdf(F_B, df_B,
df_SB_A)
p_ab = 1.0 - ss.f.cdf(F_AB, df_AB,
df_SB_A)
print('\n')
print('{0: >10s}{1: >10s}{2:
>10s}{3: >10s}{4: >10s}{5: >10s}'.
format(' ',
'SS', 'df', 'MS', 'F', 'P(>F)'))
fout.write('\n{0: >10s}{1:
>10s}{2: >10s}{3: >10s}{4: >10s}{5: >10s}\n'.
format(' ',
'SS', 'df', 'MS', 'F', 'P(>F)'))
print('{0: >10s}{1: >10.2f}{2:
>10.2f}'.
format('Btwn
Subj', SS_btwnS, df_btwnS))
fout.write('{0: >10s}{1:
>10.2f}{2: >10.2f}\n'.
format('Btwn
Subj', SS_btwnS, df_btwnS))
print('{0: >10s}{1: >10.2f}{2:
>10.2f}{3: >10.2f}{4: >10.2f}{5: >10.5f}'.
format('A',
SS_A, df_A, MS_A, F_A, p_a))
fout.write('{0: >10s}{1:
>10.2f}{2: >10.2f}{3: >10.2f}{4: >10.2f}{5: >10.5f}\n'.
format('A',
SS_A, df_A, MS_A, F_A, p_a))
print('{0: >10s}{1: >10.2f}{2:
>10.2f}{3: >10.2f}'.
format('Res(A)',
SS_SA, df_SA, MS_SA))
fout.write('{0: >10s}{1:
>10.2f}{2: >10.2f}{3: >10.2f}\n'.
format('Res(A)',
SS_SA, df_SA, MS_SA))
print('{0: >10s}{1: >10.2f}{2:
>10.2f}'.
format('Wthn
Subj', SS_wthnS, df_wthnS))
fout.write('{0: >10s}{1:
>10.2f}{2: >10.2f}\n'.
format('Wthn
Subj', SS_wthnS, df_wthnS))
print('{0: >10s}{1: >10.2f}{2:
>10.2f}{3: >10.2f}{4: >10.2f}{5: >10.5f}'.
format('B',
SS_B, df_B, MS_B, F_B, p_b))
fout.write('{0: >10s}{1:
>10.2f}{2: >10.2f}{3: >10.2f}{4: >10.2f}{5: >10.5f}\n'.
format('B',
SS_B, df_B, MS_B, F_B, p_b))
print('{0: >10s}{1: >10.2f}{2:
>10.2f}{3: >10.2f}{4: >10.2f}{5: >10.5f}'.
format('AB',
SS_AB, df_AB, MS_AB, F_AB, p_ab))
fout.write('{0: >10s}{1:
>10.2f}{2: >10.2f}{3: >10.2f}{4: >10.2f}{5: >10.5f}\n'.
format('AB',
SS_AB, df_AB, MS_AB, F_AB, p_ab))
print('{0: >10s}{1: >10.2f}{2:
>10.2f}{3: >10.2f}'.
format('Res(B)',
SS_SB_A, df_SB_A, MS_SB_A))
fout.write('{0: >10s}{1:
>10.2f}{2: >10.2f}{3: >10.2f}\n'.
format('Res(B)',
SS_SB_A, df_SB_A, MS_SB_A))
print('{0: >10s}{1: >10.2f}{2:
>10.2f}'.
format('Total',
SS_T, df_T))
fout.write('{0: >10s}{1:
>10.2f}{2: >10.2f}\n'.
format('Total',
SS_T, df_T))
fout.close()
print('\nOutput file {} was
saved.\n'.format(fout_name))
上のスクリプトを実行すると、入力データファイル名と出力用テキストファイル名の設定が求められる。
図１
データファイルはCSVファイルの形式で用意する。Excelで用意するときは図２に示すように用意して、保存ファイルのファイル拡張子として「.csv」を選べばよい。
図２
第１行目は被験者間要因Aの水準を示すものである。第１セルにAと書いて、要因Aであることを表している。第２セルからAの要因を１から水準の数a（図２ではa=3）までの整数値で表す。整数値は１からaまでの昇順である。第２行目は被験者内要因Bの水準を示すものである。第１セルにBと書いて、要因Bであることを表している。要因Aの各水準内において要因Bの水準が１から要因Bの水準の数ｂ（図２では、b=2）までの整数値が昇順に並べられている。第３行目から、各水準の組み合わせの下でのデータを設定する。図２では、例えば赤枠で囲まれたデータは、水準A1の下で水準B1とB2を繰り返し要因の水準として課された一人のデータ「５と６」、あるいは水準A3の下で繰り返し要因の水準B1とB2が課された一人のデータ「２と６」が設定されている。なお、CSV形式のファイルについては、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan, PyMC）――」丸善出版＞のpp. 33-35で説明している。拙著では、第２章においてグラフ描画の説明も行っている。
出力用テキストファイル名の設定後、計算が始まり、計算が終了すると下記のようなコンソール上への出力で終了する。
SS df
MS
F P(>F)
Btwn Subj 63.78 8.00
A
53.44
2.00
26.72
15.52 0.00425
Res(A) 10.33 6.00 1.72
Wthn Subj 22.50 9.00
B
16.06
1.00
16.06
22.23 0.00328
AB
2.11
2.00
1.06
1.46 0.30403
Res(B) 4.33 6.00 0.72
Total 86.28 17.00
Output file Results.txt was saved.
>>>
出力ファイル（上の例では、Results.txt）を開くと、図２のデータの場合、リスト２の内容になっている。
リスト２ 出力ファイル例
Data file name = Data.csv
a = 3
b = 2
N = 3
Yijk =
5.0 6.0 8.0 9.0 3.0 5.0
4.0 8.0 6.0 7.0 1.0 3.0
6.0 7.0 7.0 8.0 2.0 6.0
Y... = 5.61
Y.jk =
5.00 7.00 7.00 8.00 2.00 4.67
Y.j. =
6.00 7.50 3.33
Y..k =
4.67 6.56
SS df
MS
F P(>F)
Btwn Subj 63.78 8.00
A
53.44
2.00
26.72
15.52 0.00425
Res(A) 10.33 6.00 1.72
Wthn Subj 22.50 9.00
B 16.06 1.00 16.06 22.23 0.00328
AB
2.11
2.00
1.06
1.46 0.30403
Res(B) 4.33 6.00 0.72
Total 86.28 17.00
入力データファイルから読み込んだデータが出力されているので、データが正しく読み込まれたかどうかの確認ができる。
続いて、平均値
、
、
、
が出力されている。
最後に、分散分析表が出力されている。P(>F)はp値である。