バイオ・データ・マイニング/Rで回帰分析する
をテンプレートにして作成
開始行:
*はじめに [#gc93c04a]
ここでは、『Rによるバイオインフォマティクスデータ解析』の...
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0"...
}}
きちんと理解するために、まずは単回帰と重回帰という2種類の...
その後、多項回帰と二項回帰をやります。
*準備 [#e708ecca]
Rのインストールについては、次のページを見てください。
-[[MacでRを使う>機械学習/MacでRを使う]]
-[[WindowsでRを使う>機械学習/WindowsでRを使う]]
最初は、標準で使用できる''irisデータセット''を使います。
#geshi(rsplus){{
data(iris)
}}
このデータセットは、アヤメの種類(Species)をがく片の長さ...
長さと幅は連続値、種類はsetosa, versicolor, virginicaのい...
このデータセットには、150個のデータが含まれています。
ランダムに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...
}}
*線形回帰(単回帰) [#s1fe6159]
今から作るモデルによって説明したい変数を''目的変数''とい...
統計の分野では、目的変数を''従属変数''、説明変数を''独立...
まず、説明変数が1つの''単回帰''(single linear regression...
つまり、モデルが次のような式で表される直線をデータから求...
\[\hat{y} = \beta_0 + \beta_1 x_1\]
ここで、[math]\hat{y}[/math]がモデルによる''予測値''、[ma...
ある観測値 [math](x_i, y_i)[/math] に対して、モデルの予測...
\[ y_i - \hat{y_i} = y_i - (\beta_0 + \beta_1 x_{i, 1}) \]
線形回帰では、''残差の二乗の和''を最小化するような係数と...
つまり、単回帰では、[math]N[/math] 個の観測データに対して...
(求め方の説明は省きます。)
\[ \sum_{i=1}^N \left(y_i - (\beta_0 + \beta_1 x_{i,1})\r...
例として、がくの長さ(Petal.Length)を横軸、花びらの長さ...
#geshi(rsplus){{
plot(iris$Sepal.Length, iris$Petal.Length)
}}
#ref(iris_single_rawdata.png,nolink,50%)
どうやら、花びらの長さとがくの長さの間には正の相関(一方...
そこで、花びらの長さ(Sepal.Length)を [math]x_1[/math](...
Rでは、モデルの式を、記号''~''を使って「目的変数 ~ 説明変...
したがって、目的変数 [math]y[/math] をがくの長さ(Petal.L...
#geshi(rsplus){{
y <- iris$Petal.Length
x1 <- iris$Sepal.Length
}}
線形回帰は、標準で用意されている''lm関数''にモデル式を入...
#geshi(rsplus){{
iris.sgl <- lm(y ~ x1)
}}
学習された結果を表示してみると、次のようになります。
#geshi(rsplus){{
iris.sgl
}}
#geshi(rsplus){{
Call:
lm(formula = y ~ x1)
Coefficients:
(Intercept) x1
-7.101 1.858
}}
Coefficientsは係数を表し、Interceptは定数項の値(切片)を...
つまり、花びらの長さを [math]x_1[/math]、がくの長さを [ma...
\[ y = -7.101 + 1.858 x_1 \]
学習された単回帰モデルを図示します。
#geshi(rsplus){{
plot(x1, y)
abline(iris.sgl)
}}
#ref(iris_single_model.png,nolink,50%)
学習されたモデルを''summary関数''に渡すと、モデルに関する...
#geshi(rsplus){{
summary(iris.sgl)
}}
#geshi(rsplus){{
Call:
lm(formula = y ~ x1)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
x1 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ...
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16
}}
''Residuals''は、残差の最小値、中央値、最大値、四分位数を...
''Coefficients''は、''切片'' (Intercept) [math]\beta_0[/...
また、それぞれの推定値が 0 であるという帰無仮説に対するt...
p値が0.05を超えている場合は、その説明変数は有意水準5%で統...
''Multiple R-squared''は、''決定係数''で、このモデルがど...
単回帰では説明変数を1つしか使いませんが、説明変数を複数使...
''Adjusted R-squared''は、''自由度調整済み決定係数''で、...
説明変数の数が異なるモデルを比較するときには、自由度調整...
最後に、係数が全て 0 であるという帰無仮説に対するF検定の'...
p値が0.05を超えている場合は、このモデルは有意水準5%で統計...
次に、予測誤差を視覚化するため、それぞれのデータについて...
予測値は''predict関数''で求めることができます。
#geshi(rsplus){{
plot(y, predict(iris.sgl))
abline(a=0,b=1)
}}
#ref(iris_single_errors.png,nolink,50%)
この予測がどのくらい良いのかを表すために、''平均二乗誤差'...
\[\textit{MSE} = \frac{1}{N}\sum_{i=1}^{N} \left(y_i - (\...
#geshi(rsplus){{
mean((y - predict(iris.sgl))**2)
}}
#geshi(rsplus){{
[1] 0.743061
}}
*線形回帰(重回帰) [#bd685357]
説明変数が複数ある線形回帰を''重回帰''(multiple linear r...
重回帰のモデルは次のような式で表されます。
\[y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta...
ここで、[math]\beta[/math]と[math]\bf{x}[/math]は、それぞ...
また、単回帰のときと同様に、次式で表される二乗誤差の総和...
\[\sum_{i=1}^N \left(y_i - (\beta_0 + \beta^\mathrm{T}\bf...
ここでは、がくの長さ(Petal.Length)を目的変数 [math]y[/m...
#geshi(rsplus){{
y <- iris$Petal.Length
x1 <- iris$Sepal.Length
x2 <- iris$Sepal.Width
x3 <- iris$Petal.Width
}}
説明変数が複数あるときは、モデルを表すときにそれぞれの変...
ここでは、「y ~ x1+x2+x3」となります。
単回帰と同じように、''lm関数''を用いてモデルを学習します。
#geshi(rsplus){{
iris.mult <- lm(y ~ x1+x2+x3)
iris.mult$coefficients
}}
#geshi(rsplus){{
(Intercept) x1 x2 x3
-0.2627112 0.7291384 -0.6460124 1.4467934
}}
単回帰と同じように、学習された式は次のように書けます。
\[ y = -0.2627112 + 0.7291384 x_1 - 0.6460124 x_2 + 1.446...
単回帰と同じように、学習されたモデルを''summary関数''に渡...
#geshi(rsplus){{
summary(iris.mult)
}}
#geshi(rsplus){{
Call:
lm(formula = y ~ x1 + x2 + x3)
Residuals:
Min 1Q Median 3Q Max
-0.99333 -0.17656 -0.01004 0.18558 1.06909
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.26271 0.29741 -0.883 0.379
x1 0.72914 0.05832 12.502 <2e-16 ***
x2 -0.64601 0.06850 -9.431 <2e-16 ***
x3 1.44679 0.06761 21.399 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ...
Residual standard error: 0.319 on 146 degrees of freedom
Multiple R-squared: 0.968, Adjusted R-squared: 0.9674
F-statistic: 1473 on 3 and 146 DF, p-value: < 2.2e-16
}}
重回帰分析の場合はモデルが視覚化しにくいので、予測誤差だ...
#geshi(rsplus){{
plot(y, predict(iris.mult))
abline(a=0,b=1)
}}
#ref(iris_multiple_errors.png,nolink,50%)
単回帰と同様に、平均二乗誤差を求めます。
\[\textit{MSE} = \frac{1}{N} \sum_{i=1}^N (y_i - (\beta_0...
#geshi(rsplus){{
mean((y - predict(iris.mult))**2)
}}
#geshi(rsplus){{
[1] 0.09901965
}}
単回帰の場合よりも誤差がかなり小さくなっていることがわか...
*Lasso回帰 [#p3d71085]
''Lasso回帰''(Lasso regression)は、モデル式は重回帰と同...
\[\sum_{i=1}^N \left( y_i - (\beta_0 + \beta^\mathrm{T}\b...
ここで、[math]\lambda|\beta|[/math]を''正則化項''といい、...
Lasso回帰の正則化に用いられている[math]||\beta||_1[/math]...
\[||\beta||_1 = \sum_{i=1}^p|\beta_i|\]
つまり、L1ノルムは係数の絶対値の和を表します。
正則化パラメーター[math]\lambda[/math]の値が大きいほど0に...
Lasso回帰は''glmnetパッケージ''((Friedman, Hastie, and Ti...
そこで、glmnetパッケージをインストールして、使用します。
#geshi(rsplus){{
install.packages("glmnet")
library(glmnet)
}}
Lasso回帰を行うには、glmnetパッケージの''glmnet関数''を用...
引数には、説明変数行列と目的変数(そしてオプション)を与...
#geshi(rsplus){{
x <- cbind(x1, x2, x3)
iris.lasso <- glmnet(x, y)
}}
cbind関数は、列 (column) で結合するもので、x1, x2, x3 の...
plot関数を使うと、正則化項の大きさとそれぞれの係数の関係...
#geshi(rsplus){{
plot(iris.lasso)
}}
#ref(iris_lasso_regularization.png,nolink,50%)
黒色の折れ線が花びらの長さ(Sepal.Length)、赤色の折れ線...
縦軸はそれぞれの係数を表しています。
横軸の左端が最も正則化が強いとき、つまり、係数へのペナル...
最も正則化が最も強いときは、すべての係数が0になっています。
そこから正則化を弱めていくと、緑色の線が表すがくの幅(Pet...
その後、他の係数も0ではなくなり、正則化項が0、つまりペナ...
次に、最適な正則化パラメーター([math]\lambda[/math])の...
#geshi(rsplus){{
iris.lasso.cv = cv.glmnet(x, y)
}}
plotを用いて、[math]\lambda[/math]と平均二乗誤差の関係を...
#geshi(rsplus){{
plot(iris.lasso.cv)
}}
#ref(iris_lasso_lambda.png,nolink,50%)
今回の例では正則化パラメーターの値が小さいほうが平均二乗...
最も予測誤差が小さかったときの[math]\lambda[/math]の値は...
#geshi(rsplus){{
iris.lasso.cv$lambda.min
}}
#geshi(rsplus){{
[1] 0.008431423
}}
係数は変数betaに、定数項の値は変数a0に格納されています。
しかし、各[math]\lambda[/math]ごとの係数が格納された行列...
#geshi(rsplus){{
iris.lasso$beta[,which(iris.lasso$lambda==iris.lasso.cv$l...
}}
#geshi(rsplus){{
Sepal.Length Sepal.Width Petal.Width
0.7171234 -0.6267432 1.4504055
}}
#geshi(rsplus){{
iris.lasso$a0[which(iris.lasso$lambda==iris.lasso.cv$lamb...
}}
#geshi(rsplus){{
s57
-0.2557475
}}
重回帰と同様に表すと、学習された式を次のように書くことが...
\[y = -0.2557475 + 0.7171234 x_1 - 0.6267432 x_ 2 + 1.450...
単回帰、重回帰のときと同じように、平均二乗誤差を求め、予...
#geshi(rsplus){{
mean((y - predict(iris.lasso, x, s=iris.lasso.cv$lambda.m...
}}
#geshi(rsplus){{
[1] 0.09915377
}}
#geshi(rsplus){{
plot(y, predict(iris.lasso, x, s=iris.lasso.cv$lambda.min))
abline(a=0,b=1)
}}
#ref(iris_lasso_errors.png,nolink,50%)
*Ridge回帰 [#wa18ca12]
''Ridge回帰''は、目的関数の正則化項が[math]\frac{\lambda}...
つまり、Ridge回帰の目的関数は次のような形をしています。
\[\sum_{i=1}^N \left( y_i - (\beta_0 + \beta^\mathrm{T}\b...
Ridge回帰の正則化に用いられている[math]||\beta||_2^2[/mat...
\[||\beta||_2^2 = \sum_{i=1}^p \beta_i^2\]
Lasso回帰では正則化項のペナルティが係数の絶対値に比例して...
Ridge回帰を行うには、''glmnet関数''で学習する際にオプショ...
#geshi(rsplus){{
iris.ridge <- glmnet(x, y, alpha=0)
}}
後はLasso回帰のときと同じです。
まず、正則化項の大きさと係数の関係を表示します。
#geshi(rsplus){{
plot(iris.ridge)
}}
#ref(iris_ridge_regularization.png,nolink,50%)
Lasso回帰のグラフと比べると、Lasso回帰の方が係数を0にする...
次に、クロス・バリデーションによって最適な[math]\lambda[/...
#geshi(rsplus){{
iris.ridge.cv <- cv.glmnet(x, y, alpha=0)
iris.ridge.cv$lambda.min
}}
#geshi(rsplus){{
[1] 0.1694069
}}
#geshi(rsplus){{
iris.ridge$beta[,which(iris.ridge$lambda==iris.ridge.cv$l...
}}
#geshi(rsplus){{
Sepal.Length Sepal.Width Petal.Width
0.7882905 -0.6838069 1.2645594
}}
#geshi(rsplus){{
iris.ridge$a0[which(iris.ridge$lambda==iris.ridge.cv$lamb...
}}
#geshi(rsplus){{
s99
-0.2742464
}}
これまでと同様に表すと、学習された式は次のように書けます。
\[y = -0.2742464 + 0.7882905 x_1 - 0.6838069 x_2 + 1.2645...
最後に、平均二乗誤差を求め、予測誤差をプロットします。
#geshi(rsplus){{
mean((y - predict(iris.ridge, x, s=iris.ridge.cv$lambda.m...
}}
#geshi(rsplus){{
[1] 0.1083066
}}
#geshi(rsplus){{
plot(y, x, s=iris.ridge.cv$lambda.min))
abline(a=0,b=1)
}}
#ref(iris_ridge_errors.png,nolink,50%)
ちなみに、Lasso回帰のときと同じように[math]\lambda[/math]...
#geshi(rsplus){{
plot(iris.ridge.cv)
}}
#geshi(rsplus){{
以下にエラー axis(side = 3, at = sign.lambda * log(cvobj...
'at' と 'label' の長さが違います。 100 != 300
}}
*多項回帰 [#wd1c3279]
irisデータのように目的変数(Species)の値が離散値(カテゴ...
多項回帰を行うには、''glmnet関数''で学習する際にオプショ...
#geshi(rsplus){{
x1 <- iris$Sepal.Length
x2 <- iris$Sepal.Width
x3 <- iris$Petal.Length
x4 <- iris$Petal.Width
x <- cbind(x1, x2, x3, x4)
y <- iris$Species
iris.mnorm <- glmnet(x, y, family="multinomial")
plot(iris.mnorm)
}}
#ref(iris_multinomial_regularization.png,nolink,50%)
それから、クロス・バリデーションを用いて最適な[math]\lamb...
#geshi(rsplus){{
iris.mnorm.cv <- cv.glmnet(x, y, family="multinomial")
plot(iris.mnorm.cv)
}}
#ref(iris_multinomial_lambda.png,nolink,50%)
今回は[math]\lambda[/math]が小さすぎても大きすぎてもダメ...
最後に、予測結果を''table関数''を用いて表にします。
目的変数が離散値なので、''predict関数''にオプション''type...
#geshi(rsplus){{
table(y, predict(iris.mnorm, x, s=iris.mnorm.cv$lambda.mi...
}}
#geshi(rsplus){{
setosa versicolor virginica
setosa 50 0 0
versicolor 0 48 2
virginica 0 1 49
}}
setosaは50個全て正しく分類できていますが、versicolorのう...
*二項回帰 [#o625f55c]
目的変数の値が「はい」「いいえ」や「陽性」「陰性」など2種...
irisデータには、目的変数(Species)の値がsetosa、versicol...
まず、''subset関数''を用いてirisからSpeciesの値がsetosaか...
また、''factor関数''を用いてSpciesのレベルを再定義します。
#geshi(rsplus){{
iris2 <- subset(iris, Species=="setosa" | Species=="versi...
iris2$Species <- factor(iris2$Species, level=c("setosa","...
x1 <- iris2$Sepal.Length
x2 <- iris2$Sepal.Width
x3 <- iris2$Petal.Length
x4 <- iris2$Petal.Width
x <- cbind(x1, x2, x3, x4)
y <- iris2$Species
}}
二項回帰を行うには、''glmnet関数''で学習する際にオプショ...
後は多項回帰と同じです。
#geshi(rsplus){{
iris.binorm <- glmnet(x, y, family="binomial")
plot(iris.binorm)
}}
#ref(iris_binomial_regularization.png,nolink,50%)
#geshi(rsplus){{
iris.binorm.cv <- cv.glmnet(x, y, family="binomial")
plot(iris.binorm.cv)
}}
#ref(iris_binomial_lambda.png,nolink,50%)
#geshi(rsplus){{
table(y, predict(iris.binorm, x, s=iris.binorm.cv$lambda....
}}
#geshi(rsplus){{
setosa versicolor
setosa 50 0
versicolor 0 50
}}
こちらは全て正しく分類できました。
*糖尿病データへの応用 [#cac50579]
それでは、''diabetesデータ''に応用してみましょう。
このデータは''larsパッケージ''((Efron, Hastie, Johnstone...
#geshi(rsplus){{
install.packages("lars")
library(lars)
data(diabetes)
}}
diabetesデータには、442人の糖尿病患者から観測された年齢、...
まず、重回帰を行い、平均二乗誤差を求めます。
#geshi(rsplus){{
x <- diabetes$x
y <- diabetes$y
diabetes.mult <- lm(y ~ x)
mean ((y - predict(diabetes.mult))**2)
}}
#geshi(rsplus){{
[1] 2859.69
}}
次に、Lasso回帰を行い、正則化項と係数の関係をプロットしま...
#geshi(rsplus){{
diabetes.lasso <- glmnet(x, y)
plot(diabetes.lasso)
}}
#ref(diabetes_lasso_regularization.png,nolink,50%)
正則化が弱いときは(つまり、グラフの右の方では)、irisに...
クロス・バリデーションでLasso回帰の最適な[math]\lambda[/m...
#geshi(rsplus){{
diabetes.lasso.cv <- cv.glmnet(x, y)
diabetes.lasso.cv$lambda.min
}}
#geshi(rsplus){{
[1] 0.0735996
}}
#geshi(rsplus){{
plot(diabetes.lasso.cv)
}}
#ref(diabetes_lasso_lambda.png,nolink,50%)
Lasso回帰の平均二乗誤差を求めます。
#geshi(rsplus){{
mean((y - predict(diabetes.lasso, x, s=diabetes.lasso.cv$...
}}
#geshi(rsplus){{
[1] 2886.087
}}
最後に、Ridge回帰を行います。
#geshi(rsplus){{
diabetes.ridge <- x, y, alpha=0)
plot(diabetes.ridge)
}}
#ref(diabetes_ridge_regularization.png,nolink,50%)
クロス・バリデーションによってRidge回帰の最適な[math]\lam...
#geshi(rsplus){{
diabetes.ridge.cv <- cv.glmnet(x, y, alpha=0)
diabetes.ridge.cv$lambda.min
}}
#geshi(rsplus){{
[1] 4.516003
}}
Ridge回帰の平均二乗誤差を求めます。
#geshi(rsplus){{
mean((y - predict(diabetes.ridge, x, s=diabetes.ridge.cv$...
}}
#geshi(rsplus){{
[1] 2881.332
}}
*構造活性相関データへの応用 [#e9570b3e]
いよいよ、『Rによるバイオインフォマティクスデータ解析』の...
このデータセットは、''caretパッケージ''((Hothorn, Leisch,...
#geshi(rsplus){{
install.package("caret")
library(caret)
data(dhfr)
}}
dhfrデータセットの中を調べてみます。
#geshi(rsplus){{
dim(dhfr)
}}
#geshi(rsplus){{
[1] 325 229
}}
#geshi(rsplus){{
names(dhfr)
}}
#geshi(rsplus){{
[1] "Y" "moeGao_Abra_L" "moeG...
[4] "moeGao_Abra_acidity" "moeGao_Abra_basicity" "moeG...
[7] "moe2D_BCUT_PEOE_0" "moe2D_BCUT_PEOE_1" "moe2...
[10] "moe2D_BCUT_PEOE_3" "moe2D_BCUT_SLOGP_0" "moe2...
[13] "moe2D_BCUT_SLOGP_2" "moe2D_BCUT_SLOGP_3" "moe2...
[16] "moe2D_BCUT_SMR_1" "moe2D_BCUT_SMR_2" "moe2...
[19] "moe2D_FCharge" "moe2D_GCUT_PEOE_0" "moe2...
[22] "moe2D_GCUT_PEOE_2" "moe2D_GCUT_PEOE_3" "moe2...
[25] "moe2D_GCUT_SLOGP_1" "moe2D_GCUT_SLOGP_2" "moe2...
[28] "moe2D_GCUT_SMR_0" "moe2D_GCUT_SMR_1" "moe2...
[31] "moe2D_GCUT_SMR_3" "moe2D_Kier1" "moe2...
[34] "moe2D_Kier3" "moe2D_KierA1" "moe2...
[37] "moe2D_KierA3" "moe2D_KierFlex" "moe2...
[40] "moe2D_PEOE_PC..1" "moe2D_PEOE_RPC." "moe2...
[43] "moe2D_PEOE_VSA.0" "moe2D_PEOE_VSA.1" "moe2...
[46] "moe2D_PEOE_VSA.3" "moe2D_PEOE_VSA.4" "moe2...
[49] "moe2D_PEOE_VSA.6" "moe2D_PEOE_VSA.0.1" "moe2...
[52] "moe2D_PEOE_VSA.2.1" "moe2D_PEOE_VSA.3.1" "moe2...
[55] "moe2D_PEOE_VSA.5.1" "moe2D_PEOE_VSA.6.1" "moe2...
[58] "moe2D_PEOE_VSA_FNEG" "moe2D_PEOE_VSA_FPNEG" "moe2...
[61] "moe2D_PEOE_VSA_FPOS" "moe2D_PEOE_VSA_FPPOS" "moe2...
[64] "moe2D_PEOE_VSA_NEG" "moe2D_PEOE_VSA_PNEG" "moe2...
[67] "moe2D_PEOE_VSA_POS" "moe2D_PEOE_VSA_PPOS" "moe2...
[70] "moe2D_Q_PC..1" "moe2D_Q_RPC." "moe2...
[73] "moe2D_Q_VSA_FHYD" "moe2D_Q_VSA_FNEG" "moe2...
[76] "moe2D_Q_VSA_FPOL" "moe2D_Q_VSA_FPOS" "moe2...
[79] "moe2D_Q_VSA_HYD" "moe2D_Q_VSA_NEG" "moe2...
[82] "moe2D_Q_VSA_POL" "moe2D_Q_VSA_POS" "moe2...
[85] "moe2D_SMR" "moe2D_SMR_VSA0" "moe2...
[88] "moe2D_SMR_VSA2" "moe2D_SMR_VSA3" "moe2...
[91] "moe2D_SMR_VSA5" "moe2D_SMR_VSA6" "moe2...
[94] "moe2D_SlogP" "moe2D_SlogP_VSA0" "moe2...
[97] "moe2D_SlogP_VSA2" "moe2D_SlogP_VSA3" "moe2...
[100] "moe2D_SlogP_VSA5" "moe2D_SlogP_VSA6" "moe2...
[103] "moe2D_SlogP_VSA8" "moe2D_SlogP_VSA9" "moe2...
[106] "moe2D_VAdjEq" "moe2D_VAdjMa" "moe2...
[109] "moe2D_VDistMa" "moe2D_Weight" "moe2...
[112] "moe2D_a_ICM" "moe2D_a_acc" "moe2...
[115] "moe2D_a_base" "moe2D_a_count" "moe2...
[118] "moe2D_a_heavy" "moe2D_a_hyd" "moe2...
[121] "moe2D_a_nC" "moe2D_a_nCl" "moe2...
[124] "moe2D_a_nH" "moe2D_a_nN" "moe2...
[127] "moe2D_a_nS" "moe2D_apol" "moe2...
[130] "moe2D_b_1rotR" "moe2D_b_ar" "moe2...
[133] "moe2D_b_double" "moe2D_b_heavy" "moe2...
[136] "moe2D_b_rotR" "moe2D_b_single" "moe2...
[139] "moe2D_balabanJ" "moe2D_bpol" "moe2...
[142] "moe2D_chi0_C" "moe2D_chi0v" "moe2...
[145] "moe2D_chi1" "moe2D_chi1_C" "moe2...
[148] "moe2D_chi1v_C" "moeGao_chi2_C" "moeG...
[151] "moeGao_chi3c" "moeGao_chi3c_C" "moeG...
[154] "moeGao_chi3cv_C" "moeGao_chi3p" "moeG...
[157] "moeGao_chi3pv" "moeGao_chi3pv_C" "moeG...
[160] "moeGao_chi4c_C" "moeGao_chi4ca" "moeG...
[163] "moeGao_chi4cav" "moeGao_chi4cav_C" "moeG...
[166] "moeGao_chi4cv_C" "moeGao_chi4pc" "moeG...
[169] "moeGao_chi4pcv" "moeGao_chi4pcv_C" "moe2...
[172] "moe2D_density" "moe2D_diameter" "moe2...
[175] "moe2D_kS_aaN" "moe2D_kS_aaNH" "moe2...
[178] "moe2D_kS_aaS" "moe2D_kS_aaaC" "moe2...
[181] "moe2D_kS_aasN" "moe2D_kS_dO" "moe2...
[184] "moe2D_kS_ddssS" "moe2D_kS_dsCH" "moe2...
[187] "moe2D_kS_dssC" "moe2D_kS_sBr" "moe2...
[190] "moe2D_kS_sCl" "moe2D_kS_sF" "moe2...
[193] "moe2D_kS_sOH" "moe2D_kS_ssCH2" "moe2...
[196] "moe2D_kS_ssO" "moe2D_kS_ssS" "moe2...
[199] "moe2D_kS_sssN" "moe2D_kS_ssssC" "moe2...
[202] "moe2D_kS_tsC" "moe2D_lip_acc" "moe2...
[205] "moe2D_lip_druglike" "moe2D_lip_violation" "moe2...
[208] "moe2D_logS" "moe2D_mr" "moe2...
[211] "moe2D_opr_leadlike" "moe2D_opr_nring" "moe2...
[214] "moe2D_opr_violation" "moe2D_petitjean" "moe2...
[217] "moe2D_radius" "moe2D_rings" "moe2...
[220] "moe2D_vdw_vol" "moe2D_vsa_acc" "moe2...
[223] "moe2D_vsa_don" "moe2D_vsa_hyd" "moe2...
[226] "moe2D_vsa_pol" "moe2D_weinerPath" "moe2...
[229] "moe2D_zagreb"
}}
#geshi(rsplus){{
levels(dhfr$Y)
}}
#geshi(rsplus){{
[1] "active" "inactive"
}}
事例数が325、属性数が229です。
1個目の属性Yが目的変数で、あとの228個はすべて説明変数です。
これまでの例に比べると、説明変数がすごく多いですね。
Y、つまり目的変数の値はactiveとinactiveの2種類です。
目的変数の値が2種類のカテゴリーなので、二項回帰を行います。
#geshi(rsplus){{
dhfr.x <- as.matrix(dhfr[,-1])
dhfr.y <- dhfr[,1]
dhfr.binorm <- glmnet(dhfr.x, dhfr.y, family="binomial")
plot(dhfr.binorm)
}}
#ref(dhfr_binomial_regularization.png,nolink,50%)
説明変数の数が多いので、たくさんの線が描かれています。
意味はこれまでと同じで、一番左が最も大きく正則化して係数...
次に、クロス・バリデーションで最適な[math]\lambda[/math]...
#geshi(rsplus){{
dhfr.binorm.cv <- cv.glmnet(dhfr.x, dhfr.y, family="binom...
plot(dhfr.binorm.cv)
}}
#ref(dhfr_binomial_lambda.png,nolink,50%)
最後に、精度を求めます。
#geshi(rsplus){{
(result <- table(dhfr.y, predict(lasso, dhfr.x, s=dhfr.bi...
}}
#geshi(rsplus){{
active inactive
active 198 5
inactive 6 116
(result[1,1] + result[2,2]) / (result[1,1] + result[1,2] ...
[1] 0.9661538
}}
96.6%の正解率でした。
*まとめ [#o5e18fce]
線形回帰のモデルは次のような式で表されます。
\[y = \beta_0 + \beta^\mathrm{T}\bf{x}\]
ここで、[math]\beta[/math] と [math]\bf{x}[/math] は、そ...
次元数が [math]p=1[/math] のときを単回帰、[math]p>1[/math...
線形回帰では、二乗誤差の総和 [math]L(\beta, \bf{x})[/math...
\[L(\beta, \bf{x}) = \sum_{i=1}^N \left(y_i - (\beta_0 + ...
Lasso回帰やRidge回帰では、損失関数 [math]L(\beta, \bf{x})...
\[L(\beta, \bf{x}) + \lambda \textit{P}(\beta)\]
ここで、[math]\lambda[/math] は正則化の強さを決めるパラメ...
[math]\lambda=0[/math] のとき、正則化項が [math]0[/math] ...
また、[math]\lambda=0[/math] かつ [math]p=1[/math](次元...
glmnet関数では、[math]\lambda[/math] の値をオプションsに...
glmnet関数の正則化項 [math]P(\beta)[/math] は、次のような...
\[P(\beta) = (1 - \alpha)\frac{1}{2} ||\beta||_2^2 + \alp...
[math]\alpha[/math] は正則化の形を決めるパラメーターです。
[math]\alpha=1[/math] のときにLasso回帰、[math]\alpha=0[/...
glmnet関数では、この [math]\alpha[/math] の値をオプション...
デフォルトではalpha=1なので、alphaオプションを指定しない...
alpha=0と指定するとRidge回帰になります。
alphaの値を0と1の間に設定することもできます。
今回は、線形回帰の誤差関数に正則化項を加えたLasso回帰とRi...
また、多項回帰と二項回帰では、特に説明せずにデフォルトのL...
回帰分析についてもっと詳しく知りたい場合は、参考文献に上...
glmnet関数でロジスティック回帰、多項式回帰、ポアソン回帰...
*参考文献 [#i630eb1c]
この本の7.17節「LASSO」を参考にしています。
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0"...
}}
Lasso回帰とRidge回帰については、この本に詳しく載っていま...
この本は統計学習の教科書で、著者はglmnetパッケージの開発...
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0"...
}}
glmnetパッケージのソース、バイナリー、マニュアルがCRANか...
-[[glmnet: Lasso and elastic-net regularized generalized ...
終了行:
*はじめに [#gc93c04a]
ここでは、『Rによるバイオインフォマティクスデータ解析』の...
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0"...
}}
きちんと理解するために、まずは単回帰と重回帰という2種類の...
その後、多項回帰と二項回帰をやります。
*準備 [#e708ecca]
Rのインストールについては、次のページを見てください。
-[[MacでRを使う>機械学習/MacでRを使う]]
-[[WindowsでRを使う>機械学習/WindowsでRを使う]]
最初は、標準で使用できる''irisデータセット''を使います。
#geshi(rsplus){{
data(iris)
}}
このデータセットは、アヤメの種類(Species)をがく片の長さ...
長さと幅は連続値、種類はsetosa, versicolor, virginicaのい...
このデータセットには、150個のデータが含まれています。
ランダムに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...
}}
*線形回帰(単回帰) [#s1fe6159]
今から作るモデルによって説明したい変数を''目的変数''とい...
統計の分野では、目的変数を''従属変数''、説明変数を''独立...
まず、説明変数が1つの''単回帰''(single linear regression...
つまり、モデルが次のような式で表される直線をデータから求...
\[\hat{y} = \beta_0 + \beta_1 x_1\]
ここで、[math]\hat{y}[/math]がモデルによる''予測値''、[ma...
ある観測値 [math](x_i, y_i)[/math] に対して、モデルの予測...
\[ y_i - \hat{y_i} = y_i - (\beta_0 + \beta_1 x_{i, 1}) \]
線形回帰では、''残差の二乗の和''を最小化するような係数と...
つまり、単回帰では、[math]N[/math] 個の観測データに対して...
(求め方の説明は省きます。)
\[ \sum_{i=1}^N \left(y_i - (\beta_0 + \beta_1 x_{i,1})\r...
例として、がくの長さ(Petal.Length)を横軸、花びらの長さ...
#geshi(rsplus){{
plot(iris$Sepal.Length, iris$Petal.Length)
}}
#ref(iris_single_rawdata.png,nolink,50%)
どうやら、花びらの長さとがくの長さの間には正の相関(一方...
そこで、花びらの長さ(Sepal.Length)を [math]x_1[/math](...
Rでは、モデルの式を、記号''~''を使って「目的変数 ~ 説明変...
したがって、目的変数 [math]y[/math] をがくの長さ(Petal.L...
#geshi(rsplus){{
y <- iris$Petal.Length
x1 <- iris$Sepal.Length
}}
線形回帰は、標準で用意されている''lm関数''にモデル式を入...
#geshi(rsplus){{
iris.sgl <- lm(y ~ x1)
}}
学習された結果を表示してみると、次のようになります。
#geshi(rsplus){{
iris.sgl
}}
#geshi(rsplus){{
Call:
lm(formula = y ~ x1)
Coefficients:
(Intercept) x1
-7.101 1.858
}}
Coefficientsは係数を表し、Interceptは定数項の値(切片)を...
つまり、花びらの長さを [math]x_1[/math]、がくの長さを [ma...
\[ y = -7.101 + 1.858 x_1 \]
学習された単回帰モデルを図示します。
#geshi(rsplus){{
plot(x1, y)
abline(iris.sgl)
}}
#ref(iris_single_model.png,nolink,50%)
学習されたモデルを''summary関数''に渡すと、モデルに関する...
#geshi(rsplus){{
summary(iris.sgl)
}}
#geshi(rsplus){{
Call:
lm(formula = y ~ x1)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
x1 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ...
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16
}}
''Residuals''は、残差の最小値、中央値、最大値、四分位数を...
''Coefficients''は、''切片'' (Intercept) [math]\beta_0[/...
また、それぞれの推定値が 0 であるという帰無仮説に対するt...
p値が0.05を超えている場合は、その説明変数は有意水準5%で統...
''Multiple R-squared''は、''決定係数''で、このモデルがど...
単回帰では説明変数を1つしか使いませんが、説明変数を複数使...
''Adjusted R-squared''は、''自由度調整済み決定係数''で、...
説明変数の数が異なるモデルを比較するときには、自由度調整...
最後に、係数が全て 0 であるという帰無仮説に対するF検定の'...
p値が0.05を超えている場合は、このモデルは有意水準5%で統計...
次に、予測誤差を視覚化するため、それぞれのデータについて...
予測値は''predict関数''で求めることができます。
#geshi(rsplus){{
plot(y, predict(iris.sgl))
abline(a=0,b=1)
}}
#ref(iris_single_errors.png,nolink,50%)
この予測がどのくらい良いのかを表すために、''平均二乗誤差'...
\[\textit{MSE} = \frac{1}{N}\sum_{i=1}^{N} \left(y_i - (\...
#geshi(rsplus){{
mean((y - predict(iris.sgl))**2)
}}
#geshi(rsplus){{
[1] 0.743061
}}
*線形回帰(重回帰) [#bd685357]
説明変数が複数ある線形回帰を''重回帰''(multiple linear r...
重回帰のモデルは次のような式で表されます。
\[y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta...
ここで、[math]\beta[/math]と[math]\bf{x}[/math]は、それぞ...
また、単回帰のときと同様に、次式で表される二乗誤差の総和...
\[\sum_{i=1}^N \left(y_i - (\beta_0 + \beta^\mathrm{T}\bf...
ここでは、がくの長さ(Petal.Length)を目的変数 [math]y[/m...
#geshi(rsplus){{
y <- iris$Petal.Length
x1 <- iris$Sepal.Length
x2 <- iris$Sepal.Width
x3 <- iris$Petal.Width
}}
説明変数が複数あるときは、モデルを表すときにそれぞれの変...
ここでは、「y ~ x1+x2+x3」となります。
単回帰と同じように、''lm関数''を用いてモデルを学習します。
#geshi(rsplus){{
iris.mult <- lm(y ~ x1+x2+x3)
iris.mult$coefficients
}}
#geshi(rsplus){{
(Intercept) x1 x2 x3
-0.2627112 0.7291384 -0.6460124 1.4467934
}}
単回帰と同じように、学習された式は次のように書けます。
\[ y = -0.2627112 + 0.7291384 x_1 - 0.6460124 x_2 + 1.446...
単回帰と同じように、学習されたモデルを''summary関数''に渡...
#geshi(rsplus){{
summary(iris.mult)
}}
#geshi(rsplus){{
Call:
lm(formula = y ~ x1 + x2 + x3)
Residuals:
Min 1Q Median 3Q Max
-0.99333 -0.17656 -0.01004 0.18558 1.06909
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.26271 0.29741 -0.883 0.379
x1 0.72914 0.05832 12.502 <2e-16 ***
x2 -0.64601 0.06850 -9.431 <2e-16 ***
x3 1.44679 0.06761 21.399 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ...
Residual standard error: 0.319 on 146 degrees of freedom
Multiple R-squared: 0.968, Adjusted R-squared: 0.9674
F-statistic: 1473 on 3 and 146 DF, p-value: < 2.2e-16
}}
重回帰分析の場合はモデルが視覚化しにくいので、予測誤差だ...
#geshi(rsplus){{
plot(y, predict(iris.mult))
abline(a=0,b=1)
}}
#ref(iris_multiple_errors.png,nolink,50%)
単回帰と同様に、平均二乗誤差を求めます。
\[\textit{MSE} = \frac{1}{N} \sum_{i=1}^N (y_i - (\beta_0...
#geshi(rsplus){{
mean((y - predict(iris.mult))**2)
}}
#geshi(rsplus){{
[1] 0.09901965
}}
単回帰の場合よりも誤差がかなり小さくなっていることがわか...
*Lasso回帰 [#p3d71085]
''Lasso回帰''(Lasso regression)は、モデル式は重回帰と同...
\[\sum_{i=1}^N \left( y_i - (\beta_0 + \beta^\mathrm{T}\b...
ここで、[math]\lambda|\beta|[/math]を''正則化項''といい、...
Lasso回帰の正則化に用いられている[math]||\beta||_1[/math]...
\[||\beta||_1 = \sum_{i=1}^p|\beta_i|\]
つまり、L1ノルムは係数の絶対値の和を表します。
正則化パラメーター[math]\lambda[/math]の値が大きいほど0に...
Lasso回帰は''glmnetパッケージ''((Friedman, Hastie, and Ti...
そこで、glmnetパッケージをインストールして、使用します。
#geshi(rsplus){{
install.packages("glmnet")
library(glmnet)
}}
Lasso回帰を行うには、glmnetパッケージの''glmnet関数''を用...
引数には、説明変数行列と目的変数(そしてオプション)を与...
#geshi(rsplus){{
x <- cbind(x1, x2, x3)
iris.lasso <- glmnet(x, y)
}}
cbind関数は、列 (column) で結合するもので、x1, x2, x3 の...
plot関数を使うと、正則化項の大きさとそれぞれの係数の関係...
#geshi(rsplus){{
plot(iris.lasso)
}}
#ref(iris_lasso_regularization.png,nolink,50%)
黒色の折れ線が花びらの長さ(Sepal.Length)、赤色の折れ線...
縦軸はそれぞれの係数を表しています。
横軸の左端が最も正則化が強いとき、つまり、係数へのペナル...
最も正則化が最も強いときは、すべての係数が0になっています。
そこから正則化を弱めていくと、緑色の線が表すがくの幅(Pet...
その後、他の係数も0ではなくなり、正則化項が0、つまりペナ...
次に、最適な正則化パラメーター([math]\lambda[/math])の...
#geshi(rsplus){{
iris.lasso.cv = cv.glmnet(x, y)
}}
plotを用いて、[math]\lambda[/math]と平均二乗誤差の関係を...
#geshi(rsplus){{
plot(iris.lasso.cv)
}}
#ref(iris_lasso_lambda.png,nolink,50%)
今回の例では正則化パラメーターの値が小さいほうが平均二乗...
最も予測誤差が小さかったときの[math]\lambda[/math]の値は...
#geshi(rsplus){{
iris.lasso.cv$lambda.min
}}
#geshi(rsplus){{
[1] 0.008431423
}}
係数は変数betaに、定数項の値は変数a0に格納されています。
しかし、各[math]\lambda[/math]ごとの係数が格納された行列...
#geshi(rsplus){{
iris.lasso$beta[,which(iris.lasso$lambda==iris.lasso.cv$l...
}}
#geshi(rsplus){{
Sepal.Length Sepal.Width Petal.Width
0.7171234 -0.6267432 1.4504055
}}
#geshi(rsplus){{
iris.lasso$a0[which(iris.lasso$lambda==iris.lasso.cv$lamb...
}}
#geshi(rsplus){{
s57
-0.2557475
}}
重回帰と同様に表すと、学習された式を次のように書くことが...
\[y = -0.2557475 + 0.7171234 x_1 - 0.6267432 x_ 2 + 1.450...
単回帰、重回帰のときと同じように、平均二乗誤差を求め、予...
#geshi(rsplus){{
mean((y - predict(iris.lasso, x, s=iris.lasso.cv$lambda.m...
}}
#geshi(rsplus){{
[1] 0.09915377
}}
#geshi(rsplus){{
plot(y, predict(iris.lasso, x, s=iris.lasso.cv$lambda.min))
abline(a=0,b=1)
}}
#ref(iris_lasso_errors.png,nolink,50%)
*Ridge回帰 [#wa18ca12]
''Ridge回帰''は、目的関数の正則化項が[math]\frac{\lambda}...
つまり、Ridge回帰の目的関数は次のような形をしています。
\[\sum_{i=1}^N \left( y_i - (\beta_0 + \beta^\mathrm{T}\b...
Ridge回帰の正則化に用いられている[math]||\beta||_2^2[/mat...
\[||\beta||_2^2 = \sum_{i=1}^p \beta_i^2\]
Lasso回帰では正則化項のペナルティが係数の絶対値に比例して...
Ridge回帰を行うには、''glmnet関数''で学習する際にオプショ...
#geshi(rsplus){{
iris.ridge <- glmnet(x, y, alpha=0)
}}
後はLasso回帰のときと同じです。
まず、正則化項の大きさと係数の関係を表示します。
#geshi(rsplus){{
plot(iris.ridge)
}}
#ref(iris_ridge_regularization.png,nolink,50%)
Lasso回帰のグラフと比べると、Lasso回帰の方が係数を0にする...
次に、クロス・バリデーションによって最適な[math]\lambda[/...
#geshi(rsplus){{
iris.ridge.cv <- cv.glmnet(x, y, alpha=0)
iris.ridge.cv$lambda.min
}}
#geshi(rsplus){{
[1] 0.1694069
}}
#geshi(rsplus){{
iris.ridge$beta[,which(iris.ridge$lambda==iris.ridge.cv$l...
}}
#geshi(rsplus){{
Sepal.Length Sepal.Width Petal.Width
0.7882905 -0.6838069 1.2645594
}}
#geshi(rsplus){{
iris.ridge$a0[which(iris.ridge$lambda==iris.ridge.cv$lamb...
}}
#geshi(rsplus){{
s99
-0.2742464
}}
これまでと同様に表すと、学習された式は次のように書けます。
\[y = -0.2742464 + 0.7882905 x_1 - 0.6838069 x_2 + 1.2645...
最後に、平均二乗誤差を求め、予測誤差をプロットします。
#geshi(rsplus){{
mean((y - predict(iris.ridge, x, s=iris.ridge.cv$lambda.m...
}}
#geshi(rsplus){{
[1] 0.1083066
}}
#geshi(rsplus){{
plot(y, x, s=iris.ridge.cv$lambda.min))
abline(a=0,b=1)
}}
#ref(iris_ridge_errors.png,nolink,50%)
ちなみに、Lasso回帰のときと同じように[math]\lambda[/math]...
#geshi(rsplus){{
plot(iris.ridge.cv)
}}
#geshi(rsplus){{
以下にエラー axis(side = 3, at = sign.lambda * log(cvobj...
'at' と 'label' の長さが違います。 100 != 300
}}
*多項回帰 [#wd1c3279]
irisデータのように目的変数(Species)の値が離散値(カテゴ...
多項回帰を行うには、''glmnet関数''で学習する際にオプショ...
#geshi(rsplus){{
x1 <- iris$Sepal.Length
x2 <- iris$Sepal.Width
x3 <- iris$Petal.Length
x4 <- iris$Petal.Width
x <- cbind(x1, x2, x3, x4)
y <- iris$Species
iris.mnorm <- glmnet(x, y, family="multinomial")
plot(iris.mnorm)
}}
#ref(iris_multinomial_regularization.png,nolink,50%)
それから、クロス・バリデーションを用いて最適な[math]\lamb...
#geshi(rsplus){{
iris.mnorm.cv <- cv.glmnet(x, y, family="multinomial")
plot(iris.mnorm.cv)
}}
#ref(iris_multinomial_lambda.png,nolink,50%)
今回は[math]\lambda[/math]が小さすぎても大きすぎてもダメ...
最後に、予測結果を''table関数''を用いて表にします。
目的変数が離散値なので、''predict関数''にオプション''type...
#geshi(rsplus){{
table(y, predict(iris.mnorm, x, s=iris.mnorm.cv$lambda.mi...
}}
#geshi(rsplus){{
setosa versicolor virginica
setosa 50 0 0
versicolor 0 48 2
virginica 0 1 49
}}
setosaは50個全て正しく分類できていますが、versicolorのう...
*二項回帰 [#o625f55c]
目的変数の値が「はい」「いいえ」や「陽性」「陰性」など2種...
irisデータには、目的変数(Species)の値がsetosa、versicol...
まず、''subset関数''を用いてirisからSpeciesの値がsetosaか...
また、''factor関数''を用いてSpciesのレベルを再定義します。
#geshi(rsplus){{
iris2 <- subset(iris, Species=="setosa" | Species=="versi...
iris2$Species <- factor(iris2$Species, level=c("setosa","...
x1 <- iris2$Sepal.Length
x2 <- iris2$Sepal.Width
x3 <- iris2$Petal.Length
x4 <- iris2$Petal.Width
x <- cbind(x1, x2, x3, x4)
y <- iris2$Species
}}
二項回帰を行うには、''glmnet関数''で学習する際にオプショ...
後は多項回帰と同じです。
#geshi(rsplus){{
iris.binorm <- glmnet(x, y, family="binomial")
plot(iris.binorm)
}}
#ref(iris_binomial_regularization.png,nolink,50%)
#geshi(rsplus){{
iris.binorm.cv <- cv.glmnet(x, y, family="binomial")
plot(iris.binorm.cv)
}}
#ref(iris_binomial_lambda.png,nolink,50%)
#geshi(rsplus){{
table(y, predict(iris.binorm, x, s=iris.binorm.cv$lambda....
}}
#geshi(rsplus){{
setosa versicolor
setosa 50 0
versicolor 0 50
}}
こちらは全て正しく分類できました。
*糖尿病データへの応用 [#cac50579]
それでは、''diabetesデータ''に応用してみましょう。
このデータは''larsパッケージ''((Efron, Hastie, Johnstone...
#geshi(rsplus){{
install.packages("lars")
library(lars)
data(diabetes)
}}
diabetesデータには、442人の糖尿病患者から観測された年齢、...
まず、重回帰を行い、平均二乗誤差を求めます。
#geshi(rsplus){{
x <- diabetes$x
y <- diabetes$y
diabetes.mult <- lm(y ~ x)
mean ((y - predict(diabetes.mult))**2)
}}
#geshi(rsplus){{
[1] 2859.69
}}
次に、Lasso回帰を行い、正則化項と係数の関係をプロットしま...
#geshi(rsplus){{
diabetes.lasso <- glmnet(x, y)
plot(diabetes.lasso)
}}
#ref(diabetes_lasso_regularization.png,nolink,50%)
正則化が弱いときは(つまり、グラフの右の方では)、irisに...
クロス・バリデーションでLasso回帰の最適な[math]\lambda[/m...
#geshi(rsplus){{
diabetes.lasso.cv <- cv.glmnet(x, y)
diabetes.lasso.cv$lambda.min
}}
#geshi(rsplus){{
[1] 0.0735996
}}
#geshi(rsplus){{
plot(diabetes.lasso.cv)
}}
#ref(diabetes_lasso_lambda.png,nolink,50%)
Lasso回帰の平均二乗誤差を求めます。
#geshi(rsplus){{
mean((y - predict(diabetes.lasso, x, s=diabetes.lasso.cv$...
}}
#geshi(rsplus){{
[1] 2886.087
}}
最後に、Ridge回帰を行います。
#geshi(rsplus){{
diabetes.ridge <- x, y, alpha=0)
plot(diabetes.ridge)
}}
#ref(diabetes_ridge_regularization.png,nolink,50%)
クロス・バリデーションによってRidge回帰の最適な[math]\lam...
#geshi(rsplus){{
diabetes.ridge.cv <- cv.glmnet(x, y, alpha=0)
diabetes.ridge.cv$lambda.min
}}
#geshi(rsplus){{
[1] 4.516003
}}
Ridge回帰の平均二乗誤差を求めます。
#geshi(rsplus){{
mean((y - predict(diabetes.ridge, x, s=diabetes.ridge.cv$...
}}
#geshi(rsplus){{
[1] 2881.332
}}
*構造活性相関データへの応用 [#e9570b3e]
いよいよ、『Rによるバイオインフォマティクスデータ解析』の...
このデータセットは、''caretパッケージ''((Hothorn, Leisch,...
#geshi(rsplus){{
install.package("caret")
library(caret)
data(dhfr)
}}
dhfrデータセットの中を調べてみます。
#geshi(rsplus){{
dim(dhfr)
}}
#geshi(rsplus){{
[1] 325 229
}}
#geshi(rsplus){{
names(dhfr)
}}
#geshi(rsplus){{
[1] "Y" "moeGao_Abra_L" "moeG...
[4] "moeGao_Abra_acidity" "moeGao_Abra_basicity" "moeG...
[7] "moe2D_BCUT_PEOE_0" "moe2D_BCUT_PEOE_1" "moe2...
[10] "moe2D_BCUT_PEOE_3" "moe2D_BCUT_SLOGP_0" "moe2...
[13] "moe2D_BCUT_SLOGP_2" "moe2D_BCUT_SLOGP_3" "moe2...
[16] "moe2D_BCUT_SMR_1" "moe2D_BCUT_SMR_2" "moe2...
[19] "moe2D_FCharge" "moe2D_GCUT_PEOE_0" "moe2...
[22] "moe2D_GCUT_PEOE_2" "moe2D_GCUT_PEOE_3" "moe2...
[25] "moe2D_GCUT_SLOGP_1" "moe2D_GCUT_SLOGP_2" "moe2...
[28] "moe2D_GCUT_SMR_0" "moe2D_GCUT_SMR_1" "moe2...
[31] "moe2D_GCUT_SMR_3" "moe2D_Kier1" "moe2...
[34] "moe2D_Kier3" "moe2D_KierA1" "moe2...
[37] "moe2D_KierA3" "moe2D_KierFlex" "moe2...
[40] "moe2D_PEOE_PC..1" "moe2D_PEOE_RPC." "moe2...
[43] "moe2D_PEOE_VSA.0" "moe2D_PEOE_VSA.1" "moe2...
[46] "moe2D_PEOE_VSA.3" "moe2D_PEOE_VSA.4" "moe2...
[49] "moe2D_PEOE_VSA.6" "moe2D_PEOE_VSA.0.1" "moe2...
[52] "moe2D_PEOE_VSA.2.1" "moe2D_PEOE_VSA.3.1" "moe2...
[55] "moe2D_PEOE_VSA.5.1" "moe2D_PEOE_VSA.6.1" "moe2...
[58] "moe2D_PEOE_VSA_FNEG" "moe2D_PEOE_VSA_FPNEG" "moe2...
[61] "moe2D_PEOE_VSA_FPOS" "moe2D_PEOE_VSA_FPPOS" "moe2...
[64] "moe2D_PEOE_VSA_NEG" "moe2D_PEOE_VSA_PNEG" "moe2...
[67] "moe2D_PEOE_VSA_POS" "moe2D_PEOE_VSA_PPOS" "moe2...
[70] "moe2D_Q_PC..1" "moe2D_Q_RPC." "moe2...
[73] "moe2D_Q_VSA_FHYD" "moe2D_Q_VSA_FNEG" "moe2...
[76] "moe2D_Q_VSA_FPOL" "moe2D_Q_VSA_FPOS" "moe2...
[79] "moe2D_Q_VSA_HYD" "moe2D_Q_VSA_NEG" "moe2...
[82] "moe2D_Q_VSA_POL" "moe2D_Q_VSA_POS" "moe2...
[85] "moe2D_SMR" "moe2D_SMR_VSA0" "moe2...
[88] "moe2D_SMR_VSA2" "moe2D_SMR_VSA3" "moe2...
[91] "moe2D_SMR_VSA5" "moe2D_SMR_VSA6" "moe2...
[94] "moe2D_SlogP" "moe2D_SlogP_VSA0" "moe2...
[97] "moe2D_SlogP_VSA2" "moe2D_SlogP_VSA3" "moe2...
[100] "moe2D_SlogP_VSA5" "moe2D_SlogP_VSA6" "moe2...
[103] "moe2D_SlogP_VSA8" "moe2D_SlogP_VSA9" "moe2...
[106] "moe2D_VAdjEq" "moe2D_VAdjMa" "moe2...
[109] "moe2D_VDistMa" "moe2D_Weight" "moe2...
[112] "moe2D_a_ICM" "moe2D_a_acc" "moe2...
[115] "moe2D_a_base" "moe2D_a_count" "moe2...
[118] "moe2D_a_heavy" "moe2D_a_hyd" "moe2...
[121] "moe2D_a_nC" "moe2D_a_nCl" "moe2...
[124] "moe2D_a_nH" "moe2D_a_nN" "moe2...
[127] "moe2D_a_nS" "moe2D_apol" "moe2...
[130] "moe2D_b_1rotR" "moe2D_b_ar" "moe2...
[133] "moe2D_b_double" "moe2D_b_heavy" "moe2...
[136] "moe2D_b_rotR" "moe2D_b_single" "moe2...
[139] "moe2D_balabanJ" "moe2D_bpol" "moe2...
[142] "moe2D_chi0_C" "moe2D_chi0v" "moe2...
[145] "moe2D_chi1" "moe2D_chi1_C" "moe2...
[148] "moe2D_chi1v_C" "moeGao_chi2_C" "moeG...
[151] "moeGao_chi3c" "moeGao_chi3c_C" "moeG...
[154] "moeGao_chi3cv_C" "moeGao_chi3p" "moeG...
[157] "moeGao_chi3pv" "moeGao_chi3pv_C" "moeG...
[160] "moeGao_chi4c_C" "moeGao_chi4ca" "moeG...
[163] "moeGao_chi4cav" "moeGao_chi4cav_C" "moeG...
[166] "moeGao_chi4cv_C" "moeGao_chi4pc" "moeG...
[169] "moeGao_chi4pcv" "moeGao_chi4pcv_C" "moe2...
[172] "moe2D_density" "moe2D_diameter" "moe2...
[175] "moe2D_kS_aaN" "moe2D_kS_aaNH" "moe2...
[178] "moe2D_kS_aaS" "moe2D_kS_aaaC" "moe2...
[181] "moe2D_kS_aasN" "moe2D_kS_dO" "moe2...
[184] "moe2D_kS_ddssS" "moe2D_kS_dsCH" "moe2...
[187] "moe2D_kS_dssC" "moe2D_kS_sBr" "moe2...
[190] "moe2D_kS_sCl" "moe2D_kS_sF" "moe2...
[193] "moe2D_kS_sOH" "moe2D_kS_ssCH2" "moe2...
[196] "moe2D_kS_ssO" "moe2D_kS_ssS" "moe2...
[199] "moe2D_kS_sssN" "moe2D_kS_ssssC" "moe2...
[202] "moe2D_kS_tsC" "moe2D_lip_acc" "moe2...
[205] "moe2D_lip_druglike" "moe2D_lip_violation" "moe2...
[208] "moe2D_logS" "moe2D_mr" "moe2...
[211] "moe2D_opr_leadlike" "moe2D_opr_nring" "moe2...
[214] "moe2D_opr_violation" "moe2D_petitjean" "moe2...
[217] "moe2D_radius" "moe2D_rings" "moe2...
[220] "moe2D_vdw_vol" "moe2D_vsa_acc" "moe2...
[223] "moe2D_vsa_don" "moe2D_vsa_hyd" "moe2...
[226] "moe2D_vsa_pol" "moe2D_weinerPath" "moe2...
[229] "moe2D_zagreb"
}}
#geshi(rsplus){{
levels(dhfr$Y)
}}
#geshi(rsplus){{
[1] "active" "inactive"
}}
事例数が325、属性数が229です。
1個目の属性Yが目的変数で、あとの228個はすべて説明変数です。
これまでの例に比べると、説明変数がすごく多いですね。
Y、つまり目的変数の値はactiveとinactiveの2種類です。
目的変数の値が2種類のカテゴリーなので、二項回帰を行います。
#geshi(rsplus){{
dhfr.x <- as.matrix(dhfr[,-1])
dhfr.y <- dhfr[,1]
dhfr.binorm <- glmnet(dhfr.x, dhfr.y, family="binomial")
plot(dhfr.binorm)
}}
#ref(dhfr_binomial_regularization.png,nolink,50%)
説明変数の数が多いので、たくさんの線が描かれています。
意味はこれまでと同じで、一番左が最も大きく正則化して係数...
次に、クロス・バリデーションで最適な[math]\lambda[/math]...
#geshi(rsplus){{
dhfr.binorm.cv <- cv.glmnet(dhfr.x, dhfr.y, family="binom...
plot(dhfr.binorm.cv)
}}
#ref(dhfr_binomial_lambda.png,nolink,50%)
最後に、精度を求めます。
#geshi(rsplus){{
(result <- table(dhfr.y, predict(lasso, dhfr.x, s=dhfr.bi...
}}
#geshi(rsplus){{
active inactive
active 198 5
inactive 6 116
(result[1,1] + result[2,2]) / (result[1,1] + result[1,2] ...
[1] 0.9661538
}}
96.6%の正解率でした。
*まとめ [#o5e18fce]
線形回帰のモデルは次のような式で表されます。
\[y = \beta_0 + \beta^\mathrm{T}\bf{x}\]
ここで、[math]\beta[/math] と [math]\bf{x}[/math] は、そ...
次元数が [math]p=1[/math] のときを単回帰、[math]p>1[/math...
線形回帰では、二乗誤差の総和 [math]L(\beta, \bf{x})[/math...
\[L(\beta, \bf{x}) = \sum_{i=1}^N \left(y_i - (\beta_0 + ...
Lasso回帰やRidge回帰では、損失関数 [math]L(\beta, \bf{x})...
\[L(\beta, \bf{x}) + \lambda \textit{P}(\beta)\]
ここで、[math]\lambda[/math] は正則化の強さを決めるパラメ...
[math]\lambda=0[/math] のとき、正則化項が [math]0[/math] ...
また、[math]\lambda=0[/math] かつ [math]p=1[/math](次元...
glmnet関数では、[math]\lambda[/math] の値をオプションsに...
glmnet関数の正則化項 [math]P(\beta)[/math] は、次のような...
\[P(\beta) = (1 - \alpha)\frac{1}{2} ||\beta||_2^2 + \alp...
[math]\alpha[/math] は正則化の形を決めるパラメーターです。
[math]\alpha=1[/math] のときにLasso回帰、[math]\alpha=0[/...
glmnet関数では、この [math]\alpha[/math] の値をオプション...
デフォルトではalpha=1なので、alphaオプションを指定しない...
alpha=0と指定するとRidge回帰になります。
alphaの値を0と1の間に設定することもできます。
今回は、線形回帰の誤差関数に正則化項を加えたLasso回帰とRi...
また、多項回帰と二項回帰では、特に説明せずにデフォルトのL...
回帰分析についてもっと詳しく知りたい場合は、参考文献に上...
glmnet関数でロジスティック回帰、多項式回帰、ポアソン回帰...
*参考文献 [#i630eb1c]
この本の7.17節「LASSO」を参考にしています。
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0"...
}}
Lasso回帰とRidge回帰については、この本に詳しく載っていま...
この本は統計学習の教科書で、著者はglmnetパッケージの開発...
#html{{
<iframe style="width:120px;height:240px;" marginwidth="0"...
}}
glmnetパッケージのソース、バイナリー、マニュアルがCRANか...
-[[glmnet: Lasso and elastic-net regularized generalized ...
ページ名: