5
$\begingroup$

Is it possible to penalize coefficients toward a number other than zero in a ridge regression in R?

For example, let's say I have dependent variable Y and independent variables X1,X2,X3, and X4. Because of the multicollinear nature of the ivs, ridge regression is appropriate. But say I'm fairly certain that the coefficient of X1 is near 5, X2 is near 1, X3 is near -1, and X4 is near -5.

Is there a ridge package and method in R where I can implement penalties on the coefficients of the ivs that penalize them toward those numbers instead of 0? I'd love to see an example in R with my example data, if possible. Thank you.

$\endgroup$
2
  • $\begingroup$ For a solution that also works with logistic regression or other generalized linear models, see here $\endgroup$ Jan 4, 2019 at 12:07
  • 4
    $\begingroup$ That's known as nonzero centered ridge regresqsion and can be achieved by row-augmenting your covariate matrix X with a diagonal matrix with sqrt(lambdas) along the diagonal and adding center * sqrt(lambdas) to your outcome variable, where center is the desired target values towards which your coefficients should be shrunk. See sciencedirect.com/science/article/abs/pii/S0047259X21001603 and references therein. $\endgroup$ Jun 11, 2023 at 7:51

3 Answers 3

11
$\begingroup$

A simple way to do this is to subtract the "centering value" of the coefficient times its associated variable from the left-hand side. To go with your example,

$Y = \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4 + e$

Assume the coefficient values should be centered at (5,1,-1,-5) respectively. Then:

$Y - 5X_1 -X_2 +X_3 +5X_4 = (\beta_1-5)X_1 + (\beta_2-1)X_2 + (\beta_3+1)X_3 + (\beta_4+5)X_4 + e$

and, redefining terms, you have:

$Y^* = \beta_1^*X_1 + \beta_2^*X_2 + \beta_3^*X_3 + \beta_4^*X_4 + e$

A standard ridge regression would shrink the $\beta_i^*$ towards 0, which is equivalent to shrinking the original $\beta_i$ towards the specified centering values. to see this, consider a fully-shrunk $\beta_4^* = 0$, then $\beta_4+5 = 0$ and therefore $\beta_4 = -5$. Shrinkage accomplished!


A slightly more formal demonstration follows. The required objective function in $\beta$ is $$ (Y-X\beta)^T (Y-X\beta) + (\beta - \gamma)^T(\beta-\gamma) $$ where $\gamma$ is the vector of centring values. In the first term make the substitution $ Y-X\beta=Y-X\gamma-X\beta+X\gamma=Y^*-X\beta^*$ (the errors for a given fit are of course unchanged by this reëxpression of the model), and in the second $\beta-\gamma =\beta^*$; to obtain a zero-centered objective function in $\beta^*$: $$ (Y^* - X\beta^*)^T (Y^* - X\beta^*) + \lambda{\beta^*}^T{\beta^*}^T $$

$\endgroup$
10
  • 1
    $\begingroup$ What if you want to do it with logistic regression? $\endgroup$ Jan 4, 2019 at 11:10
  • 1
    $\begingroup$ But how does this solution allow for a given level of shrinkage towards a given centering value of the ridge regression? I don't see an lambda / regularisation parameter appear in here? $\endgroup$ Jun 20, 2023 at 14:02
  • 3
    $\begingroup$ @TomWenseleers - this reformulates the problem into a standard linear regression that can be solved using any ridge regression, Lasso, Elastic net, or other shrinkage routines with the centering being at the desired values. So, the answer to "is it possible" is clearly meant to be "yes" and "is there a ridge package" is "yes, any ridge package will do" ... as I say in the text, "a standard ridge regression would..." $\endgroup$
    – jbowman
    Jun 20, 2023 at 14:46
  • 3
    $\begingroup$ @TomWenseleers the part of the question that is asking how to write specific R code is nowadays off-topic. A question should have a focus on a statistical problem and coding examples can be part of explanations or the background of a problem but should not be the topic itself. Instructions how to write code can be found in many other places, like here: stats.meta.stackexchange.com/questions/793/… $\endgroup$ Jun 20, 2023 at 19:16
  • 2
    $\begingroup$ I think R code would be otiose for something so simple, especially as @TomWenseleers has already provided an illustration for the example given in the question. I've taken the liberty of adding some linear algebra as a footnote, but feel free of course to remove or edit it. (I know it's rather trivial, but perhaps it wasn't immediately clear to readers new to ridge regression that the two methods - yours & Tom's - give the same result.) $\endgroup$ Jun 21, 2023 at 10:17
