R biglasso 结果与 hdm 或 glmnet 不匹配
R biglasso results don't match hdm or glmnet
我一直在尝试使用 R 包 'biglasso' 来处理高维数据。但是,我得到的结果与我从 'hdm' 或 'glmnet.LASSO 函数得到的结果不匹配。 biglasso 的文档也很差。
在下面的示例中,hdm 和 glmnet 的结果非常接近但不准确,这是预期的。但是,biglasso 不会删除 'share' 变量。我已经尝试了所有不同的屏幕设置,但没有任何区别。关于如何使 biglasso 与其他人更一致的任何想法?谢谢!
编辑:对于给定的 lambda 值,结果非常相似。但是每种方法似乎 select 一个不同的 lambda.. 这对于 hdm 来说是有意义的,因为它旨在用于因果推理并且不关心样本外预测。 hdm 使用与 Belloni 等人不同的 objective 函数。 (2012),但我不确定为什么 cv.biglasso 和 cv.glmnet 会有如此大的不同。如果我 运行 biglasso 没有筛选规则,他们应该最大化相同的 objective 函数,只是在 CV 折叠中有随机差异,不是吗?
编辑 2:我编辑了下面的代码以包含 F.Privé 的代码,使 glmnet 使用类似于 biglasso 的算法,以及一些额外的代码使 biglasso 模仿 glmnet。
##########
## PREP ##
##########
## Load required libraries
library(hdm)
library(biglasso)
library(glmnet)
## Read automobile dataset
data(BLP)
df <- BLP[[1]]
## Extract outcome
Y <- scale(df$mpg)
## Rescale variables
df$price <- scale(df$price)
df$mpd <- scale(df$mpd)
df$space <- scale(df$space)
df$hpwt <- scale(df$hpwt)
df$outshr <- scale(df$outshr)
## Limit to variables I want
df <- df[,names(df) %in% c("price","mpd","space","hpwt","share","outshr","air")]
## Convert to matrix
df.mat <- data.matrix(df)
df.bm <- as.big.matrix(df.mat)
#########
## HDM ##
#########
## Set seed for reproducibility
set.seed(1233)
## Run LASSO
fit.hdm <- rlasso(x=df.mat, y=Y, post=FALSE, intercept=TRUE)
## Check results
coef(fit.hdm)
############
## GLMNET ##
############
## Set seed for reproducibility
set.seed(1233)
## LASSO with 10-fold cross-validation
fit.glmnet <- cv.glmnet(df.mat, Y, alpha=1, family="gaussian")
## Check default results
coef(fit.glmnet)
## Try to mimic results of biglasso
coef(fit.glmnet, s = "lambda.min")
##############
## BIGLASSO ##
##############
## LASSO with 10-fold cross-validation
fit.bl <- cv.biglasso(df.bm, Y, penalty="lasso", eval.metric="default",
family="gaussian", screen="None",
seed=1233, nfolds=10)
## Check default results
coef(fit.bl)
## Try to mimic results of glmnet
## Calculate threshold for CV error (minimum + 1 standard error)
thresh <- min(fit.bl$cve) + sd(fit.bl$cve)/sqrt(100)
## Identify highest lambda with CVE at or below threshold
max.lambda <- max(fit.bl$lambda[fit.bl$cve <= thresh])
## Check results for the given lambda
coef(fit.bl$fit)[,which(fit.bl$fit$lambda==max.lambda)]
选择CV后的"best" lambda基本上有两种方式:
最小化CV误差的那个(默认为{biglasso})
最简约(最高 lambda)且 CV 误差低于最小标准误差 + 1(默认为 {glmnet})的那个。
尝试coef(fit.glmnet, s = "lambda.min")
使用最小值。
此外,为了确保可重复性,请尝试设置 CV 折叠而不是一些种子。 glmnet()
中有参数foldid
,biglasso()
中有参数cv.ind
。
我一直在尝试使用 R 包 'biglasso' 来处理高维数据。但是,我得到的结果与我从 'hdm' 或 'glmnet.LASSO 函数得到的结果不匹配。 biglasso 的文档也很差。
在下面的示例中,hdm 和 glmnet 的结果非常接近但不准确,这是预期的。但是,biglasso 不会删除 'share' 变量。我已经尝试了所有不同的屏幕设置,但没有任何区别。关于如何使 biglasso 与其他人更一致的任何想法?谢谢!
编辑:对于给定的 lambda 值,结果非常相似。但是每种方法似乎 select 一个不同的 lambda.. 这对于 hdm 来说是有意义的,因为它旨在用于因果推理并且不关心样本外预测。 hdm 使用与 Belloni 等人不同的 objective 函数。 (2012),但我不确定为什么 cv.biglasso 和 cv.glmnet 会有如此大的不同。如果我 运行 biglasso 没有筛选规则,他们应该最大化相同的 objective 函数,只是在 CV 折叠中有随机差异,不是吗?
编辑 2:我编辑了下面的代码以包含 F.Privé 的代码,使 glmnet 使用类似于 biglasso 的算法,以及一些额外的代码使 biglasso 模仿 glmnet。
##########
## PREP ##
##########
## Load required libraries
library(hdm)
library(biglasso)
library(glmnet)
## Read automobile dataset
data(BLP)
df <- BLP[[1]]
## Extract outcome
Y <- scale(df$mpg)
## Rescale variables
df$price <- scale(df$price)
df$mpd <- scale(df$mpd)
df$space <- scale(df$space)
df$hpwt <- scale(df$hpwt)
df$outshr <- scale(df$outshr)
## Limit to variables I want
df <- df[,names(df) %in% c("price","mpd","space","hpwt","share","outshr","air")]
## Convert to matrix
df.mat <- data.matrix(df)
df.bm <- as.big.matrix(df.mat)
#########
## HDM ##
#########
## Set seed for reproducibility
set.seed(1233)
## Run LASSO
fit.hdm <- rlasso(x=df.mat, y=Y, post=FALSE, intercept=TRUE)
## Check results
coef(fit.hdm)
############
## GLMNET ##
############
## Set seed for reproducibility
set.seed(1233)
## LASSO with 10-fold cross-validation
fit.glmnet <- cv.glmnet(df.mat, Y, alpha=1, family="gaussian")
## Check default results
coef(fit.glmnet)
## Try to mimic results of biglasso
coef(fit.glmnet, s = "lambda.min")
##############
## BIGLASSO ##
##############
## LASSO with 10-fold cross-validation
fit.bl <- cv.biglasso(df.bm, Y, penalty="lasso", eval.metric="default",
family="gaussian", screen="None",
seed=1233, nfolds=10)
## Check default results
coef(fit.bl)
## Try to mimic results of glmnet
## Calculate threshold for CV error (minimum + 1 standard error)
thresh <- min(fit.bl$cve) + sd(fit.bl$cve)/sqrt(100)
## Identify highest lambda with CVE at or below threshold
max.lambda <- max(fit.bl$lambda[fit.bl$cve <= thresh])
## Check results for the given lambda
coef(fit.bl$fit)[,which(fit.bl$fit$lambda==max.lambda)]
选择CV后的"best" lambda基本上有两种方式:
最小化CV误差的那个(默认为{biglasso})
最简约(最高 lambda)且 CV 误差低于最小标准误差 + 1(默认为 {glmnet})的那个。
尝试coef(fit.glmnet, s = "lambda.min")
使用最小值。
此外,为了确保可重复性,请尝试设置 CV 折叠而不是一些种子。 glmnet()
中有参数foldid
,biglasso()
中有参数cv.ind
。