Plots

# Set general settings for markdown file
knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  comment = "",
  results = "hold")


# Clear environment
rm(list = ls())


# Enable/disable caching of time-consuming code chunks
knitr_cache_enabled = TRUE


# Load packages
library(dplyr)      # for data wrangling
library(ggplot2)    # for plotting
library(ggsignif)   # for adding asterisks
library(cowplot)    # for arranging plots
library(purrr)      # for calculating within-participant CIs
library(eegUtils)   # for plotting EEG topographies
library(tidyr)      # for gather function
library(sjPlot)     # for plot_model function
library(ggeffects)  # for extracting certain effects from model


# Load functions
source("./functions/summarySEwithinO.R")  # Function provided by R-cookbook: http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/
source("./functions/R_rainclouds.R")      # Function to create raincloud plots


# Define function to create common legend (from http://www.sthda.com/english/wiki/wiki.php?id_contents=7930#add-a-common-legend-for-multiple-ggplot2-graphs)
get_legend <- function(myggplot) {
  tmp      <- ggplot_gtable(ggplot_build(myggplot))
  leg      <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend   <- tmp$grobs[[leg]]
  return(legend)}


# Set figure theme
my_figure_theme <- theme_classic(base_size = 11) +
  theme(legend.position  = "bottom",
        axis.ticks.x     = element_blank(),
        plot.title       = element_text(hjust = 0.5),
        axis.title.x     = element_text(vjust = -1),
        strip.background = element_rect(color = "gray84", fill = "gray94", size = 0.3))


# Set figure colors
colors <- c("#2d5a7d", "tan1")
colors_erp <- c("royalblue1", "navy", "tan1", "sienna3")

Plot Behavioral Data (Observed)

---

This figure corresponds to Figure 2 in the manuscript.

# Load aggregated data per participant
load(file = "./saved_objects_for_plots/data_aggregated_rt.Rda")
load(file = "./saved_objects_for_plots/data_aggregated_acc.Rda")


# Define facet labels
group.labs <- c("Healthy Control Participants", "Patients With OCD")
names(group.labs) <- c("HC", "OCD")


#### Create RT plot ####

# Create raincloud plot
plot_rt <- ggplot() +
  geom_flat_violin(data = data_aggregated_rt, aes(x = gng_response_type, y = word_rt, fill = word_valence),
    position = position_nudge(x = .16, y = 0), adjust = 1.5, trim = FALSE, alpha = .7, color = NA) +
  geom_point(data = data_aggregated_rt, aes(x = as.numeric(gng_response_type) - .23, y = word_rt, color = word_valence),
    position = position_jitter(width = .05), size = 1, shape = 20, alpha = 0.8) +
  geom_boxplot(data = data_aggregated_rt, aes(x = gng_response_type, y = word_rt, fill = word_valence),
    outlier.shape = NA, alpha = 0.9, width = .3, color = "black") +
  labs(x = "Response in the Go/No-Go Task", y = "Word Categorization RT (ms)") +
  coord_cartesian(ylim = c(200, 1600)) +
  scale_y_continuous(breaks = seq(200, 1600, 200), expand = c(0, 0)) +
  scale_x_discrete(labels = c("SH" = "Slow Hit", "FH" = "Fast Hit", "FA" = "False Alarm", "IR" = "Inhibition")) +
  scale_color_manual(name = "Word Valence:", labels = c("Negative", "Positive"), values = colors,
    guide = guide_legend(override.aes = list(size = 6, shape = 15, alpha = 1), label.hjust = -2)) +
  scale_fill_manual(name  = "Word Valence:", labels = c("Negative", "Positive"), values = colors, guide = "none") +
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  my_figure_theme 


# Save and remove legend
legend <- get_legend(plot_rt)
plot_rt <- plot_rt + theme(legend.position = "none")


# Add position of asterisks
asteriks_rt <- data.frame(
  x1    = c(1.85, 2.85, 1.85, 2.85),
  x2    = c(2.15, 3.15, 2.15, 3.15),
  y1    = c(1320, 1320, 1320, 1320),
  xstar = c(2,    3,    2,    3),
  ystar = c(1330, 1330, 1330, 1330),
  lab   = c("*", "***", "**", "*"),
  group = c("HC", "HC",  "OCD", "OCD"))


# Add position of asterisks across groups
asteriks_rt2 <- data.frame(
  x1    = .370,
  x2    = .820,
  y1    = .820,
  xstar = .595,
  ystar = .830,
  lab   = "**")


