ロジスティック回帰（備忘録）

ロジスティック回帰分析ついての備忘録。（個人メモ）

# ロジスティック回帰とロジットリンク関数
確率分布に二項分布(Binomial Distribution)、リンク関数にロジットリンク関数を指定した一般化線形モデル（GLM）が、**ロジスティック回帰**。  

二項分布では、事象の発生確率 $p_i$ を指定する必要があり、$p_i$ の値の範囲 0≥$p_i$≥1 （確率なので0～100%の範囲）を表現するために適したリンク関数が、ロジットリンク関数。

## ロジスティック関数
ロジットリンク関数は、ロジスティック関数の逆関数。ロジスティック関数は、以下の式で表現される。 
$z_i$ は線形予測子 $z_i=α_1+β_1x_1+β_2x_2…$  　　・・・α：切片　$β_i$：係数

$$ p_i=logistic(x_i)=\frac{1}{1+\exp(-z_i)} $$

### ロジスティック関数を描画（シグモイドカーブ）
上記の式の $z$ を ***-10〜10*** で変化させた時の曲線を描画してみる。
$z$ がどのような値になっても 0≥𝑝𝑖≥1 に収まっていて都合が良いこと、またz=0でp=0.5となることに着目！

```
# -6～6の間で等間隔に1000個の値を取得
x = np.linspace(-10,10,1000)
# 1/(1+np.exp(-x))　がロジスティック関数
plt.plot(x, ( 1/(1+np.exp(-x))) )

plt.axvline(x=0, color='red', linestyle='dotted')
plt.axhline(y=0.5, color='red', linestyle='dotted')
plt.ylabel("確率(p)")
plt.xlabel("確率変数(z)")
plt.show()
```

下記のようなグラフが描画される。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1235367/1d86a61f-9ed1-ca9e-0541-542e52832050.png)


## ロジットリンク関数
ロジスティック関数を $z_i$= の形に変形すると以下のようになる。
$$ \log \frac{p_i}{1-p_i}=z_i $$
※ちなみに、$z_i=α_1+β_1x_1+β_2x_2…$  　　・・・α：切片　$β_i$：係数

この左辺をロジット関数という。これってつまり以下のようなこと。

$$ α_1+β_1x_1+β_2x_2… = \log \frac{p_i}{1-p_i} $$


# ロジスティック回帰によるパラメータの推定
ロジスティック回帰によるパラメータの推定。
（線形結合の式（線形予測子）の**切片**や**各係数***の推定）

尤度関数と対数尤度関数は以下で定義され、$\log L$ を最大にする推定値 $\hat{\beta_j}$ を探し出します。
$$ L(\beta_j)=\prod_i{}_n\mathrm{C}_rp_i^{y_i}(1-p_i)^{N_i-y_i} $$

両辺の対数をとると

$$ \log L({\beta_j})=\sum_i (\log {}_n \mathrm{C} _r+y_i\log(p_i)+(N_i-y_i)\log(1-p_i)) $$


# 実装
適当にデータを用意。

```
# 説明変数となるデータを1項目用意
X_train_list = np.r_[np.random.normal(3 ,1, size=50), np.random.normal(-1,1,size=50)]
# 目的変数となるデータを用意
y_train_list = np.r_[np.ones(50),np.zeros(50)]

# データフレームに変換
df = pd.DataFrame(data=[X_train_list,y_train_list]).T
df.rename(columns={0:"X",1:"Y"}, inplace=True)

# 目的変数と説明変数を分ける
X_train = df[["X"]]
y_train = df[["Y"]]

display(X_train.head())
display(y_train.head())
```

```
X
0	5.207395
1	2.692299
2	1.560036
3	3.853632
4	4.804304

Y
0	1.0
1	1.0
2	1.0
3	1.0
4	1.0
```

## sklearn.linear_model.LogisticRegressionの場合
### モデル作成
```
from sklearn.linear_model import LogisticRegression
# モデルのインスタンス化
model = LogisticRegression()
# 学習
model.fit(X_train,y_train)
```

### 予測（フラグ）

```
モデル名.predict(説明変数データ)
```

0と1のリストが返ってくる。デフォルトでは、閾値の0.5 (=50%)を超える場合には1、そうでない場合には0で設定されているとのこと。
引数で渡す説明変数データは***DataFrame***じゃないとエラーになるみたい。***Seriese***だとエラーになってしまった。

```
model.predict(X_train)
```

```
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
```

### 予測（確率）
```
モデル名.predict_proba(説明変数データ)
```
array形式で、0に分類される確率、1に分類される確率の各確率が返ってくる。

```
model.predict_proba(X_train.iloc[0:10])
```

