変数の順序は、glmnetの推定係数を変更します

Question

私はRのglmnetパッケージを使って作業しています。古い分類子を再現しようとすると問題が発生しました。説明変数が置換されていれば（例えば、逆の順序で）、cv.glmnetからの結果の係数は、不透明な設計行列を使用する係数と等しくなりません。

は、例えば、以下のデータを考慮してください

は

library(glmnet) 
set.seed(1) 

#Set initial parameters 
n <- 100 
p <- 1000 

#Simulate data 
x <- matrix(rnorm(n * p), nrow = n, ncol = p) 
colnames(x) <- as.character(1:p) 
beta <- rnorm(n = p, mean = 2, sd = 2) 
beta[rbinom(p, size = 1, prob = 0.5) == 0] <- 0 
y <- x %*% beta + rnorm(100, sd = 0.1) 

そして設計行列X及びXの置換バージョンの両方についてLASSOペナルティ（アルファ= 1）でglmnetを実行します。

#Set parameters for cross validation with cv.glmnet 
lambda <- exp(seq(-1, 1, length.out = 100)) 
alpha <- 1 
foldid <- rep(1:10, each = 10) 

#Run cross validation 
fit <- cv.glmnet(x = x, y = y, family = "gaussian", alpha = alpha, 
       lambda = lambda, 
       foldid = foldid) 

#Save coefficients 
coef1 <- as.matrix(coef(fit, s = "lambda.min")) 

#Run cross validation with rearranged design matrix 
order <- ncol(x):1 
fit2 <- cv.glmnet(x = x[,order], y = y, family = "gaussian", alpha = alpha, 
        lambda = lambda, 
        foldid = foldid) 

#Save coefficients 
coef2 <- as.matrix(coef(fit2, s = "lambda.min")) 
coef2 <- coef2[rownames(coef1),] 

次に、係数、相互検証エラー、および線形予測子を比較します。

> summary(coef2 - coef1) 
     1    
Min. :-0.2738963 
1st Qu.: 0.0000000 
Median : 0.0000000 
Mean : 0.0003739 
3rd Qu.: 0.0000000 
Max. : 0.3643040 

> min(fit$cvm) 
[1] 4584.373 
> min(fit2$cvm) 
[1] 4596.626 

> summary(cbind(1,x) %*% coef2 - cbind(1, x) %*% coef1) 
     1   
Min. :-0.5100 
1st Qu.:-0.1613 
Median : 0.0210 
Mean : 0.0000 
3rd Qu.: 0.1333 
Max. : 0.6139 

3つのすべての対策について、変数の順序のみが変更されているのに対し、モデルの違いがわかります。誰でもこれを説明できますか？

Answer 1

私は、glmnetが座標降下を使用していると考えています。ここで、変数は損失関数を最小限に抑えるために反復されます。その場合の変数の順序は、反復の順序を変更し、損失関数を最小にするために取られる経路を変更する。

Answer 2

Glmnetは、座標降下を介してLASSOの正則化パスを計算します（例：Trevor Hastieのスライド15：http://web.stanford.edu/~hastie/TALKS/glmnet.pdfを参照）。アルゴリズムが係数を反復処理するので、変数の順序が取られるパスに影響を与えます。収束しきい値および反復の最大回数に応じて、これは係数の最終値の差につながります。あなたの例の場合は、変更してみてください。

fit <- cv.glmnet(x = x, y = y, family = "gaussian", alpha = alpha, 
       lambda = lambda, 
       foldid = foldid) 

fit <- cv.glmnet(x = x, y = y, family = "gaussian", alpha = alpha, 
       lambda = lambda, 
       foldid = foldid, standardize=TRUE, thresh=1e-20, maxit=10^6) 

fit2と同じです。これには計算に1分ほどかかることがありますが、その差はごくわずかになります。

> summary(coef2 - coef1) 
     1    
Min. :-2.038e-08 
1st Qu.: 0.000e+00 
Median : 0.000e+00 
Mean : 1.050e-10 
3rd Qu.: 0.000e+00 
Max. : 3.028e-08 
> 
> min(fit$cvm) 
[1] 4598.242 
> 
> min(fit2$cvm) 
[1] 4598.242 
> 
> summary(cbind(1,x) %*% coef2 - cbind(1, x) %*% coef1) 
     1    
Min. :-5.175e-08 
1st Qu.:-1.457e-08 
Median :-2.959e-10 
Mean : 0.000e+00 
3rd Qu.: 1.503e-08 
Max. : 5.555e-08