From 10241adac7493d07e405516b224203fdf007a931 Mon Sep 17 00:00:00 2001
From: Christian Roever <christian.roever@med.uni-goettingen.de>
Date: Thu, 13 Oct 2016 08:29:34 +0200
Subject: [PATCH] fixed log-axis bug and added "prediction" argument for
forestplot() function
---
DESCRIPTION | 2 +-
R/bayesmeta.R | 125 ++++++++++++++++--------------------
man/bayesmeta-package.Rd | 2 +-
man/forestplot.bayesmeta.Rd | 32 +++++----
4 files changed, 77 insertions(+), 84 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 5933ba1..ad13260 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: bayesmeta
Type: Package
Title: Bayesian Random-Effects Meta-Analysis
Version: 1.3
-Date: 2016-10-08
+Date: 2016-10-13
Authors@R: c(person("Christian", "Roever", role=c("aut","cre"),
email="christian.roever@med.uni-goettingen.de"),
person("Tim", "Friede", role="ctb",
diff --git a/R/bayesmeta.R b/R/bayesmeta.R
index 0bcc5e9..aa75aa5 100644
--- a/R/bayesmeta.R
+++ b/R/bayesmeta.R
@@ -1534,22 +1534,23 @@ forest.bayesmeta <- function(x, xlab="effect size", refline=0, cex=1,...)
forestplot.bayesmeta <- function(x, labeltext,
exponentiate=FALSE,
+ prediction=TRUE,
shrinkage=TRUE,
digits=2,
plot=TRUE,
fn.ci_norm, col, legend,
boxsize = 0.25, ...)
-# forest plot for a "bayesmeta" object
-# based on the "forestplot" package's plotting function
#
# ARGUMENTS:
-# exponentiate : flag indicating whether to exponentiate table numbers.
+# x : a "bayesmeta" object.
+# labeltext : you may provide an alternative "labeltext" argument here
+# (see also the "forestplot()" help).
+# exponentiate : flag indicating whether to exponentiate numbers (figure and table).
+# prediction : flag indicating whether to show prediction interval.
# shrinkage : flag indicating whether to show shrinkage estimates.
# digits : number of significant digits to be shown (based on standard deviations).
-# labeltext : you may (try to) provide an alternative "labeltext" argument here
-# (see also the "forestplot" help).
# plot : flag you can use to suppress actual plotting.
-# ... : further arguments passes to the "forestplot" function.
+# ... : further arguments passed to the "forestplot" function.
#
# VALUE:
# a list with components
@@ -1563,25 +1564,31 @@ forestplot.bayesmeta <- function(x, labeltext,
stop("required 'forestplot' package not available!")
#stopifnot(require("forestplot"))
stopifnot(is.element("bayesmeta", class(x)),
- length(digits)==1, digits==round(digits), digits>=0)
+ length(digits)==1, digits==round(digits), digits>=0,
+ length(prediction)==1, is.logical(prediction),
+ length(shrinkage)==1, is.logical(shrinkage))
# the plotting data--
# the quoted estimates:
q95 <- qnorm(0.975)
ma.dat <- rbind(NA,
- cbind("mean"=exp(x$y),
- "lower"=exp(x$y - q95*x$sigma),
- "upper"=exp(x$y + q95*x$sigma)),
- exp(x$summary[c("median", "95% lower", "95% upper"),"mu"]),
- exp(x$summary[c("median", "95% lower", "95% upper"),"theta"]))
+ cbind(x$y, x$y - q95*x$sigma, x$y + q95*x$sigma),
+ x$summary[c("median", "95% lower", "95% upper"),"mu"],
+ x$summary[c("median", "95% lower", "95% upper"),"theta"])
colnames(ma.dat) <- c("estimate", "lower", "upper")
rownames(ma.dat) <- c("", x$label, "mean", "prediction")
+ if (! prediction) ma.dat <- ma.dat[-(x$k+3),]
# the shrinkage estimates:
ma.shrink <- rbind(NA,
- exp(t(x$theta)[,c("median","95% lower","95% upper")]),
+ t(x$theta)[,c("median","95% lower","95% upper")],
NA,NA)
colnames(ma.shrink) <- c("estimate", "lower", "upper")
rownames(ma.shrink) <- c("", x$label, "mean", "prediction")
- # generate data table (if not provided):
+ if (! prediction) ma.shrink <- ma.shrink[-(x$k+3),]
+ if (exponentiate) {
+ ma.dat <- exp(ma.dat)
+ ma.shrink <- exp(ma.shrink)
+ }
+ # generate data table (unless provided):
if (missing(labeltext)) {
decplaces <- function(x, signifdigits=3)
# number of decimal places (after decimal point)
@@ -1589,89 +1596,65 @@ forestplot.bayesmeta <- function(x, labeltext,
{
return(max(c(0, -(floor(log10(x))-(signifdigits-1)))))
}
+ # determine numbers of digits based on standard deviations:
if (exponentiate) {
- # determine numbers of digits based on standard deviations:
stdevs <- c(exp(x$y)*x$sigma,
- exp(x$summary["median","mu"])*x$summary["sd","mu"],
- exp(x$summary["median","theta"])*x$summary["sd","theta"])
+ exp(x$summary["median","mu"])*x$summary["sd","mu"])
+ if (prediction) stdevs <- c(stdevs, exp(x$summary["median","theta"])*x$summary["sd","theta"])
} else {
- stdevs <- c(x$sigma, x$summary["sd","mu"], x$summary["sd","theta"])
+ stdevs <- c(x$sigma, x$summary["sd","mu"])
+ if (prediction) stdevs <- c(stdevs, x$summary["sd","theta"])
}
stdevs <- abs(stdevs[is.finite(stdevs) & (stdevs != 0)])
formatstring <- paste0("%.", decplaces(stdevs, digits), "f")
- if (exponentiate) {
- labeltext <- cbind(x$labels,
- sprintf(formatstring, exp(x$y)),
- paste0("[", sprintf(formatstring, exp(x$y-q95*x$sigma)),
- ", ", sprintf(formatstring, exp(x$y+q95*x$sigma)), "]"))
- # add summary rows:
- labeltext <- rbind(labeltext,
- c("mean",
- sprintf(formatstring, exp(x$summary["median","mu"])),
- paste0("[",sprintf(formatstring, exp(x$summary["95% lower","mu"])),
- ", ", sprintf(formatstring, exp(x$summary["95% upper","mu"])), "]")),
- c("prediction",
- sprintf(formatstring, exp(x$summary["median","theta"])),
- paste0("[",sprintf(formatstring, exp(x$summary["95% lower","theta"])),
- ", ", sprintf(formatstring, exp(x$summary["95% upper","theta"])), "]")))
-
- } else {
- labeltext <- cbind(x$labels,
- sprintf(formatstring, x$y),
- paste0("[", sprintf(formatstring, x$y-q95*x$sigma),
- ", ", sprintf(formatstring, x$y+q95*x$sigma), "]"))
- # add summary rows:
- labeltext <- rbind(labeltext,
- c("mean",
- sprintf(formatstring, x$summary["median","mu"]),
- paste0("[",sprintf(formatstring, x$summary["95% lower","mu"]),
- ", ", sprintf(formatstring, x$summary["95% upper","mu"]), "]")),
- c("prediction",
- sprintf(formatstring, x$summary["median","theta"]),
- paste0("[",sprintf(formatstring, x$summary["95% lower","theta"]),
- ", ", sprintf(formatstring, x$summary["95% upper","theta"]), "]")))
+ # fill data table:
+ labeltext <- matrix(NA_character_, nrow=nrow(ma.dat), ncol=3)
+ labeltext[1,] <- c("study","estimate", "95% CI")
+ labeltext[,1] <- c("study", x$labels, "mean", "prediction")[1:nrow(ma.dat)]
+ for (i in 2:(nrow(ma.dat))) {
+ labeltext[i,2] <- sprintf(formatstring, ma.dat[i,"estimate"])
+ labeltext[i,3] <- paste0("[", sprintf(formatstring, ma.dat[i,"lower"]),
+ ", ", sprintf(formatstring, ma.dat[i,"upper"]), "]")
}
- # add header row:
- labeltext <- rbind(c("study","estimate", "95% CI"), labeltext)
}
horizl <- list(grid::gpar(col="grey"), grid::gpar(col="grey"))
names(horizl) <- as.character(c(2,x$k+2))
if (shrinkage) {
+ # default settings for estimate and shrinkage estimate plotting style:
if (missing(legend))
legend <- c("quoted estimate", "shrinkage estimate")
if (missing(col))
- col <- forestplot::fpColors(box=c("black", "grey45"),
- lines=c("black","grey45"),
- summary="grey30")
+ col <- fpColors(box=c("black", "grey45"),
+ lines=c("black","grey45"),
+ summary="grey30")
if (missing(fn.ci_norm))
- fn.ci_norm <- c(function(...) {forestplot::fpDrawPointCI(pch=15,...)},
- function(...) {forestplot::fpDrawPointCI(pch=18,...)})
+ fn.ci_norm <- c(function(...) {fpDrawPointCI(pch=15,...)},
+ function(...) {fpDrawPointCI(pch=18,...)})
mean.arg <- cbind(ma.dat[,1], ma.shrink[,1])
lower.arg <- cbind(ma.dat[,2], ma.shrink[,2])
upper.arg <- cbind(ma.dat[,3], ma.shrink[,3])
} else {
+ # default settings for estimate plotting style:
if (missing(col))
- col <- forestplot::fpColors(box="black", lines="black", summary="grey30")
+ col <- fpColors(box="black", lines="black", summary="grey30")
if (missing(fn.ci_norm))
- fn.ci_norm <- function(...) {forestplot::fpDrawPointCI(pch=15,...)}
+ fn.ci_norm <- function(...) {fpDrawPointCI(pch=15,...)}
mean.arg <- ma.dat[,1]
lower.arg <- ma.dat[,2]
upper.arg <- ma.dat[,3]
}
fp <- NULL
- if (plot) {
- fp <- forestplot::forestplot(labeltext = labeltext,
- mean = mean.arg,
- lower = lower.arg,
- upper = upper.arg,
- is.summary = c(TRUE,rep(FALSE,x$k),TRUE,TRUE),
- hrzl_lines = horizl,
- fn.ci_norm = fn.ci_norm,
- #fn.ci_sum = fpDrawBarCI,
- col = col,
- boxsize = boxsize,
- legend = legend, ...)
- }
+ if (plot) fp <- forestplot::forestplot(labeltext = labeltext,
+ mean = mean.arg,
+ lower = lower.arg,
+ upper = upper.arg,
+ is.summary = c(TRUE,rep(FALSE,x$k),TRUE,TRUE),
+ hrzl_lines = horizl,
+ fn.ci_norm = fn.ci_norm,
+ #fn.ci_sum = fpDrawBarCI,
+ col = col,
+ boxsize = boxsize,
+ legend = legend, ...)
invisible(list("data" = ma.dat[-1,],
"shrinkage" = ma.shrink[2:(x$k+1),],
"labeltext" = labeltext,
diff --git a/man/bayesmeta-package.Rd b/man/bayesmeta-package.Rd
index 4beaab3..a990e3d 100644
--- a/man/bayesmeta-package.Rd
+++ b/man/bayesmeta-package.Rd
@@ -16,7 +16,7 @@
Package: \tab bayesmeta\cr
Type: \tab Package\cr
Version: \tab 1.3\cr
- Date: \tab 2016-10-08\cr
+ Date: \tab 2016-10-13\cr
License: \tab GPL (>=2)
}
The main functionality is provided by the \code{\link{bayesmeta}()}
diff --git a/man/forestplot.bayesmeta.Rd b/man/forestplot.bayesmeta.Rd
index 2f01f13..39556ab 100644
--- a/man/forestplot.bayesmeta.Rd
+++ b/man/forestplot.bayesmeta.Rd
@@ -11,7 +11,7 @@
}
\usage{
\method{forestplot}{bayesmeta}(x, labeltext, exponentiate=FALSE,
- shrinkage=TRUE, digits=2, plot=TRUE,
+ prediction=TRUE, shrinkage=TRUE, digits=2, plot=TRUE,
fn.ci_norm, col, legend, boxsize=0.25, ...)
}
\arguments{
@@ -26,6 +26,10 @@
a logical flag indicating whether to exponentiate numbers (effect
sizes) in the table.
}
+ \item{prediction}{
+ a logical flag indicating whether to show the prediction interval
+ below the mean estimate estimate.
+ }
\item{shrinkage}{
a logical flag indicating whether to show shrinkage intervals along
with the quoted estimates.
@@ -48,7 +52,8 @@
resulting estimates (effect estimate and prediction interval,
as well as shrinkage estimates and intervals).
}
-\note{This function requires the \pkg{forestplot} package to be installed.
+\note{This function is based on the \pkg{forestplot} package's
+ \dQuote{\code{\link[forestplot]{forestplot}()}} function.
}
\author{
Christian Roever \email{christian.roever@med.uni-goettingen.de}
@@ -65,14 +70,16 @@
\seealso{
\code{\link{bayesmeta}},
\code{\link[forestplot]{forestplot}},
- \code{\link{forest.bayesmeta}}.
+ \code{\link{forest.bayesmeta}},
+ \code{\link{plot.bayesmeta}}.
}
\examples{
\dontrun{
# load data:
data("CrinsEtAl2014")
-# compute effect sizes (log-ORs):
+# compute effect sizes (log-ORs)
+# (requires "metafor" package to be installed):
require("metafor")
es.crins <- escalc(measure="OR",
ai=exp.AR.events, n1i=exp.total,
@@ -89,21 +96,24 @@ require("forestplot")
# default options:
forestplot.bayesmeta(crins.ma)
-# exponentiate values shown in table, show vertical line at OR=1:
-forestplot.bayesmeta(crins.ma, zero=1, expo=TRUE)
+# exponentiate values (shown in table and plot), show vertical line at OR=1:
+forestplot.bayesmeta(crins.ma, expo=TRUE, zero=1)
# logarithmic x-axis:
-forestplot.bayesmeta(crins.ma, xlog=TRUE)
+forestplot.bayesmeta(crins.ma, expo=TRUE, xlog=TRUE)
+
+# omit prediction interval:
+forestplot.bayesmeta(crins.ma, predict=FALSE)
# omit shrinkage intervals:
-forestplot.bayesmeta(crins.ma, zero=1, shrink=FALSE)
+forestplot.bayesmeta(crins.ma, shrink=FALSE)
# show more decimal places:
-forestplot.bayesmeta(crins.ma, zero=1, digits=3)
+forestplot.bayesmeta(crins.ma, digits=3)
# change table values:
# (here: add columns for event counts)
-fp <- forestplot.bayesmeta(crins.ma, plot=FALSE, expo=TRUE)
+fp <- forestplot.bayesmeta(crins.ma, expo=TRUE, plot=FALSE)
labtext <- fp$labeltext
labtext <- cbind(labtext[,1],
c("treatment",
@@ -116,7 +126,7 @@ labtext <- cbind(labtext[,1],
labtext[1,4] <- "OR"
print(fp$labeltext) # before
print(labtext) # after
-forestplot.bayesmeta(crins.ma, labeltext=labtext, xlog=TRUE)
+forestplot.bayesmeta(crins.ma, labeltext=labtext, expo=TRUE, xlog=TRUE)
# (see also the "forestplot" help for more arguments that you may change)
--
GitLab