R/plot.ttest.R

In nhstplot: Plot Null Hypothesis Significance Tests

#' Illustrate a one- or two-tailed t test graphically.
#'
#' This function plots the density probability distribution of a t statistic, with appropriate vertical cutlines at the t value. The p-value and the observed t value are plotted. Although largely customizable, only two arguments are required (the observed t statistic and the degrees of freedom) for a two-tailed t test. The optional argument \code{tails = "one"} plots a one-tailed test plot (the tail is on the left or right, depending on the sign of the t statistic).
#'
#' @param t A numeric value indicating the observed t statistic. Alternatively, you can pass an object of class \code{htest} created by the function \code{t.test()} or \code{cor.test()}.
#' @param df A numeric value indicating the degrees of freedom. This argument is optional if you are using an \code{htest} object as the \code{t} argument.
#' @param tails A character that indicates whether to plot a one (\code{"one"}) or two (\code{"two"}) tailed t-test (optional). By default, a two-tailed test is plotted.
#' @param blank A logical that indicates whether to hide (\code{blank = TRUE}) the test statistic value, p value and cutline. The corresponding colors are actually only made transparent when \code{blank = TRUE}, so that the output is scaled exactly the same (this is useful and especially intended for step-by-step explanations).
#' @param xmax A numeric including the maximum for the x-axis. Defaults to \code{"auto"}, which scales the plot automatically (optional).
#' @param title A character or expression indicating a custom title for the plot (optional).
#' @param xlabel A character or expression indicating a custom title for the x axis (optional).
#' @param ylabel A character or expression indicating a custom title for the y axis (optional).
#' @param fontfamily A character indicating the font family of all the titles and labels (e.g. \code{"serif"} (default), \code{"sans"}, \code{"Helvetica"}, \code{"Palatino"}, etc.) (optional).
#' @param colormiddle A character indicating the color for the "middle" area under the curve (optional).
#' @param colorsides A character indicating the color for the "side(s)" area(s) under the curve (optional).
#' @param colormiddlecurve A character indicating the color for the "middle" part of the curve (optional).
#' @param colorsidescurve A character indicating the color for the "side(s)" part of the curve (optional).
#' @param colorcut A character indicating the color for the cut line at the observed test statistic (optional).
#' @param colorplabel A character indicating the color for the label of the p-value (optional). By default, for color consistency, this color is the same as color of \code{colorright}.
#' @param theme A character indicating one of the predefined color themes. The themes are \code{"default"} (light blue and red), \code{"blackandwhite"}, \code{"whiteandred"}, \code{"blueandred"}, \code{"greenandred"} and \code{"goldandblue"}) (optional). Supersedes \code{colormiddle} and \code{colorsides} if another argument than \code{"default"} is provided.
#' @param signifdigitst A numeric indicating the number of desired significant figures reported for the t label (optional).
#' @param curvelinesize A numeric indicating the size of the curve line (optional).
#' @param cutlinesize A numeric indicating the size of the cut line(s) (optional). By default, the size of the curve line is used.
#' @param p_value_position A numeric vector of length 2, indicating the x and y coordinates of the p-value label. By default, the position is set to \code{"auto"}. Note that the absolute value is used, and the sign is ignored. The position is set to the right if the test statistic value is positive, to the left if the test statistic value is negative, and on both sides if a two tailed test is plotted.
#' @return A plot with the density of probability of t under the null hypothesis, annotated with the observed test statistic and the p-value.
#' @examples
#' #Making a t test plot with a t value of 2 and df of 10
#' plotttest(t = 2, df = 10)
#'
#' #The same plot without the t or p value
#' plotttest(2,10, blank = TRUE)
#'
#' #Plotting a one-tailed test using the "tails" parameter.
#' plotttest(t = 2, df = 10, tails = "one")
#'
#' #Using t.test() as an input
#' test <- t.test(rnorm(10), rnorm(10))
#' plotttest(test)
#'
#' #Using cor.test() as an input
#' test <- cor.test(rnorm(10), rnorm(10))
#' plotttest(test)
#'
#' @author Nils Myszkowski <nmyszkowski@pace.edu>
#' @export plotttest

