inst/doc/BayesCombo_vignette.R

In stanlazic/BayesCombo: Bayesian Evidence Combination

## ---- echo = FALSE, include=FALSE----------------------------------------
library(knitr)
opts_chunk$set(fig.align='center') 

## ----setup, echo = FALSE, warning= FALSE, message=FALSE------------------
library(BayesCombo)

## ---- echo=FALSE, out.width='90%', fig.cap='Five steps to get from a prior to a posterior probability of a hypothesis. The effect size is the parameter $\\theta$, and `H<`, `H0`, and `H>` represent the three hypotheses of negative, zero, and positive effects.'----
knitr::include_graphics("procedure.pdf")

## ---- fig.height=4, fig.width=4, fig.cap='Effect of fluoxetine (Prozac) on rats in the Forced Swim Test. Data are from Lazic [4]'----
library(labstats)
par(las=1)
plot(time.immob ~ dose, data=fluoxetine, col="royalblue",
	ylab="Time immobile (s)", ylim=c(0, 250), xlim=c(0, 250),
	xlab="Dose of fluoxetine (mg/L in drinking water)")
abline(lm(time.immob ~ dose, data=fluoxetine), col="royalblue") # add reg. line

## ------------------------------------------------------------------------
summary(lm(time.immob ~ dose, data=fluoxetine))$coef

## ---- echo=FALSE, fig.height=6.5, fig.width=3, fig.cap='Prior, Likelihood (data), and posterior distributions for an experiment. Bayes factors for the three hypotheses are calculated as the ratio of posterior to prior areas or heights.'----
# x-values for plotting
xx <- seq(-5,5, length.out = 401)

# observed data (beta=0.75, se=1)... aka likelihood
obs <- dnorm(xx, 0.75, 1) 

# prior distribution (se is for a 99% CI)
up <- dnorm(xx, 0, 1.29)


# posterior (standard Bayesian updating for normal conjugate prior)
post.b <- ((0.75/1^2) + (0/1.29^2)) / ((1/1^2) + (1/1.29^2))
post.se <- sqrt(1/((1/1^2) + (1/1.29^2)))

# posterior distribution
post <- dnorm(xx, post.b, post.se)



par(mfrow=c(3,1),
    mar=c(4,3,1,2),
    las=0)


# prior
plot(up ~ xx, type="l", ylim=c(0, 0.5), xlab="", yaxt="n")
mtext("Prior", side=4, line=1)

polygon(c(xx[1:201],0), c(up[1:201],0),
        border="darkgrey", col="lightgrey", lwd=1.5)

polygon(c(xx[201:400],0), c(up[201:400],0),
        border="royalblue", col="#348ABD", lwd=1.5)

points(up[201] ~ xx[201], pch=16)

text(x=-1, y=0.1, labels = "a", cex=1.25)
text(x=0, y=0.35, labels = "b", cex=1.25)
text(x=1, y=0.1, labels = "c", cex=1.25)


# likelihood
plot(obs ~ xx, type="l", ylim=c(0, 0.5), xlab="", yaxt="n")
mtext("Likelihood", side=4, line=1)


# posterior
plot(post  ~ xx, type="l", ylim=c(0, 0.5), xlab="", yaxt="n")
mtext("Posterior", side=4, line=1)

polygon(c(xx[1:201],0), c(post[1:201],0),
        border="darkgrey", col="lightgrey")

polygon(c(xx[201:400],0), c(post[201:400],0),
        border="royalblue", col="#348ABD")

points(post[201] ~ xx[201], pch=16)

text(x=-1+0.47, y=0.1, labels = "d", cex=1.25)
text(x=0, y=0.48, labels = "e", cex=1.25)
text(x=1, y=0.1, labels = "f", cex=1.25)

mtext("Effect size", side=1, line=2.5)

## ------------------------------------------------------------------------
x <- pph(beta = -0.252, se.beta = 0.099)
summary(x)

## ---- fig.height=4.5, fig.width=4.5, fig.cap='Prior, likelihood (data), and posterior distributions for the fluoxetine data.'----
par(las=1)
plot(x, leg.loc = "topright", xlab="Effect size (slope of regression line)")

## ---- fig.height=4.5, fig.width=4.5, fig.cap='Prior, likelihood (data), and posterior distributions for the fluoxetine data. Same as Fig. 4 but with a wider prior.'----
x2 <- pph(beta = -0.252, se.beta = 0.099, se.mult = 2)
summary(x2)

par(las=1)
plot(x2, leg.loc = "topright", xlab="Effect size (slope of regression line)")

## ---- fig.height=4.5, fig.width=4.5, fig.cap='Prior, likelihood (data), and posterior distributions for the fluoxetine data. Same as Fig. 4 but with alternative options specified.'----
x3 <- pph(beta = -0.252, se.beta = 0.099, beta0=0, se0=1.2,
          H0 = c(-0.05, 0.05), H.priors=c(0.495, 0.495, 0.01))
summary(x3)

par(las=1)
plot(x3, leg.loc = "topright", xlab="Effect size (slope of regression line)")

## ------------------------------------------------------------------------
x4 <- ev.combo(beta = c(2.3, 1.2, 0.2, 0.44),
	           se.beta = c(1.03, 0.75, 0.16, 0.28))

## ----fig.height=4, fig.width=4, fig.cap='Forest plot showing effect sizes for four hypothetical clinical studies. Positive effect sizes indicate an improvement and grey lines are the priors.'----
par(las=1)
forestplot(x4)
abline(v=0, lty=2)

## ----summaryPPH----------------------------------------------------------
summary(x4)

## ----out.width='50%', fig.height=5, fig.width=5, fig.cap='Accumulation of support for the `H>` hypothesis as studies are added. PPH = Posterior Probability of a Hypothesis.'----
par(las=1)
plot(x4, ylab="PPH", xlab="Study")

## ----out.width='90%', fig.height=5, fig.width=8, fig.cap='Posterior Probabilities for each hypothesis for four studies. Each graph shows the same information and only groups the hypothesis and studies differently.'----
par(mfrow=c(1,2))
dotchart(x4$pph.uniform, xlim=c(0,1), xlab="PPH", pch=21, bg="grey")
dotchart(t(x4$pph.uniform), xlim=c(0,1), xlab="PPH", pch=21, bg="grey")