Up

２要因分散分析

非繰り返し、固定効果

非繰り返し２要因固定効果の構造方程式（structural model）は次式で与えられる（Winer et al., 1991, p. 292）。

ここで、

である（Kirk, 1995, p. 372）。

式（１）におけるパラメータに対して次式のようにおく。

このとき、モデル（１）は次式のように書ける。

式（４）に対してStanスクリプトをリスト１のように作成した。

リスト１　繰り返し１要因固定効果モデルのStanスクリプト（TwoFctrFixed.stan）

data {

int a;

int b;

int N;

real X[N, a, b];

}

parameters {

real<lower = 0.0> sgm;

real G[a, b];

}

model {

sgm ~ exponential(0.001);

for (j in 1:a)

for (k in 1:b) {

G[j][k] ~ normal(0.0, 1000.0);

}

for (i in 1:N)

for (j in 1:a)

for (k in 1:b) {

X[i][j][k] ~ normal(G[j][k], sgm);

}

}

パラメータに対して、制約条件（２）により、主効果および交互作用のパラメータは以下のように求めることができる。

リスト１のStanスクリプトによってを求め、式（５）〜（８）によっ主効果と、交互作用を求めるPythonスクリプトを本ウェブサイトの後半に示すリスト２のように用意した。スクリプトファイルなどは、圧縮ファイルfilesTwoFctrsFixed.zipとしてまとめた。ダウンロード展開すれば、自由に利用できる。PythonにおけるStanの使い方については＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞第10章〜第12章で説明している。

リスト２では、上の式（５）〜（８）による主効果と交互作用の計算を以下のスクリプトによって行っている。

for vG in fit['G']:

Gdd = vG.mean() # G..

Gjd = vG.mean(axis = 1) # Gj.

Gdk = vG.mean(axis = 0) # G.k

f_a = Gjd - Gdd

f_b = Gdk - Gdd

f_ab = vG - Gjd.reshape((a,1)) - Gdk + Gdd

リスト２を実行すると、まず入力データファイル名の設定が求められる（図１）。

図１

入力データファイルは、図２に示す形式で用意する。

図２　データ例（data2fctrs.xlsx）

図２は、要因Factor-AがA1、A2、A3の３水準、要因Factor-BがB1、B2の２水準で、各条件におけるデータ数Nが10であるデータを表している。第１行目には要因Aの水準を書き、第２行目に要因Bの水準を書く。要因Bの水準は要因Aの水準の下に書き、水準Bの順序は、要因Aの各水準の下で同じ順序とする。第１列目には、各要因の名前と水準数、およびデータの識別値を書く。第１列目の１行目には、要因Aの名前の後に半角スラッシュ/を置いて、続けて要因Aの水準数を書く。第１列目の2行目には、要因Bの名前の後に半角スラッシュ/を置いて、続けて要因Bの水準数を書く。データ値は第３行目以降に、データの識別値に続けて該当する要因Aの水準と要因Bの水準の組み合わせの下に書く。なお、データの識別値は、図２では文字列であるが、これは通し番号でもよい。また、同じ識別値が複数回用いられてもよい。

入力データファイル名を設定して（図１では、図２のファイル名data2fctrs.xlsxが設定されている）、「Enter」キーを押すと、次に計算結果の出力テキストファイル名の設定が求められる（図１）。出力ファイル名はファイル拡張子が「.txt」である任意のファイル名を設定する。

出力ファイル名の設定後、「Enter」キーを押すと、グラフの画像保存ファイル名の設定が求められる（図１）。ファイル拡張子が「.png」である任意のファイル名を設定する。設定したファイル名の「.png」を「_a.png」などに変更した名前でファイルが保存される。

保存画像ファイル名を設定して「Enter」キーを押すと、データが読み込まれ、計算が始まる。

MCMCサンプリングが終了すると、要因Aの主効果の事後分布が表示される（図３）。

図３

図３では、３つの水準の主効果の事後分布が明瞭に分離されている。

図３のウィンドウを閉じると、次に要因Bの主効果の事後分布が表示さえっる（図４）。

図４

要因Bの主効果も明瞭に分離されている。

図４のウィンドウを閉じると、交互作用の事後分布のグラフが表示される（図５）。

図５

交互作用A1*B2とA3*B1の事後分布は０より右側にあり、交互作用A1*B1とA3*B2の事後分布は０より左側にある。交互作用A2*B1とA2*B2の事後分布は０を中ほどに含む形で分布している。

図５のウィンドウを閉じると、効果量の事後分布のグラフが表示される（図６）。

図６

