医療統計学を学ぶ大学生のブログ

医療統計学、因果推論を専攻しています。SASユーザーです。

Lasso回帰モデルにおけるSolution pathについて

f:id:NorihiroSuzuki:20210910180047p:plain


この記事は立ち位置としては、以前のこの記事の続きとなります。

norihirosuzuki.hatenablog.com

 

Lassoの概要

Lassoではまず初めに以下のような線形モデルを考えます。

f:id:NorihiroSuzuki:20210806135106p:plain

 ここでyiは目的変数、xijは説明変数、解析で推定したいパラメータは係数βjとし、iはサンプルのインデックスであると仮定します。なお誤差項εは平均0、分散σ2正規分布に独立に従うとします。

 

通常の回帰分析では、以下の式で表される残差平方和(RSS; residual sum of squares)を最小にするパラメータβjを推定します。

f:id:NorihiroSuzuki:20210806140351p:plain

 

 Tibshiraniが1996年に提案したLassoの問題は、上記のRSSの式に以下のように条件を付けたしたものです。

f:id:NorihiroSuzuki:20210806142602p:plain

この制約条件(subject to ~)はβの凸領域を与えます。n=p=2の時のこの領域を表したものがよく見るこの図です。

f:id:NorihiroSuzuki:20210806133624p:plain

 楕円の中心がy、楕円が最適化したい関数(RSS)に対する等高線、広がりを示しています。凸領域があることによって、この領域と関数が接するのが自然と軸上になり、つまり、x1=0と疎な解を与えることが分かります。

 この最適化を与えるβはtの値に影響され、tが小さければ多くのβ(パラメータ)が0に、tが大きければ0にならないパラメータが増えます。

 

ただこのTibshiraniによる問題よりもラグランジュの未定乗数法によって制約を設けず、罰則項付きの式として表されたほうがより一般的です。このどちらの表現でも問題としては等価であり、同じ解を与えるtとλは1対1に対応します。λの値が小さくなるほど0になるパラメータの値が増えます。

f:id:NorihiroSuzuki:20210910173809p:plain

 

次にλをどう選ぶかということですが、一般的な方法の一つは、k-Fold Cross Validation(k分割交差検証)です。これはデータをk個に分け1つを教師データに、k-1個を訓練データとして、最適なλの値を探索していくという方法です。

変数選択の観点からは、このλを適切に選ぶということが重要ですが、サンプルサイズnよりもパラメータpが大きい場合(n<p)の時には、やはり難しくその判断に注意が必要です。

 

この対応としては、Bootstarap法をLassoと組み合わせるStability selectionというものやSolution pathというものがあります。Stability Selectionについてはまた別にまとめるので、この記事ではSolution pathについてまとめます。

 

Solution pathについて

Solution pathを見るモチベーションとしては、正則化パラメータであるλを一組にすることで誤った結論を導いてしまうということです。Solution pathは簡単に言うと最適解として得られたλの前後で、パラメータβjの推定値がどのような挙動を見せるのかを図示したものです。

このSolution pathの作成に用いられる一般的なアルゴリズムにはLARS(Least Angle Regression Selection)や、CDA(Coordinate Descent Algorithm)があります。この記事中ではLARSについてまとめています。

 

LARSについて

LARSでは、あるλに対してその前後でλが動いても、そのλで得られる最適なパラメータのうち推定値が0でないβjの組が変わらないという仮定をまず初めに置きます。つまり、推定値が0であるパラメータはそのまま0、0でないパラメータは0でない(変化はOK)としているわけです。

この時、パラメータは値を持っている成分と、持っていない(0)の成分に分けられますが、この値を持つ成分をactive setとします。そしてこれに対応する列ベクトルをXA、パラメータをβA(λ)とします。

ここで、

f:id:NorihiroSuzuki:20210910173901p:plain

をベクトルで表記したものを微分することにより、以下の関係式が算出されます。

f:id:NorihiroSuzuki:20210910174622p:plain

XAの列数はN以下であることより、逆行列が存在するため、βA(λ)について解くと

f:id:NorihiroSuzuki:20210910175305p:plain

となりλの1次式であることから、0でないパラメータについては直線で表すことが可能ということが分かります。ちなみにこの1次式であるということが計算スピードが高速な理由です。

またこの式が成り立つのはあくまでも、あるλに対してその前後でλが動いても、そのλで得られる最適なパラメータのうち推定値が0でないβjの組が変わらないという仮定のもとであるため、この仮定が満たされない(=パラメータを持つ組が変わる)時にはそのλを区切りとしてまた、その部分で同様に上記の問題を解くことでSolition pathを作成することができます。

 

SASHELPにあるbaseballのデータを使って、Lasso回帰、Solution pathの作成を行ってみます。用いるプロシジャはGLMSELECTプロシジャで、Solution pathを出すには、plot=allか、plot=coefficientpanelとすればOKです。

baseballデータセットの中身はこんな感じ。

f:id:NorihiroSuzuki:20210910180351p:plain

 

ここで、

目的変数

  • logSalary

説明変数

  • nAtBat
  • nHits
  • nHome
  • nRuns
  • nRBI
  • nBB
  • yrMajor
  • crAtBat
  • crHits
  • crHome
  • crRuns
  • crRbi               
  • crBB
  • league
  • division
  • nOuts
  • nAssts
  • nError

とし、Lasso回帰を行った場合の選択結果とSolution pathはこんな感じ。

 

f:id:NorihiroSuzuki:20210910180714p:plain

f:id:NorihiroSuzuki:20210910180649p:plain

f:id:NorihiroSuzuki:20210910180739p:plain

結果に対する考察は行いませんが、無事にエラーなく変数選択、Solution pathの作成が行われました。

 

今回使用したSASコードは以下の通り

proc glmselect data=sashelp.baseball plots=coefficientpanel;
   class league division;
   model logSalary = nAtBat nHits nHome nRuns nRBI nBB
                  yrMajor crAtBat crHits crHome crRuns crRbi
                  crBB league division nOuts nAssts nError
                / selection=stepwise(select=SL choose=PRESS);
run;

 

参考

Lassoをはじめとする基本的なスパース推定について↓