[Genabel-commits] r1564 - pkg/MetABEL/R

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

Author: lckarssen
Date: 2014-01-26 23:28:08 +0100 (Sun, 26 Jan 2014)
New Revision: 1564

Modified:
   pkg/MetABEL/R/forestplot.R
Log:
Improved code layout of forestplot.R in MetABEL, brought it more in line with our coding standards.
Only whitespace changes, no functional changes.


Modified: pkg/MetABEL/R/forestplot.R
===================================================================
--- pkg/MetABEL/R/forestplot.R	2014-01-26 22:02:28 UTC (rev 1563)
+++ pkg/MetABEL/R/forestplot.R	2014-01-26 22:28:08 UTC (rev 1564)
@@ -1,55 +1,83 @@
-"forestplot" <- 
-function(estimate,se,labels=paste("Study",c(1:length(estimate))),CI=0.95,xexp=FALSE,...) {
-	hoff <- 3
-	del <- 10
-	mea <- !is.na(estimate)
-	estimate <- estimate[mea]
-	se <- se[mea]
-	labels <- labels[mea]
-	w2 <- 1./(se*se)
-	invsumw2 <- 1./sum(w2)
-	mestimate <- sum(estimate*w2)*invsumw2
-	mse <- sqrt(invsumw2)
-	npop <- length(estimate)
-	estimate[npop+1] <- mestimate
-	se[npop+1] <- mse
-	labels[npop+1] <- "Pooled"
-	chi2 <- round(estimate*estimate/(se*se),2)
-	p <- sprintf("%5.1e",pchisq(estimate*estimate/(se*se),1,lower=F))
-#	p[as.numeric(p)<0] <- "<1.e-16"
-	if (CI>1 || CI<0) stop("CI argument should be between 0 and 1")
-	cimultip <- qnorm(1-(1-CI)/2)
-	lower <- estimate - cimultip*se
-	upper <- estimate + cimultip*se
-	if (xexp) {
-		estimate <- exp(estimate)
-		lower <- exp(lower)
-		upper <- exp(upper)
-	}
-	cntr <- 0; if (xexp) cntr <- 1;
-	lbnd <- (-.1); if (xexp) lbnd <- 0.9
-	rbnd <- (.1); if (xexp) rbnd <- 1.1
-	minv <- min(lower)
-	minv <- minv-abs(minv/10)
-	minv <- min(lbnd,minv)
-	maxv <- max(upper)
-	maxv <- maxv+abs(maxv/10)
-	maxv <- max(rbnd,maxv)
-	hgt <- (length(estimate)+1)*del
-	if (any(is.na(estimate))) stop("estimate got NAs")
-	if (any(is.na(se))) stop("se got NAs")
-	plot(x=c(cntr,cntr),y=c(0,hgt),xlim=c(minv,maxv),ylim=c(0,hgt),typ="l",lwd=2,xlab="Beta",...)
-	for (i in c(1:(length(estimate)-1))) {
-		points(x=c(lower[i],upper[i]),y=c((i)*del,(i)*del),typ="l",lwd=2)
-		points(x=c(estimate[i]),y=c((i)*del),pch=19,cex=1)
-		text(estimate[i],i*del+1,paste(labels[i]," (Chi2=",chi2[i],", P=",p[i],")",sep=""),pos=3,cex=.7)
-	}
-	for (i in c(length(estimate))) {
-		points(x=c(lower[i],estimate[i]),y=c((i)*del,(i)*del+hoff),typ="l",lwd=2)
-		points(x=c(estimate[i],upper[i]),y=c((i)*del+hoff,(i)*del),typ="l",lwd=2)
-		points(x=c(upper[i],estimate[i]),y=c((i)*del,(i)*del-hoff),typ="l",lwd=2)
-		points(x=c(lower[i],estimate[i]),y=c((i)*del,(i)*del-hoff),typ="l",lwd=2)
-		text(estimate[i],i*del+5,paste(labels[i]," (Chi2=",chi2[i],", P=",p[i],")",sep=""),pos=3,cex=1)
-	}
-}
+"forestplot" <-
+    function(estimate, se,
+             labels=paste("Study", c(1:length(estimate))),
+             CI=0.95, xexp=FALSE, ...) {
+        hoff      <- 3
+        del       <- 10
+        mea       <- !is.na(estimate)
+        estimate  <- estimate[mea]
+        se        <- se[mea]
+        labels    <- labels[mea]
+        w2        <- 1. / (se * se)
+        invsumw2  <- 1. / sum(w2)
+        mestimate <- sum(estimate * w2) * invsumw2
+        mse       <- sqrt(invsumw2)
+        npop      <- length(estimate)
+        estimate[npop+1] <- mestimate
+        se[npop+1]     <- mse
+        labels[npop+1] <- "Pooled"
+        chi2 <- round(estimate * estimate / (se * se), 2)
+        p <- sprintf("%5.1e", pchisq(estimate * estimate / (se * se),
+                                     1, lower=F))
+        ##	p[as.numeric(p)<0] <- "<1.e-16"
 
+        if (CI > 1 || CI < 0) {
+            stop("CI argument should be between 0 and 1")
+        }
+
+        cimultip <- qnorm(1 - (1 - CI) / 2)
+        lower    <- estimate - cimultip * se
+        upper    <- estimate + cimultip * se
+
+        if (xexp) {
+            estimate <- exp(estimate)
+            lower    <- exp(lower)
+            upper    <- exp(upper)
+        }
+
+        cntr <- 0; if (xexp) cntr <- 1;
+        lbnd <- (-.1); if (xexp) lbnd <- 0.9
+        rbnd <- (.1); if (xexp) rbnd <- 1.1
+        minv <- min(lower)
+        minv <- minv - abs(minv / 10)
+        minv <- min(lbnd, minv)
+        maxv <- max(upper)
+        maxv <- maxv + abs(maxv / 10)
+        maxv <- max(rbnd, maxv)
+        hgt <- (length(estimate) + 1) * del
+
+        if (any(is.na(estimate))) stop("estimate got NAs")
+        if (any(is.na(se))) stop("se got NAs")
+
+        plot(x=c(cntr, cntr), y=c(0, hgt), xlim=c(minv, maxv),
+             ylim=c(0, hgt), typ="l", lwd=2, xlab="Beta", ...)
+
+        for (i in c(1:(length(estimate)-1))) {
+            points(x=c(lower[i], upper[i]), y=c((i) * del, (i) * del),
+                   typ="l", lwd=2)
+            points(x=c(estimate[i]), y=c((i) * del), pch=19, cex=1)
+            text(estimate[i], i * del + 1,
+                 paste(labels[i], " (Chi2=", chi2[i], ", P=", p[i],
+                       ")", sep=""),
+                 pos=3, cex=.7)
+        }
+
+        for (i in c(length(estimate))) {
+            points(x=c(lower[i], estimate[i]),
+                   y=c((i) * del, (i) * del + hoff),
+                   typ="l", lwd=2)
+            points(x=c(estimate[i], upper[i]),
+                   y=c((i) * del + hoff, (i) * del),
+                   typ="l", lwd=2)
+            points(x=c(upper[i], estimate[i]),
+                   y=c((i) * del, (i) * del - hoff),
+                   typ="l", lwd=2)
+            points(x=c(lower[i], estimate[i]),
+                   y=c((i) * del, (i) * del - hoff),
+                   typ="l", lwd=2)
+            text(estimate[i], i * del + 5,
+                 paste(labels[i], " (Chi2=", chi2[i], ", P=", p[i],
+                       ")", sep=""),
+                 pos=3, cex=1)
+        }
+    }