In epongpipat/ppi: Psychophysiological Interaction (PPI)

R Libraries

packages <- c("ppi", "tidyverse", "broom", "car", "performance", "bootPermBroom", "furrr")
xfun::pkg_attach(packages, message = F, install = T)

Data Wrangling

# mri parameters
tr <- 3
n_volumes <- 105

psy_events <- readr::read_tsv(url("https://openneuro.org/crn/datasets/ds000171/snapshots/00001/files/sub-control01:func:sub-control01_task-music_run-1_events.tsv")) %>%
  mutate(trial_type = as.factor(trial_type))

psy_contrast_table <- cbind(stimulus_vs_response = c(1, 1, -3, 1)/4,
                                          music_vs_tones = c(1, 1, 0, -2)/3,
                                          positive_music_vs_negative_music = c(-1, 1, 0, 0)/2)

upsample_factor <- 16
hrf <- create_hrf_afni("spmg1", tr, upsample_factor)
detrend_factor <- 2

phys_left_parietal <- "sub-control01_task-music_run-1_bold_space-subj_vox-32-24-38.csv"

phys_target_file <- "sub-control01_task-music_run-1_bold_space-subj_vox-48-24-36.csv"

phys_seed <- read_csv(phys_left_parietal, col_names = "seed")
psy_unlabeled_trial_type <- "response"

phys_target <- read_csv(phys_target_file, col_names = "target")

data_wrangled_left <- data_wrangling(psy_events, psy_unlabeled_trial_type, psy_contrast_table, phys_seed, detrend_factor, hrf, tr, n_volumes, upsample_factor)
data_wrangled_right <- data_wrangling(psy_events, psy_unlabeled_trial_type, psy_contrast_table, phys_target, detrend_factor, hrf, tr, n_volumes, upsample_factor)

save_data_wrangling(data_wrangled_left)

create_data_wrangling_report(data_wrangled_left)

Analysis

ROI-to-ROI

y_1d <- data_wrangled_right$phys_var$input
y_2d <- as.matrix(cbind(data_wrangled_right$phys_var$input, data_wrangled_left$phys_var$input))
y_list <- list(data_wrangled_right$phys_var$input, data_wrangled_left$phys_var$input)

X_2d <- data_wrangled_right$design_matrix
X_list <- list(data_wrangled_right$design_matrix %>% mutate(run = c(rep(1, 50), rep(2, 55))), data_wrangled_left$design_matrix %>% mutate(run = c(rep(1, 50), rep(2, 55))))

model_l1 <- model_glm_roi2roi(y_list, X_2d)
model_l1_2 <- model_glm_roi2roi(y_list, X_2d, "psy_positive_music_vs_negative_music")

x_names <- colnames(X_list[[1]])
model_lmer_1 <- model_lmer_roi2roi(y_list, X_list, glue("{x_names} + ( {x_names} | run )"))

model_l1 %>%
  unnest(vif_tol)

# extract estimate for each predictor
tidy_4d <- array(NA, c(2, 2, 8, 6))
performance_3d <- array(NA, c(2, 2, 5))
residual_3d <- array(NA, c(2, 2, n_volumes))
tolerance_3d <- array(NA, c(2, 2, 7))
n_seed <- length(model)
idx <- 1
for (i in 1:n_seed) {
  n_target <- length(model[[i]])
  for (j in 1:n_target) {
    temp_model <- model[[i]][[j]]
    tidy_4d[i,j,,] <- tidy(temp_model) %>% tidy_lm_add_r_squared(., 105) %>% select(-term) %>% as.matrix()
    performance_3d[i,j,] <- model_performance(temp_model) %>% as.matrix()
    residual_3d[i,j,] <- temp_model$residuals
    tolerance_3d[i,j,] <- 1/vif(temp_model)
    idx <- idx + 1
  }
}


# create 3d array of seed by target by statistic ----
# 1. omnibus model rsquared
# 2. omnibus model adjusted rsquared
# 3. omnibus model rmse
# 4  estimate
# 5. se
# 6. t-stat
# 7. p-value
# 8. r-squared
# 9. adjusted r-squared
model_summary <- list()
model_summary$model_rsq <- performance_3d[,,3] %>% as.data.frame()
model_summary$model_rsq_adj <- performance_3d[,,4] %>% as.data.frame()
model_summary$model_rmse <- performance_3d[,,5] %>% as.data.frame()

for (i in 1:8) {
  model_summary$estimates[[i]] <- tidy_4d[,,i,1] %>% as.data.frame()
  model_summary$se[[i]] <- tidy_4d[,,i,2] %>% as.data.frame()
  model_summary$t_stat[[i]] <- tidy_4d[,,i,3] %>% as.data.frame()
  model_summary$p_value[[i]] <- tidy_4d[,,i,4] %>% as.data.frame()
  model_summary$rsq[[i]] <- tidy_4d[,,i,5] %>% as.data.frame()
  model_summary$rsq_adj[[i]] <- tidy_4d[,,i,6] %>% as.data.frame()
}

Visualization

ROI-to-ROI

labels <- c("parietal_l", "parietal_r")

temp_data <- as.data.frame(model_summary$model_rsq)
colnames(temp_data) <- labels
temp_data <- temp_data %>% 
  mutate(seed_roi = labels) %>%
  gather(., target_roi, value, -seed_roi)

ggplot(temp_data, aes(seed_roi, target_roi, fill = value)) +
  geom_raster() +
  theme_minimal() +
  scale_fill_distiller(palette = "Greys", direction = 1)