# Add a new layer on-top to add asterisks
plot_rt <- plot_rt +
  geom_text(data = asteriks_rt,    aes(x = xstar, y = ystar, label = lab, size = 2.5),    color = "gray30") +
  geom_segment(data = asteriks_rt, aes(x = x1, xend = x2, y = y1, yend = y1), size = 0.3, color = "gray30") +
  facet_grid(. ~ group, labeller = labeller(group = group.labs))


# Add a new layer on-top to add asterisks across groups
plot_rt <- ggdraw(plot_rt) +
  geom_text(data = asteriks_rt2,    aes(x = xstar, y = ystar, label = lab, size = 2.5),    color = "gray30") +
  geom_segment(data = asteriks_rt2, aes(x = x1, xend = x2, y = y1, yend = y1), size = 0.3, color = "gray30")


#### Create accuracy plot ####

# Create raincloud plot
plot_acc <- ggplot() +
  geom_flat_violin(data = data_aggregated_acc, aes(x = gng_response_type, y = accuracy, fill = word_valence),
    position = position_nudge(x = .16, y = 0), adjust = 1.5, trim = FALSE, alpha = .7, color = NA) +
  geom_point(data = data_aggregated_acc, aes(x = as.numeric(gng_response_type) - .23, y = accuracy, color = word_valence),
    position = position_jitter(width = .05), size = 1, shape = 20, alpha = 0.8) +
  geom_boxplot(data = data_aggregated_acc, aes(x = gng_response_type, y = accuracy, fill = word_valence),
    outlier.shape = NA, alpha = 0.9, width = .3, color = "black") +
  labs(x = "Response in the Go/No-Go Task", y = "Word Categorization Accuracy (%)") +
  coord_cartesian(ylim = c(0, 105)) +
  scale_y_continuous(breaks = seq(0, 100, 10), expand = c(0, 0)) +
  scale_x_discrete(labels = c("SH" = "Slow Hit", "FH" = "Fast Hit", "FA" = "False Alarm", "IR" = "Inhibition")) +
  scale_color_manual(name = "Word Valence:", labels = c("Negative", "Positive"), values = colors) +
  scale_fill_manual(name  = "Word Valence:", labels = c("Negative", "Positive"), values = colors) +
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  annotate("rect", xmin = 0.5, xmax = 4.5, ymin = 101, ymax = 105, fill = "white") + # Add whitespace below facet
  my_figure_theme + theme(legend.position = "none")  



# Add position of asterisks
asteriks_acc <- data.frame(
  x1    = c(2.85, 2.85),
  x2    = c(3.15, 3.15),
  y1    = c(19,   19),
  xstar = c(3,    3),
  ystar = c(15,   15),
  lab   = c("***", "***"),
  group = c("HC", "OCD"))


# Add a new layer on-top to add asterisks
plot_acc <- plot_acc +
  geom_text(data = asteriks_acc,    aes(x = xstar, y = ystar, label = lab, size = 2.5),    color = "gray30") +
  geom_segment(data = asteriks_acc, aes(x = x1, xend = x2, y = y1, yend = y1), size = 0.3, color = "gray30") +
  facet_grid(. ~ group, labeller = labeller(group = group.labs))


# Arrange plots
figure_obs_val <- ggdraw() +
  draw_plot(plot_rt,  x = .00, y = .490, width = 1,  height = .51) +
  draw_plot(plot_acc, x = .00, y = .050, width = 1,  height = .42) +
  draw_plot(legend,   x = .28, y = .015, width = .5, height = .00) +
  draw_plot_label(c("A", "B"), c(0, 0), c(1, 0.47), size = 15)
figure_obs_val

Note. (A) Observed word categorization RT per preceding response type in the go/no-go task, word valence, and group. (B) Observed word categorization accuracy per preceding response type in the go/no-go task, word valence, and group. (A–B) The plots show individual data points, boxplots, and probability density plots based on raw data that were aggregated by participant. Asterisks indicate significant differences yielded by the (generalized) linear mixed model analysis. * p < .05, ** p < .01, *** p < .001.

# Save plot
ggsave(plot = figure_obs_val, "./figures/figure_2.tiff", width = 16, height = 27,
  units = "cm", dpi = 600, compression = "lzw")

Plot Behavioral Data (Predicted)

---

This figure corresponds to Figure 3 in the manuscript.

# Load (G)LMM output
LMM_rt_final   <- readRDS(file = "./saved_objects_for_plots/LMM_rt_final.rds")
GLMM_acc_final <- readRDS(file = "./saved_objects_for_plots/GLMM_acc_final.rds")


# Define facet labels
group.labs <- c("Healthy Control Participants", "Patients With OCD")
names(group.labs) <- c("HC", "OCD")


