Up

１要因分散分析

非繰り返し・固定効果

 

１要因非繰り返し固定効果の構造方程式は次式で与えられる（Winer et al., 1991, p. 84）。

式（１）においては、パラメータ値は不定であるので、をとおいてパラメータが同定可能であるようにすると、次式（２）となる。

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

 

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

 

data {

    int a;

    int N;

    real X[N, a];

}

parameters {

    real<lower = 0.0> sgm;

    real f_a[a];

}

model {

    sgm ~ exponential(0.001);

    for (j in 1:a){

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

    }

    for (i in 1:N)

        for (j in 1:a) {

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

    }

}

 

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

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

図１

 

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

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

 

１行目に、要因の水準名を書き、データ値は２行目から設定する。

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

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

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

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

図３

 

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

図４

 

効果量は次式で与えられている（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_80f4c87c00595f696cb78229e8f3641c.

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      8.81    0.02   1.09   6.94   8.07   8.73   9.44  11.24   3581    1.0

f_a[1]  70.01    0.04   2.85  64.16   68.2  69.99  71.77  75.88   4512    1.0

f_a[2]  81.46    0.04   2.82  75.99  79.56  81.47  83.37  87.02   5678    1.0

f_a[3]  58.18    0.04   2.74   52.7  56.39  58.21  60.04  63.41   5505    1.0

f_a[4]  68.72    0.04   2.84  63.04  66.87  68.74  70.55  74.21   5115    1.0

lp__   -104.2    0.05    1.7 -108.6 -105.0 -103.8 -102.9 -102.0   1371    1.0

 

Samples were drawn using NUTS at Tue Aug  4 15:48:07 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).

 

 

水準A1:

    平均値 = 70.013

    中央値 = 69.988

    Q1 = 68.200    Q3 = 71.766

    95% CI = [64.168, 75.878]

 

 

水準A2:

    平均値 = 81.459

    中央値 = 81.468

    Q1 = 79.563    Q3 = 83.370

    95% CI = [75.988, 87.021]

 

 

水準A3:

    平均値 = 58.182

    中央値 = 58.214

    Q1 = 56.390    Q3 = 60.036

    95% CI = [52.708, 63.415]

 

 

水準A4:

    平均値 = 68.717

    中央値 = 68.744

    Q1 = 66.866    Q3 = 70.549

    95% CI = [63.048, 74.200]

 

 

効果量:

    平均値 = 0.975

    中央値 = 0.971

    Q1 = 0.846    Q3 = 1.098

    95% CI = [0.625, 1.352]

 

 

 

引用文献

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

import pickle

 

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 = 'OneFctrFixed.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(a):

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

    txtout.write('\n\n水準{}:\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('水準の値', fontsize = 20)

for j in range(a):

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

plt.yticks([])

plt.legend()

plt.savefig(foutnm)

savedFiles = [foutnm]

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