[Genabel-commits] r1331 - pkg/GenABEL/R

[noreply at r-forge.r-project.org]

Author: lckarssen
Date: 2013-09-06 15:22:24 +0200 (Fri, 06 Sep 2013)
New Revision: 1331

Modified:
   pkg/GenABEL/R/estlambda.R
Log:
In GenABEL's estlambda.R: 
- update the help text and fix some typo's
- in the plot: add \chi^2 to the ylab and xlab (people may expect p-values on the axes if they send in p-value data, so this is to remove confusion.)
- some other small code layout changes.


Modified: pkg/GenABEL/R/estlambda.R
===================================================================
--- pkg/GenABEL/R/estlambda.R	2013-09-06 12:47:19 UTC (rev 1330)
+++ pkg/GenABEL/R/estlambda.R	2013-09-06 13:22:24 UTC (rev 1331)
@@ -10,16 +10,18 @@
 #'
 #' @param data A vector of reals. If all are <=1, it is assumed that this is a
 #'   vector of P-values, else it is treated as a vector of chi-squares
-#' @param plot Wether the prot should be presented
-#' @param proportion The proportion of lowest P (Chi2) to be used when
-#'   estimating the inflation factor Lambda
-#' @param method "regression", "median", or "KS" method to be used for
-#'   Lambda estimation
-#' @param filter if the test statistics with 0-value of chi-square should be
-#'   excluded prior to estimation of Lambda
-#' @param df Number of degrees of freedom
-#' @param ... arguments passed to \code{\link{plot}} function
-#' @return A list with elements \item{estimate}{Estimate of Lambda}
+#' @param plot Whether the plot should be shown or not (default).
+#' @param proportion The proportion of lowest P (or
+#'   \eqn{\chi^2}{chi^2) values to be used when estimating the inflation
+#'   factor \eqn{\lambda}{lambda}.
+#' @param method "regression" (default), "median", or "KS": method to
+#'   be used for \eqn{\lambda}{lambda} estimation.
+#' @param filter if the test statistics with 0-value of
+#'   \eqn{\chi^2}{chi^2} should be excluded prior to estimation of
+#'   \eqn{\lambda}{lambda}.
+#' @param df Number of degrees of freedom.
+#' @param ... arguments passed to the \code{\link{plot}} function.
+#' @return A list with elements \item{estimate}{Estimate of \eqn{\lambda}{lambda}}
 #'   \item{se}{Standard error of the estimate}
 #' @author Yurii Aulchenko
 #' @seealso \code{\link{ccfast}}, \code{\link{qtscore}}
@@ -35,26 +37,27 @@
 #' a <- qtscore(bt,srdta)
 #' lambda(a)
 #'
-"estlambda" <-
-                function(data,plot = FALSE, proportion=1.0, method="regression", filter=TRUE, df=1,... ) {
+"estlambda" <- function(data, plot=FALSE, proportion=1.0,
+                        method="regression", filter=TRUE, df=1,... ) {
         data <- data[which(!is.na(data))]
         if (proportion>1.0 || proportion<=0)
                 stop("proportion argument should be greater then zero and less than or equal to one")
-        ntp <- round(proportion*length(data))
-        if (ntp<1) stop("no valid measurments")
-        if (ntp==1) {
-                warning(paste("One measurment, Lambda = 1 returned"))
-                return(list(estimate=1.0,se=999.99))
+
+        ntp <- round( proportion * length(data) )
+        if ( ntp<1 ) stop("no valid measurements")
+        if ( ntp==1 ) {
+                warning(paste("One measurement, lambda = 1 returned"))
+                return(list(estimate=1.0, se=999.99))
         }
-        if (ntp<10) warning(paste("number of points is too small:",ntp))
-        if (min(data)<0) stop("data argument has values <0")
-        if (max(data)<=1) {
+        if ( ntp<10 ) warning(paste("number of points is too small:", ntp))
+        if ( min(data)<0 ) stop("data argument has values <0")
+        if ( max(data)<=1 ) {
 #		lt16 <- (data < 1.e-16)
 #		if (any(lt16)) {
 #			warning(paste("Some probabilities < 1.e-16; set to 1.e-16"))
 #			data[lt16] <- 1.e-16
 #		}
-                data <- qchisq(data,1,lower.tail=FALSE)
+                data <- qchisq(data, 1, lower.tail=FALSE)
         }
         if (filter)
         {
@@ -62,7 +65,7 @@
         }
         data <- sort(data)
         ppoi <- ppoints(data)
-        ppoi <- sort(qchisq(ppoi,df=df,lower.tail=FALSE))
+        ppoi <- sort(qchisq(ppoi, df=df, lower.tail=FALSE))
         data <- data[1:ntp]
         ppoi <- ppoi[1:ntp]
 #	s <- summary(lm(data~offset(ppoi)))$coeff
@@ -70,15 +73,15 @@
 
         out <- list()
         if (method=="regression") {
-                s <- summary(lm(data~0+ppoi))$coeff
+                s <- summary( lm(data~0+ppoi) )$coeff
                 out$estimate <- s[1,1]
                 out$se <- s[1,2]
         } else if (method=="median") {
-                out$estimate <- median(data,na.rm=TRUE)/qchisq(0.5,df)
+                out$estimate <- median(data, na.rm=TRUE)/qchisq(0.5, df)
                 out$se <- NA
         } else if (method=="KS") {
-                limits <- c(0.5,100)
-                out$estimate <- estLambdaKS(data,limits=limits,df=df)
+                limits <- c(0.5, 100)
+                out$estimate <- estLambdaKS(data, limits=limits, df=df)
                 if ( abs(out$estimate-limits[1])<1e-4 || abs(out$estimate-limits[2])<1e-4 )
                         warning("using method='KS' lambda too close to limits, use other method")
                 out$se <- NA
@@ -87,11 +90,17 @@
         }
 
         if (plot) {
-                lim <- c(0,max(data,ppoi,na.rm=T))
+                lim <- c(0, max(data, ppoi,na.rm=TRUE))
 #		plot(ppoi,data,xlim=lim,ylim=lim,xlab="Expected",ylab="Observed", ...)
-                plot(ppoi,data,xlab="Expected",ylab="Observed", ...)
-                abline(a=0,b=1)
-                abline(a=0,b=out$estimate,col="red")
+                oldmargins <- par()$mar
+                par(mar=oldmargins + 0.2)
+                plot(ppoi, data,
+                     xlab=expression("Expected " ~ chi^2),
+                     ylab=expression("Observed " ~ chi^2),
+                     ...)
+                abline(a=0, b=1)
+                abline(a=0, b=out$estimate, col="red")
+                par(mar=oldmargins)
         }
 
         out