**Nicebread » R**, and kindly contributed to R-bloggers)

One critique frequently heard about Bayesian statistics is the subjectivity of the assumed prior distribution. If one is cherry-picking a prior, of course the posterior can be tweaked, especially when only few data points are at hand. For example, see the Scholarpedia article on Bayesian statistics:

In the uncommon situation that the data are extensive and of simple structure, the prior assumptions will be unimportant and the assumed sampling model will be uncontroversial. More generally we would like to report that any conclusions are robust to reasonable changes in both prior and assumed model: this has been termed inference robustness

(David Spiegelhalter and Kenneth Rice (2009) Bayesian statistics. Scholarpedia, 4(8):5230.)

Therefore, it is suggested that …

In particular, audiences should ideally fully understand the contribution of the prior distribution to the conclusions. (ibid)

In the example of Bayes factors for t tests (Rouder, Speckman, Sun, Morey, & Iverson, 2009), the assumption that has to be defined a priori is the effect size δ expected under the H1. In the BayesFactor package for R, this can be adjusted via the r parameter. By default, it is set to 0.5, but it can be made wider (larger r’s, which means one expects larger effects) or narrower (r’s close to zero, which means one expects smaller effects in the population).

In their reanalysis of Bem’s ESP data, Wagenmakers, Wetzels, Borsboom, Kievit, and van der Maas (2011, PDF) proposed a robustness analysis for Bayes factors (BF), which simply shows the BF for a range of priors. If the conclusion is the same for a large range of priors, it could be judged to be robust.

I wrote an R function that can generate plots like this. Here’s a reproduction of Wagenmakers’ et al (2011) analysis of Bem’s data – it looks pretty identical:

^{?}View Code RSPLUS

## Bem data, two sided # provide t values, sample sizes, and the location(s) of the red dot(s) BFrobustplot( ts=c(2.51, 2.55, 2.23, 1.74, 1.92, 2.39, 2.03, 1.8, 1.31, 2.96), ns=c(100, 97, 100, 150, 100, 150, 99, 150, 200, 50), dots=1) |

You can throw in as many t values and corresponding sample sizes as you want. Furthermore, the function can compute one-sided Bayes factors as described in Wagenmakers and Morey (2013). If this approach is applied to the Bem data, the plot looks as following – everything is shifted a bit into the H1 direction:

^{?}View Code RSPLUS

# Bem data one-sided BFrobustplot( ts=c(2.51, 2.55, 2.23, 1.74, 1.92, 2.39, 2.03, 1.8, 1.31, 2.96), ns=c(100, 97, 100, 150, 100, 150, 99, 150, 200, 50), dots=1, sides="one") |

Finally, here’s the function:

^{?}View Code RSPLUS

