Up

被験者間1要因・被験者内1要因分散分析

 

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

 

被験者間要因(非繰り返し要因)A、および被験者内要因(繰り返し要因)Bのモデルを次のようにおく。

ここで、主効果A

主効果B

交互作用AB

ij番目の被験者固有の効果

jとBkとの交互作用による効果を調整した、ij番目の被験者と要因Bの水準kとの交互作用

である。

 データ値の平方和を次のように分解する。

ここで

ただし、

である。

 これらの平方和の期待値は以下のようになる。

ここで、

である。

 以上より、分散分析表は表1のように作成される。

 

表1 分散分析表

変動要因

平方和

自由度

平均平方和

F比

被験者間

 

 

残差(A)

 

被験者内

 

 

A×B

残差(B)

 

全体

 

 

 

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

 

リスト1 被験者間1要因被験者内1要因分散分析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))

 

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

図1

 

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

図2

 

第1行目は被験者間要因Aの水準を示すものである。第1セルにAと書いて、要因Aであることを表している。第2セルからAの要因を1から水準の数a(図2ではa=3)までの整数値で表す。整数値は1からaまでの昇順である。第2行目は被験者内要因Bの水準を示すものである。第1セルにBと書いて、要因Bであることを表している。要因Aの各水準内において要因Bの水準が1から要因Bの水準の数b(図2では、b=2)までの整数値が昇順に並べられている。第3行目から、各水準の組み合わせの下でのデータを設定する。図2では、例えば赤枠で囲まれたデータは、水準A1の下で水準B1B2を繰り返し要因の水準として課された一人のデータ「5と6」、あるいは水準A3の下で繰り返し要因の水準B1B2が課された一人のデータ「2と6」が設定されている。なお、CSV形式のファイルについては、拙著<岡本安晴「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析(PyStan, PyMC)――」丸善出版>のpp. 33-35で説明している。拙著では、第2章においてグラフ描画の説明も行っている。

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

 

                  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)を開くと、図2のデータの場合、リスト2の内容になっている。

 

リスト2 出力ファイル例

 

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値である。

 

Up