Functions to animate a multipanel plot of signal detection theory values in R

sdt_animation.R

# plot and animate SDT measures (assuming Gaussian densities)
# Copyright Tom Wallis, 2018.
# Shared under a CC-BY license.

library(ggplot2)
library(latex2exp)
library(animation)

# See bottom for a demo of animation

# Functions to generate plots ---------------------------------------------

line_alphas <- 0.4
distribution_plot <- function(mu_s, lambda, sigma_s=1, 
                              xlims=c(-5, 5)){
  my_colours <- ggthemes::fivethirtyeight_pal()(3)  # blue, red, green
  x <- seq(xlims[1], xlims[2], length.out = 100)
  f_n <- dnorm(x)
  f_s <- dnorm(x, mu_s, sigma_s)
  
  plt <- ggplot(data.frame(x = x, y = f_n), 
                aes(x = x, y = y)) +
    geom_line(colour = my_colours[2]) + 
    geom_line(data = data.frame(x = x, y = f_s), colour = my_colours[3]) + 
    geom_vline(xintercept = lambda, colour = my_colours[1], alpha = line_alphas) + 
    scale_x_continuous(name = TeX('Evidence $\\[\\x\\]$')) +
    scale_y_continuous(name = "", breaks = NULL) + 
    theme_classic()
  return(plt)
}

roc_plot_probability <- function(mu_s, lambda, sigma_s=1, 
                                 xlims=c(-5, 5),
                                 aucs=FALSE){
  x <- seq(xlims[1], xlims[2], length.out = 100)
  pfa <- 1 - pnorm(x)
  phit <- 1 - pnorm(x, mu_s, sigma_s)
  
  pfa_lambda <- 1 - pnorm(lambda)
  phit_lambda <- 1 - pnorm(lambda, mu_s, sigma_s)
  
  plt <- ggplot(data.frame(x = pfa, y = phit), 
                aes(x = x, y = y)) +
    geom_abline(slope = 1, linetype = "dashed", alpha = 0.5) + 
    geom_line(alpha = line_alphas) + 
    geom_point(data = data.frame(x = pfa_lambda, y = phit_lambda)) + 
    scale_x_continuous(name = TeX('\\mathrm{P_F}')) + 
    scale_y_continuous(name = TeX('\\mathrm{P_H}')) +
    theme_bw()
  
  if (aucs == TRUE){
    label <- paste0("AUC (numeric) = ", 
                    round(auc_numeric(phit, pfa), 3),
                    "\nA_z = ", 
                    round(auc_z(mu_s), 3))
    plt <- plt + 
      annotate("text", x = 0.7, y = 0.2,
               label=label)
  }
  return(plt)
}

roc_plot_gaussian <- function(mu_s, lambda, sigma_s=1, 
                              xlims=c(-5, 5)){
  x <- seq(xlims[1], xlims[2], length.out = 100)
  pfa <- 1 - pnorm(x)
  phit <- 1 - pnorm(x, mu_s, sigma_s)
  
  pfa_lambda <- 1 - pnorm(lambda)
  phit_lambda <- 1 - pnorm(lambda, mu_s, sigma_s)
  
  plt <- ggplot(data.frame(x = pfa, y = phit), 
                aes(x = qnorm(x), y = qnorm(y))) +
    geom_abline(slope = 1, linetype = "dashed", alpha = 0.5) + 
    geom_line(alpha = line_alphas) +
    geom_point(data = data.frame(x = pfa_lambda, y = phit_lambda)) + 
    scale_x_continuous(name = TeX('\\mathrm{z_F}')) + 
    scale_y_continuous(name = TeX('\\mathrm{z_H}')) +
    theme_bw()
  return(plt)
}

accuracy_plot <- function(mu_s, lambda, sigma_s=1, 
                          p_s = 0.5,
                          xlims=c(-5, 5)){
  x <- seq(xlims[1], xlims[2], length.out = 100)
  p_n <- 1 - p_s  # probability of a noise trial = 1 - p(signal)
  p_c <- p_s * (1 - pnorm(x, mu_s, sigma_s)) + p_n * pnorm(x)
  p_c_lambda <- p_s * (1 - pnorm(lambda, mu_s, sigma_s)) + p_n * pnorm(lambda)
  plt <- ggplot(data.frame(x = x, y = p_c), 
                aes(x = x, y = y)) +
    geom_line(alpha = line_alphas) + 
    geom_point(data = data.frame(x = lambda, y = p_c_lambda)) + 
    scale_x_continuous(name = TeX('Decision criterion $\\[\\lambda\\]$')) +
    scale_y_continuous(name = TeX('\\mathrm{P_C}'),
                       limits = c(0, 1),
                       breaks = seq(0, 1, by = 0.25)) + 
    annotate("text", x = 0, y = 0.2, 
             label = paste0("Proportion signal trials = ", p_s),
             size = 4) +
    theme_bw()
  return(plt)
}

