Почему регрессия glmnet ridge дает мне другой ответ, чем ручной расчет?

Я использую glmnet для расчета оценок регрессии гребня. Я получил некоторые результаты, которые сделали меня подозрительным в том, что glmnet действительно делает то, что я думаю, что делает. Чтобы проверить это, я написал простой R-скрипт, в котором я сравниваю результат регрессии гребня, выполненного execute, и результат в glmnet, разница значительна:

n    <- 1000
p.   <-  100
X.   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

beta1 <- solve(t(X)%*%X+5*diag(p),t(X)%*%Y)
beta2 <- glmnet(X,Y, alpha=0, lambda=10, intercept=FALSE, standardize=FALSE, 
                family="gaussian")$beta@x
beta1-beta2

Норма разницы обычно составляет около 20, что не может быть связано с численно различными алгоритмами, я должен делать что-то не так. Какие настройки мне нужно установить glmnet, чтобы получить тот же результат, что и для ridge?

r ridge-regression glmnet Джон
источник

Вы видели этот вопрос ?

cdeterman

Да, но я все еще не получаю тот же результат, используя нормализацию.

Джон

Не могли бы вы опубликовать свой код тогда?

Shadowtalker

У меня просто была такая же проблема! a = data.frame (a = дрожание (1:10), b = дрожание (1:10), c = дрожание (1:10), d = дрожание (1:10), e = дрожание (1:10) , f = дрожание (1:10), g = выборка (дрожание (1:10)), y = сек (10 100,10)); коэф (lm.ridge (y ~ a + b + c + d + e + f + g, a, lambda = 2,57)); coef (glmnet (as.matrix (a [, 1: 7]), a $ y, family = "gaussian", alpha = 0, lambda = 2.57 / 10)) Результаты немного отличаются и становятся намного более похожими, когда Я использую намного более высокие лямбды для glmnet.

a11msp

Интригующим. Коэффициенты, похоже, отличаются примерно в 10 раз.

Томка

Ответы:

Разница, которую вы наблюдаете, связана с дополнительным делением на количество наблюдений N, которое GLMNET использует в своей целевой функции, и неявной стандартизацией Y по стандартному отклонению выборки, как показано ниже.

\frac{1}{2 N} {| | \frac{Y}{s_{Y}} - Икс β | |}_{2}^{2} + λ | | β {| |}_{2}^{2} / 2

$\frac{1}{2N}\left\|\frac{y}{s_y}-X\beta\right\|^2_{2}+\lambda\|\beta\|^2_{2}/2$

где мы используем вместо для , $1/n$ $1/(n-1)$ $s_y$

s_{Y} знак равно \frac{\underset{я}{Σ} (Y_{я} - \bar{Y})^{2}}{N}

$s_y=\frac{\sum_i(y_i-\bar{y})^2}{n}$

Дифференцируя по отношению к бете, устанавливая уравнение на ноль,

{Икс}^{T} Икс β - \frac{{Икс}^{T} Y}{s_{Y}} + N λ β знак равно 0

$X^TX\beta-\frac{X^Ty}{s_y}+N\lambda\beta =0$

И, решив для беты, мы получаем оценку,

{\tilde{β}}_{г L M N Е T} знак равно ({Икс}^{T} Икс + N λ я_{п})^{- 1} \frac{{Икс}^{T} Y}{s_{Y}}

$\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}\frac{X^Ty}{s_y}$

Чтобы восстановить оценки (и соответствующие штрафы) по исходной метрике Y, GLMNET умножает как оценки, так и лямбды на и возвращает эти результаты пользователю, $s_y$

{\hat{β}}_{г L M N Е T} знак равно s_{Y} {\tilde{β}}_{г L M N Е T} знак равно ({Икс}^{T} Икс + N λ я_{п})^{- 1} {Икс}^{T} Y

$\hat{\beta}_{GLMNET}=s_y\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}X^Ty$

λ_{U N s T d,} знак равно s_{Y} λ

$\lambda_{unstd.}=s_y\lambda$

Сравните это решение со стандартным выводом регрессии гребня.

\hat{β} знак равно ({Икс}^{T} Икс + λ я_{п})^{- 1} {Икс}^{T} Y

$\hat{\beta}= (X^TX+\lambda I_p)^{-1}X^Ty$

Обратите внимание, что масштабируется на дополнительный коэффициент N. Кроме того, когда мы используем функцию или , штраф будет неявно масштабироваться на . То есть, когда мы используем эти функции для получения оценок коэффициентов для некоторого , мы эффективно получаем оценки для $\lambda$ predict()coef() $1/s_y$ $\lambda^*$ . $\lambda=\lambda^*/s_y$

