２要因分散分析
非繰り返し、固定効果
2025.12：CmdStanPy対応に改訂
ベイズ分析ソフトStanのPythonインターフェースCmdStanPy対応に改訂した。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
非繰り返し２要因固定効果の構造方程式（structural model）は次式で与えられる（Winer et al., 1991, p. 292）。
ここで、
である（Kirk, 1995, p. 372）。
式（１）におけるパラメータに対して次式のようにおく。
このとき、モデル（１）は次式のように書ける。
制約条件（２）は、モデル（４）における
を式（３）の要因、
、
、
、に分割するために必要である。
式（４）に対してStanスクリプトをリスト１のように作成した。
リスト１では、パラメータは、データが共通の枠組みにおいて収集されたと考え、次のハイパー事前分布を設定している。
h_mu
~ normal(0.0, 1000.0);
h_sgm
~ exponential(0.001);
for (j in 1:a)
for
(k in 1:b) {
G[j][k] ~ normal(h_mu, h_sgm);
}
リスト１ 繰り返し１要因固定効果モデルのStanスクリプト（TwoFctrFixed.stan）
data {
int a;
int b;
int N;
array[N, a, b] real X;
}
parameters {
real h_mu;
// Hyper parameter
real<lower=0.0> h_sgm; // Hyper parameter
real<lower = 0.0> sgm;
array[a, b] real G;
}
model {
h_sgm
~ exponential(0.001);
sgm
~ exponential(0.001);
for (j in 1:a)
for
(k in 1:b) {
G[j][k] ~ normal(h_mu, h_sgm);
}
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スクリプトを本ウェブサイトの後半に示すリスト２のように用意した。
該当スクリプトは、以下のようである。
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 ismpls
in range(N_Samples):
vG
= g_smpls[:,:,ismpls]
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
スクリプトファイルなどは、圧縮ファイルfilesTwoFctrsFixed.zipとしてまとめた。ダウンロード展開すれば、自由に利用できる。PythonにおけるStan（PyStan2）の使い方については＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞第10章～第12章で説明している。
リスト１のStanスクリプトファイル（TwoFctrFixed.stan）、リスト２のPythonスクリプトファイル（TwoFctrFixed.py）、および入力データファイル（例えば、図１のファイルであれば、data2fctrs.xlsx）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされた環境で、次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan) *****/filesTwoFctrsFixed$ python TwoFctrFixed.py
上のコマンドを実行すると、入力データファイル名、出力ファイル名、画像保存ファイル名の設定が求められる。
(stan) *****/filesTwoFctrsFixed$
python TwoFctrFixed.py
入力データファイル（*.xlsx）= data2fctrs.xlsx
出力ファイル（*.txt）= Results.txt
画像保存ファイル（*.png）= Figure.png
入力データファイルは、図１に示す形式で用意する。
図１ データ例（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」を画像に合わせて修正した名前でファイルが保存される。
保存画像ファイル名を設定して「Enter」キーを押すと、データが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルには、多少時間が掛かる。っコンパイルが終了すると、MCMCサンプリングが始まり、終了後、トレース図が表示される（図２）。
図２
左側に事後分布のグラフ、右側にMCMCサンプリングのトレース図が表示されている。
図２のWindowを閉じると、図３のグラフが表示される。
図３
ハイパーパラメータh_muとh_sgmに対応する正規分布normal(h_mu, h_sgm)のグラフを、MCMCサンプルからランダムにh_muとh_sgmの組を100組選んで描画したものである。赤の太線は、それらの平均である。スクリプトは、以下の通りである。
rng = np.random.default_rng()
idx_slct = rng.choice(len(fit['h_mu']), size=100,
replace=False)
xcoord = np.linspace(np.min(X), np.max(X), 1000)
curve_list = []
for idx in idx_slct:
ycoord
= ss.norm.pdf(xcoord,
loc=fit['h_mu'][idx],
scale=fit['h_sgm'][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')
図３のWindowを閉じると、図４のグラフが表示される。
図４
要因パラメータの事後分布のグラフである。ラベルに表示されている数値は、事後分布の中央値である。事後分布の代表値として、Gelman et al. (2021)は、平均値より安定している中央値を勧めている。
左上のグラフは主効果Facto-Aのグラフである。３つの水準の主効果の事後分布が明瞭に分離されている。
右上のグラフは、主効果Factor-Bのグラフである。Factor-Bの主効果も明瞭に分離されている。
下段のグラフは、交互作用のグラフである。交互作用A1*B2とA3*B1の事後分布は０より右側にあり、交互作用A1*B1とA3*B2の事後分布は０より左側にある。交互作用A2*B1とA2*B2の事後分布は０を中ほどに含む形で分布している。
図４のWindowを閉じると、図５のグラフが表示される。
図５
パラメータの事後分布の中央値のグラフである。
図５のWindowを閉じると、図６のグラフが表示される。
図６
効果量は、主効果および交互作用に対して次式で与えられている（cf. Kirk, 1995, p. 181）。
図６のWindowを閉じると、図７のグラフが表示される。
図７
事後予測値のグラフである。パラメータの事後分布のMCMCサンプルのパラメータ値に対してデータの予測値をシミュレートしている。スクリプトは、以下の通りである。
for ir in
range(a):
for ic
in range(b):
v_gs = fit[f'G[{ir+1},{ic+1}]']
smpls = []
for v_g, v_sgm in zip(v_gs, fit['sgm']): # fit[f'lmbd[{ir+1},{ic+1}]']:
v = ss.norm.rvs(v_g,
v_sgm) # poisson.rvs(vlmbd)
smpls.append(v)
ax[ir][ic].hist(smpls)
ax[ir][ic].set_title(LevelsA[ir] + ', ' + LevelsB[ic], fontsize=16)
図７のWindowを閉じると、スクリプトの実行終了である。
出力ファイルの内容は以下のようになっている。
入力データファイル: 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
Mean
MCSE StdDev ... ESS_tail ESS_bulk/s R_hat
lp__ -159.80300 0.058920 2.341240 ... 2376.83 12700.7 1.000150
h_mu 59.89800 0.160584 9.320680 ... 2245.23 31639.6 1.002230
h_sgm 21.64470 0.208937 10.679400 ... 2646.63 27848.5 0.999666
sgm
6.79533 0.011108 0.681628 ... 2623.11 30608.9 1.000230
G[1,1] 65.96140 0.028598 2.088480 ... 3162.09 42975.2 1.002180
G[1,2] 53.77100 0.029326 2.179810 ... 2943.03 44777.3 1.000290
G[2,1] 80.37860 0.029464 2.178820 ... 2958.63 44342.2 1.001620
G[2,2] 61.77910 0.027231 2.109320 ... 2645.45 48500.1 1.000530
G[3,1] 64.22800 0.032378 2.169840 ... 2985.14 36323.0 1.002520
G[3,2] 33.44920 0.029522 2.091570 ... 3049.20 40487.5 1.000060
[10 rows x 11 columns]
chain__ iter__ draw__ lp__ ... G[3,1] G[1,2] G[2,2] G[3,2]
0
1.0
1.0 1.0
-158.34811 ... 65.357755 56.213242 60.645818 34.160891
1
1.0
2.0 2.0
-159.94804 ... 64.345770 53.987089 57.322241 31.154685
2
1.0
3.0 3.0
-157.19980 ... 62.971170 53.855988 63.414587 32.239574
3
1.0
4.0 4.0
-157.56732 ... 64.642397 53.539936 59.611714 35.837484
4
1.0
5.0 5.0
-159.19429 ... 61.224465 56.562786 61.069908 32.870657
... ... ... ...
... ...
...
...
... ...
3995 4.0 996.0 3996.0 -159.31018 ...
66.642974 56.162785 60.795602 32.156635
3996 4.0 997.0 3997.0 -158.18287 ...
61.980640 52.388816 62.261442 35.550628
3997 4.0 998.0 3998.0 -158.48260 ...
66.259092 51.844354 59.427080 33.756741
3998 4.0 999.0 3999.0 -159.04392 ...
65.997956 56.272166 64.131947 32.544665
3999 4.0 1000.0 4000.0 -159.11279 ...
63.323875 51.491366 59.778435 34.519537
[4000 rows x 19 columns]
主効果：Factor-A
A1：
平均値 = -0.062
中央値 = -0.045
Q1 = -0.871 Q3 = 0.792
95% CI = [-2.589, 2.339]
A2：
平均値 = 11.151
中央値 = 11.131
Q1 = 10.345 Q3 = 11.932
95% CI = [8.704, 13.600]
A3：
平均値 = -11.089
中央値 = -11.069
Q1 = -11.867 Q3 = -10.312
95% CI = [-13.499, -8.678]
主効果：Factor-B
B1：
平均値 = 10.261
中央値 = 10.245
Q1 = 9.659 Q3 = 10.845
95% CI = [8.556, 11.973]
B2：
平均値 = -10.261
中央値 = -10.245
Q1 = -10.845 Q3 = -9.659
95% CI = [-11.973, -8.556]
交互作用：
A1・B1：
平均値 = -4.166
中央値 = -4.178
Q1 = -4.941 Q3 = -3.372
95% CI = [-6.569, -1.816]
A1・B2：
平均値 = 4.166
中央値 = 4.178
Q1 = 3.372 Q3 = 4.941
95% CI = [1.816, 6.569]
A2・B1：
平均値 = -0.962
中央値 = -0.948
Q1 = -1.791 Q3 = -0.139
95% CI = [-3.420, 1.470]
A2・B2：
平均値 = 0.962
中央値 = 0.948
Q1 = 0.139 Q3 = 1.791
95% CI = [-1.470, 3.420]
A3・B1：
平均値 = 5.128
中央値 = 5.144
Q1 = 4.319 Q3 = 5.930
95% CI = [2.696, 7.537]
A3・B2：
平均値 = -5.128
中央値 = -5.144
Q1 = -5.930 Q3 = -4.319
95% CI = [-7.537, -2.696]
効果量（主効果）：Factor-A
平均値 = 1.356
中央値 = 1.357
Q1 = 1.229 Q3 = 1.474
95% CI = [1.003, 1.716]
効果量（主効果）：Factor-B
平均値 = 1.525
中央値 = 1.519
Q1 = 1.388 Q3 = 1.660
95% CI = [1.149, 1.907]
効果量（交互作用）：Factor-A・Factor-B
平均値 = 0.588
中央値 = 0.588
Q1 = 0.497 Q3 = 0.679
95% CI = [0.323, 0.855]
Gjk's =
Gjk(A1,B1):
平均値 = 65.961
中央値 = 65.951
Q1 = 64.522 Q3 = 67.348
95% CI = [61.991, 70.191]
Gjk(A1,B2):
平均値 = 53.771
中央値 = 53.795
Q1 = 52.293 Q3 = 55.282
95% CI = [49.540, 57.968]
Gjk(A2,B1):
平均値 = 80.379
中央値 = 80.383
Q1 = 78.936 Q3 = 81.864
95% CI = [76.135, 84.669]
Gjk(A2,B2):
平均値 = 61.779
中央値 = 61.796
Q1 = 60.296 Q3 = 63.229
95% CI = [57.663, 65.835]
Gjk(A3,B1):
平均値 = 64.228
中央値 = 64.219
Q1 = 62.816 Q3 = 65.666
95% CI = [60.013, 68.536]
Gjk(A3,B2):
平均値 = 33.449
中央値 = 33.462
Q1 = 32.038 Q3 = 34.881
95% CI = [29.387, 37.488]
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.
Winer, B. J., Brown, D.R., & Michels, K.M. (1991). Statistical principles in experimental design, third ed. McGraw-Hill, Inc.
リスト２ 非繰り返し２要因固定効果モデルのPythonスクリプト例（TwoFctrFixed.py）
#
# Yasuharu
Okamoto, 2020.08, 2025.12
#
import numpy
as np
import pandas as pd
import scipy.stats
as ss
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))
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}
model = CmdStanModel(stan_file='TwoFctrFixed.stan')
fit = model.sample(data=StanData)
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['h_mu']), size=100,
replace=False)
xcoord = np.linspace(np.min(X), np.max(X), 1000)
curve_list = []
for idx in idx_slct:
ycoord
= ss.norm.pdf(xcoord,
loc=fit['h_mu'][idx],
scale=fit['h_sgm'][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.title('正規分布（h_mu、h_sgm）', fontsize=16)
plt.legend()
plt.savefig(foutnm[:-4]
+ '_h.png')
SavedFiles.append(foutnm[:-4]
+ '_h.png')
plt.show()
txtout.write('\n{}\n'.format(fit))
N_Samples = len(fit['G[1,1]'])
g_smpls = []
for i in
range(a):
v_g
= []
for j in range(b):
v_g.append(fit[f'G[{i+1},{j+1}]'])
g_smpls.append(v_g)
g_smpls = np.array(g_smpls)
plt.figure(figsize=(14,
8))
#
# 主効果（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 ismpls
in range(N_Samples):
vG
= g_smpls[:,:,ismpls]
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.subplot(2,4,(1,2))
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))
sb.kdeplot(f_as[:, j], label = LevelsA[j]+f':
{v_med:.2f}')
plt.title('主効果：{}'.format(NameA), fontsize = 20)
plt.yticks([])
plt.legend()
plt.subplot(2,4,(3,4))
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))
sb.kdeplot(f_bs[:, k], label = LevelsB[k]+f':
{v_med:.2f}')
plt.title('主効果：{}'.format(NameB), fontsize
= 20)
plt.yticks([])
plt.legend()
plt.subplot(2,4,(6,7))
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))
sb.kdeplot(f_abs[:, j, k], label
=
LevelsA[j] + '*' + LevelsB[k]
+ f': {v_med:.2f}')
plt.title('交互作用：{0}・{1}'.format(NameA, NameB), fontsize = 20)
plt.legend()
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_ab.png')
SavedFiles.append(foutnm[:-4]
+ '_ab.png')
plt.show()
g_medians = np.median(g_smpls, axis=2)
print('g_mean
=\n', g_medians)
for ir in
range(a):
plt.plot(range(b),
g_medians[ir], label=LevelsA[ir])
plt.xticks(range(b), LevelsB,
fontsize=14)
plt.legend(fontsize=12)
plt.title('Gjk(median)',
fontsize=16)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_g.png')
SavedFiles.append(foutnm[:-4]
+ '_g.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_a,
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_a))
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_a))
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_b,
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_b))
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_b))
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_ab,
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_ab))
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_ab))
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
= f'主効果：{NameA} ({v_med_a:.2f})')
sb.kdeplot(es_b, label
= f'主効果：{NameB}
({v_med_b:.2f})')
sb.kdeplot(es_ab,
label = f'交互作用
({v_med_ab:.2f})')
plt.yticks([])
plt.legend()
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 = g_smpls[j][k].mean()
v_025p, v_25p, v_med, v_75p, v_975p = \
np.percentile(g_smpls[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))
fig, ax = plt.subplots(a,
b, sharex=True, sharey=True,figsize=(12,8))
for ir in
range(a):
for ic
in range(b):
v_gs = fit[f'G[{ir+1},{ic+1}]']
smpls = []
for v_g, v_sgm in zip(v_gs, fit['sgm']): # fit[f'lmbd[{ir+1},{ic+1}]']:
v = ss.norm.rvs(v_g,
v_sgm) # poisson.rvs(vlmbd)
smpls.append(v)
ax[ir][ic].hist(smpls)
ax[ir][ic].set_title(LevelsA[ir] + ', ' + LevelsB[ic], fontsize=16)
s_size = len(fit['sgm'])
fig.suptitle(f'事後予測度数、 サンプルサイズ={s_size}',
fontsize=18)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_pred.png')
SavedFiles.append(foutnm[:-4]
+ '_pred.png')
plt.show()
txtout.close()
SavedFiles.append(txtoutnm)
print('\n\n保存ファイル:')
for nm in SavedFiles:
print(' ', nm)
print()