#### Create RT plot ####

# Extract only FA and FH from LMM
LMM_rt        <- ggpredict(LMM_rt_final, terms = c("gng_response_type", "word_valence", "group"))
LMM_rt_plot   <- LMM_rt[LMM_rt$x %in% c("FH", "FA"), ]
LMM_rt_plot$x <- droplevels(LMM_rt_plot$x)


# Plot predicted means
plot_LMM_rt <- ggplot() +
  geom_point(data = LMM_rt_plot, aes(x = x, y = predicted, color = group),
    position = position_dodge(width = .25), size = 2.5) +
  geom_errorbar(data = LMM_rt_plot, aes(x = x, ymax = conf.low, ymin = conf.high,  color = group),
    position = position_dodge(width = .25), size = .7, width = .2) +
  labs(x = "Response in the Go/No-Go Task", y = "Predicted Word Categorization RT (ms)") +
  coord_cartesian(ylim = c(500, 950)) +
  scale_y_continuous(breaks = seq(500, 950, 50), expand = c(0, 0)) +
  scale_x_discrete(labels = c("FH" = "Fast Hit", "FA" = "False Alarm")) +
  scale_color_manual(name = "Word Valence:", labels = c("Negative", "Positive"), values = colors) +
  facet_wrap(~facet, labeller = labeller(facet = group.labs)) +
  my_figure_theme


# Save and remove legend
legend <- get_legend(plot_LMM_rt)
plot_LMM_rt  <- plot_LMM_rt + theme(legend.position = "none")


# Add position of asterisks
asteriks_rt <- data.frame(
  x1    = c(.225, .415, .655, .845, .430),
  x2    = c(.275, .465, .705, .895, .880),
  y1    = c(.800, .800, .800, .800, .860),
  xstar = c(.250, .440, .680, .870, .655),
  ystar = c(.810, .810, .810, .810, .870),
  lab   = c("*",  "***", "**", "*",  "**"))


# Add a new layer on-top to add asterisks
plot_LMM_rt <- ggdraw(plot_LMM_rt) +
  geom_text(data = asteriks_rt,    aes(x = xstar, y = ystar, label = lab, size = 2.5),    color = "gray30") +
  geom_segment(data = asteriks_rt, aes(x = x1, xend = x2, y = y1, yend = y1), size = 0.3, color = "gray30")


#### Create accuracy plot ####

# Extract only FA and FH from GLMM
GLMM_acc        <- ggpredict(GLMM_acc_final, terms = c("gng_response_type", "word_valence", "group"))
GLMM_acc_plot   <- GLMM_acc[GLMM_acc$x %in% c("FH", "FA"), ]
GLMM_acc_plot$x <- droplevels(GLMM_acc_plot$x)


# Plot predicted means
plot_GLMM_acc <- ggplot() +
  geom_point(data = GLMM_acc_plot, aes(x = x, y = predicted * 100, color = group),
    position = position_dodge(width = .25), size = 2.5) +
  geom_errorbar(data = GLMM_acc_plot, aes(x = x, ymax = conf.low * 100, ymin = conf.high * 100, color = group),
    position = position_dodge(width = .25), size = .7, width = .2) +
  labs(x = "Response in the Go/No-Go Task", y = "Predicted Word Categorization Accuracy (%)") +
  coord_cartesian(ylim = c(50, 110)) +
  scale_y_continuous(breaks = seq(50, 100, 5), expand = c(0, 0)) +
  scale_x_discrete(labels = c("FH" = "Fast Hit", "FA" = "False Alarm")) +
  scale_color_manual(name = "Word Valence:", labels = c("Negative", "Positive"), values = colors) +
  facet_wrap(~facet, labeller = labeller(facet = group.labs)) +
  my_figure_theme +
  theme(legend.position = "none")


# Add position of asterisks
asteriks_acc <- data.frame(
  x1    = c(.415, .845),
  x2    = c(.465, .895),
  y1    = c(.800, .800),
  xstar = c(.440, .870),
  ystar = c(.810, .810),
  lab   = c("***", "***"))


# Add a new layer on-top with ggdraw in order to add asterisks
plot_GLMM_acc <- ggdraw(plot_GLMM_acc) +
  geom_text(data = asteriks_acc,    aes(x = xstar,  y = ystar, label = lab, size = 2.5),   color = "gray30") +
  geom_segment(data = asteriks_acc, aes(x = x1, xend = x2, y = y1, yend = y1), size = 0.3, color = "gray30")


