An interactive Bayesian Z test using the shiny package in R

This is a shiny app for performing a Bayesian Z test of

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 Z test”.
• Create a R-script file in Rstudio and call the file: ui.R.
• Save ui.R in the folder “Bayesian Z test”.
• Create a R-script file in Rstudio called: server.R.
• Save server.R in the folder “Bayesian Z test”.
• Paste the syntax below to ui.R. (scroll down for the next step).

library(shiny)

shinyUI(fluidPage(
 
 titlePanel("Bayesian multiple Z test: H0:mu=0 vs H1:mu<0 vs H2:mu>0"),
 "Made by Joris Mulder (Tilburg University, www.jorismulder.com).",
 
 sidebarLayout(
 sidebarPanel(
 sliderInput("sigma0",
 "measurement error (sigma)",
 min = 0.001,
 max = 5.000,
 value = 1),
 numericInput("samplemean",
 "sample mean",
 value = 0.8),
 sliderInput("nobs",
 "Number of observations:",
 min = 1,
 max = 100,
 value = 10),
 sliderInput("priorsd0",
 "prior standard deviation",
 min = 0.001,
 max = 10.000,
 value = 1),
 checkboxInput("viewprior", "view prior",T),
 checkboxInput("viewdata", "view likelihood",F),
 checkboxInput("viewpost", "view posterior",F),
 checkboxInput("viewConfi", "view 95%-confidence interval",F),
 checkboxInput("viewCredi", "view 95%-credibility interval",F),
 checkboxInput("H0", "H0 versus Hu",F),
 checkboxInput("H1", "H1 versus Hu",F),
 checkboxInput("H2", "H2 versus Hu",F)
 ),
 
 mainPanel(
 plotOutput("distPlot"),
 
 strong("BF of constrained hypothesis (H0, H1, or H2) versus unconstrained hypothesis Hu"),
 tableOutput("BFtu"),
 
 strong("BF Table"),
 tableOutput("BFtest"),
 
 strong("Prior probabilities hypotheses"),
 tableOutput("PriorHypo"),
 
 strong("Posterior probabilities hypotheses"),
 tableOutput("PostHypo")
 
 
 )
 )
)
)

• Paste the syntax below to server.R. (scroll down for the final step).

library(shiny)

