Frontoparietal control-default mode connectivity predicts TMS effects on cognitive control.

Code ↔ Paper

5 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 5 matches

[1] § Methods › Behavioral data processing and modeling ↔ 4_DDM_Validation.Rmd, lines 148–218 · score 0.67 · posterior predictive checks, fast dm, Validations, DDM, model
[2] § Methods › Behavioral data processing and modeling ↔ 3_DDMAnalysis.Rmd, lines 421–500 · score 0.56 · Exploratory factor, paran, parallel, psych, factor model, components
[3] § Methods › Behavioral data processing and modeling ↔ 5_Alternative_Models_Raw.Rmd, lines 473–552 · score 0.56 · Exploratory factor, paran, parallel, psych, factor model, components
[4] § Methods › Statistical analysis ↔ 7_Alternative_Model_MRI_SCRUBBED_DDM_Analysis.Rmd, lines 17–88 · score 0.56 · framewise displacement, FD, MRI, head, models
[5] § Results › Validation of factor structure ↔ 5_Alternative_Models_Raw.Rmd, lines 621–633 · score 0.52 · chi square, full model, CFA, CFI, RMSEA, SRMR

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R Markdown · 2,660 lines · 108 KB · no license · 2 matches

---
title: "Raw_Data_Analysis"
output: html_document
date: "2024-02-19"
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r}
library(tidyverse)
library(dplyr)
library(ggplot2)
```
# Factor Analysis for RT
## Load Data
**Note: The Response Time (RT) analysis is followed by an identical automated pipeline for Accuracy (ACC) starting near line ~1420.**
```{r}
# ---------------- Navon ----------------
behavioral_data_navon <- read.csv("FinalData/Navon_Behavioral_LongFormat_PRISM_final_x4_z_day.csv")
participants_without_FPCN_B_navon <- behavioral_data_navon %>%
group_by(Subj) %>%
summarise(has_FPCN_B = any(Stimulation_Site == "FPCN-B")) %>%
filter(!has_FPCN_B) %>%
pull(Subj)
cleaned_data_navon <- behavioral_data_navon %>%
filter(!Subj %in% participants_without_FPCN_B_navon) %>%
filter(RT >= 0.200 & RT <= 1.5) %>%
mutate(Stimulation_Site = factor(Stimulation_Site, levels = c("Vertex", "FPCN-B", "DAN")))
navon_raw_data <- cleaned_data_navon %>%
group_by(Subj, Stimulation_Site, Timepoint, Task_Low_High, Days) %>%
summarize(rt = mean(RT[correct == 1], na.rm=TRUE), acc = mean(correct, na.rm=TRUE), .groups = "drop")
# ---------------- Stroop ----------------
behavioral_data_stroop <- read.csv("FinalData/Stroop_Behavioral_LongFormat_PRISM_final_x3_z_day.csv")
participants_without_FPCN_B_stroop <- behavioral_data_stroop %>%
group_by(Subj) %>%
summarise(has_FPCN_B = any(Stimulation_Site == "FPCN-B")) %>%
filter(!has_FPCN_B) %>%
pull(Subj)
cleaned_data_stroop <- behavioral_data_stroop %>%
filter(!Subj %in% participants_without_FPCN_B_stroop) %>%
filter(RT >= 0.200 & RT <= 1.5) %>%
mutate(Stimulation_Site = factor(Stimulation_Site, levels = c("Vertex", "DAN", "FPCN-B")))
stroop_raw_data <- cleaned_data_stroop %>%
group_by(Subj, Stimulation_Site, Timepoint, Task_Low_High, Days) %>%
summarize(rt = mean(RT[correct == 1], na.rm=TRUE), acc = mean(correct, na.rm=TRUE), .groups = "drop")
# ---------------- NBack ----------------
behavioral_data_nback <- read.csv("FinalData/n-back_exp_results.csv")
participants_without_FPCN_B_nback <- behavioral_data_nback %>%
group_by(Subj) %>%
summarise(has_FPCN_B = any(Stimulation_Site == "FPCNB")) %>%
filter(!has_FPCN_B) %>%
pull(Subj)
cleaned_data_nback <- behavioral_data_nback %>%
filter(!Subj %in% participants_without_FPCN_B_nback) %>%
filter(RT >= 0.200 & RT <= 1.5)
nback_raw_data <- cleaned_data_nback %>%
mutate(Task_Low_High = case_when(
Condition == 0 ~ "low",
Condition == 1 ~ "medium",
Condition == 2 ~ "high"
),
Stimulation_Site = case_when(
Stimulation_Site == "vertex" ~ "Vertex",
Stimulation_Site == "FPCNB" ~ "FPCN-B",
TRUE ~ Stimulation_Site
)) %>%
group_by(Subj, Stimulation_Site, Timepoint, Task_Low_High, Days) %>%
summarize(rt = mean(RT[correct == 1], na.rm=TRUE), acc = mean(correct, na.rm=TRUE), .groups = "drop")
```
```{r}
# Filter out participants with Days z-score > 3
navon_raw_data <- navon_raw_data %>%
filter(abs((Days - mean(Days, na.rm = TRUE)) / sd(Days, na.rm = TRUE)) <= 2.5)
nback_raw_data <- nback_raw_data %>%
filter(abs((Days - mean(Days, na.rm = TRUE)) / sd(Days, na.rm = TRUE)) <= 2.5)
stroop_raw_data <- stroop_raw_data %>%
filter(abs((Days - mean(Days, na.rm = TRUE)) / sd(Days, na.rm = TRUE)) <= 2.5)
```
Quick checks on distribution of days for counterbalancing, we want to make sure that there's nothing qualitatively different about DAN
```{r}
summary_stats <- stroop_raw_data %>%
group_by(Stimulation_Site) %>%
summarise(
mean_Days = mean(Days, na.rm = TRUE),
sd_Days = sd(Days, na.rm = TRUE),
median_Days = median(Days, na.rm = TRUE)
)
# Plot histograms for each statistic in separate facets
ggplot(stroop_raw_data, aes(x = Days)) +
geom_histogram(binwidth = 1, fill = "lightblue", color = "black") +
facet_wrap(~ Stimulation_Site, scales = "free_y") +
theme_minimal() +
labs(
title = "NBack Histograms of Days",
x = "Days",
y = "Count"
)
```
```{r}
navon_wide <- navon_raw_data %>%
pivot_wider(
id_cols = c(Subj, Stimulation_Site, Timepoint),
names_from = Task_Low_High,
values_from = c(rt, acc),
names_prefix = ""
)
stroop_wide <- stroop_raw_data %>%
pivot_wider(
id_cols = c(Subj, Stimulation_Site, Timepoint),
names_from = Task_Low_High,
values_from = c(rt, acc),
names_prefix = ""
)
nback_wide <- nback_raw_data %>%
pivot_wider(
id_cols = c(Subj, Stimulation_Site, Timepoint),
names_from = Task_Low_High,
values_from = c(rt, acc),
names_prefix = ""
)
navon_pre <- subset(navon_wide, Timepoint == "pre")
stroop_pre <- subset(stroop_wide, Timepoint == "pre")
nback_pre <- subset(nback_wide, Timepoint == "pre")
navon_post <- subset(navon_wide, Timepoint == "post")
stroop_post <- subset(stroop_wide, Timepoint == "post")
nback_post <- subset(nback_wide, Timepoint == "post")
```
```{r}
# Function to rename columns with a prefix
add_prefix <- function(df, prefix) {
# Don't rename key columns (Subj, Stimulation_Site, Timepoint)
non_key_cols <- setdiff(names(df), c("Subj", "Stimulation_Site", "Timepoint"))
# Add the prefix to non-key columns
names(df)[names(df) %in% non_key_cols] <- paste(prefix, names(df)[names(df) %in% non_key_cols], sep = "_")
return(df)
}
# Add custom prefixes to each dataset
navon_wide_prefixed <- add_prefix(navon_wide, "navon")
stroop_wide_prefixed <- add_prefix(stroop_wide, "stroop")
nback_wide_prefixed <- add_prefix(nback_wide, "nback")
# Add custom prefixes to each dataset
navon_pre <- add_prefix(navon_pre, "navon")
stroop_pre <- add_prefix(stroop_pre, "stroop")
nback_pre <- add_prefix(nback_pre, "nback")
# Add custom prefixes to each dataset
navon_post <- add_prefix(navon_post, "navon")
stroop_post <- add_prefix(stroop_post, "stroop")
nback_post <- add_prefix(nback_post, "nback")
# Add custom prefixes to each dataset
navon_pre_c <- add_prefix(navon_pre, "navon_pre")
stroop_pre_c <- add_prefix(stroop_pre, "stroop_pre")
nback_pre_c <- add_prefix(nback_pre, "nback_pre")
# Add custom prefixes to each dataset
navon_post_c <- add_prefix(navon_post, "navon_post")
stroop_post_c <- add_prefix(stroop_post, "stroop_post")
nback_post_c <- add_prefix(nback_post, "nback_post")
```
## Handle Missing Data
```{r}
# Perform the join using merge
combined_data_pre <- Reduce(function(x, y) merge(x, y, by = c("Subj", "Stimulation_Site", "Timepoint"), all = TRUE),
list(navon_pre, stroop_pre, nback_pre)) # Add all your data tables
combined_data_post <- Reduce(function(x, y) merge(x, y, by = c("Subj", "Stimulation_Site", "Timepoint"), all = TRUE),
list(navon_post, stroop_post, nback_post)) # Add all your data tables
# Perform the join using merge
combined_data_pre_with_prefix <- Reduce(function(x, y) merge(x, y, by = c("Subj", "Stimulation_Site", "Timepoint"), all = TRUE),
list(navon_pre_c, stroop_pre_c, nback_pre_c)) # Add all your data tables
combined_data_post_with_prefix <- Reduce(function(x, y) merge(x, y, by = c("Subj", "Stimulation_Site", "Timepoint"), all = TRUE),
list(navon_post_c, stroop_post_c, nback_post_c)) # Add all your data tables
# Perform the join using merge
combined_data_wide <- Reduce(function(x, y) merge(x, y, by = c("Subj", "Stimulation_Site", "Timepoint"), all = TRUE),
list(navon_wide_prefixed, stroop_wide_prefixed, nback_wide_prefixed)) # Add all your data tables
# # Remove rows where Stimulation_Site equals 'DAN'
# combined_data_wide <- combined_data_wide[combined_data_wide$Stimulation_Site != "DAN", ]
```
## Combined both pre and post data into a single one
```{r}
# Remove the Timepoint column from both datasets
combined_data_pre_clean <- combined_data_pre_with_prefix %>% dplyr::select(-Timepoint)
combined_data_post_clean <- combined_data_post_with_prefix %>% dplyr::select(-Timepoint)
# Combine the datasets based on Subj and Stimulation_Site
combined_data_all <- full_join(combined_data_pre_clean, combined_data_post_clean, by = c("Subj", "Stimulation_Site"))
# View the combined dataset
head(combined_data_all)
```
# Factor Analysis for Response Time
```{r}
metric_type <- 'RT'
```
## Summary Statistics
```{r}
# Calculate Statistics for Response Times (Mean, Variance, Differences)
library(dplyr)
library(tidyr)
library(stringr)
# Filter for drift rate (v) columns
rt_cols_stats <- grep("_rt_", colnames(combined_data_all), value = TRUE)
# Create long format for statistics
stats_data <- combined_data_all %>%
dplyr::select(Subj, Stimulation_Site, all_of(rt_cols_stats)) %>%
pivot_longer(
cols = all_of(rt_cols_stats),
names_to = "full_name",
values_to = "Value"
) %>%
mutate(
Task = str_extract(full_name, "^[a-z]+"),
Timepoint = ifelse(grepl("_pre_", full_name), "Pre", "Post"),
Difficulty = str_extract(full_name, "(high|medium|low)$")
) %>%
filter(!is.na(Value))
# Calculate Summary Stats per Task
task_stats <- stats_data %>%
dplyr::select(Subj, Stimulation_Site, Task, Timepoint, Difficulty, Value) %>%
pivot_wider(names_from = Difficulty, values_from = Value) %>%
group_by(Task) %>%
summarise(
# Low Condition (Across all sites/timepoints)
Mean_Low = mean(low, na.rm = TRUE),
Var_Low = var(low, na.rm = TRUE),
# Medium Condition (Across all sites/timepoints)
Mean_Medium = mean(medium, na.rm = TRUE),
Var_Medium = var(medium, na.rm = TRUE),
# High Condition (Across all sites/timepoints)
Mean_High = mean(high, na.rm = TRUE),
Var_High = var(high, na.rm = TRUE),
# Individual Difference (High - Low)
Mean_Diff = mean(high - low, na.rm = TRUE),
Var_Diff = var(high - low, na.rm = TRUE)
)
print(task_stats)
```
```{r}
# 1. Reshape Data for Plotting
library(tidyr)
library(dplyr)
library(ggplot2)
library(stringr)
library(RColorBrewer)
# Filter for drift rate (v) columns
v_cols <- grep("_rt_", colnames(combined_data_all), value = TRUE)
# Pivot generic long format
plot_data_long <- combined_data_all %>%
dplyr::select(Subj, Stimulation_Site, all_of(v_cols)) %>%
pivot_longer(
cols = all_of(v_cols),
names_to = "full_name",
values_to = "Value"
) %>%
mutate(
# Extract Task: First word before underscore
Task = str_extract(full_name, "^[a-z]+"),
# Extract Timepoint: 'pre' or 'post'
Timepoint = ifelse(grepl("_pre_", full_name), "Pre", "Post"),
# Extract Difficulty: 'low', 'medium', 'high' at the end
Difficulty = str_extract(full_name, "(high|medium|low)$")
) %>%
mutate(
# Set factor levels for correct ordering
Difficulty = factor(str_to_title(Difficulty), levels = c("Low", "Medium", "High")),
Timepoint = factor(Timepoint, levels = c("Pre", "Post")),
Stimulation_Site = factor(Stimulation_Site, levels = c("Vertex", "FPCN-B", "DAN"), labels = c("Vertex", "L-FPN-B", "D-FPN"))
) %>%
filter(!is.na(Value))
# 2. Define Plotting Helper Function
draw_task_boxplot <- function(data, task_name) {
# Filter data for specific task
task_df <- data %>% filter(Task == task_name)
p <- ggplot(task_df, aes(x = Difficulty, y = Value)) +
# Facet by Stimulation Site to separate the groups clearly
facet_wrap(~ Stimulation_Site) +
# Boxplots:
# Fill by Stimulation Site (Consistent color within facet)
# Alpha by Timepoint (Distinguish Pre vs Post as 2 bars)
geom_boxplot(
aes(fill = Stimulation_Site, alpha = Timepoint),
position = position_dodge(width = 0.8),
width = 0.6,
outlier.shape = NA
) +
# Individual Points
geom_point(
aes(shape = Timepoint), # Shape matches Pre/Post
color = "gray25",
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.8),
size = 1.5,
alpha = 0.7
) +
# Aesthetics
scale_fill_brewer(palette = "Set1") +
scale_alpha_manual(values = c("Pre" = 0.4, "Post" = 0.9)) + # Light=Pre, Dark=Post
scale_shape_manual(values = c("Pre" = 16, "Post" = 17)) +
# Override only the alpha legend to show gray fills, keeping plot colors intact
guides(alpha = guide_legend(override.aes = list(fill = "black"))) +
theme_classic(base_size = 18) +
theme(
strip.background = element_rect(fill = "gray95", color = NA),
legend.position = "right",
plot.margin = margin(5.5, 5.5, 5.5, 15, "pt") # Increase left margin to prevent 'Response Time' cutoff
) +
labs(
title = NULL,
y = "Response Time (s)",
x = "Difficulty Condition",
fill = "Stimulation Site",
alpha = "Timepoint",
shape = "Timepoint"
)
# For N-Back, rotate labels to prevent overlap
if (task_name == "nback") {
p <- p + theme(axis.text.x = element_text(angle = 45, hjust = 1))
}
return(p)
}
# 3. Generate and Display Plots
p_navon <- draw_task_boxplot(plot_data_long, "navon")
print(p_navon)
ggsave("Figures/RT_Navon_v.png", plot = p_navon, width = 8, height = 5, dpi = 300)
p_stroop <- draw_task_boxplot(plot_data_long, "stroop")
print(p_stroop)
ggsave("Figures/RT_Stroop_v.png", plot = p_stroop, width = 8, height = 5, dpi = 300)
p_nback <- draw_task_boxplot(plot_data_long, "nback")
print(p_nback)
ggsave("Figures/RT_NBack_v.png", plot = p_nback, width = 8, height = 5, dpi = 300)
```
## Assumption Tests
We should check that our data fits the standard assumptions for a SEM (normality and linearity). Specifically, the relationships with the output should be linear, and the distributions of the input should be multivariate normal.
However, despite our inputs failing tests for normality, the MLR estimator builds robust standard errors that are able to handle non-normality.
```{r}
# Load necessary libraries
library(ggplot2)
library(MVN)
library(dplyr)
library(tidyr)
# 1) Draw histograms for navon, stroop, and nback columns in multiple subplots
# Reshape data for easy plotting
combined_data_long <- combined_data_wide %>%
dplyr::select(navon_rt_high, navon_rt_low, stroop_rt_high, stroop_rt_low, nback_rt_low, nback_rt_medium, nback_rt_high) %>%
tidyr::pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value")
# Plot histograms with each variable in a separate subplot
ggplot(combined_data_long, aes(x = Value)) +
geom_histogram(aes(fill = Variable), color = "black", bins = 20, alpha = 0.6) +
facet_wrap(~ Variable, scales = "free") +
labs(title = "Histograms for Navon, Stroop, and Nback Variables", x = "Value", y = "Frequency") +
theme_minimal() +
theme(legend.position = "none")
# 2) Multivariate normal distribution test
# Select only the relevant columns for testing
multivariate_data <- combined_data_wide %>%
dplyr::select(navon_rt_high, navon_rt_low, stroop_rt_high, stroop_rt_low, nback_rt_low, nback_rt_medium, nback_rt_high)
# Run Mardia's multivariate normality test
mardia_test <- mvn(multivariate_data, mvnTest = "mardia")
# Print results
print(mardia_test)
# Interpretation:
# If mardia_test$multivariateNormality$pValue.skew > 0.05 and mardia_test$multivariateNormality$pValue.kurt > 0.05,
# the data can be considered to follow a multivariate normal distribution.
```
### Linearity Test
We can check for linear relationships with each other. The output variable rt_general shows a general linear relationship with the otehr variables.
```{r}
# Load necessary packages
library(GGally)
library(dplyr)
# Select only the columns with 'rt_general' and the task-specific variables
plot_data <- combined_data_wide %>%
select(rt_general, navon_rt_high, navon_rt_low, stroop_rt_high, stroop_rt_low, nback_rt_low, nback_rt_medium, nback_rt_high)
# Create scatterplot matrix
p <- ggpairs(plot_data, columns = 1:ncol(plot_data),
title = "Scatterplots of rt_general with Task-specific Variables")
# Save to file with 300 DPI
ggsave("Figures/RT_Fig0_Linearity_Assumption.png", plot = p, dpi = 300, width = 12, height = 10)
```
## Run Factor Analysis (With Task Specific)
### Run Exploratory Factor Analysis with Parallel Analysis
This analysis is to determine whether the single factor or bifactor model makes more sense for our data.
For the parallel analysis, we can see when the unadjusted EV falls under the Random EV, at which point we do not want to use
that many factors. We can see in the output that 1 factor is sufficient, and 2 factors just misses the cutoff.
Therefore, we can justify our decision to use a single factor model.
#### Plot Parallel Analysis Scree Plot
```{r}
# Load necessary libraries
# install.packages(c("psych", "paran"))
library(psych)
library(paran)
library(dplyr)
# Subset the data to include only columns starting with "navon", "stroop", and "nback"
selected_data <- combined_data_wide[, grepl("^(navon|stroop|nback)_rt_", colnames(combined_data_wide))]
# selected_data$nback_rt_medium <- NULL
# Remove rows with any NA values
selected_data_complete <- na.omit(selected_data)
# Run parallel analysis on the complete dataset
paran_results <- paran(selected_data_complete, iterations = 1000, centile = 95, graph = TRUE)
# print(paran_results)
# Manually Graph Results
# Extract values from paran_results
adj_ev <- paran_results$AdjEv
ev <- paran_results$Ev
rnd_ev <- paran_results$RndEv
# Create data frame
plot_data <- data.frame(
Component = seq_along(ev),
Adjusted_EV = adj_ev,
Unadjusted_EV = ev,
Random_EV = rnd_ev
)
plot_data$Retained <- plot_data$Adjusted_EV > plot_data$Random_EV
# Create the plot
p <- ggplot(plot_data, aes(x = Component)) +
# Lines
geom_line(aes(y = Adjusted_EV, color = "Adjusted EV"), size = 1) +
geom_line(aes(y = Unadjusted_EV, color = "Unadjusted EV"), linetype = "dashed", size = 1) +
geom_line(aes(y = Random_EV, color = "Random EV"), linetype = "dotted", size = 1) +
# Points
geom_point(aes(y = Adjusted_EV, shape = Retained), size = 3, color = "black") +
geom_point(aes(y = Unadjusted_EV), size = 3, color = "red") +
geom_point(aes(y = Random_EV), size = 3, color = "blue") +
# Labels and theme
labs(
title = "Parallel Analysis Scree Plot",
x = "Number of Components",
y = "Eigenvalue",
color = "Legend",
shape = "Retained"
) +
scale_color_manual(values = c("Adjusted EV" = "black", "Unadjusted EV" = "red", "Random EV" = "blue")) +
scale_shape_manual(values = c(`TRUE` = 16, `FALSE` = 1)) +
theme_minimal(base_size = 14, base_family = "Arial") +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 28),
axis.title = element_text(size = 24),
axis.text = element_text(size = 16),
legend.title = element_text(size = 18), # Legend title font size
legend.text = element_text(size = 16), # Legend item labels
legend.box.background = element_rect(color = "black", fill = NA),
legend.box = "vertical"
)
# Save plot as high-resolution PNG
ggsave("Figures/RT_revised_parallel_analysis_plot.png", plot = p, width = 8, height = 6, dpi = 300, bg="white")
```
### Fit Model
Couple notes to keep track of:
Including the nback_rt_medium does cause the fit to get a lot worse. It's the variable with the lowest loading.
Additionally, the loading is not significant. However, removing this variable does not change the results, so I think
it's fine to proceed for now.
```{r}
library(lavaan)
set.seed(1234)
model <- '
# General factors (pre and post)
rt_general =~ navon_rt_low + navon_rt_high + stroop_rt_low + stroop_rt_high + nback_rt_low + nback_rt_high + nback_rt_medium
'
# Fit the model
fit <- cfa(model, data = combined_data_wide, estimator = "MLR", missing = "fiml", std.lv=TRUE)
# Summarize the results
summary(fit, fit.measures = TRUE, standardized = TRUE)
```
```{r}
library(lavaanPlot)
l = lavaanPlot(model = fit, coefs = TRUE)
print(l)
library(semPlot)
# Updated plotting code for top-to-bottom, black-and-white SEM diagram
custom_labels <- c("Navon V Low", "Navon V High", "Stroop V Low",
"Stroop V High", "N-Back V Low", "N-Back V High",
"N-Back V Medium", "V General")
# semPlot::semPaths(fit,
# what = "est", # Plot estimated coefficients
# edge.label.cex = 1.5, # Adjust size of edge labels
# node.label.cex = 1.5, # Adjust size of node labels
# sizeMan = 10,
# sizeLat = 10,
# sizeInt = 5,
# layout = "tree", # Layout style: tree for top-to-bottom
# rotation = 1, # Rotate diagram for top-to-bottom layout
# color = list(lat = "white", man = "white", int = "black"), # Black for edges, white for nodes
# edge.color = "black", # Black arrows
# label.color = "black",
# style = "lisrel", # Style for a clean, structured plot
# exoCov = FALSE, # Remove covariances for exogenous variables
# title = FALSE) # Apply custom variable labels) # Hide residuals
```
### Run Measurement Invariance Tests
Test configural invariance (model structure equivalence)
This step checks whether the basic factor structure (number of factors and factor loadings) is the same across groups.
The output seems to verify configural variance.
```{r}
fit_configural <- cfa(model, data = combined_data_wide, group = "Timepoint")
summary(fit_configural, fit.measures = TRUE)
```
Test metric invariance (equivalence of factor loadings)
Next, impose constraints on the factor loadings to be the same across groups.
The model does not pass the metric invariance test as clearly as it did the configural invariance test. The CFI and TLI are below ideal thresholds, and the SRMR is above the acceptable cutoff. The RMSEA is within the acceptable range, but the confidence interval suggests some variability in fit quality.
However, the chi-square test indicates that the fit of the model is not significantly worse, so you may choose to move forward with partial metric invariance testing by freeing some factor loadings if necessary or investigating specific sources of misfit.
By removing the nback_medium term, I get 'reasonable metric invariance' with slightly better fit terms. But for now, I will use the full model for completeness.
```{r}
fit_metric <- cfa(model, data = combined_data_wide, group = "Timepoint", group.equal = "loadings")
summary(fit_metric, fit.measures = TRUE)
```
### Compute Fit Metrics
We want to compute the same fit metrics as mentioned in Weigard.
We see that Omega Hierarchical for rt_general: 0.6309885, indicating that the task-general factor explains 63.1% of the data, which isn't perfect but moderate.
The Mean Lambda (λ) for rt_general: 0.456063, which is also a moderate loading, indicating that each model moderately loads onto the factor.
```{r}
library(psych)
library(lavaan)
# 1. Extract Fit Statistics
# The `summary` function with `fit.measures = TRUE` will provide key fit indices
summary(fit, fit.measures = TRUE, standardized = TRUE)
# 2. Extract standardized loadings to calculate mean λ
# Get the standardized solution
standardized_solution <- standardizedSolution(fit)
# Filter to only loadings for `rt_general` factor
rt_general_loadings <- standardized_solution[standardized_solution$lhs == "rt_general" & standardized_solution$op == "=~", "est.std"]
# Calculate mean λ for the `rt_general` factor
mean_lambda <- mean(rt_general_loadings)
cat("Mean Lambda (λ) for rt_general:", mean_lambda, "\n")
# 3. Omega Calculation for General Factor
# Use the `omega` function from the `psych` package
# First, extract the relevant columns from combined_data_wide for the general factor
rt_general_data <- combined_data_wide[, c("navon_rt_low", "navon_rt_high", "stroop_rt_low", "stroop_rt_high", "nback_rt_low", "nback_rt_high", "nback_rt_medium")]
# Run omega analysis
omega_result <- omega(rt_general_data, nfactors = 1) # Set `nfactors = 1` for one general factor
# Display omega hierarchical for rt_general
cat("Omega Hierarchical for rt_general:", omega_result$omega_h, "\n")
```
### Extract Subject Level Factors
```{r}
# Compute the factor scores for the general factor 'rt_general'
general_rt_rate <- lavPredict(fit, type = "lv") # 'lv' means latent variable
# Extract the 'rt_general' column
general_rt_rate <- as.data.frame(general_rt_rate)$rt_general
combined_data_wide$rt_general <- general_rt_rate
```
# Run Regressions
## Attach Network Metrics
```{r}
library(stringr)
library(dplyr)
# Read the connectivity file
connectivity_data <- read.csv("FinalData/bk_allnets_avg_output_final.csv")
# Extract subject ID
connectivity_data$Subj <- str_extract(connectivity_data$Row, "(?<=sub-)[A-Za-z]+(?=_big_corr)")
connectivity_data$Subj <- str_replace(connectivity_data$Subj, "^([A-Za-z])", "\\1_")
# Rename column to be "FileName"
names(connectivity_data)[names(connectivity_data) == "Row"] <- "FileName"
# Merge data
raw_metrics_full_data <- left_join(combined_data_wide, connectivity_data, by = "Subj")
```
## Attach Days Covariate
```{r}
# Load the necessary library
library(dplyr)
# Read the file
behavioral_data <- read.csv("FinalData/n-back_exp_results.csv")
# Extract unique combinations of Subj, Stimulation_Site, and Days
unique_combinations <- behavioral_data %>%
distinct(Subj, Stimulation_Site, Days, .keep_all = FALSE)
# View the results
print(unique_combinations)
library(data.table)
# Assuming unique_combinations is already a data.table. If not, convert it:
setDT(unique_combinations)
# Update Stimulation_Site values
unique_combinations[Stimulation_Site == "FPCNB", Stimulation_Site := "FPCN-B"]
unique_combinations[Stimulation_Site == "vertex", Stimulation_Site := "Vertex"]
# Merge data
raw_metrics_full_data_days <- left_join(raw_metrics_full_data, unique_combinations, by = c("Subj", "Stimulation_Site"))
```
## Look at baseline measures
```{r}
# Scale Days as well (MAKE SURE THESE METRICS WEREN'T SCALED BEFORE)
raw_metrics_full_data_days$Days = scale(raw_metrics_full_data_days$Days)
raw_metrics_full_data_days$FPCN_B.FPCN_B = scale(raw_metrics_full_data_days$FPCN_B.FPCN_B)
raw_metrics_full_data_days$DMN_Canonical.DMN_Canonical = scale(raw_metrics_full_data_days$DMN_Canonical.DMN_Canonical)
raw_metrics_full_data_days$DMN_Canonical.FPCN_B = scale(raw_metrics_full_data_days$DMN_Canonical.FPCN_B)
```
```{r}
library(lmerTest)
```
```{r}
# pre_data = filter(raw_metrics_full_data_days, Timepoint == "pre")
# model_0 <- lmer(rt_general ~ (DMN_Canonical.FPCN_B + DMN_Canonical.DMN_Canonical + FPCN_B.FPCN_B) + Days + (1 | Subj), data = pre_data, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
# summary(model_0)
pre_data = filter(raw_metrics_full_data_days, Timepoint == "pre")
model_0 <- lmer(rt_general ~ (DMN_Canonical.FPCN_B) + Days + (1 | Subj), data = pre_data, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
summary(model_0)
```
## Look at Raw Stimulation Effects without Network Moderators
```{r}
raw_metrics_full_data_days$Stimulation_Site <- factor(raw_metrics_full_data_days$Stimulation_Site)
# Set "Vertex" as the baseline (reference level)
raw_metrics_full_data_days$Stimulation_Site <- relevel(raw_metrics_full_data_days$Stimulation_Site, ref = "Vertex")
model_1a <- lmer(rt_general ~ Stimulation_Site * Timepoint + Days + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
summary(model_1a)
```
## Look at TMS Network Stimulation Effects
```{r}
# model_1 <- lmer(rt_general ~ (DMN_Canonical.FPCN_B + DMN_Canonical.DMN_Canonical + FPCN_B.FPCN_B) * Stimulation_Site * Timepoint + Days + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
# Ensure Stimulation_Site is a factor
raw_metrics_full_data_days$Stimulation_Site <- factor(raw_metrics_full_data_days$Stimulation_Site)
# Set "Vertex" as the baseline (reference level)
raw_metrics_full_data_days$Stimulation_Site <- relevel(raw_metrics_full_data_days$Stimulation_Site, ref = "Vertex")
model_1 <- lmer(rt_general ~ (DMN_Canonical.FPCN_B) * Stimulation_Site * Timepoint + Days + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
summary(model_1)
```
### Contrasts
```{r}
library(emmeans)
```
## Look at network moderation of Stimulation Effects
I'm commenting out the code to test out all three networks and only leaving in those for the FPCN-DMN anti-correlation
```{r}
# # Linear mixed-effects model to regress out 'Days'
# days_effect_rt_rate_model <- lm(rt_general ~ Days, data = raw_metrics_full_data_days)
# # Extract residuals which represent the part of 'v' not explained by 'Days'
# raw_metrics_full_data_days$rt_resid = residuals(days_effect_rt_rate_model)
#
# rt_rate_model_days_resid <- lmer(rt_resid ~ (DMN_Canonical.FPCN_B + DMN_Canonical.DMN_Canonical + FPCN_B.FPCN_B)*Stimulation_Site*Timepoint + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
# Linear mixed-effects model to regress out 'Days'
days_effect_rt_rate_model <- lm(rt_general ~ Days, data = raw_metrics_full_data_days)
# Extract residuals which represent the part of 'v' not explained by 'Days'
raw_metrics_full_data_days$rt_resid = residuals(days_effect_rt_rate_model)
rt_rate_model_days_resid <- lmer(rt_resid ~ (DMN_Canonical.FPCN_B)*Stimulation_Site*Timepoint + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
```
This is just the stimulation effects without any network moderators
```{r}
just_stim_model_resid <- lmer(rt_resid ~ Stimulation_Site*Timepoint + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
model_1_emmeans = emmeans(just_stim_model_resid, ~ Stimulation_Site * Timepoint, pbkrtest.limit = 30000, lmerTest.limit = 30000)
# Look at the pairwise comparisons for the interaction
pairs(model_1_emmeans)
# Contrast the change from pre to post for FPCN-B vs Vertex
mdl_1_emmeans = contrast(model_1_emmeans, interaction = c("revpairwise"), adjust = "none")
mdl_1_emmeans
```
Same analysis but using the full model
```{r}
just_stim_model_resid <- lmer(rt_resid ~ Stimulation_Site*Timepoint + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
model_1_just_stim_emmeans = emmeans(just_stim_model_resid, ~ Stimulation_Site, pbkrtest.limit = 30000, lmerTest.limit = 30000)
model_1_complete_emmeans = emmeans(rt_rate_model_days_resid, ~ Stimulation_Site, pbkrtest.limit = 30000, lmerTest.limit = 30000)
model_1_complete_emmeans_timepoint = emmeans(rt_rate_model_days_resid, ~ Timepoint, pbkrtest.limit = 30000, lmerTest.limit = 30000)
model_1_complete_emmeans_network = emtrends(rt_rate_model_days_resid, ~ 1, var = "DMN_Canonical.FPCN_B", pbkrtest.limit = 30000, lmerTest.limit = 30000)
# Look at the pairwise comparisons for the interaction
pairs(model_1_complete_emmeans)
# Contrast the change from pre to post for FPCN-B vs Vertex
model_1_complete_emmeans_int = emmeans(rt_rate_model_days_resid, ~ Stimulation_Site * Timepoint, pbkrtest.limit = 30000, lmerTest.limit = 30000)
mdl_1_emmeans = contrast(model_1_complete_emmeans, interaction = c("revpairwise"), adjust = "none")
mdl_1_emmeans
```
This is our main result below
```{r}
mdl_4_small = emtrends(rt_rate_model_days_resid, pairwise ~ Timepoint * Stimulation_Site, var="DMN_Canonical.FPCN_B", pbkrtest.limit = 30000, lmerTest.limit = 30000)
mdl_4_small_contrast = contrast(mdl_4_small[[1]], interaction = c("revpairwise"), adjust = "none")
mdl_4_small_contrast
```
# Draw Plots
```{r}
library(ggeffects)
library(dplyr)
library(ggplot2)
library(data.table)
# loadfonts(device = "win")
plot_metrics_stim_site = function(model, metric, metric_name, model_name){
# Remove hardcoded overrides to allow function to be generic
# metric_name = "FPCN-B and DMN\nConnectivity"
# model_name = "Response Time"
terms_vec = c(metric, "Stimulation_Site", "Timepoint")
preds <- predict_response(model, terms = terms_vec, interval="confidence", margin="mean_reference", back.transform = FALSE)
preds = as.data.table(preds)
contrast_preds = preds
preds <- preds %>%
mutate(
sem = (conf.high - conf.low) / (2 * 1.96) # approximate SEM from 95% CI
)
preds <- preds %>%
mutate(
ci.low_sem = predicted - sem,
ci.high_sem = predicted + sem
)
filtered_preds_df = preds
# --- Data Preparation for Individual Points (Added) ---
# Extract Random Intercepts
re_intercepts <- as.data.frame(ranef(model)$Subj)
re_intercepts$Subj <- rownames(re_intercepts)
intercept_col_idx <- grep("Intercept", names(re_intercepts))
re_intercepts <- re_intercepts[, c(ncol(re_intercepts), intercept_col_idx)]
names(re_intercepts) <- c("Subj", "RandomIntercept")
# Prepare raw data points with random intercept adjustment
# Note: model should be the residualized model 'rt_rate_model_days_resid'
# ensuring we use 'rt_resid' for consistency with the model response.
raw_data_plot <- raw_metrics_full_data_days %>%
mutate(across(all_of(metric), as.numeric)) %>%
left_join(re_intercepts, by = "Subj") %>%
mutate(
v_adjusted = rt_resid - RandomIntercept,
# Map Timepoint to match prediction labels if necessary
Timepoint = recode(Timepoint, "pre" = "pre", "post" = "post"),
facet = Timepoint # Use 'facet' column to match ggplot logical mapping
)
# Prepare Individual CONTRAST Data (Difference Scores)
subj_contrasts <- raw_data_plot %>%
dplyr::select(Subj, all_of(metric), Stimulation_Site, Timepoint, v_adjusted) %>%
pivot_wider(
id_cols = c(Subj, all_of(metric)),
names_from = c(Stimulation_Site, Timepoint),
values_from = v_adjusted,
names_sep = "_"
) %>%
mutate(
fpcnb_contrast = `FPCN-B_post` - `FPCN-B_pre`,
vertex_contrast = `Vertex_post` - `Vertex_pre`,
dan_contrast = `DAN_post` - `DAN_pre`,
fpcnb_vertex_contrast = fpcnb_contrast - vertex_contrast,
fpcnb_dan_contrast = fpcnb_contrast - dan_contrast,
dan_vertex_contrast = dan_contrast - vertex_contrast
)
contrast_df = dplyr::select(contrast_preds, -c(conf.low, conf.high)) %>%
pivot_wider(
names_from = c(group, facet),
values_from = c(predicted, std.error),
names_sep = "_"
)
final_contrasts <- contrast_df %>%
mutate(
fpcnb_contrast = `predicted_FPCN-B_post` - `predicted_FPCN-B_pre`,
vertex_contrast = predicted_Vertex_post - predicted_Vertex_pre,
dan_contrast = predicted_DAN_post - predicted_DAN_pre,
fpcnb_contrast.error = sqrt(`std.error_FPCN-B_post`^2 + `std.error_FPCN-B_pre`^2),
vertex_contrast.error = sqrt(std.error_Vertex_post^2 + std.error_Vertex_pre^2),
dan_contrast.error = sqrt(std.error_DAN_post^2 + std.error_DAN_pre^2),
fpcnb_vertex_contrast = fpcnb_contrast - vertex_contrast,
fpcnb_dan_contrast = fpcnb_contrast - dan_contrast,
dan_vertex_contrast = dan_contrast - vertex_contrast,
fpcnb_vertex_std.error = sqrt(fpcnb_contrast.error^2 + vertex_contrast.error^2),
fpcnb_dan_std.error = sqrt(fpcnb_contrast.error^2 + dan_contrast.error^2),
dan_vertex_std.error = sqrt(vertex_contrast.error^2 + dan_contrast.error^2)
) %>%
mutate(
dan_conf.low = dan_contrast - dan_contrast.error,
dan_conf.high = dan_contrast + dan_contrast.error,
fpcnb_conf.low = fpcnb_contrast - fpcnb_contrast.error,
fpcnb_conf.high = fpcnb_contrast + fpcnb_contrast.error,
vertex_conf.low = vertex_contrast - vertex_contrast.error,
vertex_conf.high = vertex_contrast + vertex_contrast.error,
fpcnb_vertex_conf.low = fpcnb_vertex_contrast - fpcnb_vertex_std.error,
fpcnb_vertex_conf.high = fpcnb_vertex_contrast + fpcnb_vertex_std.error,
fpcnb_dan_conf.low = fpcnb_dan_contrast - fpcnb_dan_std.error,
fpcnb_dan_conf.high = fpcnb_dan_contrast + fpcnb_dan_std.error,
dan_vertex_conf.low = dan_vertex_contrast - dan_vertex_std.error,
dan_vertex_conf.high = dan_vertex_contrast + dan_vertex_std.error
)
### Plotting
custom_theme <- theme_minimal(base_size = 24) +
theme(plot.title = element_text(size = rel(2.2), hjust = 0.5),
plot.background = element_blank(),
plot.margin = margin(t = 2, r = 1, b = 1, l = 30, unit = "pt"),
panel.background = element_rect(fill = "white"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_blank(),
axis.ticks = element_line(color = "black"),
axis.title.x = element_text(size = rel(3.0), margin = margin(t = 10), lineheight = 0.7),
axis.title.y = element_text(size = rel(3.0), margin = margin(r = 10), lineheight = 1.2),
axis.text = element_text(size = rel(2.7)),
strip.text = element_text(size = rel(2.2)),
legend.position = "right",
legend.title = element_blank(),
legend.background = element_rect(color = "black", size = .5),
legend.text = element_text(size = rel(2)),
legend.spacing.y = unit(0.5, "cm"),
panel.border = element_blank(),
text = element_text(family = "Arial"))
clean_metric_name <- gsub(" ", "", "FPCN-BandDMNConnectivity") # Kept consistent with old code behavior if needed, or use metric_name
# Ideally, use: clean_metric_name <- gsub(" ", "", metric_name)
# But assuming "FPCN-BandDMNConnectivity" was important for file naming consistency based on user prompt context "clean up",
# I will use the function argument `metric_name` logic.
clean_metric_name <- gsub("[\n ]", "", metric_name)
# --- Main Interaction Plots with Points ---
# Helper to plot site
plot_site <- function(site_name, file_suffix) {
site_preds <- filter(preds, group == site_name)
site_raw <- filter(raw_data_plot, Stimulation_Site == site_name)
p = ggplot(site_preds, aes(x = x, y = predicted, color = facet)) +
# Add Connecting Lines
geom_line(data = site_raw, aes(x = .data[[metric]], y = v_adjusted, group = Subj),
color = "gray80", alpha = 0.5) +
# Add Adjusted Individual Points
geom_point(data = site_raw, aes(x = .data[[metric]], y = v_adjusted, color = facet),
alpha = 1, size = 3) +
geom_ribbon(aes(ymin = ci.low_sem, ymax = ci.high_sem, fill = facet), linewidth = 0, alpha = 0.2) +
geom_line(size = 1.5) +
labs(y = model_name, x = metric_name) +
custom_theme
ggsave(filename = paste0("Figures/", metric_type, "_FilteredPreds_", file_suffix, "_", clean_metric_name, ".png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
}
plot_site("FPCN-B", "FPCNB")
plot_site("Vertex", "Vertex")
plot_site("DAN", "DAN")
# --- Contrast Plots ---
# Helper for contrasts
plot_contrast <- function(y_var, y_low, y_high, file_suffix, extra_theme = NULL) {
p = ggplot(final_contrasts, aes(x = x)) +
# Add Adjusted Individual Points (Contrasts)
geom_point(data = subj_contrasts, aes(x = .data[[metric]], y = .data[[y_var]]),
alpha = 0.4, size = 2.5, color = "black") +
geom_ribbon(aes(ymin = .data[[y_low]], ymax = .data[[y_high]]), linetype = "dotted", linewidth = 1, fill = "gray", alpha = 0.2, color = "gray") +
geom_line(size = 1.5, aes(y = .data[[y_var]]), color='black') +
labs(y = model_name, x = metric_name) +
custom_theme
if (!is.null(extra_theme)) {
p <- p + extra_theme
}
ggsave(filename = paste0("Figures/", metric_type, "_GeneralContrastOf_", file_suffix, "_", clean_metric_name, ".png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
}
plot_contrast("fpcnb_contrast", "fpcnb_conf.low", "fpcnb_conf.high", "FPCN")
plot_contrast("vertex_contrast", "vertex_conf.low", "vertex_conf.high", "Vertex")
plot_contrast("dan_contrast", "dan_conf.low", "dan_conf.high", "DAN")
plot_contrast("fpcnb_vertex_contrast", "fpcnb_vertex_conf.low", "fpcnb_vertex_conf.high", "FPCNVertex")
# FPCN-DAN contrast has specific margin
plot_contrast("fpcnb_dan_contrast", "fpcnb_dan_conf.low", "fpcnb_dan_conf.high", "FPCNDAN",
extra_theme = theme(plot.margin = margin(t = 20, r = 1, b = 1, l = 1, unit = "pt")))
plot_contrast("dan_vertex_contrast", "dan_vertex_conf.low", "dan_vertex_conf.high", "DANVertex")
}
```
```{r}
plot_metrics_stim_site(rt_rate_model_days_resid, "DMN_Canonical.FPCN_B", "FPCN-B and DMN\nConnectivity", "Response Time")
```
## Draw Network emmeans
```{r}
# Load necessary library
library(ggplot2)
preds <- predict_response(
model_1,
terms = c("DMN_Canonical.FPCN_B"),
interval = "confidence",
margin = "mean_reference",
back.transform = FALSE
)
preds = as.data.table(preds)
preds <- preds %>%
mutate(
sem = (conf.high - conf.low) / (2 * 1.96) # approximate SEM from 95% CI
)
preds <- preds %>%
mutate(
ci.low_sem = predicted - sem,
ci.high_sem = predicted + sem
)
preds_df = preds
# Calculate subject means for scatterplot and adjust for Random Intercepts
# This assumes the large spread is due to subject baseline differences (Random Intercepts).
# By subtracting the Random Intercept, we visualize the "Partial Residuals" - showing the
# relationship between Network and Response Time after controlling for individual baselines.
# 1. Extract Random Intercepts
re_intercepts <- as.data.frame(ranef(model_1)$Subj)
re_intercepts$Subj <- rownames(re_intercepts)
# Select only the intercept (in case there are random slopes) and rename
intercept_col_idx <- grep("Intercept", names(re_intercepts))
re_intercepts <- re_intercepts[, c(ncol(re_intercepts), intercept_col_idx)]
names(re_intercepts) <- c("Subj", "RandomIntercept")
# 2. Calculate adjusted means
subj_means <- raw_metrics_full_data_days %>%
mutate(DMN_Canonical.FPCN_B = as.numeric(DMN_Canonical.FPCN_B)) %>%
group_by(Subj) %>%
summarise(
mean_v = mean(rt_general, na.rm = TRUE),
network_val = mean(DMN_Canonical.FPCN_B, na.rm = TRUE)
) %>%
left_join(re_intercepts, by = "Subj") %>%
mutate(adjusted_mean_v = mean_v - RandomIntercept)
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"), # Center the title
plot.background = element_blank(), # White background for the plot
panel.background = element_rect(fill = "white"), # White background for the panels
panel.grid.major = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank(), # Remove minor grid lines
axis.line = element_blank(), # Define the axis lines without enclosing
axis.ticks = element_line(color = "black"), # Define ticks to make them clear
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2), # Add margin on top of x-axis title
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2), # Add margin on the right of y-axis title
axis.text = element_text(size = rel(1.8)), # Slightly larger axis text
strip.text = element_text(size = rel(2.2)), # Make facet titles bigger
legend.position = "right", # Position the legend on the right
legend.title = element_blank(), # Customize the legend title size
legend.background = element_rect(color = "black", size = .5), # Add a box around the legend
legend.text = element_text(size = rel(1.5)), # Customize the legend text size
legend.spacing.y = unit(0.5, "cm"), # Adjust the spacing between legend items
panel.border = element_blank(), # Ensure no border around the panels
text = element_text(family = "Arial")) # Apply Times New Roman (or serif font)
# Create the plot
p = ggplot(preds_df, aes(x = x, y = predicted)) +
geom_point(data = subj_means, aes(x = network_val, y = adjusted_mean_v),
color = "black", alpha = 0.6, size = 3) +
geom_line(size = 2, position = position_dodge(width = 0.5)) + # Add points
geom_ribbon(aes(ymin = ci.low_sem, ymax = ci.high_sem), linetype = "dotted", linewidth = 1, fill = "gray", alpha = 0.2) +
# scale_x_continuous(limits = c(-2, 2)) +
labs(title = "Network Connectivity",
y = "", x = "FPCN-B and DMN Connectivity") +
theme_minimal(base_size = 14) +
custom_theme +
theme(legend.position = "none")
# Use ggsave() to save the plot
ggsave(filename = paste0("Figures/RT_EMMeanNetworkDriftRate.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
```
## Draw Stimulation emmeans (No Timepoint)
```{r}
# Load necessary library
library(ggplot2)
library(emmeans)
# Create the data frame directly from the emmeans object
model_1_complete_emmeans_df <- as.data.frame(model_1_complete_emmeans)
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"), # Center the title
plot.background = element_blank(), # White background for the plot
panel.background = element_rect(fill = "white"), # White background for the panels
panel.grid.major = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank(), # Remove minor grid lines
axis.line = element_blank(), # Define the axis lines without enclosing
axis.ticks = element_line(color = "black"), # Define ticks to make them clear
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2), # Add margin on top of x-axis title
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2), # Add margin on the right of y-axis title
axis.text = element_text(size = rel(1.8)), # Slightly larger axis text
strip.text = element_text(size = rel(2.2)), # Make facet titles bigger
legend.position = "right", # Position the legend on the right
legend.title = element_blank(), # Customize the legend title size
legend.background = element_rect(color = "black", size = .5), # Add a box around the legend
legend.text = element_text(size = rel(1.5)), # Customize the legend text size
legend.spacing.y = unit(0.5, "cm"), # Adjust the spacing between legend items
panel.border = element_blank(), # Ensure no border around the panels
text = element_text(family = "Arial")) # Apply Times New Roman (or serif font)
# Convert Stimulation_Site to a factor and set the order
model_1_complete_emmeans_df$Stimulation_Site <- factor(model_1_complete_emmeans_df$Stimulation_Site, levels = c("FPCN-B", "Vertex", "DAN"))
# Create the plot
p = ggplot(model_1_complete_emmeans_df, aes(x = Stimulation_Site, y = emmean, color = Stimulation_Site)) +
geom_point(size = 6, position = position_dodge(width = 0.5)) + # Add points
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2, size = 2, position = position_dodge(width = 0.5)) + # Error bars
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Reference line
scale_color_manual(values = c("DAN" = "blue", "Vertex" = "gray", "FPCN-B" = "red")) + # Set colors for points
labs(title = "Stimulation Site",
y = "Estimated Marginal Mean of \nGeneral Response Time", x = "Stimulation Site") +
theme_minimal(base_size = 14) +
custom_theme +
theme(legend.position = "none")
# Use ggsave() to save the plot
ggsave(filename = paste0("Figures/RT_EMMeanStimSiteDriftRate_Automated.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
# Print the plot
print(p)
```
### Draw Stimulation Site separated by Timepoint
```{r}
# Load necessary libraries
library(ggplot2)
library(dplyr)
library(emmeans)
# Convert emmeans object to dataframe
# We use a new variable name to avoid overwriting the original object or the manual dataframe below
model_1_complete_emmeans_int_df <- as.data.frame(model_1_complete_emmeans_int)
# Use SE for error bars to match manual plot (Mean +/- SE)
model_1_complete_emmeans_int_df <- model_1_complete_emmeans_int_df %>%
mutate(
ci.low_sem = emmean - SE,
ci.high_sem = emmean + SE
)
# Capitalize Timepoint labels
model_1_complete_emmeans_int_df <- model_1_complete_emmeans_int_df %>%
mutate(Timepoint = recode(Timepoint, "pre" = "Pre", "post" = "Post"))
# Convert factors and set order for proper plotting
model_1_complete_emmeans_int_df$Stimulation_Site <- factor(
model_1_complete_emmeans_int_df$Stimulation_Site,
levels = c("FPCN-B", "Vertex", "DAN")
)
model_1_complete_emmeans_int_df$Timepoint <- factor(
model_1_complete_emmeans_int_df$Timepoint,
levels = c("Pre", "Post")
)
# Define custom theme
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"),
plot.background = element_blank(),
panel.background = element_rect(fill = "white"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(color = "black"),
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2),
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2),
axis.text = element_text(size = rel(1.8)),
strip.text = element_text(size = rel(2.2)),
legend.title = element_blank(),
legend.text = element_text(size = rel(1.5)),
legend.spacing.y = unit(0.5, "cm"),
panel.border = element_blank(),
text = element_text(family = "Arial"))
# Prepare raw data for plotting individual points with Random Intercept subtraction
# Extract Random Intercepts
re_intercepts <- as.data.frame(ranef(model_1)$Subj)
re_intercepts$Subj <- rownames(re_intercepts)
intercept_col_idx <- grep("Intercept", names(re_intercepts))
re_intercepts <- re_intercepts[, c(ncol(re_intercepts), intercept_col_idx)]
names(re_intercepts) <- c("Subj", "RandomIntercept")
raw_data_plot <- raw_metrics_full_data_days %>%
left_join(re_intercepts, by = "Subj") %>%
mutate(
rt_resid_centered = rt_resid - RandomIntercept,
Timepoint = recode(Timepoint, "pre" = "Pre", "post" = "Post"),
Stimulation_Site = factor(Stimulation_Site, levels = c("FPCN-B", "Vertex", "DAN")),
Timepoint = factor(Timepoint, levels = c("Pre", "Post"))
)
# Create the plot
p <- ggplot(model_1_complete_emmeans_int_df,
aes(x = Stimulation_Site, y = emmean, group = Timepoint)) +
# 1. Error bars (Bottom layer)
geom_errorbar(aes(ymin = ci.low_sem, ymax = ci.high_sem, color = Timepoint),
width = 0.2,
position = position_dodge(width = 0.5),
size = 1.2) +
# 2. Individual points (Middle layer, alpha increased for visibility)
geom_point(data = raw_data_plot, aes(y = rt_resid_centered, color = Timepoint),
position = position_jitterdodge(jitter.width = 0.25, dodge.width = 0.5),
alpha = 0.4, size = 2.5) +
# 3. Mean points (Top layer)
geom_point(size = 6, shape = 16, position = position_dodge(width = 0.5), aes(color = Timepoint)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Reference line
scale_color_manual(values = c("Pre" = "blue", "Post" = "red")) + # Different colors for error bars
labs(title = "Stimulation Site",
y = "General Response Time",
x = "Stimulation Site") +
theme_minimal(base_size = 14) +
custom_theme +
theme(legend.position = "right") # Keep the legend
# Save the plot
ggsave(filename = paste0("Figures/RT_EMMeanStimSiteAndTimepoint_Automated.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
print(p)
```
## Draw Timepoint EMMeans (automated)
```{r}
# Load necessary libraries
library(ggplot2)
library(dplyr)
library(emmeans)
# Create the data frame directly from the emmeans object
# Using a new variable name to avoid conflicts
model_1_complete_emmeans_timepoint_df <- as.data.frame(model_1_complete_emmeans_timepoint)
# Capitalize Timepoint labels to match manual plot
model_1_complete_emmeans_timepoint_df <- model_1_complete_emmeans_timepoint_df %>%
mutate(Timepoint = recode(Timepoint, "pre" = "Pre", "post" = "Post"))
# Convert Timepoint to a factor and set the order
model_1_complete_emmeans_timepoint_df$Timepoint <- factor(model_1_complete_emmeans_timepoint_df$Timepoint, levels = c("Pre", "Post"))
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"), # Center the title
plot.background = element_blank(), # White background for the plot
panel.background = element_rect(fill = "white"), # White background for the panels
panel.grid.major = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank(), # Remove minor grid lines
axis.line = element_blank(), # Define the axis lines without enclosing
axis.ticks = element_line(color = "black"), # Define ticks to make them clear
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2), # Add margin on top of x-axis title
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2), # Add margin on the right of y-axis title
axis.text = element_text(size = rel(1.8)), # Slightly larger axis text
strip.text = element_text(size = rel(2.2)), # Make facet titles bigger
legend.position = "right", # Position the legend on the right
legend.title = element_blank(), # Customize the legend title size
legend.background = element_rect(color = "black", size = .5), # Add a box around the legend
legend.text = element_text(size = rel(1.5)), # Customize the legend text size
legend.spacing.y = unit(0.5, "cm"), # Adjust the spacing between legend items
panel.border = element_blank(), # Ensure no border around the panels
text = element_text(family = "Arial")) # Apply Times New Roman (or serif font)
# Create the plot
p = ggplot(model_1_complete_emmeans_timepoint_df, aes(x = Timepoint, y = emmean, color = Timepoint)) +
geom_point(size = 6, position = position_dodge(width = 0.5)) + # Add points
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2, size = 2, position = position_dodge(width = 0.5)) + # Error bars
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Reference line
scale_color_manual(values = c("Pre" = "blue", "Post" = "red")) + # Set colors for points
labs(title = "Timepoint",
y = "Estimated Marginal Mean of \nGeneral Response Time", x = "Timepoint") +
theme_minimal(base_size = 14) +
custom_theme +
theme(legend.position = "none")
# Use ggsave() to save the plot
ggsave(filename = paste0("Figures/RT_EMMeanTimepointDriftRate_Automated.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
print(p)
```
## Draw Days Effects
```{r}
model = model_1
preds <- predict_response(model,
terms = c("Days"),
interval = "confidence",
margin = "mean_reference",
back.transform = FALSE)
preds = as.data.table(preds)
preds <- preds %>%
mutate(
sem = (conf.high - conf.low) / (2 * 1.96) # approximate SEM from 95% CI
)
preds <- preds %>%
mutate(
ci.low_sem = predicted - sem,
ci.high_sem = predicted + sem
)
filtered_preds_df = preds
# Create subject-centered data using Random Intercept subtraction
# Extract Random Intercepts
re_intercepts <- as.data.frame(ranef(model_1)$Subj)
re_intercepts$Subj <- rownames(re_intercepts)
intercept_col_idx <- grep("Intercept", names(re_intercepts))
re_intercepts <- re_intercepts[, c(ncol(re_intercepts), intercept_col_idx)]
names(re_intercepts) <- c("Subj", "RandomIntercept")
metrics_days_centered <- raw_metrics_full_data_days %>%
left_join(re_intercepts, by = "Subj") %>%
mutate(rt_general_centered = rt_general - RandomIntercept)
# Calculate intercepts and slopes for each subject based on raw data
# (Since the model only has random intercepts, we use individual OLS regressions
# to visualize the heterogeneity in slopes that exists in the raw data)
subj_trends <- metrics_days_centered %>%
group_by(Subj) %>%
summarise(
intercept = coef(lm(rt_general_centered ~ Days))[1],
slope = coef(lm(rt_general_centered ~ Days))[2]
) %>%
mutate(trend_color = ifelse(slope > 0, "green4", "red3"))
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"), # Center the title
plot.background = element_blank(), # White background for the plot
panel.background = element_rect(fill = "white"), # White background for the panels
panel.grid.major = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank(), # Remove minor grid lines
axis.line = element_blank(), # Define the axis lines without enclosing
axis.ticks = element_line(color = "black"), # Define ticks to make them clear
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2), # Add margin on top of x-axis title
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2), # Add margin on the right of y-axis title
axis.text = element_text(size = rel(1.8)), # Slightly larger axis text
strip.text = element_text(size = rel(2.2)), # Make facet titles bigger
legend.position = "right", # Position the legend on the right
legend.title = element_blank(), # Customize the legend title size
legend.background = element_rect(color = "black", size = .5), # Add a box around the legend
legend.text = element_text(size = rel(1.5)), # Customize the legend text size
legend.spacing.y = unit(0.5, "cm"), # Adjust the spacing between legend items
panel.border = element_blank(), # Ensure no border around the panels
text = element_text(family = "Arial")) # Apply Times New Roman (or serif font)
p = ggplot(filtered_preds_df, aes(x = x, y = predicted)) +
# Add individual participant trend lines (subject-centered)
geom_abline(data = subj_trends, aes(intercept = intercept, slope = slope, group = Subj, color = trend_color),
alpha = 0.2, size = 0.5) +
scale_color_identity() +
geom_ribbon(aes(ymin = ci.low_sem, ymax = ci.high_sem), linetype = "dotted", linewidth = 1, fill = "gray", alpha = 0.2) +
geom_line(size = 1.5) +
#Add confidence interval
labs(title = "Days",y = "", x = "Days") +
custom_theme
# Use ggsave() to save the plot
ggsave(filename = paste0("Figures/RT_DaysLine.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
print(p)
```
# Factor Analysis for Accuracy
```{r}
metric_type <- 'ACC'
```
## Summary Statistics
```{r}
# Calculate Statistics for Accuracys (Mean, Variance, Differences)
library(dplyr)
library(tidyr)
library(stringr)
# Filter for drift rate (v) columns
acc_cols_stats <- grep("_acc_", colnames(combined_data_all), value = TRUE)
# Create long format for statistics
stats_data <- combined_data_all %>%
dplyr::select(Subj, Stimulation_Site, all_of(acc_cols_stats)) %>%
pivot_longer(
cols = all_of(acc_cols_stats),
names_to = "full_name",
values_to = "Value"
) %>%
mutate(
Task = str_extract(full_name, "^[a-z]+"),
Timepoint = ifelse(grepl("_pre_", full_name), "Pre", "Post"),
Difficulty = str_extract(full_name, "(high|medium|low)$")
) %>%
filter(!is.na(Value))
# Calculate Summary Stats per Task
task_stats <- stats_data %>%
dplyr::select(Subj, Stimulation_Site, Task, Timepoint, Difficulty, Value) %>%
pivot_wider(names_from = Difficulty, values_from = Value) %>%
group_by(Task) %>%
summarise(
# Low Condition (Across all sites/timepoints)
Mean_Low = mean(low, na.rm = TRUE),
Var_Low = var(low, na.rm = TRUE),
# Medium Condition (Across all sites/timepoints)
Mean_Medium = mean(medium, na.rm = TRUE),
Var_Medium = var(medium, na.rm = TRUE),
# High Condition (Across all sites/timepoints)
Mean_High = mean(high, na.rm = TRUE),
Var_High = var(high, na.rm = TRUE),
# Individual Difference (High - Low)
Mean_Diff = mean(high - low, na.rm = TRUE),
Var_Diff = var(high - low, na.rm = TRUE)
)
print(task_stats)
```
```{r}
# 1. Reshape Data for Plotting
library(tidyr)
library(dplyr)
library(ggplot2)
library(stringr)
library(RColorBrewer)
# Filter for drift rate (v) columns
v_cols <- grep("_acc_", colnames(combined_data_all), value = TRUE)
# Pivot generic long format
plot_data_long <- combined_data_all %>%
dplyr::select(Subj, Stimulation_Site, all_of(v_cols)) %>%
pivot_longer(
cols = all_of(v_cols),
names_to = "full_name",
values_to = "Value"
) %>%
mutate(
# Extract Task: First word before underscore
Task = str_extract(full_name, "^[a-z]+"),
# Extract Timepoint: 'pre' or 'post'
Timepoint = ifelse(grepl("_pre_", full_name), "Pre", "Post"),
# Extract Difficulty: 'low', 'medium', 'high' at the end
Difficulty = str_extract(full_name, "(high|medium|low)$")
) %>%
mutate(
# Set factor levels for correct ordering
Difficulty = factor(str_to_title(Difficulty), levels = c("Low", "Medium", "High")),
Timepoint = factor(Timepoint, levels = c("Pre", "Post")),
Stimulation_Site = factor(Stimulation_Site, levels = c("Vertex", "FPCN-B", "DAN"), labels = c("Vertex", "L-FPN-B", "D-FPN"))
) %>%
filter(!is.na(Value))
# 2. Define Plotting Helper Function
draw_task_boxplot <- function(data, task_name) {
# Filter data for specific task
task_df <- data %>% filter(Task == task_name)
p <- ggplot(task_df, aes(x = Difficulty, y = Value)) +
# Facet by Stimulation Site to separate the groups clearly
facet_wrap(~ Stimulation_Site) +
# Boxplots:
# Fill by Stimulation Site (Consistent color within facet)
# Alpha by Timepoint (Distinguish Pre vs Post as 2 bars)
geom_boxplot(
aes(fill = Stimulation_Site, alpha = Timepoint),
position = position_dodge(width = 0.8),
width = 0.6,
outlier.shape = NA
) +
# Individual Points
geom_point(
aes(shape = Timepoint), # Shape matches Pre/Post
color = "gray25",
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.8),
size = 1.5,
alpha = 0.7
) +
# Aesthetics
scale_fill_brewer(palette = "Set1") +
scale_alpha_manual(values = c("Pre" = 0.4, "Post" = 0.9)) + # Light=Pre, Dark=Post
scale_shape_manual(values = c("Pre" = 16, "Post" = 17)) +
# Override only the alpha legend to show gray fills, keeping plot colors intact
guides(alpha = guide_legend(override.aes = list(fill = "black"))) +
theme_classic(base_size = 18) +
theme(
strip.background = element_rect(fill = "gray95", color = NA),
legend.position = "right",
plot.margin = margin(5.5, 5.5, 5.5, 15, "pt") # Increase left margin to prevent 'Accuracy' cutoff
) +
labs(
title = NULL,
y = "Accuracy (0-1)",
x = "Difficulty Condition",
fill = "Stimulation Site",
alpha = "Timepoint",
shape = "Timepoint"
)
# For N-Back, rotate labels to prevent overlap
if (task_name == "nback") {
p <- p + theme(axis.text.x = element_text(angle = 45, hjust = 1))
}
return(p)
}
# 3. Generate and Display Plots
p_navon <- draw_task_boxplot(plot_data_long, "navon")
print(p_navon)
ggsave("Figures/ACC_Navon_v.png", plot = p_navon, width = 8, height = 5, dpi = 300)
p_stroop <- draw_task_boxplot(plot_data_long, "stroop")
print(p_stroop)
ggsave("Figures/ACC_Stroop_v.png", plot = p_stroop, width = 8, height = 5, dpi = 300)
p_nback <- draw_task_boxplot(plot_data_long, "nback")
print(p_nback)
ggsave("Figures/ACC_NBack_v.png", plot = p_nback, width = 8, height = 5, dpi = 300)
```
## Assumption Tests
We should check that our data fits the standard assumptions for a SEM (normality and linearity). Specifically, the relationships with the output should be linear, and the distributions of the input should be multivariate normal.
However, despite our inputs failing tests for normality, the MLR estimator builds robust standard errors that are able to handle non-normality.
```{r}
# Load necessary libraries
library(ggplot2)
library(MVN)
library(dplyr)
library(tidyr)
# 1) Draw histograms for navon, stroop, and nback columns in multiple subplots
# Reshape data for easy plotting
combined_data_long <- combined_data_wide %>%
dplyr::select(navon_acc_high, navon_acc_low, stroop_acc_high, stroop_acc_low, nback_acc_low, nback_acc_medium, nback_acc_high) %>%
tidyr::pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value")
# Plot histograms with each variable in a separate subplot
ggplot(combined_data_long, aes(x = Value)) +
geom_histogram(aes(fill = Variable), color = "black", bins = 20, alpha = 0.6) +
facet_wrap(~ Variable, scales = "free") +
labs(title = "Histograms for Navon, Stroop, and Nback Variables", x = "Value", y = "Frequency") +
theme_minimal() +
theme(legend.position = "none")
# 2) Multivariate normal distribution test
# Select only the relevant columns for testing
multivariate_data <- combined_data_wide %>%
dplyr::select(navon_acc_high, navon_acc_low, stroop_acc_high, stroop_acc_low, nback_acc_low, nback_acc_medium, nback_acc_high)
# Run Mardia's multivariate normality test
mardia_test <- mvn(multivariate_data, mvnTest = "mardia")
# Print results
print(mardia_test)
# Interpretation:
# If mardia_test$multivariateNormality$pValue.skew > 0.05 and mardia_test$multivariateNormality$pValue.kurt > 0.05,
# the data can be considered to follow a multivariate normal distribution.
```
### Linearity Test
We can check for linear relationships with each other. The output variable acc_general shows a general linear relationship with the otehr variables.
```{r}
# Load necessary packages
library(GGally)
library(dplyr)
# Select only the columns with 'acc_general' and the task-specific variables
plot_data <- combined_data_wide %>%
select(acc_general, navon_acc_high, navon_acc_low, stroop_acc_high, stroop_acc_low, nback_acc_low, nback_acc_medium, nback_acc_high)
# Create scatterplot matrix
p <- ggpairs(plot_data, columns = 1:ncol(plot_data),
title = "Scatterplots of acc_general with Task-specific Variables")
# Save to file with 300 DPI
ggsave("Figures/ACC_Fig0_Linearity_Assumption.png", plot = p, dpi = 300, width = 12, height = 10)
```
## Run Factor Analysis (With Task Specific)
### Run Exploratory Factor Analysis with Parallel Analysis
This analysis is to determine whether the single factor or bifactor model makes more sense for our data.
For the parallel analysis, we can see when the unadjusted EV falls under the Random EV, at which point we do not want to use
that many factors. We can see in the output that 1 factor is sufficient, and 2 factors just misses the cutoff.
Therefore, we can justify our decision to use a single factor model.
#### Plot Parallel Analysis Scree Plot
```{r}
# Load necessary libraries
# install.packages(c("psych", "paran"))
library(psych)
library(paran)
library(dplyr)
# Subset the data to include only columns starting with "navon", "stroop", and "nback"
selected_data <- combined_data_wide[, grepl("^(navon|stroop|nback)_acc_", colnames(combined_data_wide))]
# selected_data$nback_acc_medium <- NULL
# Remove rows with any NA values
selected_data_complete <- na.omit(selected_data)
# Run parallel analysis on the complete dataset
paran_results <- paran(selected_data_complete, iterations = 1000, centile = 95, graph = TRUE)
# print(paran_results)
# Manually Graph Results
# Extract values from paran_results
adj_ev <- paran_results$AdjEv
ev <- paran_results$Ev
rnd_ev <- paran_results$RndEv
# Create data frame
plot_data <- data.frame(
Component = seq_along(ev),
Adjusted_EV = adj_ev,
Unadjusted_EV = ev,
Random_EV = rnd_ev
)
plot_data$Retained <- plot_data$Adjusted_EV > plot_data$Random_EV
# Create the plot
p <- ggplot(plot_data, aes(x = Component)) +
# Lines
geom_line(aes(y = Adjusted_EV, color = "Adjusted EV"), size = 1) +
geom_line(aes(y = Unadjusted_EV, color = "Unadjusted EV"), linetype = "dashed", size = 1) +
geom_line(aes(y = Random_EV, color = "Random EV"), linetype = "dotted", size = 1) +
# Points
geom_point(aes(y = Adjusted_EV, shape = Retained), size = 3, color = "black") +
geom_point(aes(y = Unadjusted_EV), size = 3, color = "red") +
geom_point(aes(y = Random_EV), size = 3, color = "blue") +
# Labels and theme
labs(
title = "Parallel Analysis Scree Plot",
x = "Number of Components",
y = "Eigenvalue",
color = "Legend",
shape = "Retained"
) +
scale_color_manual(values = c("Adjusted EV" = "black", "Unadjusted EV" = "red", "Random EV" = "blue")) +
scale_shape_manual(values = c(`TRUE` = 16, `FALSE` = 1)) +
theme_minimal(base_size = 14, base_family = "Arial") +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 28),
axis.title = element_text(size = 24),
axis.text = element_text(size = 16),
legend.title = element_text(size = 18), # Legend title font size
legend.text = element_text(size = 16), # Legend item labels
legend.box.background = element_rect(color = "black", fill = NA),
legend.box = "vertical"
)
# Save plot as high-resolution PNG
ggsave("Figures/ACC_revised_parallel_analysis_plot.png", plot = p, width = 8, height = 6, dpi = 300, bg="white")
```
### Fit Model
Couple notes to keep track of:
Including the nback_acc_medium does cause the fit to get a lot worse. It's the variable with the lowest loading.
Additionally, the loading is not significant. However, removing this variable does not change the results, so I think
it's fine to proceed for now.
```{r}
library(lavaan)
set.seed(1234)
model <- '
# General factors (pre and post)
acc_general =~ navon_acc_low + navon_acc_high + stroop_acc_low + stroop_acc_high + nback_acc_low + nback_acc_high + nback_acc_medium
'
# Fit the model
fit <- cfa(model, data = combined_data_wide, estimator = "MLR", missing = "fiml", std.lv=TRUE)
# Summarize the results
summary(fit, fit.measures = TRUE, standardized = TRUE)
```
```{r}
library(lavaanPlot)
l = lavaanPlot(model = fit, coefs = TRUE)
print(l)
library(semPlot)
# Updated plotting code for top-to-bottom, black-and-white SEM diagram
custom_labels <- c("Navon V Low", "Navon V High", "Stroop V Low",
"Stroop V High", "N-Back V Low", "N-Back V High",
"N-Back V Medium", "V General")
# semPlot::semPaths(fit,
# what = "est", # Plot estimated coefficients
# edge.label.cex = 1.5, # Adjust size of edge labels
# node.label.cex = 1.5, # Adjust size of node labels
# sizeMan = 10,
# sizeLat = 10,
# sizeInt = 5,
# layout = "tree", # Layout style: tree for top-to-bottom
# rotation = 1, # Rotate diagram for top-to-bottom layout
# color = list(lat = "white", man = "white", int = "black"), # Black for edges, white for nodes
# edge.color = "black", # Black arrows
# label.color = "black",
# style = "lisrel", # Style for a clean, structured plot
# exoCov = FALSE, # Remove covariances for exogenous variables
# title = FALSE) # Apply custom variable labels) # Hide residuals
```
### Run Measurement Invariance Tests
Test configural invariance (model structure equivalence)
This step checks whether the basic factor structure (number of factors and factor loadings) is the same across groups.
The output seems to verify configural variance.
```{r}
fit_configural <- cfa(model, data = combined_data_wide, group = "Timepoint")
summary(fit_configural, fit.measures = TRUE)
```
Test metric invariance (equivalence of factor loadings)
Next, impose constraints on the factor loadings to be the same across groups.
The model does not pass the metric invariance test as clearly as it did the configural invariance test. The CFI and TLI are below ideal thresholds, and the SRMR is above the acceptable cutoff. The RMSEA is within the acceptable range, but the confidence interval suggests some variability in fit quality.
However, the chi-square test indicates that the fit of the model is not significantly worse, so you may choose to move forward with partial metric invariance testing by freeing some factor loadings if necessary or investigating specific sources of misfit.
By removing the nback_medium term, I get 'reasonable metric invariance' with slightly better fit terms. But for now, I will use the full model for completeness.
```{r}
fit_metric <- cfa(model, data = combined_data_wide, group = "Timepoint", group.equal = "loadings")
summary(fit_metric, fit.measures = TRUE)
```
### Compute Fit Metrics
We want to compute the same fit metrics as mentioned in Weigard.
We see that Omega Hierarchical for acc_general: 0.6309885, indicating that the task-general factor explains 63.1% of the data, which isn't perfect but moderate.
The Mean Lambda (λ) for acc_general: 0.456063, which is also a moderate loading, indicating that each model moderately loads onto the factor.
```{r}
library(psych)
library(lavaan)
# 1. Extract Fit Statistics
# The `summary` function with `fit.measures = TRUE` will provide key fit indices
summary(fit, fit.measures = TRUE, standardized = TRUE)
# 2. Extract standardized loadings to calculate mean λ
# Get the standardized solution
standardized_solution <- standardizedSolution(fit)
# Filter to only loadings for `acc_general` factor
acc_general_loadings <- standardized_solution[standardized_solution$lhs == "acc_general" & standardized_solution$op == "=~", "est.std"]
# Calculate mean λ for the `acc_general` factor
mean_lambda <- mean(acc_general_loadings)
cat("Mean Lambda (λ) for acc_general:", mean_lambda, "\n")
# 3. Omega Calculation for General Factor
# Use the `omega` function from the `psych` package
# First, extract the relevant columns from combined_data_wide for the general factor
acc_general_data <- combined_data_wide[, c("navon_acc_low", "navon_acc_high", "stroop_acc_low", "stroop_acc_high", "nback_acc_low", "nback_acc_high", "nback_acc_medium")]
# Run omega analysis
omega_result <- omega(acc_general_data, nfactors = 1) # Set `nfactors = 1` for one general factor
# Display omega hierarchical for acc_general
cat("Omega Hierarchical for acc_general:", omega_result$omega_h, "\n")
```
### Extract Subject Level Factors
```{r}
# Compute the factor scores for the general factor 'acc_general'
general_acc_rate <- lavPredict(fit, type = "lv") # 'lv' means latent variable
# Extract the 'acc_general' column
general_acc_rate <- as.data.frame(general_acc_rate)$acc_general
combined_data_wide$acc_general <- general_acc_rate
```
# Run Regressions
## Attach Network Metrics
```{r}
library(stringr)
library(dplyr)
# Read the connectivity file
connectivity_data <- read.csv("FinalData/bk_allnets_avg_output_final.csv")
# Extract subject ID
connectivity_data$Subj <- str_extract(connectivity_data$Row, "(?<=sub-)[A-Za-z]+(?=_big_corr)")
connectivity_data$Subj <- str_replace(connectivity_data$Subj, "^([A-Za-z])", "\\1_")
# Rename column to be "FileName"
names(connectivity_data)[names(connectivity_data) == "Row"] <- "FileName"
# Merge data
raw_metrics_full_data <- left_join(combined_data_wide, connectivity_data, by = "Subj")
```
## Attach Days Covariate
```{r}
# Load the necessary library
library(dplyr)
# Read the file
behavioral_data <- read.csv("FinalData/n-back_exp_results.csv")
# Extract unique combinations of Subj, Stimulation_Site, and Days
unique_combinations <- behavioral_data %>%
distinct(Subj, Stimulation_Site, Days, .keep_all = FALSE)
# View the results
print(unique_combinations)
library(data.table)
# Assuming unique_combinations is already a data.table. If not, convert it:
setDT(unique_combinations)
# Update Stimulation_Site values
unique_combinations[Stimulation_Site == "FPCNB", Stimulation_Site := "FPCN-B"]
unique_combinations[Stimulation_Site == "vertex", Stimulation_Site := "Vertex"]
# Merge data
raw_metrics_full_data_days <- left_join(raw_metrics_full_data, unique_combinations, by = c("Subj", "Stimulation_Site"))
```
## Look at baseline measures
```{r}
# Scale Days as well (MAKE SURE THESE METRICS WEREN'T SCALED BEFORE)
# raw_metrics_full_data_days$Days = scale(raw_metrics_full_data_days$Days)
# raw_metrics_full_data_days$FPCN_B.FPCN_B = scale(raw_metrics_full_data_days$FPCN_B.FPCN_B)
# raw_metrics_full_data_days$DMN_Canonical.DMN_Canonical = scale(raw_metrics_full_data_days$DMN_Canonical.DMN_Canonical)
# raw_metrics_full_data_days$DMN_Canonical.FPCN_B = scale(raw_metrics_full_data_days$DMN_Canonical.FPCN_B)
```
```{r}
library(lmerTest)
```
```{r}
# pre_data = filter(raw_metrics_full_data_days, Timepoint == "pre")
# model_0 <- lmer(acc_general ~ (DMN_Canonical.FPCN_B + DMN_Canonical.DMN_Canonical + FPCN_B.FPCN_B) + Days + (1 | Subj), data = pre_data, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
# summary(model_0)
pre_data = filter(raw_metrics_full_data_days, Timepoint == "pre")
model_0 <- lmer(acc_general ~ (DMN_Canonical.FPCN_B) + Days + (1 | Subj), data = pre_data, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
summary(model_0)
```
## Look at Raw Stimulation Effects without Network Moderators
```{r}
raw_metrics_full_data_days$Stimulation_Site <- factor(raw_metrics_full_data_days$Stimulation_Site)
# Set "Vertex" as the baseline (reference level)
raw_metrics_full_data_days$Stimulation_Site <- relevel(raw_metrics_full_data_days$Stimulation_Site, ref = "Vertex")
model_1a <- lmer(acc_general ~ Stimulation_Site * Timepoint + Days + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
summary(model_1a)
```
## Look at TMS Network Stimulation Effects
```{r}
# model_1 <- lmer(acc_general ~ (DMN_Canonical.FPCN_B + DMN_Canonical.DMN_Canonical + FPCN_B.FPCN_B) * Stimulation_Site * Timepoint + Days + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
# Ensure Stimulation_Site is a factor
raw_metrics_full_data_days$Stimulation_Site <- factor(raw_metrics_full_data_days$Stimulation_Site)
# Set "Vertex" as the baseline (reference level)
raw_metrics_full_data_days$Stimulation_Site <- relevel(raw_metrics_full_data_days$Stimulation_Site, ref = "Vertex")
model_1 <- lmer(acc_general ~ (DMN_Canonical.FPCN_B) * Stimulation_Site * Timepoint + Days + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
summary(model_1)
```
### Contrasts
```{r}
library(emmeans)
```
## Look at network moderation of Stimulation Effects
I'm commenting out the code to test out all three networks and only leaving in those for the FPCN-DMN anti-correlation
```{r}
# # Linear mixed-effects model to regress out 'Days'
# days_effect_acc_rate_model <- lm(acc_general ~ Days, data = raw_metrics_full_data_days)
# # Extract residuals which represent the part of 'v' not explained by 'Days'
# raw_metrics_full_data_days$acc_resid = residuals(days_effect_acc_rate_model)
#
# acc_rate_model_days_resid <- lmer(acc_resid ~ (DMN_Canonical.FPCN_B + DMN_Canonical.DMN_Canonical + FPCN_B.FPCN_B)*Stimulation_Site*Timepoint + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
# Linear mixed-effects model to regress out 'Days'
days_effect_acc_rate_model <- lm(acc_general ~ Days, data = raw_metrics_full_data_days)
# Extract residuals which represent the part of 'v' not explained by 'Days'
raw_metrics_full_data_days$acc_resid = residuals(days_effect_acc_rate_model)
acc_rate_model_days_resid <- lmer(acc_resid ~ (DMN_Canonical.FPCN_B)*Stimulation_Site*Timepoint + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
```
This is just the stimulation effects without any network moderators
```{r}
just_stim_model_resid <- lmer(acc_resid ~ Stimulation_Site*Timepoint + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
model_1_emmeans = emmeans(just_stim_model_resid, ~ Stimulation_Site * Timepoint, pbkrtest.limit = 30000, lmerTest.limit = 30000)
# Look at the pairwise comparisons for the interaction
pairs(model_1_emmeans)
# Contrast the change from pre to post for FPCN-B vs Vertex
mdl_1_emmeans = contrast(model_1_emmeans, interaction = c("revpairwise"), adjust = "none")
mdl_1_emmeans
```
Same analysis but using the full model
```{r}
just_stim_model_resid <- lmer(acc_resid ~ Stimulation_Site*Timepoint + (1 | Subj), data = raw_metrics_full_data_days, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),REML=F)
model_1_just_stim_emmeans = emmeans(just_stim_model_resid, ~ Stimulation_Site, pbkrtest.limit = 30000, lmerTest.limit = 30000)
model_1_complete_emmeans = emmeans(acc_rate_model_days_resid, ~ Stimulation_Site, pbkrtest.limit = 30000, lmerTest.limit = 30000)
model_1_complete_emmeans_timepoint = emmeans(acc_rate_model_days_resid, ~ Timepoint, pbkrtest.limit = 30000, lmerTest.limit = 30000)
model_1_complete_emmeans_network = emtrends(acc_rate_model_days_resid, ~ 1, var = "DMN_Canonical.FPCN_B", pbkrtest.limit = 30000, lmerTest.limit = 30000)
# Look at the pairwise comparisons for the interaction
pairs(model_1_complete_emmeans)
# Contrast the change from pre to post for FPCN-B vs Vertex
model_1_complete_emmeans_int = emmeans(acc_rate_model_days_resid, ~ Stimulation_Site * Timepoint, pbkrtest.limit = 30000, lmerTest.limit = 30000)
mdl_1_emmeans = contrast(model_1_complete_emmeans, interaction = c("revpairwise"), adjust = "none")
mdl_1_emmeans
```
This is our main result below
```{r}
mdl_4_small = emtrends(acc_rate_model_days_resid, pairwise ~ Timepoint * Stimulation_Site, var="DMN_Canonical.FPCN_B", pbkrtest.limit = 30000, lmerTest.limit = 30000)
mdl_4_small_contrast = contrast(mdl_4_small[[1]], interaction = c("revpairwise"), adjust = "none")
mdl_4_small_contrast
```
# Draw Plots
```{r}
library(ggeffects)
library(dplyr)
library(ggplot2)
library(data.table)
# loadfonts(device = "win")
plot_metrics_stim_site = function(model, metric, metric_name, model_name){
# Remove hardcoded overrides to allow function to be generic
# metric_name = "FPCN-B and DMN\nConnectivity"
# model_name = "Accuracy"
terms_vec = c(metric, "Stimulation_Site", "Timepoint")
preds <- predict_response(model, terms = terms_vec, interval="confidence", margin="mean_reference", back.transform = FALSE)
preds = as.data.table(preds)
contrast_preds = preds
preds <- preds %>%
mutate(
sem = (conf.high - conf.low) / (2 * 1.96) # approximate SEM from 95% CI
)
preds <- preds %>%
mutate(
ci.low_sem = predicted - sem,
ci.high_sem = predicted + sem
)
filtered_preds_df = preds
# --- Data Preparation for Individual Points (Added) ---
# Extract Random Intercepts
re_intercepts <- as.data.frame(ranef(model)$Subj)
re_intercepts$Subj <- rownames(re_intercepts)
intercept_col_idx <- grep("Intercept", names(re_intercepts))
re_intercepts <- re_intercepts[, c(ncol(re_intercepts), intercept_col_idx)]
names(re_intercepts) <- c("Subj", "RandomIntercept")
# Prepare raw data points with random intercept adjustment
# Note: model should be the residualized model 'acc_rate_model_days_resid'
# ensuring we use 'acc_resid' for consistency with the model response.
raw_data_plot <- raw_metrics_full_data_days %>%
mutate(across(all_of(metric), as.numeric)) %>%
left_join(re_intercepts, by = "Subj") %>%
mutate(
v_adjusted = acc_resid - RandomIntercept,
# Map Timepoint to match prediction labels if necessary
Timepoint = recode(Timepoint, "pre" = "pre", "post" = "post"),
facet = Timepoint # Use 'facet' column to match ggplot logical mapping
)
# Prepare Individual CONTRAST Data (Difference Scores)
subj_contrasts <- raw_data_plot %>%
dplyr::select(Subj, all_of(metric), Stimulation_Site, Timepoint, v_adjusted) %>%
pivot_wider(
id_cols = c(Subj, all_of(metric)),
names_from = c(Stimulation_Site, Timepoint),
values_from = v_adjusted,
names_sep = "_"
) %>%
mutate(
fpcnb_contrast = `FPCN-B_post` - `FPCN-B_pre`,
vertex_contrast = `Vertex_post` - `Vertex_pre`,
dan_contrast = `DAN_post` - `DAN_pre`,
fpcnb_vertex_contrast = fpcnb_contrast - vertex_contrast,
fpcnb_dan_contrast = fpcnb_contrast - dan_contrast,
dan_vertex_contrast = dan_contrast - vertex_contrast
)
contrast_df = dplyr::select(contrast_preds, -c(conf.low, conf.high)) %>%
pivot_wider(
names_from = c(group, facet),
values_from = c(predicted, std.error),
names_sep = "_"
)
final_contrasts <- contrast_df %>%
mutate(
fpcnb_contrast = `predicted_FPCN-B_post` - `predicted_FPCN-B_pre`,
vertex_contrast = predicted_Vertex_post - predicted_Vertex_pre,
dan_contrast = predicted_DAN_post - predicted_DAN_pre,
fpcnb_contrast.error = sqrt(`std.error_FPCN-B_post`^2 + `std.error_FPCN-B_pre`^2),
vertex_contrast.error = sqrt(std.error_Vertex_post^2 + std.error_Vertex_pre^2),
dan_contrast.error = sqrt(std.error_DAN_post^2 + std.error_DAN_pre^2),
fpcnb_vertex_contrast = fpcnb_contrast - vertex_contrast,
fpcnb_dan_contrast = fpcnb_contrast - dan_contrast,
dan_vertex_contrast = dan_contrast - vertex_contrast,
fpcnb_vertex_std.error = sqrt(fpcnb_contrast.error^2 + vertex_contrast.error^2),
fpcnb_dan_std.error = sqrt(fpcnb_contrast.error^2 + dan_contrast.error^2),
dan_vertex_std.error = sqrt(vertex_contrast.error^2 + dan_contrast.error^2)
) %>%
mutate(
dan_conf.low = dan_contrast - dan_contrast.error,
dan_conf.high = dan_contrast + dan_contrast.error,
fpcnb_conf.low = fpcnb_contrast - fpcnb_contrast.error,
fpcnb_conf.high = fpcnb_contrast + fpcnb_contrast.error,
vertex_conf.low = vertex_contrast - vertex_contrast.error,
vertex_conf.high = vertex_contrast + vertex_contrast.error,
fpcnb_vertex_conf.low = fpcnb_vertex_contrast - fpcnb_vertex_std.error,
fpcnb_vertex_conf.high = fpcnb_vertex_contrast + fpcnb_vertex_std.error,
fpcnb_dan_conf.low = fpcnb_dan_contrast - fpcnb_dan_std.error,
fpcnb_dan_conf.high = fpcnb_dan_contrast + fpcnb_dan_std.error,
dan_vertex_conf.low = dan_vertex_contrast - dan_vertex_std.error,
dan_vertex_conf.high = dan_vertex_contrast + dan_vertex_std.error
)
### Plotting
custom_theme <- theme_minimal(base_size = 24) +
theme(plot.title = element_text(size = rel(2.2), hjust = 0.5),
plot.background = element_blank(),
plot.margin = margin(t = 2, r = 1, b = 1, l = 30, unit = "pt"),
panel.background = element_rect(fill = "white"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_blank(),
axis.ticks = element_line(color = "black"),
axis.title.x = element_text(size = rel(3.0), margin = margin(t = 10), lineheight = 0.7),
axis.title.y = element_text(size = rel(3.0), margin = margin(r = 10), lineheight = 1.2),
axis.text = element_text(size = rel(2.7)),
strip.text = element_text(size = rel(2.2)),
legend.position = "right",
legend.title = element_blank(),
legend.background = element_rect(color = "black", size = .5),
legend.text = element_text(size = rel(2)),
legend.spacing.y = unit(0.5, "cm"),
panel.border = element_blank(),
text = element_text(family = "Arial"))
clean_metric_name <- gsub(" ", "", "FPCN-BandDMNConnectivity") # Kept consistent with old code behavior if needed, or use metric_name
# Ideally, use: clean_metric_name <- gsub(" ", "", metric_name)
# But assuming "FPCN-BandDMNConnectivity" was important for file naming consistency based on user prompt context "clean up",
# I will use the function argument `metric_name` logic.
clean_metric_name <- gsub("[\n ]", "", metric_name)
# --- Main Interaction Plots with Points ---
# Helper to plot site
plot_site <- function(site_name, file_suffix) {
site_preds <- filter(preds, group == site_name)
site_raw <- filter(raw_data_plot, Stimulation_Site == site_name)
p = ggplot(site_preds, aes(x = x, y = predicted, color = facet)) +
# Add Connecting Lines
geom_line(data = site_raw, aes(x = .data[[metric]], y = v_adjusted, group = Subj),
color = "gray80", alpha = 0.5) +
# Add Adjusted Individual Points
geom_point(data = site_raw, aes(x = .data[[metric]], y = v_adjusted, color = facet),
alpha = 1, size = 3) +
geom_ribbon(aes(ymin = ci.low_sem, ymax = ci.high_sem, fill = facet), linewidth = 0, alpha = 0.2) +
geom_line(size = 1.5) +
labs(y = model_name, x = metric_name) +
custom_theme
ggsave(filename = paste0("Figures/", metric_type, "_FilteredPreds_", file_suffix, "_", clean_metric_name, ".png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
}
plot_site("FPCN-B", "FPCNB")
plot_site("Vertex", "Vertex")
plot_site("DAN", "DAN")
# --- Contrast Plots ---
# Helper for contrasts
plot_contrast <- function(y_var, y_low, y_high, file_suffix, extra_theme = NULL) {
p = ggplot(final_contrasts, aes(x = x)) +
# Add Adjusted Individual Points (Contrasts)
geom_point(data = subj_contrasts, aes(x = .data[[metric]], y = .data[[y_var]]),
alpha = 0.4, size = 2.5, color = "black") +
geom_ribbon(aes(ymin = .data[[y_low]], ymax = .data[[y_high]]), linetype = "dotted", linewidth = 1, fill = "gray", alpha = 0.2, color = "gray") +
geom_line(size = 1.5, aes(y = .data[[y_var]]), color='black') +
labs(y = model_name, x = metric_name) +
custom_theme
if (!is.null(extra_theme)) {
p <- p + extra_theme
}
ggsave(filename = paste0("Figures/", metric_type, "_GeneralContrastOf_", file_suffix, "_", clean_metric_name, ".png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
}
plot_contrast("fpcnb_contrast", "fpcnb_conf.low", "fpcnb_conf.high", "FPCN")
plot_contrast("vertex_contrast", "vertex_conf.low", "vertex_conf.high", "Vertex")
plot_contrast("dan_contrast", "dan_conf.low", "dan_conf.high", "DAN")
plot_contrast("fpcnb_vertex_contrast", "fpcnb_vertex_conf.low", "fpcnb_vertex_conf.high", "FPCNVertex")
# FPCN-DAN contrast has specific margin
plot_contrast("fpcnb_dan_contrast", "fpcnb_dan_conf.low", "fpcnb_dan_conf.high", "FPCNDAN",
extra_theme = theme(plot.margin = margin(t = 20, r = 1, b = 1, l = 1, unit = "pt")))
plot_contrast("dan_vertex_contrast", "dan_vertex_conf.low", "dan_vertex_conf.high", "DANVertex")
}
```
```{r}
plot_metrics_stim_site(acc_rate_model_days_resid, "DMN_Canonical.FPCN_B", "FPCN-B and DMN\nConnectivity", "Accuracy")
```
## Draw Network emmeans
```{r}
# Load necessary library
library(ggplot2)
preds <- predict_response(
model_1,
terms = c("DMN_Canonical.FPCN_B"),
interval = "confidence",
margin = "mean_reference",
back.transform = FALSE
)
preds = as.data.table(preds)
preds <- preds %>%
mutate(
sem = (conf.high - conf.low) / (2 * 1.96) # approximate SEM from 95% CI
)
preds <- preds %>%
mutate(
ci.low_sem = predicted - sem,
ci.high_sem = predicted + sem
)
preds_df = preds
# Calculate subject means for scatterplot and adjust for Random Intercepts
# This assumes the large spread is due to subject baseline differences (Random Intercepts).
# By subtracting the Random Intercept, we visualize the "Partial Residuals" - showing the
# relationship between Network and Accuracy after controlling for individual baselines.
# 1. Extract Random Intercepts
re_intercepts <- as.data.frame(ranef(model_1)$Subj)
re_intercepts$Subj <- rownames(re_intercepts)
# Select only the intercept (in case there are random slopes) and rename
intercept_col_idx <- grep("Intercept", names(re_intercepts))
re_intercepts <- re_intercepts[, c(ncol(re_intercepts), intercept_col_idx)]
names(re_intercepts) <- c("Subj", "RandomIntercept")
# 2. Calculate adjusted means
subj_means <- raw_metrics_full_data_days %>%
mutate(DMN_Canonical.FPCN_B = as.numeric(DMN_Canonical.FPCN_B)) %>%
group_by(Subj) %>%
summarise(
mean_v = mean(acc_general, na.rm = TRUE),
network_val = mean(DMN_Canonical.FPCN_B, na.rm = TRUE)
) %>%
left_join(re_intercepts, by = "Subj") %>%
mutate(adjusted_mean_v = mean_v - RandomIntercept)
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"), # Center the title
plot.background = element_blank(), # White background for the plot
panel.background = element_rect(fill = "white"), # White background for the panels
panel.grid.major = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank(), # Remove minor grid lines
axis.line = element_blank(), # Define the axis lines without enclosing
axis.ticks = element_line(color = "black"), # Define ticks to make them clear
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2), # Add margin on top of x-axis title
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2), # Add margin on the right of y-axis title
axis.text = element_text(size = rel(1.8)), # Slightly larger axis text
strip.text = element_text(size = rel(2.2)), # Make facet titles bigger
legend.position = "right", # Position the legend on the right
legend.title = element_blank(), # Customize the legend title size
legend.background = element_rect(color = "black", size = .5), # Add a box around the legend
legend.text = element_text(size = rel(1.5)), # Customize the legend text size
legend.spacing.y = unit(0.5, "cm"), # Adjust the spacing between legend items
panel.border = element_blank(), # Ensure no border around the panels
text = element_text(family = "Arial")) # Apply Times New Roman (or serif font)
# Create the plot
p = ggplot(preds_df, aes(x = x, y = predicted)) +
geom_point(data = subj_means, aes(x = network_val, y = adjusted_mean_v),
color = "black", alpha = 0.6, size = 3) +
geom_line(size = 2, position = position_dodge(width = 0.5)) + # Add points
geom_ribbon(aes(ymin = ci.low_sem, ymax = ci.high_sem), linetype = "dotted", linewidth = 1, fill = "gray", alpha = 0.2) +
# scale_x_continuous(limits = c(-2, 2)) +
labs(title = "Network Connectivity",
y = "", x = "FPCN-B and DMN Connectivity") +
theme_minimal(base_size = 14) +
custom_theme +
theme(legend.position = "none")
# Use ggsave() to save the plot
ggsave(filename = paste0("Figures/ACC_EMMeanNetworkDriftRate.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
```
## Draw Stimulation emmeans (No Timepoint)
```{r}
# Load necessary library
library(ggplot2)
library(emmeans)
# Create the data frame directly from the emmeans object
model_1_complete_emmeans_df <- as.data.frame(model_1_complete_emmeans)
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"), # Center the title
plot.background = element_blank(), # White background for the plot
panel.background = element_rect(fill = "white"), # White background for the panels
panel.grid.major = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank(), # Remove minor grid lines
axis.line = element_blank(), # Define the axis lines without enclosing
axis.ticks = element_line(color = "black"), # Define ticks to make them clear
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2), # Add margin on top of x-axis title
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2), # Add margin on the right of y-axis title
axis.text = element_text(size = rel(1.8)), # Slightly larger axis text
strip.text = element_text(size = rel(2.2)), # Make facet titles bigger
legend.position = "right", # Position the legend on the right
legend.title = element_blank(), # Customize the legend title size
legend.background = element_rect(color = "black", size = .5), # Add a box around the legend
legend.text = element_text(size = rel(1.5)), # Customize the legend text size
legend.spacing.y = unit(0.5, "cm"), # Adjust the spacing between legend items
panel.border = element_blank(), # Ensure no border around the panels
text = element_text(family = "Arial")) # Apply Times New Roman (or serif font)
# Convert Stimulation_Site to a factor and set the order
model_1_complete_emmeans_df$Stimulation_Site <- factor(model_1_complete_emmeans_df$Stimulation_Site, levels = c("FPCN-B", "Vertex", "DAN"))
# Create the plot
p = ggplot(model_1_complete_emmeans_df, aes(x = Stimulation_Site, y = emmean, color = Stimulation_Site)) +
geom_point(size = 6, position = position_dodge(width = 0.5)) + # Add points
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2, size = 2, position = position_dodge(width = 0.5)) + # Error bars
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Reference line
scale_color_manual(values = c("DAN" = "blue", "Vertex" = "gray", "FPCN-B" = "red")) + # Set colors for points
labs(title = "Stimulation Site",
y = "Estimated Marginal Mean of \nGeneral Accuracy", x = "Stimulation Site") +
theme_minimal(base_size = 14) +
custom_theme +
theme(legend.position = "none")
# Use ggsave() to save the plot
ggsave(filename = paste0("Figures/ACC_EMMeanStimSiteDriftRate_Automated.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
# Print the plot
print(p)
```
### Draw Stimulation Site separated by Timepoint
```{r}
# Load necessary libraries
library(ggplot2)
library(dplyr)
library(emmeans)
# Convert emmeans object to dataframe
# We use a new variable name to avoid overwriting the original object or the manual dataframe below
model_1_complete_emmeans_int_df <- as.data.frame(model_1_complete_emmeans_int)
# Use SE for error bars to match manual plot (Mean +/- SE)
model_1_complete_emmeans_int_df <- model_1_complete_emmeans_int_df %>%
mutate(
ci.low_sem = emmean - SE,
ci.high_sem = emmean + SE
)
# Capitalize Timepoint labels
model_1_complete_emmeans_int_df <- model_1_complete_emmeans_int_df %>%
mutate(Timepoint = recode(Timepoint, "pre" = "Pre", "post" = "Post"))
# Convert factors and set order for proper plotting
model_1_complete_emmeans_int_df$Stimulation_Site <- factor(
model_1_complete_emmeans_int_df$Stimulation_Site,
levels = c("FPCN-B", "Vertex", "DAN")
)
model_1_complete_emmeans_int_df$Timepoint <- factor(
model_1_complete_emmeans_int_df$Timepoint,
levels = c("Pre", "Post")
)
# Define custom theme
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"),
plot.background = element_blank(),
panel.background = element_rect(fill = "white"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(color = "black"),
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2),
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2),
axis.text = element_text(size = rel(1.8)),
strip.text = element_text(size = rel(2.2)),
legend.title = element_blank(),
legend.text = element_text(size = rel(1.5)),
legend.spacing.y = unit(0.5, "cm"),
panel.border = element_blank(),
text = element_text(family = "Arial"))
# Prepare raw data for plotting individual points with Random Intercept subtraction
# Extract Random Intercepts
re_intercepts <- as.data.frame(ranef(model_1)$Subj)
re_intercepts$Subj <- rownames(re_intercepts)
intercept_col_idx <- grep("Intercept", names(re_intercepts))
re_intercepts <- re_intercepts[, c(ncol(re_intercepts), intercept_col_idx)]
names(re_intercepts) <- c("Subj", "RandomIntercept")
raw_data_plot <- raw_metrics_full_data_days %>%
left_join(re_intercepts, by = "Subj") %>%
mutate(
acc_resid_centered = acc_resid - RandomIntercept,
Timepoint = recode(Timepoint, "pre" = "Pre", "post" = "Post"),
Stimulation_Site = factor(Stimulation_Site, levels = c("FPCN-B", "Vertex", "DAN")),
Timepoint = factor(Timepoint, levels = c("Pre", "Post"))
)
# Create the plot
p <- ggplot(model_1_complete_emmeans_int_df,
aes(x = Stimulation_Site, y = emmean, group = Timepoint)) +
# 1. Error bars (Bottom layer)
geom_errorbar(aes(ymin = ci.low_sem, ymax = ci.high_sem, color = Timepoint),
width = 0.2,
position = position_dodge(width = 0.5),
size = 1.2) +
# 2. Individual points (Middle layer, alpha increased for visibility)
geom_point(data = raw_data_plot, aes(y = acc_resid_centered, color = Timepoint),
position = position_jitterdodge(jitter.width = 0.25, dodge.width = 0.5),
alpha = 0.4, size = 2.5) +
# 3. Mean points (Top layer)
geom_point(size = 6, shape = 16, position = position_dodge(width = 0.5), aes(color = Timepoint)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Reference line
scale_color_manual(values = c("Pre" = "blue", "Post" = "red")) + # Different colors for error bars
labs(title = "Stimulation Site",
y = "General Accuracy",
x = "Stimulation Site") +
theme_minimal(base_size = 14) +
custom_theme +
theme(legend.position = "right") # Keep the legend
# Save the plot
ggsave(filename = paste0("Figures/ACC_EMMeanStimSiteAndTimepoint_Automated.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
print(p)
```
## Draw Timepoint EMMeans (automated)
```{r}
# Load necessary libraries
library(ggplot2)
library(dplyr)
library(emmeans)
# Create the data frame directly from the emmeans object
# Using a new variable name to avoid conflicts
model_1_complete_emmeans_timepoint_df <- as.data.frame(model_1_complete_emmeans_timepoint)
# Capitalize Timepoint labels to match manual plot
model_1_complete_emmeans_timepoint_df <- model_1_complete_emmeans_timepoint_df %>%
mutate(Timepoint = recode(Timepoint, "pre" = "Pre", "post" = "Post"))
# Convert Timepoint to a factor and set the order
model_1_complete_emmeans_timepoint_df$Timepoint <- factor(model_1_complete_emmeans_timepoint_df$Timepoint, levels = c("Pre", "Post"))
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"), # Center the title
plot.background = element_blank(), # White background for the plot
panel.background = element_rect(fill = "white"), # White background for the panels
panel.grid.major = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank(), # Remove minor grid lines
axis.line = element_blank(), # Define the axis lines without enclosing
axis.ticks = element_line(color = "black"), # Define ticks to make them clear
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2), # Add margin on top of x-axis title
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2), # Add margin on the right of y-axis title
axis.text = element_text(size = rel(1.8)), # Slightly larger axis text
strip.text = element_text(size = rel(2.2)), # Make facet titles bigger
legend.position = "right", # Position the legend on the right
legend.title = element_blank(), # Customize the legend title size
legend.background = element_rect(color = "black", size = .5), # Add a box around the legend
legend.text = element_text(size = rel(1.5)), # Customize the legend text size
legend.spacing.y = unit(0.5, "cm"), # Adjust the spacing between legend items
panel.border = element_blank(), # Ensure no border around the panels
text = element_text(family = "Arial")) # Apply Times New Roman (or serif font)
# Create the plot
p = ggplot(model_1_complete_emmeans_timepoint_df, aes(x = Timepoint, y = emmean, color = Timepoint)) +
geom_point(size = 6, position = position_dodge(width = 0.5)) + # Add points
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2, size = 2, position = position_dodge(width = 0.5)) + # Error bars
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Reference line
scale_color_manual(values = c("Pre" = "blue", "Post" = "red")) + # Set colors for points
labs(title = "Timepoint",
y = "Estimated Marginal Mean of \nGeneral Accuracy", x = "Timepoint") +
theme_minimal(base_size = 14) +
custom_theme +
theme(legend.position = "none")
# Use ggsave() to save the plot
ggsave(filename = paste0("Figures/ACC_EMMeanTimepointDriftRate_Automated.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
print(p)
```
## Draw Days Effects
```{r}
model = model_1
preds <- predict_response(model,
terms = c("Days"),
interval = "confidence",
margin = "mean_reference",
back.transform = FALSE)
preds = as.data.table(preds)
preds <- preds %>%
mutate(
sem = (conf.high - conf.low) / (2 * 1.96) # approximate SEM from 95% CI
)
preds <- preds %>%
mutate(
ci.low_sem = predicted - sem,
ci.high_sem = predicted + sem
)
filtered_preds_df = preds
# Create subject-centered data using Random Intercept subtraction
# Extract Random Intercepts
re_intercepts <- as.data.frame(ranef(model_1)$Subj)
re_intercepts$Subj <- rownames(re_intercepts)
intercept_col_idx <- grep("Intercept", names(re_intercepts))
re_intercepts <- re_intercepts[, c(ncol(re_intercepts), intercept_col_idx)]
names(re_intercepts) <- c("Subj", "RandomIntercept")
metrics_days_centered <- raw_metrics_full_data_days %>%
left_join(re_intercepts, by = "Subj") %>%
mutate(acc_general_centered = acc_general - RandomIntercept)
# Calculate intercepts and slopes for each subject based on raw data
# (Since the model only has random intercepts, we use individual OLS regressions
# to visualize the heterogeneity in slopes that exists in the raw data)
subj_trends <- metrics_days_centered %>%
group_by(Subj) %>%
summarise(
intercept = coef(lm(acc_general_centered ~ Days))[1],
slope = coef(lm(acc_general_centered ~ Days))[2]
) %>%
mutate(trend_color = ifelse(slope > 0, "green4", "red3"))
custom_theme <- theme_minimal(base_size = 20) +
theme(plot.title = element_text(size = rel(2.8), hjust = 0.5, face = "bold"), # Center the title
plot.background = element_blank(), # White background for the plot
panel.background = element_rect(fill = "white"), # White background for the panels
panel.grid.major = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank(), # Remove minor grid lines
axis.line = element_blank(), # Define the axis lines without enclosing
axis.ticks = element_line(color = "black"), # Define ticks to make them clear
axis.title.x = element_text(size = rel(2.2), margin = margin(t = 10), lineheight = 1.2), # Add margin on top of x-axis title
axis.title.y = element_text(size = rel(2.2), margin = margin(r = 10), lineheight = 1.2), # Add margin on the right of y-axis title
axis.text = element_text(size = rel(1.8)), # Slightly larger axis text
strip.text = element_text(size = rel(2.2)), # Make facet titles bigger
legend.position = "right", # Position the legend on the right
legend.title = element_blank(), # Customize the legend title size
legend.background = element_rect(color = "black", size = .5), # Add a box around the legend
legend.text = element_text(size = rel(1.5)), # Customize the legend text size
legend.spacing.y = unit(0.5, "cm"), # Adjust the spacing between legend items
panel.border = element_blank(), # Ensure no border around the panels
text = element_text(family = "Arial")) # Apply Times New Roman (or serif font)
p = ggplot(filtered_preds_df, aes(x = x, y = predicted)) +
# Add individual participant trend lines (subject-centered)
geom_abline(data = subj_trends, aes(intercept = intercept, slope = slope, group = Subj, color = trend_color),
alpha = 0.2, size = 0.5) +
scale_color_identity() +
geom_ribbon(aes(ymin = ci.low_sem, ymax = ci.high_sem), linetype = "dotted", linewidth = 1, fill = "gray", alpha = 0.2) +
geom_line(size = 1.5) +
#Add confidence interval
labs(title = "Days",y = "", x = "Days") +
custom_theme
# Use ggsave() to save the plot
ggsave(filename = paste0("Figures/ACC_DaysLine.png"), plot = p, dpi = 300, width = 12, height = 8, units = "in", bg="white")
print(p)
```
## Save Data
```{r}
save(combined_data_all, final_contrasts, subj_contrasts, raw_metrics_full_data_days, file = 'SavedOutputs/Raw_Metrics_Factor_Analysis.RData')
```

5_Alternative_Models_Raw.Rmd at commit dcc0852, no license · at the source

Overview

Authors: Brian Kim¹, John D. Medaglia^1,2

ORCID iDs: Brian Kim

Institutions: Drexel University (United States); University of Pennsylvania (United States)

Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1298

Dates: received 21 January 2026; accepted 4 June 2026; published online 14 July 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1162/imag.a.1298 · PMID 42459500 · PMCID PMC13370750 · OpenAlex W7165898388

Open access: diamond, a free copy (OpenAlex)

Status: code verified

Categories: other (modality), human (organism), cognitive (subfield)

Methods: Statistics, fMRI & imaging

Keywords: cognitive control, neuromodulation, TMS, LFPN/MFPN connectivity, DDM, individual differences, latent factor analysis

Topic: Transcranial Magnetic Stimulation Studies (Neurology, Neuroscience), according to OpenAlex

Funding: NIH (1-DP5-OD-021352-01); Starfish Neuroscience

Citations: not cited yet (Europe PMC); 74 references in the paper

Abstract

Transcranial Magnetic Stimulation (TMS) is a promising tool to probe and enhance cognitive control, yet effects are often inconsistent across individuals. These inconsistencies may arise from individual differences in Lateral Frontoparietal (Control) Network (L-FPN) and Medial Frontoparietal (Default) Network (M-FPN) interactions, essential for suppressing internal distraction and facilitating cognitive control. We tested whether baseline connectivity between the L-FPN and M-FPN moderates TMS outcomes in cognitive control tasks. We used the generalized drift rate as a task-general behavioral index of cognitive control, as it overcomes the reliability and interpretability limitations of standard difference measures. Participants completed inhibition (Stroop), working memory (n-back), and flexibility (Navon) tasks before and after intermittent theta-burst stimulation (iTBS) to the L-FPN, dorsal frontoparietal (attention) network (D-FPN), or cranial vertex. Stimulation targets were defined using individualized resting-state parcellations to maximize precision. We found that baseline connectivity moderated stimulation outcomes: individuals with more integrated L-FPN/M-FPN networks benefited most from L-FPN stimulation, whereas those with more segregated networks showed greater improvements from D-FPN stimulation. Notably, stimulation of a single site did not uniformly enhance performance, underscoring the importance of individual network profiles. These findings highlight L-FPN/M-FPN interactions as a basis for individualized TMS interventions and the utility of combining precise network mapping with robust behavioral modeling to optimize neuromodulation interventions.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 5 matches between paragraphs and lines of code.

CogNeW/project_L-FPN_M-FPN_cc_stim

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: dcc0852cf6df646a29ebb848890f5b2d251860a2, 7 July 2026

Languages: R (7)

Size: 1,587 files, 7 scripts

Software Heritage: not archived

Found in: “Data and Code Availability”

Holds: README, 7 notebooks

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: tidyverse (7 files), ggplot2 (6 files), data.table (4 files), lmerTest (4 files), emmeans (3 files), lavaan (2 files), psych (2 files), patchwork (1 file), reshape2 (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

8 files

1_Preprocess_Files.Rmd, R, 806 lines
2_Extract_DDM_Parameters.Rmd, R, 200 lines
3_DDMAnalysis.Rmd, R, 1,532 lines, 1 match
4_DDM_Validation.Rmd, R, 309 lines, 1 match
5_Alternative_Models_Raw.Rmd, R, 2,660 lines, 2 matches
6_Alternative_Model_DDM_Difference.Rmd, R, 665 lines
7_Alternative_Model_MRI_SCRUBBED_DDM_Analysis.Rmd, R, 99 lines, 1 match
README.md, Text, 42 lines

At the source: github.com/CogNeW/project_L-FPN_M-FPN_cc_stim

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
7 scripts, each with its path and the digest of its content;
5 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data and Code Availability

The R Code and raw data files required to reproduce the analyses reported are available at: https://github.com/CogNeW/project_L-FPN_M-FPN_cc_stim

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 2 authors, 7 keywords, 2 funders, 74 references.

Cite

This paper

Kim, B., & Medaglia, J. D. (2026). Frontoparietal control-default mode connectivity predicts TMS effects on cognitive control. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1298. https://doi.org/10.1162/imag.a.1298

BibTeX

@article{kim2026frontoparietal, author = {Kim, Brian and Medaglia, John D.}, title = {{Frontoparietal control-default mode connectivity predicts TMS effects on cognitive control}}, journal = {Imaging neuroscience (Cambridge, Mass.)}, year = {2026}, month = jul, volume = {4}, pages = {IMAG.a.1298}, publisher = {MIT Press}, issn = {2837-6056}, doi = {10.1162/imag.a.1298}, url = {https://doi.org/10.1162/imag.a.1298}, pmid = {42459500}, pmcid = {PMC13370750} }

RIS

TY - JOUR AU - Kim, Brian AU - Medaglia, John D. TI - Frontoparietal control-default mode connectivity predicts TMS effects on cognitive control T2 - Imaging neuroscience (Cambridge, Mass.) J2 - Imaging Neurosci (Camb) PY - 2026 DA - 2026/07/14 VL - 4 SP - IMAG.a.1298 SN - 2837-6056 PB - MIT Press DO - 10.1162/imag.a.1298 UR - https://doi.org/10.1162/imag.a.1298 LA - en ER -

CSL-JSON

{ "id": "10.1162/imag.a.1298", "type": "article-journal", "title": "Frontoparietal control-default mode connectivity predicts TMS effects on cognitive control", "container-title": "Imaging neuroscience (Cambridge, Mass.)", "author": [ { "family": "Kim", "given": "Brian" }, { "family": "Medaglia", "given": "John D." } ], "container-title-short": "Imaging Neurosci (Camb)", "volume": "4", "page": "IMAG.a.1298", "DOI": "10.1162/imag.a.1298", "PMID": "42459500", "PMCID": "PMC13370750", "ISSN": "2837-6056", "publisher": "MIT Press", "URL": "https://doi.org/10.1162/imag.a.1298", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 14 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.