因子分析: 主因子法
Python版
因子分析は、因子という潜在変数
の１次式によってデータ（変量）
の変動を説明するモデルであり、次式（１）
（１） 
によって表わされ、多くの解法が提案されている。よく用いられるものの１つとして、主因子法がある。主因子法では、標準化された変量に対して、因子と残差項
の独立性の仮定の下で、以下の繰り返しにより因子パターン
（因子負荷量）が求められる。
（Ａ）誤差項の分散を対角要素とする行列を
とおく。初期値は、
である。
（Ｂ）相関行列
を対して
の固有分解
から、固有値の大きいもの
個を選び
と
をおく。
（C）
を
と更新する。
の値が収束していると判定されれば、終了。収束していると判定されなければ、ステップ（B）に戻る。
詳しくは、岡本「心理学データ分析と測定」, pp.86-89を参照されたい。
主因子法のスクリプトをリスト１のように作成してみた。ファイルは、まとめて圧縮ファイルfapfpyfiles.zipとした。
入力データファイルは、Excelファイルとして図１a,bのように用意する。図１a,bのデータは、岡本（2014）、表5.1.1を入力したものである。
図１a
図１b
１行目に変数名を設定する。１列目は、ケースIDの番号を設定する。この値は、分析には用いられない。２列目以降に、データ値を設定する。
リスト１のスクリプトを実行すると、入力データファイル名を聞いてくる。次の例では、図１のデータファイル名Data6Subj100.xlsxが設定されている。入力データファイルは、スクリプトファイルと同じフォルダ内にあるものとする。
Input file name(*.xlsx) =
Data6Subj100.xlsx
Output file name(*.txt) = Results.txt
入力データファイル名を設定すると、出力ファイル名の設定が求められる。適当なテキストファイル名を設定する。上の例では、Results.txtが設定されている。出力ファイルは、スクリプトと同じフォルダ内に保存される。
出力テキストファイル名を設定すると、計算が始まり、固有値のスクリープロットが表示される（図２）。
図２
図２の場合、因子数は２が適当であると判断される。
図２のWindowを閉じると、因子数qの値の設定が「q = 」と求められる（図３）。
図３
図３では、２が設定されている。
因子数qを設定して、Enterキーを押すと、主因子法の計算が始まる。
因子数qの値によっては、ヘイウッドケースなど不適切な解が得られることがあるが、そのときは他の方法を試せばよい。直接最小2乗法による因子分析のウェブサイトも用意している。
主因子法の計算が終了すると、回転前の因子負荷量Aの値の布置が初期値として表示される（図４）。
図４
図４のWindowを閉じると、次に表示する布置の横軸（X軸）の表す因子の指定が求められる（図５）。
図５
因子として０以下の整数値を設定すると、初期値による布置の表示が終了する。
布置に表示する因子として０以下の整数値を設定して、初期値による布置の表示を終了すると、次にバリマックス回転による直交回転が行われる。
バリマックス回転は、次の基準を最小化するものとして、足立・山本（2024）、p.137のアルゴリズムにより行った（リスト２）。
バリマックス回転後の布置は、図６のように表示される。
図６
図６のWindowを閉じると、次の布置の表示に用いる横軸の表す因子の設定が求められる。
因子として０以下の整数値を設定すると、バリマックス回転後の布置の表示は終了する。
次に、斜交回転としてコーティミン回転が行われる。
コーティミン基準は
で与えられ、足立・山本（2024）、p.140のアルゴリズムを参考にしてリスト２に示すスクリプトを用意した。
斜交回転が終了すると、スクリプトの実行終了である（図７）。
図７
実行終了後、出力ファイルResults.txtを開くと、以下のようである。
Input fie = Data6Subj100.xlsx.
var_names =
['Jpn' 'Eng'
'Hist' 'Math' 'Phys' 'Chem']
1: [48 46 60 47 68 44]
2: [54 49 63 61 49 45]
3: [48 55 47 49 49 47]
・
・
・
98: [84 84 96 69 60 69]
99: [39 51 39 59 60 65]
100: [47 45 60 68 72 70]
mean sd
Jpn
58.750 15.213
Eng
59.310 15.591
Hist
59.310 15.059
Math
60.520 13.056
Phys
60.090 13.287
Chem
58.940 13.981
R =
Jpn Eng Hist Math Phys Chem
Jpn
1.00000 0.84066 0.87327 0.23810 0.19741 0.29807
Eng 0.84066 1.00000 0.82001 0.20751 0.12885 0.24563
Hist 0.87327 0.82001 1.00000 0.27644 0.20387 0.24637
Math 0.23810 0.20751 0.27644 1.00000 0.60998 0.67641
Phys 0.19741 0.12885 0.20387 0.60998 1.00000 0.84128
Chem 0.29807 0.24563 0.24637 0.67641 0.84128 1.00000
Eigen values of R
Cum.Sum Cum.Prop(%)
1
3.25276
3.25276
54.21
2
1.86791
5.12067
85.34
3
0.42494
5.54561
92.43
4
0.20646
5.75207
95.87
5
0.14753
5.89960
98.33
6 0.10040
6.00000
100.00
Before rotation
Factor-1 Factor-2 Communalitites
Jpn
-0.82978
-0.45147
0.89236
Eng
-0.75926
-0.46560
0.79326
Hist
-0.81036
-0.44279
0.85276
Math
-0.55047
0.44045
0.49701
Phys
-0.58987
0.65074
0.77141
Chem
-0.70145
0.65671
0.92330
Varimax rotation
Factor-1 Factor-2 Communalities
Jpn
-0.93379 0.14280
0.89236
Hist
-0.91308 0.13796
0.85276
Eng
-0.88621 0.08887
0.79326
Chem
-0.16094 0.94731
0.92330
Phys
-0.07573 0.87503
0.77141
Math
-0.17163 0.68378
0.49701
Contri
2.552
2.178
Quartimin rotation
Factor-1 Factor-2
Jpn
-0.93987
0.01669
Hist
-0.91928
0.01459
Eng
-0.89909 -0.03229
Chem
-0.02114
0.95479
Phys
0.05557
0.89213
Math
-0.07238
0.68142
Structure matrix
Factor-1 Factor-2
Jpn
-0.94451
0.27816
Hist
-0.92334
0.27034
Eng
-0.89011
0.21784
Chem
-0.28677
0.96067
Phys
-0.19263
0.87667
Math
-0.26196
0.70156
足立浩平・山本倫生（2024）主成分分析と因子分析：特異値分解を出発点として．共立出版
岡本安晴（2019）いまさら聞けないPythonでデータ分析：多変量解析、ベイズ統計分析．丸善出版
岡本安晴（2014）心理学データ分析と測定：データの見方と心の測り方．勁草書房
import pandas as pd
import numpy
as np
from mymodule
import * # arrng_pttn,
draw_map,
from rotate import * # Varimax_Orthomax,
Oblique_Quartimin
fin_nm = input('Input file name(*.xlsx) = ')
df_data = pd.read_excel(fin_nm)
# Excelファイルの読み込み
fout_nm = input('Output file name(*.txt) = ')
fout = open(fout_nm,
'w')
# 出力用テキストファイルのオープン
fout.write(f'Input fie
= {fin_nm}.\n')
var_names = np.array(df_data.keys())
caseIDs = np.array(df_data[var_names[0]])
var_names = var_names[1:]
# 変数名リスト
print('var_names
=\n', var_names)
fout.write(f'var_names
=\n{var_names}\n\n')
print('n =', len(caseIDs))
X = np.array(df_data.values)[:,1:]
# 変数値の配列（行列）
for i, v in
zip(caseIDs, X):
print(i,
v)
fout.write(f'{i}: {v}\n')
means = np.mean(X,
axis=0)
sds = np.std(X,
axis=0)
print()
print(f'{"
":10s}{"mean":>10s}{"sd":>10s}')
fout.write(f'\n\n{"
":10s}{"mean":>10s}{"sd":>10s}\n')
for i, v in enumerate(var_names):
print(f'{v:<10s}{means[i]:>10.3f}{sds[i]:>10.3f}')
fout.write(f'{v:<10s}{means[i]:>10.3f}{sds[i]:>10.3f}\n')
R = np.corrcoef(X,
rowvar=False) # 相関行列
print()
print('相関行列')
fout.write('\nR =\n')
print(f'{" ":>10s}',
end='')
fout.write(f'{" ":>10s}')
for v in var_names:
print(f'{v:>10s}',
end='')
fout.write(f'{v:>10s}')
fout.write('\n')
print()
for i, s in
enumerate(var_names):
print(f'{s:>10s}',
end='')
fout.write(f'{s:>10s}')
for j in range(len(var_names)):
print(f'{R[i][j]:>10.5f}', end='')
fout.write(f'{R[i][j]:>10.5f}')
print()
fout.write('\n')
Lmbd, V = np.linalg.eigh(R) # 相関行列の固有分解
idx = np.argsort(-Lmbd)
# 固有値の降順に並べ替えるためのインデックス
Lmbd = Lmbd[idx]
V = V.T[idx].T
n_vars = len(var_names)
cum_sum = np.cumsum(Lmbd)
print('\n相関行列の固有値')
fout.write('\n\nEigen
values of R\n')
print(f'{"Cum.Sum":>40s}{"Cum.Prop(%)":>15s}')
fout.write(f'{"Cum.Sum":>40s}{"Cum.Prop(%)":>15s}\n')
for i in
range(n_vars):
print(f'{i+1:>10d}{Lmbd[i]:>15.5f}', end ='')
print(f'{cum_sum[i]:>15.5f}{100*cum_sum[i]/cum_sum[n_vars-1]:>15.2f}')
fout.write(f'{i+1:>10d}{Lmbd[i]:>15.5f}')
fout.write(f'{cum_sum[i]:>15.5f}{100*cum_sum[i]/cum_sum[n_vars-1]:>15.2f}\n')
"""
固有値のスクリープロット
"""
plt.title('Scree plot', fontsize=16)
plt.plot(np.arange(n_vars)+1, Lmbd, ls='-')
plt.ylabel('Eigen values', fontsize=14)
plt.tight_layout()
plt.show()
q = int(input('q = ')) # 因子数の入力
if q <= 0:
q = 1
if q > n_vars:
q = n_vars
"""
主因子法
"""
Psi = np.zeros((n_vars))
RPsi = R - np.diag(Psi)
ckcnt = 0
while True:
Lmbd,
V = np.linalg.eigh(RPsi)
idx
= np.argsort(-Lmbd)
Lmbd
= Lmbd[idx]
V = V.T[idx].T
for ck in range(q):
if Lmbd[ck] < 1.0e-9:
q = ck
print(q)
break
A = V[:,:q] * (Lmbd[:q]**0.5)
R1 = A @ A.T
Psi1 = np.diag(np.identity(n_vars) - R1)
ckcnt
+= 1
ck_cnv
= np.sum(np.abs(Psi1 -
Psi))
if (ck_cnv
< 1.0e-7) or (ckcnt > 1000):
print('\nck_cnv =', ck_cnv,
' ckcnt
=', ckcnt)
print()
break
Psi = Psi1
RPsi
= R - np.diag(Psi)
"""
共通性
"""
communlty = np.sum(A
** 2, axis=1)
print('\n初期解')
fout.write('\n\nBefore
rotation\n')
print(f'{" ":10s}', end='')
fout.write(f'\n{" ":10s}')
for i in
range(q):
tstr
= f"Factor-{i+1}"
print(f'{tstr:>15s}',
end='')
fout.write(f'{tstr:>15s}')
print(f'{"Communalities":>15s}')
fout.write(f'{"Communalitites":>15s}\n')
for i in
range(n_vars):
print(f'{var_names[i]:<10s}', end='')
fout.write(f'{var_names[i]:<10s}')
for j in range(q):
print(f'{A[i][j]:>15.5f}', end='')
fout.write(f'{A[i][j]:>15.5f}')
print(f'{communlty[i]:>15.5f}')
fout.write(f'{communlty[i]:>15.5f}\n')
"""
初期値のマップ
"""
dim1 = 1
dim2 = 2 if q >= 2 else 1
while True:
draw_map(var_names, A.T[dim1-1], A.T[dim2-1],
'Before Rotation',
f'Factor-{dim1}', f'Factor-{dim2}')
dim1 = int(input('dim_x = '))
if dim1 <= 0:
break
if dim1 > q:
dim1
= q
dim2 = int(input('dim-y =
'))
if dim2 <= 0:
break
if dim2 > q:
dim2
= q
"""
バリマックス回転
"""
var_rot = Varimax_Orthomax(A)
A_vmax = var_rot.rotate()
print('\nバリマックス回転')
fout.write('\n\nVarimax
rotation\n')
#for i in
range(n_vars):
# print(var_names[i], ' ', A_vmax[i])
v_names, pttn = arrng_pttn(var_names, A_vmax)
dim1 = 1
dim2 = 2 if q >= 2 else 1
while True:
draw_map(var_names, A_vmax.T[dim1-1], A_vmax.T[dim2-1],
'Varimax rotation', f'dim-{dim1}', f'dim-{dim2}')
dim1 = int(input('dim_x = '))
if dim1 <= 0:
break
if dim1 > q:
dim1
= q
dim2 = int(input('dim-y =
'))
if dim2 <= 0:
break
if dim2 > q:
dim2
= q
comm_vmax = np.sum(pttn**2, axis=1)
contrib = np.sum(pttn**2, axis=0)
print('\nバリマックス回転後の因子パターン')
print(f'{" ":10s}', end='')
fout.write(f'\n{" ":10s}')
for i in
range(q):
vs = f'Factor-{i+1}'
print(f'{vs:>10s}',
end='')
fout.write(f'{vs:>10s}')
print(f'{"Communalities":>15s}')
fout.write(f'{"Communalities":>15s}\n')
for i in
range(n_vars):
print(f'{v_names[i]:<10s}', end='')
fout.write(f'{v_names[i]:<10s}')
for j in range(q):
print(f'{pttn[i][j]:>10.5f}',
end='')
fout.write(f'{pttn[i][j]:>10.5f}')
print(f'{comm_vmax[i]:>15.5f}')
fout.write(f'{comm_vmax[i]:>15.5f}\n')
print(f'{"Contri":<10s}',
end='')
fout.write(f'{"Contri":<10s}')
for j in range(q):
print(f'{contrib[j]:>10.3f}',
end='')
fout.write(f'{contrib[j]:>10.3f}')
print()
print()
fout.write('\n\n')
"""
コーティミン回転
"""
q_rot = Oblique_Quartimin(A)
A_qrot, Phi = q_rot.rotate()
v_names, pttn = arrng_pttn(var_names, A_qrot)
print('\nコーティミン回転後の因子パターン\n')
fout.write('\nQuartimin
rotation\n')
print(f'{" ":10s}', end='')
fout.write(f'\n{" ":10s}')
for i in
range(q):
vs = f'Factor-{i+1}'
print(f'{vs:>15s}',
end='')
fout.write(f'{vs:>15s}')
print()
fout.write('\n')
for i in
range(n_vars):
print(f'{v_names[i]:<10s}', end='')
fout.write(f'{v_names[i]:<10s}')
for j in range(q):
print(f'{pttn[i][j]:>15.5f}',
end='')
fout.write(f'{pttn[i][j]:>15.5f}')
print()
fout.write('\n')
"""
構造行列
"""
Strctr = A_qrot @
Phi
v_names, pttn = arrng_pttn(var_names, Strctr)
print('\n構造行列')
fout.write('\n\nStructure
matrix\n')
print(f'{" ":10s}', end='')
fout.write(f'{" ":10s}')
for i in
range(q):
vs = f'Factor-{i+1}'
print(f'{vs:>15s}',
end='')
fout.write(f'{vs:>15s}')
print()
fout.write('\n')
for i in
range(n_vars):
print(f'{v_names[i]:<10s}', end='')
fout.write(f'{v_names[i]:<10s}')
for j in range(q):
print(f'{pttn[i][j]:>15.5f}',
end='')
fout.write(f'{pttn[i][j]:>15.5f}')
print()
fout.write('\n')
fout.close()
print(f'\n{fout_nm}
is saved.\n')
"""
Y.Okamoto, 2025.06
参考文献
足立浩平・山本倫生（２０２４）
主成分分析と因子分析、共立出版
"""
import numpy
as np
import scipy.optimize
as so
class Varimax_Orthomax:
def __init__(self,
C):
"""
C: 負荷行列
"""
self.C = C
def rotate(self):
"""
足立・山本（2014）、p.137
"""
p = len(self.C)
Mw =
np.identity(p) - np.ones((p,p))/p
# w==1
Gmma = self.C
while True:
Gmma2 = Gmma * Gmma
S_tilde = Gmma * (Mw
@ Gmma2)
U, Delta, Vh = np.linalg.svd(self.C.T @ S_tilde)
T = U @ Vh
Gmma_new = self.C @
T
ck = np.sum((Gmma - Gmma_new)**2)
Gmma = Gmma_new
#print(ck)
if ck < 1.0e-6:
break
return Gmma # 回転後の行列
class Oblique_Quartimin():
"""
C: 負荷行列
"""
def __init__(self,
C):
self.C = C
def rotate(self):
"""
足立・山本（2014）、p.140
"""
p, m
= np.shape(self.C)
Gmma = self.C
Phi = np.identity(m)
T = np.identity(m)
def calc_tjk(gmma, phi, j, k):
class Gammas():
def __init__(self, C, Phi, j, k):
#self.c = Gmma
self.C = C
self.Phi = Phi
self.j = j
self.k = k
def f_gammas(self, t):
"""
Quartimin criterion
"""
c = self.C
Phi = self.Phi
p = np.shape(c)[0]
j = self.j
k = self.k
v = 0.0
for i in range(p):
g2_ij = (1 - 2*t*Phi[j][k] + (t**2)) * (c[i][j]**2)
g2_ik = (t**2) * (c[i][j]**2) + \
2 * t * c[i][j] * c[i][k]
+ c[i][k]**2
v += g2_ij * g2_ik
return v
calc_g = Gammas(gmma,
phi, j, k)
result = so.minimize_scalar(calc_g.f_gammas, [-10, 10])
#print(result)
return result.x
n_try = 0
while True:
Gmma0 = Gmma
for j in range(m):
for k in range(m):
if j != k:
# j < k を j != k に変更
tjk = calc_tjk(Gmma, Phi, j, k)
tjj = (1 - 2 * tjk *
Phi[j][k] + (tjk**2)) ** 0.5
T_jk = np.identity(m)
T_jk[j][j] = tjj
T_jk[j][k] = tjk
Gmma = Gmma @ T_jk
inv_T_jk = np.linalg.inv(T_jk)
Phi = inv_T_jk @ Phi @ (inv_T_jk.T)
T = T @ T_jk
n_try += 1
ck = np.sum((Gmma0 - Gmma) ** 2)
#print('ck =', ck, ' n_try =', n_try)
if (ck < 1.0e-6) or (n_try > 50):
break
else:
T = np.identity(m)
return Gmma, Phi # 回転後の負荷行列と共分散行列
リスト３ 変数の整列、布置の描画スクリプト（mymodule.py）
import copy
import numpy
as np
import matplotlib.pyplot
as plt
def arrng_pttn(names,
A):
v_nms
= copy.copy(names)
pttn
= copy.copy(A)
n = len(v_nms)
q = len(pttn[0])
vals
= np.empty(n)
for i
in range(n):
mx_p = np.argsort(-np.abs(pttn[i]))
vals[i] = np.abs(pttn[i][mx_p[0]])
* (0.1**mx_p[0])
idx
= np.argsort(-vals)
return v_nms[idx], pttn[idx]
def draw_map(names,
xcoord, ycoord, stitle, sxlabel, sylabel):
plt.figure(figsize=(5, 5))
plt.title(stitle, fontsize=14)
plt.xticks([-1,
0, 1])
plt.yticks([-1,
0, 1])
plt.xlabel(sxlabel, fontsize=12)
plt.ylabel(sylabel, fontsize=12)
plt.plot([-1,
1], [0, 0], c='k', lw=1)
plt.plot([0,
0], [-1, 1], c='k', lw=1)
n = len(names)
for i
in range(n):
plt.plot(xcoord[i], ycoord[i],
'bo', mfc='none')
for i
in range(n):
plt.text(xcoord[i], ycoord[i],
names[i], c='k')
plt.tight_layout()
plt.show()