Goal: Paired t-test on mean univariate activation (Cued-Uncued)

Preprocessing

arousal.bold.uni <- read_csv("task-arousal_maskedMeanEff.csv")
## New names:
## Rows: 1872 Columns: 4
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (2): zmap, mask dbl (2): ...1, meanEffEstimate
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
filtered_df <- arousal.bold.uni %>%
  filter(str_detect(zmap, "contrast-Cued_zmap") | str_detect(zmap, "contrast-Uncued_zmap")) %>%
  mutate(subject = str_extract(zmap, "(?<=sub-)(.*?)(?=_task)"),
    pretty_mask = map_chr(mask, ~basename(.x)),
    contrast = str_extract(zmap, "(?<=contrast-)(.*?)(?=_zmap)"))

Paired t-test

library(dplyr)
library(tidyr)
library(purrr)
library(broom)
## Warning: package 'broom' was built under R version 4.2.3
# Spread the dataframe so that each subject has one row with both Cued and Uncued values.
spread_df <- filtered_df %>% select(-c(`...1`,zmap,mask)) %>%
  spread(key = contrast, value = meanEffEstimate) 

# Perform a paired t-test for each mask and tidy the results.
t_test_results <- spread_df %>%
  group_by(pretty_mask) %>%
  do(tidy(t.test(.$Cued, .$Uncued, paired = TRUE)))

# Focus on predefined ROIs
t_test_results.filtered <-t_test_results %>% ungroup %>% filter(pretty_mask %in% c("resampled_vinsula_mask.nii.gz","resampled_vsgcv2_mask.nii.gz","resampled_vamygdala_mask.nii.gz","resampled_vorbitofrontal_mask.nii.gz"))

kable(t_test_results.filtered)
pretty_mask estimate statistic p.value parameter conf.low conf.high method alternative
resampled_vamygdala_mask.nii.gz -0.1699035 -0.8858532 0.3880560 17 -0.5745587 0.2347517 Paired t-test two.sided
resampled_vinsula_mask.nii.gz -0.4946816 -2.7680838 0.0131592 17 -0.8717248 -0.1176385 Paired t-test two.sided
resampled_vorbitofrontal_mask.nii.gz -0.2555129 -2.5601917 0.0202798 17 -0.4660772 -0.0449485 Paired t-test two.sided
resampled_vsgcv2_mask.nii.gz -0.3121361 -2.7674111 0.0131778 17 -0.5501020 -0.0741701 Paired t-test two.sided

looks like insula, ofc and sgACC is sig at p<0.05

Visualize paired t-test of mean bold per each ROI

filtered_df.2 <- filtered_df %>% select(-c(`...1`,zmap,mask)) %>%
  mutate(
    pretty_mask_name = case_when(
      str_detect(pretty_mask, "insula") ~ "Insula",
      str_detect(pretty_mask, "orbitofrontal") ~ "OFC",
      str_detect(pretty_mask, "sgc") ~ "sgACC",
      str_detect(pretty_mask, "amygdala") ~ "Amygdala",
      TRUE ~ "Other"
    )
  )

filtered_df.2 <- filtered_df.2 %>% filter(pretty_mask %in% c("resampled_vinsula_mask.nii.gz","resampled_vsgcv2_mask.nii.gz","resampled_vamygdala_mask.nii.gz","resampled_vorbitofrontal_mask.nii.gz"))

# Loading necessary libraries
library(ggstatsplot)
## Warning: package 'ggstatsplot' was built under R version 4.2.3
## You can cite this package as:
##      Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
##      Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167
# Using grouped_ggwithinstats to create the plots
grouped_ggwithinstats(
  data = filtered_df.2,
  x = contrast,
  y = meanEffEstimate,
  id.var = subject,
  grouping.var = pretty_mask_name,
  bf.message=FALSE,
  results.subtitle=FALSE)

with dotplot

# Ensure the contrast column is a factor with the correct order
filtered_df.2$contrast <- factor(filtered_df.2$contrast, levels = c("Uncued", "Cued"))

# Using ggplot to create the plot
# Using ggplot to create the plot
# Loading necessary library
library(ggplot2)

# Using ggplot to create the plot with stat_summary
plot <- ggplot(data = filtered_df.2, aes(x = contrast, y = meanEffEstimate, color = contrast)) +
  
  # Point to represent the mean
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3,
    position = position_dodge(width = 0.7)
  ) +

  # Error bars for standard errors
  stat_summary(
    fun.data = function(y) {
      y_mean <- mean(y)
      y_se <- sd(y) / sqrt(length(y))
      return(list(y = y_mean, ymin = y_mean - y_se, ymax = y_mean + y_se))
    },
    geom = "errorbar",
    width = 0.2,
    position = position_dodge(width = 0.7)
  ) +
  
  # Separate plots for each pretty_mask_name
  facet_wrap(~ pretty_mask_name, scales = "free_y") + 
  labs(title = "Within-Subject Mean and Standard Error", y = "Mean Effect Estimate") + 
  theme_apa() +
  theme(legend.position = "bottom")

print(plot)

lets see with individual points

plot <- ggplot(data = filtered_df.2, aes(x = contrast, y = meanEffEstimate, color = contrast, group = pretty_mask_name)) +
  
  # Line to connect the means
  stat_summary(
    fun.y = mean,
    geom = "line",
    aes(group = pretty_mask_name),
    color = "black",
    position = position_dodge(width = 0.7)
  ) +

  # Add individual data points with some transparency
  geom_point(aes(color = contrast), alpha = 0.5, position = position_dodge(width = 0.7), size = 1) +

  # Point to represent the mean
  stat_summary(
    fun = mean,
    shape = 15,
    geom = "point",
    size = 3,
    position = position_dodge(width = 0.7)
  ) +

  # Error bars for standard errors
  stat_summary(
    fun.data = function(y) {
      y_mean <- mean(y)
      y_se <- sd(y) / sqrt(length(y))
      return(list(y = y_mean, ymin = y_mean - y_se, ymax = y_mean + y_se))
    },
    geom = "errorbar",
    width = 0.2,
    position = position_dodge(width = 0.7)
  ) +
  
  # Separate plots for each pretty_mask_name
  facet_wrap(~ pretty_mask_name, scales = "free_y") + 
  labs(title = "Within-Subject Mean and Standard Error", y = "Mean Effect Estimate") + 
  theme_apa() +
  theme(legend.position = "bottom")
## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(plot)