Added LT-MAPIT Functions

lorinanthony · Jul 24, 2018 · 3ba58df · 3ba58df
1 parent 11dc069
commit 3ba58df
Show file tree

Hide file tree

Showing 4 changed files with 10 additions and 20 deletions.
diff --git a/OpenMP Version/LT_Simulations_OpenMP.R b/OpenMP Version/LT_Simulations_OpenMP.R
@@ -38,7 +38,7 @@ source("MAPIT_OpenMP.R"); sourceCpp("MAPIT_OpenMP.cpp")
 # (5) k = assumed disease prevelance in the population
 # (6) samp = # of samples to analyze
 
-ind = 1e6; nsnp = 1e4; pve=0.6; rho=0.5; k = 0.001; samp = 500
+ind = 1e6; nsnp = 1e4; pve=0.6; rho=0.5; k = 0.01; samp = 500
 
 ### Simulate the genotypes such that all variants have minor allele frequency (MAF) > 0.05 ###
 # NOTE: As in the paper, we center and scale each genotypic vector such that every SNP has mean 0 and standard deviation 1.
@@ -96,13 +96,8 @@ z=z_marginal+z_epi+z_err
 thresh=qnorm(1-k,mean=0,sd=1)
 
 ### Find the Number of Cases and Controls ###
-n.cases = sum(z>thresh); n.cases/length(z)
-n.controls = sum(z<=thresh); n.controls/length(z)
-
-#Bernoulli Distributed Case and Control Data
-Pheno=rep(NA,ind); 
-Pheno[z<=thresh] = rtruncnorm(n.controls,b=thresh) 
-Pheno[z>thresh] = rtruncnorm(n.cases,a=thresh)
+Pheno=rep(0,ind)
+Pheno[z>thresh] = 1
 
 ### Subsample a particular number of cases and controls ###
 cases = sample(which(z>thresh),samp/2,replace = FALSE)
@@ -132,7 +127,7 @@ cores = detectCores()
 
 ### Run LT-MAPIT ###
 ptm <- proc.time() #Start clock
-mapit = MAPIT(t(X),y,hybrid=FALSE,test="davies",cores=cores)
+mapit = MAPIT(t(X),y,hybrid=FALSE,test="davies",k=k,cores=cores)
 proc.time() - ptm #Stop clock
 
 davies.pvals = mapit$pvalues

diff --git a/OpenMP Version/MAPIT_OpenMP.R b/OpenMP Version/MAPIT_OpenMP.R
@@ -44,10 +44,10 @@ MAPIT = function(X, y, W = NULL,C = NULL,hybrid = TRUE,threshold = 0.05,test = N
 
   ### Check to See if we should run LT-MAPIT ###
   if(sum(y%in%c(0,1))==length(y)){
-    if(k = NULL){warning("The liability threshold model is going to be used but no disease prevelance is defined! Using the default 0.5.");k = 0.5}
+    if(is.null(k)){warning("The liability threshold model is going to be used but no disease prevelance is defined! Using the default 0.5.");k = 0.5}
     ### Save the labels ###
     cc = y
-
+    
     ### Find the Number of Cases and Controls ###
     n.cases = sum(cc==1)
     n.controls = sum(cc==0)

diff --git a/Standard Version/LT_Simulations.R b/Standard Version/LT_Simulations.R
@@ -96,13 +96,8 @@ z=z_marginal+z_epi+z_err
 thresh=qnorm(1-k,mean=0,sd=1)
 
 ### Find the Number of Cases and Controls ###
-n.cases = sum(z>thresh); n.cases/length(z)
-n.controls = sum(z<=thresh); n.controls/length(z)
-
-#Bernoulli Distributed Case and Control Data
-Pheno=rep(NA,ind); 
-Pheno[z<=thresh] = rtruncnorm(n.controls,b=thresh) 
-Pheno[z>thresh] = rtruncnorm(n.cases,a=thresh)
+Pheno=rep(0,ind)
+Pheno[z>thresh] = 1
 
 ### Subsample a particular number of cases and controls ###
 cases = sample(which(z>thresh),samp/2,replace = FALSE)
@@ -132,7 +127,7 @@ cores = detectCores()
 
 ### Run LT-MAPIT ###
 ptm <- proc.time() #Start clock
-mapit = MAPIT(t(X),y,hybrid=FALSE,test="davies",cores=cores)
+mapit = MAPIT(t(X),y,hybrid=FALSE,test="davies",k=k,cores=cores)
 proc.time() - ptm #Stop clock
 
 davies.pvals = mapit$pvalues

diff --git a/Standard Version/MAPIT.R b/Standard Version/MAPIT.R
@@ -39,7 +39,7 @@ MAPIT = function(X, y, W = NULL,C = NULL,hybrid = TRUE,threshold = 0.05,test = N
 
   ### Check to See if we should run LT-MAPIT ###
   if(sum(y%in%c(0,1))==length(y)){
-    if(k = NULL){warning("The liability threshold model is going to be used but no disease prevelance is defined! Using the default 0.5.");k = 0.5}
+    if(is.null(k)){warning("The liability threshold model is going to be used but no disease prevelance is defined! Using the default 0.5.");k = 0.5}
     ### Save the labels ###
     cc = y