lasso, adaptive lasso, group lasso[R]

https://www.math.mcgill.ca/yyang/comp/notes/note4code.R

library(glmnet)

このパッケージで使用されるデフォルトのモデルは、Guassian線形モデルまたは最小二乗法である。説明のためにあらかじめ作成したデータセットをロードする。ユーザは自分のデータをロードすることも、ワークスペースに保存されているデータを使用することもできる。

pacman::p_install_gh("emeryyi/gglasso")

library(gglasso)
data(bardet)
x<- bardet[["x"]]
y<- bardet[["y"]]

glmnetの最も基本的な呼び出しを使ってモデルをフィットさせる。

fit = glmnet(x, y)

"fit "はglmnetクラスのオブジェクトで、フィットされたモデルの関連情報をすべて含んでいる。我々は、ユーザが直接成分を抽出することを推奨していない。その代わりに、plot、print、coef、predictなどの様々なメソッドがオブジェクトに提供され、これらのタスクをよりエレガントに実行することができる。

plot関数を実行すれば、係数を可視化することができる：

plot(fit)

各曲線は変数に対応する。これは、λが変化する時の係数ベクトル全体の「1-norm」に対する係数のパスを示している。上の軸は、現在のλにおける非ゼロ係数の数を示しており、これはlassoの有効自由度(df)となる。ユーザーは、曲線に注釈を付けることもでき、これはplotコマンドでlabel = TRUEを設定することで可能だ。

オブジェクト名を入力するか、print機能を使えば、各ステップでのglmnetパスの概要が表示される：

print(fit)

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

    Df  %Dev   Lambda
1    0  0.00 0.097200
2    1  7.73 0.088560
3    1 14.16 0.080690
4    1 19.49 0.073520
5    1 23.91 0.066990
6    1 27.59 0.061040
7    2 31.42 0.055620
8    2 36.05 0.050680
9    4 40.59 0.046180
10   6 45.05 0.042070
...

左からゼロでない係数の数(Df)、-log(尤度)の値(%dev)、ラムダの値(Lambda)が表示される。デフォルトではglmnetは100個のλの値を要求するが、%dev%がλの値から次のλの値まで十分に変化しない場合、プログラムは早期に停止する(典型的にはパスの終わり近く)。

配列の範囲内の1つ以上のラムダで実際の係数を得ることができる：

coef0 = coef(fit,s=0.1)

関数glmnetは、ユーザーが選択できる一連のモデルを返す。多くの場合、ユーザーはソフトウェアがその中から1つを選択することを好むかもしれない。クロスバリデーションは、おそらくそのタスクのための最もシンプルで最も広く使われている方法である。

cv.glmnetは、プロットや予測などの様々なサポートメソッドとともに、クロスバリデーションを行うためのメイン関数です。ここでは、前に読み込んだサンプルデータを使っている。

cvfit = cv.glmnet(x, y)

cv.glmnetは、cv.glmnetオブジェクトを返しますが、これはここでは "cvfit"であり、クロスバリデーション適合のすべての成分を含むリストです。glmnetに関しては、選択されたλの値を見る以外には、ユーザが直接成分を抽出することは推奨しません。このパッケージは、潜在的なタスクのためによく設計された関数を提供する。

オブジェクトをプロットできる。

plot(cvfit)

クロスバリデーション曲線（赤い点線）、ラムダ配列に沿った上下の標準偏差曲線（エラーバー）が含まれる。選択された2つのラムダは縦の点線で示されている（下図参照）。

選択されたラムダと対応する係数を見ることができる。例えば、

cvfit$lambda.min

[1] 0.001226463

lambda.minは、クロスバリデーションの平均誤差を最小にするλの値である。もう1つの保存されたλはλ.1seで、これは誤差が最小の標準誤差の1つ以内に収まるような最も正則化されたモデルを与える。これを使うには、上記のlambda.minをlambda.1seに置き換えるだけでよい。

coef1 = coef(cvfit, s = "lambda.min")

係数はスパース行列形式で表現されていることに注意。その理由は、正則化パスに沿った解はスパースであることが多いため、スパースフォーマットを使用する方が時間的にも空間的にも効率的だからだ。スパースでない形式を好む場合は、出力をas.matrix()にパイプを通して出力する。

