１要因分散分析
非繰り返し・変量効果
１要因非繰り返し変量効果の構造方程式は次式で与えられる（Winer et al., 1991, p. 92）。
変量効果モデルであることは、次式
で表される。
式（１）においては、パラメータ値は不定であるので、
を
とおいてパラメータが同定可能であるようにすると、次式（３）となる。
このとき、式（２）は式（４）で表される。
上のモデルをリスト１のStanスクリプトで表す。
リスト１ １要因非繰り返し変量効果モデルのStanスクリプト（OneFctrRndm.stan）
data {
int
a;
int
N;
real X[N, a];
}
parameters {
real<lower = 0.0> sgm;
real<lower = 0.0> sgm_tau;
real mu;
real f_a[a];
}
model {
sgm
~ exponential(0.001);
sgm_tau
~ exponential(0.001);
mu ~ normal(0.0, 1000.0);
for (j in 1:a){
f_a[j] ~ normal(mu, sgm_tau);
}
for (i
in 1:N)
for
(j in 1:a) {
X[i][j] ~ normal(f_a[j], sgm);
}
}
リスト１のStanスクリプトを用いて１要因非繰り返し変量効果モデルのベイズ分析を行うスクリプトを本ウェブサイト後半にリスト２として示した。これらのスクリプトファイルなどは、圧縮ファイルfilesOneFctrRndm.zipとしてまとめた。ダウンロード展開すれば自由に利用できる。PythonにおけるStanの使い方については＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞第10章〜第12章で説明している。
リスト２のスクリプトを実行すると、まず、入力データファイル名の設定が求められる（図１）。
図１
入力データファイルは、図２のように用意する。
図２ 入力データ例（data.xlsx）
１行目に、要因の水準名を書き、データ値は２行目から設定する。
入力データファイル名を設定して「Enter」キーを押すと、次に計算結果の出力ファイル名の設定が求められる（図１）。ファイル拡張子が「.txt」である任意のテキストファイル名を設定する。
出力ファイル名の設定後、「Enter」キーを押すと、グラフの画像保存ファイル名の設定が求められる（図１）。ファイル拡張子が「.png」である任意のファイル名を設定する。設定したファイル名と設定ファイル名の「.png」を「_es.png」の２つのファイルが保存される。
保存画像ファイル名を設定して「Enter」キーを押すと、データが読み込まれ、計算が始まる。
MCMCサンプリングが終了すると、変量効果の標準偏差
の事後分布が表示される（図３）。
図３
図３のウィンドウを閉じると、次に各水準の値の事後分布のグラフが表示される（図４）。
図４
図４のウィンドウを閉じると、効果量の事後分布のグラフが表示される（図５）。
図５
効果量は次式で与えられている（cf. Kirk,
1995, p. 181）。
図５のウィンドウを閉じると、スクリプトの実行終了となる（図６）。
図６
保存されたファイル名が表示されている。これらのファイルは、スクリプトの実行終了後、開くことができる。
出力テキストファイルをPython IDLEで開こうとすると、図６のメッセージが表示されることがある。
図６
「s-jis」と設定して、「OK」ボタンをクリックすると開くことができる（図７）。
図７
出力ファイルの内容は以下のようになっている。
入力データファイル: data.xlsx
A = ['A1' 'A2' 'A3' 'A4']
a = 4
N = 10
[[66. 77. 60. 80.]
[53. 87. 59. 76.]
[63. 95. 66. 72.]
[64. 81. 57. 56.]
[76. 75. 43. 58.]
[69. 85. 55. 70.]
[77. 84. 61. 64.]
[86. 73. 61. 76.]
[75. 94. 66. 74.]
[71. 64. 54. 61.]]
Inference for Stan model:
anon_model_e631552d3878280c231e5ff128519cc8.
4 chains, each with iter=5000; warmup=3000; thin=1;
post-warmup draws per chain=2000,
total post-warmup draws=8000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sgm 8.87 0.01 1.09 7.07 8.1 8.77 9.53 11.29 5606 1.0
sgm_tau 19.7 0.97 26.8 5.31 9.46 13.61 21.07 72.38 765 1.01
mu 69.02 0.72 17.11 44.24 65.02 69.71 74.48 92.08 564 1.01
f_a[1] 69.95 0.03 2.75 64.61 68.07 69.96 71.81 75.32 6724 1.0
f_a[2] 80.8 0.03 2.78 75.28 78.98 80.82 82.67 86.2 6473 1.0
f_a[3] 58.93 0.04 2.94 53.25 56.94 58.92 60.92 64.69 5680 1.0
f_a[4] 68.75 0.03 2.79 63.11 66.94 68.77 70.63 74.24 7030 1.0
lp__ -113.9 0.06 2.24 -119.3 -115.1 -113.5 -112.2
-110.7 1506 1.0
Samples were drawn using NUTS at Wed
Aug 5 12:33: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).
変量効果の標準偏差
平均値 = 19.700
中央値 = 13.609
Q1 = 9.464 Q3 = 21.060
95% CI = [5.314, 72.159]
水準A1:
平均値 = 69.955
中央値 = 69.964
Q1 = 68.069 Q3 = 71.810
95% CI = [64.611, 75.321]
水準A2:
平均値 = 80.799
中央値 = 80.815
Q1 = 78.978 Q3 = 82.668
95% CI = [75.280, 86.196]
水準A3:
平均値 = 58.933
中央値 = 58.921
Q1 = 56.935 Q3 = 60.915
95% CI = [53.255, 64.693]
水準A4:
平均値 = 68.747
中央値 = 68.770
Q1 = 66.936 Q3 = 70.626
95% CI = [63.120, 74.236]
効果量:
平均値 = 2.260
中央値 = 1.552
Q1 = 1.072 Q3 = 2.442
95% CI = [0.576, 8.257]
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 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))
A_name = data[0]
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]
print(X)
txtout.write('\n{}\n'.format(X))
StanData = {'a': a, 'N': N, 'X': X}
sm = pystan.StanModel(file
= 'OneFctrRndm.stan')
fit = sm.sampling(data
= StanData, n_jobs = 1, iter = 5000, warmup = 3000)
print(fit)
txtout.write('\n{}\n'.format(fit))
# 日本語表記の設定
plt.rcParams['font.family']
= ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
v_mean = fit['sgm_tau'].mean()
v_025p, v_Q1, v_med,
v_Q3, v_975p = np.percentile(fit['sgm_tau'],
[2.5, 25, 50, 75, 97.5])
print('\n変量効果の標準偏差：')
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('\n\n変量効果の標準偏差\n')
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('変量効果の標準偏差の事後分布' + \
'\n平均値 =
{0:.3f} 中央値 = {1:.3f}'.format(v_mean, v_med), fontsize = 20)
sb.kdeplot(fit['sgm_tau'],
clip = (0.0, np.percentile(fit['sgm_tau'],
99)))
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_sgm_tau.png')
savedFiles = [foutnm[:-4]
+ '_sgm_tau.png']
plt.show()
for j in range(a):
print('\n水準{}:'.format(A_name[j]))
txtout.write('\n\n水準{}:\n'.format(A_name[j]))
v_mean
= fit['f_a'][:, j].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('水準の値', fontsize = 20)
for j in range(a):
sb.kdeplot(fit['f_a'][:,j], label = A_name[j])
plt.yticks([])
plt.legend()
plt.savefig(foutnm[:-4]
+ '_fa.png')
savedFiles.append(foutnm[:-4]
+ '_fa.png')
plt.show()
effcts = []
for st, sd in zip(fit['sgm_tau'], fit['sgm']):
effcts.append(st / 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('効果量' + \
'\n平均値 =
{0:.2f} 中央値 = {1:.2f}'.format(v_mean, v_med), fontsize = 20)
sb.kdeplot(effcts,
clip = (0.0, np.percentile(effcts,
99)))
plt.yticks([])
plt.tight_layout()
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