Up

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

2026.07改訂

 

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

 

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

 

 

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

水準の効果に変量モデルの考え方を適用して

 

 

とおいてみた。

 

上のモデルに基づいてStanスクリプトをリスト１のように作成した。

 

 

リスト１　Stanスクリプト（OneFacRepRndm.stan）

 

data {

    int J;

    int N;

    array[N, J] real X;

    real sd_subj;

    real sd_f;

    real mean_f;

}

parameters {

    real mu0;                   //  Hyper parameter

    real<lower=0.0001> sd0;     //  Hyper parameter

    array[N] real mu;

    array[J] real f;

    real<lower=0.0001> sgm;

}

model {

    mu0 ~ normal(mean_f, sd_f);

    sd0 - 0.0001 ~ exponential(1/sd_f);

    sgm - 0.0001 ~ uniform(0, 1000);

    for (j in 1:J) {

        f[j] ~ normal(mu0, sd0); 

    }

    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スクリプトをリスト２のように作成した。ファイルは、filesOneFacRepRndm.zipにまとめた。

 

サンプルデータとして図１に示すものを作成した。

 

図１

 

入力データファイルは、図１の形式のものであればよい。1行目に変数名を書く。１列目は、ケースIDである。データは、２行目以降に書く。

 

リスト１のStanスクリプトファイル（OneFacRepRndm.stan）、リスト２のPythonスクリプトファイル（OneFacRepRndm.py）、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリを移して、以下のようにpythonコマンドを実行する。

 

(stan) ****/filesOneFacRepRndm$ python OneFacRepRndm.py

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

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

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

 

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

ファイル名を設定すると、Stanスクリプトがコンパイルされて、MCMCサンプリングが始まる。MCMCサンプリングが終わると、トレース図が表示される（図２）。

 

図２

 

図２のWindowの右上角のX印をクリックして閉じると、図３のグラフが表示される。

 

図３

 

水準の効果値（グラフではで表されている）の事後分布のグラフである。水準１と水準４の分布がほぼ重なっており、水準３の事後分布がその左側に、水準２の事後分布がその右側に位置している。すなっわち、水準３＜（水準１、水準４）＜水準２となっている。

 

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

 

 

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

 

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.

 

 

 

リスト２　リスト１のStanスクリプトを用いるPythonスクリプト（OneFacRepRndm.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='OneFacRepRndm.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