```
array([[2.63082246e-05, 9.99973692e-01],
       [1.51340331e-02, 9.84865967e-01],
       [2.12846522e-01, 7.87153478e-01],
       [8.10623998e-04, 9.99189376e-01],
       [7.30213967e-05, 9.99926979e-01],
       [1.51565245e-02, 9.84843475e-01],
       [6.95446239e-04, 9.99304554e-01],
       [7.48031608e-02, 9.25196839e-01],
       [1.88389023e-02, 9.81161098e-01],
       [5.45899840e-04, 9.99454100e-01]])
```

ちなみに、1に分類される確率だけを出力する場合は、

```
モデル名.predict_proba(評価用の説明変数データ)[:, 1]
```

## statsmodels.apiの場合
### モデル作成
ポイントは、***sm.add_constant(説明変数データ)***をしないと、回帰式の切片が出てこないところ。

```
import statsmodels.api as sm

# 定数項（切片）のために行う
X_train = sm.add_constant(X_train)

glm_binom = sm.GLM( y_train, X_train, family=sm.families.Binomial())
res = glm_binom.fit()
print(res.summary())
```

```
                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      Y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       98
Model Family:                Binomial   Df Model:                            1
Link Function:                  logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -8.8540
Date:                Thu, 25 Mar 2021   Deviance:                       17.708
Time:                        10:54:26   Pearson chi2:                     48.1
No. Iterations:                     8                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -3.1340      1.046     -2.995      0.003      -5.185      -1.083
X              3.3071      0.952      3.476      0.001       1.442       5.172
==============================================================================
```

### 予測
1に分類される確率が返ってくる。
***predict***する前に、引数で渡す説明変数は、***sm.add_constant(説明変数データ)***をしないとエラーになります。

```
モデル.predict(説明変数データ)
```

```
res.predict(X_train)
```

```
0     0.527843
1     0.999978
2     0.999281
3     0.997063
4     0.995389
        ...   
95    0.001322
96    0.000229
97    0.004302
98    0.000940
99    0.000040
```

### AIC（赤池情報量基準）の出力

```
res.aic
```


## statsmodels.formula.apiの場合
statsmodels.apiよりも細かい指定ができるというものなのだろうか・・・

### モデル作成
formulaで　目的変数と説明変数を列挙
Rのように「全体から、この列とこの列を抜く！」という指定はできないのだろうか・・・

```
import statsmodels.formula.api as smf
# 線形予測子を定義　「1+」は切片
formula = "Y~1+X"

# リンク関数を選択
link = sm.genmod.families.links.logit

# 誤差構造(損失関数)を選択
family = sm.families.Binomial(link=link)

mod = smf.glm(formula=formula, data=df, family=family )
result = mod.fit() 
print(result.summary())
```

```
                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      Y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       98
Model Family:                Binomial   Df Model:                            1
Link Function:                  logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -8.8540
Date:                Thu, 25 Mar 2021   Deviance:                       17.708
Time:                        11:11:39   Pearson chi2:                     48.1
No. Iterations:                     8                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -3.1340      1.046     -2.995      0.003      -5.185      -1.083
X              3.3071      0.952      3.476      0.001       1.442       5.172
==============================================================================
```

### 予測
1に分類される確率が返ってくる。
***predict***する前に、引数で渡す説明変数は、***sm.add_constant(説明変数データ)***をしないとエラーになります。

```
モデル.predict(説明変数データ)
```

```
result.predict(X_train)
```

```
0     0.527843
1     0.999978
2     0.999281
3     0.997063
4     0.995389
        ...   
95    0.001322
96    0.000229
97    0.004302
98    0.000940
99    0.000040
```


### AIC（赤池情報量基準）の出力

```
result.aic
```

##解釈
ロジット関数の線形予測子を、得られた係数の値で置き換えると、以下の式になります。
例えば、切片が -19.54 $β_1$が1.95　$β_2$が2.02 だったとすると・・・

$$ \frac{p_i}{1-p_i}=\exp(\alpha+\beta_1x_i+\beta_2f_i)=\exp(\alpha)+\exp(\beta_1x_i)+\exp(\beta_2f_i)=\exp(-19.54)+\exp(1.95x_i)+\exp(2.02f_i) $$

左辺は (事象が発生する確率)/(事象が発生しない確率) のオッズになっており、例えば$x$が1増加するとオッズ比は $exp(1.95)=7.06$ 倍になる。あくまで、オッズが7.06倍になる点に注意！

## おまけ
予測結果を可視化
シグモイドカーブ描いてるね！

```
plt.scatter(df["X"],df["Y"])
plt.plot(df["X"],result.predict(df["X"]),"-o",linestyle='None',color="red")
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1235367/83e08fb0-15d3-1e94-c634-a30ee90c0bd4.png)



モデルの評価については、別で。

おしまい