glmnetパッケージのv4.0が出たので内容のまとめ
はじめに
ここでは、glmnetパッケージのバージョン4.0が出たので、その内容のまとめです。glmnetパッケージは正則化回帰を使用する際に使用していたのですが、statsパッケージのglm()に含まれるリンク関数や確率分布が使用できるようになったらしいので、そのまとめです。下記のvignettesを参考にしています。
glmnetパッケージのv4.0の変更点
正則化回帰とは、基本的にがL1正則化(Lasso)、L2正則化(Ridge)、エラスティックネット回帰などの回帰のことで、下記のような形でチューニングパラメタ\( \alpha \)とペナルティの強さをコントロールする\( \lambda \)が加った回帰分析のこと。こちらのブログが非常にわかりやすかったので、詳細はこちらを参照ください。
さて、v4.0以降のglmnetパッケージで何が変わったのかというと、statsパッケージのglm()に含まれるリンク関数や確率分布が使用できるようになり、これらの関数もplot()やpredict()が使えるようになったとのことです。また、確率分布の指定方法が変わったみたいなのです。文字で見るよりも見て動かしたほうが、変化を理解するためには早いので、早速動かしてみます。
GLMをやってみる
バージョンの確認とパッケージを済ましておきます。
# install.packages("glmnet")
packageVersion("glmnet")
[1] ‘4.0.2’
library(glmnet)
vignettesに従って、ダミーデータを生成しておきます。
set.seed(1989)
x <- matrix(rnorm(500), ncol = 5)
y <- rowSums(x[, 1:2]) + rnorm(100)
head(x)
[,1] [,2] [,3] [,4] [,5]
[1,] 1.1025783 1.73717481 0.7109263564 -0.6539329 0.2105306
[2,] 1.1178965 -0.38112662 1.7013464697 -0.1632433 -1.4645134
[3,] -1.8181019 -0.76206955 -0.4826317551 0.5241455 -0.3993483
[4,] -0.1944140 1.10435999 -0.0515429445 -1.0147559 0.5760584
[5,] -0.6131956 -0.64528426 -0.8655069038 0.6377163 -2.0518750
[6,] -0.3462673 0.06856831 0.0008600987 0.5825773 -1.0980416
head(y)
[1] 1.2335521 0.7436334 -2.4825280 2.3509169 -2.8536794
[6] 1.1825112
v4.0以降のglmnetパッケージでは、関数の呼び出し方が異なります。vignettesにもありますが、古い指定方法のほうが処理が速いようです。
# old method glmnet(x, y, family = "gaussian") # new method # glmnet(x, y, family = gaussian(link = "identity"))でももちろんOK glmnet(x, y, family = gaussian())
下記のvignettesの引用にも有るように、ここらへん変わったポイントらしくv4.0以降では、確率分布がfamilyオブジェクトであれば、どのような確率分布に対してもペナルティ付きGLMをフィットすることができます。
From v4.0 onwards, glmnet can fit penalized GLMs for any family as long as the family can be expressed as a family object.
familyオブジェクトというのはこういうこと。
family <- list(binomial(),
gaussian(),
Gamma(),
inverse.gaussian(),
poisson(),
quasi(),
quasibinomial(),
quasipoisson())
purrr::map_chr(.x = family, .f= function(x){class(rlang::enexpr(x))})
[1] "family" "family" "family" "family" "family" "family" "family" "family"
purrr::map(.x = family, .f= function(x){rlang::enexpr(x)})
[[1]]
Family: binomial
Link function: logit
[[2]]
Family: gaussian
Link function: identity
[[3]]
Family: Gamma
Link function: inverse
[[4]]
Family: inverse.gaussian
Link function: 1/mu^2
[[5]]
Family: poisson
Link function: log
[[6]]
Family: quasi
Link function: identity
[[7]]
Family: quasibinomial
Link function: logit
[[8]]
Family: quasipoisson
Link function: log
例えば、二項分布、ポアソン分布を使った回帰の場合は下記のように指定すれば動きます。statsパッケージのglm()と同じですね。
biny <- ifelse(y > 0, 1, 0) # binary data cnty <- ceiling(exp(y)) # count data # fitting binomial GLMs the new way newfit <- glmnet(x, biny, family = binomial()) # fitting Poisson GLMs the new way newfit <- glmnet(x, cnty, family = poisson())
例えば、エラスティックネットを用いたプロビット回帰、準ポアソン回帰、補対数-対数回帰、負の2項回帰は、以下のコードのようにかけます。
newfit <- glmnet(x, biny, family = binomial(link = "probit")) newfit <- glmnet(x, cnty, family = quasipoisson()) newfit <- glmnet(x, biny, family = binomial(link = "cloglog”)) library(MASS) newfit <- glmnet(x, cnty, family = negative.binomial(theta = 5))
ちなみに個人的な記憶に基づきますが、glmnet()はマトリックス型で指定する使い方しかできなかったような気がしますが…データフレームでも使えるようです。ちなみにglmnetUtilsパッケージを使えばy ~ X, data = dfという形で実行できます。
fit_mat <- glmnet(xx, y, family = gaussian(link = "identity"), alpha = 1, lambda = 0.5)
fit_mat$beta
5 x 1 sparse Matrix of class "dgCMatrix"
s0
V1 0.7817021
V2 0.4759707
V3 .
V4 .
V5 .
glmnet()の引数はこうなっています。family、alpha、lambdaを指定することで、様々な正則化GLMが実行できます。
args(glmnet)
function(
x,
y,
family = c("gaussian", "binomial", "poisson", "multinomial", "cox", "mgaussian"),
weights = NULL,
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 = NULL,
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,
...)
また、フィットさせたモデルのクラスは下記のようになりますが、重要なことは、これらはglmnetクラスを継承しているので、plot(),predict(),coef(),print()などのS3メソッドはそのまま使えます。
class(newfit) [1] "glmnetfit" "glmnet" # 下記は利用可能 fit <- glmnet(x = x, y = y, family = gaussian(link = "identity"), lambda = 0.01, alpha = 1) plot(fit) predict(fit, newx = x) coef(fit) print(fit)
また、broomパッケージのtidy()やglance()は使えます。
fit <- glmnet(x = x, y = y, family = gaussian(link = "identity"), lambda = 0.01, alpha = 1)
broom::tidy(fit)
# A tibble: 6 x 5
term step estimate lambda dev.ratio
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 1 -0.0695 0.01 0.669
2 V1 1 0.982 0.01 0.669
3 V2 1 0.978 0.01 0.669
4 V3 1 -0.0488 0.01 0.669
5 V4 1 -0.0964 0.01 0.669
6 V5 1 0.0589 0.01 0.669
broom::glance(fit)
# A tibble: 1 x 2
nulldev npasses
<dbl> <int>
1 284. 14
broom::augment(fit)
エラー: No augment method for objects of class glmnetfit
パラメタチューニング:cv.glmnet()
ここではcv.glmnet()も新しい表記で使用できるのか確認しておきます。まずは適当にダミーデータを作成します。ここでは、目的変数と関係が本当に関係がある変数が10個、ノイズとなる変数が90個、無関係な変数が900個でサンプルサイズが1000というデータを作成しています。
library(HCmodelSets)
library(tidyverse)
true_n <- 10
noise_n <- 90
sample_size <- 1000
sample_cols <- true_n + noise_n + 900
# GGP:Data generating process
dgp <- DGP(
s = true_n, # 目的変数と相関のある説明変数の数
a = noise_n, # 説明変数と相関のあるノイズ変数の数
sigStrength = 0.8, # 目的変数と説明変数の相関の強さ
rho = 0.7, # 説明変数とノイズ変数相関の強さ
n = sample_size, # サンプルサイズ
noise = 1, # 真の回帰直線の観測値の分散
var = 1, # 潜在的な説明変数の分散
d = sample_cols, # 潜在的な説明変数の数
type.response = "N", # 目的変数を正規分布から生成
intercept = 0, # 回帰直線の切片
DGP.seed = 1989) # 乱数シード
Y <- as_tibble(dgp$Y)
X <- as_tibble(dgp$X)
tmp <- X %>% bind_cols(Y)
names(tmp) <- c(paste("X", 1:sample_cols, sep = "_"), "Y")
true_val <- names(tmp[,dgp$TRUE.idx])
df <- tmp %>%
rename_at(vars(true_val), list( ~ paste0(., "_t"))) %>%
select(Y, everything())
dim(df)
[1] 1000 1001
y <- df %>% select(Y) %>% as.matrix()
X <- df %>% select(-Y) %>% as.matrix()
Lasso回帰を行うにはalpha=1に指定して実行します。正則化パラメタ\( \lambda \)を指定する必要がありますが、それを決めるためにクロスバリデーションを使って推定値と観測値の平均二乗誤差(他にも設定可)が最小となるように\( \lambda \)を決定します。本来であれば、訓練、検証、テストの3つに分けて分析するのが一般的かと思いますが、ここでは動きの確認なので、すっ飛ばしていますし、本来は説明変数を標準化してからフィットしますが、それもここでは省略します。クロスバリデーションの結果、平均二乗誤差を最小にする\( \lambda \)は 0.100665 であることがわかります。
# plot(lasso_fit_cv) lasso_fit_cv$lambda.min [1] 0.100665
最小の\( \lambda \)を使ってモデルをフィットし直します。
lasso_fit <- glmnet(x = X, y = y, family = gaussian(link = "identity"), alpha = 1, lambda = lasso_fit_cv$lambda.min) library(broom) lasso_fit %>% tidy() %>% arrange(desc(estimate)) %>% print(n = 30) # A tibble: 22 x 5 term step estimate lambda dev.ratio <chr> <dbl> <dbl> <dbl> <dbl> 1 X_572_t 1 0.883 0.101 0.979 2 X_219_t 1 0.876 0.101 0.979 3 X_828_t 1 0.786 0.101 0.979 4 X_157_t 1 0.781 0.101 0.979 5 X_660_t 1 0.718 0.101 0.979 6 X_590_t 1 0.713 0.101 0.979 7 X_643_t 1 0.679 0.101 0.979 8 X_294_t 1 0.678 0.101 0.979 9 X_13_t 1 0.677 0.101 0.979 10 X_331_t 1 0.669 0.101 0.979 11 X_577 1 0.152 0.101 0.979 12 X_149 1 0.0551 0.101 0.979 13 X_203 1 0.0416 0.101 0.979 14 X_766 1 0.0403 0.101 0.979 15 X_79 1 0.0393 0.101 0.979 16 X_722 1 0.0255 0.101 0.979 17 X_109 1 0.0252 0.101 0.979 18 (Intercept) 1 0.0228 0.101 0.979 19 X_341 1 0.0103 0.101 0.979 20 X_262 1 0.00590 0.101 0.979 21 X_171 1 0.00255 0.101 0.979 22 X_893 1 -0.0207 0.101 0.979
22個しか変数が表示されないのはglmnet()のdgCMatrixクラスの仕様です。詳細は調べてないので不明ですが、意味のある値のインデックス@iと係数@xしか持っていないので、tidy()で値を引っ張る際に、係数が0になったものが省かれて22個のみ出力されるということになります。
vals <- coef(lasso_fit) str(vals) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:22] 0 13 79 109 149 157 171 203 219 262 ... ..@ p : int [1:2] 0 22 ..@ Dim : int [1:2] 1001 1 ..@ Dimnames:List of 2 .. ..$ : chr [1:1001] "(Intercept)" "X_1" "X_2" "X_3" ... .. ..$ : chr "s0" ..@ x : num [1:22] 0.0228 0.6765 0.0393 0.0252 0.0551 ... ..@ factors : list()
