Hello R community,
I am currently working on making a Shiny app that can store an accurate test statistic for any of these three distributions (Poisson, Binomial, or Uniform). I am pretty sure the code works for the uniform distribution, but I am struggling to find code that works for any histogram produced under the rbinom and rpois functions. The part I'm needing help with is the "teststat" code.
-
For the Poisson test stat, I am having trouble when a histogram has a break with 0 versus a histogram with first break of 1.
-
For the Binomial test stat, I am having trouble making an accurate test statistic with the breaks occurring at 0.5, 1.5, 2.5, etc...
In both cases I am trying to produce a chi-squared test statistic, which is sum(((observed-expected)^2)/expected).
Any help would be appreciated. Thank you.
library(shiny)
# Define UI for application that draws a histogram
ui <- fluidPage(
titlePanel("Goodness-of-Fit Testing Hypothesis Practice"),
sidebarLayout(
sidebarPanel(
radioButtons("distribution", label = h3("Pick the type of distribution you would like to practice"), choices = c("Poisson" = "pois", "Binomial" = "binom", "Uniform" = "unif")),
actionButton(inputId = "submit_dist", label = "Press for new distribution"),
sliderInput("bins", h3("Enter the number of bins"), min = 5, max = 15, value = 10, step = 5),
numericInput("guess", "What is your calculated chi-squared test statistic?", ""),
actionButton(inputId = "submit_guess", label = "Press to submit"),
textOutput("actual"),
tags$style("body{background-color: white; color: black}")
),
mainPanel(
plotOutput("distPlot"),
textOutput("extrainfo"),
textOutput("result")
)
),
sidebarLayout(
sidebarPanel(
textOutput("forprofsmith"),
textOutput("acknowledgment"),
textOutput("credit"),
width = 3),
mainPanel(
tableOutput("table"),
textOutput("totalguesses"),
textOutput("score"),
textOutput("percentage")
)
)
)
# Define server logic required to make a histogram
server <- function(input, output){
x_val <- eventReactive(input$submit_dist, {
simulate(input$distribution)})
simulate <- function(dist, n = 250) {
switch(dist, pois = rpois(n, lambda1), binom = rbinom(n, size1, sampleprop), unif = runif(n, max = lambda1), stop("Unknown 'dist'"))
}
teststat <- eventReactive(input$submit_dist, {
findstat(input$distribution)
})
findstat <- function(dist, n = 250) {
data2 <- hist(x_val(), breaks = input$bins, plot = FALSE)
df <- data.frame(counts = data2$counts)
counter <- 2
while(counter < length(data2$breaks)){
otherexpvals[counter-1] <- data2$breaks[counter] + c(1:(data2$breaks[counter+1] - data2$breaks[counter]))
counter <- counter + 1
}
expval <- c(sum(dpois(data2$breaks[1]:data2$breaks[2], lambda1)), dpois(otherexpvals, lambda = lambda1))
unifstat <- chisq.test(df$counts, p = dunif(0:((length(df$counts))-1), max = length(df$counts)))$statistic
if(data2$breaks[2] < 1){
binomstat <- sum((df$counts[df$counts != 0] - (c(dbinom(0, size1, sampleprop), dbinom(1:length(data2$counts[data2$counts != 0])-1, size1, sampleprop)))*250)^2/((c(dbinom(0, size1, sampleprop), dbinom(1:length(data2$counts[data2$counts != 0])-1, size1, sampleprop)))*250))
}
else{
binomstat <- sum((df$counts[df$counts != 0] - (c(sum(dbinom(0:1, size1, sampleprop)), dbinom(2:length(data2$counts[data2$counts != 0]), size1, sampleprop)))*250)^2/((c(sum(dbinom(0:1, size1, sampleprop)), dbinom(2:length(data2$counts[data2$counts != 0]), size1, sampleprop)))*250))
}
switch(dist, pois = sum(((df$counts - expval*250)^2)/(expval*250)), binom = binomstat, unif = unifstat, stop("Unknown 'dist'"))
}
counter <- reactiveValues(correct = 0, guess = 0)
output$distPlot <- renderPlot({
hist(x_val(),
main = "Distribution",
xlab = "Sample Values",
ylab = "Frequency Count",
col = "darkgray",
border = "white",
labels = T,
breaks = input$bins)})
################################################
seed1 <- sample.int(1000000, 1)
set.seed(seed1)
lambda1 <- sample.int(5:15, 1)
size1 <- sample.int(10:100, 1)
sampleprop <- sample(c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9), 1)
################################################
output$table <- renderTable({
data1 <- hist(x_val(), breaks = input$bins, plot = FALSE)
data.frame(counts = data1$counts)
})
output$extrainfo <- eventReactive(input$submit_dist, {
extra(input$distribution)})
extra <- function(dist) {
switch(dist, pois = paste("The expected frequency is ", lambda1), binom = paste("The expected proportion is ", sampleprop, "and the size of each trial is", size1), unif = NULL, stop("Unknown 'dist'"))
}
output$result <- eventReactive(input$submit_guess, {
q <- input$guess
t1 <- teststat()
if(q < (t1 - 0.5)){
counter$guess <- counter$guess + 1.0
return("That estimate is too small. Try again!")
}
else if(q > (t1 + 0.5)){
counter$guess <- counter$guess + 1.0
return("That estimate is too big. Try again!")
}
else{
counter$guess <- counter$guess + 1.0
counter$correct <- counter$correct + 1.0
return("That is a correct estimation! Good job! Try again.")
}
})
output$actual <- eventReactive(input$submit_guess, {
q <- input$guess
t1 <- teststat()
if(q < (t1 - 0.5)){
paste("The actual test statistic is --.--")
}
else if(q > (t1 + 0.5)){
paste("The actual test statistic is --.--")
}
else{
paste("The true test statistic is ", (round(t1, digits = 2)))
}
})
output$totalguesses <- renderText({paste("Total guesses is ", counter$guess)})
output$score <- renderText({paste("Your score is ", counter$correct)})
output$percentage <- renderText({paste("Your percentage correct is",
(if (!isTruthy((counter$correct/counter$guess))) {
return(NULL)
}
else{
round((counter$correct/counter$guess)*100, digits = 2)
}),"%")})
output$forprofsmith <- renderText("Made for Spring 2020. ")
output$acknowledgment <- renderText("Designed by ")
output$credit <- renderText("")
}
# Run the application
shinyApp(ui = ui, server = server)