# ============================================================================= # REPRODUCIBLE R SCRIPT FOR THE JoVE FIGURES (Figures 3, 4, 5) # "Three Behavioral Paradigms to Assess Fear and Reward Learning in Relation # to Adverse Childhood Experiences" # ----------------------------------------------------------------------------- # Author : Maria Consuelo San Martin Miranda # Date : May 2026 # Project: Fondecyt Iniciacion # ----------------------------------------------------------------------------- # WHAT THIS SCRIPT DOES # --------------------- # Generates Figures 3, 4 and 5 directly from the raw Excel files, in PNG # and TIFF at 300 dpi. Every panel is rendered as a BOXPLOT (median, IQR # and individual outliers) following modern data-visualization guidance # (Weissgerber et al., 2015, PLoS Biology). The boxplot does not assume # normality and exposes the full distribution -- including the floor # effects in the CS- rating boxes -- which a bar+SE plot would hide. # # The choice of visualization is independent of the choice of inference # test. The reported tests are unchanged from `paper_analyses_R.R`: # * F3A, F4A SCR : paired-samples t-test + Bonferroni (parametric) # * F3B, F3C, F4B, F4C : paired-samples t-test + Bonferroni (parametric) # * F5A SCR : Wilcoxon signed-rank (non-parametric) # * F5B, F5C ratings : paired-samples t-test (parametric) # Brackets are shown only for SIGNIFICANT pairs (Bonferroni-corrected, # threshold .05); non-significant pairs are not drawn -- the omnibus # results are reported in the Methods/Results section. # # DATA REQUIRED # ------------- # Base_Completa_Exp_1_08_Abril.xlsx # Base_Completa_Exp_2_08_Abril.xlsx # # HOW TO RUN # ---------- # "C:/Program Files/R/R-4.5.0/bin/Rscript.exe" paper_figures_R.R # or from RStudio: # setwd("path/to/R_Script") # source("paper_figures_R.R") # # OUTPUT (saved next to the script in ../figures_v2/) # --------------------------------------------------- # Figure_3_Pavlovian_Exp1.png (and .tiff) # Figure_4_Operant_Exp1.png (and .tiff) # Figure_5_Exp2_Phases.png (and .tiff) # ============================================================================= suppressPackageStartupMessages({ pkgs <- c("readxl", "dplyr", "tidyr", "ggplot2", "ggsignif", "patchwork") for (p in pkgs) { if (!requireNamespace(p, quietly = TRUE)) install.packages(p, repos = "http://cran.us.r-project.org") library(p, character.only = TRUE) } }) # Locate raw Excel required_files <- c("Base_Completa_Exp_1_08_Abril.xlsx", "Base_Completa_Exp_2_08_Abril.xlsx") candidates <- c(normalizePath(".", winslash = "/", mustWork = FALSE), normalizePath("..", winslash = "/", mustWork = FALSE), normalizePath("../..", winslash = "/", mustWork = FALSE), normalizePath("../../..", winslash = "/", mustWork = FALSE), "C:/Users/mconsuelo/Desktop/dato_fondecytIniciacion") DATA_DIR <- NULL for (path in candidates) { if (is.na(path) || !nzchar(path) || !dir.exists(path)) next if (all(file.exists(file.path(path, required_files)))) { DATA_DIR <- path; break } } if (is.null(DATA_DIR)) stop("Could not locate Excel files.") cat("Working directory:", DATA_DIR, "\n") # Locate script directory and figure output dir get_script_dir <- function() { args <- commandArgs(trailingOnly = FALSE) m <- regmatches(args, regexpr("(?<=--file=).+", args, perl = TRUE)) if (length(m) && nzchar(m[1])) return(normalizePath(dirname(m[1]), winslash = "/", mustWork = FALSE)) f <- tryCatch(sys.frame(1)$ofile, error = function(e) NULL) if (!is.null(f) && nzchar(f)) return(normalizePath(dirname(f), winslash = "/", mustWork = FALSE)) normalizePath(".", winslash = "/", mustWork = FALSE) } SCRIPT_DIR <- get_script_dir() FIG_DIR <- file.path(dirname(SCRIPT_DIR), "figures_v2") if (!dir.exists(FIG_DIR)) dir.create(FIG_DIR, recursive = TRUE) cat("Figures will be saved in:", FIG_DIR, "\n") # ---- Load data -------------------------------------------------------------- df1 <- read_excel(file.path(DATA_DIR, "Base_Completa_Exp_1_08_Abril.xlsx")) df2 <- read_excel(file.path(DATA_DIR, "Base_Completa_Exp_2_08_Abril.xlsx")) names(df1) <- iconv(names(df1), to = "ASCII//TRANSLIT") names(df2) <- iconv(names(df2), to = "ASCII//TRANSLIT") SCR_RAW <- "CDA.SCR [muS]" df1$CDA.SCR_log <- log(pmax(df1[[SCR_RAW]], 0) + 1) df2$CDA.SCR_log <- log(pmax(df2[[SCR_RAW]], 0) + 1) # ---- Helpers --------------------------------------------------------------- fase_wide <- function(d, fases, stims, measure) { d %>% filter(Fase %in% fases, Estimulo %in% stims) %>% mutate(VAL = suppressWarnings(as.numeric(.data[[measure]]))) %>% group_by(Sujeto, Estimulo) %>% summarise(VAL = mean(VAL, na.rm = TRUE), .groups = "drop") %>% pivot_wider(names_from = Estimulo, values_from = VAL) %>% drop_na() } p_to_label <- function(p) { if (is.na(p)) return("") if (p < .001) return("***") if (p < .01) return("**") if (p < .05) return("*") "" } pair_tests <- function(wide_df, stims, x_labels = stims, test = "t") { pairs <- combn(stims, 2, simplify = FALSE) rows <- lapply(pairs, function(pp) { a <- wide_df[[pp[1]]]; b <- wide_df[[pp[2]]] if (sum(!is.na(a - b)) < 3) return(NULL) if (test == "wilcoxon") { tt <- suppressWarnings(wilcox.test(a, b, paired = TRUE)) } else { tt <- t.test(a, b, paired = TRUE) } data.frame( A_lab = x_labels[match(pp[1], stims)], B_lab = x_labels[match(pp[2], stims)], Phase = NA_character_, p_unc = tt$p.value, stringsAsFactors = FALSE ) }) out <- bind_rows(rows) out$p_bonf <- pmin(1, out$p_unc * nrow(out)) out$label <- vapply(out$p_bonf, p_to_label, character(1)) out[nzchar(out$label), , drop = FALSE] } theme_apa <- theme_classic(base_size = 11, base_family = "serif") + theme(axis.title.x = element_blank(), axis.text.x = element_text(color = "black", size = 10), axis.text.y = element_text(color = "black", size = 10), axis.line = element_line(color = "black", linewidth = 0.5), axis.ticks = element_line(color = "black"), panel.grid = element_blank(), plot.tag = element_text(size = 14, face = "bold"), plot.tag.position = c(0.02, 0.98), legend.position = "none", plot.margin = margin(8, 8, 8, 8)) # Boxplot for a single panel with K stimuli on the x-axis stim_boxplot <- function(wide_df, stims, x_labels, pairs_df, ylab, panel_tag, fill_colors, ylim = NULL) { long <- wide_df %>% select(all_of(stims)) %>% pivot_longer(everything(), names_to = "stim_id", values_to = "Value") %>% mutate(Stimulus = factor(stim_id, levels = stims, labels = x_labels)) data_top <- max(long$Value, na.rm = TRUE) if (is.null(ylim)) { span <- diff(range(long$Value, na.rm = TRUE)) span <- if (!is.finite(span) || span == 0) abs(data_top) + 1 else span y_max <- data_top + span * (0.10 + 0.12 * max(nrow(pairs_df), 1)) y_lo <- min(0, min(long$Value, na.rm = TRUE)) } else { y_lo <- ylim[1]; y_max <- ylim[2] span <- y_max - y_lo } p <- ggplot(long, aes(x = Stimulus, y = Value, fill = Stimulus)) + geom_boxplot(width = 0.6, outlier.size = 1.0, color = "black", linewidth = 0.4) + scale_fill_manual(values = setNames(fill_colors, x_labels)) + scale_y_continuous(expand = expansion(mult = c(0.02, 0)), limits = c(y_lo, y_max)) + labs(y = ylab, tag = panel_tag) + theme_apa if (nrow(pairs_df) > 0) { # Order by horizontal span so closer pairs are drawn lower ord <- order(abs(match(pairs_df$A_lab, x_labels) - match(pairs_df$B_lab, x_labels))) pairs_df <- pairs_df[ord, , drop = FALSE] base_y <- data_top + span * 0.06 step_y <- span * 0.10 y_pos <- base_y + (seq_len(nrow(pairs_df)) - 1) * step_y p <- p + ggsignif::geom_signif( annotations = pairs_df$label, y_position = y_pos, xmin = match(pairs_df$A_lab, x_labels), xmax = match(pairs_df$B_lab, x_labels), tip_length = 0.015, textsize = 3.6, vjust = -0.2, family = "serif" ) } p } # Phase x Stimulus boxplot (for F5 with 3 phases) phase_boxplot <- function(long_df, pairs_df, phase_levels, stim_levels, ylab, panel_tag, fill_colors, ylim = NULL, show_legend = FALSE) { long_df$Phase <- factor(long_df$Phase, levels = phase_levels) long_df$Stimulus <- factor(long_df$Stimulus, levels = stim_levels) data_top <- max(long_df$Value, na.rm = TRUE) if (is.null(ylim)) { span <- diff(range(long_df$Value, na.rm = TRUE)) span <- if (!is.finite(span) || span == 0) abs(data_top) + 1 else span y_max <- data_top + span * 0.20 y_lo <- min(0, min(long_df$Value, na.rm = TRUE)) } else { y_lo <- ylim[1]; y_max <- ylim[2] span <- y_max - y_lo } dodge <- position_dodge(width = 0.8) p <- ggplot(long_df, aes(x = Phase, y = Value, fill = Stimulus)) + geom_boxplot(width = 0.6, outlier.size = 1.0, color = "black", linewidth = 0.4, position = dodge) + scale_fill_manual(values = setNames(fill_colors, stim_levels)) + scale_y_continuous(expand = expansion(mult = c(0.02, 0)), limits = c(y_lo, y_max)) + labs(y = ylab, tag = panel_tag) + theme_apa if (show_legend) { p <- p + theme(legend.position = c(0.85, 0.95), legend.justification = c(1, 1), legend.title = element_blank(), legend.text = element_text(family = "serif", size = 10), legend.background = element_blank(), legend.key.height = unit(0.5, "lines")) } if (nrow(pairs_df) > 0) { bar_w <- 0.8; offset <- bar_w / 4 pos <- match(pairs_df$Phase, phase_levels) box_top <- vapply(seq_len(nrow(pairs_df)), function(i) { sub <- long_df[long_df$Phase == pairs_df$Phase[i], ] max(sub$Value, na.rm = TRUE) }, numeric(1)) p <- p + ggsignif::geom_signif( annotations = pairs_df$label, y_position = box_top + span * 0.05, xmin = pos - offset, xmax = pos + offset, tip_length = 0.012, textsize = 3.6, vjust = -0.2, family = "serif" ) } p } # ---- Build wide tables ----------------------------------------------------- stim1 <- c("EC+ Evitable","EC+ No Evitable","EC-") labels1 <- c("Avoidable\nCS+","Unavoidable\nCS+","CS-") stim2 <- c("EC+","EC-") labels2 <- c("CS+","CS-") pav1_w <- fase_wide(df1, "Pavloviano", stim1, "CDA.SCR_log") op1_w <- fase_wide(df1, "Operante", stim1, "CDA.SCR_log") pav_exp_w <- fase_wide(df1, c("Pavlov Expectativa 1","Pavlov Expectativa 2","Pavlov Expectativa 3"), stim1, "Respuestas_Preguntas") pav_ame_w <- fase_wide(df1, c("Pavlov Amenaza 1","Pavlov Amenaza 2","Pavlov Amenaza 3"), stim1, "Respuestas_Preguntas") op_exp_w <- fase_wide(df1, c("Operante Expectativa 1","Operante Expectativa 2","Operante Expectativa 3"), stim1, "Respuestas_Preguntas") op_ame_w <- fase_wide(df1, c("Operante Amenaza 1","Operante Amenaza 2","Operante Amenaza 3"), stim1, "Respuestas_Preguntas") op_aliv_w <- fase_wide(df1, c("Operante Alivio 1","Operante Alivio 2","Operante Alivio 3"), stim1, "Respuestas_Preguntas") bot_w <- df1 %>% filter(Fase == "Operante", Estimulo %in% stim1) %>% mutate(Espacios = suppressWarnings(as.numeric(`# Espacios`))) %>% group_by(Sujeto, Estimulo) %>% summarise(Espacios = mean(Espacios, na.rm = TRUE), .groups = "drop") %>% pivot_wider(names_from = Estimulo, values_from = Espacios) %>% drop_na() adq_w <- fase_wide(df2, "Pavloviano", stim2, "CDA.SCR_log") ext_w <- fase_wide(df2, "Extincion", stim2, "CDA.SCR_log") readq_w <- fase_wide(df2, "ReAdquisicion", "EC+", "CDA.SCR_log") rat_adq_exp <- fase_wide(df2, c("Pavlov Expectativa 1","Pavlov Expectativa 2"), stim2, "Respuestas_Preguntas") rat_ext_exp <- fase_wide(df2, c("Extincion Expectativa 1","Extincion Expectativa 2"), stim2, "Respuestas_Preguntas") rat_re_exp <- fase_wide(df2, c("ReAdquisicion Expectativa 1","ReAdquisicion Expectativa 2"), stim2, "Respuestas_Preguntas") rat_adq_ame <- fase_wide(df2, c("Pavlov Amenaza 1","Pavlov Amenaza 2"), stim2, "Respuestas_Preguntas") rat_ext_ame <- fase_wide(df2, c("Extincion Amenaza 1","Extincion Amenaza 2"), stim2, "Respuestas_Preguntas") rat_re_ame <- fase_wide(df2, c("ReAdquisicion Amenaza 1","ReAdquisicion Amenaza 2"), stim2, "Respuestas_Preguntas") # ---- Figure 3 (boxplot version) -------------------------------------------- COL3 <- c("#d97706","#c0392b","#2c5282") f3a <- stim_boxplot(pav1_w, stim1, labels1, pair_tests(pav1_w, stim1, labels1, "t"), ylab = expression("log SCR (log["*mu*"S + 1])"), panel_tag = "A", fill_colors = COL3) f3b <- stim_boxplot(pav_exp_w, stim1, labels1, pair_tests(pav_exp_w, stim1, labels1, "t"), ylab = "Expectancy rating (VAS 0-100)", panel_tag = "B", fill_colors = COL3, ylim = c(0, 130)) f3c <- stim_boxplot(pav_ame_w, stim1, labels1, pair_tests(pav_ame_w, stim1, labels1, "t"), ylab = "Threat rating (VAS 0-100)", panel_tag = "C", fill_colors = COL3, ylim = c(0, 130)) fig3 <- f3a + f3b + f3c + patchwork::plot_layout(ncol = 3) ggsave(file.path(FIG_DIR, "Figure_3_Pavlovian_Exp1.png"), fig3, width = 12, height = 4.2, dpi = 300, bg = "white") ggsave(file.path(FIG_DIR, "Figure_3_Pavlovian_Exp1.tiff"), fig3, width = 12, height = 4.2, dpi = 300, bg = "white", compression = "lzw") cat(" saved Figure_3 (boxplot)\n") # ---- Figure 4 (boxplot version) -------------------------------------------- f4a <- stim_boxplot(op1_w, stim1, labels1, pair_tests(op1_w, stim1, labels1, "t"), ylab = expression("log SCR (log["*mu*"S + 1])"), panel_tag = "A", fill_colors = COL3) f4b <- stim_boxplot(bot_w, stim1, labels1, pair_tests(bot_w, stim1, labels1, "t"), ylab = "Number of button presses", panel_tag = "B", fill_colors = COL3) f4c <- stim_boxplot(op_aliv_w, stim1, labels1, pair_tests(op_aliv_w, stim1, labels1, "t"), ylab = "Relief rating (VAS 0-100)", panel_tag = "C", fill_colors = COL3, ylim = c(0, 130)) fig4 <- f4a + f4b + f4c + patchwork::plot_layout(ncol = 3) ggsave(file.path(FIG_DIR, "Figure_4_Operant_Exp1.png"), fig4, width = 12, height = 4.2, dpi = 300, bg = "white") ggsave(file.path(FIG_DIR, "Figure_4_Operant_Exp1.tiff"), fig4, width = 12, height = 4.2, dpi = 300, bg = "white", compression = "lzw") cat(" saved Figure_4 (boxplot)\n") # ---- Figure 5 (boxplot version) -------------------------------------------- COL2 <- c("#c0392b","#2c5282") phase_lvls <- c("Acquisition","Extinction","Reacquisition") to_long <- function(w, phase, stims) { out <- list() for (s in stims) if (s %in% names(w)) out[[length(out)+1]] <- data.frame(Phase=phase, Stimulus=s, Value=w[[s]]) bind_rows(out) } f5a_long <- bind_rows( to_long(adq_w, "Acquisition", stim2), to_long(ext_w, "Extinction", stim2), to_long(readq_w, "Reacquisition", "EC+") ) f5a_long$Stimulus <- factor(f5a_long$Stimulus, levels = stim2, labels = labels2) f5a_pairs <- bind_rows( pair_tests(adq_w, stim2, labels2, "wilcoxon") %>% mutate(Phase = "Acquisition"), pair_tests(ext_w, stim2, labels2, "wilcoxon") %>% mutate(Phase = "Extinction") ) f5a <- phase_boxplot(f5a_long, f5a_pairs, phase_lvls, labels2, ylab = expression("log SCR (log["*mu*"S + 1])"), panel_tag = "A", fill_colors = COL2, show_legend = TRUE) f5b_long <- bind_rows( to_long(rat_adq_exp, "Acquisition", stim2), to_long(rat_ext_exp, "Extinction", stim2), to_long(rat_re_exp, "Reacquisition", stim2) ) f5b_long$Stimulus <- factor(f5b_long$Stimulus, levels = stim2, labels = labels2) f5b_pairs <- bind_rows( pair_tests(rat_adq_exp, stim2, labels2, "t") %>% mutate(Phase = "Acquisition"), pair_tests(rat_ext_exp, stim2, labels2, "t") %>% mutate(Phase = "Extinction"), pair_tests(rat_re_exp, stim2, labels2, "t") %>% mutate(Phase = "Reacquisition") ) f5b <- phase_boxplot(f5b_long, f5b_pairs, phase_lvls, labels2, ylab = "Expectancy rating (VAS 0-100)", panel_tag = "B", fill_colors = COL2, ylim = c(0, 110)) f5c_long <- bind_rows( to_long(rat_adq_ame, "Acquisition", stim2), to_long(rat_ext_ame, "Extinction", stim2), to_long(rat_re_ame, "Reacquisition", stim2) ) f5c_long$Stimulus <- factor(f5c_long$Stimulus, levels = stim2, labels = labels2) f5c_pairs <- bind_rows( pair_tests(rat_adq_ame, stim2, labels2, "t") %>% mutate(Phase = "Acquisition"), pair_tests(rat_ext_ame, stim2, labels2, "t") %>% mutate(Phase = "Extinction"), pair_tests(rat_re_ame, stim2, labels2, "t") %>% mutate(Phase = "Reacquisition") ) f5c <- phase_boxplot(f5c_long, f5c_pairs, phase_lvls, labels2, ylab = "Threat rating (VAS 0-100)", panel_tag = "C", fill_colors = COL2, ylim = c(0, 110)) fig5 <- f5a + f5b + f5c + patchwork::plot_layout(ncol = 3) ggsave(file.path(FIG_DIR, "Figure_5_Exp2_Phases.png"), fig5, width = 13, height = 4.2, dpi = 300, bg = "white") ggsave(file.path(FIG_DIR, "Figure_5_Exp2_Phases.tiff"), fig5, width = 13, height = 4.2, dpi = 300, bg = "white", compression = "lzw") cat(" saved Figure_5 (boxplot)\n") cat("\nDone.\n")