バイオ・データ・マイニング/Rで回帰分析する のバックアップソース(No.4)

• バックアップ一覧
• 差分 を表示
• 現在との差分 を表示
• バックアップ を表示
• バイオ・データ・マイニング/Rで回帰分析する へ行く。
  • 1 (2011-11-07 (月) 13:30:05)
  • 2 (2011-11-07 (月) 16:06:30)
  • 3 (2011-11-07 (月) 19:24:10)
  • 4 (2011-11-08 (火) 02:31:29)
  • 5 (2011-11-08 (火) 10:08:16)
  • 6 (2011-11-08 (火) 16:54:20)
  • 7 (2011-11-08 (火) 19:54:10)
  • 8 (2011-11-10 (木) 15:58:47)
  • 9 (2011-11-11 (金) 12:46:43)
  • 10 (2011-11-12 (土) 17:00:14)
  • 11 (2011-11-15 (火) 08:05:02)
  • 12 (2013-11-19 (火) 15:52:23)
  • 13 (2013-12-05 (木) 13:03:50)
  • 14 (2013-12-05 (木) 13:03:50)
  • 15 (2013-12-05 (木) 13:03:50)
  • 16 (2016-12-15 (木) 08:40:12)
  • 17 (2017-12-14 (木) 20:11:02)
  • 18 (2017-12-15 (金) 12:11:37)
  • 19 (2017-12-15 (金) 12:11:37)
  • 20 (2020-01-07 (火) 19:29:34)
  • 21 (2020-01-15 (水) 16:27:42)

*はじめに [#gc93c04a]

『Rによるバイオインフォマティクスデータ解析』の7.17節「LASSO」を参考にして，回帰分析をします．

#html{{
<iframe src="http://rcm-jp.amazon.co.jp/e/cm?t=tohgoroh-22&o=9&p=8&l=as1&asins=4320057082&ref=tf_til&fc1=444B4C&IS2=1&lt1=_blank&m=amazon&lc1=444B4C&bc1=FFFFFF&bg1=FFFFFF&f=ifr" style="width:120px;height:240px;" scrolling="no" marginwidth="0" marginheight="0" frameborder="0"></iframe>
}}

きちんと理解するために，まずは線形回帰をやります．
また，ついでにRidge回帰もやります．


*準備 [#e708ecca]

まずは，標準で使用できるirisデータを使います．
#geshi(sh){{
> data(iris)
}}

このデータは，ユリの種類（Species）を花びらの長さ（Sepal.Length），幅（Lepal.Width），がくの長さ（Petal.Length），幅（Petal.Width）によって分類する問題です．
長さと幅は連続値，種類はsetosa, versicolor, virginicaのいずれかをとる離散値です．

150個のデータが含まれています．
最初の10個を表示すると，次のようになります．
#geshi(sh){{
> iris[1:10,]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1           5.1         3.5          1.4         0.2  setosa
2           4.9         3.0          1.4         0.2  setosa
3           4.7         3.2          1.3         0.2  setosa
4           4.6         3.1          1.5         0.2  setosa
5           5.0         3.6          1.4         0.2  setosa
6           5.4         3.9          1.7         0.4  setosa
7           4.6         3.4          1.4         0.3  setosa
8           5.0         3.4          1.5         0.2  setosa
9           4.4         2.9          1.4         0.2  setosa
10          4.9         3.1          1.5         0.1  setosa
}}

ランダムにサンプリングした10個を表示すると，このようになります．
#geshi(sh){{
> iris[sort(sample(1:150,10)),]
    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
}}



*線形回帰（単回帰） [#s1fe6159]

まず，説明変数が1つの線形回帰（単回帰）を行います．
つまり，モデルが次のような式で表される直線をデータから求めます．
\[y = \beta_0 + \beta_1 x_1\]
ここで，[math]y[/math]が被説明変数（目的変数），[math]x[/math]が説明変数，[math]\beta_0[/math]が定数項，[math]\beta_1[/math]が（[math]x[/math]の）係数です．

線形回帰では，訓練データに対する二乗誤差を最小化するような係数および定数項を求めます．
つまり，単回帰では，[math]N[/math]個の観測データに対して，
\[\sum_{i=1}^N ||y_i - (\beta_0 + \beta_1 x_{i,1})||^2\]
を最小化するよう[math]\beta_0[/math]と[math]\beta_1[/math]を求めます．
（求め方の説明は省きます．）

ここでは，がくの長さ（Petal.Length）を目的変数，花びらの長さ（Sepal.Length）を説明変数としてモデルを学習しましょう．

この関係をプロットすると次のようになります．
#geshi(sh){{
> plot(iris$Sepal.Length, iris$Petal.Length)
}}
#ref(iris_single_rawdata.png,nolink,50%)

