第７章 関係を表すパス図
共分散構造分析
パス図の分析法として共分散構造分析（analysis of covariance structures, general structural equation technique）と呼ばれている方法がある。＜岡本安晴「心理学データ分析と測定」勁草書房＞の図7.1.4の例に対してサンプルスクリプトを作成してみた。図7.1.4を少し簡略化して図7.1.4（改）に示す。
図7.1.4（改）
図7.1.4（改）のモデルをリスト１に示すStanスクリプトで表してみた。
リスト１ 図7.1.4（改）のStanスクリプト（ファイル名SubjModel.stan）
data {
int
N;
real Jpn[N];
real Eng[N];
real Hist[N];
real Math[N];
real Phys[N];
real Chem[N];
}
parameters {
real b1;
real b2;
real FG[N];
real FL[N];
real FS[N];
real a[6];
}
transformed parameters {
real f_Jpn[N];
real f_Eng[N];
real f_Hist[N];
real f_Math[N];
real f_Phys[N];
real f_Chem[N];
for (i
in 1:N){
f_Jpn[i] = a[1] * FL[i];
f_Eng[i] = a[2] * FL[i];
f_Hist[i] = a[3] * FL[i];
f_Math[i] = a[4] * FS[i];
f_Phys[i] = a[5] * FS[i];
f_Chem[i] = a[6] * FS[i];
}
}
model {
b1 ~ uniform(0.0, 1.0);
b2 ~ uniform(-1.0, 1.0);
for (j in 1:6){
a[j]
~ uniform(0.0, 1.0);
}
for (i
in 1:N){
FG[i] ~ normal(0, 1);
FL[i] ~ normal(b1 * FG[i], sqrt(1 - b1 * b1));
FS[i] ~ normal(b2 * FG[i], sqrt(1 - b2 * b2));
Jpn[i] ~ normal(f_Jpn[i], sqrt(1
- a[1]*a[1]));
Eng[i] ~ normal(f_Eng[i], sqrt(1
- a[2]*a[2]));
Hist[i] ~ normal(f_Hist[i], sqrt(1
- a[3]*a[3]));
Math[i] ~ normal(f_Math[i], sqrt(1 - a[4]*a[4]));
Phys[i] ~ normal(f_Phys[i], sqrt(1
- a[5]*a[5]));
Chem[i] ~ normal(f_Chem[i], sqrt(1
- a[6]*a[6]));
}
}
リスト１のStanスクリプトを使って、図7.1.4（改）に対応するデータ（下図１、ファイル名Data.xlsx）を読み込み分析するPythonスクリプトを、本ウェブサイトの後半に示したリスト２のように用意した。これらのファイルは圧縮ファイルfiles7.zipとして用意したので、ダウンロードすれば自由に利用できる。
・
・
・
図１ データData.xlsx
リスト２のスクリプトのデータ読み込み部は以下のようになっている。
flnm = input('Data File Name (*.xlsx) = ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
上のスクリプトを実行すると、入力データファイル名を聞いてくる（図２）。
図２
入力データファイル名（図２の場合は、図１のデータファイル名）を入力すると、設定された名前のExcelのデータファイルが読み込まれ、PandasのDataFrame型のオブジェクトに格納されてdataによって指し示される。このDataFrame型のオブジェクトから、以下のスクリプトによってリスト１のStanスクリプト用のデータが用意される。
VarLabels = data.columns[1:].values;
print(VarLabels)
X = data.values[:,
1:]
X = ss.zscore(X)
X = X.T
Jpn = X[0]
Eng = X[1]
Hist = X[2]
Math = X[3]
Phys = X[4]
Chem = X[5]
N = len(Jpn)
データ値Xは標準得点化されている。Python用のStan（PyStan）の使い方は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第10章などを参照されたい。
Stanによるサンプリングが終了すると、まず、係数b1とb2の事後分布が表示される（図３）。
図３
図３のウィンドウを閉じると、係数a1からa6までの事後分布が表示される（図４）。
図４
図４を閉じると、因子得点fJpnの事後分布の平均値とデータ値Jpnの標準化得点との散布図が表示される（図５）。グラフの上に表示されている数値は、a1の事後分布の平均値である。
図５
図５を閉じると、因子得点fEngの事後分布の平均値とデータ値Engの標準化得点との散布図が表示される（図６）。グラフの上部の数値はa2の事後分布の平均値である。
図６
以下同様に、ウィンドウを閉じると次の散布図が表示されていく。最後に、図７の散布図が表示される。
図７
図７のウィンドウを閉じると、プログラムの終了である（図８）。
図８
コンソール画面には、各係数の平均値が出力されている。
リスト２ リスト１のStanスクリプトによる分析を行うPythonスクリプト（圧縮ファイルfiles7.zipに含まれている）
import pandas as pd
import numpy
as np
import matplotlib.pyplot
as plt
import seaborn
as sb
import scipy.stats
as ss
import pystan
import pickle
flnm = input('Data File Name (*.xlsx) = ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
VarLabels = data.columns[1:].values;
print(VarLabels)
X = data.values[:,
1:]
X = ss.zscore(X)
X = X.T
Jpn = X[0]
Eng = X[1]
Hist = X[2]
Math = X[3]
Phys = X[4]
Chem = X[5]
N = len(Jpn)
Data = {'N': N, 'Jpn':
Jpn, 'Eng': Eng, 'Hist': Hist,
'Math': Math, 'Phys': Phys,
'Chem': Chem}
def f_init():
return dict(b1
= 0.9, b2 = 0.9, a = [0.9, 0.9, 0.9, 0.9, 0.9, 0.9])
sm = pystan.StanModel(file
= 'SubjModel.stan')
fit = sm.sampling(data
= Data, init = f_init, n_jobs = 1)
print(fit)
b1_smpls = fit['b1']
b2_smpls = fit['b2']
b1 = fit['b1'].mean()
b2 = fit['b2'].mean()
print('\nb1 = {0:.5}'.format(b1))
print('b2 = {0:.5}'.format(b2))
a = fit['a'].mean(axis = 0)
for i in
range(6):
print('a[{0:}] =
{1:.5}'.format(i+1, a[i]))
sb.kdeplot(b1_smpls, label = 'b1')
sb.kdeplot(b2_smpls, label = 'b2')
plt.title('Means:\nb1 = {0:.3}, b2 = {1:.3}'.format(b1, b2), fontsize = 14)
plt.legend()
plt.show()
a_smpls = fit['a'].T
for i in
range(6):
sb.kdeplot(a_smpls[i], label =
'a{}'.format(i+1))
plt.legend()
plt.title(('Means: a1 = {0:.3}, a2 = {1:.3}\n' +
'a3 = {2:.3}, a4 =
{3:.3}, a5 = {4:.3}, a6 = {5:.3}').format(
a[0], a[1], a[2], a[3], a[4], a[5]), fontsize
= 14)
plt.show()
fJpn = fit['f_Jpn'].mean(axis
= 0)
plt.plot(fJpn, Jpn, 'bo')
plt.ylabel('Jpn(std)')
plt.xlabel('fJpn')
plt.title('a1 = {0:.5}'.format(a[0]), fontsize = 14)
plt.show()
fEng = fit['f_Eng'].mean(axis
= 0)
plt.plot(fEng, Eng, 'bo')
plt.xlabel('fEng')
plt.ylabel('Eng(std)')
plt.title('a2 = {0:.5}'.format(a[1]), fontsize = 14)
plt.show()
fHist = fit['f_Hist'].mean(axis
= 0)
plt.plot(fHist, Hist, 'bo')
plt.xlabel('fHist')
plt.ylabel('Hist(std)')
plt.title('a3 = {0:.5}'.format(a[2]), fontsize = 14)
plt.show()
fMath = fit['f_Math'].mean(axis
= 0)
plt.plot(fMath,
Math, 'bo')
plt.xlabel('fMath')
plt.ylabel('Math(std)')
plt.title('a4 = {0:.5}'.format(a[3]), fontsize = 14)
plt.show()
fPhys = fit['f_Phys'].mean(axis
= 0)
plt.plot(fPhys, Phys, 'bo')
plt.xlabel('fPhys')
plt.ylabel('Phys(std)')
plt.title('a5 = {0:.5}'.format(a[4]), fontsize = 14)
plt.show()
fChem = fit['f_Chem'].mean(axis
= 0)
plt.plot(fChem, Chem, 'bo')
plt.xlabel('fChem')
plt.ylabel('Chem(std)')
plt.title('a6 = {0:.5}'.format(a[5]), fontsize = 14)
plt.show()