[Genabel-commits] r752 - in pkg/GenABEL: . R inst/unitTests man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 15 19:50:55 CEST 2011
Author: yurii
Date: 2011-07-15 19:50:55 +0200 (Fri, 15 Jul 2011)
New Revision: 752
Added:
pkg/GenABEL/inst/unitTests/runit.polygenic.R
Modified:
pkg/GenABEL/CHANGES.LOG
pkg/GenABEL/DESCRIPTION
pkg/GenABEL/R/alleleID.R
pkg/GenABEL/R/polygenic.R
pkg/GenABEL/man/snp.data-class.Rd
Log:
Added patch of bug [#1322] + regression test. Some other small fixes coming up in checks.
Modified: pkg/GenABEL/CHANGES.LOG
===================================================================
--- pkg/GenABEL/CHANGES.LOG 2011-07-15 16:26:16 UTC (rev 751)
+++ pkg/GenABEL/CHANGES.LOG 2011-07-15 17:50:55 UTC (rev 752)
@@ -1,5 +1,12 @@
*** v. 1.6-8 (2011.07.15)
+Added patch of bug [#1322] + regression test (contributed
+by Nicola Pirastu, see
+https://lists.r-forge.r-project.org/pipermail/genabel-devel/2011-May/000276.html).
+This is still suboptimal treatment -- the covariate is just dropped;
+in say 'lm' it is 'kept' with NA as estimate, and IDs with missing
+data dropped, which is not the case with 'polygenic'.
+
Speeding up 'polygenic' by avoiding multiple inverses of
the matrix.
Modified: pkg/GenABEL/DESCRIPTION
===================================================================
--- pkg/GenABEL/DESCRIPTION 2011-07-15 16:26:16 UTC (rev 751)
+++ pkg/GenABEL/DESCRIPTION 2011-07-15 17:50:55 UTC (rev 752)
@@ -2,7 +2,7 @@
Type: Package
Title: genome-wide SNP association analysis
Version: 1.6-8
-Date: 2011-05-18
+Date: 2011-07-15
Author: GenABEL developers
Maintainer: GenABEL developers <genabel-devel at r-forge.r-project.org>
Depends: R (>= 2.10.0), methods, MASS
Modified: pkg/GenABEL/R/alleleID.R
===================================================================
--- pkg/GenABEL/R/alleleID.R 2011-07-15 16:26:16 UTC (rev 751)
+++ pkg/GenABEL/R/alleleID.R 2011-07-15 17:50:55 UTC (rev 752)
@@ -160,5 +160,5 @@
else
return(mono)
}
- return(newC)
+ return(mono)
}
Modified: pkg/GenABEL/R/polygenic.R
===================================================================
--- pkg/GenABEL/R/polygenic.R 2011-07-15 16:26:16 UTC (rev 751)
+++ pkg/GenABEL/R/polygenic.R 2011-07-15 17:50:55 UTC (rev 752)
@@ -192,6 +192,36 @@
else stop("llfun should be 'polylik' or 'polylik_eigen'")
if (!missing(data)) attach(data,pos=2,warn.conflicts=FALSE)
+
+# beging patch bug #1322 (by Nicola Pirastu)
+ if (is(formula, "formula")){
+ mf <- model.frame(formula, data, na.action = na.omit,
+ drop.unused.levels = TRUE)
+ ciccio=function(x)nlevels(as.factor(x))
+ livelli=apply(mf,2,ciccio)
+ ch=length(livelli)
+ bad=which(livelli<2)
+ if(length(bad)>0) warning(paste("Variable",names(mf)[bad],"has only one value: Dropped"))
+ livelli=livelli[livelli>1]
+ mf=mf[names(livelli)]
+ if(length(livelli)>1){
+ if(length(livelli)==2){ formula=as.formula(paste(names(mf)[1],"~",names(mf)[2]))
+
+ }else{
+ formula=(paste(names(mf)[1],"~",names(mf)[2]))
+ for(i in 3:length(livelli)){
+
+ formula=paste(formula,"+",names(mf)[i])
+
+ }
+ formula=as.formula(formula)
+ }
+ }else{
+ stop("All covariates have only one value")
+ }
+ }
+# end patch bug #1322
+
if (is(formula,"formula")) {
clafo <- "formula"
mf <- model.frame(formula,data,na.action=na.omit,drop.unused.levels=TRUE)
Added: pkg/GenABEL/inst/unitTests/runit.polygenic.R
===================================================================
--- pkg/GenABEL/inst/unitTests/runit.polygenic.R (rev 0)
+++ pkg/GenABEL/inst/unitTests/runit.polygenic.R 2011-07-15 17:50:55 UTC (rev 752)
@@ -0,0 +1,49 @@
+### --- Test setup ---
+#
+# regression tests
+#
+
+if(FALSE) {
+ ## Not really needed, but can be handy when writing tests
+ library(RUnit)
+ library(GenABEL)
+}
+
+### do not run
+#stop("SKIP THIS TEST")
+###
+
+### ---- common functions and data -----
+
+#source(paste("../inst/unitTests/shared_functions.R"))
+#source(paste(path,"/shared_functions.R",sep=""))
+
+### --- Test functions ---
+
+test.polygenic.Bug1322 <- function()
+{
+ data(ge03d2.clean)
+ df <- ge03d2.clean[1:200,autosomal(ge03d2.clean)]
+ gkin <- ibs(df,w="freq")
+# get a SNP which is informative + no missings
+ s <- summary(gtdata(df))
+ maf <- pmin(s[,"Q.2"],1.-s[,"Q.2"])
+ misSnp <- (rownames(s)[which(maf>0.2 & s[,"CallRate"]<1)])[1]
+ snpX <- as.numeric(gtdata(df[, misSnp]))
+ print(table(snpX,exclude=NULL))
+ phdata(df)$snpX <- snpX
+ pol1 <- polygenic(height~sex+age+snpX,df,kin=gkin,quiet=TRUE)
+ print(pol1$esth2)
+ nUsedIds1 <- sum(pol1$measuredIDs)
+# introduce missing indicator
+ snpInd <- rep(1,length(snpX))
+ snpInd[is.na(snpX)] <- NA
+ checkEquals(nUsedIds1,sum(!is.na(snpInd)))
+# this does did not work!
+ phdata(df)$snpInd <- snpInd
+ pol0 <- polygenic(height~sex+age+snpInd,df,kin=gkin,quiet=TRUE)
+# correct way to do comparison
+ pol0 <- polygenic(height~sex+age,df[!is.na(snpX),],kin=gkin[!is.na(snpX),!is.na(snpX)],quiet=TRUE)
+ nUsedIds0 <- sum(pol0$measuredIDs)
+ checkEquals(nUsedIds1,nUsedIds0)
+}
\ No newline at end of file
Property changes on: pkg/GenABEL/inst/unitTests/runit.polygenic.R
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Modified: pkg/GenABEL/man/snp.data-class.Rd
===================================================================
--- pkg/GenABEL/man/snp.data-class.Rd 2011-07-15 16:26:16 UTC (rev 751)
+++ pkg/GenABEL/man/snp.data-class.Rd 2011-07-15 17:50:55 UTC (rev 752)
@@ -37,6 +37,9 @@
\alias{coding}
\alias{coding,gwaa.data-method}
\alias{coding,snp.data-method}
+\alias{coding<-}
+\alias{coding<-,gwaa.data-method}
+\alias{coding<-,snp.data-method}
\alias{refallele}
\alias{refallele,gwaa.data-method}
\alias{refallele,snp.data-method}
@@ -103,6 +106,8 @@
\code{signature(object = "snp.data")}: extracts strand}
\item{coding}{\code{signature(object = "gwaa.data")},
\code{signature(object = "snp.data")}: extracts coding}
+ \item{coding<-}{\code{signature(object = "gwaa.data")},
+ \code{signature(object = "snp.data")}: assign coding}
\item{refallele}{\code{signature(object = "gwaa.data")},
\code{signature(object = "snp.data")}: extracts reference allele}
\item{effallele}{\code{signature(object = "gwaa.data")},
More information about the Genabel-commits
mailing list