R/plotting.R

Defines functions make_pelogit_fnc plot_avg_cdiff plot_avg_diff plot_avg_contour plot_avg .pac_theme

Documented in make_pelogit_fnc plot_avg plot_avg_cdiff plot_avg_contour plot_avg_diff

## SHARED CODE BEGINS HERE ##
.pac_theme <- function(base_size = 12, base_family = ""){
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.minor.y = element_blank(),
          plot.title = element_text(hjust = 0.5, vjust = 1)
    )
}
## SHARED CODE ENDS HERE ##


#' Plots average looks to interest areas.
#' 
#' \code{plot_avg} calculates the grand or conditional averages of 
#' looks to each interest area along with standard error. It then plots the results.
#' N.B.: This function will work for data with a maximum of 8 interest areas
#' and 2 conditions.
#' 
#' @export
#' @import ggplot2
#' @import dplyr
#' @import rlang
#' 
#' @param data A data table object output by either \code{\link{bin_prop}}. 
#' \code{\link{transform_to_elogit}}, or \code{\link{create_binomial}}.
#' @param type A character string indicating "proportion" or "elogit" which 
#' influences how standard error and confidence intervals are calculated.
#' @param xlim A vector of two integers specifying the limits of the x-axis.
#' @param IAColumns A named character vector specifying the desired interest 
#' area columns with custom strings for the legend.
#' @param Averaging A character string indicating how the averaging should 
#' be done. "Event" (default) will produce the overall mean in the data, while
#' "Subject" or "Item" (or, in principle, any other column name) will
#' calculate the grand mean by that factor.
#' @param Condition1 A string containing the column name corresponding to the 
#' first condition, if available. 
#' @param Condition2 A string containing the column name corresponding to the 
#' second condition, if available.
#' @param Cond1Labels A named character vector specifying the desired custom 
#' labels of the levels of the first condition. 
#' @param Cond2Labels A named character vector specifying the desired custom 
#' labels of the levels of the second condition. 
#' @param ErrorBar A logical indicating whether error bars should be
#' included in the plot.
#' @param ErrorBand A logical indicating whether error bands should be
#' included in the plot.
#' @param ErrorType A string indicating "SE" or "CI". For SE, the calculation 
#' varies for empirical logits and proportions. Further, for CI, the calculation
#' on proportions uses the Wald method. 
#' @param ConfLev A number indicating the confidence level of the CI.
#' @param CItype A string indicating "simultaneous" or "pointwise". Simultaneous
#' performs a Bonferroni correction for the interval.
#' @param VWPreTheme A logical indicating whether the theme included with the 
#' function should be applied, or ggplot2's base theme (to which any other 
#' custom theme could be added).
#' @examples
#' \dontrun{
#' library(VWPre)
#' # For plotting the grand average with the included theme and SE bars
#' plot_avg(data = dat, type = "elogit", xlim = c(0, 1000), 
#'    IAColumns = c(IA_1_ELogit = "Target", IA_2_ELogit = "Rhyme", 
#'    IA_3_ELogit = "OnsetComp", IA_4_ELogit = "Distractor"), 
#'    Averaging = "Event", Condition1 = NA, Condition2 = NA, 
#'    Cond1Labels = NA, Cond2Labels = NA,
#'    ErrorBar = TRUE, VWPreTheme = TRUE, ErrorType = "SE",
#'    ErrorBand = FALSE) 
#'    
#' # For plotting conditional averages (one condition) with the included theme
#' # and 95% simultaneous CI bars.
#' # This produces plots arranged horizontally
#' plot_avg(data = dat, type = "elogit", xlim = c(0, 1000), 
#'    IAColumns = c(IA_1_ELogit = "Target", IA_2_ELogit = "Rhyme", 
#'    IA_3_ELogit = "OnsetComp", IA_4_ELogit = "Distractor"),
#'    Averaging = "Event", Condition1 = NA, Condition2 = "talker", 
#'    Cond1Labels = NA, 
#'    Cond2Labels = c(CH1 = "Chinese 1", CH10 = "Chinese 3", CH9 = "Chinese 2", 
#'    EN3 = "English 1"), ErrorBar = TRUE, VWPreTheme = TRUE,
#'    ErrorBand = FALSE, ErrorType = "CI", ConfLev = 95, CItype = "simultaneous")
#' 
#' # For plotting conditional averages (two conditions) for one interest area
#' with the included theme and 95% simultaneous CI bands.
#' # This produces plots arranged in grid format.
#' plot_avg(data = dat, type = "elogit", xlim = c(0, 1000),
#'    IAColumns = c(IA_1_ELogit = "Target"), Averaging = "Event", 
#'    Condition1 = "talker", Condition2 = "Exp",
#'    Cond1Labels = c(CH1 = "Chinese 1", CH10 = "Chinese 3", CH9 = "Chinese 2", 
#'    EN3 = "English 1"), Cond2Labels = c(High = "H Exp", Low = "L Exp"),
#'    ErrorBar = FALSE, VWPreTheme = TRUE, ErrorBand = TRUE, 
#'    ErrorType = "CI", ConfLev = 95, CItype = "simultaneous")
#' 
#' #' # For a more complete tutorial on VWPre plotting functions:
#' vignette("SR_Plotting", package="VWPre")
#' }
#' 
plot_avg <- function(data, type = NULL, xlim = NA, IAColumns = NULL, 
                     Averaging = "Event", Condition1 = NULL, 
                     Condition2 = NULL, Cond1Labels = NA, 
                     Cond2Labels = NA, ErrorBar = TRUE, VWPreTheme = TRUE,
                     ConfLev = 95, CItype = "simultaneous",
                     ErrorBand = FALSE, ErrorType = "SE") {
  
 
	.check_for_PupilPre(type = "UseOther", suggest = "ppl_plot_avg")
	
  if(is.null(type)){
    stop("Please supply the plot type!")
  }
  
  # Set plot type and y-axis
  if (type == "proportion") {
    ylabel = "Proportion Looks"
    ylim = c(0,1)
  } else if (type == "elogit") {
    ylabel = "Empirical Logit Looks"
    ylim = c(-4,4)
  } else {
    stop("You must specify 'proportion' or 'elogit'.")
  }
  
  Columns <- IAColumns
  
  if(is.null(Columns)){
    stop("Please supply the column(s) to be plotted!")
  } else if (length(names(Columns))>8) {
    stop("You have more than 8 interest areas; you must modify this function.")
  } else {
    missing <- names(Columns)[which(!(names(Columns) %in% colnames(data)))]
    if(length(missing)>0) {
      stop(paste(utils::capture.output(print(missing)), "column(s) not present in data!"))
    }
  }
  
  Theme <- VWPreTheme
  
  {## SHARED CODE BEGINS HERE ## 
    
    if(!(Averaging %in% colnames(data))){
      stop(paste(Averaging, " column not present in data!"))
    }
    
    if(ErrorBar==TRUE && ErrorBand==TRUE){
      stop(paste("Please select either error bars OR error bands. 
                 For error bars set ErrorBar=TRUE and ErrorBand=FALSE.
                 For error bands set ErrorBar=FALSE and ErrorBand=TRUE."))
    }
    
    # Set x-axis
    if (is.na(xlim[1])) {
      xlim <- c(range(data$Time)[1], range(data$Time)[2])
      xaxis <- unique(data$Time)
    } else {
      xlim <- xlim
      data<-data[data$Time>=xlim[1] & data$Time<=xlim[2],]
      xaxis <- unique(data$Time)
    }
    
    # Check condition columns
    if(!(is.null(Condition1))) {
      if(is.na(Condition1)) {
        stop(paste("Please use NULL when not specifying a condition."))
      } else if (!(Condition1 %in% colnames(data))) {
        stop(paste(Condition1, " column not present in data!"))
      }
    }
    if(!(is.null(Condition2))) {
      if(is.na(Condition2)) {
        stop(paste("Please use NULL when not specifying a condition."))
      } else if (!(Condition2 %in% colnames(data))) {
        stop(paste(Condition2, " column not present in data!"))
      }
    }
    
    # Inform about Grand Average calculation
    message(paste0("Grand average calculated using ", Averaging, " means."))
    
    # Set columns for selection
    if(is.null(Condition1) && is.null(Condition2)){
      sel_names = quos(UQ(sym(Averaging)), Time, names(Columns))
    } else if(!is.null(Condition1) && is.null(Condition2)){
      sel_names = quos(UQ(sym(Averaging)), Time, names(Columns), UQ(sym(Condition1)))
    } else if(is.null(Condition1) && !is.null(Condition2)){
      sel_names = quos(UQ(sym(Averaging)), Time, names(Columns), UQ(sym(Condition2)))
    } else if(!is.null(Condition1) && !is.null(Condition2)){
      sel_names = quos(UQ(sym(Averaging)), Time, names(Columns), UQ(sym(Condition1)), UQ(sym(Condition2)))
    }
    
    # Set columns for gathering
    gath_col = quos(names(Columns))
    
    # Set columns for grouping
    if(is.null(Condition1) && is.null(Condition2)){
      group_names1 = quos(UQ(sym(Averaging)), IA, Time)
      group_names2 = quos(IA, Time)
    } else if(!is.null(Condition1) && is.null(Condition2)){
      group_names1 = quos(UQ(sym(Averaging)), IA, Time, UQ(sym(Condition1)))
      group_names2 = quos(IA, Time, UQ(sym(Condition1)))
    } else if(is.null(Condition1) && !is.null(Condition2)){
      group_names1 = quos(UQ(sym(Averaging)), IA, Time, UQ(sym(Condition2)))
      group_names2 = quos(IA, Time, UQ(sym(Condition2)))
    } else if(!is.null(Condition1) && !is.null(Condition2)){
      group_names1 = quos(UQ(sym(Averaging)), IA, Time, UQ(sym(Condition1)), UQ(sym(Condition2)))
      group_names2 = quos(IA, Time, UQ(sym(Condition1)), UQ(sym(Condition2)))
    }
    
    # # Make labeller for multipanel
    # if (length(Columns) > 1) {
    #   if (!is.null(Condition1) || !is.null(Condition2)) {
    #   if(is.na(Cond1Labels) && is.na(Cond2Labels)) {
    #     my_labeller <- "label_value"
    #   } else if (is.na(Cond1Labels) && !is.na(Cond2Labels)) {
    #     my_labeller <- labeller(
    #       CustCond1 = Cond1Labels
    #       )
    #   } else if (!is.na(Cond1Labels) && is.na(Cond2Labels)) {
    #     my_labeller <- labeller(
    #       CustCond2 = Cond2Labels
    #     )
    #   }
    #     else {
    #   my_labeller <- labeller(
    #     CustCond1 = Cond1Labels,
    #     CustCond2 = Cond2Labels
    #   )
    #   }
    #   }
    # }
    my_labeller <- labeller()
    
    # Prepare for averaging
    Avg <- data %>% select(!!!sel_names) %>% 
      tidyr::gather("IA", "VALUE", !!!gath_col, na.rm = FALSE, convert = FALSE) %>%
      group_by(!!!group_names1) %>% 
      summarise(VALUE = mean(VALUE, na.rm = TRUE)) %>% 
      group_by(!!!group_names2)
    
    # Execute calculation
    if(type == "elogit") {
      Avg <- Avg %>% summarise(mean = mean(VALUE, na.rm = TRUE), n=n(), se = stats::sd(VALUE, na.rm = TRUE) / sqrt(n()))
      if(CItype=="pointwise") {
        tval <- 1-(((100-ConfLev)/2)/100)
      } else if (CItype=="simultaneous") {
        tval <- 1-(((100-ConfLev)/(2*length(unique(xaxis))))/100)
      }
      Avg <- Avg %>% mutate(ci = stats::qt(tval,df=n-1)*se)
    } else if (type=="proportion") {
      Avg <- Avg %>% summarise(mean = mean(VALUE, na.rm = TRUE), n=n(), se = sqrt((mean(VALUE, na.rm = TRUE)*(1-mean(VALUE, na.rm = TRUE)))/n()))
      if(CItype=="pointwise") {
        zval <- 1-(((100-ConfLev)/2)/100)
      } else if (CItype=="simultaneous") {
        zval <- 1-(((100-ConfLev)/(2*length(unique(xaxis))))/100)
      }
      Avg <- Avg %>% mutate(ci = stats::qnorm(zval)*se)
    }
    Avg <- ungroup(Avg)
    
    # Prepare Condition columns and faceting or legend title/labels
    if (!is.null(Condition1) && is.null(Condition2)) {
      Avg <- Avg %>%
        rename(CustCond1 = UQ(sym(Condition1)))
      if (any(!is.na(Cond1Labels))) {
        lev1 <- unique(levels(Avg$CustCond1))
        for (x in 1:length(names(Cond1Labels))) {
          for(i in 1:length(lev1)) {
            if (lev1[i] == names(Cond1Labels)[x]) {
              lev1[i] <- Cond1Labels[[x]]
            }
          }
        }
        levels(Avg$CustCond1) <- lev1
      }
      Avg <- Avg %>%
        mutate(Cond = CustCond1)
      if (length(Columns)>1) {
        my_facet <- function() {
          eval(
            facet_grid(CustCond1 ~ ., labeller = my_labeller)
          )
        }
      } else {
        Cond <- Condition1
      }
    } else if (is.null(Condition1) && !is.null(Condition2)) {
      Avg <- Avg %>%
        rename(CustCond2 =  UQ(sym(Condition2))) 
      if (any(!is.na(Cond2Labels))) {
        lev2 <- unique(levels(Avg$CustCond2))
        for (x in 1:length(names(Cond2Labels))) {
          for(i in 1:length(lev2)) {
            if (lev2[i] == names(Cond2Labels)[x]) {
              lev2[i] <- Cond2Labels[[x]]
            }
          }
        }
        levels(Avg$CustCond2) <- lev2
      }
      Avg <- Avg %>%
        mutate(Cond = CustCond2)
      if (length(Columns)>1) {
        my_facet <- function() {
          eval(
            facet_grid(. ~ CustCond2, labeller = my_labeller)
          )
        }
      } else {
        Cond <- Condition2
      }
    } else if (!is.null(Condition1) && !is.null(Condition2)) {
      Avg <- Avg %>%
        rename(CustCond1 =  UQ(sym(Condition1)), CustCond2 =  UQ(sym(Condition2))) 
      if (any(!is.na(Cond1Labels))) {
        lev1 <- unique(levels(Avg$CustCond1))
        for (x in 1:length(names(Cond1Labels))) {
          for(i in 1:length(lev1)) {
            if (lev1[i] == names(Cond1Labels)[x]) {
              lev1[i] <- Cond1Labels[[x]]
            }
          }
        }
        levels(Avg$CustCond1) <- lev1
      }
      if (any(!is.na(Cond2Labels))) {
        lev2 <- unique(levels(Avg$CustCond2))
        for (x in 1:length(names(Cond2Labels))) {
          for(i in 1:length(lev2)) {
            if (lev2[i] == names(Cond2Labels)[x]) {
              lev2[i] <- Cond2Labels[[x]]
            }
          }
        }
        levels(Avg$CustCond2) <- lev2
      }
      Avg <- Avg %>%
        mutate(Cond = paste(CustCond1, CustCond2, sep = "_"))
      if (length(Columns)>1) {
        my_facet <- function() {
          eval(
            facet_grid(CustCond1 ~ CustCond2, labeller = my_labeller)
          )
        }
      } else {
        Cond <- paste(Condition1, Condition2, sep = "_by_")
      }
    }
    
    # Setting Error
    if(ErrorType=="SE"){
      Avg$error_lower <- Avg$mean - Avg$se
      Avg$error_upper <- Avg$mean + Avg$se
    } else if(ErrorType=="CI") {
      if(type=="elogit") {
        Avg$error_lower <- Avg$mean - Avg$ci
        Avg$error_upper <- Avg$mean + Avg$ci
      } else if(type=="proportion") {
        Avg$error_lower <- Avg$mean - Avg$ci
        Avg$error_lower <- ifelse(Avg$error_lower < 0, 0, Avg$error_lower)
        Avg$error_upper <- Avg$mean + Avg$ci
        Avg$error_upper <- ifelse(Avg$error_upper > 1, 1, Avg$error_upper)
      }
    }
    
    # Setting appropriate ylim
    if (type == "elogit") {
      ylim[1] = min(Avg$error_lower)
      ylim[2] = max(Avg$error_upper)
    } else if (type == "proportion") {
      ylim[1] = min(Avg$error_lower)
      ylim[2] = max(Avg$error_upper)
      if (ylim[1] > 0) {
        ylim[1] = 0
      }
      if (ylim[2] < 1) {
        ylim[2] = 1
      }
    }
    
    # Print message regarding CIs
    if (ErrorType=="CI" && (ErrorBand==TRUE | ErrorBar==TRUE)) {
      if(CItype == "pointwise") {
        message(paste0("Plot created with ", CItype, " confindence intervals (set to ", ConfLev, "%)."))
      }
      if(CItype == "simultaneous") {
        if(type=="elogit") {
          message(paste0("Plot created with ", CItype, " confindence intervals (adjusted to ", round(tval*100, 2), "% using the Bonferroni method)."))
        } else if(type=="proportion") {
          message(paste0("Plot created with ", CItype, " confindence intervals (adjusted to ", round(zval*100, 2), "% using the Bonferroni method)."))
        }
      }
    }
    
    # Set plot grouping
    if (length(Columns) > 1) {
      Col <- "IA"
    } else {
      if (!is.null(Condition1) || !is.null(Condition2)){
        Col <- "Cond"
      } else {
        Col <- NULL
      }
    }
    
    # Basic plot
    plt <- ggplot(Avg, aes_string(x = "Time", y = "mean", colour = Col)) +
      ylab(ylabel) +
      scale_x_continuous(limits = c(xlim[1], xlim[2])) +
      scale_y_continuous(limits = c(ylim[1], ylim[2]))
    
    # Basic themeing
    if (Theme == TRUE) {
      plt <- plt + .pac_theme()
    }
    
    # Determine/add faceting
    if (!is.null(Condition1) || !is.null(Condition2)) {
      if (length(Columns) > 1) {
        plt <- plt + my_facet()
      }
    }
    
    # Set errorbands according to theme and plot grouping
    if(ErrorBand==TRUE & Theme==TRUE) {
      plt <- plt + aes_string(shape = Col, colour = NULL) +
        geom_ribbon(aes_string(ymin="error_lower", ymax="error_upper", linetype=NA), fill="grey25", alpha=0.15, show.legend = FALSE) 
      if(!is.null(Col)) {
        if(Col=="IA") {
          plt <- plt + 
            scale_shape_discrete(name="Interest Area",
                                 breaks=names(Columns),
                                 labels=Columns)
        } else {
          plt <- plt + 
            scale_shape_discrete(name=Cond,
                                 breaks=unique(Avg$Cond),
                                 labels=unique(Avg$Cond))
        }
      } else {
        plt <- plt
      }
    }
    if(ErrorBand==TRUE & Theme==FALSE) {
      plt <- plt + 
        geom_ribbon(aes_string(ymin="error_lower", ymax="error_upper", linetype=NA, fill=Col), alpha=0.15, show.legend = FALSE) 
      if(!is.null(Col)) {
        if(Col=="IA") {
          plt <- plt + 
            scale_colour_hue(name="Interest Area",
                             breaks=names(Columns),
                             labels=Columns)
        } else {
          plt <- plt + 
            scale_colour_hue(name=Cond,
                             breaks=unique(Avg$Cond),
                             labels=unique(Avg$Cond))
        }
      } else {
        plt <- plt
      }
    }
    
    # Set errorbars according to theme and plot grouping
    if(ErrorBar==TRUE & Theme==TRUE) {
      plt <- plt + 
        geom_errorbar(aes_string(ymin="error_lower", ymax="error_upper"), width = .3) 
      if(!is.null(Col)) {
        if(Col=="IA") {
          plt <- plt + 
            scale_colour_grey(name="Interest Area",
                              breaks=names(Columns),
                              labels=Columns)
        } else {
          plt <- plt + 
            scale_colour_grey(name=Cond,
                              breaks=unique(Avg$Cond),
                              labels=unique(Avg$Cond))
        }
      } else {
        plt <- plt
      }
    }
    if(ErrorBar==TRUE & Theme==FALSE) {
      plt <- plt +
        geom_errorbar(aes_string(ymin="error_lower", ymax="error_upper", color=Col), width = .3)
      if(!is.null(Col)) {
        if(Col=="IA") {
          plt <- plt + 
            scale_colour_hue(name="Interest Area",
                             breaks=names(Columns),
                             labels=Columns)
        } else {
          plt <- plt + 
            scale_colour_hue(name=Cond,
                             breaks=unique(Avg$Cond),
                             labels=unique(Avg$Cond))
        }
      } else {
        plt <- plt
      }
    }
    
    # Finish plot
    plt <- plt + geom_point() + geom_line() 

  plt
  
  } ## SHARED CODE ENDS HERE ##
  
}
  

  #' Plots average contour surface of looks to a given interest area.
#' 
#' \code{plot_avg_contour} calculates the conditional average of proportions or 
#' empirical logit looks to a given interest area by Time and a specified 
#' continuous variable. It then applies a 3D smooth 
#' (derived using \code{\link[mgcv]{gam}}) over 
#' the surface and plots the results as a contour plot. 
#' 
#' @export
#' @import ggplot2
#' @import dplyr
#' @import rlang
#'  
#' @param data A data table object output by either \code{\link{bin_prop}}. 
#' \code{\link{transform_to_elogit}}, or \code{\link{create_binomial}}.
#' @param IA A string specifying the column name of the IA to use. 
#' @param type A character string indicating "proportion" or "elogit".
#' @param xlim A vector of two integers specifying the limits of the x-axis.
#' @param Var A string containing the column name corresponding to the continuous
#' variable.
#' @param Averaging A character string indicating how the averaging should 
#' be done. "Event" (default) will produce the overall mean in the data, while
#' "Subject" or "Item" (or, in principle, any other column name) will
#' calculate the grand mean by that factor.
#' @param VarLabel A string specifying the axis label to use for \code{Var}.
#' @param VWPreTheme A logical indicating whether the theme included with the 
#' function, or ggplot2's base theme (which any other custom theme could be added).
#' @param Colors A vector of two strings specifying the colrs of the contour shading - The default values represent grayscale.
#' @examples
#' \dontrun{
#' library(VWPre)
#' # For plotting a conditional contour surface...
#' plot_avg_contour(data = dat, IA = "IA_1_ELogit", type = "elogit", 
#'                Var = "Rating", VarLabel = "Accent Rating", xlim = c(0,1000), 
#'                VWPreTheme = FALSE, Colors = c("red", "white"))
#' 
#' # For a more complete tutorial on VWPre plotting functions:
#' vignette("SR_Plotting", package="VWPre")
#' }
plot_avg_contour <- function(data, IA = NULL, type = NULL, Var = NULL, 
                             Averaging = "Event",
                             VarLabel = NULL, xlim=NA, VWPreTheme=TRUE, 
                             Colors=c("gray20", "gray90")) {
  
  Column <- IA
  Theme <- VWPreTheme
  
  if(is.null(type)){
    stop("Please supply the plot type!")
  } else {
    if(type=="proportion"){
      LegendName <- "Mean\nproportion\n"
    } else {
      LegendName <- "Mean\nempirical\nlogit\n"
    } 
  }
  
  {## SHARED CODE BEGINS HERE ## 
    
    if(is.null(Column)){
      stop("Please supply the column to be plotted!")
    } else if(!(Column %in% colnames(data))) {
      stop(paste(Column, "column not present in data!"))
    }

    if(is.null(Var)){
      stop("Please supply the variable column name to be used in the contour plot!")
    } else if(!(Var %in% colnames(data))) {
      stop(paste(Var, "column not present in data!"))
    }
    
    if(is.numeric(eval(parse(text=paste0("data$", Var)))) | is.integer(eval(parse(text=paste0("data$", Var))))){
    } else {
      stop("The contour plot variable must be of class 'numeric' or class 'integer'!")
    }
    
    if(!(Averaging %in% colnames(data))){
      stop(paste(Averaging, " column not present in data!"))
    } else (
      message(paste0("Surface calculated using ", Averaging, " means."))
    )
    
    # # Check condition columns
    # if(!(is.null(Condition1))) {
    #   if(is.na(Condition1)) {
    #     stop(paste("Please use NULL when not specifying a condition."))
    #   } else if (!(Condition1 %in% colnames(data))) {
    #     stop(paste(Condition1, " column not present in data!"))
    #   }
    # }
    # if(!(is.null(Condition2))) {
    #   if(is.na(Condition2)) {
    #     stop(paste("Please use NULL when not specifying a condition."))
    #   } else if (!(Condition2 %in% colnames(data))) {
    #     stop(paste(Condition2, " column not present in data!"))
    #   }
    # }
    
    zlim <- c(-4,4)
    Colors <- Colors
    
    sel_names <- quos(UQ(sym(Averaging)), UQ(sym(Column)), Time, UQ(sym(Var)))
    group_names1 = quos(UQ(sym(Averaging)), UQ(sym(Column)), UQ(sym(Var)), Time)
    group_names2 = quos(UQ(sym(Column)), UQ(sym(Var)), Time)
    
    
    if(is.null(VarLabel)) {
      VarLabel <- Var
    }
    
    if (is.na(xlim)[1]) {
      xlim <- c(range(data$Time)[1], range(data$Time)[2])
    } else {
      xlim <- xlim
    }
    
    Avg <- data %>% select(!!!sel_names) %>% 
      group_by(!!!group_names1) %>%
      summarise(Columnmean = mean(UQ(sym(Column)), na.rm = TRUE)) %>% 
      group_by(!!!group_names2) %>% 
      summarise(mean = mean(Columnmean, na.rm = TRUE))
    
    if (type == "proportion") {
      Avg$meanon <- round(Avg$mean*100)
      Avg$meanoff <- 100 - Avg$meanon
      Avg$meanbinom <- cbind(Avg$meanon, Avg$meanoff)
      var<-list(Var)
      model <- lapply(var, function(x) {
        mgcv::gam(substitute(meanbinom ~ te(Time, i), list(i = as.name(x))), data = Avg, family = binomial)
      })
    } else {
      var<-list(Var)
      model <- lapply(var, function(x) {
        mgcv::gam(substitute(mean ~ te(Time, i), list(i = as.name(x))), data = Avg)
      })
    }
    
    model<-model[[1]]
    
    names(Avg)[names(Avg)==Var] <- "Variable"
    Varcol<-Avg$Variable
    Varcol<-unique(Varcol)
    
    # create a vector for variable A of length np
    data<-data[order(data$Event, data$Time),]
    rate <- 1000 / (data$Time[2] - data$Time[1])
    nptime<-round((xlim[2]-xlim[1])/100*(100/(1000/rate)))
    time<-seq(xlim[1],xlim[2],length.out=nptime) 
    # create a vector for variable B of length np
    npvar<-length(unique(Varcol))
    variable<-seq(range(Varcol)[1],range(Varcol)[2],length.out=npvar)
    tmp0<-expand.grid(time,variable)
    colnames(tmp0)<-c("Time", Var)
    newdat<-tmp0
    
    # compute lpmatrix
    X1<-stats::predict(model,newdat,type="lpmatrix")
    mean<-X1%*%stats::coef(model)
    predsmooth<-cbind(newdat, mean)
    names(predsmooth)[names(predsmooth)==Var] <- "Variable"
    
    if (type == "proportion") {
      predsmooth$mean <- stats::plogis(predsmooth$mean)
      zlim[1] = 0
      zlim[2] = 1
    } else {
      zlim[1] = min(predsmooth$mean) - 0.25
      zlim[2] = max(predsmooth$mean) + 0.25
    }
    
    plt <- ggplot(predsmooth) + 
      aes(x = Time, y = Variable, z = mean, fill = mean) + 
      geom_tile(aes(fill = mean)) +
      stat_contour() +
      geom_contour(color = "white") + 
      scale_fill_gradient(name = LegendName, limit=c(zlim[1], zlim[2]), low = Colors[1], high = Colors[2]) + 
      scale_x_continuous(expand=c(0,0)) +
      scale_y_continuous(expand=c(0,0)) +
      ylab(VarLabel) 
    if (Theme == TRUE) {
      plt <- plt +
        .pac_theme()
    } 
    
  plt
    
  } ## SHARED CODE ENDS HERE ##

}


#' Plots average difference between looks to two interest areas.
#' 
#' \code{plot_avg_diff} calculates the grand or conditional averages of 
#' differences between looks to two interest area along with standard error. 
#' It then plots the results.
#' 
#' @export
#' @import ggplot2
#' @import dplyr
#' @import rlang
#' 
#' @param data A data table object output by either \code{\link{bin_prop}}. 
#' \code{\link{transform_to_elogit}}, or \code{\link{create_binomial}}.
#' @param type A character string indicating "proportion" or "elogit" which 
#' influences how standard error and confidence intervals are calculated.
#' @param xlim A vector of two integers specifying the limits of the x-axis.
#' @param DiffCols A named character vector specifying the desired columns 
#' corresponding to the interest areas. 
#' @param Averaging A character string indicating how the averaging should 
#' be done. "Event" (default) will produce the overall mean in the data, while
#' "Subject" or "Item" (or, in principle, any other column name) will
#' calculate the grand mean by that factor.
#' @param Condition1 A string containing the column name corresponding to the 
#' first condition, if available. 
#' @param Condition2 A string containing the column name corresponding to the 
#' second condition, if available.
#' @param Cond1Labels A named character vector specifying the desired labels 
#' of the levels of the first condition. 
#' @param Cond2Labels A named character vector specifying the desired labels 
#' of the levels of the second condition. 
#' @param ErrorBar A logical indicating whether error bars should be
#' included in the plot.
#' @param ErrorBand A logical indicating whether error bands should be
#' included in the plot.
#' @param ErrorType A string indicating "SE" or "CI".
#' @param ConfLev A number indicating the confidence level of the CI.
#' @param CItype A string indicating "simultaneous" or "pointwise". Simultaneous
#' performs a Bonferroni correction for the interval.
#' @param VWPreTheme A logical indicating whether the theme included with the 
#' function should be applied, or ggplot2's base theme (to which any other 
#' custom theme could be added).
#' @examples
#' \dontrun{
#' library(VWPre)
#' # For plotting average differences with SE bars...
#' plot_avg_diff(data = dat, xlim = c(0, 1000), type = "proportion",
#'              DiffCols = c(IA_1_P = "Target", IA_2_P = "Rhyme"),
#'              Condition1 = NA, Condition2 = NA, Cond1Labels = NA, Cond2Labels = NA,
#'              ErrorBar = TRUE, VWPreTheme = TRUE, ErrorBand = FALSE, 
#'              ErrorType = "SE")
#'              
#' # For plotting conditional average differences (one condition) with the 
#' # included theme and 95% pointwise CI bars.
#' plot_avg_diff(data = dat, xlim = c(0, 1000), , type = "proportion",
#'              DiffCols = c(IA_1_P = "Target", IA_2_P = "Rhyme"),
#'            Condition1 = "talker", Condition2 = NA, Cond1Labels = c(CH1 = "Chinese 1", 
#'            CH10 = "Chinese 3", CH9 = "Chinese 2", EN3 = "English 1"),
#'            Cond2Labels = NA, ErrorBar = TRUE, 
#'            VWPreTheme = TRUE, ErrorBand = FALSE, 
#'              ErrorType = "CI", ConfLev = 95, CItype = "pointwise")
#'            
#' # For plotting conditional average differences (two conditions) with the 
#' # included theme and 95% simultaneous CI bands.
#' plot_avg_diff(data = dat, xlim = c(0, 1000), , type = "proportion",
#'              DiffCols = c(IA_1_P = "Target", IA_2_P = "Rhyme"),
#'            Condition1 = "talker", Condition2 = "Exp", Cond1Labels = c(CH1 = "Chinese 1", 
#'            CH10 = "Chinese 3", CH9 = "Chinese 2", EN3 = "English 1"),
#'            Cond2Labels = c(High = "H Exp", Low = "L Exp"), ErrorBar = FALSE, 
#'            VWPreTheme = TRUE, ErrorBand = TRUE, 
#'              ErrorType = "CI", ConfLev = 95, CItype = "simultaneous")
#' 
#' # For a more complete tutorial on VWPre plotting functions:
#' vignette("SR_Plotting", package="VWPre")
#' }
#' 
plot_avg_diff <- function(data, DiffCols = NULL, xlim = NA, type = NULL,
                          Averaging = "Event", 
                          Condition1 = NULL, Condition2 = NULL, Cond1Labels = NA, 
                          Cond2Labels = NA, ErrorBar = TRUE, VWPreTheme = TRUE,
                          ConfLev = 95, CItype = "simultaneous",
                          ErrorBand = FALSE, ErrorType = "SE") {
  
  .check_for_PupilPre(type = "UseOther", suggest = "ppl_plot_avg_diff")
  
  if(is.null(type)){
    stop("Please supply the plot type!")
  }
  
  if(is.null(DiffCols)){
    stop("Please supply the columns for the difference!")
  } else {
    DiffCol1 <- names(DiffCols)[1]
    if(!(DiffCol1 %in% colnames(data))){
      stop(paste(DiffCol1, "column not present in data!"))
    }
    DiffCol2 <- names(DiffCols)[2]
    if(!(DiffCol2 %in% colnames(data))){
      stop(paste(DiffCol2, "column not present in data!"))
    }
  }
  
  if(!(Averaging %in% colnames(data))){
    stop(paste(Averaging, " column not present in data!"))
  }
  
  Theme <- VWPreTheme
  
  if(ErrorBar==TRUE && ErrorBand==TRUE){
    stop(paste("Please select either error bars OR error bands. 
               For error bars set ErrorBar=TRUE and ErrorBand=FALSE.
               For error bands set ErrorBar=FALSE and ErrorBand=TRUE."))
  }
  
  # Check condition columns
  if(!(is.null(Condition1))) {
    if(is.na(Condition1)) {
      stop(paste("Please use NULL when not specifying a condition."))
    } else if (!(Condition1 %in% colnames(data))) {
      stop(paste(Condition1, " column not present in data!"))
    }
  }
  if(!(is.null(Condition2))) {
    if(is.na(Condition2)) {
      stop(paste("Please use NULL when not specifying a condition."))
    } else if (!(Condition2 %in% colnames(data))) {
      stop(paste(Condition2, " column not present in data!"))
    }
  }
  
  # Inform about Grand Average calculation
  message(paste0("Grand average of Difference calculated using ", Averaging, " means."))
  
  if (is.na(xlim[1])) {
    xlim <- c(range(data$Time)[1], range(data$Time)[2])
    xaxis <- unique(data$Time)
  } else {
    xlim <- xlim
    data<-data[data$Time>=xlim[1] & data$Time<=xlim[2],]
    xaxis <- unique(data$Time)
  }
  
  
  tmpdata <- data %>% 
    mutate(., DC1 = UQ(sym(DiffCol1)), DC2 = UQ(sym(DiffCol2)))
  
  # set column selection
  if(is.null(Condition1) && is.null(Condition2)){
    sel_names = quos(Time, DC1, DC2, UQ(sym(Averaging)))
  } else if(!is.null(Condition1) && is.null(Condition2)){
    sel_names = quos(Time, DC1, DC2, UQ(sym(Condition1)), UQ(sym(Averaging)))
  } else if(is.null(Condition1) && !is.null(Condition2)){
    sel_names = quos(Time, DC1, DC2, UQ(sym(Condition2)), UQ(sym(Averaging)))
  } else if(!is.null(Condition1) && !is.null(Condition2)){
    sel_names = quos(Time, DC1, DC2, UQ(sym(Condition1)), UQ(sym(Condition2)), UQ(sym(Averaging)))
  }
  
  # set grouping
  if(is.null(Condition1) && is.null(Condition2)){
    group_names1 = quos(UQ(sym(Averaging)), Time)
    group_names2 = quos(Time)
  } else if(!is.null(Condition1) && is.null(Condition2)){
    group_names1 = quos(UQ(sym(Averaging)), Time, UQ(sym(Condition1)))
    group_names2 = quos(Time, UQ(sym(Condition1)))
  } else if(is.null(Condition1) && !is.null(Condition2)){
    group_names1 = quos(UQ(sym(Averaging)), Time, UQ(sym(Condition2)))
    group_names2 = quos(Time, UQ(sym(Condition2)))
  } else if(!is.null(Condition1) && !is.null(Condition2)){
    group_names1 = quos(UQ(sym(Averaging)), Time, UQ(sym(Condition1)), UQ(sym(Condition2)))
    group_names2 = quos(Time, UQ(sym(Condition1)), UQ(sym(Condition2)))
  }
  
  tmpdata <- tmpdata %>% select(!!!sel_names) %>%
    group_by(!!!group_names1) %>% 
    summarise(DC1 = mean(DC1, na.rm = TRUE), DC2 = mean(DC2, na.rm = TRUE)) %>% 
    mutate(Diff = DC1 - DC2) %>% 
    group_by(!!!group_names2)
  
  if(type=="elogit") {
    tmpdata <- tmpdata %>%
      summarise(meanDiff = mean(Diff, na.rm = TRUE), DC1sd = stats::sd(DC1, na.rm = TRUE), DC2sd = stats::sd(DC2, na.rm = TRUE), n1 = n(), n2 = n()) %>%
      mutate(se = sqrt( ((DC1sd^2)/n1) + ((DC2sd^2)/n2) )) %>%
      mutate(degfree = (((((DC1sd^2)/n1) + ((DC2sd^2)/n2))^2) / (((((DC1sd^2)/n1)^2) / (n1-1)) + ((((DC2sd^2)/n2)^2) / (n2-1)))) )
    if(CItype=="pointwise") {
      tval <- 1-(((100-ConfLev)/2)/100)
    } else if (CItype=="simultaneous") {
      tval <- 1-(((100-ConfLev)/(2*length(unique(xaxis))))/100)
    }
    tmpdata <- tmpdata %>% mutate(ci = stats::qt(tval,df=degfree)*se)
  } else if (type=="proportion") {
    tmpdata <- tmpdata %>%
      summarise(meanDiff = mean(Diff, na.rm = TRUE), DC1m = mean(DC1, na.rm = TRUE), DC2m = mean(DC2, na.rm = TRUE), n1 = n(), n2 = n()) %>%
      mutate(se = sqrt(((DC1m*(1-DC1m))/n1)+((DC2m*(1-DC2m))/n2)))
    if(CItype=="pointwise") {
      zval <- 1-(((100-ConfLev)/2)/100)
    } else if (CItype=="simultaneous") {
      zval <- 1-(((100-ConfLev)/(2*length(unique(xaxis))))/100)
    }
    tmpdata <- tmpdata %>% mutate(ci = stats::qnorm(zval)*se)
  }
  # To be used for pooled variance of Elogit difference
  # if(type=="elogit") {
  #   tmpdata <- tmpdata %>%
  #     summarise(meanDiff = mean(Diff, na.rm = TRUE), DC1sd = stats::sd(DC1, na.rm = TRUE), DC2sd = stats::sd(DC2, na.rm = TRUE), n1 = n(), n2 = n()) %>%
  #     mutate(poolvar = (((n1-1)*(DC1sd^2)) + ((n2-1)*(DC2sd^2))) / (n1+n2-2)) %>%
  #     mutate(se = sqrt( ((poolvar^2)/n1) + ((poolvar^2)/n2) ))
  #   if(CItype=="pointwise") {
  #     tval <- 1-(((100-ConfLev)/2)/100)
  #   } else if (CItype=="simultaneous") {
  #     tval <- 1-(((100-ConfLev)/(2*length(unique(xaxis))))/100)
  #   }
  #   tmpdata <- tmpdata %>% mutate(ci = qt(tval,df=(n1+n2-2))*se)
  # } else if (type=="proportion") {
  #   tmpdata <- tmpdata %>%
  #     summarise(meanDiff = mean(Diff, na.rm = TRUE), DC1m = mean(DC1, na.rm = TRUE), DC2m = mean(DC2, na.rm = TRUE), n1 = n(), n2 = n()) %>%
  #     mutate(se = sqrt(((DC1m*(1-DC1m))/n1)+((DC2m*(1-DC2m))/n2)))
  #   if(CItype=="pointwise") {
  #     zval <- 1-(((100-ConfLev)/2)/100)
  #   } else if (CItype=="simultaneous") {
  #     zval <- 1-(((100-ConfLev)/(2*length(unique(xaxis))))/100)
  #   }
  #   tmpdata <- tmpdata %>% mutate(ci = qnorm(zval)*se)
  # }
  
  tmpdata <- ungroup(tmpdata)
  
  # Setting Error
  if(ErrorType=="SE"){
    tmpdata$error_lower <- tmpdata$meanDiff - tmpdata$se
    tmpdata$error_upper <- tmpdata$meanDiff + tmpdata$se
  } else if(ErrorType=="CI") {
    tmpdata$error_lower <- tmpdata$meanDiff - tmpdata$ci
    tmpdata$error_upper <- tmpdata$meanDiff + tmpdata$ci
  }
  
  # Setting ylim
  ylim <- c(0,0)
  if (type == "elogit") {
    ylim[1] = min(tmpdata$error_lower)
    ylim[2] = max(tmpdata$error_upper)
  } else if (type == "proportion") {
    ylim[1] = min(tmpdata$error_lower)
    ylim[2] = max(tmpdata$error_upper)
  }
  
  # Determine conditioning and legend labelling
  if (!is.null(Condition1) && !is.null(Condition2)) {
    tmpdata <- tmpdata %>% rename(Cond1=UQ(sym(Condition1)), Cond2=UQ(sym(Condition2)))
    if (any(!is.na(Cond1Labels))) {
      lev1 <- unique(levels(tmpdata$Cond1))
      for (x in 1:length(names(Cond1Labels))) {
        for(i in 1:length(lev1)) {
          if (lev1[i] == names(Cond1Labels)[x]) {
            lev1[i] <- Cond1Labels[[x]]
          }
        }
      }
      levels(tmpdata$Cond1) <- lev1
    }
    if (any(!is.na(Cond2Labels))) {
      lev2 <- unique(levels(tmpdata$Cond2))
      for (x in 1:length(names(Cond2Labels))) {
        for(i in 1:length(lev2)) {
          if (lev2[i] == names(Cond2Labels)[x]) {
            lev2[i] <- Cond2Labels[[x]]
          }
        }
      }
      levels(tmpdata$Cond2) <- lev2
    }
    tmpdata$Condition <- interaction(tmpdata$Cond1, tmpdata$Cond2, drop = TRUE, sep = "_")
    Cond <- TRUE
    CondName <- paste(Condition1, Condition2, sep = "_by_")
  }
  if (is.null(Condition1) && !is.null(Condition2)) {
    tmpdata <- tmpdata %>% rename(Condition=UQ(sym(Condition2)))
    if (any(!is.na(Cond2Labels))) {
      lev2 <- unique(levels(tmpdata$Condition))
      for (x in 1:length(names(Cond2Labels))) {
        for(i in 1:length(lev2)) {
          if (lev2[i] == names(Cond2Labels)[x]) {
            lev2[i] <- Cond2Labels[[x]]
          }
        }
      }
      levels(tmpdata$Condition) <- lev2
    }
    Cond <- TRUE
    CondName <- Condition2
  }
  if (!is.null(Condition1) && is.null(Condition2)) {
    tmpdata <- tmpdata %>% rename(Condition=UQ(sym(Condition1)))
    if (any(!is.na(Cond1Labels))) {
      lev1 <- unique(levels(tmpdata$Condition))
      for (x in 1:length(names(Cond1Labels))) {
        for(i in 1:length(lev1)) {
          if (lev1[i] == names(Cond1Labels)[x]) {
            lev1[i] <- Cond1Labels[[x]]
          }
        }
      }
      levels(tmpdata$Condition) <- lev1
    }
    Cond <- TRUE
    CondName <- Condition1
  }
  if (is.null(Condition1) && is.null(Condition2)) {
    Cond <- FALSE
  }
  
  # Print message regarding CIs
  if (ErrorType=="CI" && (ErrorBand==TRUE | ErrorBar==TRUE)) {
    if(CItype == "pointwise") {
      message(paste0("Plot created with ", CItype, " confindence intervals (set to ", ConfLev, "%)."))
    }
    if(CItype == "simultaneous") {
      if(type=="elogit") {
        message(paste0("Plot created with ", CItype, " confindence intervals (adjusted to ", round(tval*100, 2), "% using the Bonferroni method)."))
      } else if(type=="proportion") {
        message(paste0("Plot created with ", CItype, " confindence intervals (adjusted to ", round(zval*100, 2), "% using the Bonferroni method)."))
      }
    }
  }
  
  if (!is.null(Condition1) || !is.null(Condition2)) {
    plt <- ggplot(tmpdata, aes(x = Time, y = meanDiff, colour = Condition)) +
      ylab("Difference") +
      geom_hline(yintercept=0) +
      scale_x_continuous(limits = c(xlim[1], xlim[2])) +
      scale_y_continuous(limits = c(ylim[1], ylim[2]))
  } else {
    plt <- ggplot(tmpdata, aes(x = Time, y = meanDiff)) +
      ylab("Difference") +
      geom_hline(yintercept=0) +
      scale_x_continuous(limits = c(xlim[1], xlim[2])) +
      scale_y_continuous(limits = c(ylim[1], ylim[2]))
  }
  
  
  if (Theme == TRUE) {
    plt <- plt + .pac_theme()
  }
  
  # Set errorbands according to theme
  if(ErrorBand==TRUE & Theme==TRUE) {
    if(Cond==TRUE) {
      plt <- plt + aes(shape = Condition, colour = NULL) +
        geom_ribbon(aes(ymin=error_lower, ymax=error_upper, linetype=NA), fill="grey25", alpha=0.15, show.legend = FALSE)
      plt <- plt +
        scale_shape_discrete(name=CondName,
                             breaks=unique(tmpdata$Condition),
                             labels=unique(tmpdata$Condition))
    } else {
      plt <- plt + 
        geom_ribbon(aes(ymin=error_lower, ymax=error_upper, linetype=NA), fill="grey25", alpha=0.15, show.legend = FALSE)
      plt <- plt
    }
  }
  if(ErrorBand==TRUE & Theme==FALSE) {
    if(Cond==TRUE) {
      plt <- plt +
        geom_ribbon(aes(ymin=error_lower, ymax=error_upper, linetype=NA, fill=Condition), alpha=0.15, show.legend = FALSE)
      plt <- plt +
        scale_colour_hue(name=CondName,
                         breaks=unique(tmpdata$Condition),
                         labels=unique(tmpdata$Condition))
    } else {
      plt <- plt +
        geom_ribbon(aes(ymin=error_lower, ymax=error_upper, linetype=NA), alpha=0.15, show.legend = FALSE)
      plt <- plt
    }
  }
  
  # Set errorbars according to theme and plot grouping
  if(ErrorBar==TRUE & Theme==TRUE) {
    if(Cond==TRUE) {
      plt <- plt + 
        geom_errorbar(aes(ymin=error_lower, ymax=error_upper, colour = Condition), width = .3)
      plt <- plt +
        scale_color_grey(name=CondName,
                         breaks=unique(tmpdata$Condition),
                         labels=unique(tmpdata$Condition))
    } else {
      plt <- plt + 
        geom_errorbar(aes(ymin=error_lower, ymax=error_upper), width = .3)
      plt <- plt
    }
  }
  if(ErrorBar==TRUE & Theme==FALSE) {
    
    if(Cond==TRUE) {
      plt <- plt +
        geom_errorbar(aes(ymin=error_lower, ymax=error_upper, color=Condition), width = .3)
      plt <- plt +
        scale_colour_hue(name=CondName,
                         breaks=unique(tmpdata$Condition),
                         labels=unique(tmpdata$Condition))
    } else {
      plt <- plt +
        geom_errorbar(aes(ymin=error_lower, ymax=error_upper), width = .3)
      plt <- plt
    }
  }
  
  plt +
    geom_point() +
    geom_line() +
    ggtitle(paste("Difference between", DiffCols[[1]], "and", DiffCols[[2]], sep=" "))
  
  
  }


#' Plots average difference between two conditions.
#' 
#' \code{plot_avg_cdiff} calculates the average of differences between 
#' two specified conditions along with standard error and then plots the 
#' results.
#' 
#' @export
#' @import ggplot2
#' @import dplyr
#' @import rlang
#' 
#' @param data A data table object output by either \code{\link{bin_prop}}. 
#' \code{\link{transform_to_elogit}}, or \code{\link{create_binomial}}.
#' @param type A character string indicating "proportion" or "elogit",   
#' which influences how standard error and confidence intervals are calculated.
#' @param xlim A vector of two integers specifying the limits of the x-axis.
#' @param IAColumn A character vector specifying the desired column  
#' corresponding to the interest area. 
#' @param Averaging A character string indicating how the averaging should 
#' be done. "Subject" (default) will produce the grand mean in the data, while
#' "Item" (or, in principle, any other column name) will
#' calculate the grand mean by that factor.
#' @param Condition A list containing the column name corresponding to the 
#' condition and factor levels to be used for calculating the difference. 
#' @param CondLabels A named character vector specifying the desired labels 
#' of the levels of the condition. 
#' @param ErrorBar A logical indicating whether error bars should be
#' included in the plot.
#' @param ErrorBand A logical indicating whether error bands should be
#' included in the plot.
#' @param ErrorType A string indicating "SE" or "CI".
#' @param ConfLev A number indicating the confidence level of the CI.
#' @param CItype A string indicating "simultaneous" or "pointwise". Simultaneous
#' performs a Bonferroni correction for the interval.
#' @param VWPreTheme A logical indicating whether the theme included with the 
#' function should be applied, or ggplot2's base theme (to which any other 
#' custom theme could be added).
#' @examples
#' \dontrun{
#' library(VWPre)
#' # For plotting average difference between conditions...
#' plot_avg_cdiff(data = dat, xlim = c(0, 1000), type = "proportion",
#'              IAColumn = "IA_1_P", Condition = list(talker = c("EN3", "CH1")),
#'              CondLabels = NA, ErrorBar = TRUE, VWPreTheme = TRUE, 
#'              ErrorBand = FALSE, ErrorType = "SE")
#' 
#' # For a more complete tutorial on VWPre plotting functions:
#' vignette("SR_Plotting", package="VWPre")
#' }
#' 
plot_avg_cdiff <- function(data, IAColumn = NULL, xlim = NA, type = NULL,
                          Averaging = "Subject", 
                          Condition = NULL, CondLabels = NA, 
                          ErrorBar = TRUE, VWPreTheme = TRUE,
                          ConfLev = 95, CItype = "simultaneous",
                          ErrorBand = FALSE, ErrorType = "SE") {
  
  # Check if PupilPre is installed
  .check_for_PupilPre(type = "UseOther", suggest = "ppl_plot_avg_cdiff")
  
  if(is.null(type)){
    stop("Please supply the plot type!")
  }
  
  Column <- IAColumn
  Theme <- VWPreTheme
  
  {## SHARED CODE BEGINS HERE ## 
    
    if(is.null(Column)){
      stop("Please supply the column for which the difference will be calculated!")
    } else {
      if(!(Column %in% colnames(data))){
        stop(paste(Column, "column not present in data!"))
      }
    }
    
    # Check condition column
    if(is.null(Condition)){
      stop("Please supply a condition column along with the two levels desired for calculating the difference!")
    } 
    if(!(is.null(Condition))) {
      if (!(names(Condition)[1] %in% colnames(data))) {
        stop(paste(names(Condition)[1], " column not present in data!"))
      }
      if (length(Condition[[1]]) != 2) {
        stop(paste("Please ensure that two levels are specified for column", names(Condition)[1]))
      }
      if (Condition[[1]][1] %in% unique(data[[names(Condition)[1]]]) && 
          Condition[[1]][2] %in% unique(data[[names(Condition)[1]]])) {
      } else {
        stop(paste("One (or more) of the specified levels do not exist in column", names(Condition)[1]))
      }
    }
    
    if(is.null(Averaging)){
    } else if(!(Averaging %in% colnames(data))){
      stop(paste(Averaging, "column not present in data!"))
    }
    
    if(ErrorBar==TRUE && ErrorBand==TRUE){
      stop(paste("Please select either error bars OR error bands. 
                 For error bars set ErrorBar=TRUE and ErrorBand=FALSE.
                 For error bands set ErrorBar=FALSE and ErrorBand=TRUE."))
    }
    
    
    # Inform about Grand Average calculation
	if(is.null(Averaging)){
	message(paste0("Average of difference between ", Condition[[1]][1], " and ", Condition[[1]][2], " calculated using sample means."))
    } else {
      message(paste0("Grand average of difference between ", Condition[[1]][1], " and ", Condition[[1]][2], " calculated using ", Averaging, " means."))
    }
    
    # Prepare data
    
    if (is.na(xlim[1])) {
      xlim <- c(range(data$Time)[1], range(data$Time)[2])
      xaxis <- unique(data$Time)
    } else {
      xlim <- xlim
      data<-data[data$Time>=xlim[1] & data$Time<=xlim[2],]
      xaxis <- unique(data$Time)
    }  
    
    # Custom labels
    if(any(is.na(CondLabels))){
      main <- paste("Difference between", Condition[[1]][1], "and", Condition[[1]][2], "\n")
    } else {
      main <- paste("Difference between", 
                    CondLabels[which(names(CondLabels) == Condition[[1]][1])], "and", 
                    CondLabels[which(names(CondLabels) == Condition[[1]][2])], "\n")
      
    }
    
    tmpdata <- data %>%
      rename(Cond1 = UQ(sym(names(Condition)[1]))) %>% 
      filter(Cond1 %in% Condition[[1]]) %>% droplevels()
    tmpdata <- tmpdata %>% mutate(Cond = case_when(
      Cond1 == Condition[[1]][1] ~ "Lev1",
      Cond1 == Condition[[1]][2] ~ "Lev2"))
    tmpdata$Cond <- as.factor(tmpdata$Cond)
    
    
    # set column selection
	if(!is.null(Averaging)) {
    sel_names = quos(UQ(sym(Column)), Time, Cond, UQ(sym(Averaging)))
	} else {
	sel_names = quos(UQ(sym(Column)), Time, Cond)
	}
    
    # set grouping
	if(!is.null(Averaging)) {
    group_names1 = quos(Time, Cond, UQ(sym(Averaging)))
	}
    group_names2 = quos(Time, Cond)
    
    # Check balance of data
	if(!is.null(Averaging)) {
    chk <- tmpdata %>% group_by(UQ(sym(Averaging)), Cond) %>% summarise()
    chk <- unique(rowSums(table(chk))/2)
    if(length(chk) > 1){
      stop(paste0("Both levels of ", names(Condition)[1]," are not present for each level of", Averaging, ". Please select a different factor for \"Averaging\" and/or \"Condition\". Note that \"Averaging\" can be set to NULL so that sample means are used."))
    } else {
      if(chk == 1) {
      } else {
        stop(paste0("Both levels of ", names(Condition)[1]," are not present for each level of", Averaging, ". Please select a different factor for \"Averaging\" and/or \"Condition\". Note that \"Averaging\" can be set to NULL so that sample means are used."))
      }
    }
	}
    
    # Averaging
    tmpdata <- tmpdata %>% select(!!!sel_names) 
    
    if(!is.null(Averaging)) {
      tmpdata <- tmpdata %>% 
        group_by(!!!group_names1) %>% 
        summarise(Mean = mean(UQ(sym(Column)), na.rm = TRUE)) %>% 
        group_by(!!!group_names2) %>% 
        tidyr::spread(Cond, Mean) %>% arrange(UQ(sym(Averaging)), Time) %>% 
        mutate(Diff = Lev1 - Lev2)
    } else {
      tmpdata <- tmpdata %>% 
        group_by(!!!group_names2) %>% 
        summarise(Mean = mean(UQ(sym(Column)), na.rm = TRUE)) %>% 
        tidyr::spread(Cond, Mean) %>% 
        arrange(Time) %>% 
        mutate(Diff = Lev1 - Lev2)
    }
    
    if(type != "proportion") {
      tmpdata <- tmpdata %>%
        summarise(meanDiff = mean(Diff, na.rm = TRUE), 
                  DC1m = mean(Lev1, na.rm = TRUE), 
                  DC2m = mean(Lev2, na.rm = TRUE), 
                  DC1sd = stats::sd(Lev1, na.rm = TRUE), 
                  DC2sd = stats::sd(Lev2, na.rm = TRUE), 
                  n1 = n(), n2 = n()) %>%
        mutate(se = sqrt( ((DC1sd^2)/n1) + ((DC2sd^2)/n2) )) %>%
        mutate(degfree = (((((DC1sd^2)/n1) + ((DC2sd^2)/n2))^2) / (((((DC1sd^2)/n1)^2) / (n1-1)) + ((((DC2sd^2)/n2)^2) / (n2-1)))) )
      if(CItype=="pointwise") {
        tval <- 1-(((100-ConfLev)/2)/100)
      } else if (CItype=="simultaneous") {
        tval <- 1-(((100-ConfLev)/(2*length(unique(xaxis))))/100)
      }
      tmpdata <- tmpdata %>% mutate(ci = stats::qt(tval,df=degfree)*se)
    } else if (type=="proportion") {
      tmpdata <- tmpdata %>%
        summarise(meanDiff = mean(Diff, na.rm = TRUE), 
                  DC1m = mean(Lev1, na.rm = TRUE), 
                  DC2m = mean(Lev2, na.rm = TRUE), 
                  n1 = n(), n2 = n()) %>%
        mutate(se = sqrt(((DC1m*(1-DC1m))/n1)+((DC2m*(1-DC2m))/n2)))
      if(CItype=="pointwise") {
        zval <- 1-(((100-ConfLev)/2)/100)
      } else if (CItype=="simultaneous") {
        zval <- 1-(((100-ConfLev)/(2*length(unique(xaxis))))/100)
      }
      tmpdata <- tmpdata %>% mutate(ci = stats::qnorm(zval)*se)
    }
    
    tmpdata <- ungroup(tmpdata)
    
    # Setting Error
    if(ErrorType=="SE"){
      tmpdata$error_lower <- tmpdata$meanDiff - tmpdata$se
      tmpdata$error_upper <- tmpdata$meanDiff + tmpdata$se
    } else if(ErrorType=="CI") {
      tmpdata$error_lower <- tmpdata$meanDiff - tmpdata$ci
      tmpdata$error_upper <- tmpdata$meanDiff + tmpdata$ci
    }
    
    # Setting ylim
    ylim <- c(0,0)
    if (type != "proportion") {
      ylim[1] = min(tmpdata$error_lower)
      ylim[2] = max(tmpdata$error_upper)
    } else if (type == "proportion") {
      ylim[1] = min(tmpdata$error_lower)
      ylim[2] = max(tmpdata$error_upper)
    }
    
    # Print message regarding CIs
    if (ErrorType=="CI" && (ErrorBand==TRUE | ErrorBar==TRUE)) {
      if(CItype == "pointwise") {
        message(paste0("Plot created with ", CItype, " confindence intervals (set to ", ConfLev, "%)."))
      }
      if(CItype == "simultaneous") {
        if(type=="elogit") {
          message(paste0("Plot created with ", CItype, " confindence intervals (adjusted to ", round(tval*100, 2), "% using the Bonferroni method)."))
        } else if(type=="proportion") {
          message(paste0("Plot created with ", CItype, " confindence intervals (adjusted to ", round(zval*100, 2), "% using the Bonferroni method)."))
        }
      }
    }
    
    plt <- ggplot(tmpdata, aes(x = Time, y = meanDiff)) +
      ylab("Difference") +
      geom_hline(yintercept=0) +
      scale_x_continuous(limits = c(xlim[1], xlim[2])) +
      scale_y_continuous(limits = c(ylim[1], ylim[2]))
    
    if (Theme == TRUE) {
      plt <- plt + .pac_theme()
    }
    
    # Set errorbands according to theme
    if(ErrorBand==TRUE & Theme==TRUE) {
      plt <- plt + 
        geom_ribbon(aes(ymin=error_lower, ymax=error_upper, linetype=NA), fill="grey25", alpha=0.15, show.legend = FALSE)
      plt <- plt
    }
    
    if(ErrorBand==TRUE & Theme==FALSE) {
      plt <- plt +
        geom_ribbon(aes(ymin=error_lower, ymax=error_upper, linetype=NA), alpha=0.15, show.legend = FALSE)
      plt <- plt
    }
    
    # Set errorbars according to theme and plot grouping
    if(ErrorBar==TRUE & Theme==TRUE) {
      plt <- plt + 
        geom_errorbar(aes(ymin=error_lower, ymax=error_upper), width = .3)
      plt <- plt
    }
    
    if(ErrorBar==TRUE & Theme==FALSE) {
      plt <- plt +
        geom_errorbar(aes(ymin=error_lower, ymax=error_upper), width = .3)
      plt <- plt
    }
    
    plt +
      geom_point() +
      geom_line() +
      ggtitle(main)
    
  } ## SHARED CODE ENDS HERE ##
  
}


#' Create function for back-transforming empirical logits to proportions
#' 
#' \code{make_pelogit_fnc} creates a function that can transform empirical logit
#' values back to probability scale using the number of samples and constant
#' that were used in the original transformation. This function can then be use
#' to backtransform value predicted by a statistical model.
#' 
#' @export
#' 
#' @param ObsPerBin A positive integer indicating the number of observations 
#' used in the original transformation calculation. 
#' @param Constant A positive number used in the original transformation calculation.
#' @return A function.
#' @examples
#' \dontrun{
#' library(VWPre)
#' # Make backtransformation function
#' pelogit <- make_pelogit_fnc(20, 0.5)
#' } 
make_pelogit_fnc <- function(ObsPerBin = NULL, Constant = NULL) {
  
  .check_for_PupilPre(type = "NotAvailable")
  
  if(is.null(ObsPerBin)){
    stop("Please supply the observations per bin used in the original transformation!")
  }
  if(is.null(Constant)){
    stop("Please supply the constant used in the original transformation!")
  }
  fun <- function(predval){
    (((ObsPerBin*exp(predval) + Constant*exp(predval) - Constant) / (1 + exp(predval))) / ObsPerBin )
  }
  return(fun)
}

Try the VWPre package in your browser

Any scripts or data that you put into this service are public.

VWPre documentation built on Jan. 11, 2020, 9:28 a.m.