R/ggproptest.R

In gginference: Visualise the Results of Inferential Statistics using 'ggplot2'

#' @title 
#' Plot test of Equal or Given Proportions
#'  
#' @description
#' Visualise \code{\link{prop.test}} for testing the null that the proportions (probabilities of success) 
#' in several groups are the same, or that they equal certain given values.
#' 
#' @usage ggproptest(t, alpha=0.05,colaccept="lightsteelblue1", 
#' colreject="gray84", colstat="navyblue")
#' 
#' @param t a list result of \code{\link{prop.test}} of "htest" class
#' @param alpha alpha level for ploting distribution, when prop.test is used on more than 2 samples
#' @param colaccept color the acceptance area of the test
#' @param colreject color for the rejection area of the test
#' @param colstat color for the test statistic vline
#' 
#' @import ggplot2
#'
#' @rdname ggproptest
#' 
#' @export
#' 
#' @examples 
#' 
#' x  <- c(5, 8, 12)
#' y <- c(8, 9, 13)
#' pr_test <- prop.test(x, y)
#' pr_test
#' ggproptest(pr_test)
#' 

ggproptest <- function(t, alpha = 0.05, 
                       colaccept = "lightsteelblue1", colreject = "gray84", colstat = "navyblue"){
  df <- t$parameter 
  stat <- t$statistic
  
  textt <- t$method
  st <- strsplit(textt, split = "-")
  if (as.numeric(st[[1]][1]) <= 2) {
    level <- as.numeric(attributes(t$conf.int))
    #make a list of points for the plot
    points <- seq(-0.25, 8, length = 10000)
  } else {
    level = 1 - alpha
    #make a list of points for the plot
    points <- seq(-0.25, 20, length = 10000)
  }
  
  
  #set the upper bound
  ub <- stats::qchisq(p = level, df = df)
  #limit the points based on the upper bound
  limits <- points <= ub
  
  #make dataframes for ggplot
  dfpoly1 <- data.frame(x = c(points, ub), y = c(stats::dchisq(points, df),0))
  dfpoly2 <- rbind(dfpoly1[limits,], c(ub,0))
  
  plotpropt <- ggplot(data = data.frame(points), aes(points)) +
    stat_function(fun = stats::dchisq, args = list(df = df), col = "white") +
    geom_polygon(data = dfpoly1, aes(x, y), fill = colreject) +
    geom_polygon(data = dfpoly2, aes(x, y), fill = colaccept) +
    labs(title = "Chi-squared distribution Vs test statistic",
         subtitle = paste("based on", t$method),
         x = paste("Chi-squared distributions with", df, "degrees of freedom"),
         caption = paste("alpha=", 1 - level)) +
    scale_y_continuous(breaks = NULL) +
    ylab("") +
    #add the statistics line
    geom_vline(data = data.frame(x = stat), aes(xintercept = data.frame(x = stat)$x), colour = colstat, lwd = 1) +
    geom_text(data = data.frame(x = stat), aes(x = data.frame(x = stat)$x, y = 0, 
                                               label = paste("test statistic=", round(data.frame(x = stat)$x, 3))),
              size = 4, angle = 90, vjust = -0.4, hjust = 0, colour = colstat) +
    geom_vline(aes(xintercept = ub), linetype = 2, alpha = 0) +  
    geom_text(aes(x = ub, y = -0.02), label = round(ub,3), vjust = 0.3) +
    theme_classic()
  #plot it
  plotpropt
}