#### Arrange plots
figure_pred_val <- ggdraw() +
  draw_plot(plot_LMM_rt,   x = .00, y = .070, width = .49, height = .93) +
  draw_plot(plot_GLMM_acc, x = .52, y = .070, width = .49, height = .93) +
  draw_plot(legend,        x = .28, y = .015, width = .50, height = .00)   +
  draw_plot_label(c("A", "B"), c(0, 0.52), c(1, 1), size = 15)
figure_pred_val

Note. (A) Model-predicted mean word categorization RT per preceding response type (fast hit, false alarm) in the go/no-go task, word valence, and group, computed as partial effects from the linear mixed model (LMM). RT data were log-transformed for analysis but back-transformed to ms for visualization. (B) Model-predicted mean word categorization accuracy per preceding response type (fast hit, false alarm) in the go/no-go task, word valence, and group, computed as partial effects from the generalized linear mixed model (GLMM). (A–B) Error bars indicate 95% confidence intervals. Asterisks indicate significant differences yielded by the (G)LMM analysis. * p < .05, ** p < .01, *** p < .001.

# Save plot
ggsave(plot = figure_pred_val, "./figures/figure_3.tiff", width = 20, height = 12,
  units = "cm", dpi = 600, compression = "lzw")

Plot ERP Data

---

This figure corresponds to Figure 4 in the manuscript.

# Load data
resp_locked_data <- read.csv(file = "./data/response_locked_data_for_plots.csv", header = TRUE)


# Create relevant variables and filter data
resp_locked_data <- resp_locked_data %>%
  dplyr::mutate(
    stimulation = factor(ifelse(
      participant_id == "C_01_T1" | participant_id == "C_02_T1" | participant_id == "C_03_T2" |
      participant_id == "C_04_T2" | participant_id == "C_05_T1" | participant_id == "C_06_T2" |
      participant_id == "C_07_T1" | participant_id == "C_08_T1" | participant_id == "C_09_T2" |
      participant_id == "C_10_T2" | participant_id == "C_11_T2" | participant_id == "C_12_T2" |
      participant_id == "C_13_T1" | participant_id == "C_15_T2" | participant_id == "C_16_T1" |
      participant_id == "C_17_T2" | participant_id == "C_18_T2" | participant_id == "C_19_T1" |
      participant_id == "C_20_T1" | participant_id == "C_21_T2" | participant_id == "C_22_T2" |
      participant_id == "C_23_T1" | participant_id == "C_24_T1" | participant_id == "C_25_T1" |
      participant_id == "C_26_T1" | participant_id == "C_27_T2" | participant_id == "C_28_T2" |
      participant_id == "C_29_T1" | participant_id == "C_30_T2" | participant_id == "P_01_T1" |
      participant_id == "P_02_T2" | participant_id == "P_03_T2" | participant_id == "P_04_T1" |
      participant_id == "P_05_T1" | participant_id == "P_06_T1" | participant_id == "P_07_T2" |
      participant_id == "P_08_T2" | participant_id == "P_09_T2" | participant_id == "P_10_T1" |
      participant_id == "P_11_T2" | participant_id == "P_12_T2" | participant_id == "P_13_T1" |
      participant_id == "P_15_T1" | participant_id == "P_16_T1" | participant_id == "P_17_T2" |
      participant_id == "P_18_T1" | participant_id == "P_19_T1" | participant_id == "P_20_T2" |
      participant_id == "P_21_T2" | participant_id == "P_22_T1" | participant_id == "P_23_T1" |
      participant_id == "P_24_T2" | participant_id == "P_25_T1" | participant_id == "P_26_T2" |
      participant_id == "P_27_T1" | participant_id == "P_28_T2" | participant_id == "P_29_T2" |
      participant_id == "P_30_T1", "sham", "verum")),
    condition      = factor(ifelse(condition == 1, "correct", "incorrect"), levels = c("incorrect", "correct")),
    group          = factor(ifelse(substr(participant_id, 1, 1) == "C", "HC", "OCD")),
    session        = factor(ifelse(substr(participant_id, 6, 7) == "T1", "T1", "T2")),
    participant_id = factor(substr(participant_id, 1, 4))) %>% # to get correct no. of factor levels
  # Only keep sham session and exclude P_02 and C_02
  dplyr::filter(stimulation == "sham" & participant_id != "P_02" & participant_id != "C_02")


# Update factor levels
resp_locked_data$participant_id <- droplevels(resp_locked_data$participant_id)


# Calculate running within-participant CIs
running_ci_mfn <- resp_locked_data %>%
  split(.$time) %>%
  map(~ summarySEwithinO(
    data = .,
    measurevar  = "FCz",
    withinvars  = "condition",
    betweenvars = "group",
    idvar       = "participant_id"))