効果量は、主効果および交互作用に対して次式で与えられている（cf. Kirk, 1995, p. 181）。

図６のウィンドウを閉じると、スクリプトの実行終了となる（図７）。

図７

保存されたファイル名が表示されている。これらのファイルは、スクリプトの実行終了後、開くことができる。

出力テキストファイルをPython IDLEで開こうとすると、図８のメッセージが表示されることがある。

図８

「s-jis」と設定して、「OK」ボタンをクリックすると開くことができる（図９）。

図９

出力ファイルの内容は以下のようになっている。

入力データファイル: data2fctrs.xlsx

NameA = Factor-A

a = 3

LevelsA: ['A1', 'A2', 'A3']

NameB = Factor-B

b = 2

LevelB: ['B1' 'B2']

N = 10

X =

57.0 38.0 83.0 53.0 62.0 33.0

71.0 61.0 77.0 80.0 69.0 37.0

72.0 49.0 85.0 61.0 69.0 31.0

68.0 44.0 72.0 60.0 66.0 37.0

63.0 54.0 82.0 63.0 57.0 43.0

78.0 61.0 83.0 53.0 57.0 31.0

62.0 60.0 78.0 62.0 63.0 25.0

63.0 53.0 95.0 61.0 71.0 31.0

63.0 64.0 81.0 68.0 66.0 35.0

63.0 53.0 71.0 57.0 63.0 28.0

Inference for Stan model: anon_model_46f071e5b1f173f479a4168def849c4c.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

sgm 6.79 9.6e-3 0.68 5.61 6.3 6.73 7.22 8.26 4963 1.0

G[1,1] 66.04 0.03 2.19 61.71 64.59 66.07 67.52 70.2 6090 1.0

G[2,1] 80.72 0.03 2.15 76.45 79.29 80.72 82.14 84.88 6049 1.0

G[3,1] 64.26 0.03 2.13 60.06 62.83 64.28 65.66 68.51 6493 1.0

G[1,2] 53.72 0.03 2.17 49.41 52.23 53.73 55.22 57.93 6619 1.0

G[2,2] 61.82 0.03 2.2 57.51 60.36 61.81 63.31 66.12 5435 1.0

G[3,2] 33.1 0.03 2.11 28.99 31.72 33.12 34.48 37.23 6208 1.0

lp__ -142.3 0.05 2.0 -147.3 -143.4 -142.0 -140.9 -139.5 1482 1.0

Samples were drawn using NUTS at Sun Aug 9 15:00:45 2020.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

主効果：Factor-A

A1：

平均値 = -0.065

中央値 = -0.052

Q1 = -0.863 Q3 = 0.759

95% CI = [-2.527, 2.311]

A2：

平均値 = 11.327

中央値 = 11.338

Q1 = 10.485 Q3 = 12.167

95% CI = [8.875, 13.704]

A3：

平均値 = -11.261

中央値 = -11.257

Q1 = -12.077 Q3 = -10.464

95% CI = [-13.796, -8.782]

主効果：Factor-B

B1：

平均値 = 10.396

中央値 = 10.390

Q1 = 9.815 Q3 = 10.987

95% CI = [8.705, 12.110]

B2：

平均値 = -10.396

中央値 = -10.390

Q1 = -10.987 Q3 = -9.815

95% CI = [-12.110, -8.705]

交互作用：

A1・B1：

平均値 = -4.236

中央値 = -4.221

Q1 = -5.062 Q3 = -3.398

95% CI = [-6.779, -1.735]

A1・B2：

平均値 = 4.236

中央値 = 4.221

Q1 = 3.398 Q3 = 5.062

95% CI = [1.735, 6.779]

A2・B1：

平均値 = -0.948

中央値 = -0.948

Q1 = -1.786 Q3 = -0.128

95% CI = [-3.379, 1.517]

A2・B2：

平均値 = 0.948

中央値 = 0.948

Q1 = 0.128 Q3 = 1.786

95% CI = [-1.517, 3.379]

A3・B1：

平均値 = 5.184

中央値 = 5.202

Q1 = 4.353 Q3 = 6.031

95% CI = [2.700, 7.621]

A3・B2：

平均値 = -5.184

中央値 = -5.202

Q1 = -6.031 Q3 = -4.353

95% CI = [-7.621, -2.700]

効果量（主効果）：Factor-A

平均値 = 1.378

中央値 = 1.377

Q1 = 1.248 Q3 = 1.509

95% CI = [1.021, 1.741]

効果量（主効果）：Factor-B

平均値 = 1.546

中央値 = 1.545

Q1 = 1.409 Q3 = 1.681