予測は、フィットされたcv.glmnetオブジェクトに基づいて行うことができる。おもちゃの例を見てみよう。

predict(cvfit, newx = x[1:5,], s = "lambda.min")

  lambda.min

[1,] 8.406305 [2,] 8.353781 [3,] 8.377383 [4,] 8.294958 [5,] 8.379072

newxは新しい入力行列を表し、sは前回と同様、予測が行われるλの値（複数可）である。

adaptive lasso

ペナルティ・ファクターの例

[penalty.factor]引数により、各係数に別々のペナルティ係数を適用することができる。デフォルトは各パラメータに対して1だが、他の値を指定することもできる。特に、[penalty.factor]が0に等しい変数は、全くペナルティを受けません！v_jがj番目の変数の制約を表すとする。

ペナルティ制約は、内部的にnvarsの合計になるように再スケーリングされることに注意。

これは、変数に対する予備知識や好みがある場合に非常に便利だ。多くの場合、ある変数が非常に重要で、それを常に保持したい場合があり、これは対応する制約係数を0に設定することで実現できる： v_j は、ペナルティ係数を0に設定することで、ペナルティ係数を0にすることができる：

p.fac = rep(1, 120)
p.fac[c(31, 74, 111)] = 0
pfit = glmnet(x, y, penalty.factor = p.fac)
plot(pfit, xvar="lambda", label = TRUE)

その他の便利な引数をいくつか挙げる： [exclude]は、特定の変数をモデルにしないようにすることができる。もちろん、単純にxからこれらの変数を除外することもできるが、除外された変数をゼロに設定するだけで、係数の完全なベクトルを返すので、excludeの方が便利な場合もある。FALSEの場合、切片は強制的にゼロになる。

Adaptive lassoの例

adaptive lassoは一貫性のある第1段階を必要とする。 Zou (2006)はOLSかridge第一段lassoを推奨している。

thelasso.cv<-cv.glmnet(x,y,family = "gaussian",alpha=1) 

第1段階の係数から第2段階の重みを求める coef() はスパース行列

bhat<-as.matrix(coef(thelasso.cv,s="lambda.1se"))[-1,1] 
if(all(bhat==0)){
  bhat<-rep(.Machine$double.eps*2,length(bhat))
}

もしbhatがすべてゼロなら、すべてにゼロに極めて近いウェイトを割り当てる。これは、第2ステージのすべてにゼロのペナルティを課すことに等しい。

adaptive lasso weight

adpen<-(1/pmax(abs(bhat),.Machine$double.eps)) 

第2段階のlasso（adaptive lasso）

m_adlasso <- glmnet(x,y,family = "gaussian",alpha=1,exclude=which(bhat==0),
penalty.factor=adpen)
plot(m_adlasso)

エラスティックネットelastic-net

エラスティック・ネットの例

glmnetには、ユーザーがフィットをカスタマイズするための様々なオプションが用意されている。ここではよく使われるオプションをいくつか紹介するが、これらはglmnet関数で指定することができる。

[family="gaussian"]はglmnet関数のデフォルトのfamilyオプションだ。

[alpha]は弾性ネットの混合パラメータalphaで、alphaの範囲は[0,1]です。alpha=1はlasso（デフォルト）、alpha=0はリッジになる。

[nlambda] はシーケンス中のラムダ値の数であり、デフォルトは100ですが、[0,1]の範囲もある。

例として、α=0.2（よりリッジ回帰に近い）を設定し、オブザベーションの後半に2倍の重みを与える。 ここでの表示が長くなりすぎるのを避けるために、nlambdaを20に設定する。しかし実際には、λの値の数は100（デフォルト）以上にすることが推奨される。ほとんどの場合、アルゴリズムで使用されるウォーム・スタート warm-starts のために、余分なコストがかからず、非線形モデルではよりよい収束特性につながる。

fit = glmnet(x, y, alpha = 0.2, family="gaussian")

plot(fit, xvar = "lambda", label = TRUE)

グループlasso

グループlassoの例

install.packages("gglasso", repos = "http://cran.us.r-project.org")

https://cran.r-project.org/web/packages/gglasso/index.html

