因子分析

因子分析#

因子分析（factor analysis）はデータの分散を分解することによって何らかの共通性を抽出する手法である。

心理学などの分野において、「リーダシップ」や「社交性」のような直接観測することができない構成概念（construct）を抽出するために用いられる。

ビジネスにおいてはどういったところで応用されるのか？#

心理学以外の利用例

ビジネス関連の研究例

因子分析モデル#

1因子モデル#

\(p\)次元の観測値が\(n\)個得られているとし、\(i\)番目の\(j\)次元のデータを\(x_{ij}\)と表すことにする。\(x_{ij}\)を標準化して平均0、分散1にしたデータを\(z_{ij}\)と表すことにする。

1因子モデルは標準化データを次のように分解するモデルである。

1因子モデル（ベクトル表現）

\[ z_{ij} = a_j f_i + d_j u_{ij} \]

  • \(f_i\)は観測対象\(i\)の共通因子（common factor）で、構成概念を表す。

    • 共通因子の推定値を 共通因子スコア（common factor score） あるいは単に 因子スコア と呼ぶ。

  • \(a_j\)は観測変数\(j\)の因子負荷（factor loading）で、観測変数\(j\)が共通因子から受ける影響の程度を表す

  • \(u_{ij}\)は独自因子（unique factor）と呼ばれ、構成概念で説明しきれない要因を表す

    • 独自因子の推定値を 独自因子スコア（unique factor score） と呼ぶ。

  • \(d_j\)は観測変数\(j\)の独自係数と呼ばれ、観測変数\(j\)が独自因子から受ける影響の程度を表す

独自因子・独自係数は誤差項\(e\)としてまとめる表記

\[ z_{ij} = a_j f_i + e_{ij} \]

も存在する。独自因子と独自係数はどちらも自由に推定しようとすると解が定まらない（\(d_j\)を2倍にして\(u_{ij}\)を半分にしてもいい）ので、推定の際は通常は\(d_j=1\)に固定する。

1因子モデルをベクトルで表すと次のようになる

1因子モデル（ベクトル表現）

\[\begin{split} \begin{pmatrix} z_{1j}\\ z_{2j}\\ \vdots \\ z_{nj} \end{pmatrix} = a_j \begin{pmatrix} f_{1}\\ f_{2}\\ \vdots \\ f_{n} \end{pmatrix} + \begin{pmatrix} e_{1j}\\ e_{2j}\\ \vdots \\ e_{nj} \end{pmatrix} \\ \implies \mathbf{z}_j = a_j \mathbf{f} + \mathbf{e}_j \end{split}\]

因子スコア\(f_{i}\)が観測対象（アンケート回答者など）ごとに算出されるため、個々の観測対象のスコアリングができる。

\(m\)因子モデル#

\(m\)個の因子を扱うように拡張した場合は次のように表される。

\(m\)因子モデル

\[\begin{split} \begin{align} z_{ij} &= a_{j1} f_{i1} + \cdots + a_{jm} f_{im} + e_{ij}\\ &= \sum^m_{l=1} a_{jl} f_{il} + e_{ij}\\ \end{align} \end{split}\]

ベクトルの場合は

\[\begin{split} \begin{aligned} \mathbf{z}_j &= a_{j1} \mathbf{f}_1 + a_{j2} \mathbf{f}_2 + \cdots + a_{jm} \mathbf{f}_m + \mathbf{e}_j \\ &= \begin{pmatrix} \mathbf{f}_1 & \mathbf{f}_2 & \cdots & \mathbf{f}_m \end{pmatrix} \begin{pmatrix} a_{j1} \\ a_{j2} \\ \vdots \\ a_{jm} \end{pmatrix} + \mathbf{e}_j \\ &= \mathbf{F} \mathbf{a}_j + \mathbf{e}_j \end{aligned} \end{split}\]

と表すことができる。

因子分析モデル（行列表現）#

観測値ごとの添字をつけると\(i \times j \times l\)で3方向の次元（テンソル）になってしまうので簡略化のため観測値ではなく確率変数として扱うと、よりシンプルなモデルになる。

設問項目（item）が\(I\)個あり、\(I\)個の確率変数\(x_1,\dots,x_I\)が得られたとする。それを平均0分散1へと標準化したものを\(z_1,\dots,z_I\)とする。これらが\(T\)個の因子（trait）から影響を受けていると仮定し、\(T\)個の共通因子\(f_1,\dots,f_T\)と因子負荷量\(\lambda_{11},\lambda_{12},\dots,\lambda_{IT}\)で

