[Genabel-commits] r722 - in pkg/GenABEL: R inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 18 09:35:39 CEST 2011
Author: yurii
Date: 2011-05-18 09:35:38 +0200 (Wed, 18 May 2011)
New Revision: 722
Modified:
pkg/GenABEL/R/polygenic.R
pkg/GenABEL/R/polylik.R
pkg/GenABEL/inst/doc/index.html
pkg/GenABEL/man/polygenic.Rd
Log:
Modified: pkg/GenABEL/R/polygenic.R
===================================================================
--- pkg/GenABEL/R/polygenic.R 2011-05-18 06:12:10 UTC (rev 721)
+++ pkg/GenABEL/R/polygenic.R 2011-05-18 07:35:38 UTC (rev 722)
@@ -1,4 +1,4 @@
-#' Estimation of polygenic model
+#' Estimation of polygenic model
#'
#' This function maximises the likelihood of the data under polygenic
#' model with covariates an reports twice negative maximum likelihood estimates
@@ -81,8 +81,6 @@
#' @param maxdiffgls max difference allowed in fgls checks
#' @param patchBasedOnFGLS if FGLS checks not passed, 'patch' fixed
#' effect estimates based on FGLS expectation
-#' @param llfun function to compute likelihood (default 'polylik_eigen', also
-#' avalable -- but not recommeneded -- 'polylik')
#' @param ... Optional arguments to be passed to \code{\link{nlm}} or (\code{\link{optim}})
#' minimisation function
#'
@@ -121,7 +119,7 @@
#' Amin N, van Duijn CM, Aulchenko YS. A genomic background based method for
#' association analysis in related individuals. PLoS ONE. 2007 Dec 5;2(12):e1274.
#'
-#' @author Yurii Aulchenko, Gulnara Svischeva
+#' @author Yurii Aulchenko
#'
#' @note
#' Presence of twins may complicate your analysis. Check kinship matrix for
@@ -166,8 +164,7 @@
function(formula,kinship.matrix,data,fixh2,starth2=0.3,trait.type="gaussian",
opt.method="nlm",scaleh2=1,quiet=FALSE,
steptol=1e-8, gradtol = 1e-8, optimbou = 8,
- fglschecks=TRUE,maxnfgls=8,maxdiffgls=1e-4, patchBasedOnFGLS = TRUE,
- llfun = "polylik_eigen", ...) {
+ fglschecks=TRUE,maxnfgls=8,maxdiffgls=1e-4, patchBasedOnFGLS = TRUE, ...) {
if (!missing(data)) if (is(data,"gwaa.data"))
{
checkphengen(data)
@@ -187,9 +184,6 @@
out <- paste("opt.method argument should be one of",optargs,"\n")
stop(out)
}
- if (llfun == "polylik_eigen") llFUN <- polylik_eigen
- else if (llfun == "polylik") llFUN <- polylik
- else stop("llfun should be 'polylik' or 'polylik_eigen'")
if (!missing(data)) attach(data,pos=2,warn.conflicts=FALSE)
if (is(formula,"formula")) {
@@ -243,14 +237,7 @@
tmp <- t(relmat)
relmat[upper.tri(relmat)] <- tmp[upper.tri(tmp)]
rm(tmp);gc()
- if (llfun=="polylik") eigres <- eigen(ginv(relmat),symm=TRUE)
- else if (llfun=="polylik_eigen") {
- eigres <- eigen(relmat,symm=TRUE)
- if (any(eigres$values<0)) {
- #eigres$values <- abs(eigres$values)
- warning("some eigenvalues <=0, taking ABS for det; try option llfun='polylik'",immediate.=TRUE)
- }
- } else stop("cannot be here...")
+ eigres <- eigen(ginv(relmat),symm=TRUE)
if (!quiet) {
cat("LM estimates of fixed parameters:\n")
print(iniest);
@@ -262,21 +249,21 @@
convFGLS <- NULL;
if (!missing(fixh2)) {
- startlik <- llFUN(c(iniest,tvar),y=y,desmat=desmat,relmat=relmat,
- eigenRes=eigres,fixh2=(fixh2/scaleh2),trait.type=trait.type,
+ startlik <- polylik(c(iniest,tvar),y=y,desmat=desmat,relmat=relmat,
+ ervec=eigres$vec,fixh2=(fixh2/scaleh2),trait.type=trait.type,
scaleh2=scaleh2)
if (opt.method=="nlm") {
prnlev <- 0; if (!quiet) prnlev <- 2;
- h2an <- nlm(llFUN,p=c(iniest,tvar),y=y,desmat=desmat,relmat=relmat,eigenRes=eigres,
+ h2an <- nlm(polylik,p=c(iniest,tvar),y=y,desmat=desmat,relmat=relmat,ervec=eigres$vec,
fixh2=(fixh2/scaleh2),trait.type=trait.type,startlik=startlik,scaleh2=scaleh2,
print.level=prnlev,steptol=steptol,gradtol=gradtol,...)
-# h2an <- nlm(llFUN,p=c(iniest,tvar),y=y,desmat=desmat,relmat=relmat,eigenRes=eigres,fixh2=(fixh2/scaleh2),trait.type=trait.type,startlik=startlik,scaleh2=scaleh2,...)
+# h2an <- nlm(polylik,p=c(iniest,tvar),y=y,desmat=desmat,relmat=relmat,ervec=eigres$vec,fixh2=(fixh2/scaleh2),trait.type=trait.type,startlik=startlik,scaleh2=scaleh2,...)
} else {
lower <- c(iniest-inierr*optimbou,1.e-4)
upper <- c(iniest+inierr*optimbou,1)
cntrl <- list(); if (!quiet) cntrl <- list(trace=6,REPORT=1)
- h2an <- optim(fn=llFUN,par=c(iniest,tvar),method="L-BFGS-B",lower=lower,upper=upper,
- y=y,desmat=desmat,relmat=relmat,eigenRes=eigres,fixh2=(fixh2),trait.type=trait.type,
+ h2an <- optim(fn=polylik,par=c(iniest,tvar),method="L-BFGS-B",lower=lower,upper=upper,
+ y=y,desmat=desmat,relmat=relmat,ervec=eigres$vec,fixh2=(fixh2),trait.type=trait.type,
control=cntrl,scaleh2=1,...)
}
} else {
@@ -293,15 +280,15 @@
if (opt.method=="nlm") parsave <- c(iniest,starth2/scaleh2,tvar)
else parsave <- c(iniest,starth2,tvar)
}
- startlik<-llFUN(parsave,y=y,desmat=desmat,relmat=relmat,eigenRes=eigres,
+ startlik<-polylik(parsave,y=y,desmat=desmat,relmat=relmat,ervec=eigres$vec,
trait.type=trait.type,scaleh2=scaleh2)
if (opt.method=="nlm") {
#print(parsave)
prnlev <- 0; if (!quiet) prnlev <- 2;
- h2an <- nlm(llFUN,p=parsave,y=y,desmat=desmat,relmat=relmat,eigenRes=eigres,
+ h2an <- nlm(polylik,p=parsave,y=y,desmat=desmat,relmat=relmat,ervec=eigres$vec,
trait.type=trait.type,startlik=startlik,scaleh2=scaleh2,
print.level=prnlev,steptol=steptol,gradtol=gradtol,...)
-# h2an <- nlm(llFUN,p=c(iniest,(starth2/scaleh2),tvar),y=y,desmat=desmat,relmat=relmat,eigenRes=eigres,trait.type=trait.type,startlik=startlik,scaleh2=scaleh2,...)
+# h2an <- nlm(polylik,p=c(iniest,(starth2/scaleh2),tvar),y=y,desmat=desmat,relmat=relmat,ervec=eigres$vec,trait.type=trait.type,startlik=startlik,scaleh2=scaleh2,...)
parsave <- h2an$estimate
} else {
#print(parsave)
@@ -311,9 +298,7 @@
#print(lower)
#print(upper)
cntrl <- list(); if (!quiet) cntrl <- list(trace=6,REPORT=1)
- h2an <- optim(fn=llFUN,par=parsave,method="L-BFGS-B",lower=lower,upper=upper,
- y=y,desmat=desmat,relmat=relmat,eigenRes=eigres,trait.type=trait.type,
- control=cntrl,scaleh2=1,...)
+ h2an <- optim(fn=polylik,par=parsave,method="L-BFGS-B",lower=lower,upper=upper,y=y,desmat=desmat,relmat=relmat,ervec=eigres$vec,trait.type=trait.type,control=cntrl,scaleh2=1,...)
parsave <- h2an$par
}
#print(parsave)
Modified: pkg/GenABEL/R/polylik.R
===================================================================
--- pkg/GenABEL/R/polylik.R 2011-05-18 06:12:10 UTC (rev 721)
+++ pkg/GenABEL/R/polylik.R 2011-05-18 07:35:38 UTC (rev 722)
@@ -1,4 +1,4 @@
-"polylik" <- function(coeffs,y,desmat,relmat,eigenRes,fixh2,trait.type,
+"polylik" <- function(coeffs,y,desmat,relmat,ervec,fixh2,trait.type,
startlik=.Machine$double.xmax-2,scaleh2=1) {
#print(coeffs)
if (!missing(fixh2)) {
@@ -19,34 +19,12 @@
} else {
qt <- y - desmat %*% fixeff
}
- es <- 1./diag(t(eigenRes$vectors) %*% (sigma) %*% eigenRes$vectors)
- ginvsig <- eigenRes$vectors %*% diag(es,ncol=length(qt)) %*% t(eigenRes$vectors)
+ es <- 1./diag(t(ervec) %*% (sigma) %*% ervec)
+ ginvsig <- ervec %*% diag(es,ncol=length(qt)) %*% t(ervec)
a <- determinant(sigma,logarithm=T)
a <- a$modulus * a$sign
b <- (t(qt) %*% ginvsig %*% qt)
loglik <- a+b
- loglik[1,1]
-}
-"polylik_eigen" <- function(coeffs,y,desmat,relmat,eigenRes,fixh2,trait.type,
- startlik=.Machine$double.xmax-2,scaleh2=1) {
- if (!missing(fixh2)) {
- h2 <- fixh2
- fixeff <- coeffs[1:(length(coeffs)-1)]
- } else {
- h2 <- coeffs[(length(coeffs)-1)] # coeffs[(length(coeffs))]
- fixeff <- coeffs[1:(length(coeffs)-2)]
- }
- h2 <- h2*scaleh2
- tvar <- coeffs[(length(coeffs))]
- if (h2<1e-8 || h2>(1-1e-8) || tvar<1e-6) return(startlik+1)
- if (trait.type=="binomial") {
- expb <- exp(desmat %*% fixeff)
- qt <- y - expb/(1.+expb)
- } else {
- qt <- y - desmat %*% fixeff
- }
- y_0 <- as.vector(t(eigenRes$vectors) %*% qt)
- eigenvalueCOV <- tvar*(eigenRes$values*h2+1.-h2)
- loglik <- sum(log(abs(eigenvalueCOV))) + sum(y_0^2/eigenvalueCOV)
loglik
}
+
Modified: pkg/GenABEL/inst/doc/index.html
===================================================================
--- pkg/GenABEL/inst/doc/index.html 2011-05-18 06:12:10 UTC (rev 721)
+++ pkg/GenABEL/inst/doc/index.html 2011-05-18 07:35:38 UTC (rev 722)
@@ -1,3 +1,3 @@
-<a href="ABEL-tutorial.pdf">GenABEL tutorial</a>: Using GenABEL for genome wide association analysis
+<a href="GenABEL-tutorial.pdf">GenABEL tutorial</a>: Using GenABEL for genome wide association analysis
<br>
<a href="GenABEL-manual.pdf">GenABEL manual</a>: GenABEL reference manual
Modified: pkg/GenABEL/man/polygenic.Rd
===================================================================
--- pkg/GenABEL/man/polygenic.Rd 2011-05-18 06:12:10 UTC (rev 721)
+++ pkg/GenABEL/man/polygenic.Rd 2011-05-18 07:35:38 UTC (rev 722)
@@ -1,11 +1,10 @@
-\name{polygenic}
+\name{polygenic}
\alias{polygenic}
\title{Estimation of polygenic model...}
\usage{polygenic(formula, kinship.matrix, data, fixh2, starth2=0.3,
trait.type="gaussian", opt.method="nlm", scaleh2=1, quiet=FALSE,
steptol=1e-08, gradtol=1e-08, optimbou=8, fglschecks=TRUE,
- maxnfgls=8, maxdiffgls=1e-04, patchBasedOnFGLS=TRUE,
- llfun="polylik_eigen", ...)}
+ maxnfgls=8, maxdiffgls=1e-04, patchBasedOnFGLS=TRUE, ...)}
\description{Estimation of polygenic model}
\details{This function maximises the likelihood of the data under polygenic
model with covariates an reports twice negative maximum likelihood estimates
@@ -82,7 +81,7 @@
Amin N, van Duijn CM, Aulchenko YS. A genomic background based method for
association analysis in related individuals. PLoS ONE. 2007 Dec 5;2(12):e1274.}
-\author{Yurii Aulchenko, Gulnara Svischeva}
+\author{Yurii Aulchenko}
\note{Presence of twins may complicate your analysis. Check kinship matrix for
singularities, or rather use \code{\link{check.marker}} for identification
of twin samples. Take special care in interpretation.
@@ -131,8 +130,6 @@
\item{maxdiffgls}{max difference allowed in fgls checks}
\item{patchBasedOnFGLS}{if FGLS checks not passed, 'patch' fixed
effect estimates based on FGLS expectation}
-\item{llfun}{function to compute likelihood (default 'polylik_eigen', also
-avalable -- but not recommeneded -- 'polylik')}
\item{...}{Optional arguments to be passed to \code{\link{nlm}} or (\code{\link{optim}})
minimisation function}}
\examples{# note that procedure runs on CLEAN data
More information about the Genabel-commits
mailing list