[Genabel-commits] r663 - in pkg: DatABEL GenABEL GenABEL/R GenABEL/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 25 22:03:21 CET 2011
Author: yurii
Date: 2011-02-25 22:03:20 +0100 (Fri, 25 Feb 2011)
New Revision: 663
Added:
pkg/GenABEL/R/polygenic_hglm.R
pkg/GenABEL/man/polygenic_hglm.Rd
Modified:
pkg/DatABEL/DESCRIPTION
pkg/GenABEL/CHANGES.LOG
pkg/GenABEL/DESCRIPTION
pkg/GenABEL/NAMESPACE
pkg/GenABEL/R/GenABEL-package.R
pkg/GenABEL/R/zzz.R
pkg/GenABEL/generate_documentation.R
pkg/GenABEL/man/GenABEL-package.Rd
Log:
initial version of polygenic_hglm
Modified: pkg/DatABEL/DESCRIPTION
===================================================================
--- pkg/DatABEL/DESCRIPTION 2011-02-25 13:14:01 UTC (rev 662)
+++ pkg/DatABEL/DESCRIPTION 2011-02-25 21:03:20 UTC (rev 663)
@@ -3,7 +3,7 @@
Title: file-based access to large matrices stored on HDD in binary format
Version: 0.9-3
Date: 2011-02-09
-Author: Yurii Aulchenko, Stepan Yakovenko
+Author: Yurii Aulchenko, Stepan Yakovenko, Erik Roos, Marcel Kempenaar
Maintainer: Yurii Aulchenko <i.aoultchenko at erasmusmc.nl>
Depends: R (>= 2.4.0), methods
Description: a package providing interface to C++ FILEVECTOR
Modified: pkg/GenABEL/CHANGES.LOG
===================================================================
--- pkg/GenABEL/CHANGES.LOG 2011-02-25 13:14:01 UTC (rev 662)
+++ pkg/GenABEL/CHANGES.LOG 2011-02-25 21:03:20 UTC (rev 663)
@@ -1,3 +1,8 @@
+*** v. 1.6-6 (2011.02.25)
+
+With suggestion and guidance from Xia Shen, added procedure 'polygenic_hglm'.
+At the moment very experimental, and needs to pass checks before can be released.
+
***** v. 1.6-5 (2011.02.07)
Added '#include <cstdarg>' to iterator.cpp to solve
Modified: pkg/GenABEL/DESCRIPTION
===================================================================
--- pkg/GenABEL/DESCRIPTION 2011-02-25 13:14:01 UTC (rev 662)
+++ pkg/GenABEL/DESCRIPTION 2011-02-25 21:03:20 UTC (rev 663)
@@ -1,13 +1,13 @@
Package: GenABEL
Type: Package
Title: genome-wide SNP association analysis
-Version: 1.6-5
-Date: 2011-02-07
-Author: Yurii Aulchenko, Maksim Struchalin, et al.
+Version: 1.6-6
+Date: 2011-02-24
+Author: Yurii Aulchenko et al.
Maintainer: Yurii Aulchenko <i.aoultchenko at erasmusmc.nl>
Depends: R (>= 2.10.0), methods, MASS
-Suggests: qvalue, genetics, haplo.stats, DatABEL (>= 0.9-0)
+Suggests: qvalue, genetics, haplo.stats, DatABEL (>= 0.9-0), hglm
Description: a package for genome-wide association analysis between
- quantitative or binary traits and single-nucleiotide
+ quantitative or binary traits and single-nucleotide
polymorphisms (SNPs).
License: GPL (>= 2)
Modified: pkg/GenABEL/NAMESPACE
===================================================================
--- pkg/GenABEL/NAMESPACE 2011-02-25 13:14:01 UTC (rev 662)
+++ pkg/GenABEL/NAMESPACE 2011-02-25 21:03:20 UTC (rev 663)
@@ -74,6 +74,7 @@
plot.scan.gwaa,
plot.scan.gwaa.2D,
polygenic,
+ polygenic_hglm,
r2fast,
r2fast.old,
redundant,
Modified: pkg/GenABEL/R/GenABEL-package.R
===================================================================
--- pkg/GenABEL/R/GenABEL-package.R 2011-02-25 13:14:01 UTC (rev 662)
+++ pkg/GenABEL/R/GenABEL-package.R 2011-02-25 21:03:20 UTC (rev 663)
@@ -57,7 +57,7 @@
#' For converting of GenABEL's data to other formats, see
#' \code{\link{export.merlin}} (MERLIN and MACH formats),
#' \code{\link{export.impute}} (IMPUTE, SNPTEST and CHIAMO formats),
-#' \code{\link{export.plink}} (PLINK forat, also exports phenotypic data).
+#' \code{\link{export.plink}} (PLINK format, also exports phenotypic data).
#'
#' To load the data, see \code{\link{load.gwaa.data}}.
#'
Added: pkg/GenABEL/R/polygenic_hglm.R
===================================================================
--- pkg/GenABEL/R/polygenic_hglm.R (rev 0)
+++ pkg/GenABEL/R/polygenic_hglm.R 2011-02-25 21:03:20 UTC (rev 663)
@@ -0,0 +1,96 @@
+#' Estimation of polygenic model
+#'
+#' Estimation of polygenic model using hierarchical
+#' GLM (hglm package)
+#'
+#' @param formula Formula describing fixed effects to be used in analysis, e.g.
+#' y ~ a + b means that outcome (y) depends on two covariates, a and b.
+#' If no covariates used in analysis, skip the right-hand side of the
+#' equation.
+#' @param kinship.matrix Kinship matrix, as provided by e.g. ibs(,weight="freq"),
+#' or estimated outside of GenABEL from pedigree data.
+#' @param data An (optional) object of \code{\link{gwaa.data-class}} or a data frame with
+#' outcome and covariates
+#' @param family a description of the error distribution and link function to be
+#' used in the mean part of the model (see \code{\link{family}} for details of
+#' family functions)
+#' @param ... other parameters passed to \code{\link{hglm}} call
+#'
+#' @author Xia Shen, Yurii Aulchenko
+#'
+#' @references
+#' need to put reference here
+#'
+#' @seealso
+#' \code{\link{mmscore}},
+#' \code{\link{grammar}}
+#' \code{\link{polygenic}}
+#'
+#' @examples
+#' data(ge03d2ex.clean)
+#' df <- ge03d2ex.clean[,autosomal(ge03d2ex.clean)]
+#' gkin <- ibs(df,w="freq")
+#'
+#' # for quantitative traits
+#' h2ht <- polygenic_hglm(height ~ sex + age,kin=gkin,df)
+#' # estimate of heritability
+#' h2ht$esth2
+#' # other parameters
+#' h2ht$h2an
+#'
+#' # for binary traits
+#' h2dm <- polygenic_hglm(dm2 ~ sex + age,kin=gkin,df,family=binomial(link = 'logit'))
+#' # estimated parameters
+#' h2dm$h2an
+#'
+#'
+#' @keywords htest
+#'
+"polygenic_hglm" <- function(formula,kinship.matrix,data,family=gaussian(), ... )
+{
+ if (!require(hglm))
+ stop("this function requires 'hglm' package to be installed")
+ if (!missing(data)) if (is(data,"gwaa.data"))
+ {
+# checkphengen(data)
+ data <- phdata(data)
+ }
+ if (!missing(data))
+ if (!is(data,"data.frame"))
+ stop("data should be of gwaa.data or data.frame class")
+ allids <- data$id
+
+ relmat <- kinship.matrix
+ relmat[upper.tri(relmat)] <- t(relmat)[upper.tri(relmat)]
+ mf <- model.frame(formula,data,na.action=na.omit,drop.unused.levels=TRUE)
+ y <- model.response(mf)
+ desmat <- model.matrix(formula,mf)
+ phids <- rownames(data)[rownames(data) %in% rownames(mf)]
+ mids <- (allids %in% phids)
+ relmat <- relmat[mids,mids]
+ relmat <- relmat*2.0
+ s <- svd(relmat)
+ L <- s$u %*% diag(sqrt(s$d))
+ res_hglm <- hglm(y = y, X = desmat, Z = L, family = family, ... )
+ #sum_res_hglm <- summary(res_hglm)
+
+ out <- list()
+
+ out$measuredIDs <- mids
+ out$hglm <- res_hglm
+ out$h2an <- list()
+ tVar <- res_hglm$varRan+res_hglm$varFix
+ out$esth2 <- res_hglm$varRan/tVar
+ out$h2an$estimate <- c(res_hglm$fixef,out$esth2,tVar)
+ names(out$h2an$estimate)[(length(out$h2an$estimate)-1):(length(out$h2an$estimate))] <-
+ c("h2","totalVar")
+ out$pgresidualY <- rep(NA,length(mids))
+ out$pgresidualY[mids] <- y-res_hglm$fv
+ out$residualY <- rep(NA,length(mids))
+ out$residualY[mids] <- y - desmat %*% res_hglm$fixef
+ out$InvSigma <- ginv(tVar*out$esth2*relmat +
+ diag(tVar*(1-out$esth2),ncol=length(y),nrow=length(y)))
+
+ class(out) <- "polygenic"
+ out
+}
\ No newline at end of file
Modified: pkg/GenABEL/R/zzz.R
===================================================================
--- pkg/GenABEL/R/zzz.R 2011-02-25 13:14:01 UTC (rev 662)
+++ pkg/GenABEL/R/zzz.R 2011-02-25 21:03:20 UTC (rev 663)
@@ -1,6 +1,6 @@
.onLoad <- function(lib, pkg) {
- GenABEL.version <- "1.6-5"
- cat("GenABEL v.",GenABEL.version,"(February 07, 2011) loaded\n")
+ GenABEL.version <- "1.6-6"
+ cat("GenABEL v.",GenABEL.version,"(February 24, 2011) loaded\n")
# check for updates and news
address <- c(
Modified: pkg/GenABEL/generate_documentation.R
===================================================================
--- pkg/GenABEL/generate_documentation.R 2011-02-25 13:14:01 UTC (rev 662)
+++ pkg/GenABEL/generate_documentation.R 2011-02-25 21:03:20 UTC (rev 663)
@@ -12,6 +12,7 @@
"mach2databel.R",
#"phdata.R",
"polygenic.R",
+ "polygenic_hglm.R",
"qtscore.R"
#,
#"summary.scan.gwaa.R"
Modified: pkg/GenABEL/man/GenABEL-package.Rd
===================================================================
--- pkg/GenABEL/man/GenABEL-package.Rd 2011-02-25 13:14:01 UTC (rev 662)
+++ pkg/GenABEL/man/GenABEL-package.Rd 2011-02-25 21:03:20 UTC (rev 663)
@@ -58,7 +58,7 @@
For converting of GenABEL's data to other formats, see
\code{\link{export.merlin}} (MERLIN and MACH formats),
\code{\link{export.impute}} (IMPUTE, SNPTEST and CHIAMO formats),
-\code{\link{export.plink}} (PLINK forat, also exports phenotypic data).
+\code{\link{export.plink}} (PLINK format, also exports phenotypic data).
To load the data, see \code{\link{load.gwaa.data}}.
@@ -133,7 +133,8 @@
\code{\link{plot.check.marker}}.}
\alias{GenABEL-package}
\alias{GenABEL}
-\author{Yurii Aulchenko}
+\author{Yurii Aulchenko et al.
+(see help pages for specific functions)}
\references{If you use GenABEL package in your analysis, please cite the following work:
Aulchenko Y.S., Ripke S., Isaacs A., van Duijn C.M. GenABEL: an R package
Added: pkg/GenABEL/man/polygenic_hglm.Rd
===================================================================
--- pkg/GenABEL/man/polygenic_hglm.Rd (rev 0)
+++ pkg/GenABEL/man/polygenic_hglm.Rd 2011-02-25 21:03:20 UTC (rev 663)
@@ -0,0 +1,40 @@
+\name{polygenic_hglm}
+\alias{polygenic_hglm}
+\title{Estimation of polygenic model...}
+\usage{polygenic_hglm(formula, kinship.matrix, data, family=gaussian(), ...)}
+\description{Estimation of polygenic model}
+\details{Estimation of polygenic model using hierarchical
+GLM (hglm package)}
+\author{Xia Shen, Yurii Aulchenko}
+\references{need to put reference here}
+\seealso{\code{\link{mmscore}},
+\code{\link{grammar}}
+\code{\link{polygenic}}}
+\keyword{htest}
+\arguments{\item{formula}{Formula describing fixed effects to be used in analysis, e.g.
+y ~ a + b means that outcome (y) depends on two covariates, a and b.
+If no covariates used in analysis, skip the right-hand side of the
+equation.}
+\item{kinship.matrix}{Kinship matrix, as provided by e.g. ibs(,weight="freq"),
+or estimated outside of GenABEL from pedigree data.}
+\item{data}{An (optional) object of \code{\link{gwaa.data-class}} or a data frame with
+outcome and covariates}
+\item{family}{a description of the error distribution and link function to be
+used in the mean part of the model (see \code{\link{family}} for details of
+family functions)}
+\item{...}{other parameters passed to \code{\link{hglm}} call}}
+\examples{data(ge03d2ex.clean)
+df <- ge03d2ex.clean[,autosomal(ge03d2ex.clean)]
+gkin <- ibs(df,w="freq")
+
+# for quantitative traits
+h2ht <- polygenic_hglm(height ~ sex + age,kin=gkin,df)
+# estimate of heritability
+h2ht$esth2
+# other parameters
+h2ht$h2an
+
+# for binary traits
+h2dm <- polygenic_hglm(dm2 ~ sex + age,kin=gkin,df,family=binomial(link = 'logit'))
+# estimated parameters
+h2dm$h2an}
More information about the Genabel-commits
mailing list