12
$\begingroup$

I've read that ridge regression could be achieved by simply adding rows of data to the original data matrix, where each row is constructed using 0 for the dependent variables and the square root of $k$ or zero for the independent variables. One extra row is then added for each independent variable.

I was wondering whether it is possible to derive a proof for all cases, including for logistic regression or other GLMs.

$\endgroup$
12
$\begingroup$

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

(Often a constant is required, but not shrunken. In that case it is included in the $\beta$ and predictors -- but if you don't want to shrink it, you don't have a corresponding row for the pseudo observation. Or if you do want to shrink it, you do have a row for it. I'll write it as if it's not counted in the $p$, and not shrunken, as it's the more complicated case. The other case is a trivial change from this.)

We can write the second term as $p$ pseudo-observations if we can write each "y" and each of the corresponding $(p+1)$-vectors "x" such that

$(y_{n+j}-x_{n+j}^T\beta)^2=\lambda\beta_j^2\,,\quad j=1,\ldots,p$

But by inspection, simply let $y_{n+j}=0$, let $x_{n+j,j}=\sqrt{\lambda}$ and let all other $x_{n+j,k}=0$ (including $x_{n+j,0}=0$ typically).

Then

$(y_{n+j}-[x_{n+j,0}\beta_0+x_{n+j,1}\beta_1+x_{n+j,2}\beta_2+...+x_{n+j,p}\beta_p])^2=\lambda\beta_j^2$.

This works for linear regression. It doesn't work for logistic regression, because ordinary logistic regression doesn't minimize a sum of squared residuals.

[Ridge regression isn't the only thing that can be done via such pseudo-observation tricks -- they come up in a number of other contexts]

$\endgroup$
  • $\begingroup$ Thanks, I was already struggling with rewriting everything from logistic regression, but I simply couldn't implement the phoney data method. And I don't trust my own abilities sufficiently to be able to say that it is impossible. $\endgroup$ – Snowflake Feb 10 '15 at 12:28
  • $\begingroup$ At least I don't think it is. I'll take another look at the likelihood function. $\endgroup$ – Glen_b -Reinstate Monica Feb 10 '15 at 12:31
  • 3
    $\begingroup$ +1 Additional related regression tricks are introduced in answers at stats.stackexchange.com/a/32753 and stats.stackexchange.com/a/26187, inter alia. $\endgroup$ – whuber Feb 10 '15 at 14:15
  • $\begingroup$ GLMs are fit using iteratively reweighted least squares though, as in bwlewis.github.io/GLM, and so within each iteration one can subsitute the regular weighted least squares step with a ridge penalized weighted least squares step to get a ridge penalized GLM. In fact, in combination with adaptive ridge penalties this is used to fit L0 penalized GLMs, as in the l0ara package, see biodatamining.biomedcentral.com/articles/10.1186/… and journals.plos.org/plosone/article?id=10.1371/… $\endgroup$ – Tom Wenseleers Aug 29 at 19:25
  • $\begingroup$ @TomWenseleers thanks, yes, that makes complete sense $\endgroup$ – Glen_b -Reinstate Monica Aug 30 at 1:05
0
$\begingroup$

Generalizing this recipe to GLMs indeed is not difficult as GLMs are usually fit using iteratively reweighted least squares. Hence, within each iteration one can subsitute the regular weighted least squares step with a ridge penalized weighted least squares step to get a ridge penalized GLM. In fact, in combination with adaptive ridge penalties this recipe is used to fit L0 penalized GLMs (aka best subset, ie GLMs where the total number of nonzero coefficients are penalized). This has been implemented for example in the l0ara package, see this paper and this one for details.

It's also worth noting that the fastest closed-form way of solving a regular ridge regression is using

lmridge_solve = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  solve(crossprod(X) + diag(lambdas), crossprod(X, y))[, 1]
}

for the case where n>=p, or using

lmridge_solve_largep = function (X, Y, lambda) (t(X) %*% solve(tcrossprod(X)+lambda*diag(nrow(X)), Y))[,1]

when p>n and for a model without intercept.

This is faster than using the row augmentation recipe, i.e. doing

lmridge_rbind = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  qr.solve(rbind(X, diag(sqrt(lambdas))), c(y, rep(0, ncol(X))))
}

If you would happen to need nonnegativity constraints on your fitted coefficients then you can just do

library(nnls)

nnlmridge_solve = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  nnls(A=crossprod(X)+diag(lambdas), b=crossprod(X,Y))$x
}

which then gives a slightly more accurate result btw than

nnlmridge_rbind = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  nnls(A=rbind(X,diag(sqrt(lambdas))), b=c(Y,rep(0,ncol(X))))$x 
}

(and strictly speaking only the solution nnls(A=crossprod(X)+diag(lambdas), b=crossprod(X,Y))$x is then the correct one).

I haven't yet figured out how the nonnegativity constrained case could be further optimized for the p > n case - let me know if anyone would happen to know how to do this... [lmridge_nnls_largep = function (X, Y, lambda) t(X) %*% nnls(A=tcrossprod(X)+lambda*diag(nrow(X)), b=Y)$x doesn't work]

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.