Up

被験者内２要因分散分析

Pythonスクリプトとサンプルデータファイルは圧縮ファイルanova_two_rel.zipにまとめた。Jupyter Notebook用は圧縮ファイルanova_two_rel_jn.zipとした。これらはダウンロード・解凍して自由に利用されたい。

　被験者内２要因分散分析のモデルを以下のように設定する。

被験者内２要因、ＡおよびＢ、のモデルにおける平方和を次のように分解する。ただし、は、番目の被験者に対する要因Ａの水準、要因Ｂの水準におけるデータであり、要因Ａの水準数は、要因Ｂの水準数は、被験者数はであるとする。

ここで、

である。

　上記の各平方和の平均の期待値は次のようになる。

　以上より、分散分析表を表１のように作成する。

表１　分散分析表

変動因	平方和	自由度	平均平方和	Ｆ比
被験者間
被験者内
Ａ
残差（Ａ）
Ｂ
残差（Ｂ）
Ａ×Ｂ
残差（Ａ×Ｂ）
全体

上のモデルに基づいて、リスト１に示すスクリプトを用意した。このスクリプトのファイルとサンプルデータファイルはまとめて圧縮ファイルfiles.zipとした。

リスト１　被験者内２要因分散分析Pythonスクリプト例

import csv

import numpy as np

import scipy.stats as ss

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')

Yipp = np.zeros(N)

for i in range(N):

Yipp[i] = Y[i,:,:].mean()

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')

Yipk = np.zeros((N,b))

for i in range(N):

for k in range(b):

Yipk[i][k] = Y[i,:,k].mean()

SS_T = np.sum((Y - Yppp)**2)

SS_S = a * b * np.sum((Yipp - Yppp)**2)

SS_WS = SS_T - SS_S

SS_A = b * N * np.sum((Ypjp - Yppp)**2)

SS_AS = 0.0

for i in range(N):

for j in range(a):

SS_AS += (Yijp[i][j] - Yipp[i] - Ypjp[j] + Yppp)**2

SS_AS *= b

SS_B = a * N * np.sum((Yppk - Yppp)**2)

SS_BS = 0.0

for i in range(N):

for k in range(b):

SS_BS += (Yipk[i][k] - Yipp[i] - Yppk[k] + Yppp)**2

SS_BS *= a

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_ABS = SS_T - (SS_S + SS_A + SS_AS + SS_B + SS_BS + SS_AB)

df_T = N * a * b - 1

df_S = N - 1

df_WS = N * (a * b - 1)

df_A = a - 1

df_AS = (N - 1) * (a - 1)

df_B = b - 1

df_BS = (N - 1) * (b - 1)

df_AB = (a - 1) * (b - 1)

df_ABS = (N - 1) * (a - 1) * (b - 1)

MS_S = SS_S / df_S

MS_A = SS_A / df_A

MS_AS = SS_AS / df_AS

MS_B = SS_B / df_B

MS_BS = SS_BS / df_BS

MS_AB = SS_AB / df_AB

MS_ABS = SS_ABS / df_ABS

F_A = MS_A / MS_AS

F_B = MS_B / MS_BS

F_AB = MS_AB / MS_ABS

p_a = 1.0 - ss.f.cdf(F_A, df_A, df_AS)

p_b = 1.0 - ss.f.cdf(F_B, df_B, df_BS)

p_ab = 1.0 - ss.f.cdf(F_AB, df_AB, df_ABS)

print('\n')

fout.write('\n\n')

print('{0: >10s}{1: >10s}{2: >10s}{3: >10s}{4: >10s}{5: >10s}'.

format('SV', 'SS', 'df', 'MS', 'F', 'P(>F)'))

fout.write('{0: >10s}{1: >10s}{2: >10s}{3: >10s}{4: >10s}{5: >10s}\n'.

format('SV', 'SS', 'df', 'MS', 'F', 'P(>F)'))

print('{0: >10s}{1: >10.2f}{2: >10.2f}{3: >10.2f}'.

format('S', SS_S, df_S, MS_S))

