Home

因子分析

 

因子分析の簡単な方法として主因子法(岡本、2014pp. 85-87)、残差に制約条件を設定しない方法として直接最小2乗法(岡本、2014pp. 95-97)について考える。まず、主因子法を取り上げ、次に直接最小2乗法を扱う

 

 

主因子法

 

因子分析のモデルを以下のように設定する(岡本、2014pp. 85-87)。

ここで、は標準化されたデータ値である。は多変量の潜在変数(因子)であり、の各成分は互いに独立に平均0,標準偏差1の分布に従うものとする。残差は、平均が0で、各成分はお互いに独立であり、とも独立な分布に従うものとする。

このとき、の相関行列は次式で与えられる。

式(1)より、からを次の繰り返し法によって求める。

 

 

上の繰り返し法(主因子法)に基づくスクリプトをリストP1に示す。

リストP1のスクリプトを用いて主因子法による因子分析を行うプログラムをリストP2に示す。

スクリプトファイルとサンプルデータファイルは、まとめて圧縮ファイルFAPrg.zipとした。ダウンロード展開すれば、自由に利用できる。

リストP2の実行においては、Excel形式のファイルをPandasライブラリを用いて読み込むときにモジュールopenpyxlが必要である。

モジュールopenpyxlは、次のコマンドでインストールできる。

 

pip install openpyxl

 

リスト2のファイルを次のコマンド

 

PS ******\FAPrg> python FA.py

 

で実行すると、入力データファイル名を聞いてくる。

 

Input Excel data file (*.xlsx) =

 

入力データファイルは、Excelファイルとして用意する(図1)。

図1

 

第1行目に変数名が入力されている。データ値は2行目から設定する。1列目は、データの識別に用いるもので適当に設定する。2列目からが因子分析の対象となる。

データは、値の設定などが終われば保存してからプログラムを実行する。保存先は、リスト2のスクリプトファイルと同じフォルダにしておくと読み込みファイル名の指定が簡単である。

図1のファイル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

 

出力ファイル名を設定すると計算が始まり、固有値が端末画面に出力されるとともに、スクリープロットが表示される(図2)。

図2

 

因子数は2個が適当と思われる。

図2のフォームを閉じると、因子数を聞いてくる。

 

 

                        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 =

 

2を設定してみる。

 

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軸の設定を行うと、図3のようにパターンが表示される。

図3

 

図3のフォームを閉じると、次の描画で用いる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

 

1より小さい数を設定すると、初期解でのパターンの表示が終わり、次に進む。上の例では、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

 

上の設定で、図4のパターンが表示される。

図4

 

図4のフォームを閉じると、次の描画におけるX軸の設定が求められる。1より小さい値を設定すると、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>

 

上の例では、0を設定している。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

 

 

以下のリストP1P2のスクリプトファイルとサンプルデータファイルは、他のスクリプトファイルとともにまとめて圧縮ファイル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.")

 

 

 

直接最小2乗法

 

因子分析のモデルを次式(D1)のように設定する(de Leeuw,2004; Okamoto, 2005)。

ここで、は標準得点化されたデータ行列であり、は共通因子行列、はパターン行列、は残差行列である。足立ら(2019では、はさらに独自因子と残差に分けられているが、ここではモデル(D1)に基づいて分析を行う。モデル(D1)による分析を行うときに、残差の独立性など特別な条件は設定しない。また、ここでの分析法の場合、ヘイウッドケース(Heywood cases:共通性が1より大きくなること)が生じないことが証明されている(Okamoto, 2005, Appendix A)。

モデル(D1)のデータに対する誤差は次式(D2)で与えられる。

誤差(D2)を最小にする解は、Eckart-Young theoremにより、特異値分解で与えられる(Gifi, 1990, p.152, 518; 岡本、2014p.95)。詳しい説明が、岡本(2009pp. 96-100にある。

特異値分解を

とおく。このとき、誤差(D2)を最小にする解は次式で与えられる。

ここで、はデータ数であり、は、から特異値の大きいものから選んだものである。

また、

が成り立つので、も解である。すなわち、解を回転することができる。

 

上で説明した式(D1)による因子分析法のスクリプトをリストD1のように用意した。スクリプトファイルとサンプルデータファイルはまとめて圧縮ファイルDFAPrg.zipとした。ダウンロード展開すれば自由に利用できる。

リストD1のスクリプトを次のコマンドで実行する。

 

PS ******\DFAPrg> python DFA.py

 

実行すると、入力データファイル名を聞いてくる。

 

Input Excel data file (*.xlsx) =

 

入力データファイルは図D1のようなExcelファイルとして用意する。

D1

 

1行目に変数名を設定する。データは2行目以降に設定する。1列目は、データの識別用の値を適宜設定する。2列目以降に因子分析の対象となるデータを設定する。

データの設定が終われば保存する。保存先はプログラムのスクリプトファイルと同じフォルダにしておくと、プログラムからデータファイルの指定を行うとき簡単である。

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

 

特異値の大きさから、因子数は2個が適切であると思われる。

D2のフォームを閉じると、因子数を聞いてくる。次の例では、2と設定している。

 

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

 

上の例では、2が設定されている。

Y軸を設定すると、設定された因子でパターンが表示される(図D3)。

D3

 

第1因子は、データ全体の変動に対応しており、第2因子は文系の科目と理系の科目の区別に対応していることが分かる。

D3のフォームを閉じると、次のパターンの描画のための因子を聞いてくる。

X軸として1より小さい値を設定すると、初期値に対するパターンの描画は終了して、次にVarimax回転後のパターンの描画に進み、Varimax回転後のパターンの描画におけるX軸とY軸を聞いてくる。

次の例では、X軸に1,Y軸に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 = 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軸の因子、上の例では第1因子と第2因子、を設定すると、設定した因子でパターンが描画される(図D4)。

D4

 

第1因子が文系科目、第2因子が理系科目の因子であることが分かる。

D4のフォームを閉じると、次の描画で用いるX軸とY軸の因子の設定が求められる。

X軸として1より小さい値を設定すると、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)行列分解に基づく因子分析とその新展開.計算機統計学、3261-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.

Okamoto, Y. (2005). On the special type of principal component analysis: It should belong to factor analysis. Japan Women's University Journal (Faculty of Integrated Arts and Social Sciences),16,43-50

岡本安晴(2009統計学を学ぶための数学入門培風館

岡本安晴(2014心理学データ分析と測定勁草書房

 

 

 

リストD1 直接最小2乗法による因子分析(DFA.py

 

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

 

 

 

Home