ci_mfn <- purrr::map_df(running_ci_mfn, magrittr::extract) %>%
  dplyr::mutate(time = rep(unique(resp_locked_data$time), each = 4))
# Note. 4 refers to  no. of conditions (group x response_type = 2 x 2 = 4)


# Convert to long format for electrodes
data_topo_mfn <- resp_locked_data %>%
  # Remove channels of no interest
  dplyr::select(-IO1, -M1, -F9, -F10) %>%
  # Change from wide to long format for electrodes
  tidyr::gather(., electrode, amplitude, Fp1:O2, factor_key = TRUE)  %>%
  # Select time windows
  dplyr::filter(time >= 0 & time <= 100) %>%
  # Add electrode information
  electrode_locations(., electrode = "electrode", drop = FALSE, montage = NULL)

#### Create waveform plot ####

# Define electrode labels
label_fcz <- data.frame(time = -320, FCz = 17.5, lab = "Text") 


# Create waveform plot
plot_waveform_mfn <- ggplot(resp_locked_data, aes(time, FCz)) +
  stat_summary(fun = mean, geom = "line", size = 0.5, linetype = "solid", aes(color = condition:group)) +
  geom_ribbon(data = ci_mfn, aes(ymin = FCz - ci, ymax = FCz + ci, fill = condition:group), alpha = 0.2) +
  guides(fill = "none") +
  geom_segment(aes(x = 0,    xend = 0,   y = -8, yend = 18), linetype = "dashed", color = "gray60", size = 0.2) +
  geom_segment(aes(x = -400, xend = 800, y =  0, yend = 0),  linetype = "dashed", color = "gray60", size = 0.2) +
  geom_text(data = label_fcz, label = "FCz") +
  annotate("rect", xmin = 0, xmax = 100, ymin = -8, ymax = 18, alpha = .2, fill = "gray50") + # Shaded time area
  labs(x = "Time (ms)", y = expression(paste("Amplitude (", mu, "V)"))) +
  coord_cartesian(ylim = c(-8, 18), xlim = c(-200, 400)) +
  scale_y_continuous(breaks = seq(-8, 18, 2), expand = c(0, 0)) +
  scale_x_continuous(breaks = seq(-200, 400, 100), expand = c(0, 0)) +
  scale_color_manual(values = colors_erp, labels = c("HC: False Alarm", "OCD: False Alarm", "HC: Hit", "OCD: Hit"),
    guide = guide_legend(override.aes = list(size = 1.3))) + # Thickness legend lines
  scale_fill_manual(values = colors_erp) +
  my_figure_theme + 
  theme(axis.ticks.x = NULL, legend.title = element_blank(), legend.text = element_text(size = 9))


#### Create topography plots ####

# Define facet labels
labs_incorr <- c("HC: False Alarm", "OCD: False Alarm")
names(labs_incorr) <- c("HC", "OCD")

labs_corr <- c("HC: Hit", "OCD: Hit")
names(labs_corr) <- c("HC", "OCD")


# ERN
plot_topo_ern <- ggplot(data_topo_mfn[data_topo_mfn$condition == "incorrect", ],
  aes(x = x, y = y, fill = amplitude, z = amplitude, label = electrode)) +
  geom_topo(grid_res = 300, interp_limit = "head", chan_markers = "point", chan_size = 0.1, 
            head_size = 0.2, colour = "gray40", size = 0.2) +
  scale_fill_distiller(palette = "RdBu", limits = c(-4, 4), breaks = c(-4, 0, 4)) +
  theme_void() +
  coord_equal() +
  labs(fill = expression(paste(mu, "V"))) +
  theme(legend.position = "bottom", legend.text = element_text(size = 8, color = "gray30"),
    legend.title = element_text(size = 8, color = "gray30", vjust = 0.9), legend.key.height = unit(0.25, 'cm'),
    legend.key.width = unit(0.3, 'cm'), strip.text.x = element_text(size = 9)) +
  facet_wrap(. ~group, ncol = 2, labeller = labeller(group = labs_incorr))


# Save and remove legend
legend_ern <- get_legend(plot_topo_ern)
plot_topo_ern <- plot_topo_ern + theme(legend.position = "none")


