[Genabel-commits] r623 - in pkg/GenABEL: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Dec 7 22:42:44 CET 2010
Author: yurii
Date: 2010-12-07 22:42:43 +0100 (Tue, 07 Dec 2010)
New Revision: 623
Modified:
pkg/GenABEL/CHANGES.LOG
pkg/GenABEL/R/hom.R
pkg/GenABEL/R/polygenic.R
pkg/GenABEL/R/qtscore.R
pkg/GenABEL/man/polygenic.Rd
pkg/GenABEL/man/qtscore.Rd
Log:
* patched Roxygen documentation for qtscore
* use of 'setTxtProgressBar' in qtscore with times > 1 to indicate progress
* polygenic() added default option patchBasedOnFGLS; if no convergence based on difFGLS criterion, fixed effect betas are patched to FGLS betas before starting re-iteration
Modified: pkg/GenABEL/CHANGES.LOG
===================================================================
--- pkg/GenABEL/CHANGES.LOG 2010-11-19 05:55:19 UTC (rev 622)
+++ pkg/GenABEL/CHANGES.LOG 2010-12-07 21:42:43 UTC (rev 623)
@@ -1,5 +1,13 @@
-***** v. 1.6-5 (2010.10.16)
+***** v. 1.6-5 (2010.12.07)
+patched Roxygen documentation for qtscore
+
+use of 'setTxtProgressBar' in qtscore with times > 1 to indicate progress
+
+polygenic() added default option patchBasedOnFGLS, if no convergence based
+on difFGLS criterion, fixed effect betas are patched to FGLS betas before
+starting re-iteration
+
adding a wrapper chi2_CG -> cocohet
fixed 'demo(ge03d2)'
Modified: pkg/GenABEL/R/hom.R
===================================================================
--- pkg/GenABEL/R/hom.R 2010-11-19 05:55:19 UTC (rev 622)
+++ pkg/GenABEL/R/hom.R 2010-12-07 21:42:43 UTC (rev 623)
@@ -1,41 +1,41 @@
"hom" <-
-function(data,snpsubset,idsubset,snpfreq,n.snpfreq=1000) {
- if (is(data,"gwaa.data")) {
- data <- data at gtdata
- }
- if (!is(data,"snp.data")) {
- stop("data argument should have gwaa.data-class or snp.data-class")
- }
-
- if (!missing(snpsubset)) data <- data[,snpsubset]
- useFreq <- 0
- if (!missing(snpfreq)) {
- if (length(snpfreq) != data at nsnps) stop("snpfreq argument not equal in length to the number of SNPs in data")
- if (any(snpfreq<0) || any(snpfreq>1.)) stop("snpfreq argument: frequencies out of [0,1]")
- if (!is(snpfreq,"numeric")) stop("snpfreq argument: non-numeric class")
- useFreq <- 1
- } else {
- snpfreq <- rep(-1.,data at nsnps)
- }
- if (!missing(idsubset)) data <- data[idsubset,]
-
- totsnps <- data at nsnps
- totids <- data at nids
- if (!is(n.snpfreq,"numeric") || min(n.snpfreq)<0) stop("n.snpfreq should be positive numeric")
- if (length(n.snpfreq)==1) {
+ function(data,snpsubset,idsubset,snpfreq,n.snpfreq=1000) {
+ if (is(data,"gwaa.data")) {
+ data <- data at gtdata
+ }
+ if (!is(data,"snp.data")) {
+ stop("data argument should have gwaa.data-class or snp.data-class")
+ }
+
+ if (!missing(snpsubset)) data <- data[,snpsubset]
+ useFreq <- 0
+ if (!missing(snpfreq)) {
+ if (length(snpfreq) != data at nsnps) stop("snpfreq argument not equal in length to the number of SNPs in data")
+ if (any(snpfreq<0) || any(snpfreq>1.)) stop("snpfreq argument: frequencies out of [0,1]")
+ if (!is(snpfreq,"numeric")) stop("snpfreq argument: non-numeric class")
+ useFreq <- 1
+ } else {
+ snpfreq <- rep(-1.,data at nsnps)
+ }
+ if (!missing(idsubset)) data <- data[idsubset,]
+
+ totsnps <- data at nsnps
+ totids <- data at nids
+ if (!is(n.snpfreq,"numeric") || min(n.snpfreq)<0) stop("n.snpfreq should be positive numeric")
+ if (length(n.snpfreq)==1) {
n.snpfreq <- rep(n.snpfreq,data at nsnps)
} else {
if (length(n.snpfreq) != data at nsnps) stop("length mismatch in n.snpfreq")
}
- out <- .C("hom",as.raw(data at gtps),as.integer(data at nids),as.integer(data at nsnps),as.double(snpfreq),as.double(n.snpfreq),as.integer(useFreq),sout = double(5*data at nids), PACKAGE="GenABEL")$sout
- dim(out) <- c(data at nids,5)
- F <- (out[,3]-out[,4])/(out[,1]-out[,4])
- out <- cbind(out,F)
- out[,3] <- out[,3]/out[,1]
- out[,4] <- out[,4]/out[,1]
- out[,5] <- out[,5]/out[,2]
- colnames(out) <- c("NoMeasured","NoPoly","Hom","E(Hom)","Var","F")
- out <- as.data.frame(out,stringsAsFactors=FALSE)
- rownames(out) <- data at idnames
- out
+ out <- .C("hom",as.raw(data at gtps),as.integer(data at nids),as.integer(data at nsnps),as.double(snpfreq),as.double(n.snpfreq),as.integer(useFreq),sout = double(5*data at nids), PACKAGE="GenABEL")$sout
+ dim(out) <- c(data at nids,5)
+ F <- (out[,3]-out[,4])/(out[,1]-out[,4])
+ out <- cbind(out,F)
+ out[,3] <- out[,3]/out[,1]
+ out[,4] <- out[,4]/out[,1]
+ out[,5] <- out[,5]/out[,2]
+ colnames(out) <- c("NoMeasured","NoPoly","Hom","E(Hom)","Var","F")
+ out <- as.data.frame(out,stringsAsFactors=FALSE)
+ rownames(out) <- data at idnames
+ out
}
Modified: pkg/GenABEL/R/polygenic.R
===================================================================
--- pkg/GenABEL/R/polygenic.R 2010-11-19 05:55:19 UTC (rev 622)
+++ pkg/GenABEL/R/polygenic.R 2010-12-07 21:42:43 UTC (rev 623)
@@ -1,108 +1,110 @@
#' Estimation of polygenic model
#'
-#' This function maximises the likelihood of the data under polygenic
-#' model with covariates an reports twice negative maximum likelihood estimates
-#' and the inverse of variance-covariance matrix at the point of ML.
+#' This function maximises the likelihood of the data under polygenic
+#' model with covariates an reports twice negative maximum likelihood estimates
+#' and the inverse of variance-covariance matrix at the point of ML.
#'
-#' One of the major use of this function is to estimate residuals of the
-#' trait and the inverse of the variance-covariance matrix for
-#' further use in analysis with \code{\link{mmscore}} and
-#' \code{\link{grammar}}.
+#' One of the major use of this function is to estimate residuals of the
+#' trait and the inverse of the variance-covariance matrix for
+#' further use in analysis with \code{\link{mmscore}} and
+#' \code{\link{grammar}}.
#'
-#' Also, it can be used for a variant of GRAMMAR analysis, which
-#' allows for permutations for GW significance by use of
-#' environmental residuals as an analysis trait with \code{\link{qtscore}}.
+#' Also, it can be used for a variant of GRAMMAR analysis, which
+#' allows for permutations for GW significance by use of
+#' environmental residuals as an analysis trait with \code{\link{qtscore}}.
#'
-#' "Environmental residuals" (not to be mistaken with just "residuals") are
-#' the residual where both the effect of covariates AND the estimated
-#' polygenic effect (breeding values) are factored out. This thus
-#' provides an estimate of the trait value contributed by environment
-#' (or, turning this other way around, the part of trait not explained
-#' by covariates and by the polygene). Polygenic residuals are estimated
-#' as
+#' "Environmental residuals" (not to be mistaken with just "residuals") are
+#' the residual where both the effect of covariates AND the estimated
+#' polygenic effect (breeding values) are factored out. This thus
+#' provides an estimate of the trait value contributed by environment
+#' (or, turning this other way around, the part of trait not explained
+#' by covariates and by the polygene). Polygenic residuals are estimated
+#' as
#'
-#' \deqn{
-#' \sigma^2 V^{-1} (Y - (\hat{\mu} + \hat{\beta} C_1 + ...))
-#' }
+#' \deqn{
+#' \sigma^2 V^{-1} (Y - (\hat{\mu} + \hat{\beta} C_1 + ...))
+#' }
#'
-#' where \eqn{sigma^2} is the residual variance, \eqn{V^{-1}} is the
-#' InvSigma (inverse of the var-cov matrix at the maximum of
-#' polygenic model) and
-#' \eqn{(Y - (\hat{\mu} + \hat{\beta} C_1 + ...))} is the trait
-#' values adjusted for covariates (also at at the maximum of
-#' polygenic model likelihood).
+#' where \eqn{sigma^2} is the residual variance, \eqn{V^{-1}} is the
+#' InvSigma (inverse of the var-cov matrix at the maximum of
+#' polygenic model) and
+#' \eqn{(Y - (\hat{\mu} + \hat{\beta} C_1 + ...))} is the trait
+#' values adjusted for covariates (also at at the maximum of
+#' polygenic model likelihood).
#'
-#' It can also be used for heritability analysis.
-#' If you want to test significance of heritability,
-#' estimate the model and write down
-#' the function minimum reported at "h2an" element of the output
-#' (this is twice negative MaxLikleihood). Then do next round of
-#' estimation, but set fixh2=0. The difference between you function minima
-#' gives a test distribued as chi-squared with 1 d.f.
+#' It can also be used for heritability analysis.
+#' If you want to test significance of heritability,
+#' estimate the model and write down
+#' the function minimum reported at "h2an" element of the output
+#' (this is twice negative MaxLikleihood). Then do next round of
+#' estimation, but set fixh2=0. The difference between you function minima
+#' gives a test distribued as chi-squared with 1 d.f.
#'
-#' The way to compute the likleihood is partly based on
-#' the paper of Thompson (see refs), namely instead of
-#' taking inverse of var-cov matrix every time,
-#' eigenvectors of the inverse of G (taken only once)
-#' are used.
+#' The way to compute the likleihood is partly based on
+#' the paper of Thompson (see refs), namely instead of
+#' taking inverse of var-cov matrix every time,
+#' eigenvectors of the inverse of G (taken only once)
+#' are used.
#'
#'
-#' @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 fixh2 Optional value of heritability to be used, instead of maximisation.
-#' The uses of this option are two-fold: (a) testing significance of
-#' heritability and (b) using a priori known heritability to derive the
-#' rest of MLEs and var.-cov. matrix.
-#' @param starth2 Starting value for h2 estimate
-#' @param trait.type "gaussian" or "binomial"
-#' @param opt.method "nlm" or "optim". These two use dirrerent optimisation functions.
-#' We suggest using the default \code{\link{nlm}}, though
-#' \code{\link{optim}} may give better results in some situations
-#' @param scaleh2 Only relevant when "nlm" optimisation function is used.
-#' "scaleh2" is the heritability
-#' scaling parameter, regulating how "big" are parameter changes in h2 with the
-#' respect to changes in other parameters. As other parameters are estimated
-#' from previous regression, these are expected to change little from the
-#' initial estimate. The default value of 1000 proved to work rather well under a
-#' range of conditions.
-#' @param quiet If FALSE (default), details of optimisation process are reported.
-#' @param steptol steptal parameter of "nlm"
-#' @param gradtol gradtol parameter of "nlm"
-#' @param optimbou fixed effects boundary scale parameter for 'optim'
-#' @param fglschecks additional check for convergance on/off (convergence
-#' between estimates obtained and that from FGLS)
-#' @param maxnfgls number of fgls checks to perform
-#' @param maxdiffgls max difference allowed in fgls checks
-#' @param ... Optional arguments to be passed to \code{\link{nlm}} or (\code{\link{optim}})
-#' minimisation function
+#' @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 fixh2 Optional value of heritability to be used, instead of maximisation.
+#' The uses of this option are two-fold: (a) testing significance of
+#' heritability and (b) using a priori known heritability to derive the
+#' rest of MLEs and var.-cov. matrix.
+#' @param starth2 Starting value for h2 estimate
+#' @param trait.type "gaussian" or "binomial"
+#' @param opt.method "nlm" or "optim". These two use dirrerent optimisation functions.
+#' We suggest using the default \code{\link{nlm}}, though
+#' \code{\link{optim}} may give better results in some situations
+#' @param scaleh2 Only relevant when "nlm" optimisation function is used.
+#' "scaleh2" is the heritability
+#' scaling parameter, regulating how "big" are parameter changes in h2 with the
+#' respect to changes in other parameters. As other parameters are estimated
+#' from previous regression, these are expected to change little from the
+#' initial estimate. The default value of 1000 proved to work rather well under a
+#' range of conditions.
+#' @param quiet If FALSE (default), details of optimisation process are reported.
+#' @param steptol steptal parameter of "nlm"
+#' @param gradtol gradtol parameter of "nlm"
+#' @param optimbou fixed effects boundary scale parameter for 'optim'
+#' @param fglschecks additional check for convergance on/off (convergence
+#' between estimates obtained and that from FGLS)
+#' @param maxnfgls number of fgls checks to perform
+#' @param maxdiffgls max difference allowed in fgls checks
+#' @param patchBasedOnFGLS if FGLS checks not passed, 'patch' fixed
+#' effect estimates based on FGLS expectation
+#' @param ... Optional arguments to be passed to \code{\link{nlm}} or (\code{\link{optim}})
+#' minimisation function
#'
-#' @return
-#' A list with values
-#' \item{h2an}{A list supplied by the \code{\link{nlm}} minimisation routine.
-#' Of particular interest are elements "estimate" containing parameter
-#' maximal likelihood estimates (MLEs) (order: mean, betas for covariates,
-#' heritability, (polygenic + residual variance)). The value of
-#' twice negative maximum log-likelihood
-#' is returned as h2an\$minimum.}
-#' \item{residualY}{Residuals from analysis, based on covariate effects only;
-#' NOTE: these are NOT grammar "environmental residuals"!}
-#' \item{esth2}{Estimate (or fixed value) of heritability}
-#' \item{pgresidualY}{Environmental residuals from analysis, based on covariate effects
-#' and predicted breeding value.
-#' }
-#' \item{InvSigma}{Inverse of the variance-covariance matrix, computed at the
-#' MLEs -- these are used in \code{\link{mmscore}} and \code{\link{grammar}}
-#' functions.}
-#' \item{call}{The details of call}
-#' \item{measuredIDs}{Logical values for IDs who were used in analysis
-#' (traits and all covariates measured) == TRUE}
-#' \item{convFGLS}{was convergence achieved according to FGLS criterionE}
+#' @return
+#' A list with values
+#' \item{h2an}{A list supplied by the \code{\link{nlm}} minimisation routine.
+#' Of particular interest are elements "estimate" containing parameter
+#' maximal likelihood estimates (MLEs) (order: mean, betas for covariates,
+#' heritability, (polygenic + residual variance)). The value of
+#' twice negative maximum log-likelihood
+#' is returned as h2an\$minimum.}
+#' \item{residualY}{Residuals from analysis, based on covariate effects only;
+#' NOTE: these are NOT grammar "environmental residuals"!}
+#' \item{esth2}{Estimate (or fixed value) of heritability}
+#' \item{pgresidualY}{Environmental residuals from analysis, based on covariate effects
+#' and predicted breeding value.
+#' }
+#' \item{InvSigma}{Inverse of the variance-covariance matrix, computed at the
+#' MLEs -- these are used in \code{\link{mmscore}} and \code{\link{grammar}}
+#' functions.}
+#' \item{call}{The details of call}
+#' \item{measuredIDs}{Logical values for IDs who were used in analysis
+#' (traits and all covariates measured) == TRUE}
+#' \item{convFGLS}{was convergence achieved according to FGLS criterionE}
#'
#'
#' @references
@@ -116,16 +118,16 @@
#'
#' 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
#'
#' @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.
+#' 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.
#'
-#' If a trait (no covarites) is used, make sure that order of IDs in
-#' kinship.matrix is exactly the same as in the outcome
+#' If a trait (no covarites) is used, make sure that order of IDs in
+#' kinship.matrix is exactly the same as in the outcome
#'
#' @seealso
#' \code{\link{mmscore}},
@@ -157,7 +159,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, ...) {
+ fglschecks=TRUE,maxnfgls=8,maxdiffgls=1e-4, patchBasedOnFGLS = TRUE, ...) {
if (!missing(data)) if (is(data,"gwaa.data"))
{
checkphengen(data)
@@ -225,7 +227,7 @@
inierr <- sqrt(var(y)/length(y))
}
}
-
+
if (!missing(data)) detach(data)
tmp <- t(relmat)
relmat[upper.tri(relmat)] <- tmp[upper.tri(tmp)]
@@ -320,6 +322,10 @@
if ((h2<1e-4 || h2>(1-1e-4)) && nfgls > floor(maxnfgls/2)) {
parsave[npar-1] <- runif(1,min=0.05,max=0.95)/scaleh2
}
+ if (patchBasedOnFGLS && diffgls > maxdiffgls) {
+ if (!quiet) {cat("fixed effect betas changed to FGLS-betas for re-estimation\n")}
+ parsave[1:(npar-2)] <- betaFGLS;
+ }
} else {
nfgls <- maxnfgls+1
}
Modified: pkg/GenABEL/R/qtscore.R
===================================================================
--- pkg/GenABEL/R/qtscore.R 2010-11-19 05:55:19 UTC (rev 622)
+++ pkg/GenABEL/R/qtscore.R 2010-12-07 21:42:43 UTC (rev 623)
@@ -22,17 +22,29 @@
#' 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 data
-#' @param snpsubset
-#' @param idsubset
-#' @param strata
-#' @param trait.type
-#' @param times
-#' @param quiet
-#' @param bcast
-#' @param clambda
-#' @param propPs
-#' @param details
+#' @param data An object of \code{\link{gwaa.data-class}}
+#' @param snpsubset ndex, character or logical vector with subset of SNPs to run analysis on.
+#' If missing, all SNPs from \code{data} are used for analysis.
+#' @param idsubset ndex, character or logical vector with subset of IDs to run analysis on.
+#' If missing, all people from \code{data/cc} are used for analysis.
+#' @param strata Stratification variable. If provieded, scores are computed within strata and
+#' then added up.
+#' @param trait.type "gaussian" or "binomial" or "guess" (later option guesses trait type)
+#' @param times If more then one, the number of replicas to be used in derivation of
+#' empirical genome-wide significance. See \code{\link{emp.qtscore}}, which
+#' calls qtscore with times>1 for details
+#' @param quiet do not print warning messages
+#' @param bcast If the argument times > 1, progress is reported once in bcast replicas
+#' @param clambda If inflation facot Lambda is estimated as lower then one, this parameter
+#' controls if the original P1df (clambda=TRUE) to be reported in Pc1df,
+#' or the original 1df statistics is to be multiplied onto this "deflation"
+#' factor (clambda=FALSE). If a numeric value is provided, it is used as a correction factor.
+#' @param propPs proportion of non-corrected P-values used to estimate the inflation factor Lambda,
+#' passed directly to the \code{\link{estlambda}}
+#' @param details when FALSE, SNP and ID names are not reported in the returned object
+#' (saves some memory). This is experimental and will be not mantained anymore
+#' as soon as we achieve better memory efficiency for storage of SNP and ID names
+#' (currently default R character data type used)
#'
#' @return Object of class \code{\link{scan.gwaa-class}}
#'
@@ -159,6 +171,7 @@
if (trait.type=="binomial") bin<-1 else bin<-0
lenn <- data at gtdata@nsnps;
###out <- list()
+ if (times>1) {pb <- txtProgressBar(style = 3)}
for (j in c(1:(times+1*(times>1)))) {
if (j>1) resid <- sample(resid,replace=FALSE)
chi2 <- .C("qtscore_glob",as.raw(data at gtdata@gtps),as.double(resid),as.integer(bin),as.integer(data at gtdata@nids),as.integer(data at gtdata@nsnps), as.integer(nstra), as.integer(strata), chi2 = double(10*data at gtdata@nsnps), PACKAGE="GenABEL")$chi2
@@ -255,14 +268,17 @@
pr.2df <- pr.2df + 1*(chi2.2df < max(chi2[(lenn+1):(2*lenn)]))
pr.c1df <- pr.c1df + 1*(chi2.c1df < th1)
pr.c2df <- pr.c2df + 1*(chi2.c2df < th1)
+# if (!quiet && ((j-1)/bcast == round((j-1)/bcast))) {
+ ## cat("\b\b\b\b\b\b",round((100*(j-1)/times),digits=2),"%",sep="")
+# cat(" ",round((100*(j-1)/times),digits=2),"%",sep="")
+# flush.console()
+# }
if (!quiet && ((j-1)/bcast == round((j-1)/bcast))) {
-# cat("\b\b\b\b\b\b",round((100*(j-1)/times),digits=2),"%",sep="")
- cat(" ",round((100*(j-1)/times),digits=2),"%",sep="")
- flush.console()
+ setTxtProgressBar(pb, (j-1)/times)
}
}
}
- if (times > bcast) cat("\n")
+ if (times > bcast) {setTxtProgressBar(pb, 1.0);cat("\n")}
if (times>1) {
P1df <- pr.1df/times
Modified: pkg/GenABEL/man/polygenic.Rd
===================================================================
--- pkg/GenABEL/man/polygenic.Rd 2010-11-19 05:55:19 UTC (rev 622)
+++ pkg/GenABEL/man/polygenic.Rd 2010-12-07 21:42:43 UTC (rev 623)
@@ -4,7 +4,7 @@
\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, ...)}
+ 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
@@ -28,7 +28,7 @@
as
\deqn{
-\sigma^2 V^{-1} (Y - (\hat{\mu} + \hat{\beta} C_1 + ...))
+\sigma^2 V^{-1} (Y - (\hat{\mu} + \hat{\beta} C_1 + ...))
}
where \eqn{sigma^2} is the residual variance, \eqn{V^{-1}} is the
@@ -123,6 +123,8 @@
between estimates obtained and that from FGLS)}
\item{maxnfgls}{number of fgls checks to perform}
\item{maxdiffgls}{max difference allowed in fgls checks}
+\item{patchBasedOnFGLS}{if FGLS checks not passed, 'patch' fixed
+effect estimates based on FGLS expectation}
\item{...}{Optional arguments to be passed to \code{\link{nlm}} or (\code{\link{optim}})
minimisation function}}
\examples{# note that procedure runs on CLEAN data
Modified: pkg/GenABEL/man/qtscore.Rd
===================================================================
--- pkg/GenABEL/man/qtscore.Rd 2010-11-19 05:55:19 UTC (rev 622)
+++ pkg/GenABEL/man/qtscore.Rd 2010-12-07 21:42:43 UTC (rev 623)
@@ -1,44 +1,13 @@
\name{qtscore}
\alias{qtscore}
-\title{Fast score test for association}
-\description{
-Fast score test for association between a trait and genetic polymorphism
-}
-\usage{
-qtscore(formula,data,snpsubset,idsubset,strata,trait.type="gaussian",times=1,quiet=FALSE,bcast=10,clambda=TRUE,propPs=1.0,details=TRUE)
-}
-\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{data}{An object of \code{\link{gwaa.data-class}}}
- \item{snpsubset}{Index, character or logical vector with subset of SNPs to run analysis on.
- If missing, all SNPs from \code{data} are used for analysis.}
- \item{idsubset}{Index, character or logical vector with subset of IDs to run analysis on.
- If missing, all people from \code{data/cc} are used for analysis.}
- \item{strata}{Stratification variable. If provieded, scores are computed within strata and
- then added up.}
- \item{trait.type}{"gaussian" or "binomial" or "guess" (later option guesses trait type)}
- \item{times}{If more then one, the number of replicas to be used in derivation of
- empirical genome-wide significance. See \code{\link{emp.qtscore}}, which
- calls qtscore with times>1 for details}
- \item{quiet}{do not print warning messages}
- \item{bcast}{If the argument times > 1, progress is reported once in bcast replicas}
- \item{clambda}{If inflation facot Lambda is estimated as lower then one, this parameter
- controls if the original P1df (clambda=TRUE) to be reported in Pc1df,
- or the original 1df statistics is to be multiplied onto this "deflation"
- factor (clambda=FALSE).
- If a numeric value is provided, it is used as a correction factor.}
- \item{propPs}{proportion of non-corrected P-values used to estimate the inflation factor Lambda,
- passed directly to the \code{\link{estlambda}}}
- \item{details}{when FALSE, SNP and ID names are not reported in the returned object
- (saves some memory). This is experimental and will be not mantained anymore
- as soon as we achieve better memory efficiency for storage of SNP and ID names
- (currently default R character data type used)}
-}
-\details{
+\title{Fast score test for association...}
+\usage{qtscore(formula, data, snpsubset, idsubset, strata,
+ trait.type="gaussian", times=1, quiet=FALSE, bcast=10,
+ clambda=TRUE, propPs=1, details=TRUE)}
+\description{Fast score test for association}
+\details{Fast score test for association
+between a trait and genetic polymorphism
+
When formula contains covariates, the traits is analysed using GLM and later
residuals used when score test is computed for each of the SNPs in analysis.
Coefficients of regression are reported for the quantitative trait.
@@ -52,31 +21,50 @@
Armitage test.
This is a valid function to analyse GWA data, including X chromosome. For X chromosome,
-stratified analysis is performed (strata=sex).
-}
-\value{
- Object of class \code{\link{scan.gwaa-class}}
-}
-\references{
-Aulchenko YS, de Koning DJ, Haley C. Genomewide rapid association using mixed model
+stratified analysis is performed (strata=sex).}
+\value{Object of class \code{\link{scan.gwaa-class}}}
+\references{Aulchenko YS, de Koning DJ, Haley C. Genomewide rapid association using mixed model
and regression: a fast and simple method for genome-wide pedigree-based quantitative
trait loci association analysis. Genetics. 2007 177(1):577-85.
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.
-}
+association analysis in related individuals. PLoS ONE. 2007 Dec 5;2(12):e1274.}
\author{Yurii Aulchenko}
-%\note{}
-\seealso{
-\code{\link{mlreg}},
+\seealso{\code{\link{mlreg}},
\code{\link{mmscore}},
\code{\link{egscore}},
\code{\link{emp.qtscore}},
\code{\link{plot.scan.gwaa}},
-\code{\link{scan.gwaa-class}}
-}
-\examples{
-data(srdta)
+\code{\link{scan.gwaa-class}}}
+\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{data}{An object of \code{\link{gwaa.data-class}}}
+\item{snpsubset}{ndex, character or logical vector with subset of SNPs to run analysis on.
+If missing, all SNPs from \code{data} are used for analysis.}
+\item{idsubset}{ndex, character or logical vector with subset of IDs to run analysis on.
+If missing, all people from \code{data/cc} are used for analysis.}
+\item{strata}{Stratification variable. If provieded, scores are computed within strata and
+then added up.}
+\item{trait.type}{"gaussian" or "binomial" or "guess" (later option guesses trait type)}
+\item{times}{If more then one, the number of replicas to be used in derivation of
+empirical genome-wide significance. See \code{\link{emp.qtscore}}, which
+calls qtscore with times>1 for details}
+\item{quiet}{do not print warning messages}
+\item{bcast}{If the argument times > 1, progress is reported once in bcast replicas}
+\item{clambda}{If inflation facot Lambda is estimated as lower then one, this parameter
+controls if the original P1df (clambda=TRUE) to be reported in Pc1df,
+or the original 1df statistics is to be multiplied onto this "deflation"
+factor (clambda=FALSE). If a numeric value is provided, it is used as a correction factor.}
+\item{propPs}{proportion of non-corrected P-values used to estimate the inflation factor Lambda,
+passed directly to the \code{\link{estlambda}}}
+\item{details}{when FALSE, SNP and ID names are not reported in the returned object
+(saves some memory). This is experimental and will be not mantained anymore
+as soon as we achieve better memory efficiency for storage of SNP and ID names
+(currently default R character data type used)}}
+\examples{data(srdta)
#qtscore with stratification
a <- qtscore(qt3~sex,data=srdta)
plot(a)
@@ -93,6 +81,4 @@
add.plot(a2,col="red",cex=1.5)
# the good thing about score test is that we can do adjustment...
a2 <- qtscore(bt~age+sex,data=srdta,snps=c(1:100),trait.type="binomial")
-points(a2[,"Position"],-log10(a2[,"P1df"]),col="green")
-}
-\keyword{htest}% at least one, from doc/KEYWORDS
+points(a2[,"Position"],-log10(a2[,"P1df"]),col="green")}
More information about the Genabel-commits
mailing list