### META-ANALYSIS FUNCTIONS
#' Extract data from `rma` object
#'
#' Extract data from fitted meta-analysis of type `rma`.
#'
#' Extract effect size, sampling variance/standard error from `rma` object. Then adds a couple of more useful columns with z scores, p values (two-tailed).
#' @param .rma (rma object) The rma analysis from metafor.
#' @export
#' @return A data frame.
#' @examples
#' library(metafor); data(european_ancestry)
#' meta = rma(european_ancestry$r, sei = european_ancestry$SE_r)
#' meta_extract_data(meta)
meta_extract_data = function(.rma) {
d = tibble(
es = .rma$yi,
var = .rma$vi,
se = .rma$vi %>% sqrt,
z = es/se,
p = pnorm(z, lower.tail = F)
)
#study name
d$name = sprintf("Study %d", 1:nrow(d))
d
}
#' Plot a contingency table with ggplot2
#'
#' Makes a pretty contingency table with ggplot2 using geom_tile.
#' @param .analysis (rma object) The rma analysis from metafor.
#' @param .names (chr vector) An optional vector of names to use.
#' @param .alphabetic_sort_names (lgl sclar) Alphabetically sort names? Default yes.
#' @export
#' @examples
#' library(metafor); data(european_ancestry)
#' meta = rma(european_ancestry$r, sei = european_ancestry$SE_r)
#' GG_forest(meta, .names = european_ancestry$Author_sample)
GG_forest = function(.analysis, .names = NULL, .alphabetic_sort_names = T) {
if (!inherits(.analysis, "rma")) stop("This function only works for rma objects from the metafor package.")
#extract effect sizes, SEs and names
d = meta_extract_data(.analysis)[c("es", "se", "name")]
d$meta = "study"
#names if given
if (!is.null(.names)) {
d$name = .names
}
#make names unique if necessary
if (any(duplicated(d$name))) {
d$name %<>% str_uniquify
}
#sort?
if (.alphabetic_sort_names) {
d$name %<>% factor %>% forcats::fct_rev()
}
#extract main effect
d_meta = tibble::tibble(es = .analysis$b %>% as.vector,
se = .analysis$se,
name = "Main effect",
meta = "meta"
)
#horizontal space case
d_hline = tibble::tibble(es = .analysis$b %>% as.vector,
se = .analysis$se,
name = "",
meta = "invis"
)
#add main effect to d
d = rbind(d, d_meta, d_hline)
#make sure meta effect is in the bottom
d$name %<>% factor() %>% forcats::fct_relevel(c("Main effect", ""))
#plot
ggplot2::ggplot(d, aes(es, name, color = meta)) +
geom_point() +
geom_errorbarh(aes(xmin = es - se * 1.96,
xmax = es + se * 1.96)) +
geom_hline(yintercept = 2, linetype = "dashed") +
theme_bw() +
scale_y_discrete(name = NULL) +
scale_colour_manual(values = c("white", "black", "black"), guide = "none") +
xlab("Effect size")
}
#' Funnel plot with ggplot2
#'
#' Makes a pretty funnel plot using [ggplot2].
#' @param .analysis (rma object) The rma analysis from [metafor].
#' @param .CI (chr vector) Confidence interval to use.
#' @param .study_CI (lgl vector) Whether to plot confidence intervals for individual studies.
#' @export
#' @examples
#' library(metafor); data(european_ancestry)
#' meta = rma(european_ancestry$r, sei = european_ancestry$SE_r)
#' GG_funnel(meta)
#' GG_funnel(meta, .study_CI = T)
GG_funnel = function(.analysis, .CI = .95, .study_CI = F) {
if (!inherits(.analysis, "rma")) stop("This function only works for rma objects from the metafor package.")
#convert CI to se z
se_z = qnorm(1 - (1-.CI)/2)
#extract main effect
d_meta = tibble::tibble(es = .analysis$b %>% as.vector,
se = .analysis$se
)
#extract effect sizes and SEs
d = tibble::tibble(es = .analysis$yi,
se = sqrt(.analysis$vi),
upper = d_meta$es + se_z * se,
lower = d_meta$es - se_z * se,
outlier = !is_between(es, lower, upper)
)
#calculate funnel
d_funnel = tibble::tibble(se = seq(0, max(d$se)*1.1, length.out = 1000),
upper = d_meta$es + se * se_z,
lower = d_meta$es - se * se_z)
d_polygon = tibble::tibble(x = c(min(d_funnel$lower), d_meta$es, max(d_funnel$upper)),
y = c(max(d_funnel$se), 0, max(d_funnel$se)))
#plot
gg = ggplot2::ggplot() +
geom_line(data = d_funnel, aes(upper, se)) +
geom_line(data = d_funnel, aes(lower, se)) +
geom_polygon(data = d_polygon, aes(x, y), fill = "grey") +
geom_vline(linetype = "dashed", xintercept = d_meta$es) +
geom_point(data = d, aes(es, se, color = outlier)) +
scale_color_manual(guide = "none", values = c("black", "red")) +
scale_y_reverse() +
theme_bw() +
xlab("Effect size")
#study CIs
if (.study_CI) {
gg = gg +
geom_errorbarh(data = d, aes(xmin = es - se_z * se,
xmax = es + se_z * se,
y = se))
}
gg
}
# TIVA --------------------------------------------------------------------
#' Funnel plot with ggplot2
#'
#' Makes a pretty funnel plot using [ggplot2].
#' @param .analysis (rma object) The rma analysis from [metafor].
#' @param .CI (chr vector) Confidence interval to use.
#' @param .study_CI (lgl vector) Whether to plot confidence intervals for individual studies.
#' @export
#' @examples
#' library(metafor)
#' data(european_ancestry)
#' meta = rma(european_ancestry$r, sei = european_ancestry$SE_r)
#' meta_TIVA(meta)
meta_TIVA = function(.rma, print_plot = T) {
#extract data
d = meta_extract_data(.rma)
#test
ob_var = var(d$z)
chi_stat = var(d$z) * (nrow(d) - 1)
test_p = pchisq(chi_stat, nrow(d) - 1, lower.tail = T)
#make text
text_ = sprintf("Chi^2 test for variance <1 (TIVA)\ndf=%d\nObserved variance of z = %f\np = %f\n", nrow(d), ob_var, test_p)
#plot distribution of z's
ggplot_ = GG_denhist(d, var = "z", vline = NULL) +
GG_text(text_)
#print plot?
if (print_plot) print(ggplot_)
#return
list(observed_var = ob_var,
chi_stat = chi_stat,
p = test_p,
plot = ggplot_)
}
# meta_pcurve -------------------------------------------------------------
#' P-curve plot with ggplot2
#'
#' Makes a pretty p-curve plot using [ggplot2].
#' @param .analysis (rma object) The rma analysis from [metafor].
#' @param p_cutoff (num) P-value cutoff to use.
#' @param print_plot (lgl) Print the plot?
#' @param binom_method (chr) Which method to use to calculate binomial p-value.
#' @export
#' @examples
#' library(metafor)
#' data(european_ancestry)
#' meta = rma(european_ancestry$r, sei = european_ancestry$SE_r)
#' meta_pcurve(meta)
meta_pcurve = function(.analysis, p_cutoff = .05, print_plot = T, binom_method = "exact", text_pos = "tl") {
#get data
d = meta_extract_data(.analysis)
#get p values in less than cutoff
d %<>%
dplyr::filter(p <= p_cutoff)
#count by interval
p_bins = tibble(
bin = c(".01", ".02", ".03", ".04", ".05"),
count = c(sum(is_between(d$p, a = 0, b = .015, include_upper = F)),
sum(is_between(d$p, a = .015, b = .025, include_upper = F)),
sum(is_between(d$p, a = .025, b = .035, include_upper = F)),
sum(is_between(d$p, a = .035, b = .045, include_upper = F)),
sum(is_between(d$p, a = .045, b = 1, include_upper = T))
),
prop = count/nrow(d)
)
#CIs
p_bins$CI_lower = map_dbl(p_bins$count, ~binom::binom.confint(., n = nrow(d), methods = binom_method) %>% .[[1, "lower"]])
p_bins$CI_upper = map_dbl(p_bins$count, ~binom::binom.confint(., n = nrow(d), methods = binom_method) %>% .[[1, "upper"]])
#binomial test
binom_025 = binom.test(sum(d$p < .025), nrow(d))
#plot text
plot_text = sprintf("Binomial test for evidential value, p = %.4f\nk = %d", binom_025$p.value, nrow(d))
#plot
plot_ = ggplot(p_bins, aes(bin, prop)) +
geom_point() +
geom_line(aes(group = 1)) +
scale_y_continuous("Percentage of tests", labels = scales::percent) +
scale_x_discrete("p-value bin") +
GG_text(plot_text, text_pos = text_pos) +
theme_bw()
#print
if (print_plot) print(plot_)
#output
list(
p_binomial = binom_025,
plot = plot_
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.