When Evidence and Significance Collide
bf10_ <- function(x, se, mu, tau, alternative = "two.sided") {
bfUncorrect <- stats::dnorm(x = x, mean = mu, sd = sqrt(se^2 + tau^2))/
stats::dnorm(x = x, mean = 0, sd = se)
if (alternative == "two.sided") {
bf <- bfUncorrect
} else {
postVar <- 1/(1/se^2 + 1/tau^2)
postMean <- (x/se^2 + mu/tau^2)*postVar
if (alternative == "greater") {
bf <- bfUncorrect*pnorm(q = postMean/sqrt(postVar))/pnorm(q = mu/tau)
} else if (alternative == "less") {
bf <- bfUncorrect*pnorm(q = -postMean/sqrt(postVar))/pnorm(q = -mu/tau)
} else {
bf <- NaN
}
}
return(bf)
}
bf10 <- Vectorize(FUN = bf10_)
## make a nice contour plot
## -----------------------------------------------------------------------------
## compute contours
formatBF <- function(BF) ifelse(BF < 1, paste0("1/", round(1/BF, 2)), as.character(round(BF, 2)))
bfBreaks <- c(1/100, 1/30, 1/10, 1/3, 1, 3, 10, 30, 100)
bfLabs <- formatBF(bfBreaks)
## with ggplot
library(ggplot2)
library(colorspace)
## Significant example that lacks evidence for H1
## https://pubmed.ncbi.nlm.nih.gov/34469764/
logHR <- log(0.96)
pval <- 0.045
selogHR <- abs(logHR)/qnorm(p = 1 - pval/2)
# BF max:
1/exp(-(logHR/selogHR)^2/2)
# BF default
bf10(x = logHR, se = selogHR, mu = 0, tau = 2, alternative = "less")
## compute BF
priorM <- seq(log(0.5), log(1), length.out = 200)
## start not at zero because there numerical problems with one.sided bf
priorSD <- seq(2*.Machine$double.eps, 2, length.out = 200)
priorGrid <- expand.grid(mean = priorM, sd = priorSD)
alt <- "less"
priorGrid$bf <- bf10(x = logHR, se = selogHR, mu = priorGrid$mean,
tau = priorGrid$sd, alternative = alt)
priorGrid$bf <- ifelse(priorGrid$bf < 1/110, 1/110, priorGrid$bf)
ORbreaks <- seq(0.5, 1, 0.10)
sdBreaks <- seq(0, 2, 0.50)
contours <- c(bfBreaks)
p1 <- ggplot(data = priorGrid, aes(x = sd, y = mean)) +
geom_raster(aes(fill = bf), alpha = 0.8) +
geom_contour(aes(z = bf), breaks = contours, color = "black", lty = 2, alpha = 0.3) +
metR::geom_text_contour(data = priorGrid,
aes(x = sd, y = mean, z = bf), ## TODO add parsed label
breaks = round(contours, 2), size = 3,
skip = 0, parse = FALSE,
label.placer = metR::label_placer_random(seed = 53)) +
scale_fill_continuous_diverging(palette = "Blue-Red", trans = "log", mid = 0,
rev = FALSE, p1 = 0.5, p2 = 1.2, ## <- tweak these values
breaks = bfBreaks, limits = c(1/110, 110),
labels = bfLabs, name = bquote(BF[10] ~ " ")) +
guides(fill = guide_colorbar(barheight = 0.4, barwidth = 28)) +
scale_y_continuous(name = "Prior mean (HR)", expand = c(0, 0), breaks = log(ORbreaks),
labels = ORbreaks) +
scale_x_continuous(name = bquote("Prior standard deviation (on logHR scale)"),
expand = c(0, 0), breaks = sdBreaks) +
coord_cartesian(ylim = c(min(priorM), max(priorM))) +
geom_line(mapping = aes(x = x, y = y), data = data.frame(x = c(0, 2), y = c(log(0.70), log(0.70))), linetype = 3) +
theme_bw() +
theme(legend.position = "top",
panel.grid = ggplot2::element_blank(),
legend.text.align = 0)
## Non-significant example with evidence for H1
## https://pubmed.ncbi.nlm.nih.gov/34932079/
logHR <- log(1.63)
selogHR <- diff(log(c(0.93, 2.85)))/(2*qnorm(p = 0.975))
# BF max:
1/exp(-(logHR/selogHR)^2/2)
# BF default
bf10(x = logHR, se = selogHR, mu = 0, tau = 2, alternative = "greater")
## compute BF
priorM <- seq(log(1), log(3), length.out = 200)
## do not start at zero because there numerical problems with one.sided bf
priorSD <- seq(2*.Machine$double.eps, 2, length.out = 200)
priorGrid <- expand.grid(mean = priorM, sd = priorSD)
alt <- "greater"
priorGrid$bf <- bf10(x = logHR, se = selogHR, mu = priorGrid$mean,
tau = priorGrid$sd, alternative = alt)
sdBreaks <- seq(0, 2, 0.50)
contours <- c(bfBreaks)
ORbreaks <- c(1, 1.5, 2, 2.5, 3)
p2 <- ggplot(data = priorGrid, aes(x = sd, y = mean)) +
geom_raster(aes(fill = bf), alpha = 0.8) +
geom_contour(aes(z = bf), breaks = contours, color = "black", lty = 2, alpha = 0.3) +
metR::geom_text_contour(data = priorGrid,
aes(x = sd, y = mean, z = bf), ## TODO add parsed label
breaks = round(contours, 2), size = 3,
skip = 0, parse = FALSE,
label.placer = metR::label_placer_random(seed = 46)) +
scale_fill_continuous_diverging(palette = "Blue-Red", trans = "log", mid = 0,
rev = FALSE, p1 = 1, p2 = 1.75, ## <- tweak these values
breaks = bfBreaks, limits = c(1/110, 110),
labels = bfLabs, name = bquote(BF[10] ~ " ")) +
guides(fill = guide_colorbar(barheight = 0.4, barwidth = 28)) +
scale_y_continuous(name = "Prior mean (OR)", expand = c(0, 0), breaks = log(ORbreaks),
labels = ORbreaks) +
scale_x_continuous(name = bquote("Prior standard deviation (on logOR scale)"),
expand = c(0, 0), breaks = sdBreaks) +
coord_cartesian(ylim = c(min(priorM), max(priorM))) +
geom_line(mapping = aes(x = x, y = y), data = data.frame(x = c(0, 2), y = c(log(2), log(2))), linetype = 3) +
theme_bw() +
theme(legend.position = "top",
panel.grid = ggplot2::element_blank(),
legend.text.align = 0)
library(ggpubr)
(pCombined <- ggarrange(p1, p2, ncol = 2, common.legend = TRUE, align = "v",
labels = c("(A)", "(B)"), font.label = list(size = 12)))
ggsave(file = "figure.pdf", plot = pCombined, height = 5, width = 10)
留言
張貼留言