Up

繰り返し１要因分散分析（変量効果）

 

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

ここで、変量効果であることは次式（２）で表わされる（Kirk, 1995, p. 201）。

上のモデルを、リスト１のStanスクリプトで表す。

 

リスト１　繰り返し１要因変量効果モデルのStanスクリプト（OneFctrRndmRep.stan）

 

data {

    int a;

    int N;

    real X[N, a];

}

parameters {

    real<lower = 0.0> sgm;

    real<lower = 0.0> sgm_a;

    real f_a[a];

    real mu[N];

}

model {

    sgm ~ exponential(0.001);

    sgm_a ~ exponential(0.001);

    for (i in 1:N) {

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

    }

    for (j in 1:a) {

        f_a[j] ~ normal(0.0, sgm_a);

    }

    for (i in 1:N)

        for (j in 1:a) {

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

    }

}

 

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

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

図１

 

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

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

 

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

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

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

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

Stanスクリプトのコンパイルとビルドの後、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_2e4aa54b61eb2ef3b27eeb6d40f98959.

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.49    0.03   0.96   4.94   5.84   6.37    7.0   8.68    828    1.0

sgm_a   21.13    2.78  26.28    5.6   9.73  14.09  21.85  85.03     89   1.04

f_a[1]  -2.97    2.34  17.94  -45.5  -6.25  -0.78   4.11  19.72     59   1.07

f_a[2]   9.31    2.34  17.88 -33.05   6.02  11.23  16.37  32.13     59   1.07

f_a[3] -12.89    2.36  18.05 -56.25 -16.12 -10.52  -5.59   9.44     59   1.07

f_a[4]  -4.06    2.34  17.97 -47.72  -7.37  -1.78   3.09  18.32     59   1.07

mu[1]  189.99    2.35  18.17 167.14  182.5 187.81 193.63 231.79     60   1.07

mu[2]  193.16    2.34  18.15 169.88 185.37 191.02  197.1 234.84     60   1.07

mu[3]  222.45    2.35  18.18 199.35 214.97 220.42 226.35 264.82     60   1.07

mu[4]  204.89    2.36  18.18 181.57 197.38 202.61 208.54 248.46     59   1.07

mu[5]   193.4    2.35  18.16 170.63 185.86 191.28 197.21 236.56     60   1.06

mu[6]  236.18    2.35  18.18 212.75 228.69 233.99  239.9 280.32     60   1.07

mu[7]  198.22    2.36  18.14 175.38  190.5 195.99 202.17 241.52     59   1.07

mu[8]  210.69    2.35  18.13 187.86  203.1 208.52 214.57 252.82     60   1.07

mu[9]  195.13    2.36   18.2 171.64 187.43 193.02 199.03 238.08     60   1.07

mu[10] 205.38    2.34  18.16 182.33 197.82 203.22  209.1  246.9     60   1.07

lp__   -102.0    0.24   3.71 -110.4 -104.1 -101.5 -99.36 -96.12    250   1.01

 

Samples were drawn using NUTS at Fri Aug  7 10:15:59 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).

 

 

水準の標準偏差:

    平均値 = 21.126

    中央値 = 14.086

    Q1 = 9.736    Q3 = 21.846

    95% CI = [5.603, 84.854]

 

 

f_a水準A1:

    平均値 = -2.971

    中央値 = -0.778

    Q1 = -6.248    Q3 = 4.108

    95% CI = [-44.988, 19.703]

 

 

f_a水準A2:

    平均値 = 9.314

    中央値 = 11.229

    Q1 = 6.026    Q3 = 16.367

    95% CI = [-33.030, 32.082]

 

 

f_a水準A3:

    平均値 = -12.894

    中央値 = -10.523

    Q1 = -16.120    Q3 = -5.589

    95% CI = [-56.176, 9.434]

 

 

f_a水準A4:

    平均値 = -4.061

    中央値 = -1.779

    Q1 = -7.371    Q3 = 3.092

    95% CI = [-47.280, 18.297]

 

 

効果量:

    平均値 = 3.339

    中央値 = 2.214

    Q1 = 1.485    Q3 = 3.448

    95% CI = [0.835, 14.050]

 

 

引用文献

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 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, 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 = 'OneFctrRndmRep.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']

 

print('\n水準の標準偏差:')

txtout.write('\n\n水準の標準偏差:\n')

v_mean = fit['sgm_a'].mean()

v_025p, v_Q1, v_med, v_Q3, v_975p = np.percentile(fit['sgm_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_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('水準の標準偏差 $\sigma_a$', fontsize = 20)

for j in range(1, a):

    sb.kdeplot(fit['sgm_a'], clip = (0.0, np.percentile(fit['sgm_a'], 99)))

plt.yticks([])

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

savedFiles = [foutnm[:-4] + '_sgm_a.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()

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

    L_c, U_c = np.percentile(fit['f_a'][:, j], [1, 99])

    sb.kdeplot(fit['f_a'][:,j], clip = (L_c, U_c), 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_sgm, sd in zip(fit['sgm_a'], fit['sgm']):

    #v = np.std(a_s)

    effcts.append(a_sgm / 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)

sb.kdeplot(effcts, clip = (0.0, np.percentile(effcts, 99)))

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