plotttest <- function(t, df = t$parameter, tails = "two", blank = FALSE, xmax = "auto", title = "t Test", xlabel = "t", ylabel = "Density of probability\nunder the null hypothesis", fontfamily = "serif", colormiddle = "aliceblue", colorsides = "firebrick3", colormiddlecurve = "black", colorsidescurve = "black", colorcut = "black", colorplabel = colorsides, theme = "default", signifdigitst = 3, curvelinesize = .4, cutlinesize = curvelinesize, p_value_position = "auto") {
  x=NULL


  # If t is a t.test() object, then mine it to get t and df
  if ("htest" %in% class(t)) {
    df <- t$parameter
    if (t$alternative != "two.sided") {tails = "one"} else {tails = "two"}
    t <- t$statistic
  }

  #Unname inputs (can cause issues)
  t <- unname(t)
  df <- unname(df)

  #Create a function to restrict plotting areas to specific bounds of x
  area_range <- function(fun, min, max) {
    function(x) {
      y <- fun(x)
      y[x < min | x > max] <- NA
      return(y)
    }
  }

  # Function to format p value
  p_value_format <- function(p) {
    if (p < .001) {"< .001"} else
      if (p > .999) {"> .999"} else
        paste0("= ", substr(sprintf("%.3f", p), 2, 5))
  }

  #Store the t value provided as argument, used only for one tailed tests, to decide whether to plot left tailed or right tailed
  originalt <- t
  #Use the absolute value of the t provided for the graph
  t <- abs(t)
  #Calculate the p value
  pvalue <- stats::pt(t, df = df, lower.tail = FALSE)*2
  #Label for half the p value (two tailed)
  phalflab <- as.character(as.expression(bquote(frac(p,2)~.(p_value_format(pvalue/2)))))
  #Label for p value (one tailed)
  plab <- as.character(as.expression(bquote(p~.(p_value_format(pvalue/2)))))
  #Labels for t value and - t value (two tailed)
  tlableft <- as.character(as.expression(bquote(- group("|",t,"|") == .(signif(-t, digits=signifdigitst)))))
  tlabright <- as.character(as.expression(bquote( + group("|",t,"|") == .(signif(t, digits=signifdigitst)))))
  #Label for t value (one tailed)
  tlab <- as.character(as.expression(bquote(t == .(signif(originalt, digits=signifdigitst)))))
  #Define x axis bounds as the maximum between t*3 or 3 (this avoids only the tip of the curve to be plotted when t is small, and keeps a nice t curve shape display)
  if (xmax == "auto") {
    xbound <- max(3*t, 2)
  } else {xbound <- xmax}
  #To ensure lines plotted by stat_function are smooth
  precisionfactor <-  5000
  #To define the function to plot in stat_function
  density <- function(x) stats::dt(x, df)
  #Use the maximum density (top of the curve) to scale the y axis
  maxdensity <- density(0)

  # Set the position of the p value automatically or manually
  if (length(p_value_position) == 1) {
    #Use the density corresponding to the t value to place the label above (if this density is too high places the label lower in order to avoid the label being out above the plot)
    y_plabel <- min(density(t)+maxdensity*.1, maxdensity*.7)
    #To place the p value labels on the x axis, at the middle of the part of the curve they correspond to
    x_plabel <- t+(xbound-t)/2
  } else {
    x_plabel <- abs(p_value_position[1])
    y_plabel <- abs(p_value_position[2])
  }

  #To place t labels on the x axis where their cutline is, avoiding that they overlap if the two cutlines are too close
  x_tlabel <- max(t, .5)
  #Define the fill color of the labels as white
  colorlabelfill <- "white"
  #Theme options
  if (theme == "default") {
    colormiddle <- colormiddle
    colorsides <- colorsides
    colorplabel <- colorplabel
  } else if (theme == "blackandwhite"){
    colormiddle <- "grey96"
    colorsides <- "darkgrey"
    colorplabel <- "black"
  } else if (theme == "whiteandred") {
    colormiddle <- "grey96"
    colorsides <- "firebrick3"
    colorplabel <- "firebrick3"
  } else if (theme == "blueandred") {
    colormiddle <- "#104E8B"
    colorsides <- "firebrick3"
    colorplabel <- "firebrick3"
  } else if (theme == "greenandred") {
    colormiddle <- "seagreen"
    colorsides <- "firebrick3"
    colorplabel <- "firebrick3"
  }else if (theme == "goldandblue") {
    colormiddle <- "#FFC61E"
    colorsides <- "#00337F"
    colorplabel <- "#00337F"
  }else warning("The ",'"', "theme", '"', " argument was not recognized. See documentation for a list of available color themes. Reverting to default.")
  #To make some colors transparent when `blank` parameter is TRUE (to only plot de probability density function in that case)
  if (blank == TRUE) {
    colorsides <- grDevices::adjustcolor("white", alpha.f = 0)
    colorcut <- grDevices::adjustcolor("white", alpha.f = 0)
    colorplabel <- grDevices::adjustcolor("white", alpha.f = 0)
    colorlabelfill <- grDevices::adjustcolor("white", alpha.f = 0)
  }
  else {
    #Do nothing
  }
  #Start if else statement relative to the two vs. one tail argument
  if (tails == "two") {
          ggplot2::ggplot(data.frame(x = c(-xbound*2, xbound*2)), ggplot2::aes(x)) +
          #Axis labels
          ggplot2::labs(x=xlabel,y=ylabel, size = 10) +
          #Middle area
          ggplot2::stat_function(fun = area_range(density, -xbound, xbound), geom="area", fill=colormiddle, n=precisionfactor) +
          #Right side area
          ggplot2::stat_function(fun = area_range(density, t, xbound), geom="area", fill=colorsides, n=precisionfactor) +
          #Left side area
          ggplot2::stat_function(fun = area_range(density, -xbound, -t), geom="area", fill=colorsides, n=precisionfactor) +
          #Left side curve
          ggplot2::stat_function(fun = density, xlim = c(-xbound,-t), colour = colorsidescurve,linewidth=curvelinesize,n=precisionfactor) +
          #Right side curve
          ggplot2::stat_function(fun = density, xlim = c(t,xbound), colour = colorsidescurve,linewidth=curvelinesize,n=precisionfactor) +
          #middle curve
          ggplot2::stat_function(fun = density, xlim = c(-t,t), colour = colormiddlecurve, n=precisionfactor, linewidth=curvelinesize) +
          #Define plotting area for extraspace below the graph to place t label
          ggplot2::coord_cartesian(xlim=c(-xbound,xbound),ylim=c(-.05, maxdensity)) +
          #Left cut line
          ggplot2::geom_vline(xintercept = -t, colour = colorcut, linewidth = cutlinesize) +
          #Right cut line
          ggplot2::geom_vline(xintercept = t, colour = colorcut, linewidth = cutlinesize) +
          #Left p label
          #ggplot2::geom_label(ggplot2::aes(-x_plabel,y_plabel,label = phalflab), parse = T, fill = colorlabelfill, colour=colorplabel, family = fontfamily) +
          ggplot2::annotate(geom = "label", x = -x_plabel, y = y_plabel, label = phalflab, parse = T, fill = colorlabelfill, colour=colorplabel, family = fontfamily) +
          #Right p label
          #ggplot2::geom_label(ggplot2::aes(x_plabel,y_plabel,label = phalflab), parse = T, fill = colorlabelfill, colour=colorplabel, family = fontfamily) +
          ggplot2::annotate(geom = "label", x = x_plabel, y = y_plabel, label = phalflab, parse = T, fill = colorlabelfill, colour=colorplabel, family = fontfamily) +
          #Left t label
          #ggplot2::geom_label(ggplot2::aes(-x_tlabel,-.04,label = tlableft), fill = colorlabelfill, colour=colorcut, parse = T, family=fontfamily) +
          ggplot2::annotate(geom = "label", x = -x_tlabel, y = -.04, label = tlableft, parse = T, fill = colorlabelfill, colour=colorcut, family=fontfamily) +
          #Right t label
          #ggplot2::geom_label(ggplot2::aes(x_tlabel,-.04,label = tlabright), parse = T, fill = colorlabelfill, colour=colorcut, family=fontfamily) +
          ggplot2::annotate(geom = "label", x = x_tlabel, y = -.04, label = tlabright, parse = T, fill = colorlabelfill, colour=colorcut, family=fontfamily) +
          #Add the title
          ggplot2::ggtitle(title) +
          #Apply black and white ggplot theme to avoid grey background, etc.
          ggplot2::theme_bw() +
          #Remove gridlines and pass fontfamily argument to ggplot2
          ggplot2::theme(
            panel.grid.major = ggplot2::element_blank(),
            panel.grid.minor = ggplot2::element_blank(),
            panel.border = ggplot2::element_blank(),
            axis.title = ggplot2::element_text(family = fontfamily),
            axis.text = ggplot2::element_text(family = fontfamily),
            axis.text.x = ggplot2::element_text(family = fontfamily),
            axis.text.y = ggplot2::element_text(family = fontfamily),
            plot.title = ggplot2::element_text(family = fontfamily, hjust = .5),
            legend.text = ggplot2::element_text(family = fontfamily),
            legend.title = ggplot2::element_text(family = fontfamily))



     #If "one" is passed as the tails argument, then we distinguish between a positive and a negative t value
  } else if (tails == "one") {
    if (originalt > 0) {
          ggplot2::ggplot(data.frame(x = c(-xbound*2, xbound*2)), ggplot2::aes(x)) +
          #Axis labels
          ggplot2::labs(x=xlabel,y=ylabel, size=10) +
          #Left and middle area
          ggplot2::stat_function(fun = area_range(density, -xbound, xbound), geom="area", fill=colormiddle, n=precisionfactor) +
          #Right side area
          ggplot2::stat_function(fun = area_range(density, t, xbound), geom="area", fill=colorsides, n=precisionfactor) +
          #Left and middle curve
          ggplot2::stat_function(fun = density, xlim = c(-xbound,t), colour = colormiddlecurve,linewidth = curvelinesize,n=precisionfactor) +
          #Right side curve
          ggplot2::stat_function(fun = density, xlim = c(t,xbound), colour = colorsidescurve,linewidth = curvelinesize,n=precisionfactor) +
          #Define plotting area for extraspace below the graph to place t label
          ggplot2::coord_cartesian(xlim=c(-xbound,xbound),ylim=c(-.05, maxdensity)) +
          #Right cut line
          ggplot2::geom_vline(xintercept = t, colour = colorcut, size = cutlinesize) +
          #Right p label
          #ggplot2::geom_label(ggplot2::aes(x_plabel,y_plabel,label = plab), parse = T, fill = colorlabelfill, colour=colorplabel,family = fontfamily) +
          ggplot2::annotate(geom = "label", x = x_plabel, y = y_plabel, label = plab, parse = T, fill = colorlabelfill, colour=colorplabel,family = fontfamily) +
          #Right t label
          #ggplot2::geom_label(ggplot2::aes(x_tlabel,-.05,label = tlab), parse = T, fill = colorlabelfill, colour=colorcut, family=fontfamily) +
          ggplot2::annotate(geom = "label", x = x_tlabel, y = -.05, label = tlab, parse = T, fill = colorlabelfill, colour=colorcut, family=fontfamily) +
          #Add the title
          ggplot2::ggtitle(title) +
          #Apply black and white ggplot theme to avoid grey background, etc.
          ggplot2::theme_bw() +
          #Remove gridlines and pass fontfamily argument to ggplot2
          ggplot2::theme(
            panel.grid.major = ggplot2::element_blank(),
            panel.grid.minor = ggplot2::element_blank(),
            panel.border = ggplot2::element_blank(),
            axis.title = ggplot2::element_text(family = fontfamily),
            axis.text = ggplot2::element_text(family = fontfamily),
            axis.text.x = ggplot2::element_text(family = fontfamily),
            axis.text.y = ggplot2::element_text(family = fontfamily),
            plot.title = ggplot2::element_text(family = fontfamily, hjust = .5),
            legend.text = ggplot2::element_text(family = fontfamily),
            legend.title = ggplot2::element_text(family = fontfamily))
    } else {
    #Plot left tailed t test plot
          ggplot2::ggplot(data.frame(x = c(-xbound*2, xbound*2)), ggplot2::aes(x)) +
          #Axis labels
          ggplot2::labs(x=xlabel,y=ylabel, size=10) +
          #Middle and right area
          ggplot2::stat_function(fun = area_range(density, -xbound, xbound), geom="area", fill=colormiddle, n=precisionfactor) +
          #Left side area
          ggplot2::stat_function(fun = area_range(density, -xbound, -t), geom="area", fill=colorsides, n=precisionfactor) +
          #Left side curve
          ggplot2::stat_function(fun = density, xlim = c(-xbound,-t), colour = colorsidescurve,linewidth = curvelinesize,n=precisionfactor) +
          #Middle and right curve
          ggplot2::stat_function(fun = density, xlim = c(-t,xbound), colour = colormiddlecurve, n=precisionfactor, linewidth = curvelinesize) +
          #Define plotting area for extraspace below the graph to place t label
          ggplot2::coord_cartesian(xlim=c(-xbound,xbound),ylim=c(-.05, maxdensity)) +
          #Left cut line
          ggplot2::geom_vline(xintercept = -t, colour = colorcut, size = cutlinesize) +
          #Left p label
          #ggplot2::geom_label(ggplot2::aes(-x_plabel,y_plabel,label = plab), parse = T, fill = colorlabelfill, colour=colorplabel, family = fontfamily) +
          ggplot2::annotate(geom = "label", x = -x_plabel, y = y_plabel, label = plab, parse = T, fill = colorlabelfill, colour=colorplabel, family = fontfamily) +
          #Left t label
          #ggplot2::geom_label(ggplot2::aes(-x_tlabel,-.05,label = tlab),fill = colorlabelfill, colour=colorcut, parse = T, family=fontfamily) +
          ggplot2::annotate(geom = "label", x = -x_tlabel, y = -.04, label = tlab, parse = T, fill = colorlabelfill, colour=colorcut, family=fontfamily) +
          #Add the title
          ggplot2::ggtitle(title) +
          #Apply black and white ggplot theme to avoid grey background, etc.
          ggplot2::theme_bw() +
          #Remove gridlines and pass fontfamily argument to ggplot2
          ggplot2::theme(
            panel.grid.major = ggplot2::element_blank(),
            panel.grid.minor = ggplot2::element_blank(),
            panel.border = ggplot2::element_blank(),
            axis.title = ggplot2::element_text(family = fontfamily),
            axis.text = ggplot2::element_text(family = fontfamily),
            axis.text.x = ggplot2::element_text(family = fontfamily),
            axis.text.y = ggplot2::element_text(family = fontfamily),
            plot.title = ggplot2::element_text(family = fontfamily, hjust = .5),
            legend.text = ggplot2::element_text(family = fontfamily),
            legend.title = ggplot2::element_text(family = fontfamily))
    }
  } else
    #Stop function is invalid tails argument is provided
    warning("Please specify the number of tails as ", '"', "two", '"', " or ", '"', "one", '"', " for the ", '"', "tails", '"', " argument. Reverting to default (two-tailed).", sep="")
   }