[Genabel-commits] r1979 - pkg/MixABEL/inst/unitTests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun May 17 23:35:49 CEST 2015
Author: lckarssen
Date: 2015-05-17 23:35:49 +0200 (Sun, 17 May 2015)
New Revision: 1979
Modified:
pkg/MixABEL/inst/unitTests/runit.GWFGLS.R
Log:
Summary: Only whitespace cleanup in one of MixABEL's test files.
(No functional changes.)
Modified: pkg/MixABEL/inst/unitTests/runit.GWFGLS.R
===================================================================
--- pkg/MixABEL/inst/unitTests/runit.GWFGLS.R 2015-05-15 17:17:37 UTC (rev 1978)
+++ pkg/MixABEL/inst/unitTests/runit.GWFGLS.R 2015-05-17 21:35:49 UTC (rev 1979)
@@ -22,44 +22,44 @@
test.GWFGLS_all_data_types <- function()
{
-
+
unlink("tmp*")
-
-
+
+
data(ge03d2.clean)
-
+
idprop <- runif(1,min=0.2,max=0.5)
IDS <- sort(sample(1:nids(ge03d2.clean),floor(nids(ge03d2.clean)*idprop)))
dfd <- ge03d2.clean[IDS,autosomal(ge03d2.clean)]
gkin <- ibs(dfd,w="freq")
-
+
snpprop <- runif(1,min=0.005,max=0.01)
SNPS <- sort(sample(1:nsnps(dfd),floor(nsnps(dfd)*snpprop)))
dfd <- dfd[,SNPS]
s <- summary(dfd at gtdata)
maf <- pmin(s$Q.2,1-s$Q.2)
dfd <- dfd[,(maf>0.05)]
-
+
modelh2 <- 0.8
covars <- c("sex","age")
s2 <- 12^2
betas_cov <- c(170,12,0.01)
-
+
X <- as.matrix(cbind(rep(1,length(IDS)),phdata(dfd)[,covars]))
grel <- gkin
grel[upper.tri(grel)] <- t(grel)[upper.tri(grel)]
grel <- 2*grel
- Y <- as.vector(rmvnorm(1,mean=(X %*% betas_cov),sigma=s2*(modelh2*grel+diag(dim(grel)[1])*(1-modelh2))))
+ Y <- as.vector(rmvnorm(1,mean=(X %*% betas_cov),sigma=s2*(modelh2*grel+diag(dim(grel)[1])*(1-modelh2))))
#+ 10*as.numeric(dfd at gtdata[,10])
length(Y)
Y[2] <- NA
dfd at phdata$Y <- Y
-
+
mdl <- Y~age+sex
h2 <- polygenic(mdl,data=dfd,kin=gkin,gradtol=1e-5,steptol=1e-5,maxdiffgls=1e-2)
h2$h2an
mdl <- Y~age+sex+SNP
-
+
# run tests with different types of genodata
gtOld <- dfd at gtdata
gtReal <- as.double(gtOld)
@@ -72,97 +72,97 @@
for (model in c("additive","dominantB","recessiveB","overdominant","genotypic")) {
print(c(varcov,verbosity,residuals,test,model))
aIR <- GWFGLS(mdl,data=phdata(dfd),genodata = gtReal,
- model.SNP=model, test = test, residuals = residuals,
- verbosity = verbosity, varcov = varcov,
+ model.SNP=model, test = test, residuals = residuals,
+ verbosity = verbosity, varcov = varcov,
include.means = include.means)
- aIO <- GWFGLS(mdl,data=phdata(dfd),genodata = gtOld,
- model.SNP=model, test = test, residuals = residuals,
- verbosity = verbosity, varcov = varcov,
+ aIO <- GWFGLS(mdl,data=phdata(dfd),genodata = gtOld,
+ model.SNP=model, test = test, residuals = residuals,
+ verbosity = verbosity, varcov = varcov,
include.means = include.means)
- aIO1<- GWFGLS(mdl,data=dfd,
- model.SNP=model, test = test, residuals = residuals,
- verbosity = verbosity, varcov = varcov,
+ aIO1<- GWFGLS(mdl,data=dfd,
+ model.SNP=model, test = test, residuals = residuals,
+ verbosity = verbosity, varcov = varcov,
include.means = include.means)
- aIN <- GWFGLS(mdl,data=phdata(dfd),genodata = gtNew,
- model.SNP=model, test = test, residuals = residuals,
- verbosity = verbosity, varcov = varcov,
+ aIN <- GWFGLS(mdl,data=phdata(dfd),genodata = gtNew,
+ model.SNP=model, test = test, residuals = residuals,
+ verbosity = verbosity, varcov = varcov,
include.means = include.means)
checkEquals(aIR,aIO,tolerance=5*.Machine$double.eps^0.5)
checkEquals(aIR,aIO1,tolerance=5*.Machine$double.eps^0.5)
checkEquals(aIR,aIN,tolerance=5*.Machine$double.eps^0.5)
aWR <- GWFGLS(mdl,data=phdata(dfd),W=h2$InvSigma,genodata = gtReal,
- model.SNP=model, test = test, residuals = residuals,
- verbosity = verbosity, varcov = varcov,
+ model.SNP=model, test = test, residuals = residuals,
+ verbosity = verbosity, varcov = varcov,
include.means = include.means)
aWO <- GWFGLS(mdl,data=phdata(dfd),W=h2$InvSigma,genodata = gtOld,
- model.SNP=model, test = test, residuals = residuals,
- verbosity = verbosity, varcov = varcov,
+ model.SNP=model, test = test, residuals = residuals,
+ verbosity = verbosity, varcov = varcov,
include.means = include.means)
aWO1<- GWFGLS(mdl,data=dfd,W=h2$InvSigma,
- model.SNP=model, test = test, residuals = residuals,
- verbosity = verbosity, varcov = varcov,
+ model.SNP=model, test = test, residuals = residuals,
+ verbosity = verbosity, varcov = varcov,
include.means = include.means)
aWN <- GWFGLS(mdl,data=phdata(dfd),W=h2$InvSigma,genodata = gtNew,
- model.SNP=model, test = test, residuals = residuals,
- verbosity = verbosity, varcov = varcov,
+ model.SNP=model, test = test, residuals = residuals,
+ verbosity = verbosity, varcov = varcov,
include.means = include.means)
checkEquals(aWR,aWO,tolerance=5*.Machine$double.eps^0.5)
checkEquals(aWR,aWO1,tolerance=5*.Machine$double.eps^0.5)
checkEquals(aWR,aWN,tolerance=5*.Machine$double.eps^0.5)
}
-
-
+
+
rm(list=ls());gc()
-
+
unlink("tmp*")
}
test.GWFGLS_equal_to_LM <- function()
{
-
+
unlink("tmp*")
-
-
+
+
data(ge03d2.clean)
-
+
idprop <- runif(1,min=0.25,max=0.5)
IDS <- sort(sample(1:nids(ge03d2.clean),floor(nids(ge03d2.clean)*idprop)))
dfd <- ge03d2.clean[IDS,autosomal(ge03d2.clean)]
gkin <- ibs(dfd,w="freq")
-
+
snpprop <- runif(1,min=0.005,max=0.01)
SNPS <- sort(sample(1:nsnps(dfd),floor(nsnps(dfd)*snpprop)))
dfd <- dfd[,SNPS]
smr <- summary(dfd at gtdata)
maf <- pmin(smr$Q.2,1-smr$Q.2)
dfd <- dfd[,(maf>0.01 & smr$P.11>2 & smr$P.12>2 & smr$P.22>2)]
-
+
modelh2 <- 0.8
covars <- c("sex","age")
s2 <- 12^2
betas_cov <- c(170,12,0.01)
-
+
X <- as.matrix(cbind(rep(1,length(IDS)),phdata(dfd)[,covars]))
grel <- gkin
grel[upper.tri(grel)] <- t(grel)[upper.tri(grel)]
grel <- 2*grel
- Y <- as.vector(rmvnorm(1,mean=(X %*% betas_cov),sigma=s2*(modelh2*grel+diag(dim(grel)[1])*(1-modelh2))))
+ Y <- as.vector(rmvnorm(1,mean=(X %*% betas_cov),sigma=s2*(modelh2*grel+diag(dim(grel)[1])*(1-modelh2))))
#+ 10*as.numeric(dfd at gtdata[,10])
length(Y)
Y[2] <- NA
dfd at phdata$Y <- Y
-
+
mdl <- Y~age+sex
h2 <- polygenic(mdl,data=dfd,kin=gkin,gradtol=1e-5,steptol=1e-5,maxdiffgls=1e-2)
h2$h2an
mdl <- Y~age+sex+SNP
-
+
# check equivalence to LM
for (model in c("additive","dominantB","recessiveB","overdominant","genotypic")) {
print(model)
aIR <- GWFGLS(mdl,data=dfd,
- model.SNP=model, test = "wald", residuals = FALSE,
+ model.SNP=model, test = "wald", residuals = FALSE,
verbosity = 2, varcov = TRUE, with.lm = TRUE)
passed <- rep(TRUE,dim(dfd at gtdata)[2])
for (i in 1:dim(dfd at gtdata)[2]) {
@@ -174,22 +174,22 @@
aIR$LM <- aIR$LM[passed]
checkEqualsNumeric(aIR$LM,aIR$FGLS[, "P-value"],tolerance=5*.Machine$double.eps^0.5)
}
-
+
rm(list=ls());gc()
-
+
unlink("tmp*")
}
test.GWFGLS_Old_vs_New <- function()
{
-
+
# library("RUnit")
# library("GenABEL")
# library("DatABEL")
# library("MixABEL")
# library("MASS")
# library("mvtnorm")
-
+
data(ge03d2.clean)
propNids <- runif(1,min=0.25,max=0.5)
IDS <- sort(sample(1:nids(ge03d2.clean),round(propNids*nids(ge03d2.clean))))
@@ -199,12 +199,12 @@
maf <- pmin(s$Q.2,1-s$Q.2)
df <- df[,(maf>0.05)]
gkin <- ibs(df[,autosomal(df)],w="freq")
-
+
h2mod <- 0.5
covars <- c("sex","age")
s2 <- 12^2
betas_cov <- c(170,12,0.01)
-
+
X <- as.matrix(cbind(rep(1,NIDS),phdata(ge03d2.clean)[1:NIDS,covars]))
grel <- gkin
grel[upper.tri(grel)] <- t(grel)[upper.tri(grel)]
@@ -212,14 +212,14 @@
Y <- as.vector(rmvnorm(1,mean=(X %*% betas_cov),sigma=s2*(h2mod*grel+diag(dim(grel)[1])*(1-h2mod))))
Y[2] <- NA
df at phdata$Y <- Y
-
+
h2 <- polygenic(Y~sex+age,data=df,kin=gkin,gradtol=1e-5,steptol=1e-5,maxdiffgls=1e-2)
h2$h2an
-
+
propSnps <- runif(1,min=0.005,max=0.01)
SNPS <- sort(sample(1:nsnps(df),round(propSnps*nsnps(df))))
NSNPS <- length(SNPS)
-
+
dfn <- df[,SNPS]
smr <- summary(dfn at gtdata)
## actually should automatically do that!
@@ -231,14 +231,14 @@
#dfn <- dfn[,which(maf<0.01 | smr$P.11 < 3 | smr$P.12 < 3 & smr$P.22 < 3)]
nsnps(dfn)
#summary(dfn at gtdata)
-
+
# check that iterator produces NAs with mono-snps
resIteratorNew <- GWFGLS(Y~age+sex*SNP,data=dfn at phdata,W=h2$InvSigma,old=FALSE,ver=2,varcov=TRUE,
include.means=FALSE,genodata=matrix(rep(1,length(Y)*2),ncol=2))
print(resIteratorNew)
checkTrue(all(is.na(resIteratorNew[1,4:dim(resIteratorNew)[2]])))
checkTrue(all(is.na(resIteratorNew[2,4:dim(resIteratorNew)[2]])))
-
+
t0 <- proc.time()
resIteratorOld <- GWFGLS(Y~age+sex*SNP,data=dfn,W=h2$InvSigma,old=TRUE)
proc.time() - t0
@@ -250,8 +250,8 @@
resIteratorNew <- resIteratorNew[!newNA,]
resIteratorOld <- resIteratorOld[!newNA,]
checkEquals(resIteratorOld,resIteratorNew)
-
-
+
+
for (varcov in c(TRUE,FALSE))
for (include.means in c(TRUE,FALSE))
for (verbosity in c(1,2))
@@ -278,5 +278,5 @@
checkEquals(resIteratorOld,resIteratorNew,tolerance=5*(.Machine$double.eps^0.5))
#}
}
-
+
}
More information about the Genabel-commits
mailing list