group-lasso ペナルティ付き最小二乗法、ロジスティック回帰、Huberized SVM、二乗SVMの解経路を効率的に計算するための統一アルゴリズム、blockwise-majorization-descent (BMD)。 このパッケージは、Yang, Y. and Zou, H. (2015) DOI: <doi:10.1007/s11222-014-9498-5> の実装だ。

library(gglasso)

bardetデータセットの読み込み

data(bardet)

20のグループを定義する

group1 <- rep(1:20,each=5)

ペナルティ付き最小二乗"group lasso"にフィットさせる。

m1 <- gglasso(x=bardet$x,y=bardet$y,group=group1,loss="ls")

plot(m1)

ペナルティ付きグループLassoでls回帰を用いた5-fold クロスバリデーション

cv <- cv.gglasso(x=bardet$x, y=bardet$y, group=group1, loss="ls",
pred.loss="L2", lambda.factor=0.05, nfolds=5)

plot(cv)

colon データセットを読み込む

data(colon)

グループ指標を定義する

group2 <- rep(1:20,each=5)

グループlassoペナルティ付きロジスティック回帰にフィット

m2 <- gglasso(x=colon$x,y=colon$y,group=group2,loss="logit")

plot(m2)

グループlassoペナルティ付きロジスティック回帰を用いた5-foldクロスバリデーション

cv2 <- cv.gglasso(x=colon$x, y=colon$y, group=group2, loss="logit",
pred.loss="misclass", lambda.factor=0.05, nfolds=5)

plot(cv2)

λ=λ.1seにおける係数

pre = coef(cv$gglasso.fit, s = cv$lambda.1se)

シミュレーション・デモ-1

ロジスティック回帰と最小二乗法におけるlasso、adaptive lasso、scad、mcpの比較

Rライブラリをロード
glmnet: LASSO、adaptive LASSO、elastic net ペナルティ付き最小2乗法およびロジスティック回帰 (すべてポアソン、多項式、コックスモデルをサポート)
ncvreg: SCAD および MCP ペナルティ化最小2乗法およびロジスティック回帰用

library(glmnet)
library(ncvreg)
library(MASS)

PART I 最小二乗法

n=100
p=200

真のベータ

truebeta <- c(4,4,4,-6*sqrt(2),4/3,rep(0,p-5))

誤差分散

sigma2 <- 0.3

XjとXkの間の共分散はcov(X_j, X_k) = rho^|i-j|共分散行列である。

covmat <- function(rho, p) {
  rho^(abs(outer(seq(p), seq(p), "-")))
}

rho = 0.1、Xの共分散行列を生成する。

sigma <- covmat(0.1,p)

X ~ N(0, sigma)
epsilon ~ N(0, sigma2)
真のモデル： y = x*truebeta + epsilon

x <- mvrnorm(n,rep(0,p),sigma)
epsilon <- rnorm(n,0,sd=sqrt(sigma2))
y <- x %*% truebeta + epsilon

lassoをフィットさせ、5-fold CVでラムダを選択する。

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, family = "gaussian")
plot(cvfit)

CVで選ばれたラムダ

cvfit$lambda.min

[1] 0.0821079

解決パスをプロットする

plot(cvfit$glmnet.fit,label=TRUE,xvar="lambda")

λ.minのプロット

abline(v=log(cvfit$lambda.min),lty=1)

λ.1seをプロットする

abline(v=log(cvfit$lambda.1se),lty=2)

model_compare <- matrix(NA, nrow=5, ncol=p,
                dimnames=list(c("true model",
                "lasso","adaptive lasso","mcp","scad"),
                paste("V",seq(p),sep="")))

真のモデルを保存する

model_compare[1, ] <- truebeta

lassoによって推定された係数

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, 
family = "gaussian")
tmp <- cvfit$glmnet.fit$beta
lasso_beta <- as.matrix(tmp[,cvfit$lambda==cvfit$lambda.min])
model_compare[2, ] <- lasso_beta

重みを計算し、adaptive lassoにフィットさせる

weight = 1/(lasso_beta)^2

推定係数がゼロの場合、対応するウェイトがInfとなり、数値誤差を防ぐため、Infを大きな数値（例えば1e6）に変換する。

weight[weight==Inf] = 1e6

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, 
        family = "gaussian", 
        nfolds = 5, penalty.factor=weight)

