こんにちは、やまもとです。
心理学系の論文を読んでいると、よくリッカート式のアンケート結果に対する因子分析の結果が掲載されています。因子分析には、探索的因子分析と確証的因子分析の２種類があります。探索的因子分析は、どのような因子が得られるか分からずに行う因子分析です。確証的因子分析は、理論的に因子構造が予測されていて、それを確認するために行う因子分析のことです。
一方、因子分析の計算では、因子数が予め知られていなければなりません。これは、因子分析では、因子数が入力パラメータだからです。確証的因子分析の場合は、理論的に因子数も予測されているため、このことは問題になりません。しかし、探索的因子分析では、因子数も因子構造も分からない状況で分析を始めるため、このことが問題になります。
因子数を予測する方法として、相関行列を見て予測する方法や、スクリープロットを見て予測する方法があります。しかし、これらの方法は人が目で見て判断する方法なので、分析者の能力に依存するところが大きくなります。しかも、予測された因子数で因子分析を行ったとしても、内的一貫性などの検証の結果、再分析が必要になることもあります。そのため、因子分析は、最終的に良い結果が得られるまで、属人的で試行錯誤と手間隙と時間のかかる作業が必要です。
しかし、人が行うべき価値ある作業は、分析結果の解釈です。
そこで、この時間のかかる作業を一部自動化するクラスExploratoryFactorAnalysisを作ってみました。
目次
使い方
テストデータ
テストデータは、下記で公開されているBig5+知性のデータを使います。Excelファイルでダウンロードして、CSV形式で保存しておきます。
実行方法
ExploratoryFactorAnalysisクラスのを使った分析方法は、下記のようなシンプルなものです。
import pandas as pd
d = pd.read_csv("data08-01.csv",index_col=0) #データ読み込み
efa = ExploratoryFactorAnalysis() #インスタンス作成
nf = efa.explore(d) # 因子数推定
fs = efa.analyze(d,nf) # 因子分析
print(efa) # 分析結果表示
- explore関数の役割
- 最適化によって、全ての因子について内的一貫性が閾値より大きくなる因子数を推定します。
- 最適化の反復を繰り返すため、データによっては計算時間がかかる場合があります。
- analyze関数の役割
- 実際に因子分析を行い、分析結果として必要な値を一通り計算します。
- 計算結果を、インスタンスefaのメンバーに格納しておきます。
実行結果
efa.explore(d)を実行すると、次のようなログを出力します。
predicted max_factors = 2
try nfactors = 2
found optimum nfactors = 2
Predicted nfactors = 2
テストデータだと、因子数＝２を１回目で当ててしまっているので、ありがたみがありませんね。こ
のような幸運なデータでなければ、”try nfactors = “が繰り返し出力され、現在どの因子数を試しているのかが分かります。
関数analyze()を実行した後に、print(efa)を実行すると、次のように分析結果が出力されます。
predicted max_factors = 2
try nfactors = 2
found optimum nfactors = 2
Predicted nfactors = 2
'fa' does not have a menber 'psi_'
Likert's scale : (1, 7)
Cutoff Loading : 0.5
Cutoff Alpha : 0.7
Rotation : promax
Num. of Factors: 2
Loadings:
F1 F2
素直さ 0.935606 -0.093399
知性 0.660452 0.041540
信頼性 0.652540 0.030894
積極性 -0.112861 0.867653
外向性 0.074178 0.701946
社交性 0.032645 0.662783
Cronbach's alpha:
F1 F2
alpha 0.785032 0.783893
Split-half rho:
F1 F2
rho 0.722344 0.737691
I-T Correlations:
F1 F2
素直さ 0.881199 NaN
知性 0.822455 NaN
信頼性 0.805127 NaN
積極性 NaN 0.859822
外向性 NaN 0.831353
社交性 NaN 0.815611
Rotation Matrix:
F1 F2
F1 0.627145 -0.583059
F2 0.871851 0.901929
Factor Correlations:
F1 F2
F1 1.000000 0.204861
F2 0.204861 1.000000
Communalities:
comm.
外向性 0.498231
社交性 0.440347
積極性 0.765560
知性 0.437923
信頼性 0.426763
素直さ 0.884081
Uniquenesses:
uniq.
外向性 0.501769
社交性 0.559653
積極性 0.234440
知性 0.562077
信頼性 0.573237
素直さ 0.115919
因子負荷量（Loadings）は、論文で見るような形式（因子が共通する質問毎に降順にソートしてある形式）で出力されていることが分かります。因子負荷量は、Cutoff Loadingの値を下限値として、それ以上の値を持つ質問を並べ替えています。
内的一貫性（Cronbach’s alpha）は、因子の信頼性を図るための１つの指標です。0.8以上であれば全く問題ありませんが、0.7でも大丈夫と判断される場合もあります。因子数推定では、Cutoff Alphaの値を下限値として推定していました。
参考：「因子分析」の解釈の仕方
因子分析の解説は、下記が参考になりました。
- http://www.u.tsukuba.ac.jp/~hirai.akiyo.ft/forstudents/eigokyouikugaku7.files/M4_9(N.K.).pdf
- http://www.u.tsukuba.ac.jp/~hirai.akiyo.ft/forstudents/HP%20title/ibunka2014/20141008%20Ch.9.pdf
コード解説
コード全体は最後に記載することにして、少しコードの解説をしておきたいと思います。自分も忘れてしまいそうなので！
使用パッケージ
ほとんどは標準的なパッケージですが、因子分析にはfactor-analyzerを使用しています。
- NumPy
- PANDAS
- MatPlotLib
- japanize-matplotlib
- factor-analyzer
因子分析は、scikit-learnにも実装されていますが、論文でよく使われるpromax回転が未実装だったので利用を諦めました。
Pythonでの因子分析のやり方は、下記を参考にしました。
scikit-learnを使った因子分析のやり方
factor-analyzerを使った因子分析のやり方
自作関数
いくつか自分で作ったメンバー関数があるので、それらを解説しておきます。
クロンバックαの計算
計算は、統計ウェブのクロンバックα係数に基づいて、下記の公式を使用しています。
これを、コードにすると次のようになります。pandas.DataFrame型で、因子で使用する質問のみの回答データdataが入力されることを想定しています。
def cronbach_alpha(data, ddof=False):
#ddof=True:不偏分散、False:標本分散
M = data.shape[1]
if M > 1 :
x = data.sum(axis=1) # 行方向への和
sigma_x = x.var(ddof=ddof)
sigma_i = 0e0
for i in range(M):
sigma_i += data.iloc[:,i].var(ddof=ddof)
cronbach_alpha = M / (M-1) * ( 1 - sigma_i/sigma_x )
else:
cronbach_alpha = float(M)
return cronbach_alpha
このコードは、公式通りにプログラミングしただけですが、例外としてm=0の場合とm=1の場合を条件分岐させておく必要があります。理由と対処は、以下の通りです。
- m=0の場合
- 公式からα=0は明らかで計算が無駄になる。
- そこで、計算せずにα=m=0を返却する。
- m=1の場合
- 0÷0が発生し公式では計算できない。
- そこで、0÷0=1として公式を変形し、計算せずにα=m=1を返却する。
標準化の計算
特にリッカート式のアンケートの場合、ローデータだと各値=1,2,3,4,…と大きすぎるため、ローデータのまま因子分析するとおかしな結果になります。そのため、データを標準化あるいは規格化しないといけません。
そこで、下記のようなコードを実装しました。
def standarize(data,ddof=False):
stddat = data - data.mean() / data.std(ddof=ddof)
stddat.index = data.index
stddat.columns = data.columns
return stddat
参考
因子負荷量からクロンバックαを計算する内部関数
クロンバックαを計算するには、ある１つの因子について
- 因子負荷量の大きい質問群を選ぶ（因子が影響している質問を選ぶ）
- 選択した質問群のみの回答データ（標準化済）を作る
- 選択された回答データからクロンバックαを計算する
という処理が因子の数だけ必要です。
ただし、因子負荷量は負に大きい場合もあります。その場合、回答データを逆転して意味を逆にする必要があります。
この処理は、Excelでも実現することができますが、かなり手間がかかります。そのため、自動化しておくと便利です。
この処理を自動化するコードは、次のようなになりました。
def __loading_to_alpha(self,loading,data):
# リスト型だと、pandas.DataFrameに変換するとき、行と列が逆になってしまうため、辞書型にした。
item_dict = {}
# 因子負荷量の絶対値が閾値より大きい質問を取り出す
for i in range(len(loading)):
if loading[i] >= self.tol_loading_ :
item_dict[data.columns[i]] = data.values[:,i]
elif loading[i] <= -self.tol_loading_ : #逆転処理
item_dict[data.columns[i]] = self.__reverse(data.values[:,i])
#クロンバックαの計算
alpha = 0e0
if len(item_dict) > 1 : # １質問の場合は、共通因子はなかったと見なす
alpha = self.cronbach_alpha(pd.DataFrame(item_dict),ddof=self.ddof_)
#print(item_dict.keys()," alpha="+str(alpha))
return alpha
実装上の注意点は２つあります。
- 因子を構成する質問を一時的に辞書item_dictに格納
- 辞書型なのは、pandas.DataFrame変換時に列として認識させるためです。リスト型だと、行として認識されます。
- 因子を構成する質問が１問の場合は、α=0にして因子化に失敗させる
- １質問で１因子が構成された場合、それは共通因子とは呼べないため、除外しています。
因子負荷量のソート
因子分析の作業の中で、意外と手間がかかるのは、因子負荷量の並べ直し作業です。
因子負荷量は、論文に載せるとき、因子毎に降順に並べ直します。ただし、この時も、
- 正負の因子負荷量を混ぜてソートする必要がある
- 全ての因子について１つ１つソートする必要がある
といった要求があり、これを手作業で行うと時間がかかります。
そのため、これも自動化しておくと後々楽になります。
この処理をコードに直すと次のようになります。
def sort_loadings(loadings,tol=0.4):
i=1
n=len(loadings.columns)
names=[]
while i <= n:
name="A"+str(i)
names.append(name) # 絶対値のカラム名を作成
i+=1
abs_loadings = loadings.abs() # 絶対値テーブル
abs_loadings.columns = names # カラム名をせっと
tmp_loadings = loadings.copy() # 元データを破壊しないために、deep copyで一時変数を用意
tmp_loadings = tmp_loadings.join(abs_loadings) # 一時変数と絶対値のテーブルを結合
results=pd.DataFrame() # 空のDataFrameを作成
ilabel=n # データが、F1,...,Fn,A1,...,Anの並びなので、A1から開始するためにnにする
for label in names: # A1,...,Anのループ
x = tmp_loadings.sort_values(label,ascending=False) # Axカラムで降順にソート
j=0
while j < len(x.index) and x.iloc[j,ilabel] >= tol: # 因子負荷量がtol以上の要素数をカウント
j += 1
results = results.append(x.iloc[0:j,0:n]) # 因子負荷量>=tolの部分だけをコピー
tmp_loadings = x.iloc[j:,:] # 因子負荷量>=tolの要素を削除して次へ
ilabel += 1
return results
因子負荷行列に、その絶対値の行列を結合して、絶対値のデータで因子ごとに降順にソートしていく処理をしています。
最適化アルゴリズム
因子数を推定するアルゴリズムはとても単純で、スクリープロットで最大因子数を取得したら、それを１つずつ減らしながら因子分析を実行し、全因子で内的一貫性があれば止める、というものです。アルゴリズムにすると、次のようになります。
- 相関行列の固有値を計算し、１以上の固有値の個数をmax_factorsとする
- 因子数を初期化：nfactors = max_factors
- 因子分析を実行し、因子負荷行列を取得：loadings
- 因子番号を初期化：ifactor = nfactors – 1
- 第ifactor因子のクロンバックα係数を計算：alpha
- α係数が閾値tol未満であれば、8.へ（内的一貫性がないため）
- 次の因子番号へ変更し、5.へ戻る：ifactor = ifactor – 1
- ifactor<0の場合、10.へ
- ifactor>=0の場合、nfactorsを変更して、3.へ戻る：nfactors=nfactors-1
- nfactorsを戻す
内側のループで、因子番号ifactorを最大値からデクリメントしているのがコツです。これは、第１因子よりも第N因子の方が内的一貫性が崩れることが多いためです。そのため、第１因子からチェックしていくと、失敗が判明するまでに時間がかかります。早く失敗がわかれば、すぐ次へ行けるため、第N因子からチェックするようにしています。
ソースコード
下記に、ソースコードを載せておきます。使用する場合は、コピペして実行してください。（2023年2月12日更新）
GitHub : https://github.com/jyam45/stat_analytics
ExploratoryFactorAnalysisクラス
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
from factor_analyzer import FactorAnalyzer
class ExploratoryFactorAnalysis:
def __init__(self,
likert=(1,7),tol_loading=0.5,tol_alpha=0.7,rotation='promax',
ddof=False,method='minres',bounds=(0.005,1),impute='median'):
self.likert_ = likert # Likert's scale for reversing (default:(1,7))
self.tol_loading_ = tol_loading # Cutoff value for factor loadings (default:0.5)
self.tol_alpha_ = tol_alpha # Cutoff value for Cronbach's alpha (default:0.7)
self.rotation_ = rotation # [None|'varimax'|'promax'|'oblimin'|'oblimax'|'quartimin'|'quatimax'|'eqamax']
self.ddof_ = ddof # (default:False) 不偏分散
self.method_ = method # ['minres'|'ml'|'principal']
self.bounds_ = bounds # L-BFGS-M bounds
self.impute_ = impute # ['median'|'drop'|'mean']
self.nfactors_ = 0 # number of factors
self.max_factors_ = 0 # maximum number of factors
self.scree_eigvals_ = None # Eigen values for scree plot
self.factors_ = None # factor values
self.alphas_ = None # Cronbach's alpha
self.loadings_ = None # The factor loading matrix
self.corr_ = None # The original correlation matrix
self.rotation_matrix_ = None # The rotation matrix
self.structure_ = None # The structure loading matrix for promax
self.factor_corr_ = None # The factor correlations matrix (psi)
self.communalities_ = None #
self.eigenvalues_ = None
self.factor_variance_ = None
self.uniquenesses_ = None
self.sorted_loadings_ = None
self.itcorr_ = None
self.rhos_ = None
#クーロンバックα係数
@staticmethod
def cronbach_alpha(data, ddof=False):
#ddof=True:不偏分散、False:標本分散
M = data.shape[1]
if M > 1 :
x = data.sum(axis=1) # 行方向への和
#分散
sigma_x = x.var(ddof=ddof)
sigma_i = 0e0
for i in range(M):
sigma_i += data.iloc[:,i].var(ddof=ddof)
cronbach_alpha = M / (M-1) * ( 1 - sigma_i/sigma_x )
else:
cronbach_alpha = float(M)
return cronbach_alpha
# I-T相関
# https://bellcurve.jp/statistics/glossary/7396.html
@staticmethod
def item_total_corr(data):
m = len(data.columns)
x = data.copy() # deep copy
x = x.join(pd.DataFrame(data.sum(axis=1),columns=["I-T corr."])) # 行方向の和を追加する
c = x.corr()
z = pd.DataFrame(c.loc[:,"I-T corr."].iloc[:m,])
z = z.sort_values("I-T corr.",ascending=False)
return z
# 折半法(split-half method)
# https://bellcurve.jp/statistics/glossary/7447.html
@staticmethod
def split_half_rho(data):
x = {}
x["even"] = data.iloc[:,0::2].sum(axis=1)
x["odd" ] = data.iloc[:,1::2].sum(axis=1)
z = pd.DataFrame(data=x)
c = z.corr()
r = c.loc["even","odd"]
rho = 2*r /(1+r)
return rho
#標準化
#https://note.nkmk.me/python-list-ndarray-dataframe-normalize-standardize/
@staticmethod
def standarize(data,ddof=False):
stddat = data - data.mean() / data.std(ddof=ddof)
stddat.index = data.index
stddat.columns = data.columns
return stddat
#逆転（内部関数）
def __reverse(self,value):
#if self.rtype is "likert" :
return (self.likert_[0]+self.likert_[1]) - value
#因子負荷量から項目データを選択する（内部関数）
def __select_items_by_loading(self,loading,data):
# 格納用辞書：pandas.DataFrameに変換するときに、列データとするため辞書に
item_dict = {}
# 因子負荷量の絶対値が閾値より大きい質問を取り出す
for i in range(len(loading)):
if loading[i] >= self.tol_loading_ :
item_dict[data.columns[i]] = data.values[:,i]
elif loading[i] <= -self.tol_loading_ :
#item_dict[data.columns[i]] = self.__reverse(data.values[:,i])
item_dict[data.columns[i]] = -data.values[:,I] # dataは標準化済みと仮定する
return pd.DataFrame(item_dict)
#因子負荷量からクロンバックαを計算する（内部関数）
def __loading_to_alpha(self,loading,data):
items = self.__select_items_by_loading(loading,data)
#クロンバックαの計算
alpha = 0e0
if len(items.columns) > 1 : # １質問の場合は、共通因子はなかったと見なす
alpha = self.cronbach_alpha(items,ddof=self.ddof_)
#print(item_dict.keys()," alpha="+str(alpha))
return alpha
#I-T相関を計算する
def __loading_to_itcorr(self,loading,data):
items = self.__select_items_by_loading(loading,data)
itcorr= self.item_total_corr(items)
return itcorr
#折半法を計算する
def __loading_to_rho(self,loading,data):
items = self.__select_items_by_loading(loading,data)
rho = self.split_half_rho(items)
return rho
#因子間相関を計算する
def __loading_to_corr(self,loadings,data):
m = len(data.index)
nf = len(loadings.columns)
fnames = loadings.columns
x = pd.DataFrame(np.zeros(m*nf).reshape(m,nf),index=data.index,columns=fnames)
for i in range(nf):
picks = self.__select_items_by_loading(loadings.iloc[:,i],data)
x[fnames[i]] = picks.mean(axis=1)
return x.corr()
#複数因子に含まれる項目を見つける
def __find_drop_items(self,loadings):
#print(loadings)
a = (abs(loadings) >= self.tol_loading_).sum(axis=1) # 複数因子に含まれる項目確認
df = pd.DataFrame(a,index=loadings.index,columns=["count"])
#print(df)
x = df[ df["count"] > 1 ]
return x.index
#ソート済因子負荷量
@staticmethod
def sort_loadings(loadings,tol=0.4):
i=1
n=len(loadings.columns)
#print(n)
names=[]
while i <= n:
name="A"+str(i)
names.append(name)
i+=1
abs_loadings = loadings.abs()
abs_loadings.columns = names
tmp_loadings = loadings.copy() # deep copy
tmp_loadings = tmp_loadings.join(abs_loadings)
#print(tld)
results=pd.DataFrame()
ilabel=n
for label in names:
x = tmp_loadings.sort_values(label,ascending=False)
#print(label)
#print(x)
# explore
j=0
while j < len(x.index) and x.iloc[j,ilabel] >= tol:
j += 1
#print(x.iloc[0:j,0:n])
results = results.append(x.iloc[0:j,0:n])
tmp_loadings = x.iloc[j:,:]
ilabel += 1
#print(results)
return results
#因子数推定
def explore( self, data ):
#標準化
stddat = self.standarize(data,ddof=self.ddof_)
#最大因子数推定
ev = pd.DataFrame(np.linalg.eigvals(stddat.corr()))
max_factors = 0
for e in ev[0]:
if e > 1.0 :
max_factors += 1
print('predicted max_factors = '+str(max_factors))
#因子数探索ループ
nfactors = max_factors
while nfactors > 0 :
print('try nfactors = '+str(nfactors))
#因子分析
fa = FactorAnalyzer(nfactors,rotation=self.rotation_,method=self.method_,
bounds=self.bounds_,impute=self.impute_)
fa.fit(stddat) # 因子分析（特異値分解SVD）
#因子負荷量(因子×質問）
#loadings = pd.DataFrame(fa.components_.T)
#loadings = pd.DataFrame(fa.loadings_)
loadings = pd.DataFrame(fa.loadings_,index=stddat.columns)
#信頼性＆妥当性（内的一貫性）のチェック
ifactor = nfactors - 1
while ifactor >= 0:
alpha = self.__loading_to_alpha(loadings.iloc[:,ifactor],stddat)
if alpha < self.tol_alpha_ : # 閾値未満なら一貫性がないため、因子数が妥当ではない。
break
ifactor -= 1
#探索終了条件：全因子の妥当性が確認された場合、探索を打ち切る
if ifactor < 0 :
#print('found optimum nfactors = '+str(nfactors))
#break
drop_items = self.__find_drop_items(loadings) #重複項目の探索
#print(stddat)
if len(drop_items) > 0 :
print("Dropped Items : ",drop_items)
stddat = stddat.drop(drop_items,axis=1)
nfactors = max_factors+1 # 再探索する
else:
print('found optimum nfactors = '+str(nfactors))
break
#探索継続条件：因子が確認できなかった場合、因子数を減らして次へ。
nfactors -= 1
#推定因子数（もし0なら失敗）
if nfactors > 0:
print("Predicted nfactors = "+str(nfactors))
else:
print("Prediction error : No factors in this data, because nfactors = 0")
#インスタンス変数に格納
self.max_factors_ = max_factors
self.nfactors_ = nfactors
self.scree_eigvals_ = ev
return nfactors # 成功 nfactors>0, 失敗 nfactors=0
def scree_plot(self):
# 基準線(固有値1)
ev_1 = np.ones(len(self.scree_eigvals_))
# 変数
plt.plot(self.scree_eigvals_, 's-') # 主成分分析による固有値
plt.plot(ev_1, 's-') # ダミーデータ
# 軸名
plt.xlabel("因子の数")
plt.ylabel("固有値")
plt.grid()
plt.show()
#因子分析の解説
#
def analyze(self,data,nfactors):
#エラー処理
if nfactors < 1 :
raise ValueError('NFACTORS in analyze() must be more than zero!')
#仮因子名生成 "Fn"
factnames = [] # 因子名配列
i=0
while i < nfactors: #因子ラベル生成
factname = 'F' + str(i+1)
factnames.append(factname)
i+=1
#標準化
stddat = self.standarize(data,ddof=self.ddof_)
#因子分析（MINRES法を使用）
fa = FactorAnalyzer(nfactors,rotation=self.rotation_,method=self.method_,
bounds=self.bounds_,impute=self.impute_)
fa.fit(stddat)
#重複項目の削除 & 再分析
drop_items = self.__find_drop_items(pd.DataFrame(fa.loadings_,index=data.columns,columns=factnames))
if len(drop_items) > 0 :
print("Dropped Items : ",drop_items)
stddat = stddat.drop(drop_items,axis=1)
fa.fit(stddat)
#因子負荷量（質問×因子）
self.loadings_ = pd.DataFrame(fa.loadings_,index=data.columns,columns=factnames)
self.sorted_loadings_ = self.sort_loadings(self.loadings_,tol=self.tol_loading_)
#内的一貫性
self.alphas_ = pd.DataFrame(np.zeros(nfactors).reshape(1,nfactors),index=["alpha"],columns=factnames)
for i in range(nfactors):
alpha = self.__loading_to_alpha(self.loadings_.iloc[:,I],stddat)
self.alphas_.iloc[0,i] = alpha
#I-T相関
self.itcorr_ = pd.DataFrame()
for i in range(nfactors):
itcorr = self.__loading_to_itcorr(self.loadings_.iloc[:,I],stddat)
itcorr.columns = [factnames[i]]
self.itcorr_ = self.itcorr_.append(itcorr)
#split-half
self.rhos_ = pd.DataFrame(np.zeros(nfactors).reshape(1,nfactors),index=["rho"],columns=factnames)
for i in range(nfactors):
rho = self.__loading_to_rho(self.loadings_.iloc[:,I],stddat)
self.rhos_.iloc[0,i] = rho
#因子得点
self.factors_ = pd.DataFrame(fa.transform(stddat),index=data.index,columns=factnames) # 因子得点に変換
#共通性（共通分散の値のこと、複数の観測変数に影響する度合い）
x = fa.get_communalities()
if x is not None :
self.communalities_ = pd.DataFrame(x,index=data.columns,columns=['comm.'])
#print(self.communalities_)
#オリジナル相関の固有値、因子相関の固有値
x = fa.get_eigenvalues()
if x is not None :
self.eigenvalues_ = pd.DataFrame(x,index = ['original','factor'],columns = data.columns)
#print(self.eigenvalues_)
#因子寄与、因子寄与率、累積因子寄与率
x = fa.get_factor_variance()
if x is not None :
self.factor_variance_ = pd.DataFrame(x,index=['var.','prop.','cumu.'],columns = factnames)
#print(self.factor_variance_)
#独自性（独自分散の値のこと、単数の観測変数に影響する度合い）
x = fa.get_uniquenesses()
if x is not None :
self.uniquenesses_ = pd.DataFrame(x,index = data.columns,columns = ['uniq.'])
#print(self.uniquenesses_)
#相関行列（項目間の相関）
self.corr_ = pd.DataFrame(fa.corr_,index=data.columns,columns=data.columns)
#print(self.corr_)
#回転行列（斜交回転の場合だけ計算され、軸の回転を表す）
if hasattr(fa,'rotation_matrix_') :
self.rotation_matrix_ = pd.DataFrame(fa.rotation_matrix_,index=factnames,columns=factnames)
#print(self.rotation_matrix_)
else:
print("'fa' does not have a menber 'rotation_matrix_'")
#構造行列（promaxの場合だけ計算される）
if hasattr(fa,'structure_'):
self.structure_ = pd.DataFrame(fa.structure_,index=data.columns,columns=factnames)
#print(self.structure_)
else:
print("'fa' does not have a menber 'structure_'")
#因子相関行列（斜交回転の場合だけ計算される、promaxで計算されないのはバグ？）
if hasattr(fa,'psi_'):
self.factor_corr_ = pd.DataFrame(fa.psi_)
#print(self.factor_corr_)
else:
print("Computed factor's correlations instead because 'fa' does not have a member 'psi_'")
self.factor_corr_ = self.__loading_to_corr(self.loadings_,stddat)
return self.factors_
# CFAのための構造方程式を計算する
# formating : 出力形式 [lavaan|sem] (default:lavaan)
# model_type : CFAモデル型 [order1|order2|bifactor|multilayer](default:order1)
def cfa_model(self,formating="lavaan",model_type="order1",gfactor="GF"):
model = ""
if formating == "lavaan":
left = "~"
right= "=~"
both = "~~"
elif formating == "sem":
left = "<-"
right= "->"
both = "<>"
else:
print(f"cfa_model: irregal option formating={formating}. ['lavaan','sem']")
factors = self.loadings_.columns
items = self.loadings_.index
#潜在変数の定義
for factor in factors:
for item in items:
loading = self.sorted_loadings_.loc[item,factor]
if abs(loading) > self.tol_loading_:
model += f"{factor} {right} {item}\n"
#model += f"{factor} {both} 1*{factor}\n"
#CFAモデル
if model_type == "order1":
# 通常のCFAモデル
nfactors = len(factors)
for j in range(1,nfactors):
jfactor = factors[j]
for i in range(0,j):
ifactor = factors[i]
model += f"{ifactor} {both} {jfactor}\n"
elif model_type == "order2":
# ２次因子モデル
nfactors = len(factors)
for j in range(0,nfactors):
jfactor = factors[j]
model += f"{gfactor} {right} {jfactor}\n"
#model += f"DEFINE(latent) {gfactor}\n"
elif model_type == "bifactor":
# bifactorモデル
for item in items:
model += f"{gfactor} {right} {item}\n"
# 潜在因子の独立性
for factor in factors:
model += f"{gfactor} {both} 0*{factor}\n"
nfactors = len(factors)
for j in range(1,nfactors):
jfactor = factors[j]
for i in range(0,j):
ifactor = factors[i]
model += f"{ifactor} {both} 0*{jfactor}\n"
#model += f"DEFINE(latent) {gfactor}\n"
#model += f"{gfactor} {both} 1*{gfactor}\n"
elif model_type == "multilayer":
# 階層因子モデル
nfactors = len(factors)
for j in range(0,nfactors):
jfactor = factors[j]
model += f"{gfactor} {right} {jfactor}\n"
for item in items:
model += f"{gfactor} {right} {item}\n"
#model += f"DEFINE(latent) {gfactor}\n"
else:
print("cfa_model: invalide option model_type={model_type}. ['order1','order2','bifactor','multilayer']")
return model
def __str__(self):
x = "\n"
x += "Likert's scale : " + str(self.likert_) + "\n"
x += "Cutoff Loading : " + str(self.tol_loading_) + "\n"
x += "Cutoff Alpha : " + str(self.tol_alpha_) + "\n"
x += "Rotation : " + self.rotation_ + "\n"
x += "Num. of Factors: " + str(self.nfactors_) + "\n"
x += "\n"
if self.loadings_ is not None:
pd.set_option('display.max_columns', len(self.loadings_.columns))
pd.set_option('display.max_rows', len(self.loadings_.index))
if self.sorted_loadings_ is not None:
x += "Loadings:\n" + str(self.sorted_loadings_) + "\n\n"
if self.alphas_ is not None:
x += "Cronbach's alpha:\n" + str(self.alphas_) + "\n\n"
if self.rhos_ is not None:
x += "Split-half rho:\n" + str(self.rhos_) + "\n\n"
if self.itcorr_ is not None:
x += "I-T Correlations:\n" + str(self.itcorr_) + "\n\n"
if self.rotation_matrix_ is not None:
x += "Rotation Matrix:\n" + str(self.rotation_matrix_) + "\n\n"
if self.factor_corr_ is not None:
x += "Factor Correlations:\n" + str(self.factor_corr_) + "\n\n"
if self.factor_variance_ is not None:
x += "Contributions:\n" + str(self.factor_variance_) + "\n\n"
if self.communalities_ is not None:
x += "Communalities:\n" + str(self.communalities_) + "\n\n"
if self.uniquenesses_ is not None:
x += "Uniquenesses:\n" + str(self.uniquenesses_) + "\n\n"
return x
実行コードサンプル
import pandas as pd
import factor_analysis as fa
#import sem_analysis as sa
import semopy
d = pd.read_csv("data.csv",index_col=0) #データ読み込み
efa = fa.ExploratoryFactorAnalysis() #インスタンス作成
nf = efa.explore(d) # 因子数推定
fs = efa.analyze(d,nf) # 因子分析
print(efa) # 分析結果表示
print("CFA model:")
cfa = efa.cfa_model(formating="lavaan")
print(cfa)
#cfi = StructEquationModel(cfa)
mod = semopy.Model(cfa)
res = mod.fit(d)
print(res)
ins = mod.inspect()
pd.set_option('display.max_columns',len(ins.columns))
pd.set_option('display.max_rows',len(ins.index))
print(ins)
gfi = semopy.calc_stats(mod)
pd.set_option('display.max_columns',len(gfi.index))
pd.set_option('display.max_rows',len(gfi.columns))
print(gfi.T)
print("Graphviz:")
g = semopy.semplot(mod,"dat.png",plot_covs=True,std_ests=True,show=True)
print(g)
あとは、フロントエンドを用意すれば、データのコピペだけで因子分析できますね。
