【Pythonで統計モデル】因子分析3：尺度構成のための因子分析
はじめに
株式会社GA technologiesの福中です。今、僕は、Advanced Innovation Strategy Center（以下AISC）という部署でChief Data Scientistとして働いています（詳しい自己紹介はAISC WEBサイトのプロフィールをご覧ください）。前回に引き続き、【Pythonで統計モデル】というシリーズの記事を書いていこうと思います。前回は、実践的な多因子モデルを構築し、Pythonで本格的な探索的因子分析ができることを確認しました。
第3回目である今回は、心理学分野でよく研究されている「尺度構成」について、Pythonを使った因子分析で行う方法を見ていきます。
注意点
本記事は統計学およびPythonの初心者向け解説記事ではありません。普段からR等を使って統計モデルによる分析や研究を行っている研究者の方が、改めてPythonで因子分析を行う必要が出た際のリファレンスとしての利用を想定しています。したがって、因子分析に不慣れな方は、先に以下で紹介するテキストをはじめとした初心者用の教科書で勉強するとよいと思います。
テキスト
本記事では東京図書から出版されている因子分析入門―Ｒで学ぶ最新データ解析―の分析結果を、Pythonを使って再現していきます。分析用データは東京図書のホームページからダウンロードできますので、興味のある方はぜひ試してみてください。
環境
本記事の分析は、OS: Windows 11、Python 3.11.2の環境下で行っています。また、必要なパッケージは以下となります。インストール済みの方は割愛してください。
# packageのインストール
pip install numpy pandas matplotlib japanize_matplotlib
pip install pingouin
pip install factor_analyzerなおPythonにおける因子分析のパッケージであるfactor_analyzerのドキュメントは以下となります。本記事ではここから重要な部分をいくつかピックアップしてご紹介します。網羅的な解説はしませんので、適宜、必要に応じてご参照ください。
項目分析
それでは尺度構成を行う際の基本である項目分析から始めていきましょう。今回使用するデータは、テキスト第2章で利用されている心理的ストレス反応を測定するために開発された「心理的ストレス反応データ」です。このデータは全37項目、すべて5件法で構成されています。まずは以下のようにしてパッケージとデータの読み込みを行いましょう。
# packageの読み込み
import numpy as np
import pandas as pd
import pingouin as pg
import matplotlib.pyplot as plt
import japanize_matplotlib
from factor_analyzer import FactorAnalyzer# 心理的ストレスデータの読み込み
df = pd.read_csv("dat/SR2.csv", encoding='shift_jis')
df.head()
次にテキストP.37 の表2.1で示されているような各設問の平均と標準偏差および項目全体相関を算出してみましょう。各項目の平均と標準偏差については、pandasのdescribeメソッドを利用すると簡単に算出できます。
# 平均と標準偏差の算出
statistic = df.describe().loc[('mean','std'),].T
print(statistic) mean std
C1 2.850962 1.297597
C2 3.319712 1.242305
D1 3.564904 1.143075
A1 2.906250 1.279442
C3 2.850962 1.261821
F1 2.769231 1.347775
C4 2.865385 1.238403
C5 2.901442 1.273416
F2 2.855769 1.268094
A2 2.536058 1.323976
F3 2.367788 1.143582
A3 2.971154 1.333221
F4 3.144231 1.254722
D2 2.850962 1.232843
F5 2.242788 1.204888
A4 2.853365 1.283871
D3 3.466346 1.221819
B1 2.199519 1.301010
A5 2.432692 1.222397
F6 3.134615 1.332143
D4 2.605769 1.197735
D5 2.290865 1.169703
E1 3.144231 1.352682
B2 2.206731 1.261307
...
F10 1.848558 1.043002
B5 1.701923 1.063115
C6 3.278846 1.354736
C7 1.819712 1.010163次に項目全体相関ですが、これを算出するためには最初に行ごとの合計点、すなわち検査協力者ごとの心理的ストレス反応得点（テスト得点）を計算しておかないといけません。スクリプトは以下のように記述すればよいでしょう。
# 心理的ストレス反応得点（合計点）の算出
df['sum'] = df.sum(axis=1)
df.head()
相関係数の算出にはpandasのcorrメソッドを利用します。以下では心理ストレス反応得点（テスト得点）と各列との相関係数を計算していますが、テスト得点同士の相関は1となり不要なので削除していることに注意してください。
# 項目全体相関の算出
statistic['R'] = df.corr()['sum'].drop('sum')
print(statistic) mean std R
C1 2.850962 1.297597 0.443808
C2 3.319712 1.242305 0.386667
D1 3.564904 1.143075 0.453415
A1 2.906250 1.279442 0.602114
C3 2.850962 1.261821 0.529553
F1 2.769231 1.347775 0.574909
C4 2.865385 1.238403 0.510091
C5 2.901442 1.273416 0.540531
F2 2.855769 1.268094 0.717014
A2 2.536058 1.323976 0.590936
F3 2.367788 1.143582 0.703411
A3 2.971154 1.333221 0.654158
F4 3.144231 1.254722 0.598650
D2 2.850962 1.232843 0.532062
F5 2.242788 1.204888 0.564427
A4 2.853365 1.283871 0.466820
D3 3.466346 1.221819 0.523403
B1 2.199519 1.301010 0.496591
A5 2.432692 1.222397 0.626397
F6 3.134615 1.332143 0.650997
D4 2.605769 1.197735 0.622168
D5 2.290865 1.169703 0.621183
E1 3.144231 1.352682 0.604516
B2 2.206731 1.261307 0.535553
...
F10 1.848558 1.043002 0.618347
B5 1.701923 1.063115 0.520113
C6 3.278846 1.354736 0.473334
C7 1.819712 1.010163 0.542233このようにテキストP.37の表2.1の値が小数第2位の精度で一致していることが確認できます。
信頼性係数の算出
信頼性係数にはいろいろありますが、ここでは内的整合性の指標として一般的によく利用されているクロンバックの$${\alpha}$$係数とマクドナルドの$${\omega}$$係数を計算してみましょう。ただしPythonで計算する場合、$${\omega}$$係数を算出する便利なパッケージやメソッドはありませんので、自分で式を書く必要があることにご注意ください。
まず$${\alpha}$$係数についてみてみましょう。$${\alpha}$$係数の定義式は以下となります。
$$
\alpha = \frac{p}{p-1}\left(1-\sum_{i=1}^{p}\frac{s_i^2}{s_x^2}\right)
$$
ここで$${p}$$は項目の数、$${s_i^2}$$は各項目の分散、$${s_x^2}$$はテスト得点（本例の場合は心理的ストレス反応得点）の分散を表しています。この定義式に従い、以下のようにスクリプトを書いて計算してみましょう。
#α信頼性係数
p = statistic.shape[0]
s_i = statistic['std']
s_x = df['sum'].std()
(p/(p-1))*(1-((s_i**2).sum()/s_x**2))0.9375393314732781このように小数点第3位の精度でテキストの値と一致していることがわかります。パッケージを利用する場合、pingouinの中にあるcronbach_alpha()というメソッドを利用します。
pg.cronbach_alpha(df.drop('sum', axis=1))(0.9375393314732781, array([0.929, 0.946]))同じく小数点第3位の精度でテキストの値と一致していることがわかります。なお、[0.929, 0.946]の値は95％信頼区間を表しています。
このクロンバックの$${\alpha}$$係数ですが、古典的テストモデルの枠組みではタウ等化測定の仮定の下での信頼性係数になります。しかし、このタウ等化測定は「真の得点からの影響がすべて等しい」という少々現実的でない仮定になっています。
そこで、このタウ等化測定の制約を外し、同族測定の仮定のみで定義した信頼性係数がマクドナルドの$${\omega}$$係数になります。このあたりの説明をすると長くなるので今回は省略しますが、興味のある方はテキスト（あるいはテスト理論に関する専門書）をご参照ください。
テキストで紹介されている$${\omega}$$係数の定義式は以下となります。
$$
\omega = 1 - \sum_{i=1}^p \frac{d_i^2}{\sigma_x^2}
$$
ここで$${p}$$は項目の数、$${d_i^2}$$は因子分析における独自分散（独自性と観測変数の分散との積）、$${\sigma_x^2}$$は合計点（本例の場合は心理的ストレス反応得点）の分散を表しています。
実はこの$${\omega}$$係数は、時代が進むとともにいろいろな定義式が登場しました。提案したマクドナルド自身、1978年に最初に発表した定義式を、1999年には修正して再発表しています。またマクドナルド以外にも「レイコフの$${\omega}$$」など、様々なものが提案されています。本ブログではテキストに従い、1978年版のマクドナルドの$${\omega}$$係数を紹介しています。各$${\omega}$$係数の違いについて知りたい方はこちらをご参照ください。
上記の定義式からわかるように、マクドナルドの$${\omega}$$係数を計算するためには、先に因子分析を行う必要があります。そこで以下のようにして因子分析を実行してみましょう。
# 1因子モデルの因子分析の実行
fa = FactorAnalyzer(n_factors=1, rotation=None, method='ml')
fa.fit(df.drop('sum', axis=1))
loadings_df = pd.DataFrame(fa.loadings_, columns=["因子負荷量"])
loadings_df.index = df.drop('sum', axis=1).columns
loadings_df["共通性"] = fa.get_communalities()
loadings_df["独自性"] = fa.get_uniquenesses()
print(loadings_df) 因子負荷量 共通性 独自性
C1 0.392358 0.153945 0.846055
C2 0.340761 0.116118 0.883882
D1 0.439358 0.193035 0.806965
A1 0.610949 0.373259 0.626741
C3 0.494746 0.244773 0.755227
F1 0.566297 0.320692 0.679308
C4 0.467481 0.218538 0.781462
C5 0.493469 0.243511 0.756489
F2 0.730103 0.533050 0.466950
A2 0.594008 0.352846 0.647154
F3 0.723024 0.522764 0.477236
A3 0.660557 0.436335 0.563665
F4 0.588137 0.345905 0.654095
D2 0.517245 0.267542 0.732458
F5 0.560323 0.313962 0.686038
A4 0.452547 0.204799 0.795201
D3 0.503824 0.253839 0.746161
B1 0.450343 0.202809 0.797191
A5 0.624694 0.390243 0.609757
F6 0.635495 0.403854 0.596146
D4 0.600883 0.361060 0.638940
D5 0.600702 0.360843 0.639157
E1 0.582242 0.339006 0.660994
B2 0.496796 0.246806 0.753194
...
F10 0.625776 0.391595 0.608405
B5 0.481063 0.231421 0.768579
C6 0.439857 0.193474 0.806526
C7 0.509617 0.259709 0.740291次に1因子性が妥当かを確かめるためにスクリープロットを描いてみましょう。
# スクリーテスト
eigen = fa.get_eigenvalues()
left = list(range(1, len(eigen[0])+1))
plt.plot(left, eigen[0], marker="o")
plt.ylabel("固有値")
このようにスクリーテストの結果を見ても1因子性が妥当であることが確認できます。そこで、上記で紹介したマクドナルドの$${\omega}$$係数の定義式に従って、信頼性係数を算出してみましょう。
#ω信頼性係数
sum_d2 = loadings_df["独自性"].dot(statistic['std']**2)
sigma_x2 = (df['sum'].std())**2
1-(sum_d2/sigma_x2)0.9384119865276651このように、値がテキストと一致していることが確認できます。
最後に
いかがでしたでしょうか？第1回目～第3回目の記事で、Pythonを使っても基本的な因子分析ならRと同等に実行できることが分かったと思います。ではもっと特殊な場合の因子分析ではどうでしょうか？やはりどこまでPythonでできるのかが気になります。
次回の記事（因子分析4）では変数に順序カテゴリカルデータが含まれる場合の因子分析のやり方、すなわちカテゴリカル因子分析もPythonを使って実行できるということを見ていきたいと思います。
参考文献
McDonald, R. P. (1978). Generalizability in Factorable Domains: "Domain Validity and Generalizability". Educational and iPsychological Measurement, 38, 75-79.
McDonald, R. P. (1999). Test Theory: A unified treatment. Psychology Press.
この記事が気に入ったらサポートをしてみませんか？