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

はじめに †

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

きちんと理解するために，まずは単回帰と重回帰という2種類の線形回帰をやり，その後にlasso回帰とRidge回帰をやります．その後，多項回帰と二項回帰をやります．

準備 †

Rのインストールについては，次のページを見てください．

最初は，標準で使用できるirisデータセットを使います．

data(iris)

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

このデータセットには，150個のデータが含まれています．ランダムに10個のデータを選択して，見てみましょう．

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

↑

線形回帰（単回帰） †

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

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

例として，がくの長さ（Petal.Length）を横軸，花びらの長さ（Sepal.Length）を縦軸にとって，150個のデータをプロットしてみます．

    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

どうやら，花びらの長さとがくの長さの間には正の相関（一方が大きいと他方も大きいという関係）がありそうですね．そこで，花びらの長さ（Sepal.Length）を [math]x_1[/math]（説明変数），がくの長さ（Petal.Length）を [math]y[/math]（目的変数）として，単回帰を行いましょう．

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

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

線形回帰は，標準で用意されているlm関数にモデル式を入力して行います．

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

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

iris.sgl <- lm(y ~ x1)

Coefficientsは係数を表し，Interceptは定数項の値（切片）を表します．つまり，花びらの長さを [math]x_1[/math]，がくの長さを [math]y[/math] とすると，学習された式は次のように書けます． \[ y = -7.101 + 1.858 x_1 \]

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

iris.sgl

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

Call:
lm(formula = y ~ x1)

Coefficients:
(Intercept)           x1  
     -7.101        1.858

この予測がどのくらい良いのかを表すために，平均二乗誤差（実際の値と予測値の差の二乗の平均）を求めます． \[\textit{MSE} = \frac{1}{N}\sum_{i=1}^{N} \left(y_i - (\beta_0 + \beta_1 x_{i,1})\right)^2\]

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

↑

線形回帰（重回帰） †

説明変数が複数ある線形回帰を重回帰（multiple linear regression）といいます．

重回帰のモデルは次のような式で表されます． \[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]を要素とするベクトルです．

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

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

summary(iris.sgl)

説明変数が複数あるときは，モデルを表すときにそれぞれの変数を記号+で連結して表します．ここでは，「y ~ x1+x2+x3」となります．

単回帰と同じように，lm関数を用いてモデルを学習します．

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

単回帰と同じように，学習された式は次のように書けます． \[ y = -0.2627112 + 0.7291384 x_1 - 0.6460124 x_2 + 1.4467934 x_3 \]

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

plot(y, predict(iris.sgl))
abline(a=0,b=1)

単回帰と同様に，平均二乗誤差を求めます． \[\textit{MSE} = \frac{1}{N} \sum_{i=1}^N (y_i - (\beta_0 + \beta^\mathrm{T}\bf{x}_i))^2\]

mean((y - predict(iris.sgl))**2)

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

↑

lasso回帰 †

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

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

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

[1] 0.743061

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

y <- iris$Petal.Length
x1 <- iris$Sepal.Length
x2 <- iris$Sepal.Width
x3 <- iris$Petal.Width

cbind関数は、列 (column) で結合するもので、x1, x2, x3 のそれぞれを列とした行列を作っています。

plot関数を使うと，正則化項の大きさとそれぞれの係数の関係をプロットすることができます．

iris.mult <- lm(y ~ x1+x2+x3)
iris.mult$coefficients

黒色の折れ線が花びらの長さ（Sepal.Length），赤色の折れ線が花びらの幅（Sepal.Width），緑色の折れ線ががくの幅（Petal.Width）を表しています．

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

最も正則化が最も強いときは，すべての係数が0になっています．そこから正則化を弱めていくと，緑色の線が表すがくの幅（Petal.Width）の係数だけが大きくなっていきます．その後，他の係数も0ではなくなり，正則化項が0，つまりペナルティなしのときに重回帰のときと同じ係数になっています．

次に，最適な正則化パラメーター（[math]\lambda[/math]）の値を交差検定（クロス・バリデーション）を用いて求めます．

(Intercept)          x1          x2          x3 
 -0.2627112   0.7291384  -0.6460124   1.4467934

plotを用いて，[math]\lambda[/math]と平均二乗誤差の関係を表示すると次のようになります．

summary(iris.mult)

今回の例では正則化パラメーターの値が小さいほうが平均二乗誤差が小さくなりましたが，いつもこうなるわけではありません．

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

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

係数は変数betaに，定数項の値は変数a0に格納されています．しかし，各[math]\lambda[/math]ごとの係数が格納された行列やベクトルになっていますので，which関数を用いて交差検定の予測誤差が最も小さかったときの[math]\lambda[/math]のインデックスを調べ，そのときの係数と定数項の値だけを表示します．