95% CI = [1.168, 1.942]

効果量（交互作用）：Factor-A・Factor-B

平均値 = 0.595

中央値 = 0.594

Q1 = 0.504 Q3 = 0.687

95% CI = [0.315, 0.873]

Gjk's =

Gjk(A1,B1):

平均値 = 66.038

中央値 = 66.071

Q1 = 64.593 Q3 = 67.519

95% CI = [61.721, 70.186]

Gjk(A1,B2):

平均値 = 53.717

中央値 = 53.734

Q1 = 52.229 Q3 = 55.216

95% CI = [49.425, 57.921]

Gjk(A2,B1):

平均値 = 80.717

中央値 = 80.724

Q1 = 79.290 Q3 = 82.138

95% CI = [76.454, 84.869]

Gjk(A2,B2):

平均値 = 61.822

中央値 = 61.811

Q1 = 60.357 Q3 = 63.310

95% CI = [57.521, 66.117]

Gjk(A3,B1):

平均値 = 64.263

中央値 = 64.277

Q1 = 62.830 Q3 = 65.662

95% CI = [60.062, 68.505]

Gjk(A3,B2):

平均値 = 33.101

中央値 = 33.119

Q1 = 31.716 Q3 = 34.480

95% CI = [28.996, 37.225]

引用文献

Winer, B. J., Brown, D.R., & Michels, K.M. (1991). Statistical principles in experimental design, third ed. McGraw-Hill, Inc.

Kirk, R.E. (1995). Experimental design: Procedures for the behavioral sciences, third ed. Brooks/Cole Publishing Company.

リスト２　非繰り返し２要因固定効果モデルのPythonスクリプト例

#

# Yasuharu Okamoto, 2020.08

#

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import pystan

import seaborn as sb

finnm = input('入力データファイル（*.xlsx）= ')

txtoutnm = input('出力ファイル（*.txt）= ')

foutnm = input('画像保存ファイル（*.png）= ')

txtout = open(txtoutnm, 'w')

#

# Excelファイルの読み込み

#　　　モジュール xlrd がインストール済みであること

#

xlsx_fl = pd.ExcelFile(finnm) # ExcelファイルからExcelFile型オブジェクトを生成

data_xlsx = pd.read_excel(xlsx_fl, header = None) # ExcelFile型オブジェクトからDataFrame型オブジェクトを生成

data = data_xlsx.values # Dataframe型オブジェクトからnumpyの配列オブジェクトを生成

print('\n', data)

txtout.write('\n入力データファイル: {}\n'.format(finnm))

NameA = data[0][0].split(sep = '/')[0]

a = int(data[0][0].split(sep = '/')[1])

NameB = data[1][0].split(sep = '/')[0]

b = int(data[1][0].split(sep = '/')[1])

LevelsA = [data[0][1 + j*b] for j in range(a)]

LevelsB = data[1][1:b+1]

print('\nNameA = ', NameA)

print('a = ', a)

print('LevelsA: ', LevelsA)

print('NameB = ', NameB)

print('b = ', b)

print('LevelB: ', LevelsB)

txtout.write('\nNameA = {}\n'.format(NameA))

txtout.write('a = {}\n'.format(a))

txtout.write('LevelsA: {}\n'.format(LevelsA))

txtout.write('\nNameB = {}\n'.format(NameB))

txtout.write('b = {}\n'.format(b))

txtout.write('LevelB: {}\n'.format(LevelsB))

N = len(data) - 2

print('N = ', N)

txtout.write('\nN = {}\n'.format(N))

X = np.empty((N, a, b))

for i in range(N):

for j in range(a):

for k in range(b):

X[i][j][k] = data[i + 2][1 + j * b + k]

print('\nX =')

txtout.write('\nX =\n')

for i in range(N):

for j in range(a):

for k in range(b):

print(' ', X[i][j][k], end = '')

txtout.write(' {}'.format(X[i][j][k]))

print()

txtout.write('\n')

StanData = {'a': a, 'b': b, 'N': N, 'X': X}

sm = pystan.StanModel(file = 'TwoFctrFixed.stan')

fit = sm.sampling(data = StanData, n_jobs = 1)

print(fit)

txtout.write('\n{}\n'.format(fit))

N_Samples = len(fit['G'])

#

# 主効果（f_a, f_b）と交互作用（f_ab）

#

f_as = np.empty((N_Samples, a))

f_bs = np.empty((N_Samples, b))

f_abs = np.empty((N_Samples, a, b))

n_s = 0

for vG in fit['G']:

Gdd = vG.mean() # G..

Gjd = vG.mean(axis = 1) # Gj.

