Pythonで潜在クラス分析ができる！Stepmixの使い方

# はじめに
カテゴリカルな変数で記述されるデータがあるとき、どのようなクラス（グループ）があるかを確認したい場合、一つの方法として **潜在クラス分析** (LCA) と呼ばれる統計手法があります。
以前はこの手法を実行する場合、RではpoLCAなどのライブラリがあり比較的簡単に実行できましたが、Pythonではライブラリがありませんでした（少なくとも私は知らなかった）。しかし、2023年に `Stepmix` というライブラリが公開されて、Pythonでも実行できるようになりました[1]。
アンケートなどでカテゴリカルな変数で構成されるデータがある場合、これからはPythonでも潜在クラス分析のライブラリが利用できることは便利だと思います。この記事では、簡単な例を使って、Pythonでの潜在クラス分析の実施方法を紹介します。

# **1. 潜在クラス分析（LCA）とは**
潜在クラス分析（Latent Class Analysis, LCA）は、観測データから潜在的なクラス（グループ）を推定する統計手法です。LCAはカテゴリカルデータの背後にある構造を明らかにし、マーケティング、心理学、医療統計などで広く活用されています。

## 潜在クラス分析の概要

潜在クラス分析（Latent Class Analysis: LCA）とは、観測されたデータに基づいて、観測されなかった潜在的なクラスターを推定する統計的手法です。
アンケート調査などで得られた回答データから、回答者の背後に存在する可能性のあるグループ（潜在クラス）を特定するために使われます。

## 潜在クラス分析とクラスタリングの違い

潜在クラス分析と似た手法として、クラスタリングがあります。
どちらも与えられたデータに基づいてグループ分けを行うという点で共通していますが、いくつかの重要な違いがあります。

| 項目 | 潜在クラス分析 | クラスタリング |
| --- | --- | --- |
| モデル | 確率モデル | 距離に基づくモデル |
| データの仮定 | 確率分布に従う | 特になし |
| カテゴリカル変数 | 扱える | 扱えない場合が多い |
| 目的 | 潜在クラスの推定 | 類似度の高いデータのグループ化 |
| メリット | 統計的基盤がしっかりしている | 計算が比較的容易 |
| デメリット | 計算に時間がかかる場合がある | 結果の解釈が難しい場合がある |

潜在クラス分析は、各データがどの潜在クラスに属するかを確率的に推定します。
一方、クラスタリングは、データ間の距離や類似度に基づいてデータをグループ化します。
また、k-meansなど多くのクラスタリング手法は、カテゴリカル変数を直接扱うことができません。
潜在クラス分析は、カテゴリカル変数を自然に扱うことができます。