Adaptive Lassoによって推定された係数

tmp <- cvfit$glmnet.fit$beta
adaptive_lasso_beta <- as.matrix(tmp[,cvfit$lambda==cvfit$lambda.min])
model_compare[3, ] <- adaptive_lasso_beta

mcpによって推定された係数

cvfit <- cv.ncvreg(X = x, y = y, penalty = "MCP", family = "gaussian")
mcp_beta <- cvfit$fit$beta[, cvfit$min]
model_compare[4, ] <- mcp_beta[-1]

scadによって推定された係数

cvfit <- cv.ncvreg(X = x, y = y, penalty = "SCAD", family = "gaussian")
scad_beta <- cvfit$fit$beta[, cvfit$min]
model_compare[5, ] <- scad_beta[-1]

4つの方法で推定された係数を比較すると、lasso over-selected、adaptive lasso、scad、mcpが問題を解決していることがわかる。

model_compare

               V1           V2 V3        V4       V5 V6 V7 V8 V9 V10
true model      4  4.000000000  4 -8.485281 1.333333  0  0  0  0   0
lasso           0 -0.022023796  0  0.000000 0.000000  0  0  0  0   0
adaptive lasso  0  0.000000000  0  0.000000 0.000000  0  0  0  0   0
mcp             0  0.000000000  0  0.000000 0.000000  0  0  0  0   0
scad            0 -0.004545409  0  0.000000 0.000000  0  0  0  0   0
                       V11         V12        V13       V14       V15 V16
true model     0.000000000  0.00000000  0.0000000 0.0000000  0.000000   0
lasso          0.006311609  0.00000000  0.0000000 0.0000000  0.000000   0
adaptive lasso 0.000000000  0.00000000  0.0000000 0.0000000  0.000000   0
mcp            0.000000000 -0.08207367 -0.1647326 0.2088095 -1.000567   0
scad           0.000000000  0.00000000  0.0000000 0.0000000 -1.033549   0
                      V17         V18 V19 V20 V21          V22         V23
true model     0.00000000  0.00000000   0   0   0 0.0000000000  0.00000000
lasso          0.00000000 -0.01762063   0   0   0 0.0001952964 -0.05328850
adaptive lasso 0.00000000  0.00000000   0   0   0 0.0000000000  0.00000000
mcp            0.03504008  0.00000000   0   0   0 0.0000000000  0.00000000
scad           0.00000000  0.00000000   0   0   0 0.0000000000 -0.08488739
                      V24        V25        V26 V27        V28 V29        V30
true model      0.0000000  0.0000000 0.00000000   0  0.0000000   0  0.0000000
lasso           0.0000000 -0.5648782 0.06002556   0  0.0000000   0 -0.3371098
adaptive lasso  0.0000000 -1.0114037 0.00000000   0  0.0000000   0  0.0000000
mcp            -0.4379709  0.0000000 0.00000000   0 -0.1994319   0  0.0000000
scad            0.0000000  0.0000000 0.00000000   0  0.0000000   0  0.0000000
                      V31 V32        V33 V34 V35       V36 V37 V38 V39 V40 V41
true model      0.0000000   0  0.0000000   0   0 0.0000000   0   0   0   0   0
lasso          -0.1165544   0  0.0000000   0   0 0.0482393   0   0   0   0   0
adaptive lasso  0.0000000   0  0.0000000   0   0 0.0000000   0   0   0   0   0
mcp            -0.4561359   0 -0.1259705   0   0 0.1201394   0   0   0   0   0
scad            0.0000000   0  0.0000000   0   0 0.0000000   0   0   0   0   0
               V42 V43         V44 V45        V46       V47 V48 V49        V50
true model       0   0 0.000000000   0  0.0000000 0.0000000   0   0  0.0000000
lasso            0   0 0.009982214   0 -0.2897772 0.0000000   0   0  0.0000000
adaptive lasso   0   0 0.000000000   0 -0.6616183 0.0000000   0   0  0.0000000
mcp              0   0 0.000000000   0 -0.4648836 0.1156758   0   0 -0.1438189
scad             0   0 0.000000000   0 -0.1163973 0.0000000   0   0  0.0000000
                      V51 V52        V53        V54 V55 V56 V57         V58 V59