fout.write('{0: >10s}{1: >10.2f}{2: >10.2f}{3: >10.2f}\n'.

format('S', SS_S, df_S, MS_S))

print('{0: >10s}{1: >10.2f}{2: >10.2f}'.

format('SSw.s', SS_WS, df_WS))

fout.write('{0: >10s}{1: >10.2f}{2: >10.2f}\n'.

format('SSw.s', SS_WS, df_WS))

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('AS', SS_AS, df_AS, MS_AS))

fout.write('{0: >10s}{1: >10.2f}{2: >10.2f}{3: >10.2f}\n'.

format('AS', SS_AS, df_AS, MS_AS))

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}'.

format('BS', SS_BS, df_BS, MS_BS))

fout.write('{0: >10s}{1: >10.2f}{2: >10.2f}{3: >10.2f}\n'.

format('BS', SS_BS, df_BS, MS_BS))

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('ABS', SS_ABS, df_ABS, MS_ABS))

fout.write('{0: >10s}{1: >10.2f}{2: >10.2f}{3: >10.2f}\n'.

format('ABS', SS_ABS, df_ABS, MS_ABS))

print('{0: >10s}{1: >10.2f}{2: >10.2f}'.

format('SStotal', SS_T, df_T))

fout.write('{0: >10s}{1: >10.2f}{2: >10.2f}\n'.

format('SStotal', 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）までの整数値が昇順に並べられている。第３行目から、各水準の組み合わせの下でのデータを設定する。図２では、例えば赤枠で囲まれたデータは、要因Aの各水準A1、A2、A3と要因Bの各水準B1、B2の組み合わせの下での一人の被験者のデータである。要因Aと要因Bは被験者内要因（繰り返し要因）なので、各被験者は水準Aと水準Bの全ての組み合わせ条件のデータがある。赤線で囲まれたデータ「37、・・・、59」はある一人の被験者のデータである。なお、CSV形式のファイルについては、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan, PyMC）――」丸善出版＞のpp. 33-35で説明している。拙著では、第２章においてグラフ描画の説明も行っている。

出力用テキストファイル名の設定後、計算が始まり、計算が終了すると下記のようなコンソール上への出力で終了する。

SV SS df MS F P(>F)

S 49.14 5.00 9.83

SSw.s 1515.17 30.00

A 770.39 2.00 385.19 30.50 0.00006

AS 126.28 10.00 12.63

B 476.69 1.00 476.69 24.88 0.00415

BS 95.81 5.00 19.16

AB 0.06 2.00 0.03 0.01 0.99398

ABS 45.94 10.00 4.59

SStotal 1564.31 35.00

Output file Results.txt was saved.

>>>

出力ファイル（上の例では、Results.txt）を開くと、図２のデータの場合、リスト２の内容になっている。

リスト２　出力ファイル例

Data file name = Data.csv

a = 3

b = 2

N = 6

Yijk =

37.0 41.0 46.0 53.0 49.0 59.0

42.0 47.0 48.0 56.0 55.0 59.0

45.0 46.0 48.0 53.0 53.0 53.0

37.0 50.0 47.0 55.0 48.0 57.0

40.0 56.0 46.0 55.0 47.0 59.0

38.0 42.0 49.0 56.0 53.0 62.0

Y... = 49.64

Y.jk =

39.83 47.00 47.33 54.67 50.83 58.17

Y.j. =

43.42 51.00 54.50

Y..k =

46.00 53.28

SV SS df MS F P(>F)

S 49.14 5.00 9.83

SSw.s 1515.17 30.00

A 770.39 2.00 385.19 30.50 0.00006

AS 126.28 10.00 12.63

B 476.69 1.00 476.69 24.88 0.00415

BS 95.81 5.00 19.16

AB 0.06 2.00 0.03 0.01 0.99398

ABS 45.94 10.00 4.59

SStotal 1564.31 35.00

入力データファイルから読み込んだデータが出力されているので、データが正しく読み込まれたかどうかの確認ができる。

続いて、平均値、、、が出力されている。

最後に、分散分析表が出力されている。P(>F)はp値である。

Up