# CRN
plot_topo_crn <- ggplot(data_topo_mfn[data_topo_mfn$condition == "correct", ],
  aes(x = x, y = y, fill = amplitude, z = amplitude, label = electrode)) +
  geom_topo(grid_res = 300, interp_limit = "head", chan_markers = "point", chan_size = 0.1, 
            head_size = 0.2, colour = "gray40", size = 0.2) +
  scale_fill_distiller(palette = "RdBu", limits = c(-4, 8.4), breaks = c(-4, 0, 4, 8)) +
  theme_void() +
  coord_equal() +
  labs(fill = expression(paste(mu, "V"))) +
  theme(legend.position = "bottom", legend.text = element_text(size = 8, color = "gray30"),
    legend.title = element_text(size = 8, color = "gray30", vjust = 0.9), legend.key.height = unit(0.25, 'cm'),
    legend.key.width = unit(0.3, 'cm'), strip.text.x = element_text(size = 9)) +
  facet_wrap(. ~group, ncol = 2, labeller = labeller(group = labs_corr))


# Save and remove legend
legend_crn <- get_legend(plot_topo_crn)
plot_topo_crn <- plot_topo_crn + theme(legend.position = "none")


# Overlay waveform and topopgraphy plots
aligned_plots1_mfn <- align_plots(plot_waveform_mfn, plot_topo_ern, axis = "tr")
aligned_plots2_mfn <- align_plots(plot_waveform_mfn, plot_topo_crn, axis = "tr")


# Arrange plots
figure_ern_crn <- ggdraw() +
  draw_plot(aligned_plots1_mfn[[1]], x = 0,   y = .32, width = .99, height = .64) +
  draw_plot(aligned_plots1_mfn[[2]], x = .09, y = .05, width = .42, height = .21) +
  draw_plot(aligned_plots2_mfn[[2]], x = .55, y = .05, width = .42, height = .21) +
  draw_plot(legend_ern,              x = .07, y = 0,   width = .42, height = .05) +
  draw_plot(legend_crn,              x = .53, y = 0,   width = .42, height = .05) +
  draw_plot_label(c("A", "B"), c(0, 0), c(1, 0.285), size = 15)
figure_ern_crn

Note. (A) Grand average response-locked waveforms at electrode site FCz with 95% confidence intervals in patients with OCD and HC participants for correct (hit) and incorrect (false alarm) responses in the go/no-go task. The gray shaded time interval was used to quantify ERN and CRN amplitude. (B) Scalp topographies of ERN and CRN in both groups in the time window 0–100 ms after false alarms or hits, respectively.

# Save plot
ggsave(plot = figure_ern_crn, "./figures/figure_4.tiff", width = 13.5, height = 16,
  units = "cm", dpi = 600, compression = "lzw")

Plot Effects of OCD Characteristics

---

OCD Symptom Severity (OCI-R)

This figure corresponds to Figure 5 in the manuscript.

# Load LMM output and correlation data
LMM_rt_oci_ocd  <- readRDS(file = "./saved_objects_for_plots/LMM_rt_oci_ocd.rds")
load(file = "./saved_objects_for_plots/df_corr.Rda")


# Create plot of estimated interaction effects
plot_oci_1 <- plot_model(LMM_rt_oci_ocd,
                         type   = "pred",
                         terms  = c("oci_centered", "pos_neg", "FA_FH"),
                         ci.lvl = .95,
                         transform = "exp",
                         line.size = 1) +
  labs(title = NULL, y = "Word Categorization RT (ms)", x = "OCI-R Score (Group-Mean-Centered)") +
  coord_cartesian(ylim = c(500, 900), xlim = c(-13, 20)) +
  scale_y_continuous(breaks = seq(500, 900, 50), expand = c(0, 0)) +
  scale_x_continuous(breaks = seq(-10, 20, 5), expand = c(0, 0)) +
  scale_color_manual(name = "Word Valence:", values = colors, labels = c("Negative", "Positive")) +
  scale_fill_manual(name  = "Word Valence:", values = colors, labels = c("Negative", "Positive")) +
  my_figure_theme +
  theme(strip.text.x = element_text(size = 11), axis.ticks.x = NULL) +
  guides(fill = guide_legend(override.aes = list(color = colors, size = 0.65)))


# Change facet labels
plot_oci_1$data$facet <- ifelse(plot_oci_1$data$facet == "FA_FH = -0.5", "Fast Hits", "False Alarms")


# Create scatterplot
plot_oci_2 <- ggplot(data = df_corr[df_corr$group == "OCD", ], aes(x = oci, y = priming, fill = group)) +
  geom_point(color = "#003192", size = 1) +
  geom_smooth(method = 'lm', color = "#003192") +
  labs(title = NULL, x = "OCI-R Score", y = "Overall Priming Effect (ms)") +
  coord_cartesian(ylim = c(-150, 400), xlim = c(12.9, 45.5)) +
  scale_y_continuous(breaks = seq(-150, 400, 50), expand = c(0, 0)) +
  scale_x_continuous(breaks = seq(15, 45, 5), expand = c(0, 0)) +
  scale_fill_manual(values = "#003192") +
  my_figure_theme +
  theme(legend.position = "none", axis.ticks.x = NULL) +
  annotate("text", x = 45, y = -120, label = paste("list(italic(r)(26) ==", -.41, ", italic(p) == .029)"),
    parse = TRUE, hjust = 1, size = 3, color = "gray20")