shinyServer(
 function(input, output){
 
 output$distPlot = renderPlot({
 seq1 = seq(-4,4,length=1e3+1)
 sigma0 = input$sigma0
 samplemean = input$samplemean
 nobs = input$nobs
 priorsd0 = input$priorsd0
 
 priorvar = priorsd0**2
 datavar = sigma0**2/nobs
 postvar = 1/(1/datavar + 1/priorvar)
 postmean = (samplemean/datavar)/(1/datavar+1/priorvar)
 
 ymax = dnorm(0,mean=0,sd=sqrt(postvar))
 
 seqL = seq(min(seq1),0,length=1e3)
 seqH = seq(0,max(seq1),length=1e3)
 
 plot(-100,-100,ylab="",xlab="mu",lwd=4,ylim=c(0,ymax*1.1),xlim=c(min(seq1),max(seq1)))
 legend(2.5,ymax*1.1,legend=c("prior","likelihood","posterior"),lwd=rep(2,3),col=c(2,1,4),lty=c(2,4,1))
 
 if(input$viewprior>0){
 lines(seq1,dnorm(seq1,mean=0,sd=priorsd0),lwd=2,col=2,lty=2)
 
 if(input$H0>0){
 dprior = dnorm(0,mean=0,sd=priorsd0)
 lines(c(0,0),c(0,dprior),lty=3,lwd=2,col=2)
 lines(c(min(seq1),0),rep(dprior,2),lty=3,lwd=2,col=2)
 text(-4,dprior+.05,round(dprior,2),col=2)
 points(0,dprior,col="white",lwd=2,pch=20)
 points(0,dprior,col=2,lwd=2,pch=21)
 }
 
 if(input$H1>0){
 polygon(c(seqL,rev(seqL)),c(dnorm(seqL,mean=0,sd=priorsd0),rep(0,1e3)), col = rgb(1,0,0,.6), border = "NA")
 prob1 = pnorm(0,mean=0,sd=priorsd0)
 densprior0 = dnorm(0,mean=0,sd=priorsd0)
 text(-.4,densprior0*.25,round(prob1,2),col=2,lwd=2)
 }
 
 if(input$H2>0){
 polygon(c(seqH,rev(seqH)),c(dnorm(seqH,mean=0,sd=priorsd0),rep(0,1e3)), col = rgb(1,0,0,.3), border = "NA")
 prob2 = 1-pnorm(0,mean=0,sd=priorsd0)
 densprior0 = dnorm(0,mean=0,sd=priorsd0)
 text(.4,densprior0*.25,round(prob2,2),col=2,lwd=2)
 }
 }
 if(input$viewdata>0){
 lines(seq1,dnorm(seq1,mean=samplemean,sd=sqrt(datavar)),lwd=2,col=1,lty=4)
 
 if(input$viewConfi>0){
 UB = samplemean + 1.96*sqrt(datavar)
 LB = samplemean - 1.96*sqrt(datavar)
 conf95 = seq(UB,LB,length=1e3)
 points(conf95,rep(0,1e3),pch=20)
 text(LB,.1,round(LB,2),col=1)
 text(UB,.1,round(UB,2),col=1)
 }
 }
 if(input$viewpost>0){
 lines(seq1,dnorm(seq1,mean=postmean,sd=sqrt(postvar)),lwd=2,col=4)
 if(input$H0>0){
 dpost = dnorm(0,mean=postmean,sd=sqrt(postvar))
 lines(c(0,0),c(0,dpost),lty=3,lwd=2,col=4)
 lines(c(min(seq1),0),rep(dpost,2),lty=3,lwd=2,col=4)
 text(-4,dpost+.05,round(dpost,2),col=4)
 points(0,dpost,col="white",lwd=2,pch=20)
 points(0,dpost,col=4,lwd=2,pch=21)
 }
 
 if(input$H1>0){
 polygon(c(seqL,rev(seqL)),c(dnorm(seqL,mean=postmean,sd=sqrt(postvar)),rep(0,1e3)), col = rgb(0,0,1,.6), border = "NA")
 prob1 = pnorm(0,mean=postmean,sd=sqrt(postvar))
 denspost0 = dnorm(0,mean=postmean,sd=sqrt(postvar))
 text(-.4,denspost0*.5,round(prob1,2),col=4,lwd=2)
 }
 
 if(input$H2>0){
 polygon(c(seqH,rev(seqH)),c(dnorm(seqH,mean=postmean,sd=sqrt(postvar)),rep(0,1e3)), col = rgb(0,0,1,.3), border = "NA")
 prob2 = 1-pnorm(0,mean=postmean,sd=sqrt(postvar))
 denspost0 = dnorm(0,mean=postmean,sd=sqrt(postvar))
 text(.4,denspost0*.5,round(prob2,2),col=4,lwd=2)
 }
 
 if(input$viewCredi>0){
 UB = postmean + 1.96*sqrt(postvar)
 LB = postmean - 1.96*sqrt(postvar)
 conf95 = seq(UB,LB,length=1e3)
 points(conf95,rep(0,1e3),pch=20,col=4)
 text(LB,.2,round(LB,2),col=4)
 text(UB,.2,round(UB,2),col=4)
 }
 }
 }
)
 
 output$BFtu <- renderTable({
 
 sigma0 = input$sigma0
 samplemean = input$samplemean
 nobs = input$nobs
 priorsd0 = input$priorsd0
 
 priorvar = priorsd0**2
 datavar = sigma0**2/nobs
 postvar = 1/(1/datavar + 1/priorvar)
 postmean = (samplemean/datavar)/(1/datavar+1/priorvar)
 
 BF0u = dnorm(0,mean=postmean,sd=sqrt(postvar))/dnorm(0,mean=0,sd=sqrt(priorvar))
 BF1u = pnorm(0,mean=postmean,sd=sqrt(postvar))/pnorm(0,mean=0,sd=sqrt(priorvar))
 BF2u = (1-pnorm(0,mean=postmean,sd=sqrt(postvar)))/(1-pnorm(0,mean=0,sd=sqrt(priorvar)))
 
 results <- data.frame(Hu=c(BF0u,BF1u,BF2u),row.names=c("H0","H1","H2"))
 },
 digits = 3,include.rownames=T, sanitize.rownames.function = function(x) paste('<b>',x,'</b>', sep ='')
 )
 
 output$BFtest <- renderTable({
 
 sigma0 = input$sigma0
 samplemean = input$samplemean
 nobs = input$nobs
 priorsd0 = input$priorsd0
 
 priorvar = priorsd0**2
 datavar = sigma0**2/nobs
 postvar = 1/(1/datavar + 1/priorvar)
 postmean = (samplemean/datavar)/(1/datavar+1/priorvar)
 
 BF0u = dnorm(0,mean=postmean,sd=sqrt(postvar))/dnorm(0,mean=0,sd=sqrt(priorvar))
 BF1u = pnorm(0,mean=postmean,sd=sqrt(postvar))/pnorm(0,mean=0,sd=sqrt(priorvar))
 BF2u = (1-pnorm(0,mean=postmean,sd=sqrt(postvar)))/(1-pnorm(0,mean=0,sd=sqrt(priorvar)))
 
 BFvec = c(BF0u,BF1u,BF2u)
 
 BFmat = t(matrix(c(BF0u/BFvec,BF1u/BFvec,BF2u/BFvec),ncol=3))
 
 results <- data.frame(H0=BFmat[,1],H1=BFmat[,2],H2=BFmat[,3],row.names=c("H0","H1","H2"))
 },
 digits = 3, sanitize.rownames.function = function(x) paste('<b>',x,'</b>', sep ='')
 )
 
 output$PriorHypo <- renderTable({
 results <- data.frame(H0=.333,H1=.333,H2=.333)
 },
 digits = 3,include.rownames=F
 )
 
 output$PostHypo <- renderTable({
 
 sigma0 = input$sigma0
 samplemean = input$samplemean
 nobs = input$nobs
 priorsd0 = input$priorsd0
 
 priorvar = priorsd0**2
 datavar = sigma0**2/nobs
 postvar = 1/(1/datavar + 1/priorvar)
 postmean = (samplemean/datavar)/(1/datavar+1/priorvar)
 
 BF0u = dnorm(0,mean=postmean,sd=sqrt(postvar))/dnorm(0,mean=0,sd=sqrt(priorvar))
 BF1u = pnorm(0,mean=postmean,sd=sqrt(postvar))/pnorm(0,mean=0,sd=sqrt(priorvar))
 BF2u = (1-pnorm(0,mean=postmean,sd=sqrt(postvar)))/(1-pnorm(0,mean=0,sd=sqrt(priorvar)))
 
 BFvec = c(BF0u,BF1u,BF2u)
 
 results <- data.frame(H0=BF0u/sum(BFvec),H1=BF1u/sum(BFvec),H2=BF2u/sum(BFvec))
 },
 digits = 3,include.rownames=F
 )
 }
)

• Press “Run App” in Rstudio. Enjoy :)