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

はじめに †

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

LASSO回帰の前に，線形回帰をやります．また，Ridge回帰もやります．

準備 †

まずは，標準で使用できるirisデータを使います．

data(iris)

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

150個のデータが含まれています．最初の10個を表示すると，次のようになります．

iris[sort(sample(1:150,10)),]

ランダムにサンプリングした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

↑

線形回帰（単回帰） †

まず，説明変数が1つの線形回帰（単回帰）を行います．

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

この関係をプロットすると次のようになります．

plot(iris$Sepal.Length, iris$Petal.Length)

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

線形回帰は，標準で用意されているlm関数で行います．

y <- iris$Petal.Length
x1 <- iris$Sepal.Length

学習された結果を表示してみると，次のようになります．

iris.sgl <- lm(y ~ x1)

Coefficientsは係数を表し，Interceptは定数項（[mathy[/math]切片）の値を表します．つまり，花びらの長さを[math]x[/math]，がくの長さを[math]y[/math]とすると，学習された式は[math]y = 1.858433 x - 7.101443[/math]です．

学習された線形回帰モデルを図示します．

iris.sgl

すると，次のようなグラフがプロットされます．

単回帰だとモデルが視覚化しやすいです．

次に，予測誤差を視覚化するため，それぞれのデータについて，横軸に実際のがくの長さ（iris$Petal.Length），縦軸にその予測値をとってグラフを表示します．予測値はpredict関数で求めることができます．

Call:
lm(formula = y ~ x1)

Coefficients:
(Intercept)           x1  
     -7.101        1.858

↑

線形回帰（重回帰） †

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

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

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

単回帰と同じように，lm関数を用いて線形回帰を行います．

plot(x1, y)
abline(iris.sgl)

重回帰分析の場合はモデルが視覚化しにくいので，予測誤差だけ視覚化します．

plot(iris$Petal.Length, predict())

abline(a=0,b=1)

単回帰の場合よりも誤差が小さくなっていることがわかります．

↑

LASSO回帰 †

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

summary(iris.sgl)

私の環境ではこのような警告が出ますが，これは「utilsパッケージのhead関数とtail関数をgolmnetパッケージのhead関数とtail関数で置き換えましたよ」という意味なので，気にしなくていいです．

LASSO回帰は，glmnetパッケージのglmnet関数を用いて行います．引数には，説明変数行列，目的変数，そしてオプションを与えます．

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 ‘ ’ 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

ここで，