久しぶりの更新です（いつも言っています）。
背景
データサイエンス入門シリーズの「スパース回帰分析とパターン認識」を読んでいたら大変面白かったので、いつものように glmnet の中身を見てみることにしました。
なお私は業務でLasso/Ridgeを使った経験があまりないため理解が間違っているかもしれませんが、その点あらかじめご了承ください。
こちらの本です。良い本です。
glmnet の実行結果
前回の GAM の時と同様に、まずは glmnet でどのような結果を得ることができるのか確認してみましょう。「スパース回帰分析とパターン認識」（以下、教科書）P12 コード1.2を（少し改変して）実行してみます。
なおこれらのコードはこちらからダウンロードすることもできます。
環境は以下のような感じです。
> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] compiler_3.6.0 tools_3.6.0 grid_3.6.0 lattice_0.20-38
library(glmnet)
library(plotmo)
x <- scale(LifeCycleSavings[, 2:5])
y <- LifeCycleSavings[, 1] - mean(LifeCycleSavings[, 1])
lasso <- glmnet(x, y, family = "gaussian", alpha = 1)
ridge <- glmnet(x, y, family = "gaussian", alpha = 0)
png("./Image/glmnet_dive_01_01.png", width = 600, height = 400)
plot_glmnet(lasso, xvar = "lambda", label = TRUE)
dev.off()
png("./Image/glmnet_dive_01_02.png", width = 600, height = 400)
plot_glmnet(ridge, xvar = "lambda", label = TRUE)
dev.off()
結果の解釈などについて詳しくは教科書を見て頂くとして、 glmnet は目的関数に回帰係数の規模に応じた罰則を設けることで、回帰係数を0に向かって縮小させながらフィッティングを行います。
またグラフのように罰則の大きさを色々と動かすことで各変数への回帰係数がどのように変化するかを評価することができます。
このグラフでは左から右に向かって罰則が強くなりますが、それにつれてLasso/Ridgeの両方とも回帰係数が0に向かって小さくなっている（縮小している）ことがわかります。
なお Lasso では回帰係数が0に収束している一方、 Ridge では微小ながら最後まで係数が0とならずに残っていることがわかりますが（グラフ上部の Degrees of Freedom が 4 のままとなっています）、 Lasso のように一部の回帰係数を正確に 0 と推定することが可能な手法をスパース推定と呼びます。
glmnet の実装
それでは glmnet という関数がどのように実装されているのか見ていきましょう。
まずはいつものように全体を眺め、見通しをよくします。
function (x, y, family = c("gaussian", "binomial", "poisson",
"multinomial", "cox", "mgaussian"), weights, offset = NULL,
alpha = 1, nlambda = 100, lambda.min.ratio = ifelse(nobs <
nvars, 0.01, 1e-04), lambda = NULL, standardize = TRUE,
intercept = TRUE, thresh = 1e-07, dfmax = nvars + 1, pmax = min(dfmax *
2 + 20, nvars), exclude, penalty.factor = rep(1, nvars),
lower.limits = -Inf, upper.limits = Inf, maxit = 1e+05, type.gaussian = ifelse(nvars <
500, "covariance", "naive"), type.logistic = c("Newton",
"modified.Newton"), standardize.response = FALSE, type.multinomial = c("ungrouped",
"grouped"), relax = FALSE, trace.it = 0, ...)
{
family = match.arg(family)
if (alpha > 1) {
warning("alpha >1; set to 1")
alpha = 1
}
if (alpha < 0) {
warning("alpha<0; set to 0")
alpha = 0
}
alpha = as.double(alpha)
this.call = match.call()
nlam = as.integer(nlambda)
y = drop(y)
np = dim(x)
if (is.null(np) | (np[2] <= 1))
stop("x should be a matrix with 2 or more columns")
nobs = as.integer(np[1])
if (missing(weights))
weights = rep(1, nobs)
else if (length(weights) != nobs)
stop(paste("number of elements in weights (", length(weights),
") not equal to the number of rows of x (", nobs,
")", sep = ""))
nvars = as.integer(np[2])
dimy = dim(y)
nrowy = ifelse(is.null(dimy), length(y), dimy[1])
if (nrowy != nobs)
stop(paste("number of observations in y (", nrowy, ") not equal to the number of rows of x (",
nobs, ")", sep = ""))
vnames = colnames(x)
if (is.null(vnames))
vnames = paste("V", seq(nvars), sep = "")
ne = as.integer(dfmax)
nx = as.integer(pmax)
if (missing(exclude))
exclude = integer(0)
if (any(penalty.factor == Inf)) {
exclude = c(exclude, seq(nvars)[penalty.factor == Inf])
exclude = sort(unique(exclude))
}
if (length(exclude) > 0) {
jd = match(exclude, seq(nvars), 0)
if (!all(jd > 0))
stop("Some excluded variables out of range")
penalty.factor[jd] = 1
jd = as.integer(c(length(jd), jd))
}
else jd = as.integer(0)
vp = as.double(penalty.factor)
internal.parms = glmnet.control()
if (internal.parms$itrace)
trace.it = 1
else {
if (trace.it) {
glmnet.control(itrace = 1)
on.exit(glmnet.control(itrace = 0))
}
}
if (any(lower.limits > 0)) {
stop("Lower limits should be non-positive")
}
if (any(upper.limits < 0)) {
stop("Upper limits should be non-negative")
}
lower.limits[lower.limits == -Inf] = -internal.parms$big
upper.limits[upper.limits == Inf] = internal.parms$big
if (length(lower.limits) < nvars) {
if (length(lower.limits) == 1)
lower.limits = rep(lower.limits, nvars)
else stop("Require length 1 or nvars lower.limits")
}
else lower.limits = lower.limits[seq(nvars)]
if (length(upper.limits) < nvars) {
if (length(upper.limits) == 1)
upper.limits = rep(upper.limits, nvars)
else stop("Require length 1 or nvars upper.limits")
}
else upper.limits = upper.limits[seq(nvars)]
cl = rbind(lower.limits, upper.limits)
if (any(cl == 0)) {
fdev = glmnet.control()$fdev
if (fdev != 0) {
glmnet.control(fdev = 0)
on.exit(glmnet.control(fdev = fdev))
}
}
storage.mode(cl) = "double"
isd = as.integer(standardize)
intr = as.integer(intercept)
if (!missing(intercept) && family == "cox")
warning("Cox model has no intercept")
jsd = as.integer(standardize.response)
thresh = as.double(thresh)
if (is.null(lambda)) {
if (lambda.min.ratio >= 1)
stop("lambda.min.ratio should be less than 1")
flmin = as.double(lambda.min.ratio)
ulam = double(1)
}
else {
flmin = as.double(1)
if (any(lambda < 0))
stop("lambdas should be non-negative")
ulam = as.double(rev(sort(lambda)))
nlam = as.integer(length(lambda))
}
is.sparse = FALSE
ix = jx = NULL
if (inherits(x, "sparseMatrix")) {
is.sparse = TRUE
x = as(x, "CsparseMatrix")
x = as(x, "dgCMatrix")
ix = as.integer(x@p + 1)
jx = as.integer(x@i + 1)
x = as.double(x@x)
}
if (trace.it) {
if (relax)
cat("Training Fit\n")
pb <- createPB(min = 0, max = nlam, initial = 0, style = 3)
}
kopt = switch(match.arg(type.logistic), Newton = 0, modified.Newton = 1)
if (family == "multinomial") {
type.multinomial = match.arg(type.multinomial)
if (type.multinomial == "grouped")
kopt = 2
}
kopt = as.integer(kopt)
fit = switch(family, gaussian = elnet(x, is.sparse, ix, jx,
y, weights, offset, type.gaussian, alpha, nobs, nvars,
jd, vp, cl, ne, nx, nlam, flmin, ulam, thresh, isd, intr,
vnames, maxit), poisson = fishnet(x, is.sparse, ix, jx,
y, weights, offset, alpha, nobs, nvars, jd, vp, cl, ne,
nx, nlam, flmin, ulam, thresh, isd, intr, vnames, maxit),
binomial = lognet(x, is.sparse, ix, jx, y, weights, offset,
alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, flmin,
ulam, thresh, isd, intr, vnames, maxit, kopt, family),
multinomial = lognet(x, is.sparse, ix, jx, y, weights,
offset, alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam,
flmin, ulam, thresh, isd, intr, vnames, maxit, kopt,
family), cox = coxnet(x, is.sparse, ix, jx, y, weights,
offset, alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam,
flmin, ulam, thresh, isd, vnames, maxit), mgaussian = mrelnet(x,
is.sparse, ix, jx, y, weights, offset, alpha, nobs,
nvars, jd, vp, cl, ne, nx, nlam, flmin, ulam, thresh,
isd, jsd, intr, vnames, maxit))
if (trace.it) {
utils::setTxtProgressBar(pb, nlam)
close(pb)
}
if (is.null(lambda))
fit$lambda = fix.lam(fit$lambda)
fit$call = this.call
fit$nobs = nobs
class(fit) = c(class(fit), "glmnet")
if (relax)
relax.glmnet(fit, x = x, y = y, weights = weights, offset = offset,
lower.limits = lower.limits, upper.limits = upper.limits,
check.args = FALSE, ...)
else fit
}
glmnet では以上のように、
- パラメータの設定、前処理、エラーチェック
- フィッティング
- 後処理
といったステップで処理が進んでおり、これは過去にみてきた glm や gam と同様ですね。
それでは各ステップを細かく見ていきましょう。
1. パラメータの設定、前処理、エラーチェック
まずはパラメータの設定や前処理に関わる部分ですが、はじめに family の指定が問題ないかをチェックします。
family = match.arg(family)
glmnet で使用可能な family は glm とは異なっており、Gamma / inverse.gaussian / quasi- が使えない代わりに、 multinomial / cox / mgaussian が使えるようになっています。
ここで multinomial は多項分布、mgaussian は多変量正規分布を意味するようです。
family のチェックには match.arg 関数が使われています。
この関数の挙動を理解するのは少し難しいのですが、こちらのブログが参考になります。
続いて alpha をチェックします：
if (alpha > 1) {
warning("alpha >1; set to 1")
alpha = 1
}
if (alpha < 0) {
warning("alpha<0; set to 0")
alpha = 0
}
alpha = as.double(alpha)
glmnet においてこの alpha は、回帰係数のL1およびL2ノルムそれぞれに対する罰則の割合をコントロールします。
より具体的には、 glmnet では罰則項は以下によって定義されます（https://cran.r-project.org/web/packages/glmnet/glmnet.pdf の P19より）：
冒頭のコードでは alpha = 1 または alpha = 0 としましたが、上の式から alpha = 1 のときにL2ノルムに対する罰則が消えてL1ノルムのみが残り（Lasso）、逆に alpha = 0 とするとL1ノルムに対する罰則が消えてL2ノルムが残る（Ridge）ことがわかります。
また alpha を (0, 1) とすると両者がそれぞれの割合でブレンドされます。
なお、ここでL2ノルムに対する罰則が1/2になっている理由はよくわかりませんでした。
glmnet の help で引用されているこちらの論文ではすでに $(1-\alpha)1/2||\beta||^2_2$ として定義されています。
またscikit-learnでも同様にL2ノルムに対しては0.5を乗じているようです（https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html）。
誰か理由を教えてください。
続いて match.call() を用いて引数の指定を正式なものに直します：
this.call = match.call()
これだけだと何を言っているかちょっとわからないと思いますので、以下の例で確認してみましょう：
myfun <- function(abc, def, ghi) {
return(abc + 2*def + 3*ghi)
}
上のように引数として abc 、 def 、 ghi を取る関数を定義します。
このとき R では、引数の指定がない場合には順番通りに入力されます：
> myfun(1, 2, 3)
[1] 14
一部の引数のみ指定がある場合では指定された引数だけがその通りに入力され、残りは順番通りに割り当てられるようです。
> myfun(def = 3, 4, 5)
[1] 25
ところでこの引数の指定は、一意に決まれば指定は省略することができます：
> myfun(d = 3, 4, 5)
[1] 25
一方、例えば以下のような呼び出しでは g から始まる引数が２つあるため一意に決まらず、エラーとなってしまいます。
> myfun2 <- function(abc, def, ghi, gjk) {
+ return(abc + 2*def + 3*ghi + 4*gjk)
+ }
> myfun2(g = 3, 4, 5, 6)
myfun2(g = 3, 4, 5, 6) でエラー: 引数 1 が複数の仮引数に一致します
では match.call を使って関数を呼び出すとどうなるかと言うと：
> match.call(myfun, call("myfun", 1, def = 3, ghi = 5))
myfun(abc = 1, def = 3, ghi = 5)
この通り、各引数に対して何を割り当てたかを得ることができます。
便利ですね。
さらに続いて、 nlambda の指定です。
ここでは $\lambda$ （罰則の大きさ）そのものではなく、検証する $\lambda$ の数（nubmer of lambda）を指定します（デフォルトは100）。
nlam = as.integer(nlambda)
ここからは y 、 x および weight のチェックです：
y = drop(y)
np = dim(x)
if (is.null(np) | (np[2] <= 1))
stop("x should be a matrix with 2 or more columns")
nobs = as.integer(np[1])
if (missing(weights))
weights = rep(1, nobs)
else if (length(weights) != nobs)
stop(paste("number of elements in weights (", length(weights),
") not equal to the number of rows of x (", nobs,
")", sep = ""))
nvars = as.integer(np[2])
dimy = dim(y)
nrowy = ifelse(is.null(dimy), length(y), dimy[1])
if (nrowy != nobs)
stop(paste("number of observations in y (", nrowy, ") not equal to the number of rows of x (",
nobs, ")", sep = ""))
vnames = colnames(x)
if (is.null(vnames))
vnames = paste("V", seq(nvars), sep = "")
y に対する drop ですが、これは length が 1 であるような冗長な次元を落とす関数です。
続いて x の行数が weight や y と合わない場合にエラーを返しています。
以下ではモデルに含める変数や非ゼロとする変数などを指定します
（ nx(=pmax) の方はちょっと理解がアヤシイので help の説明を書いておきます）：
ne = as.integer(dfmax)
nx = as.integer(pmax)
if (missing(exclude))
exclude = integer(0)
次に変数ごとに異なるペナルティを与えるために penalty.factor を指定します。
この数値が lambda に乗じられるため、例えば特定の変数に対して penalty.factor = 0 としておけば罰則を与えないようにすることが可能となります（結果として常にモデルに採用されるようになる）：
if (any(penalty.factor == Inf)) {
exclude = c(exclude, seq(nvars)[penalty.factor == Inf])
exclude = sort(unique(exclude))
}
if (length(exclude) > 0) {
jd = match(exclude, seq(nvars), 0)
if (!all(jd > 0))
stop("Some excluded variables out of range")
penalty.factor[jd] = 1
jd = as.integer(c(length(jd), jd))
}
else jd = as.integer(0)
vp = as.double(penalty.factor)
これはせっかくなので実際にやってみましょう。
冒頭のコードを持ってきて、以下のように lambda を適当に設定してみます。
x <- scale(LifeCycleSavings[, 2:5])
y <- LifeCycleSavings[, 1] - mean(LifeCycleSavings[, 1])
> coef(glmnet(x, y, family = "gaussian", alpha = 1, lambda = 0.3))
5 x 1 sparse Matrix of class "dgCMatrix"
s0
(Intercept) 1.182354e-15
pop15 -1.691002e+00
pop75 .
dpi .
ddpi 9.816514e-01
このとき、2・3番目の変数である pop75 と dpi は 0 と推定されてしまいました。
ここでこれらの変数の penalty.factor を 0 とすると
> coef(glmnet(x, y, family = "gaussian", alpha = 1, lambda = 0.3,
+ penalty.factor = c(1, 0, 0, 1)))
5 x 1 sparse Matrix of class "dgCMatrix"
s0
(Intercept) 9.523943e-16
pop15 -7.827680e-01
pop75 8.127991e-01
dpi -1.560908e-01
ddpi 6.812498e-01
ちゃんと推定されるようになっています。
逆に pop15 の penalty.factor を大きくすると
> coef(glmnet(x, y, family = "gaussian", alpha = 1, lambda = 0.3,
+ penalty.factor = c(2, 0, 0, 1)))
5 x 1 sparse Matrix of class "dgCMatrix"
s0
(Intercept) 7.266786e-16
pop15 .
pop75 1.374655e+00
dpi 2.586151e-02
ddpi 9.300500e-01
このようにモデルから除外されることになります。
さらに penalty.factor = Inf とすると、その変数は exclude として扱われるようになります。
続いて glmnet.control で持っているパラメータを渡します。
internal.parms = glmnet.control()
if (internal.parms$itrace)
trace.it = 1
else {
if (trace.it) {
glmnet.control(itrace = 1)
on.exit(glmnet.control(itrace = 0))
}
}
次に、回帰係数に対する上限・下限を設定します。
なお下限は non-positive 、上限は non-negative しか設定できないようですね。
if (any(lower.limits > 0)) {
stop("Lower limits should be non-positive")
}
if (any(upper.limits < 0)) {
stop("Upper limits should be non-negative")
}
lower.limits[lower.limits == -Inf] = -internal.parms$big
upper.limits[upper.limits == Inf] = internal.parms$big
if (length(lower.limits) < nvars) {
if (length(lower.limits) == 1)
lower.limits = rep(lower.limits, nvars)
else stop("Require length 1 or nvars lower.limits")
}
else lower.limits = lower.limits[seq(nvars)]
if (length(upper.limits) < nvars) {
if (length(upper.limits) == 1)
upper.limits = rep(upper.limits, nvars)
else stop("Require length 1 or nvars upper.limits")
}
else upper.limits = upper.limits[seq(nvars)]
cl = rbind(lower.limits, upper.limits)
if (any(cl == 0)) {
fdev = glmnet.control()$fdev
if (fdev != 0) {
glmnet.control(fdev = 0)
on.exit(glmnet.control(fdev = fdev))
}
}
storage.mode(cl) = "double"
標準化と切片に対する指定です。
標準化の処理そのものは以降の関数の内部で実行されるため、ここでは指定のみを行います。
isd = as.integer(standardize)
intr = as.integer(intercept)
if (!missing(intercept) && family == "cox")
warning("Cox model has no intercept")
jsd = as.integer(standardize.response)
収束を判定する閾値を指定します。
thresh = as.double(thresh)
次に、 lambda に関する指定となりますが、 flmin および ulam の使われ方がよく理解できなかったため、これらの説明は省略します。
なお help にもありますが、通常は lambda には単一の値ではなく、候補となる値のベクトルを与えます。
Avoid supplying a single value for lambda (for predictions after CV use predict() instead).
if (is.null(lambda)) {
if (lambda.min.ratio >= 1)
stop("lambda.min.ratio should be less than 1")
flmin = as.double(lambda.min.ratio)
ulam = double(1)
}
else {
flmin = as.double(1)
if (any(lambda < 0))
stop("lambdas should be non-negative")
ulam = as.double(rev(sort(lambda)))
nlam = as.integer(length(lambda))
}
次に疎行列の指定です。
入力 X が疎行列である場合、dgCMatrix 形式に変換されます。
ここで dgCMatrix とは列方向の志向性を持つ疎行列の形式です。
is.sparse = FALSE
ix = jx = NULL
if (inherits(x, "sparseMatrix")) {
is.sparse = TRUE
x = as(x, "CsparseMatrix")
x = as(x, "dgCMatrix")
ix = as.integer(x@p + 1)
jx = as.integer(x@i + 1)
x = as.double(x@x)
}
ここも、せっかくなので疎行列における数値の格納方法についても見ておきましょう。
以下のように疎行列を作成します：
set.seed(1234)
i <- c(1, 5, 18)
j <- c(4, 13, 19)
n <- rnorm(3)
m <- matrix(0, 20, 20)
for (k in 1:length(n)) {
m[i[k], j[k]] <- n[k]
}
s_m <- as(m, "dgCMatrix")
ここで s_m は行列 m を疎行列として扱ったものです。
str() で確認すると、 s_m には
@ i ：非ゼロの要素の入っていた行番号( 0-index であることに注意)@ p ：各列における非ゼロの要素の個数を積み上げたもの@ Dim ：行列の次元@ Dimnames ：行列の各次元の名前@ x ：非ゼロの要素の数値@ factors ：（これはちょっとわかりませんでした）
が格納されています。
> str(s_m)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:3] 0 4 17
..@ p : int [1:21] 0 0 0 0 1 1 1 1 1 1 ...
..@ Dim : int [1:2] 20 20
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ x : num [1:3] -1.207 0.277 1.084
..@ factors : list()
ここで @ i には非ゼロである各要素の行番号が入るため行列 m を作ったときの行番号の指定 i に対応しますが、0-index であるため数字は1つずつ小さくなっています。
> print(i- 1)
[1] 0 4 17
> print(s_m@i)
[1] 0 4 17
ちょっとわかりにくいのが @ p で、ここには各列における非ゼロの要素の個数の累積が格納され、列数に対応します（ただし最初に 0 が追加されるため、列数 + 1 の長さになります）。
今回の例では行列の列数が 20 なので、length が 21 となります。
> length(s_m@p)
[1] 21
このベクトルには非ゼロの要素の個数の累積が入っているため、差分を取ると元の行列で非ゼロの要素が入っていた列を得ることができます。
> diff(s_m@p)
[1] 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0
列番号を指定した j と比較してみましょう：
> which(diff(s_m@p) == 1)
[1] 4 13 19
> j
[1] 4 13 19
合っていますね。
続く処理では、 ix には各列における非ゼロの要素の累積個数(+1)を
、 jx には行番号を代入しています。
また x には元の疎行列における非ゼロの要素の値そのものをベクトルとして入力しており、説明変数の行列が疎行列であった場合、この時点で行列ではなくベクトルとして扱われることになります。
次に、プログレスバーの指定です（出せるんですね）。
if (trace.it) {
if (relax)
cat("Training Fit\n")
pb <- createPB(min = 0, max = nlam, initial = 0, style = 3)
}
そして最後に最適化の手法についての指定です。
family が `binomial または multinomial の場合、 glmnet の引数である type.logistic および type.multinomial が評価され、（後の工程で）それに応じて呼ばれる関数が変わります。
具体的には lognet2m 、 lognetn および multlognetn のどれが選ばれるかが決まります。
これは別の機会に解説します（予定です）。
kopt = switch(match.arg(type.logistic), Newton = 0, modified.Newton = 1)
if (family == "multinomial") {
type.multinomial = match.arg(type.multinomial)
if (type.multinomial == "grouped")
kopt = 2
}
kopt = as.integer(kopt)
最初の方で family のチェックに使われ、ここでも使われている match.arg ですが、せっかくなので挙動を確認しておきましょう：
myfun <- function(a = "aaa", type.logistic = c("Newton", "modified.Newton")) {
kopt <- switch(match.arg(type.logistic), Newton = 0, modified.Newton = 1)
kopt
}
上のような関数を定義し、以下のように呼び出すと、結果はそれぞれ 0, 0, 1 となります。
> myfun()
[1] 0
> myfun(type.logistic = "Newton")
[1] 0
> myfun(type.logistic = "modified.Newton")
[1] 1
2. フィッティング
以上でパラメータの設定や前処理が終わりましたので次はフィッティングです。
といってもここでは family に応じて呼び出す関数を変えているだけなので、詳細は一旦スキップしましょう。
fit = switch(family,
gaussian = elnet(x, is.sparse, ix, jx,
y, weights, offset, type.gaussian, alpha, nobs, nvars,
jd, vp, cl, ne, nx, nlam, flmin, ulam, thresh, isd, intr,
vnames, maxit),
poisson = fishnet(x, is.sparse, ix, jx,
y, weights, offset, alpha, nobs, nvars, jd, vp, cl, ne,
nx, nlam, flmin, ulam, thresh, isd, intr, vnames, maxit),
binomial = lognet(x, is.sparse, ix, jx, y, weights, offset,
alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, flmin,
ulam, thresh, isd, intr, vnames, maxit, kopt, family),
multinomial = lognet(x, is.sparse, ix, jx, y, weights,
offset, alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam,
flmin, ulam, thresh, isd, intr, vnames, maxit, kopt,
family),
cox = coxnet(x, is.sparse, ix, jx, y, weights,
offset, alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam,
flmin, ulam, thresh, isd, vnames, maxit),
mgaussian = mrelnet(x,
is.sparse, ix, jx, y, weights, offset, alpha, nobs,
nvars, jd, vp, cl, ne, nx, nlam, flmin, ulam, thresh,
isd, jsd, intr, vnames, maxit))
if (trace.it) {
utils::setTxtProgressBar(pb, nlam)
close(pb)
}
なおここでそれぞれの関数に渡されている引数を比較すると以下のようになります（一部はよくわかりませんでした）：
| 引数 |
説明 |
elnet |
fishnet |
lognet |
coxnet |
mrelnet |
| x |
説明変数の行列 |
〇 |
〇 |
〇 |
〇 |
〇 |
| is.sparse |
疎行列であるかの指定 |
〇 |
〇 |
〇 |
〇 |
〇 |
| ix |
疎行列における非ゼロの要素の累積個数 |
〇 |
〇 |
〇 |
〇 |
〇 |
| jx |
疎行列における非ゼロの要素の行番号 |
〇 |
〇 |
〇 |
〇 |
〇 |
| y |
目的変数の行列 |
〇 |
〇 |
〇 |
〇 |
〇 |
| weights |
観測値に対する重み |
〇 |
〇 |
〇 |
〇 |
〇 |
| offset |
オフセット |
〇 |
〇 |
〇 |
〇 |
〇 |
| type.gaussian |
1:covariance, 2:naïve |
〇 |
- |
- |
- |
- |
| alpha |
L1とL2に対する重みの調整パラメータ |
〇 |
〇 |
〇 |
〇 |
〇 |
| nobs |
レコード数 |
〇 |
〇 |
〇 |
〇 |
〇 |
| nvars |
説明変数の数 |
〇 |
〇 |
〇 |
〇 |
〇 |
| jd |
? |
〇 |
〇 |
〇 |
〇 |
〇 |
| vp |
各変数に対する罰則の重み（penalty.factor） |
〇 |
〇 |
〇 |
〇 |
〇 |
| cl |
? |
〇 |
〇 |
〇 |
〇 |
〇 |
| ne |
モデルに含まれる変数の上限。ne = dfmax = nvars + 1 |
〇 |
〇 |
〇 |
〇 |
〇 |
| nx |
非ゼロとする変数の個数の上限？ |
〇 |
〇 |
〇 |
〇 |
〇 |
| nlam |
lambdaの数 |
〇 |
〇 |
〇 |
〇 |
〇 |
| flmin |
? |
〇 |
〇 |
〇 |
〇 |
〇 |
| ulam |
? |
〇 |
〇 |
〇 |
〇 |
〇 |
| thresh |
収束判定の閾値 |
〇 |
〇 |
〇 |
〇 |
〇 |
| isd |
standardizeするかの指定 |
〇 |
〇 |
〇 |
〇 |
〇 |
| jsd |
? |
- |
- |
- |
- |
〇 |
| intr |
切片（Intercept）を含めるかの指定 |
〇 |
〇 |
〇 |
- |
〇 |
| vnames |
変数名 |
〇 |
〇 |
〇 |
〇 |
〇 |
| maxit |
反復回数の上限 |
〇 |
〇 |
〇 |
〇 |
〇 |
| kopt |
最適化の手法 |
- |
- |
〇 |
- |
- |
| family |
family |
- |
- |
〇 |
- |
- |
3. 後処理
最後に後処理です。
if (is.null(lambda))
fit$lambda = fix.lam(fit$lambda)
fit$call = this.call
fit$nobs = nobs
class(fit) = c(class(fit), "glmnet")
if (relax)
relax.glmnet(fit, x = x, y = y, weights = weights, offset = offset,
lower.limits = lower.limits, upper.limits = upper.limits,
check.args = FALSE, ...)
else fit
この後処理で目立つ工程としては relax の部分でしょう。
ここで relax は help によると、
If relax=TRUE a duplicate sequence of models is produced, where each active set in the elastic-net path is refit without regularization. The result of this is a matching "glmnet" object which is stored on the original object in a component named "relaxed", and is part of the glmnet output.
ということで、glmnet によって変数選択された結果を用いて、罰則なしで再度フィッティングを行うオプションのようです。
これも実際にやってみるのが早いと思いますので、以下のように実行してみます：
lasso_02 <- glmnet(x, y, family = "gaussian", relax = T)
すると、先程の結果（ lasso ）に、 lasso_02$relaxed という結果が追加されていることがわかりますが、内容は lasso とほとんど同じです。
> str(lasso)
List of 12
$ a0 : Named num [1:68] 6.11e-16 6.71e-16 7.26e-16 7.76e-16 8.22e-16 ...
..- attr(*, "names")= chr [1:68] "s0" "s1" "s2" "s3" ...
$ beta :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. ..@ i : int [1:216] 0 0 0 0 0 3 0 3 0 3 ...
.. ..@ p : int [1:69] 0 0 1 2 3 4 6 8 10 12 ...
.. ..@ Dim : int [1:2] 4 68
.. ..@ Dimnames:List of 2
.. .. ..$ : chr [1:4] "pop15" "pop75" "dpi" "ddpi"
.. .. ..$ : chr [1:68] "s0" "s1" "s2" "s3" ...
.. ..@ x : num [1:216] -0.181 -0.347 -0.497 -0.634 -0.757 ...
.. ..@ factors : list()
$ df : int [1:68] 0 1 1 1 1 2 2 2 2 2 ...
$ dim : int [1:2] 4 68
$ lambda : num [1:68] 2.02 1.84 1.68 1.53 1.39 ...
$ dev.ratio: num [1:68] 0 0.0352 0.0645 0.0888 0.1089 ...
$ nulldev : num 984
$ npasses : int 562
$ jerr : int 0
$ offset : logi FALSE
$ call : language glmnet(x = x, y = y, family = "gaussian", alpha = 1)
$ nobs : int 50
- attr(*, "class")= chr [1:2] "elnet" "glmnet"
> str(lasso_02$relaxed)
List of 12
$ a0 : Named num [1:68] 6.11e-16 1.29e-15 1.29e-15 1.29e-15 1.29e-15 ...
..- attr(*, "names")= chr [1:68] "s0" "s1" "s2" "s3" ...
$ beta :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. ..@ i : int [1:216] 0 0 0 0 0 3 0 3 0 3 ...
.. ..@ p : int [1:69] 0 0 1 2 3 4 6 8 10 12 ...
.. ..@ Dim : int [1:2] 4 68
.. ..@ Dimnames:List of 2
.. .. ..$ : chr [1:4] "pop15" "pop75" "dpi" "ddpi"
.. .. ..$ : chr [1:68] "s0" "s1" "s2" "s3" ...
.. ..@ x : num [1:216] -2.04 -2.04 -2.04 -2.04 -1.98 ...
.. ..@ factors : list()
$ df : int [1:68] 0 1 1 1 1 2 2 2 2 2 ...
$ dim : int [1:2] 4 68
$ lambda : num [1:68] 2.02 1.84 1.68 1.53 1.39 ...
$ dev.ratio: num [1:68] 0 0.208 0.208 0.208 0.208 ...
$ nulldev : num 984
$ npasses : int 562
$ jerr : int 0
$ offset : logi FALSE
$ call : language glmnet(x = x, y = y, family = "gaussian", relax = T)
$ nobs : int 50
- attr(*, "class")= chr [1:2] "elnet" "glmnet"
ここで lasso_02$relaxed の中身を少し見てみると、例えば beta には以下のような数値が入っています。
> lasso_02$relaxed$beta[, 1:6]
4 x 6 sparse Matrix of class "dgCMatrix"
s0 s1 s2 s3 s4 s5
pop15 . -2.040996 -2.040996 -2.040996 -2.040996 -1.980216
pop75 . . . . . .
dpi . . . . . .
ddpi . . . . . 1.270865
これは何かと言うと、少しずつ罰則の重みを変えたことで変数が選択された状態で通常の線形回帰を当てはめたときの回帰係数となっています。
例えば lasso_02$relaxed$beta[, 6] には、変数として選択された pop15 と ddpi それぞれの回帰係数が入っています。
実際に lm の結果と一致するか見てみましょう：
> coef(lm(y ~ x[, c(1, 4)]))
(Intercept) x[, c(1, 4)]pop15 x[, c(1, 4)]ddpi
1.364331e-15 -1.980216e+00 1.270865e+00
合っていますね。
ところで切片の推定値が入っている lasso_02$relaxed$a0 の値は少し異なるようです：
> lasso_02$relaxed$a0[6]
s5
1.28119e-15
なんでやろか。。
もしかしたら標準化の違いかとも思いましたがそれでもないようで、この理由はわかりませんでした。
lasso_03 <- glmnet(x, y, family = "gaussian", relax = T, standardize = F)
> lasso_03$relaxed$a0[6]
s5
1.28119e-15
glmnet() の実装は以上となります。
次回はフィッティングの部分で呼ばれている elnet を詳しく見ていきましょう。
なお gam のときとは違い、 glmnet では library をインストールしてもソースコードは付いてきませんでしたので、こちらを参考に fortran のソースコードを取得しました。
ではまた次回。
