An interactive Bayesian test of multiple hypotheses using the shiny package in R
This is a shiny app for performing a Bayesian test of multiple hypotheses on a proportion p:
- H0: p = .5
- H1: p < .5
- H2: p > .5
To run the app please go through the following steps.
- Install R.
- Install Rstudio.
- Via Tools -> install packages -> shiny.
- In the console, write: “library(shiny)”.
- Create folder named “Bayesian proportion test”.
- Create a R-script file in Rstudio and call the file: ui.R.
- Save ui.R in the folder “Bayesian proportion test”.
- Create a R-script file in Rstudio called: server.R.
- Save server.R in the folder “Bayesian proportion test”.
Paste the syntax below to ui.R. (scroll down for the next step).
library(shiny)
shinyUI(fluidPage(
# Application title
titlePanel("Bayesian multiple hypothesis test of a proportion: H0:p=.5 vs H1:p.5"),
# Sidebar with a slider input
sidebarLayout(
sidebarPanel(
sliderInput("alpha0",
"prior alpha",
min = 0.001,
max = 20,
value = 1),
sliderInput("beta0",
"prior beta",
min = .001,
max = 20,
value = 1),
sliderInput("nobs",
"Number of experiments:",
min = 1,
max = 100,
value = 60),
sliderInput("ncorrect",
"Number of correct answers:",
min = 1,
max = 100,
value = 38),
checkboxInput("H0", "H0 versus Hu",F),
checkboxInput("H1", "H1 versus Hu",F),
checkboxInput("H2", "H2 versus Hu",F)
),
# Show a plot of the generated distribution
mainPanel(
plotOutput("distPlot"),
textOutput("BF0u"),
textOutput("BF1u"),
textOutput("BF2u"),
textOutput("PrHP0"),
textOutput("PrHP1"),
textOutput("PrHP2"),
textOutput("PHP0"),
textOutput("PHP1"),
textOutput("PHP2")
)
)
))
- Paste the syntax below to server.R. (scroll down for the final step).
library(shiny)
shinyServer(
function(input, output) {
output$distPlot = renderPlot({
seq1 = seq(0,1,length=1e3)
alpha0 = input$alpha0
beta0 = input$beta0
nobs = input$nobs
ncorrect = input$ncorrect
alpha1 = ncorrect+alpha0
beta1 = nobs-ncorrect+beta0
ymax = max(dbeta(seq1,shape1=alpha1,shape2=beta1),dbeta(seq1,shape1=alpha0,shape2=beta0))
seqL = seq(0,.5,length=1e3)
seqH = seq(.5,1,length=1e3)
plot(seq1,dbeta(seq1,shape1=alpha1,shape2=beta1),"l",ylab="density",xlab="theta",lwd=4,ylim=c(0,ymax+.5),xlim=c(0,1))
lines(seq1,dbeta(seq1,shape1=alpha0,shape2=beta0),lwd=2,col=2)
if(input$H1>0){
polygon(c(seqL,rev(seqL)),c(dbeta(seqL,shape1=alpha0,shape2=beta0),rep(0,1e3)), col = rgb(1,0,0,.6), border = "NA",ylim=c(0,ymax+.5),xlim=c(0,1))
polygon(c(seqL,rev(seqL)),c(dbeta(seqL,shape1=alpha1,shape2=beta1),rep(0,1e3)), col = rgb(0,0,0,.8), border = "NA",ylim=c(0,ymax+.5),xlim=c(0,1))
}
if(input$H2>0){
polygon(c(rev(seqH),seqH),c(rep(0,1e3),dbeta(seqH,shape1=alpha0,shape2=beta0)), col = rgb(1,0,0,.15), border = "NA")
polygon(c(rev(seqH),seqH),c(rep(0,1e3),dbeta(seqH,shape1=alpha1,shape2=beta1)), col = rgb(0,0,0,.2), border = "NA")
}
legend(0,ymax,c("unconstrained posterior","unconstrained prior"),lty=c(1,1),col=c(1,2),lwd=c(4,2))
dpost = dbeta(.5,shape1=alpha1,shape2=beta1)
dprior = dbeta(.5,shape1=alpha0,shape2=beta0)
if(input$H0>0){
lines(c(.5,.5),c(0,max(dprior,dpost)),lty=2,lwd=2)
lines(c(0,.5),c(dprior,dprior),lty=2,lwd=2)
lines(c(0,.5),c(dpost,dpost),lty=2,lwd=2)
text(0,dprior+.2,round(dprior,2))
text(0,dpost+.2,round(dpost,2))
}
})
output$BF0u <- renderText({
alpha0 = input$alpha0
beta0 = input$beta0
nobs = input$nobs
ncorrect = input$ncorrect
alpha1 = ncorrect+alpha0
beta1 = nobs-ncorrect+beta0
dpost = dbeta(.5,shape1=alpha1,shape2=beta1)
dprior = dbeta(.5,shape1=alpha0,shape2=beta0)
BF0u = round(dpost/dprior,2)
dpost = round(dpost,3)
dprior = round(dprior,3)
paste("BF(H0,Hu) = ", dpost, " / ", dprior, " = ",BF0u)
})
output$BF1u <- renderText({
alpha0 = input$alpha0
beta0 = input$beta0
nobs = input$nobs
ncorrect = input$ncorrect
alpha1 = ncorrect+alpha0
beta1 = nobs-ncorrect+beta0
ppost = pbeta(.5,shape1=alpha1,shape2=beta1)
pprior = pbeta(.5,shape1=alpha0,shape2=beta0)
BF1u = round(ppost/pprior,2)
ppost = round(ppost,3)
pprior = round(pprior,3)
paste("BF(H1,Hu) = ", ppost, " / ", pprior, " = ",BF1u)
})
output$BF2u <- renderText({
alpha0 = input$alpha0
beta0 = input$beta0
nobs = input$nobs
ncorrect = input$ncorrect
alpha1 = ncorrect+alpha0
beta1 = nobs-ncorrect+beta0
ppost = pbeta(.5,shape1=alpha1,shape2=beta1)
pprior = pbeta(.5,shape1=alpha0,shape2=beta0)
BF2u = round((1-ppost)/(1-pprior),2)
ppost = round(ppost,3)
pprior = round(pprior,3)
paste("BF(H2,Hu) = ", (1-ppost), " / ", (1-pprior), " = ",BF2u)
})
output$PHP0 <- renderText({
alpha0 = input$alpha0
beta0 = input$beta0
nobs = input$nobs
ncorrect = input$ncorrect
alpha1 = ncorrect+alpha0
beta1 = nobs-ncorrect+beta0
ppost = pbeta(.5,shape1=alpha1,shape2=beta1)
pprior = pbeta(.5,shape1=alpha0,shape2=beta0)
dpost = dbeta(.5,shape1=alpha1,shape2=beta1)
dprior = dbeta(.5,shape1=alpha0,shape2=beta0)
BF0u = dpost/dprior
BF1u = ppost/pprior
BF2u = (1-ppost)/(1-pprior)
PHP0 = BF0u/(BF0u+BF1u+BF2u)
PHP1 = BF1u/(BF0u+BF1u+BF2u)
PHP2 = BF2u/(BF0u+BF1u+BF2u)
paste("Pr(H0|data) = ",round(PHP0,3))
})
output$PHP1 <- renderText({
alpha0 = input$alpha0
beta0 = input$beta0
nobs = input$nobs
ncorrect = input$ncorrect
alpha1 = ncorrect+alpha0
beta1 = nobs-ncorrect+beta0
ppost = pbeta(.5,shape1=alpha1,shape2=beta1)
pprior = pbeta(.5,shape1=alpha0,shape2=beta0)
dpost = dbeta(.5,shape1=alpha1,shape2=beta1)
dprior = dbeta(.5,shape1=alpha0,shape2=beta0)
BF0u = dpost/dprior
BF1u = ppost/pprior
BF2u = (1-ppost)/(1-pprior)
PHP0 = BF0u/(BF0u+BF1u+BF2u)
PHP1 = BF1u/(BF0u+BF1u+BF2u)
PHP2 = BF2u/(BF0u+BF1u+BF2u)
paste("Pr(H1|data) = ",round(PHP1,3))
})
output$PHP2 <- renderText({
alpha0 = input$alpha0
beta0 = input$beta0
nobs = input$nobs
ncorrect = input$ncorrect
alpha1 = ncorrect+alpha0
beta1 = nobs-ncorrect+beta0
ppost = pbeta(.5,shape1=alpha1,shape2=beta1)
pprior = pbeta(.5,shape1=alpha0,shape2=beta0)
dpost = dbeta(.5,shape1=alpha1,shape2=beta1)
dprior = dbeta(.5,shape1=alpha0,shape2=beta0)
BF0u = dpost/dprior
BF1u = ppost/pprior
BF2u = (1-ppost)/(1-pprior)
PHP0 = BF0u/(BF0u+BF1u+BF2u)
PHP1 = BF1u/(BF0u+BF1u+BF2u)
PHP2 = BF2u/(BF0u+BF1u+BF2u)
paste("Pr(H2|data) = ",round(PHP2,3))
})
output$PrHP0 <- renderText({
paste("Pr(H0) = ",.333)
})
output$PrHP1 <- renderText({
paste("Pr(H1) = ",.333)
})
output$PrHP2 <- renderText({
paste("Pr(H2) = ",.333)
})
})
- Press “Run App” in Rstudio. Enjoy :)
joris mulder 