Ridge regression 2

Simulated data

y <- c(10, 5, 8)
y

[1] 10  5  8

set.seed(40)
n <- 3  # individuals
m <- 2  # markers 
W <- matrix(rbinom(n = n * m, size = 2, prob = 0.5), nrow = n, ncol = m)
W

     [,1] [,2]
[1,]    1    0
[2,]    2    0
[3,]    1    1

Ordinary least squares

# SNP effects
a_ols <- solve(t(W) %*% W) %*% t(W) %*% y
a_ols

     [,1]
[1,]    4
[2,]    4

Ridge regression with lambda = 0.1

lambda = 0.1

# Directly obtain ridge solutions
a_ridge1 <- solve(t(W) %*% W + diag(lambda, m)) %*% t(W) %*% y
a_ridge1

         [,1]
[1,] 3.992995
[2,] 3.642732

# Indirectly obtain ridge solutions from OLS
a_ridge2 <- solve(diag(m) + solve(t(W) %*% W) * lambda) %*% a_ols
a_ridge2

         [,1]
[1,] 3.992995
[2,] 3.642732

Ridge regression with lambda = 10

lambda = 10

# Directly obtain ridge solutions
a_ridge3 <- solve(t(W) %*% W + diag(lambda, m)) %*% t(W) %*% y
a_ridge3

          [,1]
[1,] 1.7142857
[2,] 0.5714286

# Indirectly obtain ridge solutions from OLS
a_ridge4 <- solve(diag(m) + solve(t(W) %*% W) * lambda) %*% a_ols
a_ridge4

          [,1]
[1,] 1.7142857
[2,] 0.5714286

glmnet R package

# install.packages('glmnet')
library(glmnet)

library(BGLR)  # use the wheat dataset
data(wheat)
`?`(wheat)
y <- wheat.Y[, 2]  # use the second phenotype
W <- wheat.X
dim(W)

[1]  599 1279

# lambda = 0.2
a_ridge5 <- glmnet(W, y, alpha = 0, lambda = 0.2)  # use alpha = 0 for ridge
names(a_ridge5)

 [1] "a0"        "beta"      "df"        "dim"       "lambda"   
 [6] "dev.ratio" "nulldev"   "npasses"   "jerr"      "offset"   
[11] "call"      "nobs"     

a_ridge5$a0

       s0 
-1.014927 

head(a_ridge5$beta)  # the first six marker effects

6 x 1 sparse Matrix of class "dgCMatrix"
                  s0
wPt.0538 -0.01779275
wPt.8463 -0.06768257
wPt.6348 -0.02959608
wPt.9992  0.01144462
wPt.2838 -0.03363664
wPt.8266 -0.04012267

# lambda = 10
a_ridge6 <- glmnet(W, y, alpha = 0, lambda = 10)
a_ridge6$a0

        s0 
-0.4921524 

head(a_ridge6$beta)  # the first six marker effects

6 x 1 sparse Matrix of class "dgCMatrix"
                    s0
wPt.0538  0.0010426896
wPt.8463  0.0094178218
wPt.6348  0.0025946282
wPt.9992  0.0004541006
wPt.2838  0.0015773713
wPt.8266 -0.0037416161

# lambda = 1000
a_ridge7 <- glmnet(W, y, alpha = 0, lambda = 1000)
a_ridge7$a0

         s0 
-0.03698497 

head(a_ridge7$beta)  # the first six marker effects

6 x 1 sparse Matrix of class "dgCMatrix"
                   s0
wPt.0538 2.180747e-04
wPt.8463 3.633714e-04
wPt.6348 1.920957e-04
wPt.9992 2.709440e-04
wPt.2838 2.832555e-04
wPt.8266 1.151529e-05