[Genabel-commits] r751 - in pkg/GenABEL: . R inst/unitTests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 15 18:26:16 CEST 2011
Author: yurii
Date: 2011-07-15 18:26:16 +0200 (Fri, 15 Jul 2011)
New Revision: 751
Modified:
pkg/GenABEL/CHANGES.LOG
pkg/GenABEL/R/polygenic.R
pkg/GenABEL/inst/unitTests/report.html
pkg/GenABEL/inst/unitTests/report.txt
pkg/GenABEL/inst/unitTests/reportSummary.txt
Log:
Speeding up 'polygenic'; fixed bug in 'polygenic' causing the offset in names(...$residualY) and names(...$pgresidualY)
Modified: pkg/GenABEL/CHANGES.LOG
===================================================================
--- pkg/GenABEL/CHANGES.LOG 2011-07-15 14:23:17 UTC (rev 750)
+++ pkg/GenABEL/CHANGES.LOG 2011-07-15 16:26:16 UTC (rev 751)
@@ -1,5 +1,12 @@
*** v. 1.6-8 (2011.07.15)
+Speeding up 'polygenic' by avoiding multiple inverses of
+the matrix.
+
+Fixed bug in 'polygenic' causing the offset in
+names(...$residualY) and names(...$pgresidualY) in
+case of missing observations.
+
Updated the merge.snp.data procedure to allow meaningful
merging of mono- and poly- codings. Take care -- this is
new and not tested much yet! You best check few results
Modified: pkg/GenABEL/R/polygenic.R
===================================================================
--- pkg/GenABEL/R/polygenic.R 2011-07-15 14:23:17 UTC (rev 750)
+++ pkg/GenABEL/R/polygenic.R 2011-07-15 16:26:16 UTC (rev 751)
@@ -204,7 +204,8 @@
sglm <- summary(rglm)
iniest <- sglm$coef[,1]
inierr <- sglm$coef[,2]
- phids <- rownames(data)[rownames(data) %in% rownames(mf)]
+ allids <- rownames(data)
+ phids <- allids[rownames(data) %in% rownames(mf)]
#print("bb")
#print(phids)
relmat <- kinship.matrix[phids,phids]*2.0
@@ -220,7 +221,8 @@
y <- formula
if (length(y) != dim(kinship.matrix)[1]) stop("dimension of outcome and kinship.matrix do not match")
mids <- (!is.na(y))
- phids <- data$id[mids]
+ allids <- data$id
+ phids <- allids[mids]
y <- y[mids]
relmat <- kinship.matrix[mids,mids]*2.0
sdy <- sd(y)
@@ -322,7 +324,13 @@
if (fglschecks && missing(fixh2)) {
npar <- length(parsave)
h2 <- parsave[npar-1]*scaleh2
- iSigma <- ginv(h2*relmat + (1-h2)*diag(x=1,ncol=length(y),nrow=length(y)))
+#
+# iSigma <- ginv(h2*relmat + (1-h2)*diag(x=1,ncol=length(y),nrow=length(y)))
+# start new
+ if (llfun=="polylik") eigres <- eigen(relmat,symm=TRUE) # ensure eigres contain eigen of RelMat (not Inv(RelMat))
+ es <- (eigres$value*h2+1.-h2)*parsave[npar]*sdy*sdy
+ iSigma <- (eigres$vec) %*% diag(1./es,ncol=length(es)) %*% t(eigres$vec)
+# END new
LHest <- parsave[1:(npar-2)]
betaFGLS <- as.vector(ginv(t(desmat) %*% iSigma %*% desmat) %*%
(t(desmat) %*% iSigma %*% y))
@@ -414,30 +422,65 @@
}
out$esth2 <- h2
tvar <- h2an$estimate[npar]
- ervec <- eigres$vec
- sigma <- h2*tvar*relmat + (1-h2)*tvar*diag(x=1,ncol=length(y),nrow=length(y))
+# old variant
+# time0 <- proc.time()
+# ervec <- eigres$vec
+# sigma <- h2*tvar*relmat + (1-h2)*tvar*diag(x=1,ncol=length(y),nrow=length(y))
# es <- 1./diag(t(ervec) %*% (sigma) %*% ervec)
# ginvsig <- ervec %*% diag(es,ncol=length(y)) %*% t(ervec)
- out$InvSigma <- ginv(sigma) #ginvsig
+# out$InvSigma <- ginv(sigma) #ginvsig
+# print(proc.time()-time0)
+# end old variant
+# old variant'
+# time0 <- proc.time()
+# ervec <- eigres$vec
+# sigma <- h2*tvar*relmat + (1-h2)*tvar*diag(x=1,ncol=length(y),nrow=length(y))
+# es <- 1./diag(t(ervec) %*% (sigma) %*% ervec)
+# ginvsig <- ervec %*% diag(es,ncol=length(y)) %*% t(ervec)
+# out$InvSigma <- ginv(sigma) #ginvsig
+# print(proc.time()-time0)
+# end old variant'
+# new implementation of InvSigma
+ if (fglschecks && missing(fixh2)) {
+ ginvsig <- iSigma # already computed it in FGLS checks
+ } else {
+ if (llfun=="polylik") eigres <- eigen(relmat,symm=TRUE) # ensure eigres contain eigen of RelMat (not Inv(RelMat))
+ es <- tvar*(eigres$value*h2+1.-h2)
+ #print(es[1:5])
+ #print(eigres$vec[1:5,1:5])
+ #print(((diag(1./es,ncol=length(es))))[1:5,1:5])
+ ginvsig <- (eigres$vec) %*% diag(1./es,ncol=length(es)) %*% t(eigres$vec)
+ #print(ginvsig[1:5,1:5])
+ }
+ out$InvSigma <- ginvsig #ginvsig
+# END new implementation of InvSigma
rownames(out$InvSigma) <- phids
colnames(out$InvSigma) <- phids
- pgres <- as.vector((1.-h2) * tvar * (out$InvSigma %*% out$residualY))
+ InvSigma_x_residualY <- (out$InvSigma %*% out$residualY)
+ pgres <- as.vector((1.-h2) * tvar * InvSigma_x_residualY)
out$measuredIDs <- mids
# need to fix -- now only measured names, while length is >
# names(out$measuredIDs) <- phids
out$pgresidualY <- rep(NA,length(mids))
out$pgresidualY[mids] <- pgres
- names(out$pgresidualY) <- phids
+ names(out$pgresidualY) <- allids #phids
resY <- out$residualY
out$residualY <- rep(NA,length(mids))
out$residualY[mids] <- resY
- names(out$residualY) <- phids
+ names(out$residualY) <- allids #phids
out$call <- match.call()
out$convFGLS <- convFGLS
- a <- determinant(sigma,logarithm=T)
- a <- a$modulus * a$sign
- b <- (t(resY) %*% out$InvSigma %*% resY)
+# old variant
+# time0 <- proc.time()
+# a <- determinant(sigma,logarithm=T)
+# a <- a$modulus * a$sign
+# print(proc.time()-time0)
+# end old variant
+# new variant
+ a <- sum(log(abs(es))) # logarithm of determinant of sigma as sum of logarithm of eigenvalues
+# end new variant
+ b <- (t(resY) %*% InvSigma_x_residualY)
loglik <- a+b
out$h2an$minimum <- as.vector(loglik)
Modified: pkg/GenABEL/inst/unitTests/report.html
===================================================================
--- pkg/GenABEL/inst/unitTests/report.html 2011-07-15 14:23:17 UTC (rev 750)
+++ pkg/GenABEL/inst/unitTests/report.html 2011-07-15 16:26:16 UTC (rev 751)
@@ -1,8 +1,8 @@
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
"http://www.w3.org/TR/html4/transitional.dtd">
-<html><head><title>RUNIT TEST PROTOCOL--Fri Jul 15 16:20:40 2011</title>
+<html><head><title>RUNIT TEST PROTOCOL--Fri Jul 15 18:24:27 2011</title>
</head>
-<body><h1 TRUE>RUNIT TEST PROTOCOL--Fri Jul 15 16:20:40 2011</h1>
+<body><h1 TRUE>RUNIT TEST PROTOCOL--Fri Jul 15 18:24:27 2011</h1>
<p>Number of test functions: 10</p>
<p>Number of errors: 0</p>
<p>Number of failures: 0</p>
@@ -23,7 +23,7 @@
<hr>
<h3 TRUE>Details</h3>
<p><a name="GenABEL unit testing"><h5 TRUE>Test Suite: GenABEL unit testing</h5>
-</a>Test function regexp: ^test.+<br/>Test file regexp: ^runit.+\.[rR]$<br/>Involved directory:<br/>/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests<br/><ul><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.convert.snp.R">Test file: runit.convert.snp.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.convert.snp.R_test.convert.snp">test.convert.snp: (4 checks) ... OK (0.98 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.convert.snp.ped.R">Test file: runit.convert.snp.ped.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.convert.snp.ped.R_test.convert.snp.ped">test.convert.snp.ped: (0 checks) ... OK (0.02 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.descriptives.trait.R">Test file: runit.descriptives.trait.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.descriptives.trait.R_test.descriptives.trait">test.descriptives.trait: (1 checks) ... OK (0.38 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.findRelatives.R">Test file: runit.findRelatives.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.findRelatives.R_test.findRelatives">test.findRelatives: (10 checks) ... OK (77.47 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.impute2xxx.R">Test file: runit.impute2xxx.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.impute2xxx.R_test.impute2databel">test.impute2databel: (23 checks) ... OK (0.47 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.impute2xxx_large.R">Test file: runit.impute2xxx_large.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.impute2xxx_large.R_test.impute2xxx_large">test.impute2xxx_large: (0 checks) ... OK (0 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.iterator.R">Test file: runit.iterator.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.iterator.R_test.qtscore">test.qtscore: (0 checks) ... OK (0 seconds)<br/></a></li><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.iterator.R_test.summary_snp_data">test.summary_snp_data: (3 checks) ... OK (4.41 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.mach2databel.R">Test file: runit.mach2databel.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.mach2databel.R_test.mach2databel">test.mach2databel: (8 checks) ... OK (0.12 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.polylik.R">Test file: runit.polylik.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.polylik.R_test.polylik">test.polylik: (6 checks) ... OK (7.73 seconds)<br/></a></li></ul></li></ul><hr>
+</a>Test function regexp: ^test.+<br/>Test file regexp: ^runit.+\.[rR]$<br/>Involved directory:<br/>/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests<br/><ul><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.convert.snp.R">Test file: runit.convert.snp.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.convert.snp.R_test.convert.snp">test.convert.snp: (4 checks) ... OK (0.93 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.convert.snp.ped.R">Test file: runit.convert.snp.ped.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.convert.snp.ped.R_test.convert.snp.ped">test.convert.snp.ped: (0 checks) ... OK (0.05 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.descriptives.trait.R">Test file: runit.descriptives.trait.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.descriptives.trait.R_test.descriptives.trait">test.descriptives.trait: (1 checks) ... OK (0.42 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.findRelatives.R">Test file: runit.findRelatives.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.findRelatives.R_test.findRelatives">test.findRelatives: (10 checks) ... OK (51.72 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.impute2xxx.R">Test file: runit.impute2xxx.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.impute2xxx.R_test.impute2databel">test.impute2databel: (23 checks) ... OK (0.44 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.impute2xxx_large.R">Test file: runit.impute2xxx_large.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.impute2xxx_large.R_test.impute2xxx_large">test.impute2xxx_large: (0 checks) ... OK (0 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.iterator.R">Test file: runit.iterator.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.iterator.R_test.qtscore">test.qtscore: (0 checks) ... OK (0 seconds)<br/></a></li><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.iterator.R_test.summary_snp_data">test.summary_snp_data: (3 checks) ... OK (4.03 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.mach2databel.R">Test file: runit.mach2databel.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.mach2databel.R_test.mach2databel">test.mach2databel: (8 checks) ... OK (0.13 seconds)<br/></a></li></ul></li><li><a href="/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.polylik.R">Test file: runit.polylik.R</a><ul><li><a name="GenABEL unit testing__Users_yuryaulchenko_eclipse_workspace_genabel_pkg_GenABEL_tests_.._inst_unitTests_runit.polylik.R_test.polylik">test.polylik: (6 checks) ... OK (7.72 seconds)<br/></a></li></ul></li></ul><hr>
<table border="0" width="80%" >
<tr><th>Name</th>
<th>Value</th>
Modified: pkg/GenABEL/inst/unitTests/report.txt
===================================================================
--- pkg/GenABEL/inst/unitTests/report.txt 2011-07-15 14:23:17 UTC (rev 750)
+++ pkg/GenABEL/inst/unitTests/report.txt 2011-07-15 16:26:16 UTC (rev 751)
@@ -1,4 +1,4 @@
-RUNIT TEST PROTOCOL -- Fri Jul 15 16:20:40 2011
+RUNIT TEST PROTOCOL -- Fri Jul 15 18:24:27 2011
***********************************************
Number of test functions: 10
Number of errors: 0
@@ -19,29 +19,29 @@
/Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests
---------------------------
Test file: /Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.convert.snp.R
-test.convert.snp: (4 checks) ... OK (0.98 seconds)
+test.convert.snp: (4 checks) ... OK (0.93 seconds)
---------------------------
Test file: /Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.convert.snp.ped.R
-test.convert.snp.ped: (0 checks) ... OK (0.02 seconds)
+test.convert.snp.ped: (0 checks) ... OK (0.05 seconds)
---------------------------
Test file: /Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.descriptives.trait.R
-test.descriptives.trait: (1 checks) ... OK (0.38 seconds)
+test.descriptives.trait: (1 checks) ... OK (0.42 seconds)
---------------------------
Test file: /Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.findRelatives.R
-test.findRelatives: (10 checks) ... OK (77.47 seconds)
+test.findRelatives: (10 checks) ... OK (51.72 seconds)
---------------------------
Test file: /Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.impute2xxx.R
-test.impute2databel: (23 checks) ... OK (0.47 seconds)
+test.impute2databel: (23 checks) ... OK (0.44 seconds)
---------------------------
Test file: /Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.impute2xxx_large.R
test.impute2xxx_large: (0 checks) ... OK (0 seconds)
---------------------------
Test file: /Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.iterator.R
test.qtscore: (0 checks) ... OK (0 seconds)
-test.summary_snp_data: (3 checks) ... OK (4.41 seconds)
+test.summary_snp_data: (3 checks) ... OK (4.03 seconds)
---------------------------
Test file: /Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.mach2databel.R
-test.mach2databel: (8 checks) ... OK (0.12 seconds)
+test.mach2databel: (8 checks) ... OK (0.13 seconds)
---------------------------
Test file: /Users/yuryaulchenko/eclipse_workspace/genabel/pkg/GenABEL/tests/../inst/unitTests/runit.polylik.R
-test.polylik: (6 checks) ... OK (7.73 seconds)
+test.polylik: (6 checks) ... OK (7.72 seconds)
Modified: pkg/GenABEL/inst/unitTests/reportSummary.txt
===================================================================
--- pkg/GenABEL/inst/unitTests/reportSummary.txt 2011-07-15 14:23:17 UTC (rev 750)
+++ pkg/GenABEL/inst/unitTests/reportSummary.txt 2011-07-15 16:26:16 UTC (rev 751)
@@ -1,4 +1,4 @@
-RUNIT TEST PROTOCOL -- Fri Jul 15 16:20:40 2011
+RUNIT TEST PROTOCOL -- Fri Jul 15 18:24:27 2011
***********************************************
Number of test functions: 10
Number of errors: 0
More information about the Genabel-commits
mailing list