Rでは，モデルの式を，記号~を使って「目的変数 ~ 説明変数」というように表します．
したがって，目的変数をがくの長さ（Petal.Length），説明変数を花びらの長さ（Sepal.Length）とするとき，モデルの式は「iris$Petal.Length ~ iris$Sepal.Length」となります．

線形回帰は，標準で用意されているlm関数で行います．
#geshi(sh){{
> Single <- lm(iris$Petal.Length ~ iris$Sepal.Length)
}}

学習された結果を表示してみると，次のようになります．
#geshi(sh){{
> Single$coefficients
      (Intercept) iris$Sepal.Length 
        -7.101443          1.858433 
}}
Coefficientsは係数を表し，Interceptは定数項（[math]y[/math]切片）の値を表します．
つまり，花びらの長さを[math]x[/math]，がくの長さを[math]y[/math]とすると，学習された式は次のように書けます．
\[y =  -7.101443 + 1.858433 x_1\]

学習された線形回帰モデルを図示します．
#geshi(sh){{
> plot(iris$Sepal.Length, iris$Petal.Length)
> abline(Single)
}}

すると，次のようなグラフがプロットされます．
#ref(iris_single_model.png,nolink,50%)
単回帰だとモデルが視覚化しやすいです．

次に，予測誤差を視覚化するため，それぞれのデータについて，横軸に実際のがくの長さ（Petal.Length），縦軸にその予測値をとってグラフを表示します．
予測値はpredict関数で求めることができます．
#geshi(sh){{
> plot(iris$Petal.Length, predict(Single))
> abline(a=0,b=1)
}}
#ref(iris_single_errors.png,nolink,50%)

この予測がどのくらい良いのかを表すために，平均二乗誤差（実際の値と予測値の差の二乗の平均）を求めます．
#geshi(sh){{
> summary((iris$Petal.Length - predict(Single))**2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00007 0.11160 0.35950 0.74310 0.97580 6.22600 
}}
Meanのところが平均二乗誤差です．



*線形回帰（重回帰） [#bd685357]

説明変数が複数ある回帰を重回帰といいます．

重回帰のモデルは次のような式で表されます．
\[y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p\ = \beta_0 + \beta^\mathrm{T}\bf{x}\]
ここで，[math]\beta[/math]と[math]\bf{x}[/math]は，それぞれ[math]\beta_1, \dots, \beta_p[/math]と[math]x_1, \dots, x_p[/math]を要素とするベクトルです．

また，単回帰のときと同様に，二乗誤差
\[\sum_{i=1}^N || y_i - (\beta_0 + \beta^\mathrm{T}\bf{x}_i)||^2\]
を最小化するような[math]\beta_0, \dots, \beta_p[/math]を求めます．

ここでは，がくの長さ（Petal.Length）を目的変数，花びらの長さ（Sepal.Length），幅（Sepal.Width），がくの幅（Petal.Width）を説明変数として重回帰を行いましょう．

説明変数が複数あるときは，モデルを表すときにそれぞれの変数を記号+で連結して表します．
ここでは，「iris$Petal.Length ~ iris$Sepal.Length+iris$Sepal.Width+iris$Petal.Width」となります．

単回帰と同じように，lm関数を用いて線形回帰を行います．
#geshi(sh){{
> Multi <- lm(iris$Petal.Length ~ iris$Sepal.Length+iris$Sepal.Width+iris$Petal.Width)
> Multi$coefficients
      (Intercept) iris$Sepal.Length  iris$Sepal.Width  iris$Petal.Width 
       -0.2627112         0.7291384        -0.6460124         1.4467934 
}}
これも，花びらの長さ（Sepal.Length）を[math]x_1[/math]，幅（Sepal.Width）を[math]x_2[/math]，がくの幅（Petal.Width）を[math]x_3[/math]，がくの長さ（Petal.Length）を[math]y[/math]とすると，学習された式は次のように書けます．
\[y = -0.2627112 + 0.7291384 x_1 - 0.6460124 x_2 + 1.4467934 x_3\]


重回帰分析の場合はモデルが視覚化しにくいので，予測誤差だけ視覚化します．
#geshi(sh){{
> plot(iris$Petal.Length, predict(Multi))
> abline(a=0,b=1)
}}
#ref(iris_multiple_errors.png,nolink,50%)

