An interactive 2-sided p-value test of a proportion using the shiny package in R

The 2-sided test of a proportion ‘theta’ can be written as

• H0: theta = theta0 (i.e., the proportion is equal to a null value)
• H1: theta != theta0 (i.e., the proportion is unequal to a null value)

The shiny app computes the 2-tailed p-value for this test based on a binomial sampling plan (i.e., perform a certain experiment N number of times out of which we observe x successes) and a negative binomial sampling plan (i.e., conduct an experiment until x successes have been observed for which we needed N number of experiments in total).

• Install R.
• Install Rstudio.
• Via Tools -> install packages -> shiny.
• In the console, write: “library(shiny)”.
• Create folder named “P-value sampling plan”.
• Create a R-script file in Rstudio and call the file: ui.R.
• Save ui.R in the folder “P-value sampling plan”.
• Create a R-script file in Rstudio called: server.R.
• Save server.R in the folder “P-value sampling plan”.
• Paste the syntax below to ui.R. (scroll down for the next step).

#astro posterior

library(shiny)

shinyUI(fluidPage(
  
  # Application title
  titlePanel("Classical test of H0: theta=.5 versus H1: theta!=.5"),
  
  # Sidebar with a slider input
  sidebarLayout(
    sidebarPanel(
      sliderInput("nobs",
                  "Number of experiments:",
                  min = 1,
                  max = 100,
                  value=60),
      sliderInput("ncorrect",
                  "Number of correct answers:",
                  min = 1,
                  max = 100,
                  value = 38),
      sliderInput("theta0",
                  "H0:theta=",
                  min = .001,
                  max = .99,
                  value = .5)
#      checkboxInput("pvalue", "View 2-sided p-value",F)
    ),
    
    # Show a plot of the generated distribution
    mainPanel(
      plotOutput("distPlot1"),
      textOutput("pvalue1"),
      plotOutput("distPlot2"),
      textOutput("pvalue2")
    )
  )
))

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

library(shiny)

shinyServer(function(input, output) {
    
  output$distPlot1 = renderPlot({
    
    theta0 = input$theta0
    N = input$nobs
    x = input$ncorrect
    xright = x
    if(xN/2){xleft=N-x}
    probs = dbinom(0:N,size=N,prob=theta0)
    pvalue1 = sum(probs[1:(xleft+1)])
    paste("The p-value is ",round(pvalue1*2,3)) 
    
    text(x=xleft+2,y=max(probs)/10,round(pvalue1,3),pos=2)
    text(x=xright+10,y=max(probs)/10,round(pvalue1,3),pos=4)
    

  })

output$pvalue1 N/2){xleft=N-x}
  probs = dbinom(0:N,size=N,prob=input$theta0)
  pvalue1 = sum(probs[1:(xleft+1)])
  paste0("The p-value is ",round(pvalue1*2,3),".")  
})

output$distPlot2 = renderPlot({
  
  dnegbinom = function(n,xx,prob){
    #stop until x successes. total number of experiments is n.
    #prob is probability of success per experiment
    exp(xx*log(prob)+(n-xx)*log(1-prob)+lfactorial(n-1)-lfactorial(xx-1)-lfactorial(n-xx))
  }
  
  theta0 = input$theta0
  N = input$nobs
  x = input$ncorrect
  
  maxx = x/theta0*2
  probs = dnegbinom(n=x:maxx,x=x,prob=theta0)
  barplot(probs,names.arg=x:maxx,xlim=c(0,maxx),ylim=c(0,max(probs)))
  if(x/N>theta0){
    par(new=T)
    #at least as extreme as N is N and less
    barplot(c(probs[1:(N-x+1)],rep(0,maxx-N+x)),xlim=c(0,maxx),ylim=c(0,max(probs)),col=2)
    pleft = sum(probs[1:(N-x+1)])
    text(N-x-8,max(probs)/10,round(pleft,3))
    pxleft = 1-x/N
    for(nn in (N+1):maxx){
      if(x/nn<=pxleft){
        Nright = nn
        break
      }
    }
    par(new=T)
    #at least as extreme as N is Nleft and more
    barplot(c(rep(0,(Nright-x)),probs[(Nright-x+1):(maxx-x+1)]),xlim=c(0,maxx),ylim=c(0,max(probs)),col=2)
    pright = sum(probs[(Nright-x+1):(maxx-x+1)])
    text(N+x,max(probs)/10,round(pright,3))
  }
  if(x/N<theta0){
    par(new=T)
    barplot(c(rep(0,(N-x-1)),probs[(N-x):(maxx-x+1)]),xlim=c(0,maxx),ylim=c(0,max(probs)),col=2)
    pright = sum(probs[(N-x):(maxx-x+1)])
    text(N+x,max(probs)/10,round(pright,3))
    pxright = 1-x/N
    for(nn in (x+1):N){
      if(x/nn<=pxright){
        Nleft = nn-1
        break
      }
    }
    par(new=T)
    barplot(c(probs[1:(Nleft-x)],rep(0,maxx-Nleft+x)),xlim=c(0,maxx),ylim=c(0,max(probs)),col=2)
    pleft = sum(probs[1:(Nleft-x)])
    text(N-x-18,max(probs)/10,round(pleft,3))
  }
  
})


output$pvalue2 theta0){
#    par(new=T)
    #at least as extreme as N is N and less
#    barplot(c(probs[1:(N-x+1)],rep(0,maxx-N+x)),xlim=c(0,maxx),ylim=c(0,max(probs)),col=2)
    pleft = sum(probs[1:(N-x+1)])
#    text(N-x-8,max(probs)/10,round(pleft,3))
    pxleft = 1-x/N
    for(nn in (N+1):maxx){
      if(x/nn<=pxleft){
        Nright = nn
        break
      }
    }
#    par(new=T)
    #at least as extreme as N is Nleft and more
#    barplot(c(rep(0,(Nright-x)),probs[(Nright-x+1):(maxx-x+1)]),xlim=c(0,maxx),ylim=c(0,max(probs)),col=2)
    pright = sum(probs[(Nright-x+1):(maxx-x+1)])
#    text(N+x,max(probs)/10,round(pright,3))
    pvalue2 = pleft+pright
  }
  if(x/N<theta0){
#    par(new=T)
#    barplot(c(rep(0,(N-x-1)),probs[(N-x):(maxx-x+1)]),xlim=c(0,maxx),ylim=c(0,max(probs)),col=2)
    pright = sum(probs[(N-x):(maxx-x+1)])
#    text(N+x,max(probs)/10,round(pright,3))
    pxright = 1-x/N
    for(nn in (x+1):N){
      if(x/nn<=pxright){
        Nleft = nn-1
        break
      }
    }
    par(new=T)
#    barplot(c(probs[1:(Nleft-x)],rep(0,maxx-Nleft+x)),xlim=c(0,maxx),ylim=c(0,max(probs)),col=2)
    pleft = sum(probs[1:(Nleft-x)])
#    text(N-x-18,max(probs)/10,round(pleft,3))
    pvalue2 = pleft+pright
  }
  
  paste0("The p-value is ",round(pvalue2,3),".")  
})

})

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