Gdk = vG.mean(axis = 0) # G.k

f_a = Gjd - Gdd

f_b = Gdk - Gdd

f_ab = vG - Gjd.reshape((a,1)) - Gdk + Gdd

f_as[n_s] = f_a

f_bs[n_s] = f_b

f_abs[n_s] = f_ab

n_s += 1

# 日本語表記の設定

plt.rcParams['font.family'] = ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']

print('\n主効果：{}'.format(NameA))

txtout.write('\n主効果：{}\n'.format(NameA))

for j in range(a):

print('\n{}：'.format(LevelsA[j]))

txtout.write('\n{}：\n'.format(LevelsA[j]))

v_mean = f_as[:, j].mean()

v_025p, v_25p, v_med, v_75p, v_975p = \

np.percentile(f_as[:, j], [2.5, 25, 50, 75, 97.5])

print(' 平均値 = {0:.3f}'.format(v_mean))

print(' 中央値 = {0:.3f}'.format(v_med))

print(' Q1 = {0:.3f} Q3 = {1:.3f}'.format(v_25p, v_75p))

print(' 95% CI = [{0:.3f}, {1:.3f}]'.format(v_025p, v_975p))

txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))

txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))

txtout.write(' Q1 = {0:.3f} Q3 = {1:.3f}\n'.format(v_25p, v_75p))

txtout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(v_025p, v_975p))

plt.title('主効果：{}'.format(NameA), fontsize = 20)

for j in range(a):

sb.kdeplot(f_as[:, j], label = LevelsA[j])

plt.savefig(foutnm[:-4] + '_a.png')

plt.yticks([])

SavedFiles = [foutnm[:-4] + '_a.png']

plt.show()

print('\n主効果：{}'.format(NameB))

txtout.write('\n主効果：{}\n'.format(NameB))

for k in range(b):

print('\n{}：'.format(LevelsB[k]))

txtout.write('\n{}：\n'.format(LevelsB[k]))

v_mean = f_bs[:, k].mean()

v_025p, v_25p, v_med, v_75p, v_975p = \

np.percentile(f_bs[:, k], [2.5, 25, 50, 75, 97.5])

print(' 平均値 = {0:.3f}'.format(v_mean))

print(' 中央値 = {0:.3f}'.format(v_med))

print(' Q1 = {0:.3f} Q3 = {1:.3f}'.format(v_25p, v_75p))

print(' 95% CI = [{0:.3f}, {1:.3f}]'.format(v_025p, v_975p))

txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))

txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))

txtout.write(' Q1 = {0:.3f} Q3 = {1:.3f}\n'.format(v_25p, v_75p))

txtout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(v_025p, v_975p))

plt.title('主効果：{}'.format(NameB), fontsize = 20)

for k in range(b):

sb.kdeplot(f_bs[:, k], label = LevelsB[k])

plt.yticks([])

plt.savefig(foutnm[:-4] + '_b.png')

SavedFiles.append(foutnm[:-4] + '_b.png')

plt.show()

print('\n交互作用：')

txtout.write('\n交互作用：\n')

for j in range(a):

for k in range(b):

print('\n{0}・{1}：'.format(LevelsA[j], LevelsB[k]))

txtout.write('\n{0}・{1}：\n'.format(LevelsA[j], LevelsB[k]))

v_mean = f_abs[:, j, k].mean()

v_025p, v_25p, v_med, v_75p, v_975p = \

np.percentile(f_abs[:, j, k], [2.5, 25, 50, 75, 97.5])

print(' 平均値 = {0:.3f}'.format(v_mean))

print(' 中央値 = {0:.3f}'.format(v_med))

print(' Q1 = {0:.3f} Q3 = {1:.3f}'.format(v_25p, v_75p))

print(' 95% CI = [{0:.3f}, {1:.3f}]'.format(v_025p, v_975p))

txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))

txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))

txtout.write(' Q1 = {0:.3f} Q3 = {1:.3f}\n'.format(v_25p, v_75p))

txtout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(v_025p, v_975p))

plt.title('交互作用：{0}・{1}'.format(NameA, NameB), fontsize = 20)

for j in range(a):

for k in range(b):

sb.kdeplot(f_abs[:, j, k], label = LevelsA[j] + '*' + LevelsB[k])

plt.yticks([])

plt.savefig(foutnm[:-4] + '_ab.png')

SavedFiles.append(foutnm[:-4] + '_ab.png')

plt.show()

es_a = []

es_b = []

es_ab = []

for a_s, b_s, ab_s, sgm_s in zip(f_as, f_bs, f_abs, fit['sgm']):