単回帰と同様に，平均二乗誤差を求めます．
#geshi(sh){{
> summary((iris$Petal.Length - predict(Multi))**2)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000007 0.0058040 0.0328400 0.0990200 0.1324000 1.1430000 
}}
単回帰の場合よりも誤差がかなり小さくなっていることがわかります．



*LASSO回帰 [#p3d71085]

LASSO回帰は，モデル式は重回帰と同じですが，最小化する目的関数が次のような形をしています．
\[\sum_{i=1}^N || y_i - (\beta_0 + \beta^\mathrm{T}\bf{x}_i)||^2 + \lambda |\beta|\]
ここで，[math]\lambda|\beta|[/math]を正則化項といい，[math]\lambda[/math]を正則化パラメーターといいます．

LASSO回帰の正則化をL1ノルムといい，次のように計算します．
\[|\beta| = \sum_{i=1}^p|\beta_i|\]
つまり，L1ノルムは係数の絶対値の和を表します．
正則化パラメーター[math]\lambda[/math]の値が大きいほど0に近い係数を好むようになります．


LASSO回帰はglmnetパッケージに含まれています．
glmnetパッケージは
そこで，glmnetパッケージをインストールして，使用します．

#geshi(sh){{
> install.packages("glmnet")
> library(glmnet)
 要求されたパッケージ Matrix をロード中です 
 要求されたパッケージ lattice をロード中です 

 次のパッケージを付け加えます: 'Matrix' 

The following object(s) are masked from 'package:base':

    det

Loaded glmnet 1.7.1

 警告メッセージ： 
1:  以前のインポート ‘head’ を置き換えています（‘utils’のロード中）  
2:  以前のインポート ‘tail’ を置き換えています（‘utils’のロード中）  
>
}}
私の環境ではこのような警告が出ますが，これは「utilsパッケージのhead関数とtail関数をglmnetパッケージのhead関数とtail関数で置き換えましたよ」という意味なので，気にしなくていいです．

LASSO回帰は，glmnetパッケージのglmnet関数を用いて行います．
引数には，説明変数行列と目的変数（そしてオプション）を与えます．
#geshi(sh){{
> LASSO <- glmnet(as.matrix(iris[,-5][,-3]), iris$Petal.Length)
}}
iris[,-5][,-3]は，irisデータの5列目と3列目を除いたもので，これをas.matrix関数を用いて行列として与えています．

plot関数を使うと，正則化項の大きさとそれぞれの係数の関係をプロットすることができます．
#geshi(sh){{
> plot(LASSO)
}}
#ref(iris_lasso_regularization.png,nolink,50%)
黒色の折れ線が花びらの長さ（Sepal.Length），赤色の折れ線が花びらの幅（Sepal.Width），緑色の折れ線ががくの幅（Petal.Width）を表しています．

縦軸はそれぞれの係数を表しています．
横軸の左端が最も正則化が強いとき，つまり，係数が大きいときのペナルティを最も大きくしたときであり，右端が最も正則化が弱いとき，つまり，係数の大きさに対するペナルティなしのときです．

最も正則化が強いときは，すべての係数が0になっています．
また，ペナルティなしのときの係数は，重回帰のときと同じになっています．

次に，最適な正則化パラメーター（[math]\lambda[/math]）の値を交差検定（クロス・バリデーション）を用いて求めます．
#geshi(sh){{
> LASSO.cv = cv.glmnet(as.matrix(iris[,-5][,-3]), iris$Petal.Length)
}}

plotを用いて，[math]\lambda[/math]と平均二乗誤差の関係を表示すると次のようになります．
#geshi(sh){{
> plot(LASSO.cv)
}}
#ref(iris_lasso_lambda.png,nolink,50%)
今回の例では正則化パラメーターの値が小さいほうが平均二乗誤差が小さくなりましたが，いつもこうなるわけではありません．

最も予測誤差が小さかったときの[math]\lambda[/math]の値は，クロス・バリデーションによって得られたオブジェクトの変数lambda.minに格納されています．
#geshi(sh){{
> LASSO.cv$lambda.min
[1] 0.008431423
}}

係数は変数betaに，定数項の値は変数a0に格納されています．
しかし，各[math]\lambda[/math]ごとの係数が格納された行列やベクトルになっていますので，which関数を用いて交差検定の予測誤差が最も小さかったときの[math]\lambda[/math]のインデックスを調べ，そのときの係数と定数項の値だけを表示します．
#geshi(sh){{
> LASSO$beta[,which(LASSO$lambda == LASSO.cv$lambda.min)]
Sepal.Length  Sepal.Width  Petal.Width 
   0.7171234   -0.6267432    1.4504055 