# Arrange plots
figure_oci <- ggdraw() +
  draw_plot(plot_oci_1, x = .00, y = .000, width = .60, height = 1) +
  draw_plot(plot_oci_2, x = .61, y = .115, width = .39, height = .83) +
  draw_plot_label(c("A", "B"), c(0, .61), c(1, 1), size = 15)
figure_oci

Note. (A) Model-predicted interaction effect between Obsessive-Compulsive Inventory-Revised (OCI-R) score, preceding response type (fast hit, false alarm) in the go/no-go task, and word valence on word categorization response time (RT) in patients with OCD, computed as partial effects from the linear mixed model. RT data were log-transformed for analysis but were back-transformed to ms for data visualization. (B) Negative correlation between the OCI-R score and the overall priming effect across participants. The overall priming effect refers to the priming effect after false alarms and fast hits and was calculated by subtracting the RT in congruent conditions (positive words after fast hits and negative words after false alarms) from the RT in incongruent conditions (positive words after false alarms and negative words after fast hits) for each participant (Aarts et al., 2012). (A–B) Shaded bands represent 95% confidence intervals.

# Save plot
ggsave(plot = figure_oci, "./figures/figure_5.tiff", width = 20, height = 12,
  units = "cm", dpi = 600, compression = "lzw")

Trait Anxiety (STAI)

This figure corresponds to Figure S1 in the supplemental material.

# Load LMM output
LMM_rt_stai_ocd <- readRDS(file = "./saved_objects_for_plots/LMM_rt_stai_ocd.rds")


# Create plot of estimated interaction effects
plot_stai_1 <- plot_model(LMM_rt_stai_ocd,
                          type   = "pred",
                          terms  = c("stai_centered", "pos_neg", "FA_FH"),
                          ci.lvl = .95,
                          transform = "exp",
                          line.size = 1) +
  labs(title = NULL, y = "Word Categorization RT (ms)", x = "STAI Trait Score (Group-Mean-Centered)") +
  coord_cartesian(ylim = c(450, 900), xlim = c(-18, 25)) +
  scale_y_continuous(breaks = seq(450, 900, 50), expand = c(0, 0)) +
  scale_x_continuous(breaks = seq(-15, 25, 5), expand = c(0, 0)) +
  scale_color_manual(name = "Word Valence:", values = colors, labels = c("Negative", "Positive")) +
  scale_fill_manual(name  = "Word Valence:", values = colors, labels = c("Negative", "Positive")) +
  my_figure_theme +
  theme(strip.text.x = element_text(size = 11), axis.ticks.x = NULL) +
  guides(fill = guide_legend(override.aes = list(color = colors, size = 0.65)))


# Change facet labels
plot_stai_1$data$facet <- ifelse(plot_stai_1$data$facet == "FA_FH = -0.5", "Fast Hits", "False Alarms")


# Create scatterplot
plot_stai_2 <- ggplot(data = df_corr[df_corr$group == "OCD", ], aes(x = stai, y = priming, fill = group)) +
  geom_point(color = "#003192", size = 1) +
  geom_smooth(method = 'lm', color = "#003192") +
  labs(title = NULL, x = "STAI Trait Score", y = "Overall Priming Effect (ms)") +
  coord_cartesian(ylim = c(-150, 400), xlim = c(35.8, 80)) +
  scale_y_continuous(breaks = seq(-150, 400, 50), expand = c(0, 0)) +
  scale_x_continuous(breaks = seq(40, 80, 10), expand = c(0, 0)) +
  scale_fill_manual(values = "#003192") +
  my_figure_theme +
  theme(legend.position = "none", axis.ticks.x = NULL) +
  annotate("text", x = 80, y = -115, label = paste("list(italic(r)(26) ==", -.34, ", italic(p) == .078)"),
    parse = TRUE, hjust = 1, size = 3, color = "gray20")


# Arrange plots
figure_stai <- ggdraw() +
  draw_plot(plot_stai_1, x = .00, y = .000, width = .60, height = 1) +
  draw_plot(plot_stai_2, x = .61, y = .115, width = .39, height = 0.83) +
  draw_plot_label(c("A", "B"), c(0, .61), c(1, 1), size = 15)
figure_stai

