## Ordinary least squares

set.seed(40)
n <- 3  # individuals
m <- 5  # markers
W <- matrix(rbinom(n = n * m, size = 2, prob = 0.5), nrow = n, ncol = m)
W
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    2    0    1    0    0
[3,]    1    1    1    1    2
y <- c(10, 5, 8)
y
[1] 10  5  8
# SNP effects
try(a <- solve(t(W) %*% W) %*% t(W) %*% y)

# determinant
det(t(W) %*% W)
[1] 0

## The role of $$\lambda$$

lambda <- 0.1
diag(lambda, m)
     [,1] [,2] [,3] [,4] [,5]
[1,]  0.1  0.0  0.0  0.0  0.0
[2,]  0.0  0.1  0.0  0.0  0.0
[3,]  0.0  0.0  0.1  0.0  0.0
[4,]  0.0  0.0  0.0  0.1  0.0
[5,]  0.0  0.0  0.0  0.0  0.1
# W'W
t(W) %*% W
     [,1] [,2] [,3] [,4] [,5]
[1,]    6    1    3    1    2
[2,]    1    1    1    1    2
[3,]    3    1    2    1    2
[4,]    1    1    1    1    2
[5,]    2    2    2    2    4
# SNP effects
a <- solve(t(W) %*% W + diag(lambda, m)) %*% t(W) %*% y
a
          [,1]
[1,]  7.367792
[2,]  1.531720
[3,] -8.711284
[4,]  1.531720
[5,]  3.063440
# determinant
det(t(W) %*% W + diag(lambda, m))
[1] 0.10041

## Scalar form

# marker 1
a[1, ]
[1] 7.367792
# marker 2
a[2, ]
[1] 1.53172
# marker 1
(W[, 1] %*% (y - W[, 2:5] %*% matrix(a[-1, ])))/(sum(W[, 1]^2) + diag(lambda,
1))
         [,1]
[1,] 7.367792
# marker 2
(W[, 2] %*% (y - W[, -2] %*% matrix(a[-2, ])))/(sum(W[, 2]^2) + diag(lambda,
1))
        [,1]
[1,] 1.53172

## Marker specific shrinkage

# marker 1
sum(W[, 1]^2)/(sum(W[, 1]^2) + diag(lambda, 1))
          [,1]
[1,] 0.9836066
# marker 2
sum(W[, 2]^2)/(sum(W[, 2]^2) + diag(lambda, 1))
          [,1]
[1,] 0.9090909
# marker 3
sum(W[, 3]^2)/(sum(W[, 3]^2) + diag(lambda, 1))
         [,1]
[1,] 0.952381

February 1, 2018