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