統計学の代表的な手法を実践する (2) 〜 ロジスティック回帰

[昨日](http://qiita.com/ynakayama/items/bd0b65716279ae2df957)の続きで統計学の実践的な手法を追っていきます。今日はロジスティック回帰についてです。

## 参考文献

書籍をお持ちの方は[最強である](http://www.amazon.co.jp/dp/4478022216) p189 [実践編](http://www.amazon.co.jp/dp/4478028230) p221 あたりを参照してください。 (後者のほうが詳しいです)

またロジスティック回帰の説明としては

概要については
ロジスティック回帰の考え方・使い方 - TokyoR #33
http://www.slideshare.net/horihorio/tokyo-r33-logi

説明については
ロジスティック回帰と変数選択 (三重大学・奥村先生)
http://oku.edu.mie-u.ac.jp/~okumura/stat/140921.html

また実務における例としては
回帰分析とその応用③　～ロジスティック回帰分析
http://www.gixo.jp/blog/2492

マーケターのためのデータマイニング・ヒッチハイクガイド
第10回：ロジスティック回帰分析
http://jpn.teradata.jp/library/ma/ins_1310a.html

あたりが詳しくわかりやすいと思います。

## ロジスティック回帰の使いどころ

ロジスティック回帰は true か false という二値論理に関するアウトカムを分析するための回帰分析であると言われています。 ([前回の表](http://qiita.com/ynakayama/items/bd0b65716279ae2df957)を参照)

統計の分野では一般的な方法ですが、機械学習の分野でも [PRML 4 章線形識別モデル 4.3.2 ロジスティック回帰 (p204)](http://www.amazon.co.jp/dp/4621061224)として登場します。こちらでは名前は回帰ですが、実際は分類のためのモデルであると強調されています。クラス C1 の事後確率が特徴ベクトル Φ の線形関数のロジスティックシグモイド関数として書けることを利用したモデルです。

以前に挙げた[線形回帰](http://qiita.com/ynakayama/items/d6e8513791b43f18d28b)と比較すると、ごく単純な線形回帰である単回帰分析では直線に回帰するのに対し、あてはめのイメージが対数オッズにもとづいた曲線となるところに特徴があります。

業務上扱いたいアウトカムは「この人は一度でも来店したことがあるかどうか」「この会員は解約するのかどうか」というふうに 0 か 1 かで表現できる論理だったりすることが多いです。

たとえば[前回](http://qiita.com/ynakayama/items/bd0b65716279ae2df957)も登場したケースコントロール調査のような例で、似たような条件なのに病気になった人とならなかった人の間でどのような違いがあるのか、このようなことをより深く調べていこうといったときに、重回帰分析以上のあてはまりを見せることがあります。

## 実践してみる

単純な例では奥村先生のページの[ロジスティック回帰](http://oku.edu.mie-u.ac.jp/~okumura/stat/logistic.html)の例がわかりやすいです。

y = 入試の合否 x1 = 内申点 x2 = 模試の偏差値としたときに

|入試の合否(y)|内申点(x1)|偏差値(x2)|
|:--|:--|:--|
|0|3.6|60.1|
|1|4.1|52.0|
|0|3.7|62.5|
|0|4.9|60.6|
|1|4.4|54.1|
|0|3.6|63.6|
|1|4.9|68.0|
|0|4.8|38.2|
|1|4.1|59.5|
|0|4.3|47.3|

これをロジスティック回帰すると

```r
# データを読み込む
data1 <- read.csv("data1.csv", header=TRUE)
# ロジスティック回帰
mylogit = glm(y ~ x1 + x2, data=data1, 
family=binomial(link="logit"))
# サマリーを見る
summary(mylogit)
#=>
# Call:
# glm(formula = y ~ x1 + x2, family = binomial(link ="logit"), data = data1)
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -1.4754  -0.8584  -0.8007   1.1905   1.5719  

# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept) -9.44589    9.12237  -1.035    0.300
# x1           1.27158    1.49423   0.851    0.395
# x2           0.06424    0.08739   0.735    0.462

# (Dispersion parameter for binomial family taken to be 1)

#     Null deviance: 13.460  on 9  degrees of freedom
# Residual deviance: 12.345  on 7  degrees of freedom
# AIC: 18.345

# Number of Fisher Scoring iterations: 4

# フィットされた値を見る
fitted(mylogit)
#=>
#         1         2         3         4         5         6         7         8         9        10 
# 0.2675020 0.2907201 0.3260736 0.6632579 0.4072105 0.3137831 0.7600997 0.2914609 0.3988916 0.2810008 
```

つまりは以下の式で予測される、またフィットされた値もわかるというわけです。

```math
logit(\pi) = -9.44589 + 1.27158x_{1} + 0.06424x_{2}
```

※ ただし有意ではありません

## 可視化する

カリフォルニア大学ロサンゼルス校のページに掲載されている[こちら](http://www.ats.ucla.edu/stat/r/dae/logit.htm)の例に沿って R で可視化した例を紹介します。ここでは[仮のデータ](http://www.ats.ucla.edu/stat/data/binary.csv)に対してロジスティック回帰をする説明が掲載されています。

```r
# データ読み込み
mydata <- read.csv("binary.csv")
mydata$rank <- factor(mydata$rank)
# ロジスティック回帰
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")
# サマリーを見る
summary(mylogit)
#=>
# Call:
# glm(formula = admit ~ gre + gpa + rank, family = "binomial", data = mydata)

# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -1.6268  -0.8662  -0.6388   1.1490   2.0790  

# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept) -3.989979   1.139951  -3.500 0.000465 ***
# gre          0.002264   0.001094   2.070 0.038465 *  
# gpa          0.804038   0.331819   2.423 0.015388 *  
# rank2       -0.675443   0.316490  -2.134 0.032829 *  
# rank3       -1.340204   0.345306  -3.881 0.000104 ***
# rank4       -1.551464   0.417832  -3.713 0.000205 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# (Dispersion parameter for binomial family taken to be 1)

#     Null deviance: 499.98  on 399  degrees of freedom
# Residual deviance: 458.52  on 394  degrees of freedom
# AIC: 470.52

# Number of Fisher Scoring iterations: 4

newdata1 <- with(mydata, data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))
newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response")
newdata2 <- with(mydata, data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100),
    4), gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))
newdata3 <- cbind(newdata2, predict(mylogit, newdata = newdata2, type = "link",
    se = TRUE))
newdata3 <- within(newdata3, {
    PredictedProb <- plogis(fit)
    LL <- plogis(fit - (1.96 * se.fit))
    UL <- plogis(fit + (1.96 * se.fit))
})
# グラフィックライブラリを使う
library("ggplot2")
ggplot(newdata3, aes(x = gre, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
    ymax = UL, fill = rank), alpha = 0.2) + geom_line(aes(colour = rank),
    size = 1)
# 結果を画像ファイルに落とす
ggsave("image.png")
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/38501/8f2e7da7-f78f-e17f-3d25-ebeb0a0e22c2.png)

## まとめ

ひとまず今回はロジスティック回帰の概要と R で分析する例を紹介しました。ロジスティック回帰の話は次回に続きます。