true model      0.0000000   0 0.00000000 0.00000000   0   0   0  0.00000000   0
lasso           0.0000000   0 0.05634673 0.13775797   0   0   0  0.00000000   0
adaptive lasso  0.0000000   0 0.00000000 0.09460430   0   0   0  0.00000000   0
mcp             0.0000000   0 0.07499241 0.11837706   0   0   0 -0.05640483   0
scad           -0.1008567   0 0.00000000 0.08972644   0   0   0  0.00000000   0
               V60 V61 V62        V63 V64       V65         V66 V67 V68 V69 V70
true model       0   0   0 0.00000000   0 0.0000000  0.00000000   0   0   0   0
lasso            0   0   0 0.00000000   0 0.2563789 -0.12041928   0   0   0   0
adaptive lasso   0   0   0 0.00000000   0 0.2957797  0.00000000   0   0   0   0
mcp              0   0   0 0.02528918   0 0.3411462  0.00000000   0   0   0   0
scad             0   0   0 0.00000000   0 0.1248035 -0.03840615   0   0   0   0
               V71 V72        V73 V74 V75         V76       V77        V78 V79
true model       0   0 0.00000000   0   0  0.00000000 0.0000000 0.00000000   0
lasso            0   0 0.05161909   0   0  0.00000000 0.0000000 0.04393271   0
adaptive lasso   0   0 0.00000000   0   0  0.00000000 0.0000000 0.00000000   0
mcp              0   0 0.07055282   0   0  0.00000000 0.0652188 0.12140394   0
scad             0   0 0.00000000   0   0 -0.06290252 0.0000000 0.00000000   0
               V80 V81 V82 V83 V84 V85        V86          V87 V88       V89
true model       0   0   0   0   0   0 0.00000000  0.000000000   0 0.0000000
lasso            0   0   0   0   0   0 0.00000000 -0.001347216   0 0.2044193
adaptive lasso   0   0   0   0   0   0 0.00000000  0.000000000   0 0.1554728
mcp              0   0   0   0   0   0 0.04738824  0.000000000   0 0.4817615
scad             0   0   0   0   0   0 0.00000000  0.000000000   0 0.0000000
               V90 V91 V92 V93        V94 V95 V96 V97 V98 V99       V100 V101
true model       0   0   0   0 0.00000000   0   0   0   0   0 0.00000000    0
lasso            0   0   0   0 0.07206646   0   0   0   0   0 0.06736159    0
adaptive lasso   0   0   0   0 0.00000000   0   0   0   0   0 0.00000000    0
mcp              0   0   0   0 0.44369829   0   0   0   0   0 0.14673570    0
scad             0   0   0   0 0.00000000   0   0   0   0   0 0.00000000    0
                       V102 V103 V104 V105 V106 V107 V108 V109 V110        V111
true model      0.000000000    0    0    0    0    0    0    0    0 0.000000000
lasso          -0.022023796    0    0    0    0    0    0    0    0 0.006311609
adaptive lasso  0.000000000    0    0    0    0    0    0    0    0 0.000000000
mcp             0.000000000    0    0    0    0    0    0    0    0 0.000000000
scad           -0.004545409    0    0    0    0    0    0    0    0 0.000000000
                      V112       V113      V114      V115 V116       V117
true model      0.00000000  0.0000000 0.0000000  0.000000    0 0.00000000
lasso           0.00000000  0.0000000 0.0000000  0.000000    0 0.00000000
adaptive lasso  0.00000000  0.0000000 0.0000000  0.000000    0 0.00000000
mcp            -0.08207367 -0.1647326 0.2088095 -1.000567    0 0.03504008
scad            0.00000000  0.0000000 0.0000000 -1.033549    0 0.00000000
                      V118 V119 V120 V121         V122        V123       V124
true model      0.00000000    0    0    0 0.0000000000  0.00000000  0.0000000
lasso          -0.01762063    0    0    0 0.0001952964 -0.05328850  0.0000000
adaptive lasso  0.00000000    0    0    0 0.0000000000  0.00000000  0.0000000
mcp             0.00000000    0    0    0 0.0000000000  0.00000000 -0.4379709
scad            0.00000000    0    0    0 0.0000000000 -0.08488739  0.0000000
                     V125       V126 V127       V128 V129       V130       V131
