Exploring the robustness of Bayes Factors: A convenient plotting function

August 22, 2013
By

(This article was first published on 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)

 

Bayes factor robustness analysis, two-sided (cf. Wagenmakers et al., 2011)

Bayes factor robustness analysis, two-sided (cf. Wagenmakers et al., 2011)

 

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")

 

Bayes factor robustness analysis, one-sided (cf. Wagenmakers et al., 2011; Wagenmakers and Morey, 2013)

Bayes factor robustness analysis, one-sided (cf. Wagenmakers et al., 2011; Wagenmakers and Morey, 2013)

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]

To leave a comment for the author, please follow the link and comment on his blog: Nicebread » R.

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



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.