Up

繰り返し１要因分散分析（固定効果）

繰り返し１要因固定効果のstrictly additive modelは次式で与えられる（Winer et al., 1991, p. 261）。

ここで、

である。

式（１）においては、パラメータ値は不定であるので、以下のようにパラメータを再設定する。

このとき、式（１）は以下のようになる。

式（２）、（４）より、

である。

式（５）、（６）のモデルをリスト１のStanスクリプトで表す。

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

// Yasuharu Okamoto, 2020.08

data {

int a;

int N;

real X[N, a];

}

parameters {

real<lower = 0.0> sgm;

real tau[a - 1];

real mu[N];

}

transformed parameters {

real f_a[a];

real v = 0.0;

for (j in 2:a) {

v += tau[j - 1];

}

f_a[1] = -v / a;

for (j in 2:a) {

f_a[j] = tau[j - 1] + f_a[1];

}

}

model {

sgm ~ exponential(0.001);

for (i in 1:N) {

mu[i] ~ normal(0.0, 1000.0);

}

for (j in 1:(a-1)){

tau[j] ~ normal(0.0, 1000.0);

}

for (i in 1:N)

for (j in 1:a) {

if (j == 1) {

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

}

else {

X[i][j] ~ normal(mu[i] + tau[j-1], sgm);

}

}

}

リスト１のStanスクリプトを用いて繰り返し１要因固定効果モデルのベイズ分析を行うスクリプトを、本ウェブサイト後半にリスト２として示した。これらのスクリプトファイルなどは、圧縮ファイルfilesOneFctrFixedRep.zipとしてまとめた。ダウンロード展開すれば自由に利用できる。PythonにおけるStanの使い方については＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞第10章〜第12章で説明している。

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

図１

入力データファイルは、図２のように用意する。

図２　入力データ例（datarep.xlsx）

１行目に、変数名と要因の水準名を書き、データ値は２行目から設定する。1列目は、実験参加者など各データ組の識別値であるが、同じ値を複数回用いてもよい。文字列でも通し番号でもよい。２列目からそれぞれの水準のデータ値を書く。

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

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

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

MCMCサンプリングが終了すると、各水準の値としての事後分布のグラフが表示される（図３）。

図３

図３のウィンドウを閉じると、次に、各水準の値としての事後分布が表示される（図４）。

図４

図４では、はf_aで表されている。

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

図５

効果量は次式で与えられている（cf. Kirk, 1995, p. 181）。

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

図６

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

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

図６

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

図７

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

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

A = ['A1' 'A2' 'A3' 'A4']

a = 4

N = 10

[[182. 199. 185. 183.]

[186. 212. 184. 180.]

[222. 233. 207. 217.]

[207. 209. 186. 207.]

[187. 203. 181. 192.]

[235. 245. 227. 227.]

[197. 206. 189. 190.]

[205. 228. 197. 202.]

[195. 198. 181. 196.]

[204. 213. 180. 214.]]

Inference for Stan model: anon_model_e59af728c281d219c4b020b880acf617.

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.45 0.02 0.93 4.94 5.79 6.36 6.97 8.52 2519 1.0

tau[1] 12.54 0.07 2.93 6.76 10.56 12.48 14.58 18.26 1774 1.0

tau[2] -10.36 0.07 2.88 -16.13 -12.21 -10.35 -8.46 -4.62 1610 1.0

tau[3] -1.24 0.07 2.87 -6.88 -3.12 -1.28 0.68 4.41 1866 1.0

mu[1] 186.95 0.07 3.77 179.59 184.48 186.9 189.39 194.45 2769 1.0

mu[2] 190.24 0.07 3.67 182.88 187.81 190.33 192.73 197.4 2819 1.0

mu[3] 219.56 0.07 3.61 212.47 217.19 219.49 221.94 226.75 2612 1.0

mu[4] 202.02 0.07 3.68 194.81 199.59 202.06 204.45 209.1 2517 1.0

mu[5] 190.47 0.07 3.72 183.2 187.95 190.38 192.99 197.89 2724 1.0

mu[6] 233.23 0.07 3.66 226.09 230.83 233.24 235.62 240.59 2722 1.0

mu[7] 195.27 0.07 3.63 188.29 192.77 195.26 197.69 202.35 2802 1.0

mu[8] 207.77 0.07 3.73 200.09 205.31 207.85 210.28 214.79 2598 1.0

mu[9] 192.23 0.07 3.71 185.01 189.81 192.27 194.73 199.54 2567 1.0

mu[10] 202.56 0.07 3.71 195.26 200.03 202.67 204.95 209.77 2815 1.0

f_a[1] -0.23 0.05 1.77 -3.69 -1.42 -0.21 0.95 3.23 1242 1.01

f_a[2] 12.3 0.02 1.79 8.79 11.14 12.28 13.47 15.79 5185 1.0

