Up

因子分析

直接最小２乗法

　因子分析は、因子という潜在変数の１次式によってデータ（変量）の変動を説明するモデルであり、次式（１）

（１）

によって表わされる。多くの解法が提案されているが、最小２乗法の基準により、直接、因子負荷行列を求める方法（岡本、2014、pp.95-97）を以下に試みた。

直接最小２乗法は、共通性（communality）が１より大きくなるというヘイウッドケースが無く、因子得点も直接求められる。

一般に、行列の、ランクの行列による最小2乗近似は、特異値分解の特異値が大きいものから個とったもので与えられる。これはエッカート・ヤングの定理（Eckart and Young, 1936）と呼ばれているが、それより以前に発表されている(e.g. Schmidt, 1906)とGifi (1990), p.518は説明している。

いま、式（１）の左辺データは標準得点化されているとし、その特異分解を

とおく。このとき、ランクの最小２乗近似は、特異値の大きいものから個を選んで

で与えられる。式（２）と式（１）を比べ、因子が平均0、分散1の標準化得点であることに注意すると

で与えられることがわかる。ここで、はデータ数である。また、より、である。

主成分分析でも形式は同じ式が使われるが、主成分分析の場合は、データの正射影による圧縮であるのに対して、上記（３）はデータに対する因子モデル（１）による因子得点と因子負荷の推定値である。データを効率よく圧縮したもの（主成分分析）とデータをよく表す推定値（因子分析）が数値として一致することは自然なことであると言えよう。

式（１）より、次式を得る。

両辺の（h,h）成分に注目して、次式を得る。

ここで、

より、

が導かれる。すなわち、共通性は１を超えることはない。

正則行列に対して、次式が成り立つ。

このは、回転と呼ばれている。

直接最小2乗法のスクリプトをリスト１のように作成してみた。ファイルは、まとめて圧縮ファイルfa_direct_files.zipとした。

入力データファイルは、Excelファイルとして図１a,bのように用意する。図１a,bのデータは、岡本（2014）、表5.1.1を入力したものである。

図１a

図１b

１行目に変数名を設定する。１列目は、ケースIDの番号を設定する。ケースIDの番号は、整数であれば任意であり、昇順でなくてもランダムな順序でもよいが、因子得点の出力において用いられるが、分析には用いられない。２列目以降に、データ値を設定する。

リスト１のスクリプトを実行すると、入力データファイル名を聞いてくる。次の例では、図１のデータファイル名Data6Subj100.xlsxが設定されている。入力データファイルは、スクリプトファイルと同じフォルダ内にあるものとする。

Input file name(*.xlsx) = Data6Subj100.xlsx

Output file name(*.txt) = Results.txt

入力データファイル名を設定すると、出力ファイル名の設定が求められる。適当なテキストファイル名を設定する。上の例では、Results.txtが設定されている。出力ファイルは、スクリプトと同じフォルダ内に保存される。

出力テキストファイル名を設定すると、計算が始まり、相関行列の固有値のスクリープロットが表示される（図２）。

図２

図２の場合、因子数は２が適当であると判断される。

図２のWindowを閉じると、因子数qの値の設定が「q = 」と求められる（図３）。

図３

図３では、２が設定されている。

因子数qを設定して、Enterキーを押すと、直接最小2乗法による計算が始まる。

計算が終了すると、回転前の因子負荷量Aの値の布置が初期値として表示される（図４）。

図４

図４のWindowを閉じると、次に表示する布置の横軸（X軸）の表す因子の指定が求められる（図５）。

図５

因子として０以下の整数値を設定すると、初期値による布置の表示が終了する。

布置に表示する因子として０以下の整数値を設定して、初期値による布置の表示を終了すると、次にバリマックス回転による直交回転が行われる。

バリマックス回転は、次の基準を最小化するものとして、足立・山本（2024）、p.137のアルゴリズムにより行った（リスト２）。

バリマックス回転後の布置は、図６のように表示される。

図６

図６のWindowを閉じると、次の布置の表示に用いる横軸の表す因子の設定が求められる。

因子として０以下の整数値を設定すると、バリマックス回転後の布置の表示は終了する。

次に、斜交回転としてコーティミン回転が行われる。

コーティミン基準は

で与えられ、足立・山本（2024）、p.140のアルゴリズムを参考にしてリスト２に示すスクリプトを用意した。

斜交回転が終了すると、スクリプトの実行終了である（図７）。

図７

回転前の初期解の因子得点が、fi.csvに、バリマックス回転後の因子得点がfv.csvに、コーティミン回転後の因子得点がfq.csvに保存されていることが示されている。これらのCSVファイルは、Excelで開くことができる。Excelで開かれたCSVファイルは、拡張子を.xlsxに変更してExcelファイルとして保存することができる。

実行終了後、出力ファイル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.82210 0.48694 0.91297

Eng 0.77575 0.52727 0.87981

Hist 0.81441 0.48464 0.89813

Math 0.64182 -0.53731 0.70064

Phys 0.62907 -0.67123 0.84628

Chem 0.71006 -0.61535 0.88285

Varimax rotation

Factor-1 Factor-2 Communalities

Jpn 0.94506 -0.14083 0.91297

Hist 0.93763 -0.13776 0.89813

Eng 0.93454 -0.08029 0.87981

Chem 0.16272 -0.92540 0.88285

Phys 0.06463 -0.91766 0.84628

Math 0.15901 -0.82180 0.70064

Contri 2.702 2.419

Quartimin rotation

Correlation between factors...

Factor-1 Factor-2

comp.1 1.00000 -0.26029

comp.2 -0.26029 1.00000

Factor loadings

Factor-1 Factor-2

Jpn 0.94954 -0.02195

Eng 0.94735 0.03888

Hist 0.94236 -0.01976

Phys -0.06404 -0.93452

Chem 0.03549 -0.92974

Math 0.04640 -0.82376

Structure matrix

Factor-1 Factor-2

Jpn 0.95526 -0.26911

Hist 0.94751 -0.26505

Eng 0.93723 -0.20770

Chem 0.27749 -0.93897

Phys 0.17921 -0.91785

Math 0.26082 -0.83584

参考文献

足立浩平・山本倫生（2024）主成分分析と因子分析：特異値分解を出発点として．共立出版

Gifi, A. (1990). Nonlinear multivariate analysis. John Wiley & Sons.

岡本安晴（2019）いまさら聞けないPythonでデータ分析：多変量解析、ベイズ統計分析．丸善出版

岡本安晴（2014）心理学データ分析と測定：データの見方と心の測り方．勁草書房

リスト１　直接最小2乗法のスクリプト（fa_direct_main.py）

import scipy.stats as ss

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')

N = len(caseIDs)

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')

X = ss.zscore(X, axis=0) # 標準得点化

print(np.mean(X, axis=0))

print(np.std(X, axis=0))

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]

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

"""

直接最小2乗法

"""

U, Lmbd, Vh = np.linalg.svd(X)

print('Vh.T[:q] =\n', Vh.T[:q])

print('np.diag(Lmbd[:q]) =\n', np.diag(Lmbd[:q]))

A = Vh[:q].T @ np.diag(Lmbd[:q]) / (N**0.5)

"""

共通性

"""

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')

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')

リスト２　回転のスクリプト（rotate.py）

"""

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()

Home