因子分析
因子分析の簡単な方法として主因子法（岡本、2014、pp. 85-87）、残差に制約条件を設定しない方法として直接最小２乗法（岡本、2014、pp. 95-97）について考える。まず、主因子法を取り上げ、次に直接最小２乗法を扱う。
因子分析のモデルを以下のように設定する（岡本、2014、pp. 85-87）。
ここで、
は標準化されたデータ値である。
は多変量の潜在変数（因子）であり、
の各成分は互いに独立に平均０，標準偏差１の分布に従うものとする。残差
は、平均が０で、各成分はお互いに独立であり、とも独立な分布に従うものとする。
このとき、の相関行列は次式で与えられる。
式（１）より、
から
を次の繰り返し法によって求める。
上の繰り返し法（主因子法）に基づくスクリプトをリストP1に示す。
リストP1のスクリプトを用いて主因子法による因子分析を行うプログラムをリストP2に示す。
スクリプトファイルとサンプルデータファイルは、まとめて圧縮ファイルFAPrg.zipとした。ダウンロード展開すれば、自由に利用できる。
リストP2の実行においては、Excel形式のファイルをPandasライブラリを用いて読み込むときにモジュールopenpyxlが必要である。
モジュールopenpyxlは、次のコマンドでインストールできる。
pip install
openpyxl
リスト２のファイルを次のコマンド
PS ******\FAPrg> python FA.py
で実行すると、入力データファイル名を聞いてくる。
Input Excel
data file (*.xlsx) =
入力データファイルは、Excelファイルとして用意する（図１）。
図１
第１行目に変数名が入力されている。データ値は２行目から設定する。１列目は、データの識別に用いるもので適当に設定する。２列目からが因子分析の対象となる。
データは、値の設定などが終われば保存してからプログラムを実行する。保存先は、リスト２のスクリプトファイルと同じフォルダにしておくと読み込みファイル名の指定が簡単である。
図１のファイルDataSubjects.xlsxを指定するときは、以下のようになる。
Input Excel
data file (*.xlsx) = DataSubjects.xlsx
上のファイル名の指定では、ファイルはスクリプトファイルと同じフォルダにあるとして処理される。
入力データファイルを指定すると、次に出力ファイル名の設定が求められる。出力ファイル名は任意である。次の例では、出力ファイル名Results.txtが設定されている。
Input Excel data file (*.xlsx) =
DataSubjects.xlsx
Output text data file (*.txt) =
Results.txt
出力ファイル名を設定すると計算が始まり、固有値が端末画面に出力されるとともに、スクリープロットが表示される（図２）。
図２
因子数は２個が適当と思われる。
図２のフォームを閉じると、因子数を聞いてくる。
Lambda cum.sqr
%
Lambda-1 =
3.36473 3.36473 56.08
Lambda-2 =
1.97245 5.33718 88.95
Lambda-3 =
0.28630 5.62349 93.72
Lambda-4 =
0.21584 5.83932 97.32
Lambda-5 =
0.09666 5.93599 98.93
Lambda-6 =
0.06401 6.00000 100.00
Number of factors =
２を設定してみる。
Number of
factors = 2
因子数を設定すると、パターン行列の描画における横軸（X軸）と縦軸（Y軸）の設定が求められる。
Set the factor number (from 1 to 2) to
plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 1
ydim = 2
Y軸の設定を行うと、図３のようにパターンが表示される。
図３
図３のフォームを閉じると、次の描画で用いるX軸の指定が求められる。
Set the factor number (from 1 to 2) to
plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 0
１より小さい数を設定すると、初期解でのパターンの表示が終わり、次に進む。上の例では、０が設定されているので、次のVarimax解に進む。。
Varimax回転後のパターンの表示におけるX軸とY軸の指定が求められる。
Set the factor number (from 1 to 2) to
plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 0
Set the factor number (from 1 to 2) to
plot the pattern for the Varimax criterion.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 1
ydim = 2
上の設定で、図４のパターンが表示される。
図４
図４のフォームを閉じると、次の描画におけるX軸の設定が求められる。１より小さい値を設定すると、Varimax回転後のパターンの描画は終了して次に進む。
Set the factor number (from 1 to 2) to
plot the pattern for the Varimax criterion.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 1
ydim = 2
Set the factor number (from 1 to 2) to
plot the pattern for the Varimax criterion.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 0
Results.txt was saved.
PS ******FAPrg>
上の例では、０を設定している。Varimax回転の次はPromax回転が行われるが、斜交回転であるので直交座標によるパターンの描画表示は行われない。Promax回転の結果は出力ファイルに保存される。
Promax回転が終わると、プログラムの終了である。上の端末出力例では、プログラムの終了時に出力ファイルがResults.txtであることが表示されている。。
上の場合の出力ファイルResults.txtの内容は以下のようになっている。
Input data file = DataSubjects.xlsx
R =
Jpn
Eng
Hist
Math Phys Chem
Jpn
1.00000 0.82892 0.89625 0.37343 0.00028 0.14148
Eng
0.82892 1.00000 0.80985 0.39106 0.18418 0.34114
Hist 0.89625 0.80985 1.00000 0.36743 0.12012 0.19902
Math 0.37343 0.39106 0.36743 1.00000 0.70281 0.83197
Phys 0.00028 0.18418 0.12012 0.70281 1.00000 0.83240
Chem 0.14148 0.34114 0.19902 0.83197 0.83240 1.00000
Lambda...
Lambd cum.sqr
%
Lambda-1 =
3.36473 3.36473 56.08
Lambda-2 =
1.97245 5.33718 88.95
Lambda-3 =
0.28630 5.62349 93.72
Lambda-4 =
0.21584 5.83932 97.32
Lambda-5 =
0.09666 5.93599 98.93
Lambda-6 =
0.06401 6.00000 100.00
Pattern matrix for the initial factors
fctr.1
fctr.2
communality
Jpn -0.76881 -0.61370
0.96770
Eng -0.77771 -0.38441
0.75261
Hist -0.76900 -0.51095
0.85243
Math -0.76204 0.41462
0.75261
Phys -0.55778 0.65470
0.73976
Chem -0.72637 0.66435
0.96897
Varimax rotation was applied...
Pattern for the Varimax rotated
factors =
fctr.1 fctr.2
Jpn -0.98298 0.03795
Eng -0.83987 0.21731
Hist -0.91597 0.11584
Math -0.30581 0.81184
Phys 0.00569 0.86007
Chem -0.11560 0.97755
Sorted Pattern for Varimax rotation =
fctr.1 fctr.2
Jpn -0.98298 0.03795
Hist -0.91597 0.11584
Eng -0.83987 0.21731
Chem -0.11560 0.97755
Phys 0.00569 0.86007
Math -0.30581 0.81184
Promax rotation was applied...
Pattern for the Promax rotated factors
=
fctr.1 fctr.2
Jpn -1.00561 -0.09742
Eng -0.83239 0.10710
Hist -0.92538 -0.00793
Math -0.19685 0.79330
Phys 0.13056 0.88620
Chem 0.02283 0.99030
Correlation matrix of the factors...
fctr.1 fctr.2
fctr.1 1.00000 -0.27068
fctr.2 -0.27068 1.00000
Sorted Pattern for Promax rotation =
fctr.1 fctr.2
Jpn -1.00561 -0.09742
Hist -0.92538 -0.00793
Eng -0.83239 0.10710
Chem 0.02283 0.99030
Phys 0.13056 0.88620
Math -0.19685 0.79330
Structure Matrix
fctr.1 fctr.2
Jpn -0.97924 0.17478
Eng -0.86138 0.33242
Hist -0.92324 0.24256
Math -0.41158 0.84658
Phys -0.10932 0.85086
Chem -0.24523 0.98412
以下のリストP1とP2のスクリプトファイルとサンプルデータファイルは、他のスクリプトファイルとともにまとめて圧縮ファイルFAPrg.zipとした。ダウンロード展開すれば自由に利用できる。
リストP1 主因子法のスクリプトファイル（init_fa.py）
import numpy
import matplotlib.pyplot as plt
from math import *
import mymodule as mm
import copy
class InitFA:
"""
Calculation of the Initial Factors
"""
def __init__(self, data,
f_out):
self.data = data
self.f_out
= f_out
def calc_fa(self):
print('data =\n', self.data)
n_vars = len(self.data[0]) - 1
var_names = self.data[0][1:]
X =
self.data[1:][:,1:]
print('X = \n', X)
X =
numpy.array(X, dtype = 'float')
R =
numpy.corrcoef(X.T)
print('R =\n', R)
self.f_out.write("\n\nR = \n")
self.f_out.write("{0:10s}".format(" "))
for
i in range(n_vars):
self.f_out.write("{0:10s}".format(var_names[i]))
self.f_out.write('\n')
for
i in range(n_vars):
self.f_out.write("{0:10s}".format(var_names[i]))
for j in range(n_vars):
self.f_out.write("{0:<10.5f}".format(R[i, j]))
self.f_out.write('\n')
Lmbd2, V = numpy.linalg.eig(R)
Lmbd
= []
for
i in range(0, len(Lmbd2)):
if Lmbd2[i] > 0.0:
Lmbd.append(Lmbd2[i])
else:
Lmbd.append(0.0)
n_eign = len(Lmbd)
ListLV = []
for
i in range(0, n_eign):
ListLV.append([])
ListLV[i].append(Lmbd[i])
tempL = []
for j in range(0, n_vars):
tempL.append(V[j,i])
ListLV[i].append(tempL)
ListLV = sorted(ListLV, key = lambda x: -x[0])
self.f_out.write("\nLambda...\n")
sum_L2 = 0.0
for
i in range(0, n_eign):
sum_L2 += ListLV[i][0]
cum_sum = 0.0
self.f_out.write("%30s"%"Lambd ")
self.f_out.write("%10s"%"cum.sqr")
self.f_out.write("%10s\n"%"%")
print("\n")
print(("%30s"%"Lambda") +
("%10s"%"cum.sqr") + ("%10s"%"%") +
"\n")
for
i in range(0, n_eign):
self.f_out.write("%20s"%("Lambda-%d = "%(i + 1)))
self.f_out.write("%10.5f"%ListLV[i][0])
cum_sum += ListLV[i][0]
self.f_out.write("%10.5f"%cum_sum)
self.f_out.write("%10.2f\n"%(100.0 * cum_sum / sum_L2))
print("%20s"%("Lambda-%d
= "%(i + 1)), end="")
print("%10.5f"%ListLV[i][0], end = "")
print("%10.5f"%cum_sum, end = "")
print("%10.2f"%(100.0 * cum_sum / sum_L2))
plt.title('Scree Plot', fontsize = 20)
plt.plot([v + 1 for v in range(n_eign)],
[v[0] for v in ListLV])
plt.plot([v + 1 for v in range(n_eign)],
[v[0] for v in ListLV], 'bo')
plt.ylabel('Eigen value', fontsize = 16)
plt.show()
s0 =
input("\nNumber of factors = ")
n_comp = int(s0)
if
(n_comp > n_vars):
print("Number of factors must be less than or equal to the number
of variables")
n_comp = n_vars
elif
(n_comp < 1):
print("Number of
factors should be positive integer.")
n_comp = 1
if
(n_comp > n_eign):
print("Number of factors should be less than or equal to the rank
of R")
n_comp = n_eign
Pat
= None
aPsy
= []
aPsy1 = []
for
i in range(n_vars):
aPsy.append([])
aPsy1.append([])
for j in range(n_vars):
aPsy1[i].append(0.0)
aPsy[i].append(0.0)
Psy
= numpy.zeros((n_vars, n_vars)) #
numpy.matrix(aPsy)
Psy1
= numpy.zeros((n_vars, n_vars)) #
numpy.matrix(aPsy1)
n_iter = 0
sum0
= 1.0e9
while True:
Psy1 = copy.deepcopy(Psy)
RPsy = R - Psy
vLmbd, Up = numpy.linalg.eigh(RPsy)
ck_idx = numpy.abs(vLmbd) > 1.0e-7
vLmbd = vLmbd[ck_idx]
Up = Up[:, ck_idx]
ck_idx = numpy.argsort(-vLmbd)
vLmbd = vLmbd[ck_idx]
Up = Up[:,ck_idx]
rnk = len(vLmbd)
if rnk > n_comp:
vLmbd = vLmbd[:n_comp]
Up = Up[:,:n_comp]
A = Up @ numpy.diag(vLmbd**0.5)
A = numpy.matrix(A)
AAp = A * numpy.transpose(A)
for i in range(n_vars):
Psy[i, i] = 1.0 - AAp[i, i]
sum = 0.0
for i in range(n_vars):
sum += abs(Psy[i, i] - Psy1[i, i])
#print('sum =', sum)
if (abs(sum - sum0) < 1.0e-7) or (sum > sum0):
Pat = A
break
sum0 = sum
n_iter += 1
#print('n_iter =', n_iter)
commu = []
for
i in range(n_vars):
sum = 0.0
for j in range(n_comp):
sum += A[i, j] ** 2.0
commu.append(sum)
self.f_out.write("\n\nPattern matrix for the initial
factors\n")
self.f_out.write("%10s"%" ")
for
i in range(0, n_comp):
self.f_out.write("%10s"%("fctr.%d"%(i
+ 1)))
self.f_out.write("%15s\n"%"communality")
for
i in range(0, n_vars):
self.f_out.write("%10s"%var_names[i])
for j in range(0, n_comp):
self.f_out.write("%10.5f"%Pat[i, j])
self.f_out.write("%15.5f\n"%commu[i])
return var_names, n_vars, n_comp, Pat
リストP2 リストP1のスクリプトを用いて因子分析をおこなうスクリプトファイル（FA.py）
import numpy
import pandas as pd
import math
import mymodule as mm
from init_fa import InitFA
from rotate import Varimax, Promax
from PlotComp import PlotPattern
from ArrangePat import ArrangePtn
"""
Principal Factor Analysis
"""
#
# Excelデータファイルの読み込み
#
s = input("Input Excel data file
(*.xlsx) = ")
xlsx_fl = pd.ExcelFile(s)
data_xlsx =
pd.read_excel(xlsx_fl, header =
None)
data = data_xlsx.values
#
# 計算結果出力用テキストファイル
#
s_out = input("Output text data
file (*.txt) = ")
f_out = open(s_out, "w")
f_out.write("Input data file =
" + s + "\n\n")
#
# 初期値
#
calc_init_fa = InitFA(data, f_out)
var_names, n_vars, n_comp, Pat =
calc_init_fa.calc_fa()
"""
var_names:
A list of variable names
n_vars: The
number of variables
n_com:
The number of factors
Pat:
The pattern matrix of the principal components
"""
PlotPtn = PlotPattern( Pat, var_names,
n_vars, n_comp, "Initial Pattern Matrix" )
while True:
print("\nSet the factor
number (from 1 to %d) to plot the pattern."%n_comp)
print("If you do not
want to plot the pattern, set the number less than 1.")
xdim = int(input("xdim
= "))
if xdim < 1:
break
ydim = int(input("ydim
= "))
PlotPtn.Plot(xdim, ydim)
rot = Varimax(Pat, n_vars, n_comp)
PatV = rot.rotateV()
f_out.write("\n\nVarimax rotation
was applied...\n\n")
f_out.write("\nPattern for the
Varimax rotated factors =\n")
f_out.write("%10s"%"
")
for j in range(0, n_comp):
f_out.write("%10s"%("fctr.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_vars):
f_out.write("%10s"%var_names[i])
for j in range(0, n_comp):
f_out.write("%10.5f"%PatV[i, j])
f_out.write("\n")
ArgPtn = ArrangePtn(var_names, PatV,
n_vars, n_comp)
PtnArgd = ArgPtn.arrange()
f_out.write("\n\nSorted Pattern
for Varimax rotation =\n")
f_out.write("%10s"%"
")
for j in range(0, n_comp):
f_out.write("%10s"%("fctr.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_vars):
f_out.write("%10s"%PtnArgd[i][0])
for j in range(0, n_comp):
f_out.write("%10.5f"%PtnArgd[i][1][j])
f_out.write("\n")
PlotPtn = PlotPattern( PatV,
var_names, n_vars, n_comp, "Pattern Matrix(Varimax)" )
while True:
print("\nSet the factor
number (from 1 to %d) to plot the pattern for the Varimax
criterion."%n_comp)
print("If you do not
want to plot the pattern, set the number less than 1.")
xdim = int(input("xdim
= "))
if xdim < 1:
break
ydim = int(input("ydim
= "))
PlotPtn.Plot(xdim, ydim)
obrot = Promax(PatV, n_vars, n_comp,
4.0)
Patob, Phi = obrot.rotateP();
f_out.write("\n\nPromax rotation
was applied...\n\n")
f_out.write("\nPattern for the
Promax rotated factors =\n")
f_out.write("%10s"%"
")
for j in range(0, n_comp):
f_out.write("%10s"%("fctr.%d"%(j
+ 1)))
f_out.write("\n")
for i in range(0, n_vars):
f_out.write("%10s"%var_names[i])
for j in range(0, n_comp):
f_out.write("%10.5f"%Patob[i, j])
f_out.write("\n")
f_out.write("\n\nCorrelation
matrix of the factors...\n")
f_out.write("%10s"%"
")
for j in range(0, n_comp):
f_out.write("%10s"%("fctr.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_comp):
f_out.write("%10s"%("fctr.%d"%(i + 1)))
for j in range(0, n_comp):
f_out.write("%10.5f"%Phi[i, j])
f_out.write("\n")
ArgPtn = ArrangePtn(var_names, Patob,
n_vars, n_comp)
PtnArgd = ArgPtn.arrange()
f_out.write("\n\nSorted Pattern
for Promax rotation =\n")
f_out.write('{0:10s}'.format(' '))
for i in range(n_comp):
str = 'fctr.{}'.format(i +
1)
f_out.write('{0:>10s}'.format(str))
f_out.write('\n')
for i in range(0, n_vars):
f_out.write("%10s"%PtnArgd[i][0])
for j in range(0, n_comp):
f_out.write("%10.5f"%PtnArgd[i][1][j])
f_out.write("\n")
S = Patob * Phi
f_out.write('\nStructure Matrix\n')
f_out.write('{0:10s}'.format(' '))
for i in range(n_comp):
str = 'fctr.{}'.format(i +
1)
f_out.write('{0:>10s}'.format(str))
f_out.write('\n')
for i in range(n_vars):
f_out.write('{0:>10s}'.format(var_names[i]))
for j in range(n_comp):
f_out.write('{0:10.5f}'.format(S[i, j]))
f_out.write('\n')
f_out.close()
print(s_out + " was saved.")
因子分析のモデルを次式（D1）のように設定する（de Leeuw,2004; Okamoto, 2005）。
ここで、
は標準得点化されたデータ行列であり、
は共通因子行列、はパターン行列、
は残差行列である。足立ら（2019）では、はさらに独自因子と残差に分けられているが、ここではモデル（D1）に基づいて分析を行う。モデル（D1）による分析を行うときに、残差の独立性など特別な条件は設定しない。また、ここでの分析法の場合、ヘイウッドケース（Heywood cases：共通性が１より大きくなること）が生じないことが証明されている（Okamoto, 2005, Appendix A）。
モデル（D1）のデータに対する誤差は次式（D2）で与えられる。
誤差（D2）を最小にする解は、Eckart-Young theoremにより、特異値分解で与えられる（Gifi, 1990, p.152, 518; 岡本、2014、p.95）。詳しい説明が、岡本（2009）、pp. 96-100にある。
特異値分解を
とおく。このとき、誤差（D2）を最小にする解は次式で与えられる。
ここで、
はデータ数であり、
、
、
は、
、
、
から特異値の大きいものから選んだものである。
また、
が成り立つので、
と
も解である。すなわち、解を回転することができる。
上で説明した式（D1）による因子分析法のスクリプトをリストD1のように用意した。スクリプトファイルとサンプルデータファイルはまとめて圧縮ファイルDFAPrg.zipとした。ダウンロード展開すれば自由に利用できる。
リストD1のスクリプトを次のコマンドで実行する。
PS ******\DFAPrg>
python DFA.py
実行すると、入力データファイル名を聞いてくる。
Input Excel
data file (*.xlsx) =
入力データファイルは図D1のようなExcelファイルとして用意する。
図D1
１行目に変数名を設定する。データは２行目以降に設定する。１列目は、データの識別用の値を適宜設定する。２列目以降に因子分析の対象となるデータを設定する。
データの設定が終われば保存する。保存先はプログラムのスクリプトファイルと同じフォルダにしておくと、プログラムからデータファイルの指定を行うとき簡単である。
図D1のファイルDataSubjects.xlsxを指定するときは、次のようにすればよい。
Input Excel
data file (*.xlsx) = DataSubjects.xlsx
入力データファイルを設定すると、次に出力テキストファイル名の設定が求められる。
Output Text
file(*.txt) =
適当なテキストファイル名を設定する。このテキストファイルに因子分析の結果が出力される。
Output Text
file(*.txt) = Results.txt
上の例では、Results.txtが設定されている。
出力ファイル名が設定されると計算が始まり、特異値が表示される。
Singular values =
[8.20333268 6.280843 2.39291599 2.07768208 1.39040963
1.13150139]
特異値のスクリープロットも表示される（図D2）
図D2
特異値の大きさから、因子数は２個が適切であると思われる。
図D2のフォームを閉じると、因子数を聞いてくる。次の例では、２と設定している。
Number of
factors = 2
因子数を設定すると、パターンの表示における横軸（X軸）の因子の設定が求められる。
Set the factor number (from 1 to 2) to
plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim =
X軸を設定すると、縦軸（Y軸）の設定が求められる。
xdim = 1
ydim = 2
上の例では、２が設定されている。
Y軸を設定すると、設定された因子でパターンが表示される（図D3）。
図D3
第１因子は、データ全体の変動に対応しており、第２因子は文系の科目と理系の科目の区別に対応していることが分かる。
図D3のフォームを閉じると、次のパターンの描画のための因子を聞いてくる。
X軸として１より小さい値を設定すると、初期値に対するパターンの描画は終了して、次にVarimax回転後のパターンの描画に進み、Varimax回転後のパターンの描画におけるX軸とY軸を聞いてくる。
次の例では、X軸に１，Y軸に２が設定されている。
Set the factor number (from 1 to 2) to
plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 1
ydim = 2
Set the factor number (from 1 to 2) to
plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 0
Set the factor number (from 1 to 2) to
plot the pattern for the Varimax criterion.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 1
ydim = 2
X軸とY軸の因子、上の例では第１因子と第２因子、を設定すると、設定した因子でパターンが描画される（図D4）。
図D4
第１因子が文系科目、第２因子が理系科目の因子であることが分かる。
図D4のフォームを閉じると、次の描画で用いるX軸とY軸の因子の設定が求められる。
X軸として１より小さい値を設定すると、Varimax回転後のパターンの描画が終了する。
Varimax回転後のパターンの描画が終了すると、次にPromax回転が行われるが、斜交回転であるのでパターンの描画は行われない。
Promax回転が終了すると、出力テキストファイルが保存されたことのメッセージが表示されて、プログラムの終了となる。
Set the factor number (from 1 to 2) to
plot the pattern for the Varimax criterion.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 0
Results.txt was saved.
PS ******\FA\DFAPrg>
上の出力例では、テキストファイルResults.txtが保存されたことが表示されている。
テキストファイルの内容は、以下のようである。
Input data = Results.txt
data =
[['ID' '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]
・
・
・
[18 62 73 63 60 64 67]
[19 66 60 61 64 43 63]
[20 70 66 71 50 55 52]]
Pattern matrix and communalities for
the initial factors
fctr.1
fctr.2
communality
Jpn 0.75184 0.61221 0.94007
Eng 0.81468 0.43846
0.85595
Hist 0.78026 0.54381
0.90454
Math 0.80603 -0.44283
0.84579
Phys 0.59683 -0.71379
0.86569
Chem 0.72185 -0.63568 0.92515
Varimax rotation was applied...
Pattern for the Varimax rotated
factors
fctr.1 fctr.2
Jpn 0.96874 -0.04003
Eng 0.90061 -0.21177
Hist 0.94468 -0.11008
Math 0.30993 -0.86587
Phys -0.02632 -0.93005
Chem 0.11907 -0.95445
Sorted Pattern for Varimax rotation
fctr.1 fctr.2
Jpn 0.96874 -0.04003
Hist 0.94468 -0.11008
Eng 0.90061 -0.21177
Chem 0.11907 -0.95445
Phys -0.02632 -0.93005
Math 0.30993 -0.86587
Promax rotation was applied...
Pattern for the Promax rotated factors
fctr.1 fctr.2
Jpn 0.98905 0.08941
Eng 0.89541 -0.09616
Hist 0.95468 0.01422
Math 0.19873 -0.84757
Phys -0.15534 -0.95880
Chem -0.00944 -0.96427
Correlation matrix of the factors
fctr.1 fctr.2
fctr.1 1.00000 -0.26097
fctr.2 -0.26097 1.00000
Sorted Pattern for Promax rotation
fctr.1 fctr.2
Jpn 0.98905 0.08941
Hist 0.95468 0.01422
Eng 0.89541 -0.09616
Chem -0.00944 -0.96427
Phys -0.15534 -0.95880
Math 0.19873 -0.84757
Structure Matrix
fctr.1 fctr.2
Jpn 0.96572 -0.16870
Eng 0.92051 -0.32983
Hist 0.95097 -0.23492
Math 0.41992 -0.89943
Phys 0.09488 -0.91826
Chem 0.24220 -0.96180
読み込まれたデータ値が表示されているので、正しくデータ設定されたか確認することができる。
足立浩平・伊藤真道・宇野浩平（2019）行列分解に基づく因子分析とその新展開．計算機統計学、32、61-77．
De Leeuw, J. (2004). Least squares optimal scaling of partially observed linear systems. In Recent Developments on Structural Equations Models: Theory and Applications (van Montfort, K., Oud, J & Satorra, A. (Eds.)), Kluwer Academic Publishers, 121-134.
Gifi, A. (1990). Nonlinear multivariate analysis. John Wiley & Sons.
岡本安晴（2009）統計学を学ぶための数学入門［下］．培風館．
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from rotate import Varimax, Promax
from PlotComp import PlotPattern
from ArrangePat import ArrangePtn
"""
Direct Least Squares Method
for Factor Analysis
Yasuharu Okamoto, 2023.11
"""
#
# Excelデータファイルの読み込み
#
s = input("Input Excel data file
(*.xlsx) = ")
xlsx_fl = pd.ExcelFile(s)
data_xlsx =
pd.read_excel(xlsx_fl, header =
None)
data = data_xlsx.values
f_out_nm = input('Output Text
file(*.txt) = ')
f_out = open(f_out_nm, 'w')
print('data =\n', data)
f_out.write(f'\nInput data =
{f_out_nm}\n')
f_out.write(f'\ndata =\n{data}\n')
var_names = data[0]
data_id = data.T[0][1:]
X = data[1:][:,1:]
X = np.array(X, dtype='float')
n, n_vars = np.shape(X)
print('var_names =\n', var_names)
print('data_id =\n', data_id)
print('n =', n)
print('n_vars =', n_vars)
print('X =\n', X)
var_means = np.mean(X, axis=0)
print('var_means =\n', var_means)
var_sds = np.std(X, axis=0)
print('var_sds =\n', var_sds)
X = (X - var_means) / var_sds
print('X-means =\n', np.sum(X,axis=0)
/ n)
print('X-sds =\n', np.sum(X**2,
axis=0) / n)
U, S, Vh = np.linalg.svd(X,
full_matrices = False) # 特異値分解
print('U =\n', U)
print('S =\n', S)
print('Vh =\n', Vh)
ck_idx = S > 1.0e-7
print('ck_idx =\n', ck_idx)
U = U[:, ck_idx]
print('U =\n', U)
S = S[ck_idx]
print('S =\n', S)
Vh = Vh[ck_idx]
print('Vh =\n', Vh)
print('Singular values =\n', S)
plt.plot(np.arange(len(S)), S)
plt.xticks(np.arange(len(S)),
np.arange(1, len(S) + 1))
plt.ylabel('Singular value', fontsize
= 12)
plt.show()
p = int(input('Number of factors = '))
if p > len(X[0]):
p = len(X[0])
A = Vh.T[:,:p] @ np.diag(S[:p]) /
(n**0.5)
print('A =\n', A)
PlotPtn = PlotPattern(A,
var_names[1:], n_vars, p, "Initial Pattern Matrix" )
while True:
print("\nSet the factor
number (from 1 to %d) to plot the pattern."%p)
print("If you do not
want to plot the pattern, set the number less than 1.")
xdim = int(input("xdim
= "))
if xdim < 1:
break
ydim = int(input("ydim
= "))
PlotPtn.Plot(xdim, ydim)
commu = np.sum(A**2, axis=1)
f_out.write("\n\nPattern matrix
and communalities for the initial factors\n")
f_out.write("\n%10s"%"
")
for i in range(0, p):
f_out.write("%10s"%("fctr.%d"%(i + 1)))
f_out.write("%15s\n"%"communality")
for i in range(0, n_vars):
f_out.write("%10s"%var_names[i+1])
for j in range(0, p):
f_out.write("%10.5f"%A[i, j])
f_out.write("%15.5f\n"%commu[i])
rot = Varimax(A, n_vars, p) # Varimax回転
PatV = rot.rotateV()
PlotPtn = PlotPattern( PatV,
var_names[1:], n_vars, p, "Pattern Matrix(Varimax)" )
while True:
print("\nSet the factor
number (from 1 to %d) to plot the pattern for the Varimax criterion."%p)
print("If you do not
want to plot the pattern, set the number less than 1.")
xdim = int(input("xdim
= "))
if xdim < 1:
break
ydim = int(input("ydim
= "))
PlotPtn.Plot(xdim, ydim)
f_out.write("\n\nVarimax rotation
was applied...\n\n")
f_out.write("\nPattern for the
Varimax rotated factors\n\n")
f_out.write("%10s"%"
")
for j in range(0, p):
f_out.write("%10s"%("fctr.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_vars):
f_out.write("%10s"%var_names[i+1])
for j in range(0, p):
f_out.write("%10.5f"%PatV[i, j])
f_out.write("\n")
ArgPtn = ArrangePtn(var_names[1:],
PatV, n_vars, p)
PtnArgd = ArgPtn.arrange()
f_out.write("\n\nSorted Pattern
for Varimax rotation\n\n")
f_out.write("%10s"%"
")
for j in range(0, p):
f_out.write("%10s"%("fctr.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_vars):
f_out.write("%10s"%PtnArgd[i][0])
for j in range(0, p):
f_out.write("%10.5f"%PtnArgd[i][1][j])
f_out.write("\n")
obrot = Promax(PatV, n_vars, p,
4.0) # Promax回転
Patob, Phi = obrot.rotateP();
f_out.write("\n\nPromax rotation
was applied...\n\n")
f_out.write("\nPattern for the
Promax rotated factors\n\n")
f_out.write("%10s"%"
")
for j in range(0, p):
f_out.write("%10s"%("fctr.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_vars):
f_out.write("%10s"%var_names[i+1])
for j in range(0, p):
f_out.write("%10.5f"%Patob[i, j])
f_out.write("\n")
f_out.write("\n\nCorrelation
matrix of the factors\n\n")
f_out.write("%10s"%"
")
for j in range(0, p):
f_out.write("%10s"%("fctr.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, p):
f_out.write("%10s"%("fctr.%d"%(i + 1)))
for j in range(0, p):
f_out.write("%10.5f"%Phi[i, j])
f_out.write("\n")
ArgPtn = ArrangePtn(var_names[1:], Patob,
n_vars, p)
PtnArgd = ArgPtn.arrange()
f_out.write("\n\nSorted Pattern
for Promax rotation\n\n")
f_out.write('{0:10s}'.format(' '))
for i in range(p):
str = 'fctr.{}'.format(i +
1)
f_out.write('{0:>10s}'.format(str))
f_out.write('\n')
for i in range(0, n_vars):
f_out.write("%10s"%PtnArgd[i][0])
for j in range(0, p):
f_out.write("%10.5f"%PtnArgd[i][1][j])
f_out.write("\n")
S = Patob * Phi
f_out.write('\n\nStructure
Matrix\n\n')
f_out.write('{0:10s}'.format(' '))
for i in range(p):
str = 'fctr.{}'.format(i +
1)
f_out.write('{0:>10s}'.format(str))
f_out.write('\n')
for i in range(n_vars):
f_out.write('{0:>10s}'.format(var_names[i+1]))
for j in range(p):
f_out.write('{0:10.5f}'.format(S[i, j]))
f_out.write('\n')
f_out.close()
print(f'\n{f_out_nm} was saved.\n')