5
$\begingroup$

That's known as nonzero centered ridge regression and can be achieved by row-augmenting your covariate matrix X with a diagonal matrix with sqrt(lambdas) along the diagonal and concatenating prior.mean * sqrt(lambdas) to your outcome variable, where prior.mean is the desired target values towards which your coefficients should be shrunk and lambdas is the vector with the desired L2 ridge penalty on your coefficients; and then performing an ordinary least-squares regression of the (row-augmented) target vector on the (row-augmented) matrix X. See https://www.sciencedirect.com/science/article/abs/pii/S0047259X21001603 and references therein. You can also compare with this post for the same recipe for a regular zero-centered ridge regression.

For your example:

First, generate some example data:

set.seed(123)

# Number of observations
n <- 100

# Generate some correlated variables
X1 <- rnorm(n, mean = 5)
X2 <- X1 + rnorm(n, mean = 1)
X3 <- X1 + rnorm(n, mean = -1)
X4 <- X1 + rnorm(n, mean = -5)

# Make a covariate matrix X
X <- cbind(X1, X2, X3, X4)

# Generate a response y
y <- 5*X1 + X2 - X3 - 5*X4 + rnorm(n)

# This is our center (the coefficients we believe are "true", ie the mean of our implicit Gaussian prior)
prior.mean <- c(5, 1, -1, -5)

Next, we can apply a non-zero centered ridge regression. Note that lambda is a regularization parameter and would typically be chosen via cross-validation:

# Regularization parameter 
# note: you can use a vector of lambda values if you would like
# to penalize different coefficients to different extents
# (e.g. penalty on intercept is usually set to zero)
lambda <- 100

# Number of features
p <- ncol(X)

# Augment X
X_augmented <- rbind(X, diag(sqrt(lambda), p))

# Augment y
y_augmented <- c(y, prior.mean * sqrt(lambda))

# Solve nonzero centered ridge regression as this least squares problem
beta <- solve(t(X_augmented) %*% X_augmented, t(X_augmented) %*% y_augmented)

So in other words all we need to do to solve this nonzero centered least-square ridge regression is solve a least square problem but with some extra observations added. We can do this just using a regular linear model (lm), so no need even to use any fancy ridge regression package like ridge or glmnet :

lm(y_augmented ~ X_augmented)

For example if we would set lambda to zero we would then get the standard (zero intercept) least square result

         [,1]
X1  5.1977104
X2  0.8903286
X3 -1.0616111
X4 -5.1220296

but with a lambda of 100 we would get values shrunken more towards the center of your Gaussian prior (nonzero centered ridge regression is in a Bayesian sense equivalent to the MAP estimate with a nonzero-centered Gaussian prior, while standard ridge regression would be equivalent to using a zero-centered Gaussian prior).

         [,1]
X1  5.0477868
X2  0.9773907
X3 -1.0079250
X4 -5.0728437

Adjusting the response to shrink towards your prior mean comes down to the same thing:

y_star <- y - X %*% prior.mean
X_aug <- rbind(X, diag(sqrt(lambda), ncol(X)))
y_aug_star <- c(y_star, rep(0, ncol(X)))
beta_star <- solve(t(X_aug) %*% X_aug, t(X_aug) %*% y_aug_star)
beta <- beta_star + matrix(prior.mean, nrow = length(prior.mean), ncol = 1)

Some of the maths behind this: the objective function for standard (zero-centered) Ridge Regression is:

argmin(β) [(y - Xβ)ᵀ(y - Xβ) + λβᵀβ]

This corresponds to a zero-centered Gaussian prior for β.

To consider a non-zero centered Ridge Regression, we adjust the objective function to center at a given target γ:

argmin(β) [(y - Xβ)ᵀ(y - Xβ) + λ(β - γ)ᵀ(β - γ)]

Now, let's reformulate this as a least squares problem:

If we define z = [y; γ*sqrt(λ)] and W = [X; sqrt(λ)I], where I is the identity matrix, the problem becomes:

argmin(β) [(z - Wβ)ᵀ(z - Wβ)]

To see this, we can expand out the terms:

