/
forest_plot.R
132 lines (121 loc) · 5.21 KB
/
forest_plot.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#' Forest Plot
#'
#' Creates a forest plot for \code{\link{gMAP}} analysis objects.
#'
#' @param x \code{\link{gMAP}} object.
#' @param prob confidence interval width and probability mass of credible intervals.
#' @param est can be set to one of \code{both} (default), \code{MAP}, \code{Mean} or \code{none}. Controls which model estimates are to be included.
#' @param model controls which estimates are displayed per study. Either \code{stratified} (default), \code{both} or \code{meta}.
#' @param point_est shown point estimate. Either \code{median} (default) or \code{mean}.
#' @param size controls point and linesize.
#' @param alpha transparency of reference line. Setting \code{alpha=0}
#' suppresses the reference line.
#'
#' @details The function creates a forest plot suitable for
#' \code{\link{gMAP}} analyses. Note that the Meta-Analytic-Predictive
#' prior is included by default in the plot as opposed to only showing
#' the estimated model mean. See the examples below to obtain standard
#' forest plots.
#'
#' Also note that the plot internally flips the x and
#' y-axis. Therefore, if you want to manipulate the x-axis, you have
#' to give commands affecting the y-axis (see examples).
#'
#' @template plot-help
#'
#' @return The function returns a \pkg{ggplot2} plot object.
#'
#' @seealso \code{\link{gMAP}}
#'
#' @examples
#' # we consider the example AS MAP analysis
#' example(AS)
#'
#' # default forest plot for a gMAP analysis
#' forest_plot(map_AS)
#'
#' # standard forest plot (only stratified estimate and Mean)
#' forest_plot(map_AS, est=c("Mean"), model="stratified")
#'
#' # to further customize these plots, first load bayesplot and ggplot2
#' library(bayesplot)
#' library(ggplot2)
#'
#' # to make plots with red colors, big fonts for presentations, suppress
#' # the x axis label and add another title (with a subtitle)
#' color_scheme_set("red")
#' theme_set(theme_default(base_size=16))
#' forest_plot(map_AS, size=2) +
#' yaxis_title(FALSE) +
#' ggtitle("Ankylosing Spondylitis Forest Plot",
#' subtitle="Control Group Response Rate")
#'
#' # the defaults are set with
#' color_scheme_set("blue")
#' theme_set(theme_default(base_size=12))
#'
#' @export
forest_plot <- function(x,
prob=0.95,
est = c("both", "MAP", "Mean", "none"),
model = c("stratified", "both", "meta"),
point_est = c("median", "mean"),
size=1.25,
alpha=0.5) {
assert_number(prob, lower=0, upper=1)
assert_that(inherits(x, "gMAP"))
assert_that(x$has_intercept)
est <- match.arg(est)
low <- (1-prob)/2
up <- 1-low
strat <- as.data.frame(x$est_strat(1-prob))
strat <- cbind(strat[1:2], median=strat$mean, strat[3:4])
names(strat)[3:4] <- c("low", "up")
fit <- as.data.frame(fitted(x, type="response", probs=c(0.5, low, up)))
est <- match.arg(est)
model <- match.arg(model)
point_est <- match.arg(point_est)
if(est == "both") est <- c("MAP", "Mean")
if(model == "both") model <- c("stratified", "meta")
pred_est <- as.data.frame(do.call(rbind, summary(x, probs=c(0.5, low, up), type="response")[c("theta.pred", "theta")]))
pred_est <- transform(pred_est, study=c("MAP", "Mean") , model="meta")
pred_est <- pred_est[c("MAP", "Mean") %in% est,]
names(pred_est)[1:5] <- names(strat) <- names(fit) <- c("mean", "sem", "median", "low", "up")
comb <- rbind(if("stratified" %in% model) transform(strat, study=rownames(strat), model="stratified"),
if("meta" %in% model) transform(fit, study=rownames(strat), model="meta"),
pred_est
)
comb <- within(comb, {
model <- factor(model, levels=c("meta", "stratified"))
study <- factor(study, levels=rev(c(rownames(strat), "Mean", "MAP")))
})
opts <- list(position=position_dodge(width=0.3), size=size)
xlab_str <- switch(x$family$family,
gaussian="Response",
binomial="Response Rate",
poisson="Counting Rate")
graph <- ggplot(comb, aes(x=.data$study, y=.data[[point_est]], ymin=.data$low, ymax=.data$up, linetype=.data$model, color=.data$model))
if(any(c("MAP", "Mean") %in% est)) {
ref_line <- est[est %in% c("Mean", "MAP")][1]
ref_data <- subset(pred_est, study == ref_line)
no_ref <- sum(est %in% c("Mean", "MAP"))
graph <- graph + geom_rect(ymin=-Inf, ymax=Inf, xmin=0, xmax=no_ref + 0.5,
fill=get_color("l"),
color=get_color("l"), show.legend=FALSE) +
geom_hline(yintercept=ref_data[1,point_est],
color=get_color("mh"),
alpha=alpha,
size=size)
}
graph <- graph +
scale_color_manual("Model", values=get_color(c("mh", "m"))) +
do.call(geom_pointrange, opts) +
ylab(xlab_str) +
scale_linetype_discrete("Model") +
theme(axis.line.y=element_blank(), axis.ticks.y=element_blank()) +
bayesplot::bayesplot_theme_get() +
bayesplot::xaxis_title(FALSE) +
coord_flip() +
bayesplot::legend_none()
graph
}