\[ z_i = \lambda_{i1} f_1 + \cdots + \lambda_{iT} f_T + \varepsilon_i \quad (i = 1,\dots, I) \]

と説明する線形モデルが因子分析モデルとなる。行列とベクトルで表すと次のようになる

因子分析モデル

\[ \mathbf{z} = \mathbf{\Lambda} \mathbf{f} + \boldsymbol{\varepsilon} \]

  • \(\mathbf{z}\)：標準化した観測変数（\(I\)次元ベクトル）

  • \(\mathbf{\Lambda}\)：因子負荷行列（\(I\times T\)次元行列）

  • \(\mathbf{f}\)：共通因子ベクトル（\(T\)次元）

  • \(\boldsymbol{\varepsilon}\)：独自因子（誤差）ベクトル（\(I\)次元）

また、共通因子\(\mathbf{f}\)と独自因子\(\mathbf{\varepsilon}\)には次の仮定がおかれる：

  • \(\mathbf{f}\)と\(\mathbf{\varepsilon}\)は統計的に独立

  • \(\operatorname{E}[\mathbf{f}] = 0, \quad \operatorname{Var}[\mathbf{f}] = \mathbf{\Phi}\)

  • \(\operatorname{E}[\boldsymbol{\varepsilon}] = 0, \quad \operatorname{Var}[\boldsymbol{\varepsilon}] = \mathbf{\Psi}\)

    • （\(\mathbf{\Psi}\)は対角行列、つまり 異なる変数に対する独自因子は無相関とする）

この仮定により、観測変数の分散共分散行列を考えると

\[ \operatorname{Var}[\mathbf{z}] =\mathbf{\Sigma} =\mathbf{\Lambda \Phi \Lambda^\top} + \mathbf{\Psi} \]

となる。つまり、 観測された変数の分散共分散行列をモデルのパラメータで説明するモデル であることがわかる。

\(\mathbf{z}\)を中心化ではなく標準化している場合、分散共分散行列ではなく相関係数行列に対して上記の構造化を考えることになる。

（参考：因子分析 - Wikipedia）

パラメータの推定#

因子分析モデルのパラメータは観測変数の相関行列をもとにして最尤推定法などで推定される。

\(\mathbf{\Lambda}, \mathbf{\Psi}\)の推定（最尤推定法による例）#

標本共分散行列\(\mathbf{S}\)を使って、尤度関数は

\[ \mathcal{L}(\mathbf{\Lambda}, \mathbf{\Psi})=-\frac{n}{2}(\log |\mathbf{\Sigma}|+\operatorname{tr}(\mathbf{S} \mathbf{\Sigma}^{-1})) \]

と表される。これをEMアルゴリズムなどで最大化して\(\mathbf{\Lambda}, \mathbf{\Psi}\)の推定値を得る。

\(\mathbf{f}\)の推定（バートレット法による例）#

最小二乗法にちょっと似てる数式で、独自因子の分散の逆行列で重み付け（誤差が大きい場合はスコアへの寄与が小さくなる）を行いつつ推定する

\[ \hat{\mathbf{f}}=\left(\boldsymbol{\Lambda}^{\top} \mathbf{\Psi}^{-1} \boldsymbol{\Lambda}\right)^{-1} \boldsymbol{\Lambda}^{\top} \mathbf{\Psi}^{-1} \mathbf{x} \]

実装例#

import pandas as pd

# テストの点数データ
scores = pd.DataFrame([
    {"英語": 98, "国語": 95, "数学": 75},
    {"英語": 90, "国語": 75, "数学": 70},
    {"英語": 60, "国語": 57, "数学": 80},
    {"英語": 50, "国語": 70, "数学": 60},
    {"英語": 30, "国語": 50, "数学": 40},
])
Z = (scores - scores.mean()) / scores.std()

from statsmodels.multivariate.factor import Factor

fa = Factor(Z, n_factor=1, method="ml").fit()
fa.summary()

/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/statsmodels/multivariate/factor.py:417: UserWarning: Fitting did not converge
  warnings.warn("Fitting did not converge")
Communality
英語 国語 数学
0.9977 0.7723 0.5347
Pre-rotated loadings
factor 0
英語 0.9995
国語 0.8719
数学 0.7180 