es_a.append(np.std(a_s) / sgm_s)

es_b.append(np.std(b_s) / sgm_s)

es_ab.append(np.std(ab_s) / sgm_s)

es_a = np.array(es_a)

es_b = np.array(es_b)

es_ab = np.array(es_ab)

print('\n効果量（主効果）：{}'.format(NameA))

txtout.write('\n\n効果量（主効果）：{}\n'.format(NameA))

v_mean = es_a.mean()

v_025p, v_25p, v_med, v_75p, v_975p = \

np.percentile(es_a, [2.5, 25, 50, 75, 97.5])

print(' 平均値 = {0:.3f}'.format(v_mean))

print(' 中央値 = {0:.3f}'.format(v_med))

print(' Q1 = {0:.3f} Q3 = {1:.3f}'.format(v_25p, v_75p))

print(' 95% CI = [{0:.3f}, {1:.3f}]'.format(v_025p, v_975p))

txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))

txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))

txtout.write(' Q1 = {0:.3f} Q3 = {1:.3f}\n'.format(v_25p, v_75p))

txtout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(v_025p, v_975p))

print('\n効果量（主効果）：{}'.format(NameB))

txtout.write('\n効果量（主効果）：{}\n'.format(NameB))

v_mean = es_b.mean()

v_025p, v_25p, v_med, v_75p, v_975p = \

np.percentile(es_b, [2.5, 25, 50, 75, 97.5])

print(' 平均値 = {0:.3f}'.format(v_mean))

print(' 中央値 = {0:.3f}'.format(v_med))

print(' Q1 = {0:.3f} Q3 = {1:.3f}'.format(v_25p, v_75p))

print(' 95% CI = [{0:.3f}, {1:.3f}]'.format(v_025p, v_975p))

txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))

txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))

txtout.write(' Q1 = {0:.3f} Q3 = {1:.3f}\n'.format(v_25p, v_75p))

txtout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(v_025p, v_975p))

print('\n効果量（交互作用）：{0}・{1}'.format(NameA, NameB))

txtout.write('\n効果量（交互作用）：{0}・{1}\n'.format(NameA, NameB))

v_mean = es_ab.mean()

v_025p, v_25p, v_med, v_75p, v_975p = \

np.percentile(es_ab, [2.5, 25, 50, 75, 97.5])

print(' 平均値 = {0:.3f}'.format(v_mean))

print(' 中央値 = {0:.3f}'.format(v_med))

print(' Q1 = {0:.3f} Q3 = {1:.3f}'.format(v_25p, v_75p))

print(' 95% CI = [{0:.3f}, {1:.3f}]'.format(v_025p, v_975p))

txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))

txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))

txtout.write(' Q1 = {0:.3f} Q3 = {1:.3f}\n'.format(v_25p, v_75p))

txtout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(v_025p, v_975p))

plt.title('効果量', fontsize = 20)

sb.kdeplot(es_a, label = '主効果：{}'.format(NameA))

sb.kdeplot(es_b, label = '主効果：{}'.format(NameB))

sb.kdeplot(es_ab, label = '交互作用')

plt.yticks([])

plt.savefig(foutnm[:-4] + '_es.png')

SavedFiles.append(foutnm[:-4] + '_es.png')

plt.show()

print("\n\nGjk's =")

txtout.write("\n\nGjk's =\n")

for j in range(a):

for k in range(b):

print('\nGjk({0},{1}):'.format(LevelsA[j], LevelsB[k]))

txtout.write('\nGjk({0},{1}):\n'.format(LevelsA[j], LevelsB[k]))

v_mean = fit['G'][:, j, k].mean()

v_025p, v_25p, v_med, v_75p, v_975p = \

np.percentile(fit['G'][:, j, k], [2.5, 25, 50, 75, 97.5])

print(' 平均値 = {0:.3f}'.format(v_mean))

print(' 中央値 = {0:.3f}'.format(v_med))

print(' Q1 = {0:.3f} Q3 = {1:.3f}'.format(v_25p, v_75p))

print(' 95% CI = [{0:.3f}, {1:.3f}]'.format(v_025p, v_975p))

txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))

txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))

txtout.write(' Q1 = {0:.3f} Q3 = {1:.3f}\n'.format(v_25p, v_75p))

txtout.write(' 95% CI = [{0:.3f}, {1:.3f}]\n'.format(v_025p, v_975p))

txtout.close()

SavedFiles.append(txtoutnm)

print('\n\n保存ファイル:')

for nm in SavedFiles:

print(' ', nm)

print()

Up