クラスタリングの手法でカテゴリカル変数を扱う場合は、K-Prototypeという方法があります。この方法については拙著の[こちらの記事](https://qiita.com/shinji_komine/items/b634146bd873d8d876e5)をご参照ください。

# **2. StepMixによるLCAの実装**
Pythonの`StepMix`は、潜在クラス分析を行うためのライブラリで、EMアルゴリズムを用いてパラメータ推定を行います。

## **2.1 データの準備**
今回の分析では、Rで潜在クラス分析を行う `poLCA` ライブラリに入っている `cheating` というデータセットを使います。同じデータセットを使うことで、後でRの結果を再現できているかを念のために確認します。

`cheating` データセットは、Rの`poLCA`パッケージに含まれているサンプルデータで、319人の大学生を対象に、カンニング行為に関する4つの質問への二値回答と、各学生の学業成績（GPA）を記録しています。  

| 変数名   | 内容                                                                                     | 値の範囲                                                                                 |
|-----------|------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------|
| `LIEEXAM` | 試験を受けるのを避けるために嘘をついたことがあるか                                         | 1: いいえ<br>2: はい                                                                     |
| `LIEPAPER`| レポートの提出期限を延ばすために嘘をついたことがあるか                                     | 1: いいえ<br>2: はい                                                                     |
| `FRAUD`   | 自分のものとして提出するためにレポートを購入した、または試験前に試験のコピーを入手したことがあるか | 1: いいえ<br>2: はい                                                                     |
| `COPYEXAM`| 試験中に近くの人の答えをコピーしたことがあるか                                             | 1: いいえ<br>2: はい                                                                     |
| `GPA`     | 学業成績（GPA）                                                                           | 1: 2.99以下<br>2: 3.00-3.25<br>3: 3.26-3.50<br>4: 3.51-3.75<br>5: 3.76-4.00<br>欠損値あり |

なお、`GPA`のデータは、カンニングを報告しなかった4人の学生については欠損しています。

私はRでこのデータをcsvファイルに出力してから、Pythonへ読み込むことにしました。

```python
# cheatデータの読み込み
df_cheat = pd.read_csv('cheating.csv')
df_cheat = df_cheat.iloc[:,1:]
print(df_cheat.shape)
display(df_cheat.head())
```

データには複数のカテゴリカル変数が含まれ、それらをもとに潜在クラスを推定します。

## **2.2 StepMixによるモデル適用**
StepMixを使って、潜在クラス分析を実施します。

StepMixでは、ベルヌーイ分布に従う二値項目は0,1である必要があるので、変換しておきます。
```python
# ライブラリの読み込み
import pandas as pd
from stepmix.stepmix import StepMix

# 説明変数のリスト作成
exp_vars = df_cheat.columns.tolist()
exp_vars.remove('GPA')

# Stepmixでは、ベルヌーイ分布に従う二値項目は0,1である必要があるので、変換
for exp_var in exp_vars:
    df_cheat[exp_var] = df_cheat[exp_var] - 1
```
次に `StepMix` クラスを使って、潜在クラスのモデルを作成します。潜在クラスの数は2つと仮定します。
`n_components=2` として、潜在クラス数を2と設定します。

```python
# モデルの設定
model = StepMix(n_components=2, n_steps=3, measurement='bernoulli',  random_state=1)

# モデルの学習
model.fit(df_cheat[exp_vars])
```

EMアルゴリズムによってパラメータ推定が行われ、各観測データがどのクラスに属するかの確率が得られます。


## **2.3 結果の確認**

結果は、get_parameters_df()で確認できます。
```python
# 自由度
print(f'degree of freedom : {model.n_parameters}')

# 推定されたパラメータ
model.get_parameters_df()
```

![StepMix_get_param.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/225870/d7a40674-8f63-49ce-baec-14e15e0dc5e7.png)

一番下の２行が2つの潜在クラスに属する事前確率です。0と1の2つの潜在クラスのうち、潜在クラス1に属する確率が0.839と高くなっています。
上から16個のパラメータは、それぞれの潜在クラス0と1に属しているとき、yesあるいはnoと回答する確率が出ています。
パラメータは 2+16=18 と合計18個計算されています。ただし、確率は合計すると1という制約条件があるため、自由度は9になっています。

StepMixでは、モデルの適合度を評価するために、AIC（赤池情報量基準）、BIC（ベイズ情報量基準）、対数尤度（log-likelihood）を確認できます。

```python
# 情報基準統計量の確認
print(f'aic : {model.aic(df_cheat[exp_vars])}')
print(f'bic : {model.bic(df_cheat[exp_vars])}')

# 対数尤度の確認
print(f'average log-likelihood : { model.score(df_cheat[exp_vars])}' )
print(f'log-likelihood : { model.score(df_cheat[exp_vars]) * df_cheat.shape[0]}')
```
![StepMix_AIC.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/225870/b0ed6502-f7eb-4552-bee2-b02500a828d6.png)

## 2.4 **モデルにより推論された潜在クラスの確認**

作ったモデルで、データが潜在クラス0と1のどちらに割り当てられたかを、4つの説明変数のyes、noごとに件数を確認します。

```python
df_cheat2 = df_cheat.copy()
df_cheat2['pred_class'] = model.predict(df_cheat[exp_vars])

pd.pivot_table(data=df_cheat2, index=exp_vars, columns=['pred_class'], values='GPA', aggfunc='count', fill_value=0)
```
![StepMix_pivot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/225870/248f7e2e-a62d-450f-b52b-a3a577a42ecc.png)

LIEEXAMとLIEPAPERで「いいえ」と回答した人が潜在クラス1、これ以外の人が潜在クラス0と予測されています。きれいに分かれていますね。


## **3. RのpoLCAとの比較**
Rの `poLCA` ライブラリを使って、同じデータを使って、同じ潜在クラス分析を実施しました。

### **3.1 poLCAの実装（R）**
Rでの実装は以下のようになります。

```r
library("poLCA")

data(cheating)
dim(cheating); head(cheating)

# LCAモデルの設定
f <- cbind(LIEEXAM,LIEPAPER,FRAUD,COPYEXAM)~1
ch2 <- poLCA(f,data,nclass=2)

# 結果の表示
summary(lca_model)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/225870/96b422fc-3fb7-43b8-9d6f-6f6f6041a43e.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/225870/a1bf29fd-e616-4d2d-8dd1-aae6e1b97bfb.png)

推定されたパラメータ、AICとBICの情報量基準統計量、対数尤度がStepMixとpoLCAで一致しています。当たり前ですけど、実績のあるライブラリと一致したので、StepMixの使い方はこれで正しいはずです。これでPythonでも安心して潜在クラス分析を実行できます。


## **まとめ**
- 潜在クラス分析（LCA）は、カテゴリカルデータを基に潜在クラスを推定する手法です
- Pythonの`StepMix`を用いることで、LCAを実施できる
- モデルの推定結果、対数尤度、AICをStepMixで確認できる
- 結果は、無事RのpoLCAに一致することを確認できました


# 参考
[1] [StepMix: A Python Package for Pseudo-Likelihood Estimation of Generalized Mixture Models with External Variables, Sacha Morin et al, 2023](https://arxiv.org/abs/2304.03853)