Up

被験者間２要因分散分析

被験者間２要因の分散分析を行うPythonスクリプトを、＜岡本安晴「データを使いこなすための統計入門」22世紀アート、pp. 115-128＞に説明されている式に基づいて作成した。スクリプトとサンプルデータは圧縮ファイルanova_two_ind.zipとしてまとめた。Jupyter Notebook用のものは圧縮ファイルanova_two_ind_jn.zipとした。それぞれダウウンロード・解凍して自由に利用されたい。

被験者間２要因分散分析の基本式は、全分散を、以下のように平方和に分解するものである。

ここで、

これより、平均平方和は次式で与えられる。

自由度dfは次式で与えられる。

このとき、F値とF分布のパラメータは次の通りである。

主効果A：

主効果B：

交互作用AB：

上の式に基づいて、リスト１に示すスクリプトを用意した。このスクリプトのファイルとサンプルデータファイルはまとめて圧縮ファイル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...')

fout.write('\nData...\n')

for v in data:

for w in v:

print('{0: >7s}'.format(w), end = '')

fout.write('{0: >7s}'.format(w))

print()

fout.write('\n')

ck_cnt = data[0][1:]

cnts = np.bincount(ck_cnt)

print(cnts)

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)

print('b = ', b)

print('N = ', 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])

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

print()

Yppp = Y.mean()

print('Y... = ', Yppp)

fout.write('\nY... = {}\n'.format(Yppp))

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_A = b * N * np.sum((Ypjp - Yppp)**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_S_AB = SS_T - SS_A - SS_B - SS_AB

dfT = a * b * N

dfA = a - 1

dfB = b - 1

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

dfS_AB = a * b * (N-1)

MS_A = SS_A / dfA

MS_B = SS_B / dfB

MS_AB = SS_AB / dfAB

MS_S_AB = SS_S_AB / dfS_AB

F_A = MS_A / MS_S_AB

F_B = MS_B / MS_S_AB

F_AB = MS_AB / MS_S_AB

import scipy.stats as ss

p_a = 1.0 - ss.f.cdf(F_A, dfA, dfS_AB)

p_b = 1.0 - ss.f.cdf(F_B, dfB, dfS_AB)

p_ab = 1.0 - ss.f.cdf(F_AB, dfAB, dfS_AB)

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}{3: >10.2f}{4: >10.2f}{5: >10.5f}'.

format('A', SS_A, dfA, 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, dfA, MS_A, F_A, p_a))

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

format('B', SS_B, dfB, 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, dfB, 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, dfAB, 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, dfAB, MS_AB, F_AB, p_ab))

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

format('S/AB', SS_S_AB, dfS_AB, MS_S_AB))

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

format('S/AB', SS_S_AB, dfS_AB, MS_S_AB))

fout.close()

print('\nOutput file {} was saved.\n'.format(fout_name))

上のスクリプトを実行すると、入力データファイル名と出力用テキストファイル名の設定が求められる。

図１

データファイルはCSVファイルの形式で用意する。Excelで用意するときは図２に示すように用意して、保存ファイルのファイル拡張子として「.csv」を選べばよい。

図２

第１行目は要因Aの水準を示すものである。第１セルにAと書いて、要因Aであることを表している。第２セルからAの要因を１から水準の数a（図２ではa=2）までの整数値で表す。整数値は１からaまでの昇順である。第２行目は要因Bの水準を示すものである。第１セルにBと書いて、要因Bであることを表している。要因Aの各水準内において要因Bの水準が１から要因Bの水準の数ｂ（図２では、b=2）までの整数値が昇順に並べられている。第３行目から、各水準の組み合わせの下でのデータを設定する。図２では、要因A1,B1のデータとして、63、69、72が設定されおり、要因A2,B2のデータとして66、71、73が設定されている。なお、CSV形式のファイルについては、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan, PyMC）――」丸善出版＞のpp. 33-35で説明している。拙著では、第２章においてグラフ描画の説明も行っている。

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

SS df MS F P(>F)

A 1452.00 1.00 1452.00 72.60 0.00003

B 1200.00 1.00 1200.00 60.00 0.00006

AB 12.00 1.00 12.00 0.60 0.46086

S/AB 160.00 8.00 20.00

Output file Results.txt was saved.

>>>

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

リスト２　出力ファイル例

Data file name = Data.csv

Data...

A 1 1 2 2

B 1 2 1 2

63 46 87 66

69 49 92 71

72 55 97 73

Y... = 70.0

Y.jk =

68.00 50.00 92.00 70.00

Y.j. =

59.00 81.00

Y..k =

80.00 60.00

SS df MS F P(>F)

A 1452.00 1.00 1452.00 72.60 0.00003

B 1200.00 1.00 1200.00 60.00 0.00006

AB 12.00 1.00 12.00 0.60 0.46086

S/AB 160.00 8.00 20.00

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

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

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

Up