На основании этих наблюдений, наказание используется в GLMNET должна быть расширена на коэффициент . $s_y/N$

set.seed(123)

n    <- 1000
p   <-  100
X   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

beta1 <- solve(t(X)%*%X+10*diag(p),t(X)%*%(Y))[,1]

fit_glmnet <- glmnet(X,Y, alpha=0, standardize = F, intercept = FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

           [,1]        [,2]
[1,]  0.23793862  0.23793862
[2,]  1.81859695  1.81859695
[3,] -0.06000195 -0.06000195
[4,] -0.04958695 -0.04958695
[5,]  0.41870613  0.41870613
[6,]  1.30244151  1.30244151
[7,]  0.06566168  0.06566168
[8,]  0.44634038  0.44634038
[9,]  0.86477108  0.86477108
[10,] -2.47535340 -2.47535340

Результаты обобщаются на включение перехвата и стандартизированных переменных X. Мы модифицируем стандартизированную X-матрицу, включив в нее столбец единиц и диагональную матрицу, чтобы иметь дополнительный нулевой элемент в позиции [1,1] (т.е. не штрафовать за перехват). Затем вы можете отстранить оценки от оценки их соответствующих стандартных отклонений (снова убедитесь, что вы используете 1 / n при вычислении стандартного отклонения).

{\hat{β}}_{J} знак равно \frac{\tilde{β_{J}}}{s_{{Икс}_{J}}}

$\hat\beta_{j}=\frac{\tilde{\beta_j}}{s_{x_j}}$

{\hat{β}}_{0} знак равно \tilde{β_{0}} - {\bar{Икс}}^{T} \hat{β}

$\hat\beta_{0}=\tilde{\beta_0}-\bar{x}^T\hat{\beta}$

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)
X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}
X_scaled_ones <- cbind(rep(1,n), X_scaled)

beta3 <- solve(t(X_scaled_ones)%*%X_scaled_ones+1000*diag(x = c(0, rep(1,p))),t(X_scaled_ones)%*%(Y))[,1]
beta3 <- c(beta3[1] - crossprod(mean_x,beta3[-1]/sd_x), beta3[-1]/sd_x)

fit_glmnet2 <- glmnet(X,Y, alpha=0, thresh = 1e-20)
beta4 <- as.vector(coef(fit_glmnet2, s = sd_y*1000/n, exact = TRUE))

cbind(beta3[1:10], beta4[1:10])
             [,1]        [,2]
 [1,]  0.24534485  0.24534485
 [2,]  0.17661130  0.17661130
 [3,]  0.86993230  0.86993230
 [4,] -0.12449217 -0.12449217
 [5,] -0.06410361 -0.06410361
 [6,]  0.17568987  0.17568987
 [7,]  0.59773230  0.59773230
 [8,]  0.06594704  0.06594704
 [9,]  0.22860655  0.22860655
[10,]  0.33254206  0.33254206

Добавлен код для отображения стандартизированного X без перехвата:

set.seed(123)

n <- 1000
p <-  100
X <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)

X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}

beta1 <- solve(t(X_scaled)%*%X_scaled+10*diag(p),t(X_scaled)%*%(Y))[,1]