argmin(β) [(z - Wβ)ᵀ(z - Wβ)]

= argmin(β) [zᵀz - zᵀWβ - βᵀWᵀz + βᵀWᵀWβ]

= argmin(β) [(y - Xβ)ᵀ(y - Xβ) + γᵀγλ - 2γᵀβλ + βᵀβλ]

= argmin(β) [(y - Xβ)ᵀ(y - Xβ) + λ(β - γ)ᵀ(β - γ)]

Therefore, solving the least squares problem with z = [y; γ*sqrt(λ)] and W = [X; sqrt(λ)I] corresponds to solving the non-zero centered Ridge Regression problem. The regularization term now encourages the coefficients to be close to γ instead of zero.

In a Bayesian context, the standard Ridge regression minimizes the negative log posterior, which is equivalent to minimizing:

[(y - Xβ)ᵀ(y - Xβ) + λβᵀβ]

This can be viewed as the sum of the negative log likelihood [(y - Xβ)ᵀ(y - Xβ)] and the negative log prior λβᵀβ, where λβᵀβ corresponds to the log of a Gaussian distribution with mean zero and variance 1/λ.

If you move to a non-zero centered Ridge regression, you're effectively changing the prior to a Gaussian distribution that is centered around a non-zero mean. Specifically, for a center γ, the new objective function becomes:

[(y - Xβ)ᵀ(y - Xβ) + λ(β - γ)ᵀ(β - γ)]

Which is equivalent to minimizing the sum of the negative log likelihood [(y - Xβ)ᵀ(y - Xβ)] and the negative log prior λ(β - γ)ᵀ(β - γ). Here, λ(β - γ)ᵀ(β - γ) corresponds to the log of a Gaussian distribution with mean γ and variance 1/λ.

(Full disclosure: this is a ChatGPT4 based answer, but I verified everything for correctness, i.e. I tested the code & checked the maths)

$\endgroup$
10
  • 2
    $\begingroup$ Please don't generate answers with chatGPT especially when there is already a perfectly sufficient answer. This answers is verbose and it could be made more intuitive by shortening it. I would suggest to fully write out the matrix equation from the example. $\endgroup$ Jun 13, 2023 at 16:16
  • 4
    $\begingroup$ None of the other answers gave working R code so I thought this was useful... And none made reference to nonzero centered ridge regression, which is what this technique is known as. $\endgroup$ Jun 13, 2023 at 16:27
  • 1
    $\begingroup$ @SextusEmpiricus Also the question asked explicitly about "Ridge Regression in R" so working R code I think is quite essential here. The Bayesian approach suggested above would work, but I believe the OP was interested in a faster ML solution. $\endgroup$ Jun 13, 2023 at 16:37
  • 1
    $\begingroup$ @SextusEmpiricus The fact that the OP did not flag any of the other answers as the correct one to me also suggests he wasn't happy with any of the other answers... $\endgroup$ Jun 15, 2023 at 8:33
  • 1
    $\begingroup$ "weighted least squares"? - it's just ordinary least squares. $\endgroup$ Jun 20, 2023 at 12:48
2
$\begingroup$

Ridge regression sets a normal prior centered at zero for each parameter. To penalize coefficients towards different values, just center the priors around your target instead of around 0. I think you can accomplish this using the bayesGLM function of the arm package with the following parameters:

  • prior.mean : set this to vector containing the values you want to use as the centers for your priors (the values you want to penalize "towards")
  • prior.df : set this to Inf to get normal priors. Possibly not necessary. There should be a way to control your prior variance, but it's not clear to me how this would be accomplished from looking at the documentation.

http://cran.r-project.org/web/packages/arm/arm.pdf

$\endgroup$
2
  • $\begingroup$ Hmm, I'm sure I'm missing something but I tried this and it gave me the same coefficients as a normal lm() regression $\endgroup$
    – James
    Nov 19, 2013 at 5:42
  • $\begingroup$ You would also still need to set the SD of the priors, which to mimick ridge regression (interpreted as the MAP solution of a Bayesian LM with Gaussian priors then) with a given lambda would still be a function of the noise and your data & the desired level of regularisation lambda. I think the OP was rather interested in a ML solution, in which case my answer is better & more convient I think (and computationally much faster)... $\endgroup$ Jun 15, 2023 at 8:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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