R/ggvartest.R

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

#' @title 
#' F test plot
#'  
#' @description
#' Visualise F test to compare two variances
#' 
#' @usage ggvartest(t, colaccept = "lightsteelblue1", 
#' colreject = "gray84", colstat = "navyblue")
#' 
#' @param t a list result of \code{\link{var.test}} of "htest" class
#' @param colaccept color the acceptance area of the test, see  \code{\link{colors}}
#' @param colreject color for the rejection area of the test
#' @param colstat color for the test statistic vline
#' 
#' 
#' @import ggplot2
#' @importFrom rlang .data
#' 
#' @rdname ggvartest
#' 
#' @export
#' 
#' @examples 
#' 
#' x <- rnorm(50, mean = 0, sd = 2)
#' y <- rnorm(30, mean = 1, sd = 1)
#' var_test <- var.test(x, y)
#' var_test
#' ggvartest(var_test)
#' 

ggvartest <- function(t, colaccept = "lightsteelblue1", 
                      colreject = "gray84", colstat = "navyblue"){
  
  #1-alpha,is needed for F distribution
  level <- as.numeric(attributes(t$conf.int))
  
  #degrees of freedom for variance test
  if (t$alternative == "two.sided") {
    if (t$statistic > 1) {
      df1 <- t$parameter[1]
      df2 <- t$parameter[2]
      t$statistic <- t$statistic
    }
    else {
      df1 <- t$parameter[2]
      df2 <- t$parameter[1]
      t$statistic <- 1/t$statistic
    }
  } else if (t$alternative == "greater") {
    df1 <- t$parameter[1]
    df2 <- t$parameter[2]
  } else {
    df1 <- t$parameter[2]
    df2 <- t$parameter[1]
  }
  
  #sequence of points for F-distribution
  points <- seq(-0.25, 7, length = 10000)
  
  #upper bound for all alternatives
  ub <- stats::qf(p = level, df1 = df1, df2 = df2)
  limits <- points <= ub
  
  #make dataframes for ggplot
  dfpoly1 <- data.frame(x = c(points, ub), y = c(stats::df(points, df1, df2), 0))
  dfpoly2 <- rbind(dfpoly1[limits,], c(ub,0))
  
  #plot it
  ggplot(data = data.frame(points), aes(points)) +
    stat_function(
      fun = stats::df, 
      n = 101, 
      args = list(df1 = df1, df2 = df2), 
      colour = "white") +
    
    geom_polygon(
      data = dfpoly1,
      aes(x, y), 
      fill = colreject) +
    
    geom_polygon(
      data = dfpoly2, 
      aes(x, y), 
      fill = colaccept) +
    
    labs(
      title = "F distribution Vs test statistic",
      subtitle = paste("Alternative hypothesis:", t$alternative),
      x = paste("F distributions with", df1, "&", df2, "degrees of freedom"),
      caption = paste("alpha=", 1 - level)) +
    
    scale_y_continuous(breaks = NULL) +
    
    ylab("") +
    
    #add the statistics line
    geom_vline(
      data = data.frame(x = t$statistic),
      aes(xintercept = data.frame(x = t$statistic)$x),
      colour = colstat, lwd = 1) +
    
    geom_text(
      data = data.frame(x = t$statistic), 
      aes(x = data.frame(x = t$statistic)$x,
          y = 0, 
          label = paste("test statistic=", round(data.frame(x = t$statistic)$x, 3))),
      size = 4, angle = 90, vjust = -0.4, hjust = 0, colour = colstat) +
    
    theme_classic() +
    
    #add the F-score
    geom_vline(aes(xintercept = ub), linetype = 2, alpha = 0) +  
    
    geom_text(aes(x = ub, y = -0.025), label = round(ub,3), vjust = 0.3)
  
}