*はじめに [#g88a8135]
ロジスティック回帰は、教師信号 [math]y[/math] が 1 または 0 の2値分類問題（カテゴリーが「はい」と「いいえ」のカテゴリー分析）において、[math]y = 1[/math] となる確率を表す関数を回帰分析によって求める手法です。
『Rによるバイオインフォマティクスデータ解析』には載っていません。
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0" marginheight="0" scrolling="no" frameborder="0" src="https://rcm-fe.amazon-adsystem.com/e/cm?ref=tf_til&t=tohgorohmatsu-22&m=amazon&o=9&p=8&l=as1&IS2=1&detail=1&asins=4320057082&linkId=f6db5311dcfd9a82f5dee20859d39574&bc1=ffffff<1=_blank&fc1=444b4c&lc1=444b4c&bg1=ffffff&f=ifr"></iframe>
}}
*準備 [#e5c35924]
Rのインストールについては，次のページを見てください．
-[[MacでRを使う>機械学習/MacでRを使う]]
-[[WindowsでRを使う>機械学習/WindowsでRを使う]]
ここでは、標準で使用できる''irisデータセット''を使います。
#geshi(rsplus){{
data(iris)
}}
このデータセットは、アヤメの種類（Species）をがく片の長さ（Sepal.Length）、幅（Lepal.Width）、花びらの長さ（Petal.Length）、幅（Petal.Width）によって分類する問題です。
長さと幅は連続値，種類はsetosa, versicolor, virginicaのいずれかをとる離散値です。
このデータセットには、setosa, versicolor, virginicaという3種類のアヤメについて、それぞれ50個ずつ、合計150個のデータが含まれています。
ランダムに10個のデータを選択して、見てみましょう。
#geshi(rsplus){{
iris[sort(sample(1:150,10)),]
}}
#geshi(rsplus){{
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
4 4.6 3.1 1.5 0.2 setosa
22 5.1 3.7 1.5 0.4 setosa
65 5.6 2.9 3.6 1.3 versicolor
97 5.7 2.9 4.2 1.3 versicolor
100 5.7 2.8 4.1 1.3 versicolor
108 7.3 2.9 6.3 1.8 virginica
116 6.4 3.2 5.3 2.3 virginica
122 5.6 2.8 4.9 2.0 virginica
136 7.7 3.0 6.1 2.3 virginica
146 6.7 3.0 5.2 2.3 virginica
}}
2値分類にするため、Speciesがsetosaなら 1、そうでないなら 0 という新しい列 Setosa を作ります。
これを、versicolor, virginica についても行います。
#geshi(rsplus){{
iris$Setosa <- ifelse(iris$Species=="setosa", 1, 0)
iris$Versicolor <- ifelse(iris$Species=="versicolor", 1, 0)
iris$Virginica <- ifelse(iris$Species=="virginica", 1, 0)
}}
もう一度、ランダムに10個のデータを選択して、追加した列を確認します。
#geshi(rsplus){{
Sepal.Length Sepal.Width Petal.Length Petal.Width Species Setosa Versicolor Virginica
21 5.4 3.4 1.7 0.2 setosa 1 0 0
43 4.4 3.2 1.3 0.2 setosa 1 0 0
50 5.0 3.3 1.4 0.2 setosa 1 0 0
54 5.5 2.3 4.0 1.3 versicolor 0 1 0
73 6.3 2.5 4.9 1.5 versicolor 0 1 0
76 6.6 3.0 4.4 1.4 versicolor 0 1 0
84 6.0 2.7 5.1 1.6 versicolor 0 1 0
114 5.7 2.5 5.0 2.0 virginica 0 0 1
139 6.0 3.0 4.8 1.8 virginica 0 0 1
150 5.9 3.0 5.1 1.8 virginica 0 0 1
}}
これをランダム・サンプリングして80%を訓練データとし、のこりの20％をテスト・データとします。
#geshi(rsplus){{
idx <- sample(nrow(iris), as.integer(nrow(iris)*0.8))
iris.train <- iris[idx,]
iris.test <- iris[-idx,]
}}
*ロジスティック回帰 [#u50975a4]
説明のため、まずは説明変数をPetal.Length、目的変数をVirginicaとします。
#geshi(rsplus){{
x.train <- iris.train$Petal.Length
y.train <- iris.train$Virginica
plot(x.train, y.train)
}}
#ref(iris_petal_length_virginica.png,nolink)
このデータに対して線形回帰でモデルを作成すると、うまくフィットしません。
#geshi(rsplus){{
model <- lm(y.train ~ x.train)
abline(model)
}}
#ref(iris_petal_length_virginica_lm.png,nolink)
そこで、目的変数の値が 1 になる確率を [math]p[/math] とし、この確率を
$$ p = \frac{1}{1+\exp(-(\beta_0+\beta_1 x))} $$
という関数で近似することを考えます。
この関数は、次のような 0 から 1 までの値をとり、S字の形をしています。
#ref(sigmoid.png,nolink)
この式を変形すると、次のようになります。
$$ \log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 x$$
左辺の [math]\frac{p}{1-p}[/math] を''オッズ'' (odds) といい、[math]\log\left(\frac{p}{1-p}\right)[/math] を''ロジット関数''といいます。
このロジット関数を線形回帰で求めるのが''ロジスティック回帰''です。
*ロジスティック回帰を実行する [#z45c2dde]
ロジスティック回帰を行うには、''glm''関数を使います。
glm関数には、第１引数にモデル式、第２引数にデータフレーム、第3引数に ''family=binomial(link="logit")'' を指定します。
#geshi(rsplus){{
model <- glm(Virginica ~ Petal.Length, data=iris.train, family=binomial(link="logit"))
}}
学習したモデルのサマリーを出力してみます。
#geshi(rsplus){{
summary(model)
}}
#geshi(rsplus){{
Call:
glm(formula = Virginica ~ Petal.Length, family = binomial(link = "logit"),
data = iris.train)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.02963 -0.03626 0.00000 0.02848 2.51612
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -40.965 11.058 -3.705 0.000212 ***
Petal.Length 8.409 2.267 3.709 0.000208 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 154.11 on 119 degrees of freedom
Residual deviance: 31.59 on 118 degrees of freedom
AIC: 35.59
Number of Fisher Scoring iterations: 9
}}
CoefficientsのEstimateがロジスティック回帰で求められたロジット関数の係数を表しており、Interceptが定数項 [math]\beta_0[/math] です。
したがって、目的変数の値が 1 になる確率は次のように表されます。
$$ p = \frac{1}{1 + \exp(-(-40.965 + 8.409 x))} $$
これをプロットしてみます。
#geshi(rsplus){{
plot(x.train, y.train)
curve(1/(1+exp(-(-40.965+8.409*x))), add=TRUE)
}}
#ref(iris_petal_length_virginica_logistic.png,nolink)
訓練データに対する確率は、''fitted''関数で出力できます。
#geshi(rsplus){{
fitted(model)
}}
#geshi(rsplus){{
26 85 58 100 86 105 48 66 135 129 90
1.129131e-12 4.219567e-02 1.825599e-06 1.522187e-03 4.219567e-02 9.995945e-01 2.100456e-13 1.864666e-02 9.978238e-01 9.978238e-01 6.570957e-04
118 11 67 19 127 2 84 102 130 68 41
9.999998e-01 4.870000e-13 4.219567e-02 2.617940e-12 3.544564e-01 2.100456e-13 8.725072e-01 8.725072e-01 9.995945e-01 1.522187e-03 2.220446e-16...
}}
インデックスをランダムに選んだままでソートしていないので、順番はバラバラになっています。
テストデータに対して予測をするときは、''predict.glm''関数を用います。
predict.glm関数の第１引数には学習したモデル、第２引数にはテストデータのデータフレーム、第３引数には''type="response"''を指定します。
#geshi(rsplus){{
predict.glm(model, iris.test, type="response")
}}
#geshi(rsplus){{
5 7 18 22 24 30 33 36 38 39 43
2.100456e-13 2.100456e-13 2.100456e-13 4.870000e-13 2.617940e-12 1.129131e-12 4.870000e-13 2.220446e-16 2.100456e-13 2.220446e-16 2.220446e-16
47 55 59 62 70 72 76 91 93 94 103
1.129131e-12 9.267631e-02 9.267631e-02 3.522188e-03 2.835147e-04 6.570957e-04 1.864666e-02 1.864666e-02 6.570957e-04 1.825599e-06 9.998250e-01
108 115 124 126 131 141 145 149
9.999939e-01 8.725072e-01 5.600666e-01 9.999245e-01 9.999675e-01 9.978238e-01 9.990602e-01 9.884120e-01
}}
*モデルの評価 [#c6e49a5d]
ロジスティック回帰で求めたモデルを評価するには、''ROC曲線''の下側の面積 (''AUC'': Area Under the Curve) を使います。
ROC曲線とAUCを求めるには''pROC''パッケージを使います。
#geshi(rsplus){{
install.packages("pROC")
library(pROC)
}}
テストデータの目的変数とテストデータに対する予測値を変数に代入しておき、''roc''関数でROC曲線とAUCを求めます。
#geshi(rsplus){{
y.test <- iris.test$Virginica
p.test <- predict.glm(model, iris.test, type="response")
roc(y.test, p.test)
}}
#geshi(rsplus){{
Call:
roc.default(response = y.test, predictor = p.test)
Data: p.test in 19 controls (y.test 0) < 11 cases (y.test 1).
Area under the curve: 0.9928
}}
ROC曲線をプロットするには、roc関数で生成されるオブジェクトをplot関数に渡します。
#geshi(rsplus){{
roc <- roc(y.test, p.test)
plot(roc)
}}
#ref(roc.png,nolink)
*まとめ [#z31308ae]
ロジスティック回帰は、目的変数が 1 になる確率のロジット関数を線形回帰によって求める手法で、2値分類（カテゴリーが「はい」または「いいえ」のカテゴリー分析）に用いる手法です。
*参考文献 [#g8f00b3a]
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0" marginheight="0" scrolling="no" frameborder="0" src="https://rcm-fe.amazon-adsystem.com/e/cm?ref=tf_til&t=tohgorohmatsu-22&m=amazon&o=9&p=8&l=as1&IS2=1&detail=1&asins=4320057082&linkId=f6db5311dcfd9a82f5dee20859d39574&bc1=ffffff<1=_blank&fc1=444b4c&lc1=444b4c&bg1=ffffff&f=ifr"></iframe>
<iframe style="width:120px;height:240px;" marginwidth="0" marginheight="0" scrolling="no" frameborder="0" src="https://rcm-fe.amazon-adsystem.com/e/cm?ref=qf_sp_asin_til&t=tohgorohmatsu-22&m=amazon&o=9&p=8&l=as1&IS2=1&detail=1&asins=4621061224&linkId=632f1721f5fa1b024003cfb1b98ce366&bc1=ffffff<1=_blank&fc1=444b4c&lc1=444b4c&bg1=ffffff&f=ifr"></iframe>
<iframe style="width:120px;height:240px;" marginwidth="0" marginheight="0" scrolling="no" frameborder="0" src="https://rcm-fe.amazon-adsystem.com/e/cm?ref=qf_sp_asin_til&t=tohgorohmatsu-22&m=amazon&o=9&p=8&l=as1&IS2=1&detail=1&asins=4621061240&linkId=24ebe891b6a9641e04a983f2259eb3ab&bc1=ffffff<1=_blank&fc1=444b4c&lc1=444b4c&bg1=ffffff&f=ifr"></iframe>
<iframe style="width:120px;height:240px;" marginwidth="0" marginheight="0" scrolling="no" frameborder="0" src="https://rcm-fe.amazon-adsystem.com/e/cm?ref=qf_sp_asin_til&t=tohgorohmatsu-22&m=amazon&o=9&p=8&l=as1&IS2=1&detail=1&asins=4320121341&linkId=f7fea767d9538f35cc900e57d2e4a1b1&bc1=ffffff<1=_blank&fc1=444b4c&lc1=444b4c&bg1=ffffff&f=ifr"></iframe>
}}
![[PukiWiki]](https://xn--p8ja5bwe1i.jp:443/wiki/image/banner.png "[PukiWiki]")