plot(y, predict(iris.mult))
abline(a=0,b=1)

重回帰と同様に表すと，学習された式を次のように書くことができます． \[y = -0.2557475 + 0.7171234 x_1 - 0.6267432 x_ 2 + 1.4504055 x_3\]

単回帰，重回帰のときと同じように，平均二乗誤差を求め，予測誤差を視覚化します．

mean((y - predict(iris.mult))**2)

↑

Ridge回帰 †

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

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

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

Ridge回帰を行うには，glmnet関数で学習する際にオプションでalpha=0を指定します．

[1] 0.09901965

後はlasso回帰のときと同じです．

まず，正則化項の大きさと係数の関係を表示します．

install.packages("glmnet")
library(glmnet)

lasso回帰のグラフと比べると，lasso回帰の方が係数を0にする傾向が強いことがわかります．

次に，クロス・バリデーションによって最適な[math]\lambda[/math]を求め，そのときの係数と定数項を表示します．

x <- cbind(x1, x2, x3)
iris.lasso <- glmnet(x, y)

これまでと同様に表すと，学習された式は次のように書けます． \[y = -0.2742464 + 0.7882905 x_1 - 0.6838069 x_2 + 1.2645594 x_3\]

最後に，平均二乗誤差を求め，予測誤差をプロットします．

plot(iris.lasso)

ちなみに，lasso回帰のときと同じように[math]\lambda[/math]と平均二乗誤差の関係をプロットしようとしましたが，できませんでした．

iris.lasso.cv = cv.glmnet(x, y)

↑

多項回帰 †

irisデータのように目的変数（Species）の値が離散値（カテゴリー）のときは，多項回帰（Multinomial regression）を用います．

多項回帰を行うには，glmnet関数で学習する際にオプションでfamily="multinomial"を指定します．

plot(iris.lasso.cv)

それから，クロス・バリデーションを用いて最適な[math]\lambda[/math]を求めます．

iris.lasso.cv$lambda.min

今回は[math]\lambda[/math]が小さすぎても大きすぎてもダメです．

最後に，予測結果をtable関数を用いて表にします．目的変数が離散値なので，predict関数にオプションtype="class"を指定します．

[1] 0.008431423

setosaは50個全て正しく分類できていますが，versicolorのうち2つをvirginicaと誤って分類し，virginicaの1つをversicolorと誤って分類しています．

↑

二項回帰 †

目的変数の値が「はい」「いいえ」や「陽性」「陰性」など2種類の離散値（カテゴリー）のときは，二項回帰（binomial regression）を用います．

irisデータには，目的変数（Species）の値がsetosa，versicolor，virginicaであるデータが50個ずつ含まれていますから，ここからsetosaとversicolorだけを取り出して新しいデータ・フレームを作成し，二項回帰をやってみます．

まず，subset関数を用いてirisからSpeciesの値がsetosaかversicolorであるデータだけを取り出します．また，factor関数を用いてSpciesのレベルを再定義します．

iris.lasso$beta[,which(iris.lasso$lambda==iris.lasso.cv$lambda.min)]

二項回帰を行うには，glmnet関数で学習する際にオプションでfamily="binomial"を指定します．後は多項回帰と同じです．

Sepal.Length  Sepal.Width  Petal.Width 
   0.7171234   -0.6267432    1.4504055

iris.lasso$a0[which(iris.lasso$lambda==iris.lasso.cv$lambda.min)]

       s57 
-0.2557475

こちらは全て正しく分類できました．

↑

糖尿病データへの応用 †

それでは，diabetesデータに応用してみましょう．

このデータはlarsパッケージ*2に含まれていますので，まずはlarsパッケージをインストールします．

mean((y - predict(iris.lasso, x, s=iris.lasso.cv$lambda.min))**2)

diabetesデータには，442人の糖尿病患者から観測された年齢，性別，BMI，平均血圧，6種類の血清測定値から成る10次元の説明変数行列xと1年後の糖尿病の進行度を表す目的変数ベクトルyが格納されています．

まず，重回帰を行い，平均二乗誤差を求めます．

[1] 0.09915377

次に，lasso回帰を行い，正則化項と係数の関係をプロットします．

plot(y, predict(iris.lasso, x, s=iris.lasso.cv$lambda.min))
abline(a=0,b=1)

正則化が弱いときは（つまり，グラフの右の方では），irisに比べて係数がかなり大きいことがわかります．

クロス・バリデーションでlasso回帰の最適な[math]\lambda[/math]を求めます．

iris.ridge <- glmnet(x, y, alpha=0)