> LASSO$a0[which(LASSO$lambda == LASSO.cv$lambda.min)]
       s57 
-0.2557475 
}}
重回帰と同様に表すと，学習された式を次のように書くことができます．
\[y = -0.2557475 + 0.7171234 x_1 - 0.6267432 x_ 2 + 1.4504055 x_3\]

単回帰，重回帰のときと同じように，平均二乗誤差を求め，予測誤差を視覚化します．
#geshi(sh){{
> summary((iris$Petal.Length - predict(LASSO, as.matrix(iris[,-5][,-3]), s=LASSO.cv$lambda.min))**2)
       1            
 Min.   :1.499e-05  
 1st Qu.:5.382e-03  
 Median :3.340e-02  
 Mean   :9.915e-02  
 3rd Qu.:1.349e-01  
 Max.   :1.167e+00  
> plot(iris$Petal.Length, predict(LASSO, as.matrix(iris[,-5][,-3]), s=LASSO.cv$lambda.min))
> abline(a=0,b=1)
}}
#ref(iris_lasso_errors.png,nolink,50%)




*Ridge回帰 [#wa18ca12]

Ridge回帰は，目的関数の正則化項が[math]\lambda||\beta||_2^2[/math]となっている回帰です．
つまり，Ridge回帰の目的関数は次のような形をしています．
\[\sum_{i=1}^N || y_i - (\beta_0 + \beta^\mathrm{T}\bf{x}_i)||^2 + \lambda ||\beta||_2^2\]

Ridge回帰の正則化をL2ノルムといい，次のように計算します．
\[||\beta||_2^2 = \sum_{i=1}^p \beta_i^2\]

LASSO回帰では正則化項のペナルティが係数の絶対値に比例して増加しましたが，Ridge回帰では係数の2乗に比例して増加します．


Ridge回帰を行うには，glmnet関数で学習する際にオプションでalpha=0を指定します．
#geshi(sh){{
> Ridge <- glmnet(as.matrix(iris[,-5][,-3]), iris$Petal.Length, alpha=0)
}}
後はLASSO回帰のときと同じです．

まず，正則化項の大きさと係数の関係を表示します．
#geshi(sh){{
> plot(Ridge)
}}
#ref(iris_ridge_regularization.png,nolink,50%)
LASSO回帰のグラフと比べると，LASSO回帰の方が係数を0にする傾向が強いことがわかります．

次に，クロス・バリデーションによって最適な[math]\lambda[/math]を求め，そのときの係数と定数項を表示します．
#geshi(sh){{
> Ridge.cv <- cv.glmnet(as.matrix(iris[,-5][,-3]), iris$Petal.Length, alpha=0)
> Ridge.cv$lambda.min
[1] 0.1694069
> Ridge$beta[,which(Ridge$lambda == Ridge.cv$lambda.min)]
Sepal.Length  Sepal.Width  Petal.Width 
   0.7882905   -0.6838069    1.2645594 
> Ridge$a0[which(Ridge$lambda == Ridge.cv$lambda.min)]
       s99 
-0.2742464 
}}
これまでと同様に表すと，学習された式は次のように書けます．
\[y = -0.2742464 + 0.7882905 x_1 - 0.6838069 x_2 + 1.2645594 x_3\]

最後に，平均二乗誤差を求め，予測誤差をプロットします．
#geshi(sh){{
> summary((iris$Petal.Length - predict(Ridge, as.matrix(iris[,-5][,-3]), s=Ridge.cv$lambda.min))**2)
       1            
 Min.   :1.296e-06  
 1st Qu.:9.084e-03  
 Median :4.576e-02  
 Mean   :1.083e-01  
 3rd Qu.:1.463e-01  
 Max.   :1.152e+00  
> plot(iris$Petal.Length, predict(Ridge, as.matrix(iris[,-5][,-3]), s=Ridge.cv$lambda.min))
> abline(a=0,b=1)
}}
#ref(iris_ridge_errors.png,nolink,50%)

ちなみに，LASSO回帰のときと同じように[math]\lambda[/math]と平均二乗誤差の関係をプロットしようとしましたが，できませんでした．
#geshi(sh){{
> plot(Ridge.cv)
 以下にエラー axis(side = 3, at = sign.lambda * log(cvobj$lambda), labels = paste(cvobj$nz),  : 
   'at' と 'label' の長さが違います。 100 != 300 
}}