bias_plot <- function(mu_s, lambda, sigma_s=1, xlims=c(-5, 5)){
  x <- seq(xlims[1], xlims[2], length.out = 100)
  log_beta_x <- log_beta(x, mu_s, sigma_s)
  log_beta_lambda <- log_beta(lambda, mu_s, sigma_s)
  plt <- ggplot(data.frame(x = x, y = log_beta_x), 
                aes(x = x, y = y)) +
    geom_line(alpha = line_alphas) + 
    geom_point(data = data.frame(x = lambda, y = log_beta_lambda)) + 
    geom_hline(yintercept = 0, linetype = "dashed") + 
    annotate("text", x = 0, y = 0.3, label = "'No' bias") + 
    annotate("text", x = 0, y = -0.3, label = "'Yes' bias") + 
    scale_x_continuous(name = TeX('Decision criterion $\\[\\lambda\\]$')) +    
    scale_y_continuous(name = TeX('Response bias \\[$\\log \\beta$\\]')) + 
    theme_bw()
  return(plt)  
}

sensitivity_plot <- function(mu_s, sigma_s=1, xlims=c(-5, 5)){
  # create a vector of mu_s for plot:
  mu_vec <- seq(0, 4, length.out = 100)
  if (sigma_s == 1){
    # equal variance gaussian model
    auc_mu_s <- auc_z(mu_s)
    auc_vec <- auc_z(mu_vec)
  } else {
    x <- seq(xlims[1], xlims[2], length.out = 100)
    calc_auc_vec <- function(mu, x, sigma_s){
      pfa <- 1 - pnorm(x)
      phit <- 1 - pnorm(x, mu, sigma_s)
      auc <- auc_numeric(phit, pfa)
      return(auc)
    }
    auc_mu_s <- calc_auc_vec(mu_s, x, sigma_s)
    auc_vec <- sapply(mu_vec, calc_auc_vec, x=x, sigma_s=sigma_s)
  }
  
  plt <- ggplot(data.frame(x = mu_vec, y = auc_vec), 
                aes(x = x, y = y)) +
    geom_line(alpha = line_alphas) + 
    geom_point(data = data.frame(x = mu_s, y = auc_mu_s)) + 
    scale_x_continuous(name = TeX('$\\mu_s$')) +
    scale_y_continuous(name = "AUC") + 
    theme_bw()
  return(plt)  
}

summary_plot <- function(mu_s, lambda, sigma_s=1, 
                         p_s = 0.5,
                         xlims=c(-5, 5)){
  
  dists <- distribution_plot(mu_s = mu_s, lambda = lambda, 
                             sigma_s = sigma_s, 
                             xlims=xlims)
  
  roc_1 <- roc_plot_probability(mu_s = mu_s, lambda = lambda, 
                                sigma_s = sigma_s, 
                                xlims=xlims)
  
  roc_2 <- roc_plot_gaussian(mu_s = mu_s, lambda = lambda, 
                             sigma_s = sigma_s, 
                             xlims=xlims)
  
  accuracy <- accuracy_plot(mu_s = mu_s, lambda = lambda, 
                            sigma_s = sigma_s, 
                            p_s = p_s,
                            xlims=xlims)  
  
  bias <- bias_plot(mu_s = mu_s, lambda = lambda, 
                    sigma_s = sigma_s, 
                    xlims=xlims)
  
  sensitivity <- sensitivity_plot(mu_s = mu_s,
                                  sigma_s = sigma_s, 
                                  xlims=xlims)
  
  plt <- gridExtra::grid.arrange(dists, roc_1, roc_2,
                                 sensitivity, bias, accuracy, 
                                 ncol = 3)
  return(plt)  
}

log_beta <- function(x, mu_s, sigma_s){
  return(0.5 * (x^2 - ((x - mu_s)^2) / sigma_s^2) - log(sigma_s))
}


auc_numeric <- function(hr, far){
  # a simple numerical computation of AUC.
  # adapted from http://blog.revolutionanalytics.com/2016/11/calculating-auc.html
  # check that inputs are sorted inputs:
  hr <- sort(hr)
  far <- sort(far)
  
  dfar <- c(diff(far), 0)
  dhr <- c(diff(hr), 0)
  return(sum(hr * dfar) + sum(dhr * dfar)/2)
}

auc_z <- function(mu_s){
  # assumes equal-variance Gaussian model
  return(pnorm(mu_s / sqrt(2)))
}


# Demo animation ----------------------------------------------------------

# demo with equal variance Gaussian model; dprime = 1
mu_s <- 1
sigma_s <- 1
p_s <- 0.5

# do animation:
x <- seq(-2, 4, length.out = 50)  # you could change this to animate another parameter
saveGIF({
  for (i in 1:length(x)){
    plt <- summary_plot(mu_s = mu_s, 
                        lambda = x[i], 
                        sigma_s = sigma_s,
                        p_s = p_s,
                        xlims = c(-2, 4))
  }
}, movie.name = "sdt_animation.gif", interval = 0.1, nmax = length(x), 
ani.width = 800, ani.height = 400)