fa.uniqueness

array([0.00228581, 0.22772564, 0.46534024])

# 共通因子スコア
fa.factor_scoring(Z)

array([[ 1.15171807],
       [ 0.85998399],
       [-0.19932969],
       [-0.54814633],
       [-1.26422604]])

主成分分析との違い#

主成分分析も分散を分解して共通した成分を取り出す手法であり、因子分析と似ているが、誤差を分離するかどうかが異なる。

因子分析は観測変数を共通した部分と誤差とに分解した上で共通因子を抽出するのに対し、主成分分析は誤差の分離を行わずに共通した部分を取り出そうとする。

妥当性と信頼性#

信頼性係数#

真の得点の分散の、テスト得点の分散に対する比を、テスト得点\(y\)の 信頼性係数 (reliability coefficient) という。

\[ \rho=\frac{\sigma_t^2}{\sigma_y^2} =\frac{\sigma_t^2}{\sigma_t^2+\sigma_e^2}=1-\frac{\sigma_e^2}{\sigma_y^2} \]

信頼性係数は測定の誤差分散が小さいほど1に近い値をとり、誤差分散が大きいとき0に近い値をとる。

クロンバックの\(\alpha\)係数#

信頼性係数\(\rho\)は母集団レベルの話であり、標本の情報から推定するには別の値を使う。 信頼性係数の推定方法はいくつかあるが、「項目を超えた一貫性」に関する検討を行うクロンバックのα係数が代表的。

「項目を超えた一貫性」とは「テストに含まれる項目群を、同様の別の項目群に置き換えても結果は安定しているか」という考え方で、そういう意味での信頼性を近似的に与える指標がクロンバックのα係数である。

クロンバックの\(\alpha\)
\[ \alpha=\frac{J}{J-1}\left(1-\frac{\sum_{j=1}^J s_j^2}{s_y^2}\right) \]

  • \(J\)：項目数

  • \(s_y^2\)：テスト得点の分散

  • \(s_j^2\)：項目\(j\)の分散

\(\alpha\)係数は

  • 項目間の相関が高いとき

  • 項目数が多いとき

に高い値をとる。項目が多いだけでも高くなってしまうので注意が必要で参考程度に用いられる。

測定の標準誤差（SEM）は信頼性係数の推定値\(\hat\rho\)を用いて推定される。

測定の標準誤差（SEM）の推定量
\[ \operatorname{SEM} = s_y \sqrt{1-\hat{\rho}} \]

Note

こうしたCTTの信頼性係数や測定の標準誤差の注意点は

  • 1つのテスト & 1つの受験者集団ごとに算出される

  • 異なる受験者集団に対して測定の標準誤差を求めると、異なる値が得られる（標本依存性）

Messickの妥当性概念#

Messickは妥当性検証のための側面を6つ挙げた。

側面（facet）

内容（簡潔な説明）

① 内容的妥当性（Content）

測定内容が理論的構成概念をどれだけ網羅しているか

② 構成概念妥当性（Substantive）

測定が理論に基づき、認知的に妥当か

③ 構造的妥当性（Structural）

得点構造が想定した構成概念と整合しているか（因子構造など）

④ 外的妥当性（External）

他の指標との関係が理論通りか（収束的・弁別的妥当性）

⑤ 結果的妥当性（Consequential）

測定の利用が社会的に望ましい結果をもたらすか（公正性・偏りなど）

⑥ 実証的妥当性（Generalizability）

測定結果が一般化可能か（状況・集団を超えて妥当か）

上記6つの側面にわたる証拠が備えられていることが望ましい。しかし妥当性検証とは尺度が利用されることで明らかになる証拠を積み重ねるプロセスであるため、最初から全部を備えることは不可能。

探索的因子分析#

探索的因子分析（exploratory factor analysis: EFA）は、因子と観測変数の関係についての仮説や制約を置かずに、観測変数のみから相関係数を計算し、観測変数間に相関関係をもたらす因子を推定する方法。

たとえばモデルの因子の数は固有値から探索して決める。例えば「1より大きい固有値をもつ因子を採用する」とする（Guttman基準）

サンプルデータを生成して因子分析を実行してみる#

Hide code cell source

 
import numpy as np
import pandas as pd

