Irisデータセットとstatsmodelsによるロジスティック回帰の実装

この記事では，インターネット上に公開されている [Iris plants dataset](https://archive.ics.uci.edu/dataset/53/iris) (以下Iris データセット) を例にロジスティック回帰分析をstatsmodels (Python) で実装する方法を説明します．

## 概要
Iris (アヤメ) データセットは，セトサ (Iris-Setosa)，バーシクル (Iris-Versicolor)，バージニカ (Iris-virginica) の3種類のアヤメの花の「がく片 (sepal)」と「花弁 (petal)」の「幅」と「長さ」が収められたデータセットです．

このデータを用いて，がく片の長さ ($x_1$)，がく片の幅 ($x_2$)，花弁の長さ ($x_3$)，花弁の幅 ($x_4$) からアヤメの花がバージニカ種かそれ以外か ($y$) を予測するためのモデルを多重ロジスティック回帰分析を用いて構築します．
ここで
* $y$ は反応変数で，バージニカ種のときは$y = 1$，そうでないときは$y = 0$です．
* $x_1, x_2, x_3, x_4$ は説明変数です．
* $x$ の前の数値は係数です．

この例で，アヤメの花がバージニカ種である確率を$p$とすると，ロジスティック回帰モデルでは$p$は次の式で表されます．
$$
p = \frac{1}{1 + \exp[-(\beta_0 + \beta_1 x_1 + \cdots + \beta_4 x_4)]}
$$
ここで$\beta_0$は切片です．この式は次のように書き換えることもできます．
$$
\log(p/(1-p)) = \beta_0 + \beta_1 x_1 + \cdots + \beta_4 x_4
$$

## セットアップ
はじめに，本説明で用いるモジュールをあらかじめインポートしておきます．

```python
import pandas as pd
import statsmodels.api as sm
```

## データをロードする
`scikit-learn`の`load_iris()`関数を用いてデータをロードします．`pandas`の`DataFrame`または`Series`型でロードするために，`as_frame`を`True`に設定しています．

```python
from sklearn.datasets import load_iris
iris = load_iris(as_frame=True)
iris.keys()
```
出力: 
```
dict_keys(['data', 'target', 'frame', 'target_names', 'DESCR', 'feature_names', 'filename', 'data_module'])
```

`iris.data`に観察したアヤメの花のデータが、`iris.target`にアヤメの花の分類データが格納されています．アヤメの花がセトサのときは0，バーシクルのときは1，バージニカのときは2で表現されています (`iris.target_names`のインデックスと対応しています)．
```python
iris_data = iris.data
iris_target = iris.target
```

## モデルを構築するためにデータを整形する
反応変数を`endog`, 説明変数を`exog`とします．`iris_target`には3種それぞれに数値が与えられているので，これを「バージニカ種のときは1，それ以外は0」となるように変換します．statsmodels の統計モデルにデータを渡すとき，デフォルトでは切片は含まれていないので，`add_constant()` 関数で切片項を追加する必要があります．
```python
endog = (iris_target == 2)
exog = sm.add_constant(iris_data)
```

## 単ロジスティック回帰分析
上記で定義した`logit_res`に対して`summary()`メソッドを呼び出すことで，各説明変数に対する係数，標準偏差，Wald統計検定量，P値，95%信頼区間などの統計報告を得ることができます．

詳細は，[statsmodels.discrete.discrete_model.LogitResults.summary - statsmodels 0.14.0](https://www.statsmodels.org/stable/generated/statsmodels.discrete.discrete_model.LogitResults.summary.html)をご覧ください．

まずはそれぞれの特徴量について単ロジスティック回帰分析を行い、$P \leq 0.1$である有意な特徴量について探索します。

```python
from io import StringIO

print("有意な特徴量:\n")
# const以外のcolumnをそれぞれattrとして反復する
for attr in exog.columns[1:]:
    logit_mod = sm.Logit(endog, exog[["const", attr]])
    logit_res = logit_mod.fit()
    summary = logit_res.summary()

    table = pd.read_csv(StringIO(summary.tables[1].as_csv()), index_col=0)
    p = table.at[attr, "P>|z| "]
    if p <= 0.1:
        print(f"{attr} (P = {p})")
```
出力: 
```
有意な特徴量:

Optimization terminated successfully.
         Current function value: 0.391152
         Iterations 7
sepal length (cm) (P = 0.0)
Optimization terminated successfully.
         Current function value: 0.627128
         Iterations 5
sepal width (cm) (P = 0.1)
Optimization terminated successfully.
         Current function value: 0.111440
         Iterations 11
petal length (cm) (P = 0.0)
Optimization terminated successfully.
         Current function value: 0.111403
         Iterations 10
petal width (cm) (P = 0.0)
```

## 重ロジスティック回帰分析
単ロジスティック回帰分析の結果をもとに、重ロジスティック回帰モデルを作成します。

```python
logit_mod = sm.Logit(endog, exog)
logit_res = logit_mod.fit()
summary = logit_res.summary()
print(summary)
```
出力: 
```
Optimization terminated successfully.
         Current function value: 0.039662
         Iterations 14
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                 target   No. Observations:                  150
Model:                          Logit   Df Residuals:                      145
Method:                           MLE   Df Model:                            4
Date:                Fri, 22 Sep 2023   Pseudo R-squ.:                  0.9377
Time:                        22:50:38   Log-Likelihood:                -5.9493
converged:                       True   LL-Null:                       -95.477
Covariance Type:            nonrobust   LLR p-value:                 1.189e-37
=====================================================================================
                        coef    std err          z      P>|z|      [0.025      0.975]
-------------------------------------------------------------------------------------
const               -42.6378     25.708     -1.659      0.097     -93.024       7.748
sepal length (cm)    -2.4652      2.394     -1.030      0.303      -7.158       2.228
sepal width (cm)     -6.6809      4.480     -1.491      0.136     -15.461       2.099
petal length (cm)     9.4294      4.737      1.990      0.047       0.145      18.714
petal width (cm)     18.2861      9.743      1.877      0.061      -0.809      37.381
=====================================================================================

Possibly complete quasi-separation: A fraction 0.73 of observations can be
perfectly predicted. This might indicate that there is complete
quasi-separation. In this case some parameters will not be identified.
```

## 参考
- [scikit-learn、Keras、TensorFlowによる実践機械学習 第2版](https://www.oreilly.co.jp/books/9784873119281/)
- [statsmodels 0.14.0](https://www.statsmodels.org/stable/index.html)