true model      0.0000000 0.00000000    0  0.0000000    0  0.0000000  0.0000000
lasso          -0.5648782 0.06002556    0  0.0000000    0 -0.3371098 -0.1165544
adaptive lasso -1.0114037 0.00000000    0  0.0000000    0  0.0000000  0.0000000
mcp             0.0000000 0.00000000    0 -0.1994319    0  0.0000000 -0.4561359
scad            0.0000000 0.00000000    0  0.0000000    0  0.0000000  0.0000000
               V132       V133 V134 V135      V136 V137 V138 V139 V140 V141
true model        0  0.0000000    0    0 0.0000000    0    0    0    0    0
lasso             0  0.0000000    0    0 0.0482393    0    0    0    0    0
adaptive lasso    0  0.0000000    0    0 0.0000000    0    0    0    0    0
mcp               0 -0.1259705    0    0 0.1201394    0    0    0    0    0
scad              0  0.0000000    0    0 0.0000000    0    0    0    0    0
               V142 V143        V144 V145       V146      V147 V148 V149
true model        0    0 0.000000000    0  0.0000000 0.0000000    0    0
lasso             0    0 0.009982214    0 -0.2897772 0.0000000    0    0
adaptive lasso    0    0 0.000000000    0 -0.6616183 0.0000000    0    0
mcp               0    0 0.000000000    0 -0.4648836 0.1156758    0    0
scad              0    0 0.000000000    0 -0.1163973 0.0000000    0    0
                     V150       V151 V152       V153       V154 V155 V156 V157
true model      0.0000000  0.0000000    0 0.00000000 0.00000000    0    0    0
lasso           0.0000000  0.0000000    0 0.05634673 0.13775797    0    0    0
adaptive lasso  0.0000000  0.0000000    0 0.00000000 0.09460430    0    0    0
mcp            -0.1438189  0.0000000    0 0.07499241 0.11837706    0    0    0
scad            0.0000000 -0.1008567    0 0.00000000 0.08972644    0    0    0
                      V158 V159 V160 V161 V162       V163 V164      V165
true model      0.00000000    0    0    0    0 0.00000000    0 0.0000000
lasso           0.00000000    0    0    0    0 0.00000000    0 0.2563789
adaptive lasso  0.00000000    0    0    0    0 0.00000000    0 0.2957797
mcp            -0.05640483    0    0    0    0 0.02528918    0 0.3411462
scad            0.00000000    0    0    0    0 0.00000000    0 0.1248035
                      V166 V167 V168 V169 V170 V171 V172       V173 V174 V175
true model      0.00000000    0    0    0    0    0    0 0.00000000    0    0
lasso          -0.12041928    0    0    0    0    0    0 0.05161909    0    0
adaptive lasso  0.00000000    0    0    0    0    0    0 0.00000000    0    0
mcp             0.00000000    0    0    0    0    0    0 0.07055282    0    0
scad           -0.03840615    0    0    0    0    0    0 0.00000000    0    0
                      V176      V177       V178 V179 V180 V181 V182 V183 V184
true model      0.00000000 0.0000000 0.00000000    0    0    0    0    0    0
lasso           0.00000000 0.0000000 0.04393271    0    0    0    0    0    0
adaptive lasso  0.00000000 0.0000000 0.00000000    0    0    0    0    0    0
mcp             0.00000000 0.0652188 0.12140394    0    0    0    0    0    0
scad           -0.06290252 0.0000000 0.00000000    0    0    0    0    0    0
               V185       V186         V187 V188      V189 V190 V191 V192 V193
true model        0 0.00000000  0.000000000    0 0.0000000    0    0    0    0
lasso             0 0.00000000 -0.001347216    0 0.2044193    0    0    0    0
adaptive lasso    0 0.00000000  0.000000000    0 0.1554728    0    0    0    0
mcp               0 0.04738824  0.000000000    0 0.4817615    0    0    0    0
scad              0 0.00000000  0.000000000    0 0.0000000    0    0    0    0
                     V194 V195 V196 V197 V198 V199       V200