fit_glmnet <- glmnet(X_scaled,Y, alpha=0, standardize = F, intercept = 
FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

             [,1]        [,2]
 [1,]  0.23560948  0.23560948
 [2,]  1.83469846  1.83469846
 [3,] -0.05827086 -0.05827086
 [4,] -0.04927314 -0.04927314
 [5,]  0.41871870  0.41871870
 [6,]  1.28969361  1.28969361
 [7,]  0.06552927  0.06552927
 [8,]  0.44576008  0.44576008
 [9,]  0.90156795  0.90156795
[10,] -2.43163420 -2.43163420

skijunkie
источник

+6. Добро пожаловать в CV и спасибо за ответ на этот старый вопрос в такой ясной форме.

говорит амеба, восстанови Монику

β

$\beta$

\tilde{β}

$\tilde{\beta}$

Я также заметил, что для второй части, где вы сказали «Результаты обобщаются на включение перехвата и стандартизированных переменных X»; для этой части, если вы исключите перехват, то после тех же самых вычислений результаты glmnet станут отличными от вычислений вручную.

user1769197

β

$\beta$

gaussianglmnet()

\begin{matrix} (1) & \frac{1}{2 N} Σ_{я знак равно 1}^{N} (Y_{я} - β_{0} - {Икс}_{я}^{T} β)^{2} + λ Σ_{J знак равно 1}^{п} (α | β_{J} | + (1 - α) β_{J}^{2} / 2), \end{matrix}

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda\sum_{j=1}^p(\alpha|\beta_j| +(1-\alpha)\beta_j^2/2). \tag{1}$

glmnet(x, y, alpha=1) $x$ $\lambda$

\frac{1}{2 N} Σ_{я знак равно 1}^{N} (Y_{я} - β_{0} - {Икс}_{я}^{T} β)^{2} + λ Σ_{J знак равно 1}^{п} | β_{J} |,

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \sum_{j=1}^p |\beta_j|.$ glmnet_2.0-13glmnet(x, y, alpha=0)

λ

$\lambda$

\frac{1}{2 N} Σ_{я знак равно 1}^{N} (Y_{я} - β_{0} - {Икс}_{я}^{T} β)^{2} + λ \frac{1}{2 s_{Y}} Σ_{J знак равно 1}^{п} β_{J}^{2},

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \frac{1}{2s_y} \sum_{j=1}^p \beta_j^2.$

s_{y}

$s_y$

y

$y$

λ / s_{y}

$\lambda/s_y$ .

$y$ $y_0$

\begin{matrix} (2) & \frac{1}{2 N} Σ_{я знак равно 1}^{N} (Y_{0 я} - {Икс}_{я}^{T} γ)^{2} + η Σ_{J знак равно 1}^{п} (α | γ_{J} | + (1 - α) γ_{J}^{2} / 2), \end{matrix}

$\frac{1}{2n} \sum_{i=1}^n (y_{0i}-x_i^T\gamma)^2 +\eta \sum_{j=1}^p(\alpha|\gamma_j| +(1-\alpha)\gamma_j^2/2), \tag{2}$ что эффективно сводит к минимуму

\frac{1}{2 N s_{Y}^{2}} Σ_{я знак равно 1}^{N} (Y_{я} - β_{0} - {Икс}_{я}^{T} β)^{2} + η \frac{α}{s_{Y}} Σ_{J знак равно 1}^{п} | β_{J} | + η \frac{1 - α}{2 s_{Y}^{2}} Σ_{J знак равно 1}^{п} β_{J}^{2},

$\frac{1}{2n s_y^2} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta \frac{\alpha}{s_y} \sum_{j=1}^p |\beta_j| +\eta \frac{1-\alpha}{2s_y^2} \sum_{j=1}^p \beta_j^2,$ или эквивалентно, чтобы минимизировать

\frac{1}{2 N} Σ_{я знак равно 1}^{N} (Y_{я} - β_{0} - {Икс}_{я}^{T} β)^{2} + η s_{Y} α Σ_{J знак равно 1}^{п} | β_{J} | + η (1 - α) Σ_{J знак равно 1}^{п} β_{J}^{2} / 2.

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta s_y \alpha \sum_{j=1}^p |\beta_j| +\eta (1-\alpha) \sum_{j=1}^p \beta_j^2/2.$

Для лассо ( $\alpha=1$ ), масштабирование $\eta$ вернуться, чтобы сообщить о казни как $\eta s_y$ имеет смысл. Тогда для всех $\alpha$ , $\eta s_y$ должен быть указан как штраф, чтобы сохранить непрерывность результатов по всему $\alpha$ , Это, вероятно, является причиной проблемы выше. Отчасти это связано с использованием (2) для решения (1). Только тогда, когда $\alpha=0$ или $\alpha=1$ существует некоторая эквивалентность между задачами (1) и (2) (т. е. соответствие между $\lambda$ в (1) и $\eta$ в (2)). Для любого другого $\alpha\in(0,1)$ задачи (1) и (2) - это две разные проблемы оптимизации, и нет однозначного соответствия между $\lambda$ в (1) и $\eta$ в (2).

Чун Ли
источник

Я не вижу, чем ваш ответ отличается от предыдущего. Не могли бы вы объяснить, пожалуйста?

Firebug

@Firebug Я хотел бы пролить свет на то, почему функция сообщает о лямбде таким образом, который кажется неестественным, если смотреть исключительно с точки зрения регрессии гребня, но имеет смысл (или должен быть таким), если смотреть с точки зрения всего спектра. в том числе и хребет и лассо.

Чун Ли