def generate_factor_analysis_data(n_samples=1000, n_variables=3, noise_std=0.1, random_state=42):
    """因子分析のサンプルデータを生成する"""
    np.random.seed(random_state)
    n_factors = 1 # 1因子モデル

    # 因子負荷量行列（loadings）: n_variables × n_factors
    # loadings = np.random.uniform(0.5, 0.9, size=(n_variables, n_factors))
    loadings = np.array([
        [0.9],
        [0.5],
        [0.1],
    ])

    # 潜在因子スコア（factors）: n_samples × n_factors
    factors = np.random.normal(size=(n_samples, n_factors))

    # ノイズ項（ユニーク性）： n_samples × n_variables
    noise = np.random.normal(scale=noise_std, size=(n_samples, n_variables))

    # 観測変数: n_samples × n_variables
    data = factors @ loadings.T + noise

    # データフレームに変換
    columns = [f"X{i+1}" for i in range(n_variables)]
    df = pd.DataFrame(data, columns=columns)
    return df, factors, loadings

# サンプルデータ生成
df, true_factors, true_loadings = generate_factor_analysis_data()

display(pd.DataFrame(true_factors).head().style.set_caption("真の因子スコアの最初の5行"))
display(pd.DataFrame(true_loadings).style.set_caption("因子負荷量"))
display(df.head().style.set_caption("サンプルデータの最初の5行"))
真の因子スコアの最初の5行
  0
0 0.496714
1 -0.138264
2 0.647689
3 1.523030
4 -0.234153
因子負荷量
  0
0 0.900000
1 0.500000
2 0.100000
サンプルデータの最初の5行
  X1 X2 X3
0 0.586978 0.340820 0.055634
1 -0.189132 0.000690 0.025522
2 0.672439 0.387361 0.169724
3 1.317203 0.893254 0.172063
4 -0.003212 -0.185995 0.150181

statsmodelsパッケージでの推定#

from statsmodels.multivariate.factor import Factor
fa = Factor(endog=df, n_factor=1, method="ml").fit()
fa.rotate("promax")
fa.summary()

/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/statsmodels/multivariate/factor_rotation/_gpa_rotation.py:101: RuntimeWarning: divide by zero encountered in log10
  table.append([i_try, f, np.log10(s), al])
Communality
X1 X2 X3
0.9821 0.9655 0.4711
Pre-rotated loadings
factor 0
X1 0.9910
X2 0.9826
X3 0.6863
promax rotated loadings
factor 0
X1 0.9910
X2 0.9826
X3 0.6863 

factors_ = fa.factor_scoring(df)
display(pd.DataFrame(factors_).head().style.set_caption("因子スコアの最初の5行"))
因子スコアの最初の5行
  0
0 0.653819
1 -0.160928
2 0.762085
3 1.589657
4 -0.136692

factor-analyzerパッケージでの推定#

from factor_analyzer import FactorAnalyzer
fa = FactorAnalyzer(n_factors=1, rotation='promax', method='ml', is_corr_matrix=False)
fa.fit(df)

# 因子負荷量（loadings）
loadings = pd.DataFrame(fa.loadings_, index=df.columns, columns=["第1因子"])
loadings.style.set_caption("因子負荷量")

