ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

罰則付き回帰モデル（LASSO回帰, Ridge回帰, Elastic Net）で過学習を抑える（glmnetパッケージ）[R]

variable selection prediction/diagnostic validation R neco

回帰モデルの説明変数は多いほど良いという訳ではない。
15組の(x,y)を2次式で回帰した場合（赤・点線）と10次式で回帰した場合（青・実線）を比べると、後者の方が15のデータからのズレが小さいが、新しいデータには前者の方が上手く当てはまる。

図1：2次式（赤・点線）と10次式（青・実線）を当てはめた

このようにデータに過度に適合しすぎて返って使い物にならなくなることを、過剰適合（overfitting）とか過学習という。

罰則付き回帰モデル（penalized regression model）で過学習を抑える

xが少し変化しただけでyが大きく変化するモデル、つまり、xの係数βの絶対値が大きすぎるモデルは過学習を起こしている可能性があるので、|β|の大きさに応じてペナルティーを与えることで適正なモデルを探そうという話。 正則化回帰モデル（regularized regression model）ともいう。

下のように、「モデルからのズレ（=損失関数）」と「ペナルティー（=正則化項）」の和（=目的関数）を小さくするモデルを選ぶ。

glmnet( )のヘルプには

The objective function for "gaussian" is
　　1/2 RSS/nobs + λ*penalty,
and for the other models it is
　　-loglik/nobs + λ*penalty.

とあるので、正規分布モデルでは残差平方和（residual sum of squares, RSS）、その他のモデルでは対数尤度を「ズレの指標」に使っているみたい。

LASSO回帰, Ridge回帰, Elastic Netは正則化項（ペナルティー）が違う

ラッソ回帰（LASSO*1 regression）ではパラメータの絶対値の和 \Sigma |\beta| が、リッジ回帰（Ridge regression）ではパラメータの2乗和 \Sigma \beta ^2 が正則化項として用いられる *2。 Elastic Netは下のように両者を混合したもの（凸結合*3したもの）。

\displaystyle \begin{align} (penalty) = \sum_{j=1}^p \left\{ \alpha | \beta_j | + (1-\alpha) \beta_j^2 \right\}, 0\leq \alpha \leq 1 \end{align}

α=1のときはLASSO回帰、α=0のときはRidge回帰になる。
Lasso回帰では影響の小さい変数が除かれていく（=係数βが0になる）ので変数選択が可能だが、Ridge回帰では基本的に全ての変数が残るので変数選択には使えない。

glmnet( )でペナルティーの重みλを色々変えて係数βを推定する

glmnetパッケージのQuickStartExampleデータを使う。

library(glmnet)
data(QuickStartExample)
x <- QuickStartExample$x
y <- QuickStartExample$y 

> str(QuickStartExample)
List of 2
 $ x: num [1:100, 1:20] 0.274 2.245 -0.125 -0.544 -1.459 ...
 $ y: num [1:100, 1] -1.275 1.843 0.459 0.564 1.873 ...

従属変数yは連続変数、説明変数xは20個の連続変数・100観察分のmatrixになっている。

説明変数xと従属変数y以外に、glmnet( )の引数で使いそうなものは次のとおり。

  • family：モデルの種類。
    • "gaussian"：正規分布モデル（デフォルト）
    • "binomial"：二項分布モデル
    • "poisson"：ポアソン分布
    • "multinomial"：多項分布モデル
    • "cox"：Coxモデル
    • "mgaussian"：多変量正規分布モデル（アウトカムが複数の連続変数で構成される）
  • alpha：0〜1の値を指定する。1ならLasso回帰、0ならRidge回帰になる。デフォルトは1。
  • nlambda：計算するλの個数。デフォルトは100。
  • standardize：モデルを当てはめる前にxを標準化するかどうか*4。デフォルトはTRUE。
  • lower.limits, upper.limits：βの下限・上限。デフォルトは-Inf〜Inf。
  • penalty.factor：変数毎に罰則の大きさを変えたい時に指定する。
fit <- glmnet(x,y)

結果fitをprint( )に渡すと、試されたλ毎の結果が表示される。

> print(fit)

Call:  glmnet(x = x, y = y) 

   Df  %Dev  Lambda
1   0  0.00 1.63100
2   2  5.53 1.48600
3   2 14.59 1.35400
4   2 22.11 1.23400
5   2 28.36 1.12400
~~~
  • Lambda：重み付けパラメータλ
  • DF：係数が0になっていない説明変数の個数
  • %Dev：モデルで説明されるデビアンスの割合