f_a[3] -10.6 0.03 1.77 -14.31 -11.7 -10.61 -9.43 -7.06 4594 1.0

f_a[4] -1.48 0.02 1.78 -5.03 -2.65 -1.5 -0.32 2.17 5569 1.0

v 0.93 0.2 7.07 -12.92 -3.81 0.84 5.69 14.74 1242 1.01

lp__ -92.08 0.09 3.3 -99.51 -93.95 -91.66 -89.66 -86.89 1273 1.0

Samples were drawn using NUTS at Thu Aug 6 14:13:28 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).

tau-水準A1:

0に固定

tau-水準A2:

平均値 = 12.538

中央値 = 12.483

Q1 = 10.558 Q3 = 14.579

95% CI = [6.762, 18.259]

tau-水準A3:

平均値 = -10.362

中央値 = -10.348

Q1 = -12.213 Q3 = -8.459

95% CI = [-16.116, -4.628]

tau-水準A4:

平均値 = -1.242

中央値 = -1.283

Q1 = -3.118 Q3 = 0.678

95% CI = [-6.881, 4.405]

f_a水準A1:

平均値 = -0.234

中央値 = -0.210

Q1 = -1.423 Q3 = 0.951

95% CI = [-3.683, 3.226]

f_a水準A2:

平均値 = 12.305

中央値 = 12.280

Q1 = 11.139 Q3 = 13.466

95% CI = [8.793, 15.783]

f_a水準A3:

平均値 = -10.595

中央値 = -10.612

Q1 = -11.702 Q3 = -9.430

95% CI = [-14.306, -7.058]

f_a水準A4:

平均値 = -1.476

中央値 = -1.504

Q1 = -2.645 Q3 = -0.322

95% CI = [-5.027, 2.164]

効果量:

平均値 = 1.308

中央値 = 1.303

Q1 = 1.148 Q3 = 1.462

95% CI = [0.866, 1.773]

引用文献

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スクリプト例

import numpy as np

import scipy.stats as ss

import scipy.optimize as so

import pandas as pd

import matplotlib.pyplot as plt

import pystan

import seaborn as ss

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

A_name = data[0, 1:]

a = len(A_name)

N = len(data) - 1

X = np.empty((N, a), dtype = float)

print('A = ', A_name)

print('a = ', a)

print('N = ', N)

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

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

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

for i in range(N):

for j in range(a):

X[i][j] = data[i + 1][j + 1]

print(X)

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

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

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

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

print(fit)

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

# 日本語表記の設定

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

for j in range(0, a):

if j == 0:

print('\ntau-{}:'.format(A_name[j]))

print(' 0に固定')

txtout.write('\n\ntau-水準{}:\n'.format(A_name[j]))

txtout.write(' 0に固定\n')

else:

print('\ntau-{}:'.format(A_name[j]))

txtout.write('\n\ntau-水準{}:\n'.format(A_name[j]))

v_mean = fit['tau'][:, j-1].mean()

#A_mean.append(v_mean)

v_025p, v_Q1, v_med, v_Q3, v_975p = np.percentile(fit['tau'][:, j-1],

[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_Q1, v_Q3))

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_Q1, v_Q3))

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

plt.figure(figsize = (8,5))

plt.title('水準の値 $\\tau$' + \

'\n$\\tau$({}) = 0, 固定'.format(A_name[0]), fontsize = 20)

for j in range(1, a):

ss.kdeplot(fit['tau'][:,j-1], label = A_name[j])

plt.yticks([])

plt.legend()

plt.tight_layout()

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

savedFiles = [foutnm[:-4] + '_tau.png']

plt.show()

for j in range(a):

print('\nf_a水準{}:'.format(A_name[j]))

txtout.write('\n\nf_a水準{}:\n'.format(A_name[j]))

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

#A_mean.append(v_mean)

v_025p, v_Q1, v_med, v_Q3, v_975p = np.percentile(fit['f_a'][:, 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_Q1, v_Q3))

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_Q1, v_Q3))

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

plt.figure(figsize = (8,5))

plt.title('水準の値 f_a', fontsize = 20)

for j in range(a):

ss.kdeplot(fit['f_a'][:,j], label = A_name[j])

plt.yticks([])

plt.legend()

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

savedFiles.append(foutnm[:-4] + '_f_a.png')

plt.show()

effcts = []

for a_s, sd in zip(fit['f_a'], fit['sgm']):

v = np.std(a_s)

effcts.append(v / sd)

effcts = np.array(effcts)

print('\n効果量:')

txtout.write('\n\n効果量:\n')

v_mean = effcts.mean()

v_025p, v_Q1, v_med, v_Q3, v_975p = np.percentile(effcts, [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_Q1, v_Q3))

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_Q1, v_Q3))

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

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

ss.kdeplot(effcts)

plt.ylabel([])

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

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

plt.show()

txtout.close()

savedFiles.append(txtoutnm)

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

for nm in savedFiles:

print(' ', nm)

print()

Up