【Pythonで統計モデル】因子分析1：とりあえず1因子モデルを構築してみる
はじめに
株式会社GA technologiesの福中です。今、僕は、Advanced Innovation Strategy Center（以下AISC）という部署でChief Data Scientistとして働いています（詳しい自己紹介はAISC WEBサイトのプロフィールをご覧ください）。今回からしばらくの間、【Pythonで統計モデル】というシリーズの記事を書いていこうと思います。このシリーズではさまざまな統計モデルをPythonで実行する方法について紹介します。
統計分析ソフトの変遷
2000年代初頭の話です。このころは統計分析といえば、多くの人がSPSSやSASを使っていました。もちろんMathematicaやMatlabのような行列言語もあり、それらを利用している人もいましたが、どちらかといえば特殊な数値計算が必要なモデルを研究する人が主に利用していたと記憶しています。大部分の統計ユーザーは利便性の高いSPSSやSASを利用していました。僕が所属していた早稲田大学の豊田研究室でも主にSASを使っていました。これは、豊田研究室の専門の1つが構造方程式モデリング（Structural Equation Modeling: SEM）であり、この分析を実行するためにはSASのCALIS プロシジャがとても優秀だったためです。
時が進み、2005年から2010年にかけて、統計解析ソフトRが一気に浸透しました。統計学会はほぼR一色になったと言っても過言ではありません。Rをプログラミング言語としてみたときは、その設計思想に賛否両論ありましたが、それでも統計分析を行う際の汎用性の高さと拡張性の高さが群を抜いており、何よりプログラミング言語に不慣れな初学者でも扱いやすいという大きなメリットがあったため、急速に普及したのでしょう。
このころも当然Pythonはありましたが、こと統計学者の間で使われることはあまりなかったと記憶しています。僕もPythonの勉強はしましたが、当時はRで解析することがほとんどでした。なぜならPythonは統計ライブラリが貧弱だったことと、何より出力される結果が統計学者には不十分だったからです。その後、機械学習（AI）全盛の時代になり、Pythonが一気に普及しましたが、統計業界でRを使い続けている人が多かったのはこのような理由のためでしょう。
なぜこの記事を書こうと思ったのか？
さて、本題に入ります。上記のように統計モデルを扱う上では不十分な点が多かったPythonですが、現在ではどうなっているのだろうと疑問に思い、調査してみました。「分析をするだけならRでも構わないのでは？」と思われるかもしれません。なぜここまでPythonにこだわるかというと、AISCでは分析をするだけではなく、そのモデルをAWSのLambda等を使ってAPIに実装し、それをエンジンにしたプロトタイプの開発まで1人で行う必要があるからです。アプリを開発する上ではどうしてもRよりPythonに軍配が上がります。また、AISC内の共通言語としてPythonを採用しており、そこに部分的にRを利用すると引継ぎの際に苦労するという理由もあります。したがって、統計モデルを構築する際にもPythonで行うメリットは十分にあると考えました。そこで本記事では、よく利用される統計モデルをPythonで実行することができるのかを調査した結果を書いていきたいと思います。
注意点
本記事は統計学およびPythonの初心者向け解説記事ではありません。普段からR等を使って統計モデルによる分析や研究を行っている研究者の方が、改めてPythonで因子分析を行う必要が出た際のリファレンスとしての利用を想定しています。したがって、因子分析に不慣れな方は、先に以下で紹介するテキストをはじめとした初心者用の教科書で勉強するとよいと思います。
テキスト
本記事では東京図書から出版されている因子分析入門―Ｒで学ぶ最新データ解析―の分析結果を、Pythonを使って再現していきます。分析用データは東京図書のホームページからダウンロードできますので、興味のある方はぜひ試してみてください。
環境
本記事の分析は、OS: Windows 11、Python 3.11.2の環境下で行っています。また、必要なパッケージは以下となります。インストール済みの方は割愛してください。
# packageのインストール
pip install numpy pandas
pip install scikit-learn
pip install factor_analyzer1因子モデル
それでは手始めに、「第1章 速習・因子分析」の最初に出てくる「理科の試験データ」の例を使って、1因子モデルをPythonで表現する方法について学習しましょう。このデータは物理・化学・生物・地学の4教科で構成されており、100点満点で採点された10人の多変量データとなっています。まずは以下のようにしてパッケージとデータの読み込みを行いましょう。
# packageの読み込み
import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from factor_analyzer import FactorAnalyzer# データの読み込み
# データ自体は東京図書のHPからダウンロード可能です。
df = pd.read_csv("dat/理科試験.csv", encoding='shift_jis')
df = df.iloc[:, 1:5]
df.head()ここで使用するStandardScalerは多変量データ（Pandasのデータフレーム）を（平均0、分散1に）標準化するためのもの、FactorAnalyzerが探索的因子分析（exploratory factor analysis: EFA）を行うためのパッケージになります。FactorAnalyzerの詳しい仕様に関しては以下のドキュメントをご覧ください。本記事ではいくつか抜粋してご紹介します。
StandardScalerを使って実際に標準化すると以下のようになります。
# データを平均0、分散1に標準化する
scaler = StandardScaler()
scaler.fit(df)
df_std = pd.DataFrame(scaler.transform(df), columns=df.columns)
df_std.head()
テキストP.4の表1.2に一致していることがわかる
次に1因子モデルで因子分析を行ってみましょう。
# 探索的因子分析を実行する
fa = FactorAnalyzer(n_factors=1, rotation=None, method='ml')
fa.fit(df_std)"n_factors=" で因子数の指定を行います。今回は1因子モデルなので "n_factors=1" にしています。同じく1因子モデルなので、回転は行わず、"rotation=None" （初期解）を指定しています。また、"method=" では推定法を指定でき、今回は最尤推定法（maximum likelihood estimation）を指定しています。なお、推定法に関して、デフォルトでは"method='minres'" となっています。ただし、ミンレス法は結局のところ最小二乗法を非効率的に計算しているだけなので、母数推定には現代統計学のデファクトスタンダードである最尤推定法を選択するようにしましょう。なので、特に理由がない限り"ml" を指定するように注意してください。
参考までにFactorAnalyzerで使用できる回転法と推定法の一覧を以下に掲載しておきます。
None(初期解)
varimax (orthogonal rotation)
promax (oblique rotation、デフォルト)
oblimin (oblique rotation)
oblimax (orthogonal rotation)
quartimin (oblique rotation)
quartimax (orthogonal rotation)
equamax (orthogonal rotation)
minres (ミンレス法、デフォルト)
ml (最尤推定法)
principal (主成分法)
結果の出力をします。ここでは因子負荷量と共通性、独自性を出力しています。
loadings_df = pd.DataFrame(fa.loadings_, columns=["因子負荷量"])
loadings_df.index = df.columns
loadings_df["共通性"] = fa.get_communalities()
loadings_df["独自性"] = fa.get_uniquenesses()
print(loadings_df) 因子負荷量 共通性 独自性
物理 -0.937183 0.878311 0.121689
化学 -0.894242 0.799669 0.200331
生物 -0.864184 0.746814 0.253186
地学 -0.502930 0.252939 0.747061テキストの結果と照らし合わせるとわかるように、小数点第2位の値まで一致しています。ただし、因子負荷量の符号がテキストと逆になっています（テキストでは全部正の値）。しかし、因子分析において因子負荷量の符号が反転することはよくあり、結果に問題ありません。本結果のように因子負荷量が全部負の値だと、因子の名称は「理科学力」というよりは「理科の苦手さ」と命名する方が自然ではあるでしょう。しかし、解釈しにくくなるので、符号を反転して解釈し、「理科学力」と命名しても問題ありません。
最後に因子得点を計算してみましょう。FactorAnalyzerで因子得点を計算する場合、transform() というメソッドを利用します。引数にはデータセットを指定するだけです。
# 因子得点の算出
fa.transform(df_std)array([[-1.70034289],
[-1.0741958 ],
[-0.86126037],
[-0.32192505],
[-0.20332637],
[ 0.26423881],
[ 0.49441201],
[ 0.898387 ],
[ 0.96345161],
[ 1.54056105]])算出された値がテキストと異なりますが、これは因子得点の算出ロジックがテキストとは異なるからです。テキストの因子スコアはRのデフォルトであるサーストンの因子スコアで計算されているのに対して、Pythonのtransformメソッドではバートレットの方法で算出されています。因子スコアを計算する方法には、サーストンやバートレットの方法以外にも数多くの方法が提案されています。しかしながらtransformメソッドは計算方法を指定する引数が存在しないため、バートレットの方法以外で算出したい場合は自分で式を書かないといけないことにご注意ください。
本例の場合、計算結果の値は異なりますが、その順序性は一致しています。また因子得点自体、絶対評価するものではありませんし、どの方法論を使って算出すべきかは一概に決められないので特に問題ありません。また、符号が反転していることにも注意してください。これは因子負荷量がすべて負の値として推定されていたためこのような結果となっています。つまり、因子名を「理科の苦手さ」とした場合の順序になっているということです。これも「理科学力」として解釈したい場合はすべての因子得点の符号を反転させて解釈するとよいでしょう。
この記事が気に入ったらサポートをしてみませんか？