λが1.631より大きくなると全てのβが0になるのでモデルで説明されるdevianceが0になったということ。

結果fitをplot( )に渡すと、λを変化させたときの各係数の推移をプロットする。
指定する引数は、次のとおり。

  • xvar：横軸に何をとるか。
    • "norm"とするとL1正則（=罰則の大きさ）
    • "lambda"とするとλの対数
    • "dev"とするとモデルで説明されているdevianceの割合
  • label：=TRUEとすると変数名も表示する

図2：対数λに対して係数βの推移をプロットした

λが小さいとき（グラフ左）は罰則が相対的に小さいので全ての変数が残っているが、λを大きくしていく（グラフ右）と影響の小さい変数の係数が0に落ちていく。

図3：L1正則に対して係数βの推移をプロットした

変数がたくさん残っている（グラフ右）うちはL1正則は大きいが、λを大きくしていくにつれて変数が脱落していき（グラフ左）L1正則は小さくなっていく。

cv.glmnet( )で重み付けパラメータλを変えながらモデルの交差検証を行う

交差検証（cross-validation, CV）はモデルの汎化性能（= 新しく測定されるデータをどれくらい誤差なく予測できるか）を評価する方法。 K分割交差検証（K-fold cross-validation）ではデータをK個に分割し、そのうちのK-1個のデータで作ったモデルを残りの1個のデータで評価する。

cv.glmnet( )で次の引数を指定する。

  • family：glmnetと同じ。
  • nfold：分けるグループ数。
  • foldid：グループの分け方を揃えたいときに指定する。
  • type.measure：交差検証に使う損失関数。指定できる種類・デフォルトは、familyで指定するモデルの種類によって決まる（下の表を参照。■はdefault）。
    • "default"：モデルのfamilyに応じたデフォルトの損失関数を使用する。
    • "mse"：平均2乗誤差（mean square error）
    • "deviance"：デビアンス
    • "class"：誤分類率（misclassification error）
    • "auc"：ROC曲線の曲線下面積（area under the curve）
    • "mae"：平均絶対誤差（mean absolute error）
    • "C"：Harrel's C

図4：モデルの種類別利用可能な損失関数

結果の渡し方はglmnet( )と同じ。

> cvfit <- cv.glmnet(x, y)
> print(cvfit)

Call:  cv.glmnet(x = x, y = y) 

Measure: Mean-Squared Error 

     Lambda Index Measure     SE Nonzero
min 0.06897    35   1.129 0.1276      10
1se 0.14517    27   1.247 0.1676       8

minでは損失関数が最小となるとなるときのλなどが表示される。
1seで表示されるλは何かというと、

largest value of lambda such that error is within 1 standard error of the minimum.

とのこと。CV errorの最小値に、そのときのCV errorの1標準誤差分を足したCV errorを与えるλを使うことを推奨する「one standard error rule」なるものがあるらしい*5。

任意のλのときの係数βを取り出したいときは、coef( )を使う。指定したいλは（なぜか）s=で指定する*6。
CV errorを最小にするλを与えたいときは"lambda.min"、one standar error ruleに基づいたλを与えたいときは"lambda.1se"を指定する（デフォルトはこちら）。

> coef(cvfit, s="lambda.min")
21 x 1 sparse Matrix of class "dgCMatrix"
                      s1
(Intercept)  0.147927056
V1           1.337393911
V2           .          
V3           0.704086480
~~~~~  
V18          .          
V19          .          
V20         -1.061382443

ちなみにglmnet()オブジェクトを使って次のようにしても同じ。

fit <- glmnet(x, y)
cvfit <- cv.glmnet(x, y)
coef(fit, s=cvfit$lambda.min

最後に交差検証のプロットを描く。

plot(cvfit)

図5：対数λに対して交差検証の結果をプロットした

プロットで表示される縦線は、lambda.minとlambda.1seに対応している。

おわりに

  • 目にする機会が増えたので勉強してみましたが、自分で使う機会があるかどうか...
  • 換毛期です。茶トラは明るい毛色なので大変です。

ランキング参加中
統計

参考資料

*1:least absolute shrinkage and selection operator

*2:正則化項について分かりやすく解説されています。qiita.com

*3:convex combination：和が1となる非負係数の線型結合

*4:標準化の詳細はこちら：A deep dive into glmnet:&nbsp;standardizestatisticaloddsandends.wordpress.com

*5:1SE ruleに関する議論：stats.stackexchange.com

*6:理由は"In case we want to allow one to specify the model size in other ways in the future."だからとのこと