新しい「LASSO」たち（2）
状況によって回帰係数を変えられる「Pliable Lasso」
もう一つ、新しいLassoのアルゴリズムを紹介しましょう。
Lassoを代表とするような回帰分析においては、
\( \hat{y} = \sum_{j = 1}^p X_j \hat{\beta}_j \)
というように、各説明変数\(X_j\)に対する回帰係数\(\hat{\beta}_j\)を推定していました。
一方、現実の問題では、この回帰係数を状況によって変えたいという要望が存在します。例えば以下のような場合です。
「ある薬は、男性にはよく効くが、女性にはあまり効かない」
数学的に言い換えれば、この薬を表す説明変数を\(X_j\)とすると、\(X_j\)に対する回帰係数\(\hat{\beta}_j\)を、男性に投与するときは大きく、女性に投与するときは小さく推定したいということになります。
これを可能にするアルゴリズムが、Lassoの開発者であるRobert Tibshirani教授によって開発された「Pliable Lasso」です[2]。
注
[2] Tibshirani, Robert, and Jerome Friedman. "A pliable lasso." Journal of Computational and Graphical Statistics 29.1 (2020): 215-225.
Pliable Lassoでは、説明変数\(X\)の他に、修正変数\(Z\)という変数を導入します。\(Z\)は先程の例でいう性別などを表す変数です。
そして、回帰モデルを
\( \hat{y} = \sum_{j = 1}^p X_j \circ (\hat{\beta}_j \mathbf{1} + Z \hat{\theta}_j) \)
と考えます。ただし、\(\mathbf{1}\)は\(1\)を縦に\(p\)個並べたベクトル、\(\circ\)はベクトルの要素ごとの積を表します。
このモデルは、修正変数\(Z\)の値によって、説明変数\(X_j\)にかかる回帰係数の値が変化するということを表しています。
Pliable Lassoでは、上のモデルに、さらに「\(\hat{\beta}_j = 0\)のときは、\(\hat{\theta}_j = 0\)である」という階層制約を加え、修正変数\(Z\)が単体で入力\(y\)の主要因になってしまうことを防ぎます。
このような回帰係数\(\beta\)および\(\theta\)は、次の関数を最小化することにより推定することができます。
\( J(\beta, \theta) = \frac{1}{2} \| y - \hat{y} \|_2^2 + (1 - \alpha) \lambda \sum_{j = 1}^p (\| (\beta_j, \theta_j)\|_2^2 + \|\theta_j\|_2^2) + \alpha \lambda \sum_{j, k} |\theta_{jk}| \)
ここで、第2項は先程の階層制約を表していて、第3はの回帰係数\(\hat{\theta}\)のスパース性を表しています。
それでは、実際に、Pliable Lassoをプログラムで動かしてみましょう。まず、Pliable Lassoのパッケージpliableをインストールします。
# Install package for Pliable Lasso
if(!require("pliable")){
url <- "https://cran.r-project.org/src/contrib/Archive/pliable/pliable_1.1.tar.gz"
pkgFile <- "pliable_1.1.tar.gz"
download.file(url = url, destfile = pkgFile)
install.packages(pkgs=pkgFile, type="source", repos=NULL)
unlink(pkgFile)
}
今回は、各要素を正規分布でランダムに初期化した行列\(X\)と、各要素を\(0\)または\(1\)でランダムに初期化した行列\(Z\)を用いて、\(y\)を
\( y = 2 X_1 - 2 X_2 + (2 + 2 Z_1) X_3 + (2 - 4 Z_2) X_4 + 0.5 \varepsilon \)
と生成したデータを用います。ただし、\(\varepsilon\)は標準正規分布に従うノイズとします。
このデータを使ってCross Validationを用いたPliable Lassoで回帰係数の推定を行ってみます。
fit_pliable <- pliable(x, z, y) cv_fit_pliable <- cv.pliable(fit_pliable, x, z, y, verbose=FALSE) fit_pliable <- pliable(x, z, y, lambda=cv_fit_pliable$lambda.min) print(coef(fit_pliable)) # $betaz # [1] 0.01100791 -0.06641263 0.10949380 0.03993577 # # $beta # [1] 1.96552143 -2.01533988 1.95736209 1.88119097 -0.08082027 -0.03948929 # [7] 0.00000000 -0.07799262 -0.02051081 -0.02314566 # # $theta # [,1] [,2] [,3] [,4] # [1,] 0.087065566 -0.000131529 0.000000000 -0.041261783 # [2,] 0.000000000 0.000000000 0.000000000 0.000000000 # [3,] 1.941422377 -0.034798144 0.000000000 -0.002797550 # [4,] -0.005704073 -3.955412349 0.000000000 0.026515585 # [5,] -0.031601711 0.000000000 0.069946806 -0.076491350 # [6,] -0.007395559 0.000000000 0.000000000 0.004435835 # [7,] 0.000000000 0.000000000 0.000000000 0.000000000 # [8,] 0.000000000 -0.031844665 0.058509824 0.000000000 # [9,] 0.009148531 0.000000000 -0.006090111 -0.001607628 # [10,] 0.000000000 0.016534991 -0.006355646 0.000000000)
出力結果と、上で\(y\)を生成した式を見比べてみると、回帰係数\(\hat{\beta}\)および\(\hat{\theta}\)の値がほぼ正しく推定できていることがわかります。
また、詳しい読者の方ならお気づきかもしれませんが、Pliable Lassoと似たモデルとして、交互作用モデルというものがあります。
これは、2種類の説明変数\(X\)と\(Z\)の組み合わせで表されるような要因が考えられる際に、両者をかけ合わせた交互作用項\(XZ\)を新たに説明変数に加えるようなモデルのことを指します。
実は、Pliable Lassoのモデルの式の
\( \hat{y} = \sum_{j = 1}^p X_j \circ (\hat{\beta}_j \mathbf{1} + Z \hat{\theta}_j) \)
は、次のように変形することで、\(X\)と\(Z\)の交互作用として表すことができます。
\( \hat{y} = X \hat{\beta} + \sum_{j = 1}^p (X_j \circ Z) \hat{\theta}_j \)
つまり、Pliable Lassoは、２種類の説明変数\(X\)と\(Z\)の交互作用を考慮した交互作用モデルに「\(\hat{\beta}_j = 0\)のときは、\(\hat{\theta}_j = 0\)である」という階層制約を加えたモデルと言えます。
そうすることにより、「\(Z\)の値によって、\(X_j\)にかかる回帰係数の値を修正する」という、交互作用モデルよりもより解釈性の高いモデルを実現しているのです。
最後に
本連載では、データ量が少なくても分析が実現できる「スパースモデリング」について、その歴史から最先端の世界までをご紹介してきました。
株式会社HACARUSでは、スパースモデリングを利用した産業・医療向けAIサービスを提供しています。本連載を通してスパースモデリングに興味を持っていただけた読者の皆様は、一度当社HPもご覧になっていただければ幸いです。
お読みいただき、ありがとうございました。
