１要因分散分析
非繰り返し・固定効果
2025.12： CmdStanPy対応に改訂
ベイズ分析のためのソフトStanのPythonインターフェースCmdStanPy用に改訂した。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
１要因非繰り返し固定効果の構造方程式は次式で与えられる（Winer et al., 1991, p. 84）。
式（１）においては、パラメータ値は不定であるので、
を
とおいてパラメータが同定可能であるようにすると、次式（２）となる。
式（２）のモデルをリスト１のStanスクリプトで表す。
リスト１ １要因非繰り返し固定効果モデルのStanスクリプト（OneFctrFixed.stan）
data {
int a;
int N;
array[N, a] real X;
}
parameters {
real<lower=1.0e-9>
sgm0;
// Hyper parameter
real mu0;
// Hyper parameter
real<lower = 1.0e-9>
sgm;
array[a] real f_a;
}
model {
mu0 ~ normal(0.0,
1000.0);
// Hyper prior distribution
sgm0 - 1.0e-9 ~
exponential(0.01); // Hyper prior distribution
for (j in 1:a){
f_a[j] ~ normal(mu0,
sgm0); // Prior distribution
}
sgm - 1.0e-9 ~
exponential(0.01);
// Prior distribution
for (i in 1:N)
for
(j in 1:a) {
X[i][j] ~ normal(f_a[j], sgm);
}
}
データが共通の枠組みで収集されていることを踏まえて、水準の効果f_a[j]を共通のハイパーパラメータmu0とsgm0を持つ事前分布からのサンプルであるとしている。
（３） f_a[j] ~ normal(mu0, sgm0);
上式のモデル（３）は、f_a[j]に対して事前分布を設定するものであるが、水準値f_a[j]が正規分布normal(mu0, sgm0)に従うrandom effectであると見ることもできる。すなわち、fixed effectの階層モデル（３）は、random effectモデルを表しているとも考えられ、ベイズモデルではfixed effectモデルとrandom effectモデルは事前分布という確率モデルにより統合されると考えられる。
リスト１のStanスクリプトを用いて１要因非繰り返し固定効果モデルのベイズ分析を行うスクリプトを本ウェブサイト後半にリスト２として示した。これらのスクリプトファイルなどは、圧縮ファイルfilesOneFctrFixed.zipとしてまとめた。ダウンロード展開すれば自由に利用できる。PythonにおけるStan（PyStan2）の使い方については＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞第10章〜第12章で説明している。
リスト１のStanスクリプトファイル（OneFctrFixed.stan）、リスト２のPythonスクリプトファイル（OneFctrFixed.py）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされている環境において次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan)
*****/filesOneFctrFixed$ python OneFctrFixed.py
上のコマンドを実行すると、入力データファイル名、出力ファイル名、画像保存ファイル名の設定が求められる。
(stan) *****/filesOneFctrFixed$ python
OneFctrFixed.py
入力データファイル（*.xlsx）= data.xlsx
出力ファイル（*.txt）= Results.txt
画像保存ファイル（*.png）= Figure.png
上の例では、入力データファイル名として図１のファイル名が設定されている。
図１ 入力データ例（data.xlsx）
入力データファイルは、Excelファイルとして、１行目に、要因の水準名を書き、データ値は２行目から設定する。
入力データファイル名を設定して「Enter」キーを押すと、次に計算結果の出力ファイル名の設定が求められる。ファイル拡張子が「.txt」である任意のテキストファイル名を設定する。
出力ファイル名の設定後、「Enter」キーを押すと、グラフの画像保存ファイル名の設定が求められる。ファイル拡張子が「.png」である任意のファイル名を設定する。設定したファイル名の「.png」を画像ファイルに合わせて変更したファイル名でそれぞれのファイルが保存される。
保存画像ファイル名を設定して「Enter」キーを押すと、データが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルに多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、終了後、トレース図が表示される（図２）。
図２
左側に事後分布のグラフ、右側にMCMCサンプリングのトレース図が表示されている。
図２のWindowを閉じると、図３のグラフが表示される。
図３
水準の値f_a[j]のハイパー分布（３）の正規分布のグラフを、MCMCサンプルからハイパーパラメータmu0とsgm0を100対ランダムに選び描画したものである。赤の太い曲線は、それらの平均値である。スクリプトを以下に示す。
idx_slct = rng.choice(len(fit['mu0']),
size=100, replace=False)
xcoord = np.linspace(0, 150, 1000)
curve_list = []
for idx in idx_slct:
ycoord = ss.norm.pdf(xcoord,
loc=fit['mu0'][idx], scale=fit['sgm0'][idx])
plt.plot(xcoord,ycoord,
alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(xcoord, mn_curve, c='r',
label='Mean values')
plt.legend()
plt.title('Normal Distributions with
Hyperparameters', fontsize=16)
plt.tight_layout()
plt.savefig(foutnm[:-4] +
'_h_normal.png')
savedFiles.append(foutnm[:-4] +
'_h_normal.png')
plt.show()
図３のWindowを閉じると、図４のグラフが表示される。
図４
水準の値f_a[j]の事後分布のグラフである。グラフのラベルに表示されている値は、事後分布の中央値である。事後分布の代表値として、Gelman et al.(2021)は、平均値より安定しているとして中央値を勧めている。
図4のWindowを閉じると、図５のグラフが表示される。
図５
効果量の事後分布のグラフである。効果量は次式で与えられている（Kirk, 1995； Myers & Well, 2003）。
図５のWindowを閉じると、図６のグラフが表示される。
図６
データの事後予測値のグラフである。スクリプトを以下に示す。
for j in range(a):
v_smpls = []
for va, vs in
zip(fit[f'f_a[{j+1}]'], fit ['sgm']):
v =
ss.norm.rvs(va, vs)
v_smpls.append(v)
sb.kdeplot(v_smpls,
label=A_name[j])
plt.title('Predictive Posterior
Values', fontsize=16)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig(foutnm[:-4] + '_pred.png')
savedFiles.append(foutnm[:-4] +
'_pred.png')
plt.show()
図６のWindowを閉じると、スクリプトの実行終了である。
出力ファイルの内容は以下のようになっている。
入力データファイル: 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.]]
Mean
MCSE StdDev MAD ...
ESS_bulk ESS_tail ESS_bulk/s R_hat
lp__ -113.84600 0.061333 2.14481 1.93476 ... 1353.62 1499.370 7647.60 1.001060
sgm0 16.04620 0.404799 12.87410 6.30120 ... 1463.15 1199.230 8266.38 1.002210
mu0 70.14560 0.298139 9.59263 6.42210 ... 1719.94 827.311 9717.18 1.001980
sgm
8.84134 0.020753 1.10384 1.07073 ... 2951.92 2442.930 16677.50 1.002180
f_a[1] 69.92040 0.045344 2.72009 2.65213 ... 3637.73 2737.370 20552.20 1.001610
f_a[2] 80.65930 0.049099 2.81496 2.82273 ... 3369.45 2666.280 19036.40 1.000450
f_a[3] 58.94630 0.054253 2.89181 2.87354 ... 2882.78 2145.210 16286.90 1.001610
f_a[4] 68.76780 0.045951 2.69093 2.67248 ... 3456.75 2561.910 19529.70 0.999549
[8 rows x 11 columns]
水準A1:
平均値 = 69.920
中央値 = 69.901
Q1 = 68.135 Q3 = 71.731
95% CI = [64.517, 75.161]
水準A2:
平均値 = 80.659
中央値 = 80.636
Q1 = 78.750 Q3 = 82.558
95% CI = [75.179, 85.987]
水準A3:
平均値 = 58.946
中央値 = 58.974
Q1 = 56.976 Q3 = 60.848
95% CI = [53.329, 64.793]
水準A4:
平均値 = 68.768
中央値 = 68.811
Q1 = 66.967 Q3 = 70.571
95% CI = [63.476, 74.017]
効果量:
平均値 = 0.911
中央値 = 0.911
Q1 = 0.778 Q3 = 1.043
95% CI = [0.510, 1.292]
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
Kirk, R.E. (1995). Experimental design: Procedures for the behavioral sciences, third ed. Brooks/Cole Publishing Company.
Myers, J. L. & Well, A. D. (2003). Research design and statistical analysis, second ed. Lawrence Erlbaum Associates, Publishers.
Winer, B. J., Brown, D.R., & Michels, K.M. (1991). Statistical principles in experimental design, third ed. McGraw-Hill, Inc.
リスト２ １要因非繰り返し固定効果モデルのPythonスクリプト例（OneFctrFixed.py）
import numpy as np
import scipy.stats as ss
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
from cmdstanpy import CmdStanModel
import seaborn as sb
import arviz as az
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}
model =
CmdStanModel(stan_file='OneFctrFixed.stan')
fit = model.sample(data=StanData,
adapt_delta=0.99)
print(fit.diagnose())
print(fit.summary())
txtout.write(f'{fit.summary()}')
infdata = az.from_cmdstanpy(fit) # Arviz InfereceData型へ変換
az.plot_trace(infdata)
plt.tight_layout()
plt.savefig(foutnm[:-4] +
'_trace.png')
savedFiles = [foutnm[:-4] +
'_trace.png']
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型への変換
print(fit.keys())
rng = np.random.default_rng()
idx_slct = rng.choice(len(fit['mu0']),
size=100, replace=False)
xcoord = np.linspace(0, 150, 1000)
curve_list = []
for idx in idx_slct:
ycoord = ss.norm.pdf(xcoord,
loc=fit['mu0'][idx], scale=fit['sgm0'][idx])
plt.plot(xcoord,ycoord,
alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(xcoord, mn_curve, c='r',
label='Mean values')
plt.legend()
plt.title('Normal Distributions with
Hyperparameters', fontsize=16)
plt.tight_layout()
plt.savefig(foutnm[:-4] +
'_h_normal.png')
savedFiles.append(foutnm[:-4] +
'_h_normal.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'f_a[{j+1}]'].mean()
v_025p, v_Q1, v_med, v_Q3,
v_975p = \
np.percentile(fit[f'f_a[{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('水準の値', fontsize = 20)
for j in range(a):
v_med =
np.median(fit[f'f_a[{j+1}]'])
sb.kdeplot(fit[f'f_a[{j+1}]'],
label = A_name[j] + f': {v_med:.2f}')
plt.yticks([])
plt.legend()
plt.savefig(foutnm)
savedFiles.append(foutnm)
plt.show()
a_smpls = []
for i in range(a):
a_smpls.append(fit[f'f_a[{i+1}]'])
a_smpls = np.array(a_smpls)
print('a_smpls shape =',
np.shape(a_smpls))
print(a_smpls)
effcts = []
for i in range(len(a_smpls[0])):
v =
np.std(a_smpls[:,i])
sd = fit['sgm'][i]
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('効果量' + '\n' +
f'Q1={v_Q1:.2f},
Med={v_med:.2f},
Q3={v_Q3:.2f}', fontsize = 20)
sb.kdeplot(effcts)
plt.ylabel([])
plt.tight_layout()
plt.savefig(foutnm[:-4] + '_es.png')
savedFiles.append(foutnm[:-4] +
'_es.png')
plt.show()
for j in range(a):
v_smpls = []
for va, vs in
zip(fit[f'f_a[{j+1}]'], fit ['sgm']):
v =
ss.norm.rvs(va, vs)
v_smpls.append(v)
sb.kdeplot(v_smpls,
label=A_name[j])
plt.title('Predictive Posterior
Values', fontsize=16)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig(foutnm[:-4] + '_pred.png')
savedFiles.append(foutnm[:-4] +
'_pred.png')
plt.show()
txtout.close()
savedFiles.append(txtoutnm)
print('\n保存ファイル:')
for nm in savedFiles:
print(' ', nm)
print()