Note. (A) Model-predicted interaction effect between State-Trait Anxiety Inventory (STAI) trait score, preceding response type (fast hit, false alarm) in the go/no-go task, and word valence on word categorization response time (RT) in patients with OCD, computed as partial effects from the linear mixed model. RT data were log-transformed for analysis but were back-transformed to ms for data visualization. (B) Negative correlation between the STAI trait score and the overall priming effect across participants. The overall priming effect refers to the priming effect after false alarms and fast hits and was calculated by subtracting the RT in congruent conditions (positive words after fast hits and negative words after false alarms) from the RT in incongruent conditions (positive words after false alarms and negative words after fast hits) for each participant (Aarts et al., 2012). (A–B) Shaded bands represent 95% confidence intervals.

# Save plot
ggsave(plot = figure_stai, "./figures/figure_S1.tiff", width = 20, height = 12,
  units = "cm", dpi = 600, compression = "lzw")

References

---

Aarts, K., De Houwer, J., & Pourtois, G. (2012). Evidence for the automatic evaluation of self-generated actions. Cognition, 124(2), 117–127. https://doi.org/10.1016/j.cognition.2012.05.009

Session Info

---

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] lmerTest_3.1-3  lme4_1.1-25     Matrix_1.5-1    ggeffects_1.1.3
 [5] sjPlot_2.8.11   tidyr_1.2.1     eegUtils_0.6.1  purrr_0.3.5    
 [9] cowplot_1.1.1   ggsignif_0.6.4  ggplot2_3.3.6   dplyr_1.0.10   

loaded via a namespace (and not attached):
  [1] nlme_3.1-140        matrixStats_0.58.0  RColorBrewer_1.1-3 
  [4] insight_0.18.4      httr_1.4.4          numDeriv_2016.8-1.1
  [7] tools_3.6.1         backports_1.4.1     bslib_0.4.0        
 [10] utf8_1.2.2          R6_2.5.1            sjlabelled_1.2.0   
 [13] DBI_1.1.3           lazyeval_0.2.2      mgcv_1.8-28        
 [16] colorspace_2.0-3    withr_2.5.0         tidyselect_1.2.0   
 [19] emmeans_1.6.0       compiler_3.6.1      textshaping_0.3.6  
 [22] performance_0.10.0  cli_3.4.1           plotly_4.10.0      
 [25] labeling_0.4.2      bayestestR_0.13.0   sass_0.4.0         
 [28] scales_1.2.1        mvtnorm_1.1-3       systemfonts_1.0.4  
 [31] stringr_1.4.1       digest_0.6.30       minqa_1.2.4        
 [34] rmarkdown_2.16      pkgconfig_2.0.3     htmltools_0.5.3    
 [37] parallelly_1.32.1   highr_0.9           fastmap_1.1.0      
 [40] htmlwidgets_1.5.4   rlang_1.0.6         rstudioapi_0.14    
 [43] shiny_1.7.2         farver_2.1.1        jquerylib_0.1.4    
 [46] generics_0.1.3      jsonlite_1.8.3      magrittr_2.0.3     
 [49] parameters_0.19.0   Rcpp_1.0.9          munsell_0.5.0      
 [52] fansi_1.0.3         abind_1.4-5         lifecycle_1.0.3    
 [55] stringi_1.6.1       yaml_2.2.1          MASS_7.3-51.4      
 [58] plyr_1.8.7          grid_3.6.1          parallel_3.6.1     
 [61] listenv_0.8.0       promises_1.2.0.1    sjmisc_2.8.9       
 [64] miniUI_0.1.1.1      lattice_0.20-45     splines_3.6.1      
 [67] sjstats_0.18.1      knitr_1.40          pillar_1.8.1       
 [70] boot_1.3-22         estimability_1.4.1  effectsize_0.7.0.5 
 [73] future.apply_1.9.1  codetools_0.2-16    glue_1.6.2         
 [76] evaluate_0.16       modelr_0.1.9        data.table_1.14.4  
 [79] renv_0.12.0         nloptr_1.2.2.2      vctrs_0.5.0        
 [82] httpuv_1.6.1        gtable_0.3.1        future_1.28.0      
 [85] assertthat_0.2.1    datawizard_0.6.2    cachem_1.0.4       
 [88] xfun_0.33           mime_0.12           xtable_1.8-4       
 [91] broom_1.0.1         pracma_2.4.2        later_1.2.0        
 [94] ragg_1.2.4          viridisLite_0.4.1   signal_0.7-6       
 [97] tibble_3.1.8        statmod_1.4.37      globals_0.16.1     
[100] ellipsis_0.3.2