Nothing
R/plotting.R
In VWPre: Tools for Preprocessing Visual World Data
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 Nov. 30, 2020, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.