n = 100
p = 1000
m = 10
# Generate design matrix (genotype matrix)
set.seed(1)
f = runif( p, .1, .2 ) # simulate minor allele frequency
x = matrix( nrow = n, ncol = p )
colnames(x) = 1:p
for(j in 1:p)
x[,j] = rbinom( n, 2, f[j] )
# Generate true effect sizes
causal_snps = sample( 1:p, m )
beta = rep( 0, p )
set.seed(1)
beta[causal_snps] = rnorm(m, mean = 0, sd = 2 )
# Generate continuous (phenotype) data
y.cont = x %*% beta + rnorm(n, 0, 1)
# Generate binary (phenotype) data
prob = exp(x %*% beta)/(1 + exp(x %*% beta))
y.bin = sapply(1:n, function(i)rbinom(1,1,prob[i]) )
# Fix scale parameter tau
tau = 0.022
# GWASinlps analysis
cor_xy = c(cor(x,y.cont))
names(cor_xy) = colnames(x)
nlps_cont = nlpsLM(y.cont, x, cor_xy=cor_xy, prior="mom",
tau=tau, k0=2, rxx=0.3, niter=10000, verbose=TRUE)
nlps_cont
library(fastglm)
mode(x) = "double" #needed for fastglm() function below
mmle_xy = apply( x, 2, function(z) coef( fastglm(y=y.bin,
x=cbind(1,matrix(z)), family = binomial(link = "logit")) )[2] )
nlps_bin = nlpsGLM(y.bin, x, mmle_xy=mmle_xy, prior="mom",
tau=tau, k0=2, rxx=0.3, niter=10000, verbose=TRUE)
nlps_bin
Run the code above in your browser using DataLab