Up

繰り返し１要因分散分析（固定効果）

2026.07改訂

 

分散分析は、分散（平方和）の分割に基づく分析法であるが、目的は分散の分割に基づいて要因の効果の有無を調べることである。要因の効果は、要因の各水準の値が示されると分かり易い。各水準と従属変数の関係を確率モデルで表して、水準の値の事後分布を求めてみた。

 

ケースの水準における値が次の確率分布に従うとする。

 

 

ケースパラメータがであり、水準パラメータがである。の事前分布を平均値０の正規分布として原点を設定する。

このモデルに対するStanスクリプトをリスト１のように用意した。ファイルは、filesOneFctrFixedRep.zipにまとめた。

 

 

リスト１　１要因分散分析のStanスクリプト（OneFctrFixedRep.stan）

 

data {

    int J;

    int N;

    array[N, J] real X;

    real sd_subj;

    real sd_f;

    real mean_f;

}

parameters {

    array[N] real mu;

    array[J] real f;

    real<lower=0.0001> sgm;

}

model {

    sgm - 0.0001 ~ uniform(0, 1000);

    for (j in 1:J) {

        f[j] ~ normal(mean_f, sd_f); 

    }

    for (i in 1:N) {

        mu[i] ~  normal(0, sd_subj); 

        for (j in 1:J) {

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

        }

    }

}

 

 

リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプトをリスト２のように作成した。

 

リスト２用の入力データのサンプルを図１のように用意した。

 

図１

 

図１の形式に従うものであれば他のものでもよい。１行目に水準名、１列目にケース名を置き、水準名とケース名の交点のセルに該当するデータを書く。

 

リスト１とリスト２のスクリプトファイルと図１のデータファイル（他のExcelファイルでもよい）を同じフォルダに置く。そのフォルダにカレントディレクトリを移して、次のpythonコマンドを実行するとファイル名の設定が求められる。

 

(stan) ****/filesOneFctrFixedRep$ python OneFctrFixedRep.py

入力データファイル（*.xlsx）= datarep.xlsx

出力ファイル（*.txt）= results.txt

画像保存ファイル（*.png）= Fig.png

 

 

上の例では、入力データファイル名は図１のファイル名が設定されている。出力ファイル名と画像保存ファイル名は、ファイル拡張子がそれぞれ「*.txt」、「*.png」の任意のファイル名である。

ファイル名の設定が終わると、ファイルが読み込まれ、Stanスクリプトがコンパイルされる。コンパイル後、MCMCサンプリングが始まり、終了後、トレース図が表示される（図２）。

 

図２

 

図２のWindowを閉じると、図３のグラフが表示される。

 

図３

 

水準の効果（グラフではとなっている）の事後分布のグラフである。水準１（a1）と水準４（a4）が事後分布が重なっており、水準３（a3）の事後分布はそれらより左側、水準２（a2）のグラフはそれらより右側に位置している。すなわち、水準のパラメータ値は、水準３＜（水準１、水準４）＜水準２の関係にあることが分かる。

 

図３のWindowを閉じると、スクリプトの実行終了である。

 

 

なお、現在（2026.07.09）matplotlib 3.11で実行するとエラーが出る。matplotlib 3.10に換えると大丈夫である。matplotlib 3.10への交換は、次のコマンド

 

conda install matplotlib=3.10

 

でできた。

 

 

 

文献

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スクリプト例（OneFctrFixedRep.py）

 

import numpy as np

import scipy.stats as ss

import scipy.optimize as so

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, 1:]

a = len(A_name)

N = len(data) - 1

 

X = np.array(data[1:,1:], 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))

 

print(X)

txtout.write('\n{}\n'.format(X))

 

min_X = np.min(X)

max_X = np.max(X)

 

sd_subj = np.mean(np.std(X, axis=0))

print('sd_subj =', sd_subj)

 

mean_f = np.mean(X)

print('mean_f =', mean_f)

sd_f = np.mean(np.std(X, axis=1))

print('np.sd =', np.std(X, axis=1))

print('sd_f =', sd_f)

 

StanData = {'J': a, 'N': N, 'X': X,

            'sd_subj':sd_subj, 'sd_f':sd_f, 'mean_f':mean_f}

   

print('a,N', a, N,'\nX =\n', X)

 

model = CmdStanModel(stan_file='OneFctrFixedRep.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_dist(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())

 

for j in range(a):

    sb.kdeplot(fit[f"f[{j+1}]"], label=f"a{j+1}")

plt.xlabel('a', fontsize=16)

plt.legend(fontsize=14)

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

savedFiles += [foutnm[:-4] + '_a.png']

plt.show()

 

txtout.close()

savedFiles.append(txtoutnm)

 

print('\n保存ファイル:')

for nm in savedFiles:

    print('    ', nm)

print()

 

 

 

Up