## This source code is licensed under the FreeBSD license ## (c) 2013 Felix Schönbrodt #' @title Plots a comparison of a sequence of priors for t test Bayes factors #' #' @details #' #' #' @param ts A vector of t values #' @param ns A vector of corresponding sample sizes #' @param rs The sequence of rs that should be tested #' @param labels Names for the studies (displayed in the facet headings) #' @param dots Values of r's which should be marked with a red dot #' @param plot If TRUE, a ggplot is returned. If false, a data frame with the computed Bayes factors is returned #' @param sides If set to "two" (default), a two-sided Bayes factor is computed. If set to "one", a one-sided Bayes factor is computed. In this case, it is assumed that positive t values correspond to results in the predicted direction and negative t values to results in the unpredicted direction. For details, see Wagenmakers, E. J., & Morey, R. D. (2013). Simple relation between one-sided and two-sided Bayesian point-null hypothesis tests. #' @param nrow Number of rows of the faceted plot. #' #' @references #' #' Rouder, J. N., Speckman, P. L., Sun, D., Morey, R. D., & Iverson, G. (2009). Bayesian t-tests for accepting and rejecting the null hypothesis. Psychonomic Bulletin and Review, 16, 225-237. #' Wagenmakers, E.-J., & Morey, R. D. (2013). Simple relation between one-sided and two-sided Bayesian point-null hypothesis tests. Manuscript submitted for publication #' Wagenmakers, E.-J., Wetzels, R., Borsboom, D., Kievit, R. & van der Maas, H. L. J. (2011). Yes, psychologists must change the way they analyze their data: Clarifications for Bem, Utts, & Johnson (2011) BFrobustplot <- function(ts, ns, rs=seq(0, 2, length.out=200), labels=c(), dots=1, plot=TRUE, sides="two", nrow=2, xticks=3) { library(BayesFactor) # compute one-sided p-values from ts and ns ps <- pt(ts, df=ns-1, lower.tail = FALSE) # one-sided test # add the dots location to the sequences of r's rs <- c(rs, dots) res <- data.frame() for (r in rs) { # first: calculate two-sided BF B_e0 <- c() for (i in 1:length(ts)) B_e0 <- c(B_e0, exp(ttest.tstat(t = ts[i], n1 = ns[i], rscale=r)$bf)) # second: calculate one-sided BF B_r0 <- c() for (i in 1:length(ts)) { if (ts[i] > 0) { # correct direction B_r0 <- c(B_r0, (2 - 2*ps[i])*B_e0[i]) } else { # wrong direction B_r0 <- c(B_r0, (1 - ps[i])*2*B_e0[i]) } } res0 <- data.frame(t=ts, n=ns, BF_two=B_e0, BF_one=B_r0, r=r) if (length(labels) > 0) { res0$labels <- labels res0$heading <- factor(1:length(labels), labels=paste0(labels, "\n(t = ", ts, ", df = ", ns-1, ")"), ordered=TRUE) } else { res0$heading <- factor(1:length(ts), labels=paste0("t = ", ts, ", df = ", ns-1), ordered=TRUE) } res <- rbind(res, res0) } # define the measure to be plotted: one- or two-sided? # log-transform and reverse the sign in order to get the direction of BF as plotted in Wagenmakers, E., Wetzels, R., & Borsboom, D. (n.d.). Yes, Psychologists Must Change the Way They Analyze Their Data: Clarifications for Bem, Utts, and Johnson (2011). Journal of Personality. res$BF <- 1/res[, paste0("BF_", sides)] if (plot==TRUE) { library(ggplot2) p1 <- ggplot(res, aes(x=r, y=log(BF))) + geom_line() + facet_wrap(~heading, nrow=nrow) + theme_bw() + ylab("log(BF)") p1 <- p1 + geom_hline(yintercept=c(c(-log(c(30, 10, 3)), log(c(3, 10, 30)))), linetype="dotted", color="darkgrey") p1 <- p1 + geom_hline(yintercept=log(1), linetype="dashed", color="darkgreen") # add the dots p1 <- p1 + geom_point(data=res[res$r %in% dots,], aes(x=r, y=log(BF)), color="red", size=2) # add annotation p1 <- p1 + annotate("text", x=max(rs)*1.8, y=-2.85, label="Strong~H[1]", hjust=1, vjust=.5, size=3, color="black", parse=TRUE) p1 <- p1 + annotate("text", x=max(rs)*1.8, y=-1.7 , label="Moderate~H[1]", hjust=1, vjust=.5, size=3, color="black", parse=TRUE) p1 <- p1 + annotate("text", x=max(rs)*1.8, y=-.55 , label="Anectodal~H[1]", hjust=1, vjust=.5, size=3, color="black", parse=TRUE) p1 <- p1 + annotate("text", x=max(rs)*1.8, y=2.86 , label="Strong~H[0]", hjust=1, vjust=.5, size=3, color="black", parse=TRUE) p1 <- p1 + annotate("text", x=max(rs)*1.8, y=1.7 , label="Moderate~H[0]", hjust=1, vjust=.5, size=3, color="black", parse=TRUE) p1 <- p1 + annotate("text", x=max(rs)*1.8, y=.55 , label="Anectodal~H[0]", hjust=1, vjust=.5, vjust=.5, size=3, color="black", parse=TRUE) # set scale ticks p1 <- p1 + scale_y_continuous(breaks=c(c(-log(c(30, 10, 3)), 0, log(c(3, 10, 30)))), labels=c("-log(30)", "-log(10)", "-log(3)", "log(1)", "log(3)", "log(10)", "log(30)")) p1 <- p1 + scale_x_continuous(breaks=seq(min(rs), max(rs), length.out=xticks)) return(p1) } else { return(res) } } |

*References*

Rouder, J. N., Speckman, P. L., Sun, D., Morey, R. D., & Iverson, G. (2009). Bayesian t-tests for accepting and rejecting the null hypothesis. *Psychonomic Bulletin and Review, 16*, 225-237. [for a PDF, see bottom of this page]

Wagenmakers, E.-J., & Morey, R. D. (2013). Simple relation between one-sided and two-sided Bayesian point-null hypothesis tests. Manuscript submitted for publication (website)

Wagenmakers, E.-J., Wetzels, R., Borsboom, D., Kievit, R. & van der Maas, H. L. J. (2011). Yes, psychologists must change the way they analyze their data: Clarifications for Bem, Utts, & Johnson (2011) [PDF]

**leave a comment**for the author, please follow the link and comment on their blog:

**Nicebread » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...