バイオ・データ・マイニング/Rでロジスティック回帰を使う
をテンプレートにして作成
開始行:
*はじめに [#g88a8135]
ロジスティック回帰は、教師信号 [math]y[/math] が 1 または...
『Rによるバイオインフォマティクスデータ解析』には載ってい...
*準備 [#e5c35924]
Rのインストールについては，次のページを見てください．
-[[MacでRを使う>機械学習/MacでRを使う]]
-[[WindowsでRを使う>機械学習/WindowsでRを使う]]
ここでは、標準で使用できる''irisデータセット''を使います。
#geshi(rsplus){{
data(iris)
}}
このデータセットは、アヤメの種類（Species）をがく片の長さ...
長さと幅は連続値，種類はsetosa, versicolor, virginicaのい...
このデータセットには、setosa, versicolor, virginicaという...
ランダムに10個のデータを選択して、見てみましょう。
#geshi(rsplus){{
iris[sort(sample(1:150,10)),]
}}
#geshi(rsplus){{
Sepal.Length Sepal.Width Petal.Length Petal.Width ...
4 4.6 3.1 1.5 0.2 ...
22 5.1 3.7 1.5 0.4 ...
65 5.6 2.9 3.6 1.3 ver...
97 5.7 2.9 4.2 1.3 ver...
100 5.7 2.8 4.1 1.3 ver...
108 7.3 2.9 6.3 1.8 vi...
116 6.4 3.2 5.3 2.3 vi...
122 5.6 2.8 4.9 2.0 vi...
136 7.7 3.0 6.1 2.3 vi...
146 6.7 3.0 5.2 2.3 vi...
}}
2値分類にするため、Speciesがsetosaなら 1、そうでないなら ...
これを、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 ...
21 5.4 3.4 1.7 0.2 ...
43 4.4 3.2 1.3 0.2 ...
50 5.0 3.3 1.4 0.2 ...
54 5.5 2.3 4.0 1.3 ver...
73 6.3 2.5 4.9 1.5 ver...
76 6.6 3.0 4.4 1.4 ver...
84 6.0 2.7 5.1 1.6 ver...
114 5.7 2.5 5.0 2.0 vi...
139 6.0 3.0 4.8 1.8 vi...
150 5.9 3.0 5.1 1.8 vi...
}}
これをランダム・サンプリングして80%を訓練データとし、のこ...
#geshi(rsplus){{
idx <- sample(nrow(iris), as.integer(nrow(iris)*0.8))
iris.train <- iris[idx,]
iris.test <- iris[-idx,]
}}
*ロジスティック回帰 [#u50975a4]
説明のため、まずは説明変数をPetal.Length、目的変数をVirgi...
#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) と...
このロジット関数を線形回帰で求めるのが''ロジスティック回...
*ロジスティック回帰を実行する [#z45c2dde]
ロジスティック回帰を行うには、''glm''関数を使います。
glm関数には、第１引数にモデル式、第２引数にデータフレーム...
#geshi(rsplus){{
model <- glm(Virginica ~ Petal.Length, data=iris.train, f...
}}
学習したモデルのサマリーを出力してみます。
#geshi(rsplus){{
summary(model)
}}
#geshi(rsplus){{
Call:
glm(formula = Virginica ~ Petal.Length, family = binomial...
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 ...
(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がロジスティック回帰で求められたロ...
したがって、目的変数の値が 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 ...
1.129131e-12 4.219567e-02 1.825599e-06 1.522187e-03 4.219...
118 11 67 19 ...
9.999998e-01 4.870000e-13 4.219567e-02 2.617940e-12 3.544...
}}
インデックスをランダムに選んだままでソートしていないので...
テストデータに対して予測をするときは、''predict.glm''関数...
predict.glm関数の第１引数には学習したモデル、第２引数には...
#geshi(rsplus){{
predict.glm(model, iris.test, type="response")
}}
#geshi(rsplus){{
5 7 18 22 ...
2.100456e-13 2.100456e-13 2.100456e-13 4.870000e-13 2.617...
47 55 59 62 ...
1.129131e-12 9.267631e-02 9.267631e-02 3.522188e-03 2.835...
108 115 124 126 ...
9.999939e-01 8.725072e-01 5.600666e-01 9.999245e-01 9.999...
}}
*モデルの評価 [#c6e49a5d]
ロジスティック回帰で求めたモデルを評価するには、''ROC曲線...
ROC曲線とAUCを求めるには''pROC''パッケージを使います。
#geshi(rsplus){{
install.packages("pROC")
library(pROC)
}}
テストデータの目的変数とテストデータに対する予測値を変数...
#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...
Area under the curve: 0.9928
}}
ROC曲線をプロットするには、roc関数で生成されるオブジェク...
#geshi(rsplus){{
roc <- roc(y.test, p.test)
plot(roc)
}}
#ref(roc.png,nolink)
*まとめ [#z31308ae]
ロジスティック回帰は、目的変数が 1 になる確率のロジット関...
*参考文献 [#g8f00b3a]
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0"...
<iframe style="width:120px;height:240px;" marginwidth="0"...
}}
終了行:
*はじめに [#g88a8135]
ロジスティック回帰は、教師信号 [math]y[/math] が 1 または...
『Rによるバイオインフォマティクスデータ解析』には載ってい...
*準備 [#e5c35924]
Rのインストールについては，次のページを見てください．
-[[MacでRを使う>機械学習/MacでRを使う]]
-[[WindowsでRを使う>機械学習/WindowsでRを使う]]
ここでは、標準で使用できる''irisデータセット''を使います。
#geshi(rsplus){{
data(iris)
}}
このデータセットは、アヤメの種類（Species）をがく片の長さ...
長さと幅は連続値，種類はsetosa, versicolor, virginicaのい...
このデータセットには、setosa, versicolor, virginicaという...
ランダムに10個のデータを選択して、見てみましょう。
#geshi(rsplus){{
iris[sort(sample(1:150,10)),]
}}
#geshi(rsplus){{
Sepal.Length Sepal.Width Petal.Length Petal.Width ...
4 4.6 3.1 1.5 0.2 ...
22 5.1 3.7 1.5 0.4 ...
65 5.6 2.9 3.6 1.3 ver...
97 5.7 2.9 4.2 1.3 ver...
100 5.7 2.8 4.1 1.3 ver...
108 7.3 2.9 6.3 1.8 vi...
116 6.4 3.2 5.3 2.3 vi...
122 5.6 2.8 4.9 2.0 vi...
136 7.7 3.0 6.1 2.3 vi...
146 6.7 3.0 5.2 2.3 vi...
}}
2値分類にするため、Speciesがsetosaなら 1、そうでないなら ...
これを、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 ...
21 5.4 3.4 1.7 0.2 ...
43 4.4 3.2 1.3 0.2 ...
50 5.0 3.3 1.4 0.2 ...
54 5.5 2.3 4.0 1.3 ver...
73 6.3 2.5 4.9 1.5 ver...
76 6.6 3.0 4.4 1.4 ver...
84 6.0 2.7 5.1 1.6 ver...
114 5.7 2.5 5.0 2.0 vi...
139 6.0 3.0 4.8 1.8 vi...
150 5.9 3.0 5.1 1.8 vi...
}}
これをランダム・サンプリングして80%を訓練データとし、のこ...
#geshi(rsplus){{
idx <- sample(nrow(iris), as.integer(nrow(iris)*0.8))
iris.train <- iris[idx,]
iris.test <- iris[-idx,]
}}
*ロジスティック回帰 [#u50975a4]
説明のため、まずは説明変数をPetal.Length、目的変数をVirgi...
#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) と...
このロジット関数を線形回帰で求めるのが''ロジスティック回...
*ロジスティック回帰を実行する [#z45c2dde]
ロジスティック回帰を行うには、''glm''関数を使います。
glm関数には、第１引数にモデル式、第２引数にデータフレーム...
#geshi(rsplus){{
model <- glm(Virginica ~ Petal.Length, data=iris.train, f...
}}
学習したモデルのサマリーを出力してみます。
#geshi(rsplus){{
summary(model)
}}
#geshi(rsplus){{
Call:
glm(formula = Virginica ~ Petal.Length, family = binomial...
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 ...
(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がロジスティック回帰で求められたロ...
したがって、目的変数の値が 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 ...
1.129131e-12 4.219567e-02 1.825599e-06 1.522187e-03 4.219...
118 11 67 19 ...
9.999998e-01 4.870000e-13 4.219567e-02 2.617940e-12 3.544...
}}
インデックスをランダムに選んだままでソートしていないので...
テストデータに対して予測をするときは、''predict.glm''関数...
predict.glm関数の第１引数には学習したモデル、第２引数には...
#geshi(rsplus){{
predict.glm(model, iris.test, type="response")
}}
#geshi(rsplus){{
5 7 18 22 ...
2.100456e-13 2.100456e-13 2.100456e-13 4.870000e-13 2.617...
47 55 59 62 ...
1.129131e-12 9.267631e-02 9.267631e-02 3.522188e-03 2.835...
108 115 124 126 ...
9.999939e-01 8.725072e-01 5.600666e-01 9.999245e-01 9.999...
}}
*モデルの評価 [#c6e49a5d]
ロジスティック回帰で求めたモデルを評価するには、''ROC曲線...
ROC曲線とAUCを求めるには''pROC''パッケージを使います。
#geshi(rsplus){{
install.packages("pROC")
library(pROC)
}}
テストデータの目的変数とテストデータに対する予測値を変数...
#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...
Area under the curve: 0.9928
}}
ROC曲線をプロットするには、roc関数で生成されるオブジェク...
#geshi(rsplus){{
roc <- roc(y.test, p.test)
plot(roc)
}}
#ref(roc.png,nolink)
*まとめ [#z31308ae]
ロジスティック回帰は、目的変数が 1 になる確率のロジット関...
*参考文献 [#g8f00b3a]
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0"...
<iframe style="width:120px;height:240px;" marginwidth="0"...
}}
ページ名:
![[PukiWiki]](https://xn--p8ja5bwe1i.jp:443/wiki/image/banner.png "[PukiWiki]")