Skip to content
Snippets Groups Projects
Commit 10241ada authored by Christian Roever's avatar Christian Roever
Browse files

fixed log-axis bug and added "prediction" argument for forestplot() function

parent e7d829c7
Branches forestplot
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
+ 1
− 1
View file @ 10241ada
......@@ -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",
......
R/bayesmeta.R
+ 54
− 71
View file @ 10241ada
......@@ -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,
......
man/bayesmeta-package.Rd
+ 1
− 1
View file @ 10241ada
......@@ -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}()}
......
man/forestplot.bayesmeta.Rd
+ 21
− 11
View file @ 10241ada
......@@ -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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment