Polygenic model with known lambda = sigma2e / sigma2a

polygenic.R

n <- 1000
m <- 200
X <- matrix(rbinom(m*n, 2, 0.5), n)
p <- apply(X, 2, mean)/2
pv <- apply(X, 2, sd)
Xp <- t((t(X) - 2*p) / pv)
A <- Xp %*% t(Xp) / m
A[1:10,1:10]


lambda <- 0.8/0.2

eff <- rnorm(m)
y <- drop(X %*% eff)
y <- (y - mean(y)) / sd(y) * sqrt(0.2) + rnorm(n, sd=sqrt(0.8))
mean(y)
sd(y)
var(y)

eigenA <- eigen(A)
eig.vec <- eigenA$vectors
eig.val <- eigenA$values
M <- diag(1/sqrt(eig.val*lambda+1)) %*% t(eig.vec)
Y <- M %*% y
X.Data <- data.frame(M %*% X)
lm1 <- lm(Y~.-1,data=X.Data)

anova(lm1)