*応用 [#cac50579]

diabetes（糖尿病）データに応用してみましょう．

このデータはlarsパッケージ((EFRON, HASTIE,  JOHNSTONE, and TIBSHIRANI: Least Angle Regression, The Annals of Statistics, Vol. 32, No. 2, pp. 407–499, 2004.))に含まれていますので，まずはlarsパッケージをインストールします．
#geshi(sh){{
> install.packages("lars")
> library(lars)
> data(diabetes)
}}
diabetesデータには，442個の訓練データが含まれており，10次元の説明変数を表す行列xと目的変数を表すベクトルyが格納されています．

まず，重回帰を行い，平均二乗誤差を求めます．
#geshi(sh){{
> Multi <- lm(diabetes$y ~ diabetes$x)
> summary((diabetes$y - predict(Multi))**2)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
    0.041   331.200  1484.000  2860.000  3928.000 24280.000 
}}


次に，LASSO回帰を行い，正則化項と係数の関係をプロットします．
#geshi(sh){{
> LASSO <- glmnet(diabetes$x, diabetes$y)
> plot(LASSO)
}}
#ref(diabetes_lasso_regularization.png,nolink,50%)
正則化が弱いときは（つまり，グラフの右の方では），irisに比べて係数がかなり大きいことがわかります．

クロス・バリデーションでLASSO回帰の最適な[math]\lambda[/math]を求めます．
#geshi(sh){{
> LASSO.cv <- cv.glmnet(diabetes$x, diabetes$y)
> LASSO.cv$lambda.min
[1] 0.0735996
> plot(LASSO.cv)
}}
#ref(diabetes_lasso_lambda.png,nolink,50%)

LASSO回帰の平均二乗誤差を求めます．
#geshi(sh){{
> summary((diabetes$y - predict(LASSO, diabetes$x, s=LASSO.cv$lambda.min))**2)
       1            
 Min.   :1.924e-05  
 1st Qu.:3.149e+02  
 Median :1.504e+03  
 Mean   :2.862e+03  
 3rd Qu.:3.853e+03  
 Max.   :2.392e+04  
}}

最後に，Ridge回帰を行います．
#geshi(sh){{
> Ridge <- glmnet(diabetes$x, diabetes$y, alpha=0)
> plot(Ridge)
}}
#ref(diabetes_ridge_regularization.png,nolink,50%)

クロス・バリデーションによってRidge回帰の最適な[math]\lambda[/math]を求めます．
#geshi(sh){{
> Ridge.cv <- cv.glmnet(diabetes$x, diabetes$y, alpha=0)
> Ridge.cv$lambda.min
[1] 4.516003
}}

Ridge回帰の平均二乗誤差を求めます．
#geshi(sh){{
> summary((diabetes$y - predict(Ridge, diabetes$x, s=Ridge.cv$lambda.min))**2)
       1            
 Min.   :1.750e-02  
 1st Qu.:3.142e+02  
 Median :1.514e+03  
 Mean   :2.881e+03  
 3rd Qu.:3.872e+03  
 Max.   :2.353e+04  
}}


*まとめ [#o5e18fce]
LASSO回帰やRidge回帰は，誤差関数に正則化項を加えることによってモデルに制約を加えるものです．




*参考文献 [#i630eb1c]

7.17節「LASSO」を理解するためにこのページを作成しました．
#html{{
<iframe src="http://rcm-jp.amazon.co.jp/e/cm?t=tohgoroh-22&o=9&p=8&l=as1&asins=4320057082&ref=tf_til&fc1=444B4C&IS2=1&lt1=_blank&m=amazon&lc1=444B4C&bc1=FFFFFF&bg1=FFFFFF&f=ifr" style="width:120px;height:240px;" scrolling="no" marginwidth="0" marginheight="0" frameborder="0"></iframe>
}}

LASSO回帰とRidge回帰については，この本に詳しく載っています．
#html{{
<iframe src="http://rcm-jp.amazon.co.jp/e/cm?lt1=_blank&bc1=FFFFFF&IS2=1&bg1=FFFFFF&fc1=000000&lc1=444B4C&t=tohgoroh-22&o=9&p=8&l=as4&m=amazon&f=ifr&ref=ss_til&asins=0387848576" style="width:120px;height:240px;" scrolling="no" marginwidth="0" marginheight="0" frameborder="0"></iframe>
}}




*注 [#a434ba10]

この記事はまだ書きかけです．