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

追加された行はこの色です。
削除された行はこの色です。
*はじめに [#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>
}}

LASSO回帰の前に，線形回帰と最小二乗法をやります．
LASSO回帰の前に，線形回帰をやります．
また，Ridge回帰もやります．


*準備 [#e708ecca]

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

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

150個のデータが含まれています．
最初の10個を表示すると，次のようになります．
#geshi(sh){{
> data(iris)
> names(iris)
[1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"     
> 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つの線形回帰（単回帰）を行います．

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

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

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

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

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

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

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

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


*最小二乗法 [#aac422ee]


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

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

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

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

まずは，線形モデルによる回帰を行います．
単回帰と同じように，lm関数を用いて線形回帰を行います．
#geshi(sh){{
> Linear <- lm(iris$Petal.Length ~ iris$Sepal.Length+iris$Sepal.Width+iris$Petal.Width)
> Linear$coefficients
      (Intercept) iris$Sepal.Length  iris$Sepal.Width  iris$Petal.Width 
       -0.2627112         0.7291384        -0.6460124         1.4467934 
}}

重回帰分析の場合はモデルが視覚化しにくいので，予測誤差だけ視覚化します．
> plot(iris$Petal.Length, predict())
> abline(a=0,b=1)
#ref(iris_multiple_errors.png,nolink,50%)
単回帰の場合よりも誤差が小さくなっていることがわかります．

ここでは，線形回帰，LASSO回帰，Ridge回帰を行います．


*LASSO回帰 [#p3d71085]

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関数をgolmnetパッケージのhead関数とtail関数で置き換えましたよ」という意味なので，気にしなくていいです．

データは，標準で使えるiris，glmnet
LASSO回帰は，glmnetパッケージのglmnet関数を用いて行います．
引数には，説明変数行列，目的変数，そしてオプションを与えます．
#geshi(sh){{
> LASSO <- glmnet(matrix(iris$Sepal.Length), iris$Petal.Length, alpha=1)
}}
ここで，


この記事はまだ書きかけです．
バイオ・データ・マイニング/Rで回帰分析する のバックアップ差分(No.2)

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