lasso回帰の平均二乗誤差を求めます．

plot(iris.ridge)

最後に，Ridge回帰を行います．

iris.ridge.cv <- cv.glmnet(x, y, alpha=0)
iris.ridge.cv$lambda.min

クロス・バリデーションによってRidge回帰の最適な[math]\lambda[/math]を求めます．

[1] 0.1694069

Ridge回帰の平均二乗誤差を求めます．

iris.ridge$beta[,which(iris.ridge$lambda==iris.ridge.cv$lambda.min)]

↑

構造活性相関データへの応用 †

いよいよ，『Rによるバイオインフォマティクスデータ解析』の7.17節で用いられている構造活性相関のデータであるdhfrデータセットを対象に回帰を行います．

このデータセットは，caretパッケージ*3に含まれていますので，インストールしてから使用します．

Sepal.Length  Sepal.Width  Petal.Width 
   0.7882905   -0.6838069    1.2645594

dhfrデータセットの中を調べてみます．

iris.ridge$a0[which(iris.ridge$lambda==iris.ridge.cv$lambda.min)]

事例数が325，属性数が229です． 1個目の属性Yが目的変数で，あとの228個はすべて説明変数です．これまでの例に比べると，説明変数がすごく多いですね． Y，つまり目的変数の値はactiveとinactiveの2種類です．

目的変数の値が2種類のカテゴリーなので，二項回帰を行います．

       s99 
-0.2742464

説明変数の数が多いので，たくさんの線が描かれています．意味はこれまでと同じで，一番左が最も大きく正則化して係数を小さくしたとき，一番右が正則化をまったく行っていないときです．

次に，クロス・バリデーションで最適な[math]\lambda[/math]を求めます．

mean((y - predict(iris.ridge, x, s=iris.ridge.cv$lambda.min))**2)

最後に，精度を求めます．

[1] 0.1083066

96.6%の正解率でした．

↑

まとめ †

線形回帰のモデルは次のような式で表されます． \[y = \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] を要素とするベクトルです．次元数が [math]p=1[/math] のときを単回帰，[math]p>1[/math] のときを重回帰といいます．

線形回帰では，二乗誤差の総和 [math]L(\beta, \bf{x})[/math] を最小化するような [math]\beta[/math] を求めます． \[L(\beta, \bf{x}) = \sum_{i=1}^N \left(y_i - (\beta_0 + \beta^\mathrm{T}\bf{x}_i)\right)^2\]

lasso回帰やRidge回帰では，損失関数 [math]L(\beta, \bf{x})[/math] に正則化項 [math]\lambda P(\beta)[/math] を加えたもの最小化するような [math]\beta[/math] を求めます． \[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]（次元数が1）のとき，正則化なしの単回帰になります． glmnet関数では，[math]\lambda[/math] の値をオプションsに設定します．

glmnet関数の正則化項 [math]P(\beta)[/math] は，次のような一般的な形をしています． \[P(\beta) = (1 - \alpha)\frac{1}{2} ||\beta||_2^2 + \alpha ||\beta||_1\]

[math]\alpha[/math] は正則化の形を決めるパラメーターです． [math]\alpha=1[/math] のときにlasso回帰，[math]\alpha=0[/math] のときにRidge回帰になります． glmnet関数では，この [math]\alpha[/math] の値をオプションalphaに設定します．デフォルトではalpha=1なので，alphaオプションを指定しないときはlasso回帰になります． alpha=0と指定するとRidge回帰になります． alphaの値を0と1の間に設定することもできます．

今回は，線形回帰の誤差関数に正則化項を加えたlasso回帰とRidge回帰をやりましたが，ロジスティック回帰，多項式回帰，ポアソン回帰，Cox回帰の誤差関数に正則化項を加えたlasso回帰やRidge回帰もできます．また，多項回帰と二項回帰では，特に説明せずにデフォルトのlasso回帰を用いましたが，Ridge回帰を用いた多項回帰や二項回帰もできます．

回帰分析についてもっと詳しく知りたい場合は，参考文献に上げた『The Elements of Statistical Learning』を読むといいでしょう． glmnet関数でロジスティック回帰，多項式回帰，ポアソン回帰，Cox回帰を使いたい場合は，glmnetパッケージのマニュアルを読みましょう．

↑

参考文献 †

この本の7.17節「LASSO」を参考にしています．

lasso回帰とRidge回帰については，この本に詳しく載っています．この本は統計学習の教科書で，著者はglmnetパッケージの開発者です．

glmnetパッケージのソース，バイナリー，マニュアルがCRANからダウンロードできます．

glmnet: Lasso and elastic-net regularized generalized linear models | CRAN

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