/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/sklearn/utils/deprecation.py:151: FutureWarning: 'force_all_finite' was renamed to 'ensure_all_finite' in 1.6 and will be removed in 1.8.
  warnings.warn(
/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/factor_analyzer/factor_analyzer.py:663: UserWarning: No rotation will be performed when the number of factors equals 1.
  warnings.warn(
因子負荷量
  第1因子
X1 -0.991009
X2 -0.982624
X3 -0.686349 
# 相関行列を渡す場合
from factor_analyzer import FactorAnalyzer
fa = FactorAnalyzer(n_factors=1, rotation='promax', method='ml', is_corr_matrix=True)
fa.fit(df.corr())

# 因子負荷量（loadings）
loadings = pd.DataFrame(fa.loadings_, index=df.columns, columns=["第1因子"])
loadings.style.set_caption("因子負荷量")

/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/sklearn/utils/deprecation.py:151: FutureWarning: 'force_all_finite' was renamed to 'ensure_all_finite' in 1.6 and will be removed in 1.8.
  warnings.warn(
/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/factor_analyzer/factor_analyzer.py:663: UserWarning: No rotation will be performed when the number of factors equals 1.
  warnings.warn(
因子負荷量
  第1因子
X1 0.991009
X2 0.982624
X3 0.686349

固有値を高い順に並べる スクリープロット

Hide code cell source

 
# スクリープロットで因子数を探索
import matplotlib.pyplot as plt
import matplotlib_fontja

eigen = fa.get_eigenvalues()
left = list(range(1, len(eigen[0])+1))
fig, ax = plt.subplots()
ax.plot(left, eigen[0], marker="o")
ax.set(xlabel="因子数（固有値の数）", ylabel="固有値", title="スクリープロット")
ax.axhline(y=1, color="gray", linestyle="--", alpha=0.5)
fig.show()
../_images/04f99c49c164655d947261be43f447917feffe84506bc8ec12dd2ca4bc6497f1.png

scikit-learnの因子分析#

# 2因子モデルで因子分析を実行（scikit-learn）
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_std = scaler.fit_transform(df)

from sklearn.decomposition import FactorAnalysis
fa = FactorAnalysis(n_components=1, random_state=0, rotation="varimax")
fa.fit(X_std)

# 因子スコア
factors_ = fa.transform(X_std)

# 因子負荷量（loadings）
loadings_ = pd.DataFrame(fa.components_.T, index=df.columns, columns=["第1因子"])  # shape: (n_variables, n_factors)
loadings_.style.set_caption("因子負荷量")
因子負荷量
  第1因子
X1 0.987525
X2 0.986081
X3 0.686483 
plt.scatter(true_factors, factors_)

<matplotlib.collections.PathCollection at 0x7f0cdc628130>
../_images/0cdefad8e764c7e069456ed63c83ce42f861b03a95448cf4fef73435591170c9.png

確認的因子分析#

確認的因子分析（confirmatory factor analysis: CFA）は、因子の数や観測変数との関係性などについての仮説（モデル）をあらかじめ立てておき、その仮説の正しさを検証するためにモデルをデータにあてはめていく。

その妥当性はデータの分散をモデルがどれだけ説明できたか（適合度）の指標などによって確認される。

import semopy

カテゴリカル因子分析#

通常の因子分析では、観測変数も連続値（正規分布を仮定）であることを前提にしている。しかし実際の調査データでは、

  • 「1（全くそう思わない）」〜「5（非常にそう思う）」のリッカート尺度

  • 「ある／ない」「はい／いいえ」の2値データ

  • 評価カテゴリ（例：低・中・高）

といった順序尺度のカテゴリカル変数が含まれることがある。

因子分析では相関行列をもとに計算できるため、カテゴリカル変数であっても相関関係へと変換してから扱えば通常の計算手続きに持ち込むことができる。

カテゴリカル変数の相関関係には「観測値はカテゴリカル変数だが、その背後には連続値の潜在変数があり、ある閾値をこえるかどうかで離散的な観測値になっている」という考え方で統計モデルを組んで相関係数を算出するタイプのものが存在する。この仮定は因子分析モデルの仮定（離散値の背後には潜在的な因子スコアが存在する）とも整合的なのでそういった相関係数が使われる。具体的には、

  • 連続変数 × カテゴリカル変数 の相関 → ポリシリアル（polyserial）相関係数

  • カテゴリカル変数 × カテゴリカル変数 の相関 → ポリコリック（polychoric）相関係数

が使われる。Pythonだと ordinalcorr パッケージで計算できる。

実装例#

先程のサンプルデータを離散化してみる

Hide code cell source

 
def discretize_to_categorical(df, n_bins=4, labels=None):
    """
    各列をカテゴリカルに変換（等頻度ビンに分割）
    
    Parameters:
        df: pd.DataFrame - 入力の連続値データ
        n_bins: int - ビンの数（カテゴリ数）
        labels: list or None - カテゴリラベル。Noneなら自動で0～(n_bins-1)

    Returns:
        pd.DataFrame - カテゴリカル変数データ
    """
    if labels is None:
        labels = list(range(n_bins))

    df_cat = df.copy()
    for col in df.columns:
        df_cat[col] = pd.qcut(df[col], q=n_bins, labels=labels)
    return df_cat

# カテゴリカルに変換
df_cat = discretize_to_categorical(df)
df_cat.head().style.set_caption("サンプルデータの最初の5行")
サンプルデータの最初の5行
  X1 X2 X3
0 2 2 2
1 1 1 2
2 3 3 3
3 3 3 3
4 1 1 3

相関行列に変換する

from ordinalcorr import hetcor
corr_matrix = hetcor(df_cat)
corr_matrix
X1 X2 X3
X1 1.000000 0.970696 0.679743
X2 0.970696 1.000000 0.682401
X3 0.679743 0.682401 1.000000

FactorAnalyzerでis_corr_matrix=Trueを指定して相関行列を入力データとして渡すようにする

from factor_analyzer import FactorAnalyzer
fa = FactorAnalyzer(n_factors=1, rotation='promax', method='ml', is_corr_matrix=True)
fa.fit(corr_matrix)

# 因子負荷量（loadings）
loadings_cat = pd.DataFrame(fa.loadings_, index=df.columns, columns=["第1因子"])
loadings_cat.style.set_caption("因子負荷量（カテゴリカル）")

/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/sklearn/utils/deprecation.py:151: FutureWarning: 'force_all_finite' was renamed to 'ensure_all_finite' in 1.6 and will be removed in 1.8.
  warnings.warn(
/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/factor_analyzer/factor_analyzer.py:663: UserWarning: No rotation will be performed when the number of factors equals 1.
  warnings.warn(
因子負荷量（カテゴリカル）
  第1因子
X1 -0.983318
X2 -0.987164
X3 -0.691274 
factors_cat = fa.transform(df)
plt.scatter(true_factors, factors_cat)

/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/sklearn/utils/deprecation.py:151: FutureWarning: 'force_all_finite' was renamed to 'ensure_all_finite' in 1.6 and will be removed in 1.8.
  warnings.warn(
/home/runner/work/notes/notes/.venv/lib/python3.10/site-packages/factor_analyzer/factor_analyzer.py:766: UserWarning: Could not find original mean and standard deviation; usingthe mean and standard deviation from the current data set.
  warnings.warn(

<matplotlib.collections.PathCollection at 0x7f0cd8798490>
../_images/e7bb3f4507130cb24aef4718af17ee6425e9ded2f5c4f00e355be8066c02923c.png

カテゴリカル変数にする前の連続変数に対する推定結果とある程度近い推定値となっていることがわかる

(loadings_cat - loadings).abs().style.set_caption("連続変数の場合との推定値の絶対誤差")
連続変数の場合との推定値の絶対誤差
  第1因子
X1 1.974327
X2 1.969787
X3 1.377623

Hide code cell source

 
# スクリープロットで因子数を探索
eigen = fa.get_eigenvalues()
left = list(range(1, len(eigen[0])+1))
fig, ax = plt.subplots()
ax.plot(left, eigen[0], marker="o")
ax.set(xlabel="因子数（固有値の数）", ylabel="固有値", title="スクリープロット")
ax.axhline(y=1, color="gray", linestyle="--", alpha=0.5)
fig.show()
../_images/d341fadfcb7378b833e31636f61f24f87a9788cb669beb9fd79ef2d5be64e51c.png

例2：実際のアンケートデータ#

bfi データ（Big Five Inventory）は、心理学でよく使われる**5因子性格モデル（OCEAN）**に基づいたアンケートデータ

25個の性格項目（各5問 × 5因子）

  • A1 ～ A5：Agreeableness（協調性）

  • C1 ～ C5：Conscientiousness（誠実性）

  • E1 ～ E5：Extraversion（外向性）

  • N1 ～ N5：Neuroticism（神経症傾向）

  • O1 ～ O5：Openness（開放性）

import statsmodels.api as sm

# psychパッケージの Big Five Inventory（性格5因子） データ
df = sm.datasets.get_rdataset("bfi", package="psych").data
df.head()
A1 A2 A3 A4 A5 C1 C2 C3 C4 C5 ... N4 N5 O1 O2 O3 O4 O5 gender education age
rownames
61617 2.0 4.0 3.0 4.0 4.0 2.0 3.0 3.0 4.0 4.0 ... 2.0 3.0 3.0 6 3.0 4.0 3.0 1 NaN 16
61618 2.0 4.0 5.0 2.0 5.0 5.0 4.0 4.0 3.0 4.0 ... 5.0 5.0 4.0 2 4.0 3.0 3.0 2 NaN 18
61620 5.0 4.0 5.0 4.0 4.0 4.0 5.0 4.0 2.0 5.0 ... 2.0 3.0 4.0 2 5.0 5.0 2.0 2 NaN 17
61621 4.0 4.0 6.0 5.0 5.0 4.0 4.0 3.0 5.0 5.0 ... 4.0 1.0 3.0 3 4.0 3.0 5.0 2 NaN 17
61622 2.0 3.0 3.0 4.0 5.0 4.0 4.0 5.0 3.0 2.0 ... 4.0 3.0 3.0 3 4.0 3.0 3.0 1 NaN 17

5 rows × 28 columns

from statsmodels.multivariate.factor import Factor
fa = Factor(endog=df.iloc[:, 0:25], n_factor=5, method="ml").fit()
fa.rotate("varimax")
fa.get_loadings_frame()
  factor 0 factor 1 factor 2 factor 3 factor 4
E2 -0.670905 0.186767 -0.094653 0.206146 -0.097616
E4 0.581940 -0.106014 0.067780 -0.421312 0.000625
E1 -0.564269 -0.003799 0.040074 0.183339 -0.110895
E5 0.452352 0.073793 0.304980 -0.190193 0.262146
E3 0.425604 0.010923 0.060131 -0.369124 0.348025
N1 0.028616 0.844284 -0.046915 0.104798 -0.050977
N2 -0.021902 0.807156 -0.023921 0.102075 0.010221
N3 -0.159165 0.697670 -0.082174 -0.055294 0.019294
N4 -0.432049 0.516176 -0.187678 -0.012137 0.068812
N5 -0.245634 0.483898 -0.059355 -0.141936 -0.132161
C4 -0.105730 0.210326 -0.654757 0.026594 -0.072408
C2 -0.017581 0.057436 0.623406 -0.153290 0.125050
C5 -0.221359 0.252726 -0.569024 0.062313 0.047009
C3 0.009739 -0.041651 0.549671 -0.136527 -0.012136
C1 0.034580 -0.007740 0.536762 -0.083552 0.209220
A3 0.197586 -0.038162 0.086844 -0.691715 0.085198
A2 0.115756 -0.024847 0.124441 -0.623642 0.072931
A5 0.291354 -0.166135 0.057848 -0.601124 0.104497
A4 0.145934 -0.093663 0.213993 -0.471772 -0.102302
A1 0.079856 0.159325 0.013873 0.368451 -0.051519
O3 0.211458 0.004002 0.074508 -0.179411 0.631765
O1 0.136184 -0.020774 0.113180 -0.101021 0.532319
O5 0.016270 0.090103 -0.094123 -0.027126 -0.506041
O2 -0.003365 0.163857 -0.130711 -0.121798 -0.443965
O4 -0.281219 0.150597 -0.022332 -0.132553 0.358567 
# シンプルにしてみる
df = df.iloc[:, 0:10].dropna()

from statsmodels.multivariate.factor import Factor
fa = Factor(endog=df, n_factor=2, method="ml").fit()
fa.rotate("varimax")
fa.get_loadings_frame()
  factor 0 factor 1
A3 -0.757357 0.079160
A2 -0.652163 0.112862
A5 -0.617412 0.111762
A4 -0.466231 0.210560
A1 0.381575 0.009641
C4 0.106303 -0.657542
C2 -0.101881 0.621765
C5 0.146665 -0.574452
C1 -0.056792 0.555669
C3 -0.122728 0.548047 
df.corr().style.format("{:.3f}").background_gradient()
  A1 A2 A3 A4 A5 C1 C2 C3 C4 C5
A1 1.000 -0.343 -0.266 -0.146 -0.183 0.020 0.015 -0.021 0.121 0.045
A2 -0.343 1.000 0.490 0.341 0.388 0.097 0.134 0.191 -0.145 -0.125
A3 -0.266 0.490 1.000 0.365 0.506 0.100 0.138 0.124 -0.124 -0.152
A4 -0.146 0.341 0.365 1.000 0.313 0.097 0.231 0.133 -0.157 -0.245
A5 -0.183 0.388 0.506 0.313 1.000 0.128 0.115 0.134 -0.125 -0.168
C1 0.020 0.097 0.100 0.097 0.128 1.000 0.437 0.316 -0.350 -0.259
C2 0.015 0.134 0.138 0.231 0.115 0.437 1.000 0.367 -0.391 -0.304
C3 -0.021 0.191 0.124 0.133 0.134 0.316 0.367 1.000 -0.351 -0.347
C4 0.121 -0.145 -0.124 -0.157 -0.125 -0.350 -0.391 -0.351 1.000 0.482
C5 0.045 -0.125 -0.152 -0.245 -0.168 -0.259 -0.304 -0.347 0.482 1.000 
# 相関係数をとってみる
from ordinalcorr import hetcor
corr = hetcor(df)
corr.style.format("{:.3f}").background_gradient()
  A1 A2 A3 A4 A5 C1 C2 C3 C4 C5
A1 1.000 -0.410 -0.324 -0.173 -0.230 -0.006 0.007 -0.024 0.145 0.055
A2 -0.410 1.000 0.561 0.396 0.446 0.119 0.157 0.224 -0.191 -0.159
A3 -0.324 0.561 1.000 0.414 0.574 0.128 0.155 0.154 -0.168 -0.186
A4 -0.173 0.396 0.414 1.000 0.361 0.117 0.269 0.169 -0.199 -0.286
A5 -0.230 0.446 0.574 0.361 1.000 0.165 0.141 0.154 -0.164 -0.199
C1 -0.006 0.119 0.128 0.117 0.165 1.000 0.491 0.353 -0.411 -0.302
C2 0.007 0.157 0.155 0.269 0.141 0.491 1.000 0.408 -0.444 -0.337
C3 -0.024 0.224 0.154 0.169 0.154 0.353 0.408 1.000 -0.394 -0.389
C4 0.145 -0.191 -0.168 -0.199 -0.164 -0.411 -0.444 -0.394 1.000 0.535
C5 0.055 -0.159 -0.186 -0.286 -0.199 -0.302 -0.337 -0.389 0.535 1.000 
fa = Factor(corr=corr, n_factor=2, method="ml").fit()
fa.rotate("varimax")
fa.get_loadings_frame()
  factor 0 factor 1
var2 -0.802822 0.094423
var1 -0.700661 0.129877
var4 -0.659263 0.130041
var3 -0.494776 0.237330
var0 0.435057 0.006267
var8 0.137793 -0.697534
var6 -0.105498 0.657520
var9 0.165186 -0.603882
var5 -0.072780 0.598489
var7 -0.136864 0.574023 
# カテゴリの数を減らしてみる
df_cat = df.copy()
for col in df.columns:
    df_cat[col] = pd.cut(df[col], bins=3, labels=range(3)).astype(int)

# 相関係数をとってみる
from ordinalcorr import hetcor
corr = hetcor(df_cat, n_unique=50)
corr.style.format("{:.3f}").background_gradient()
  A1 A2 A3 A4 A5 C1 C2 C3 C4 C5
A1 1.000 -0.416 -0.298 -0.194 -0.215 0.035 0.016 -0.028 0.150 0.057
A2 -0.416 1.000 0.574 0.448 0.454 0.098 0.174 0.238 -0.161 -0.155
A3 -0.298 0.574 1.000 0.454 0.589 0.122 0.172 0.172 -0.138 -0.192
A4 -0.194 0.448 0.454 1.000 0.362 0.112 0.287 0.146 -0.172 -0.280
A5 -0.215 0.454 0.589 0.362 1.000 0.151 0.107 0.150 -0.131 -0.191
C1 0.035 0.098 0.122 0.112 0.151 1.000 0.499 0.374 -0.423 -0.291
C2 0.016 0.174 0.172 0.287 0.107 0.499 1.000 0.435 -0.423 -0.338
C3 -0.028 0.238 0.172 0.146 0.150 0.374 0.435 1.000 -0.387 -0.399
C4 0.150 -0.161 -0.138 -0.172 -0.131 -0.423 -0.423 -0.387 1.000 0.537
C5 0.057 -0.155 -0.192 -0.280 -0.191 -0.291 -0.338 -0.399 0.537 1.000 
fa = Factor(corr=corr, n_factor=2, method="ml").fit()
fa.rotate("varimax")
fa.get_loadings_frame()
  factor 0 factor 1
var2 -0.806967 0.095949
var1 -0.719831 0.121336
var4 -0.663792 0.098485
var3 -0.544785 0.214299
var0 0.415463 0.014423
var8 0.107723 -0.683504
var6 -0.123003 0.662741
var5 -0.056390 0.619970
var9 0.171923 -0.593443
var7 -0.151973 0.590409

参考#

Contents