Frontoparietal physics and multiple demand differentiation analyses
Sam Maione
06-10-24
For the main analyses in this script: (1) Computes the possible interaction between matching tasks and preferred function (e.g. spatial working memory voxels match the task they were taken from (spWM) and ‘prefer’ the hard condition compared to the easy condition), (2) computes the overlap in the top 10% of activated voxels from each estimated network, and (3) compares behavioral accuracy from each task.
For the exploratory analyses in this script, we (1) evaluate the effect size as a function of voxel selection size (section 0.0.6) and (2) evaluate the effect size and (3) overlap within multiple demand parcels not used in the main analyses (non-frontoparietal MD ROIs; section 0.0.7).
0.0.1 Read in data
neural_data = readRDS(here("inputs/neural_data.rds"))
dice_data = readRDS(here("inputs/dice_data.rds"))
behavioral_data = readRDS(here("inputs/behavioral_data.rds"))
voxel_selections = readRDS(here("inputs/neural_data_by_voxel-selection.rds"))
all_MD = readRDS(here("inputs/additional_ROIs.rds"))
all_MD = all_MD[!grepl("_bilateral", all_MD$roi), ]
manyregions = read.csv(here("inputs/manyregions_info.csv"), header=TRUE)
dice_data_additional = readRDS(here("inputs/dice_data_additional_ROIs.rds"))0.0.2 Is there functional overlap in the top 10% of MD and IP frontoparietal voxels?
With this analysis, we will study whether MD ROIs are engaged in physical reasoning, and whether IP ROIs are engaged in multiple demand. To conduct this analysis, we will extract and average the responses of both sets of ROIs to the physics task (physics > social), and the MD task (hard > easy), per subject. For each ROI (MD and IP ROIs will be analyzed separately), we will model its responses as predicted by (1) which task was used to extract its responses (MD or physics), (2) the extraction condition (either the preferred or dispreferred condition for each task; preferred conditions are “hard” and “physics”; dispreferred conditions are “easy” and “social”), and (3) their interaction.
This analysis allows us to test, in each ROI, whether the responses are domain-specific (leading to an interaction, such that the MD ROIs are more responsive to the hard > easy than the physics > social contrast, and the IP ROIs are more responsive to the physics > social than the hard > easy contrast), or domain-general (leading to a main effect of extraction condition, preferred > dispreferred, regardless of task). Each region will be classified as domain-specific (a significant interaction as described above), or domain-general (no interaction, and a significant main effect of condition, as described above).
functional_dataframe = data.frame()
for(voxel_amount in unique(neural_data$top_voxel_selection_method)){
for(fROI in unique(neural_data$roi)){
for(selection in unique(neural_data$selection_contrast)){
interaction_subset = neural_data[which(neural_data$roi == fROI & neural_data$selection_contrast == selection & neural_data$top_voxel_selection_method == voxel_amount),]
interaction_model = lmer(scale(mean_topvoxels_extracted_cope) ~ pref_dispref * matching_tasks + (1|subjectID), data = interaction_subset)
pref_sig = summary(lsmeans(interaction_model, pairwise ~ pref_dispref | matching_tasks, adjust = "none")$contrasts)
new_row = (c(
fROI,
voxel_amount,
selection,
#interaction
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',3],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',4],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',5],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',1],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',2],
#mismatch
pref_sig[1, "df"],
pref_sig[1, "t.ratio"],
pref_sig[1, "p.value"],
pref_sig[1, "estimate"],
pref_sig[1, "SE"],
#match
pref_sig[2, "df"],
pref_sig[2, "t.ratio"],
pref_sig[2, "p.value"],
pref_sig[2, "estimate"],
pref_sig[2, "SE"]))
functional_dataframe = rbind(functional_dataframe, new_row)
}
}
}
names(functional_dataframe) = c('fROI',
'voxel_percentage',
'selection_contrast',
'interaction_df',
'interaction_t',
'interaction_p',
'interation_effectsize',
'interation_effectsize_SE',
'mismatch_pref_df',
'mismatch_pref_t',
'mismatch_pref_p',
'mismatch_pref_std_estimate',
'mismatch_pref_SE',
'match_pref_df',
'match_pref_t',
'match_pref_p',
'match_pref_std_estimate',
'match_pref_SE')
#convert all to numbers, where available
functional_dataframe[] <- lapply(functional_dataframe, function(x) if (all(!is.na(as.numeric(x)), na.rm = TRUE)) as.numeric(x) else x)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercionfunctional_dataframe = functional_dataframe %>%
mutate(interaction_p_star = as.factor(
case_when(interaction_p < .001 ~ "***",
interaction_p < .01 ~ "**",
interaction_p < .05 ~ "*",
interaction_p < .1 ~ "~",
TRUE ~ " ")
))
functional_dataframe = functional_dataframe %>%
mutate(mismatch_pref_p_star = as.factor(
case_when(mismatch_pref_p < .001 ~ "***",
mismatch_pref_p < .01 ~ "**",
mismatch_pref_p < .05 ~ "*",
mismatch_pref_p < .1 ~ "~",
TRUE ~ " ")
))
functional_dataframe = functional_dataframe %>%
mutate(match_pref_p_star = as.factor(
case_when(match_pref_p < .001 ~ "***",
match_pref_p < .01 ~ "**",
match_pref_p < .05 ~ "*",
match_pref_p < .1 ~ "~",
TRUE ~ " ")
))
functional_dataframe = functional_dataframe %>%
mutate(across(where(is.numeric), ~round(., 3)))
functional_dataframe = functional_dataframe %>%
mutate(fROI = factor(fROI, levels = c('PPC_L', 'MPC_L', 'APC_L', 'PC_L', 'PPC_R', 'MPC_R', 'APC_R', 'PC_R'))) %>%
arrange(fROI)
voxelwise_extraction = functional_dataframe
DT::datatable(voxelwise_extraction %>%
arrange(voxel_percentage, selection_contrast),
options = list(scrollX = TRUE, pageLength = 32))mini_graph_data =
summarySEwithin(
data = neural_data,
measurevar = "mean_topvoxels_extracted_cope",
withinvars = c("task", "roi", "selection_contrast", "extraction_contrast"),
idvar = "subjectID"
)## Automatically converting the following non-factors to factors: task, roi## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------mini_graph_data$selection_contrast = recode_factor(mini_graph_data$selection_contrast, 'SWM' = 'WM ROIs', 'IP' = 'Physics ROIs')
mini_graph_data$roi = factor(mini_graph_data$roi, levels = c('PPC_L', 'MPC_L', 'APC_L', 'PC_L', 'PPC_R', 'MPC_R', 'APC_R', 'PC_R'))
mini_graph_data$selection_contrast = factor(mini_graph_data$selection_contrast, levels = c('WM ROIs', 'Physics ROIs'))
mini_graph = ggplot(mini_graph_data, aes(x = extraction_contrast, y = mean_topvoxels_extracted_cope_norm, fill = extraction_contrast))+
geom_bar(stat="identity", position=position_dodge(.7), width = 0.7) + theme_bw(12) +
geom_errorbar(width = .3, position=position_dodge(.7), aes(ymin = mean_topvoxels_extracted_cope_norm-se, ymax=mean_topvoxels_extracted_cope_norm+se)) +
guides(fill = guide_legend(title = "Condition ( > Rest)")) + theme(legend.position = "bottom") +
xlab('Functional ROI') +
ylab('Amplitude of Response') +
scale_fill_discrete(labels = c('Easy', 'Hard', 'Physics', 'Social')) +
scale_fill_manual(values=c('#ffab40','#FCD299', '#4a86e8', '#A4DBE8')) +
coord_cartesian(ylim = c(0, 15)) +
scale_y_continuous(position = "right") +
facet_grid(selection_contrast ~ roi, switch = "y") +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
theme(panel.grid = element_blank())## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.mini_graph
0.0.3 Is there voxel overlap (Dice’s coefficient) in top 10% of MD and IP frontoparietal voxels?
We will also take the same ROIs and compute, for each person, for each fROI, Dice’s coefficient. Dice’s coefficient is a measure of the extent of overlap between two regions. If there is at least moderate (.4 and above) overlap between all IP and MD fROIs, this would support the hypothesis that the two networks have shared function. If there is low overlap (.4 and below) between all IP and MD fROIs, this would support the two regions being functionally distinct. If some fROIs show at least moderate overlap while others show low overlap, this would support the hypothesis that some regions within these networks share function, whereas other regions are functionally distinct.
dice_summary =
summarySEwithin(
data = dice_data,
measurevar = "overlap_dices_coef",
withinvars = c("roi", "combined", "top_voxel_percent"),
idvar = "participantID")## Automatically converting the following non-factors to factors: roi, combined, top_voxel_percentdice_summary_10 = dice_summary[which(dice_summary$top_voxel_percent == '10' & dice_summary$combined == "phys_gt_soc vs. hard_gt_easy"),]
dice_participant = dice_data[which(dice_data$top_voxel_percent == '10' & dice_data$combined == "phys_gt_soc vs. hard_gt_easy"),]
dice_graph = ggplot(dice_summary_10, aes(x = roi, y = overlap_dices_coef_norm, fill = combined))+
geom_bar(stat="identity", position="dodge", width = 0.7, color="black", fill="white") +
theme(legend.position="none") +
ggtitle("Dice's Coefficients per Region (Hard > Easy vs. Physics > Social)") + xlab('Functional ROI') +
geom_errorbar(width = .3, position=position_dodge(.7), aes(ymin = overlap_dices_coef_norm-se, ymax=overlap_dices_coef_norm+se)) +
ylab("Dice's Coefficient")
participant_dice_graph = dice_graph + geom_point(data = dice_summary, aes(x = roi, y = overlap_dices_coef, fill = combined), alpha = 0.25)
participant_dice_graph
0.0.4 How does the voxel overlap (Dice’s coefficient) change with voxel selection in the frontoparietal ROIs?
Instead of only choosing the top 10% of voxels in the IP and MD localized regions, does the overlap grow stronger as the percentage of voxels used grows?
dice_summary =
summarySEwithin(
data = dice_data,
measurevar = "overlap_dices_coef",
withinvars = c("roi", "combined", "top_voxel_percent"),
idvar = "participantID")## Automatically converting the following non-factors to factors: roi, combined, top_voxel_percentdice_summary$top_voxel_percent = as.numeric(dice_summary$top_voxel_percent)
ggplot(dice_summary, aes(x = top_voxel_percent, y = overlap_dices_coef_norm, color = roi)) +
geom_point(shape = 1) + #geom_smooth(method = lm, se = FALSE) +
scale_x_continuous()
0.0.5 Are these tasks (DOTS and spWM) difficulty matched?
One possible explanation for the observed MD regions not showing function in the DOTS task is that the DOTS task could be too easy for the MD network to be involved (since MD activity scales with task difficulty). However, if DOTS accuracy and (hard) spWM accuracy are similar, this would provide evidence against this alternative explanation (which does seem to be the case).
behavioral_data = subset(behavioral_data, select = -X)
behavioral_data$Condition = factor(behavioral_data$Condition, levels = c('DOTS_Physics_Accuracy', 'DOTS_Social_Accuracy', 'spWM_Hard_Accuracy', 'spWM_Easy_Accuracy'))
behavioral_data$Condition = relevel(behavioral_data$Condition, ref = "DOTS_Physics_Accuracy")
behavioral_model = lmer(Percentage_accuracy ~ Condition + (1|Subject), data = behavioral_data)
lsmeans(behavioral_model, pairwise ~ Condition, adjust = "none")$contrasts## contrast estimate SE df t.ratio
## DOTS_Physics_Accuracy - DOTS_Social_Accuracy -0.01084 0.0182 186 -0.597
## DOTS_Physics_Accuracy - spWM_Hard_Accuracy -0.00221 0.0182 186 -0.122
## DOTS_Physics_Accuracy - spWM_Easy_Accuracy -0.14646 0.0182 186 -8.067
## DOTS_Social_Accuracy - spWM_Hard_Accuracy 0.00864 0.0182 186 0.476
## DOTS_Social_Accuracy - spWM_Easy_Accuracy -0.13561 0.0182 186 -7.470
## spWM_Hard_Accuracy - spWM_Easy_Accuracy -0.14425 0.0182 186 -7.945
## p.value
## 0.5510
## 0.9034
## <.0001
## 0.6348
## <.0001
## <.0001
##
## Degrees-of-freedom method: kenward-rogerDT::datatable(behavioral_data %>%
arrange(Subject),
options = list(scrollX = TRUE, pageLength = 16))0.0.6 How does the effect size of the interaction change with the number of voxels used?
The stronger the interaction, the more preference for one function over the other. If this effect remains relatively constant throughout voxel selection sizes, this provides greater evidence for this being a true effect.
functional_dataframe = data.frame()
for(voxel_amount in unique(voxel_selections$top_voxel_selection_method)){
for(fROI in unique(voxel_selections$roi)){
for(selection in unique(voxel_selections$selection_contrast)){
interaction_subset = voxel_selections[which(voxel_selections$roi == fROI & voxel_selections$selection_contrast == selection & voxel_selections$top_voxel_selection_method == voxel_amount),]
interaction_subset$pref_dispref = relevel(interaction_subset$pref_dispref, ref = "Preferred")
interaction_model = lmer(scale(mean_topvoxels_extracted_cope) ~ pref_dispref * matching_tasks + (1|subjectID), data = interaction_subset)
pref_sig = summary(lsmeans(interaction_model, pairwise ~ pref_dispref | matching_tasks, adjust = "none")$contrasts)
new_row = (c(
fROI,
voxel_amount,
selection,
#interaction
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',3],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',4],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',5],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',1],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',2],
#mismatch
pref_sig[1, "df"],
pref_sig[1, "t.ratio"],
pref_sig[1, "p.value"],
pref_sig[1, "estimate"],
pref_sig[1, "SE"],
#match
pref_sig[2, "df"],
pref_sig[2, "t.ratio"],
pref_sig[2, "p.value"],
pref_sig[2, "estimate"],
pref_sig[2, "SE"]))
functional_dataframe = rbind(functional_dataframe, new_row)
}
}
}
names(functional_dataframe) = c('fROI',
'voxel_percentage',
'selection_contrast',
'interaction_df',
'interaction_t',
'interaction_p',
'interation_effectsize',
'interation_effectsize_SE',
'mismatch_pref_df',
'mismatch_pref_t',
'mismatch_pref_p',
'mismatch_pref_std_estimate',
'mismatch_pref_SE',
'match_pref_df',
'match_pref_t',
'match_pref_p',
'match_pref_std_estimate',
'match_pref_SE')
#convert all to numbers, where available
functional_dataframe[] <- lapply(functional_dataframe, function(x) if (all(!is.na(as.numeric(x)), na.rm = TRUE)) as.numeric(x) else x)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercionfunctional_dataframe = functional_dataframe %>%
mutate(interaction_p_star = as.factor(
case_when(interaction_p < .001 ~ "***",
interaction_p < .01 ~ "**",
interaction_p < .05 ~ "*",
interaction_p < .1 ~ "~",
TRUE ~ " ")
))
functional_dataframe = functional_dataframe %>%
mutate(mismatch_pref_p_star = as.factor(
case_when(mismatch_pref_p < .001 ~ "***",
mismatch_pref_p < .01 ~ "**",
mismatch_pref_p < .05 ~ "*",
mismatch_pref_p < .1 ~ "~",
TRUE ~ " ")
))
functional_dataframe = functional_dataframe %>%
mutate(match_pref_p_star = as.factor(
case_when(match_pref_p < .001 ~ "***",
match_pref_p < .01 ~ "**",
match_pref_p < .05 ~ "*",
match_pref_p < .1 ~ "~",
TRUE ~ " ")
))
functional_dataframe = functional_dataframe %>%
mutate(across(where(is.numeric), ~round(., 3)))
functional_dataframe = functional_dataframe %>%
mutate(fROI = factor(fROI, levels = c('PPC_L', 'MPC_L', 'APC_L', 'PC_L', 'PPC_R', 'MPC_R', 'APC_R', 'PC_R'))) %>%
arrange(fROI)
voxelwise_extraction = functional_dataframe
DT::datatable(voxelwise_extraction %>%
arrange(voxel_percentage, selection_contrast),
options = list(scrollX = TRUE, pageLength = 32))functional_dataframe = voxelwise_extraction
functional_dataframe$voxel_percentage = as.numeric(functional_dataframe$voxel_percentage)
effectdf = functional_dataframe %>% pivot_longer(cols = c('mismatch_pref_std_estimate', 'match_pref_std_estimate'),
names_to = 'matching',
values_to = 'std_estimate')
effectdf$extraction_task = 'blank'
effectdf[which(effectdf$matching == 'match_pref_std_estimate' & effectdf$selection_contrast == 'IP'),]$extraction_task = 'DOTS'
effectdf[which(effectdf$matching == 'mismatch_pref_std_estimate' & effectdf$selection_contrast == 'IP'),]$extraction_task = 'spWM'
effectdf[which(effectdf$matching == 'match_pref_std_estimate' & effectdf$selection_contrast == 'SWM'),]$extraction_task = 'spWM'
effectdf[which(effectdf$matching == 'mismatch_pref_std_estimate' & effectdf$selection_contrast == 'SWM'),]$extraction_task = 'DOTS'
sep_by_vox_pairwise = ggplot(effectdf, aes(x = voxel_percentage, y = std_estimate, color = extraction_task, group = extraction_task))+
geom_line(aes(group = matching)) + geom_point() +
facet_grid(. ~ selection_contrast + fROI) +
scale_x_continuous(limits = c(1, 50)) +
xlab('Functional ROI') +
ylab('Standardized Pairwise (Prefer > Disprefer) Effect') + theme_bw() +
guides(color = guide_legend(title = "Extraction Task")) + theme(legend.position = "bottom") +
scale_color_manual(values=c('#4a86e8', '#ffab40'))
sep_by_vox_pairwise## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 32 rows containing missing values or values outside the scale range
## (`geom_point()`).
#I think the takeaway is that IP selections, when mismatched task, are greater than SWM selections with mismatched tasks0.0.7 Is there functional overlap in the MD and IP outside of the originally selected frontoparietal voxels?
The greater the effect size for preferred task > dispreferred condition (e.g. phys & hard > soc & easy), the greater evidence that a given ROI supports both intuitive physics AND spatial working memory function. It appears that there is greater disparity between the spWM extraction effect size in SWM voxels and DOTS extraction effect size in SWM voxels, when compared to the disparity between the spWM extraction effect size and DOTS extraction effect size in IP voxels. This suggests, when taken in conjunction with the statistics listed in the next dataframe, that the IP network–even outside the frontoparietal regions–may not be selective, unlike the MD network, which seems to be mostly selective for SWM across the whole brain.
functional_dataframe = data.frame()
for(voxel_amount in unique(all_MD$top_voxel_selection_method)){
for(fROI in unique(all_MD$roi)){
for(selection in unique(all_MD$selection_contrast)){
interaction_subset = all_MD[which(all_MD$roi == fROI & all_MD$selection_contrast == selection & all_MD$top_voxel_selection_method == voxel_amount),]
interaction_subset$pref_dispref = relevel(interaction_subset$pref_dispref, ref = "Preferred")
interaction_model = lmer(scale(mean_topvoxels_extracted_cope) ~ pref_dispref * matching_tasks + (1|subjectID), data = interaction_subset)
pref_sig = summary(lsmeans(interaction_model, pairwise ~ pref_dispref | matching_tasks, adjust = "none")$contrasts)
new_row = (c(
fROI,
voxel_amount,
selection,
#interaction
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',3],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',4],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',5],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',1],
summary(interaction_model)$coefficients['pref_dispref1:matching_tasks1',2],
#mismatch
pref_sig[1, "df"],
pref_sig[1, "t.ratio"],
pref_sig[1, "p.value"],
pref_sig[1, "estimate"],
pref_sig[1, "SE"],
#match
pref_sig[2, "df"],
pref_sig[2, "t.ratio"],
pref_sig[2, "p.value"],
pref_sig[2, "estimate"],
pref_sig[2, "SE"]))
functional_dataframe = rbind(functional_dataframe, new_row)
}
}
}
names(functional_dataframe) = c('fROI',
'voxel_percentage',
'selection_contrast',
'interaction_df',
'interaction_t',
'interaction_p',
'interation_effectsize',
'interation_effectsize_SE',
'mismatch_pref_df',
'mismatch_pref_t',
'mismatch_pref_p',
'mismatch_pref_std_estimate',
'mismatch_pref_SE',
'match_pref_df',
'match_pref_t',
'match_pref_p',
'match_pref_std_estimate',
'match_pref_SE')
#convert all to numbers, where available
functional_dataframe[] <- lapply(functional_dataframe, function(x) if (all(!is.na(as.numeric(x)), na.rm = TRUE)) as.numeric(x) else x)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercionfunctional_dataframe = functional_dataframe %>%
mutate(interaction_p_star = as.factor(
case_when(interaction_p < .001 ~ "***",
interaction_p < .01 ~ "**",
interaction_p < .05 ~ "*",
interaction_p < .1 ~ "~",
TRUE ~ " ")
))
functional_dataframe = functional_dataframe %>%
mutate(mismatch_pref_p_star = as.factor(
case_when(mismatch_pref_p < .001 ~ "***",
mismatch_pref_p < .01 ~ "**",
mismatch_pref_p < .05 ~ "*",
mismatch_pref_p < .1 ~ "~",
TRUE ~ " ")
))
functional_dataframe = functional_dataframe %>%
mutate(match_pref_p_star = as.factor(
case_when(match_pref_p < .001 ~ "***",
match_pref_p < .01 ~ "**",
match_pref_p < .05 ~ "*",
match_pref_p < .1 ~ "~",
TRUE ~ " ")
))
functional_dataframe = functional_dataframe %>%
mutate(across(where(is.numeric), ~round(., 3)))
roi_names_functional_dataframe_no_dice = merge(functional_dataframe, manyregions, by.x = "fROI", by.y = "ROI_name", all.x = TRUE)
roi_names_functional_dataframe_no_dice = transform(roi_names_functional_dataframe_no_dice, fROI = ifelse(is.na(ROI_name_final), fROI, ROI_name_final))
roi_names_functional_dataframe_no_dice = subset(roi_names_functional_dataframe_no_dice, select = -c(ROI_name_final, ROI_category, bilateral, focal_region, manyregions_region, old_ROI, parcel_overlaps_with))
roi_names_functional_dataframe_no_dice$fROI = ifelse(roi_names_functional_dataframe_no_dice$fROI %in% c("midParietal_L",
"midParietal_R",
"postParietal_L",
"postParietal_R",
"antParietal_L",
"antParietal_R",
"precentral_A_preCG_L",
"precentral_B_IFGop_L",
"precentral_A_preCG_R",
"precentral_B_IFGop_R"),
paste(roi_names_functional_dataframe_no_dice$fROI, "_*", sep = ""),
roi_names_functional_dataframe_no_dice$fROI)
indices_to_drop = grepl("_\\*$", roi_names_functional_dataframe_no_dice$fROI)
roi_names_functional_dataframe_no_dice = roi_names_functional_dataframe_no_dice[!indices_to_drop, ]
roi_names_functional_dataframe_no_dice$voxel_percentage = as.character(roi_names_functional_dataframe_no_dice$voxel_percentage)
roi_names_functional_dataframe_no_dice$voxel_percentage = as.integer(roi_names_functional_dataframe_no_dice$voxel_percentage)mini_graph_data =
summarySEwithin(
data = all_MD,
measurevar = "mean_topvoxels_extracted_cope",
withinvars = c("task", "roi", "selection_contrast", "extraction_contrast"),
idvar = "subjectID"
)## Automatically converting the following non-factors to factors: task, roimini_graph_data$selection_contrast = recode_factor(mini_graph_data$selection_contrast, 'SWM' = 'WM ROIs', 'IP' = 'Physics ROIs')
#mini_graph_data$roi = factor(mini_graph_data$roi, levels = c('PPC_L', 'MPC_L', 'APC_L', 'PC_L', 'PPC_R', 'MPC_R', 'APC_R', 'PC_R'))
mini_graph_data$selection_contrast = factor(mini_graph_data$selection_contrast, levels = c('WM ROIs', 'Physics ROIs'))
mini_graph = ggplot(mini_graph_data, aes(x = extraction_contrast, y = mean_topvoxels_extracted_cope_norm, fill = extraction_contrast))+
geom_bar(stat="identity", position=position_dodge(.7), width = 0.7) + theme_bw(12) +
geom_errorbar(width = .3, position=position_dodge(.7), aes(ymin = mean_topvoxels_extracted_cope_norm-se, ymax=mean_topvoxels_extracted_cope_norm+se)) +
guides(fill = guide_legend(title = "Condition ( > Rest)")) + theme(legend.position = "bottom") +
xlab('Functional ROI') +
ylab('Amplitude of Response') +
scale_fill_discrete(labels = c('Easy', 'Hard', 'Physics', 'Social')) +
scale_fill_manual(values=c('#ffab40','#FCD299', '#4a86e8', '#A4DBE8')) +
coord_cartesian(ylim = c(0, 15)) +
scale_y_continuous(position = "right") +
facet_grid(selection_contrast ~ roi, switch = "y") +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
theme(panel.grid = element_blank())## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.mini_graph
roi_names_functional_dataframe_no_dice$voxel_percentage = as.numeric(roi_names_functional_dataframe_no_dice$voxel_percentage)
effectdf = roi_names_functional_dataframe_no_dice %>% pivot_longer(cols = c('mismatch_pref_std_estimate', 'match_pref_std_estimate'),
names_to = 'matching',
values_to = 'std_estimate')
effectdf$extraction_task = 'blank'
effectdf[which(effectdf$matching == 'match_pref_std_estimate' & effectdf$selection_contrast == 'IP'),]$extraction_task = 'DOTS'
effectdf[which(effectdf$matching == 'mismatch_pref_std_estimate' & effectdf$selection_contrast == 'IP'),]$extraction_task = 'spWM'
effectdf[which(effectdf$matching == 'match_pref_std_estimate' & effectdf$selection_contrast == 'SWM'),]$extraction_task = 'spWM'
effectdf[which(effectdf$matching == 'mismatch_pref_std_estimate' & effectdf$selection_contrast == 'SWM'),]$extraction_task = 'DOTS'
sep_by_vox_pairwise = ggplot(effectdf, aes(x = voxel_percentage, y = std_estimate, color = extraction_task, group = extraction_task))+
geom_line(aes(group = matching)) + geom_point() +
facet_grid(. ~ selection_contrast + fROI) +
scale_x_continuous(limits = c(1, 50)) +
xlab('Functional ROI') +
ylab('Standardized Pairwise (Prefer > Disprefer) Effect') + theme_bw() +
guides(color = guide_legend(title = "Extraction Task")) + theme(legend.position = "bottom") +
scale_color_manual(values=c('#4a86e8', '#ffab40'))
sep_by_vox_pairwise
DT::datatable(roi_names_functional_dataframe_no_dice %>%
arrange(voxel_percentage, selection_contrast),
options = list(scrollX = TRUE, pageLength = 32))0.0.8 How does the overlap in MD and IP voxels change, for ROIs outside original frontoparietal selections?
As the percentage of voxels selected for increases, overlap in terms of Dice’s coefficient increases, for ROIs outside of the main frontoparietal selection.
dice_data_additional_phys_hard = dice_data_additional[which(dice_data_additional$combined == "phys_gt_soc vs. hard_gt_easy"), ]
roi_names_functional_dataframe = merge(dice_data_additional_phys_hard, manyregions, by.x = c("roi"), by.y = c("ROI_name"), all.x = TRUE)
roi_names_functional_dataframe = transform(roi_names_functional_dataframe, roi = ifelse(is.na(ROI_name_final), roi, ROI_name_final))
roi_names_functional_dataframe = subset(roi_names_functional_dataframe, select = -c(ROI_name_final, ROI_category, bilateral, focal_region, manyregions_region, old_ROI, parcel_overlaps_with))
roi_names_functional_dataframe$roi = ifelse(roi_names_functional_dataframe$roi %in% c("midParietal_L", "midParietal_R", "postParietal_L", "postParietal_R", "antParietal_L",
"antParietal_R", "precentral_A_preCG_L", "precentral_B_IFGop_L", "precentral_A_preCG_R",
"precentral_B_IFGop_R"),
paste(roi_names_functional_dataframe$roi, "_*", sep = ""),
roi_names_functional_dataframe$roi)
roi_names_functional_dataframe = merge(roi_names_functional_dataframe, roi_names_functional_dataframe_no_dice, by.x = c("roi", "top_voxel_percent"), by.y = c("fROI", "voxel_percentage"), all = TRUE)
indices_to_drop = grepl("_\\*$", roi_names_functional_dataframe$roi)
roi_names_functional_dataframe_new = roi_names_functional_dataframe[!indices_to_drop, ]
roi_names_functional_dataframe_new$top_voxel_percent = as.character(roi_names_functional_dataframe_new$top_voxel_percent)
roi_names_functional_dataframe_new$top_voxel_percent = as.integer(roi_names_functional_dataframe_new$top_voxel_percent)
ggplot(roi_names_functional_dataframe_new, aes(x = top_voxel_percent, y = overlap_dices_coef_norm, color = roi)) +
geom_point(shape = 1) + #geom_smooth(method = lm, se = FALSE) +
scale_x_continuous()