true model     0.00000000    0    0    0    0    0 0.00000000
lasso          0.07206646    0    0    0    0    0 0.06736159
adaptive lasso 0.00000000    0    0    0    0    0 0.00000000
mcp            0.44369829    0    0    0    0    0 0.14673570
scad           0.00000000    0    0    0    0    0 0.00000000

シミュレーション・デモ-2

PART II ロジスティック回帰

データを生成する

n <- 200 
p <- 8

真のベータ

truebeta <- c(6,3.5,0,5,rep(0,p-4)) 
truebeta

真のモデルからxとyを生成 真のモデル： P(y=1|x) = exp(xtruebeta)/(exp(xtruebeta)+1)

x <- matrix(rnorm(n*p), n, p)
feta <- x %*% truebeta 
fprob <- ifelse(feta < 0, exp(feta)/(1+exp(feta)), 1/(1 + exp(-feta)))
y <- rbinom(n, 1, fprob)


model_compare <- matrix(NA, nrow=5, ncol=p,
                        dimnames=list(c("true model","lasso",
                        "adaptive lasso","mcp","scad"),
                        paste("V",seq(p),sep="")))

真のモデルを保存する

model_compare[1, ] <- truebeta

lassoのケース
cv.glmfit ラッソモデルをフィットし、ラムダの選択にクロス検証を使用する。

family = "binomial", ロジスティック回帰
family = "gaussian", 最小二乗法

alphaはL1ペナルティ項の度合いを制御する、
alpha = 1, lasso、
alpha = 0, ridge、
alpha = (0,1), 弾性ネット

nfolds = 5, 5回クロスバリデーション(CV)

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, family = "binomial", nfolds = 5)

CV結果をプロットする 左の縦線(λ.min)は、最小の偏差を与えるλに対応する。
右の縦線(λ.1se)は、1標準偏差ルールのλに対応する。

plot(cvfit)

解決パスをプロットする

plot(cvfit$glmnet.fit,label=TRUE,xvar="lambda")

plot(cvfit$glmnet.fit,label=TRUE,xvar="lambda")

λ.minをプロットする

abline(v=log(cvfit$lambda.min),lty=1)

λ.1seをプロットする

abline(v=log(cvfit$lambda.1se),lty=2)

各列はラムダ値の推定値を表す。

tmp <- cvfit$glmnet.fit$beta
tmp

CVが選択したλ.minに対応するベータ値

lasso_beta <- as.matrix(tmp[,cvfit$lambda==cvfit$lambda.min])
model_compare[2, ] <- lasso_beta

重みを計算し、Adaptive Lassoをフィットさせる

weight = 1/(lasso_beta)^2

推定係数がゼロの場合、対応するウェイトがInfとなり、数値誤差を防ぐため、Infを大きな数値（例えば1e6）に変換する。

weight[weight==Inf] = 1e6

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, family = "binomial", 
                        nfolds = 5, penalty.factor=weight)

Adaptive Lassoによって推定された係数

tmp <- cvfit$glmnet.fit$beta
adaptive_lasso_beta <- as.matrix(tmp[,cvfit$lambda==cvfit$lambda.min])
model_compare[3, ] <- adaptive_lasso_beta

mcpによって推定された係数

cvfit <- cv.ncvreg(X = x, y = y, penalty = "MCP", family = "binomial")
mcp_beta <- cvfit$fit$beta[, cvfit$min]
model_compare[4, ] <- mcp_beta[-1]

scadによって推定された係数

cvfit <- cv.ncvreg(X = x, y = y, penalty = "SCAD", family = "binomial")
scad_beta <- cvfit$fit$beta[, cvfit$min]
model_compare[5, ] <- scad_beta[-1]

4つの方法から推定された係数を比較する。

model_compare

                     V1       V2 V3       V4 V5          V6         V7
true model           NA       NA NA       NA NA          NA         NA
lasso          4.467528 2.375300  0 3.096532  0 -0.07315415 -0.3499353
adaptive lasso 4.664295 2.480856  0 3.335227  0  0.00000000  0.0000000
mcp            4.459921 2.533341  0 3.366028  0  0.00000000  0.0000000
scad           4.808539 2.744814  0 3.615375  0  0.00000000  0.0000000
                      V8
true model            NA
lasso          0.3801377
adaptive lasso 0.0000000
mcp            0.0000000
scad           0.0000000