Genomic best linear unbiased prediction

Data

library(BGLR)
data(mice)
y <- mice.pheno$Obesity.BMI
W0 <- mice.X
dim(W0)

[1]  1814 10346

p <- colSums(W0)/(2 * nrow(W0))
maf <- ifelse(p > 0.5, 1 - p, p)  # or pmin(p, 1-p)
maf.index <- which(maf < 0.1)
p0 <- p[-maf.index]
W0 <- W0[, -maf.index]
dim(W0)

[1] 1814 9340

Genomic relationship matrix

W <- scale(W0, center = TRUE, scale = FALSE)
G <- tcrossprod(W)/sum(2 * p0 * (1 - p0))

Genomic best linear unbiased prediction

sigma2g <- 0.0007119091
sigma2e <- 0.002941943

I <- diag(nrow(G))  # diagonal matrix
Ginv <- solve(G)  # inverse of G
ghat <- solve(I + Ginv * (sigma2e/sigma2g)) %*% matrix(y)  # GBLUP

head(ghat)

                   [,1]
A048005080 -0.004752057
A048006063  0.030142424
A048006555  0.003771144
A048007096  0.018250869
A048010273 -0.019232979
A048010371  0.002255862

tail(ghat)

                   [,1]
A084291611  0.001567918
A084291618  0.034358666
A084291628  0.018702081
A084291785  0.010346188
A084291787 -0.011264509
A084292044  0.024042050