OSCR

Fusiform face area development correlates with development in higher-order social brain regions.

Code ↔ Paper

10 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 10 matches
  1. [1] § Materials and methods › fMRI data analysis › Motion treatment ↔ scripts/05.motion_exclusions/mark_motion_exclusions.py, lines 142–217 · score 0.80 · composite motion, framewise displacement, standardised DVARS, rapidart, thresholds, workflow
  2. [2] § Results › Development of FFA face response and functional connectivity ↔ scripts/FRIC_DevelopmentalChange.Rmd, lines 1614–1713 · score 0.76 · left STS, left amygdala, left MMPFC, right STS, right amygdala, right MMPFC
  3. [3] § Results › Development of FFA face response and functional connectivity ↔ scripts/FRIC_FunctionalMaturityFFA.Rmd, lines 1209–1300 · score 0.75 · left STS, left amygdala, left MMPFC, right STS, right amygdala, right MMPFC
  4. [4] § Materials and methods › fMRI data analysis › Region of interest (ROI) definition ↔ scripts/06.first_level/define_fROIs.py, lines 113–201 · score 0.71 · ROI definition, fROI, search space, ranking, maps, voxel
  5. [5] § Materials and methods › Statistical analyses › Developmental change in functional responses of FFA, MMPFC, amygdala and STS ↔ scripts/FRIC_DevelopmentalChange.Rmd, lines 697–838 · score 0.58 · scene events, face events, S01, S12, interaction, hemisphere
  6. [6] § Materials and methods › Statistical analyses › Associations between functional maturity of FFA and its functional connectivity to MMPFC, amygdala and STS ↔ scripts/FRIC_FunctionalMaturityFFA.Rmd, lines 304–429 · score 0.56 · preregistered linear mixed, independent variable, functional maturity, reversed, right FFA, interaction
  7. [7] § Materials and methods › fMRI data analysis › Inter-region correlations (i.e., functional connectivity) ↔ scripts/FRIC_MRIvariables.Rmd, lines 717–776 · score 0.53 · inter region correlations, Pearson correlations, ROIs, STS, MMPFC, FFA
  8. [8] § Materials and methods › fMRI data analysis › Timecourse extraction ↔ scripts/06.first_level/firstlevel_pipeline.py, lines 205–257 · score 0.52 · aCompCor, outlier volumes, regressed, filtered
  9. [9] § Materials and methods › fMRI data analysis › Region of interest (ROI) definition ↔ scripts/06.first_level/extract_stats.py, lines 20–63 · score 0.51 · localiser task, beliefs, fROI, split, events, model
  10. [10] § Materials and methods › fMRI data analysis › Functional maturity ↔ scripts/FRIC_MRIvariables.Rmd, lines 209–309 · score 0.51 · Pearson correlation, fROIs, scored, child, adult, STS

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 authors' code

R Markdown · 2,474 lines · 81 KB · no license · 2 matches

  1. ---
  2. title: "FRIC_DevelopmentalChange"
  3. author: "LJS"
  4. date: "2025-01-14"
  5. output: html_document
  6. ---
  7. # Libraries
  8. ```{r setup, include=FALSE}
  9. library("rlang") #toolbox
  10. library("readr") #toolbox
  11. library("dplyr") #data manipulation
  12. library("purrr") #data manipulation
  13. library("tidyr") #data manipulation
  14. library("here") #toolbox
  15. library("tidyverse") #data manipulation
  16. library("ggpubr") #plotting
  17. library("tidytext") #data manipulation
  18. library("knitr") #rmarkdown
  19. library("htmltools") #rmarkdown
  20. library("markdown") #rmarkdown
  21. library("httpuv") #rmarkdown
  22. library("NCmisc") #toolbox
  23. library("stringr") #data manipulation
  24. library("reshape") #data manipulation
  25. library("reshape2") #data manipulation
  26. library("lme4") #regression models
  27. library("lmerTest") #regression models
  28. library("broom.mixed") #save regression results
  29. library("car") #VIF
  30. library("VIM") #visualise missing data
  31. library("gridExtra") #combine plots
  32. library("patchwork") #combine plots
  33. ```
  34. # Face Responses in Childhood (FRIC) project
  35. This script characterises age-related changes in functional maturity, response magnitude to face/scene events and functional connectivity of FFA, MMPFC, STS and amygdala, as well as age-related changes in lateralisation of face response in FFA.
  36. Structure:
  37. - Section I: Functional maturity measure (fm) i.e., similarity between children’s and adults’ timecourses’
  38. - Section II: Magnitude of responses of each ROI at face/scene events
  39. - Section III: Age effects on functional connectivity (fx) between regions
  40. - Section IV: Age effects on lateralisation of face response (non-preregistered)
  41. ## Section I: Functional maturity measure (fm) i.e., similarity between children’s and adults’ timecourses’
  42. ### Read and tidy up data
  43. ```{r}
  44. # Covariates
  45. participants <- read_tsv(here("raw_data/covariates", "participants.tsv"), col_names=TRUE)
  46. colnames(participants)[which(names(participants) == "participant_id")] <- "ID"
  47. ## Remove the word "pixar" and trim extra whitespace
  48. participants$ID <- gsub("[sub-pixar-]","", participants$ID)
  49. participants$ID <- trimws(participants$ID)
  50. # Motion
  51. motion <- read_tsv(here("raw_data/motion", "outlier_info.tsv"), col_names=TRUE)
  52. colnames(motion)[which(names(motion) == "subject")] <- "ID"
  53. ## Remove the word "pixar" and trim extra whitespace
  54. motion$ID <- gsub("[sub-pixar-]","", motion$ID)
  55. motion$ID <- trimws(motion$ID)
  56. # Fm
  57. child_fm <- read_csv(here("processed_data", "child_MRIvariables.csv"), col_names=TRUE)
  58. child_fm <- child_fm[,c(1:9)]
  59. ## Reshape the dataset to long format
  60. long_data <- child_fm %>%
  61. pivot_longer(cols = starts_with("zscored_fm_"),
  62. names_to = c("hemisphere", "region"),
  63. names_pattern = "zscored_fm_(l|r)(.*)", # Capture hemisphere and region
  64. values_to = "fm") %>%
  65. mutate(hemisphere = ifelse(hemisphere == "l", "left", "right")) # Convert l/r to left/right
  66. # Merge age, motion and fm data
  67. fm_data <- merge(long_data, participants[,c(1:2)], by="ID")
  68. fm_data <- merge(fm_data, motion[,c(4,6)], by="ID")
  69. colnames(fm_data) <- c("ID", "hemisphere", "region", "fm", "age", "meanFD")
  70. # Ensure categorical variables are treated as factors
  71. fm_data$ID <- as.factor(fm_data$ID)
  72. fm_data$hemisphere <- as.factor(fm_data$hemisphere)
  73. fm_data$region <- as.factor(fm_data$region)
  74. # Scale predictors
  75. fm_data[, c("age", "meanFD")] <- scale(fm_data[, c("age", "meanFD")])
  76. ```
  77. ### Models
  78. With interaction effects
  79. ```{r}
  80. # FFA
  81. fm_data_FFA <- subset(fm_data, region == "FFA")
  82. model_FFA <- lmer(fm ~ age*hemisphere + meanFD + (1 | ID), data = fm_data_FFA)
  83. summary_FFA <- tidy(model_FFA)
  84. # MMPFC
  85. fm_data_MMPFC <- subset(fm_data, region == "MMPFC")
  86. model_MMPFC <- lmer(fm ~ age*hemisphere + meanFD + (1 | ID), data = fm_data_MMPFC)
  87. summary_MMPFC <- tidy(model_MMPFC)
  88. # Amygdala
  89. fm_data_amygdala <- subset(fm_data, region == "amygdala")
  90. model_amygdala <- lmer(fm ~ age*hemisphere + meanFD + (1 | ID), data = fm_data_amygdala)
  91. summary_amygdala <- tidy(model_amygdala)
  92. # STS
  93. fm_data_STS <- subset(fm_data, region == "STS")
  94. model_STS <- lmer(fm ~ age*hemisphere + meanFD + (1 | ID), data = fm_data_STS)
  95. summary_STS <- tidy(model_STS)
  96. # Summary
  97. all_summaries <- bind_rows(
  98. summary_FFA %>% mutate(region = "FFA"),
  99. summary_MMPFC %>% mutate(region = "MMPFC"),
  100. summary_amygdala %>% mutate(region = "amygdala"),
  101. summary_STS %>% mutate(region = "STS")
  102. )
  103. view(all_summaries)
  104. ```
  105. Since interaction effects are non-sig, repeat models without interaction effects
  106. ```{r}
  107. # FFA
  108. fm_data_FFA <- subset(fm_data, region == "FFA")
  109. model_FFA <- lmer(fm ~ age + hemisphere + meanFD + (1 | ID), data = fm_data_FFA)
  110. summary_FFA <- tidy(model_FFA)
  111. # MMPFC
  112. fm_data_MMPFC <- subset(fm_data, region == "MMPFC")
  113. model_MMPFC <- lmer(fm ~ age + hemisphere + meanFD + (1 | ID), data = fm_data_MMPFC)
  114. summary_MMPFC <- tidy(model_MMPFC)
  115. # Amygdala
  116. fm_data_amygdala <- subset(fm_data, region == "amygdala")
  117. model_amygdala <- lmer(fm ~ age + hemisphere + meanFD + (1 | ID), data = fm_data_amygdala)
  118. summary_amygdala <- tidy(model_amygdala)
  119. # STS
  120. fm_data_STS <- subset(fm_data, region == "STS")
  121. model_STS <- lmer(fm ~ age + hemisphere + meanFD + (1 | ID), data = fm_data_STS)
  122. summary_STS <- tidy(model_STS)
  123. # Save results
  124. all_summaries <- bind_rows(
  125. summary_FFA %>% mutate(region = "FFA"),
  126. summary_MMPFC %>% mutate(region = "MMPFC"),
  127. summary_amygdala %>% mutate(region = "amygdala"),
  128. summary_STS %>% mutate(region = "STS")
  129. ) %>%
  130. mutate(
  131. estimate = round(estimate, 2),
  132. std.error = round(std.error, 2),
  133. statistic = round(statistic, 2),
  134. df = round(df, 0),
  135. p.value = if_else(
  136. p.value < 0.001,
  137. format(p.value, scientific = TRUE, digits = 4),
  138. sprintf("%.3f", p.value)
  139. )
  140. )
  141. write.csv(all_summaries, here("results/RQ1", "Fm_age.csv"), row.names = FALSE)
  142. ```
  143. ### Model plots
  144. Age and functional maturity
  145. ```{r}
  146. # Data wrangling
  147. data_plot <- fm_data
  148. data_plot$age <- NULL
  149. data_plot$meanFD <- NULL
  150. data_plot <- merge(data_plot, participants[,c(1:2)], by="ID")
  151. colnames(data_plot) <- c("ID", "hemisphere", "region", "fm", "age")
  152. data_plot_FFA <- subset(data_plot, region == "FFA")
  153. data_plot_STS <- subset(data_plot, region == "STS")
  154. data_plot_MMPFC <- subset(data_plot, region == "MMPFC")
  155. data_plot_amygdala <- subset(data_plot, region == "amygdala")
  156. # Colours
  157. custom_colours_amygdala <- c("#E1AA9F", "#D34835")
  158. custom_colours_MMPFC <- c("#87BACF", "#578CAD")
  159. custom_colours_FFA <- c("#C39EDA", "#652975")
  160. custom_colours_STS <- c("#FFF865", "#FFCE1B")
  161. # FFA plot
  162. plot_FFA <- ggplot(data_plot_FFA, aes(x = age, y = fm, color = hemisphere)) +
  163. geom_point(alpha = 0.8, size = 2.5, shape = 19) +
  164. geom_smooth(method = "lm", se = FALSE, linetype = "solid", linewidth = 2) +
  165. scale_color_manual(values = custom_colours_FFA) +
  166. labs(title = "",
  167. x = "Age",
  168. y = "",
  169. color = "FFA Hemisphere") +
  170. theme_minimal() +
  171. theme(
  172. legend.position = "bottom",
  173. legend.text = element_text(size = 12),
  174. legend.title = element_text(size = 12, face = "bold"),
  175. axis.text = element_text(size = 12),
  176. axis.title = element_text(size = 14, face = "bold"),
  177. strip.text = element_text(size = 14, face = "bold"),
  178. strip.background = element_blank(),
  179. axis.line.x = element_blank(),
  180. axis.text.y = element_blank(),
  181. plot.title = element_text(hjust = 0.5, size = 14, face = "bold")) +
  182. coord_cartesian(ylim = c(-3, 3)) +
  183. annotate(
  184. "text",
  185. x = Inf, y = -Inf,
  186. label = "Age effect: \nβ(SE)=0.26(0.07), p<0.001",
  187. hjust = 1.05,
  188. vjust = -0.5,
  189. size = 4,
  190. color = "black"
  191. )
  192. # MMPFC plot
  193. plot_MMPFC <- ggplot(data_plot_MMPFC, aes(x = age, y = fm, color = hemisphere)) +
  194. geom_point(alpha = 0.8, size = 2.8, shape = 19) +
  195. geom_smooth(method = "lm", se = FALSE, linetype = "solid", linewidth = 2) +
  196. scale_color_manual(values = custom_colours_MMPFC) +
  197. labs(title = "",
  198. x = "Age",
  199. y = "Functional maturity",
  200. color = "MMPFC Hemisphere") +
  201. theme_minimal() +
  202. theme(
  203. legend.position = "bottom",
  204. legend.text = element_text(size = 12),
  205. legend.title = element_text(size = 12, face = "bold"),
  206. axis.text = element_text(size = 12),
  207. axis.title = element_text(size = 14, face = "bold"),
  208. strip.text = element_text(size = 14, face = "bold"),
  209. strip.background = element_blank(),
  210. axis.line.x = element_blank(),
  211. plot.title = element_text(hjust = 0.5, size = 14, face = "bold")) +
  212. coord_cartesian(ylim = c(-3, 3)) +
  213. annotate(
  214. "text",
  215. x = Inf, y = -Inf,
  216. label = "Age effect: \nβ(SE)=0.35(0.07), p<0.001",
  217. hjust = 1.05,
  218. vjust = -0.5,
  219. size = 4,
  220. color = "black"
  221. )
  222. # Amygdala plot
  223. plot_amygdala <- ggplot(data_plot_amygdala, aes(x = age, y = fm, color = hemisphere)) +
  224. geom_point(alpha = 0.8, size = 2.8, shape = 19) +
  225. geom_smooth(method = "lm", se = FALSE, linetype = "solid", linewidth = 2) +
  226. scale_color_manual(values = custom_colours_amygdala) +
  227. labs(
  228. title = "",
  229. x = "Age",
  230. y = "",
  231. color = "Amygdala Hemisphere"
  232. ) +
  233. theme_minimal() +
  234. theme(
  235. legend.position = "bottom",
  236. legend.text = element_text(size = 12),
  237. legend.title = element_text(size = 12, face = "bold"),
  238. axis.text = element_text(size = 12),
  239. axis.title = element_text(size = 14, face = "bold"),
  240. strip.text = element_text(size = 14, face = "bold"),
  241. strip.background = element_blank(),
  242. axis.line.x = element_blank(),
  243. axis.text.y = element_blank(),
  244. plot.title = element_text(hjust = 0.5, size = 14, face = "bold")
  245. ) +
  246. coord_cartesian(ylim = c(-3, 3)) +
  247. annotate(
  248. "text",
  249. x = Inf, y = -Inf,
  250. label = "Age effect: \nβ(SE)=0.20(0.08), p=0.016",
  251. hjust = 1.05,
  252. vjust = -0.5,
  253. size = 4,
  254. color = "black"
  255. )
  256. # STS plot
  257. plot_STS <- ggplot(data_plot_STS, aes(x = age, y = fm, color = hemisphere)) +
  258. geom_point(alpha = 0.8, size = 2.5, shape = 19) +
  259. geom_smooth(method = "lm", se = FALSE, linetype = "solid", linewidth = 2) +
  260. scale_color_manual(values = custom_colours_STS) +
  261. labs(title = "",
  262. x = "Age",
  263. y = "",
  264. color = "STS Hemisphere") +
  265. theme_minimal() +
  266. theme(
  267. legend.position = "bottom",
  268. legend.text = element_text(size = 12),
  269. legend.title = element_text(size = 12, face = "bold"),
  270. axis.text = element_text(size = 12),
  271. axis.title = element_text(size = 14, face = "bold"),
  272. strip.text = element_text(size = 14, face = "bold"),
  273. strip.background = element_blank(),
  274. axis.line.x = element_blank(),
  275. axis.text.y = element_blank(),
  276. plot.title = element_text(hjust = 0.5, size = 14, face = "bold")) +
  277. coord_cartesian(ylim = c(-3, 3)) +
  278. annotate(
  279. "text",
  280. x = Inf, y = -Inf,
  281. label = "Age effect: \nβ(SE)=0.24(0.07), p<0.001",
  282. hjust = 1.05,
  283. vjust = -0.5,
  284. size = 4,
  285. color = "black"
  286. )
  287. # Combine and save plots
  288. fm_age <- (plot_MMPFC | plot_amygdala | plot_STS | plot_FFA ) +
  289. plot_annotation(title = "Functional maturity across Age") &
  290. theme(plot.title = element_text(hjust = 0.5, size = 16, face = "bold"))
  291. fm_age
  292. ggsave(here("results/figures", "Fm_age.png"), plot = fm_age, width = 13, height = 5, units = "in", dpi = 300)
  293. ```
  294. ## Section II: Magnitude of responses of each ROI at face/scene events
  295. ### Read data
  296. ```{r}
  297. # Events magnitude
  298. child_eventsmagnitude <- read_csv(here("processed_data", "child_MRIvariables.csv"), col_names=TRUE)
  299. child_eventsmagnitude <- child_eventsmagnitude[,c(1, 10:201)]
  300. # Reshape the dataset to long format
  301. long_data <- child_eventsmagnitude %>%
  302. pivot_longer(
  303. cols = -ID, # Exclude the ID column
  304. names_to = c("region", "event"), # Split names into region and event
  305. names_pattern = "([lr]?[^_]+)_(F\\d{2}|S\\d{2})", # Match region and event
  306. values_to = "value" # The name of the new column for values
  307. ) %>%
  308. mutate(
  309. hemisphere = if_else(grepl("^l", region), "left", "right"), # Add hemisphere based on region prefix
  310. region = gsub("^[lr]", "", region) # Remove hemisphere prefix from region
  311. ) %>%
  312. pivot_wider(
  313. names_from = "event", # Make "event" column headers
  314. values_from = "value" # Fill the values for each event
  315. )
  316. long_data <- long_data %>%
  317. select(ID, hemisphere, region, F01, F02, F03, F04, F05, F06, F07, F08, F09, F10, F11, F12,
  318. S01, S02, S03, S04, S05, S06, S07, S08, S09, S10, S11, S12)
  319. # Merge age, motion and fm data
  320. eventsmagnitude_data <- merge(long_data, participants[,c(1:2)], by="ID")
  321. eventsmagnitude_data <- merge(eventsmagnitude_data, motion[,c(4,6)], by="ID")
  322. colnames(eventsmagnitude_data) <- c("ID", "hemisphere", "region", "F01", "F02", "F03", "F04", "F05", "F06", "F07", "F08", "F09", "F10", "F11", "F12", "S01", "S02", "S03", "S04", "S05", "S06", "S07", "S08", "S09", "S10", "S11", "S12", "age", "meanFD")
  323. # Scale predictors
  324. eventsmagnitude_data[, c("age", "meanFD")] <- scale(eventsmagnitude_data[, c("age", "meanFD")])
  325. # Ensure categorical variables are treated as factors
  326. eventsmagnitude_data$ID <- as.factor(eventsmagnitude_data$ID)
  327. eventsmagnitude_data$hemisphere <- as.factor(eventsmagnitude_data$hemisphere)
  328. eventsmagnitude_data$region <- as.factor(eventsmagnitude_data$region)
  329. ```
  330. ### Models
  331. FFA
  332. ```{r}
  333. # Read data
  334. eventsmagnitude_data_FFA <- subset(eventsmagnitude_data, region == "FFA")
  335. ## Face events
  336. ### Models with interaction effects
  337. events <- paste0("F", sprintf("%02d", 1:12))
  338. model_summaries <- list()
  339. for (event in events) {
  340. #### Create the formula dynamically
  341. formula <- as.formula(paste(event, "~ age*hemisphere + meanFD + (1 | ID)"))
  342. #### Fit the model
  343. model <- lmer(formula, data = eventsmagnitude_data_FFA)
  344. #### Summarise the model
  345. summary_tidy <- tidy(model)
  346. #### Store the summary in the list with the event name as the key
  347. model_summaries[[event]] <- summary_tidy
  348. }
  349. summaries_FFA_face_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  350. model_summaries[[event]] %>% mutate(event = event) # Add the event name to each summary
  351. }))
  352. summaries_FFA_face_events$region <- "FFA"
  353. summaries_FFA_face_events <- summaries_FFA_face_events %>%
  354. mutate(sig = case_when(
  355. p.value < 0.001 ~ "***",
  356. p.value < 0.01 ~ "**",
  357. p.value < 0.05 ~ "*",
  358. TRUE ~ ""
  359. ))
  360. ### Interaction effects are non-sig for all events, repeat models without interaction effects for all events
  361. model_summaries <- list()
  362. for (event in events) {
  363. #### Create the formula dynamically
  364. formula <- as.formula(paste(event, "~ age + hemisphere + meanFD + (1 | ID)"))
  365. #### Fit the model
  366. model <- lmer(formula, data = eventsmagnitude_data_FFA)
  367. #### Summarise the model
  368. summary_tidy <- tidy(model)
  369. #### Store the summary in the list with the event name as the key
  370. model_summaries[[event]] <- summary_tidy
  371. }
  372. summaries_FFA_face_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  373. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  374. }))
  375. summaries_FFA_face_events$region <- "FFA"
  376. summaries_FFA_face_events$p.value <- as.numeric(summaries_FFA_face_events$p.value)
  377. summaries_FFA_face_events <- summaries_FFA_face_events %>%
  378. mutate(sig = case_when(
  379. p.value < 0.004 ~ "significant after correction",
  380. TRUE ~ ""
  381. ))
  382. ## Scene events
  383. ### Models with interaction effects
  384. events <- paste0("S", sprintf("%02d", 1:12)) #generates S01, S02, ..., S12
  385. model_summaries <- list()
  386. for (event in events) {
  387. #### Create the formula dynamically
  388. formula <- as.formula(paste(event, "~ age*hemisphere + meanFD + (1 | ID)"))
  389. #### Fit the model
  390. model <- lmer(formula, data = eventsmagnitude_data_FFA)
  391. #### Summarise the model
  392. summary_tidy <- tidy(model)
  393. #### Store the summary in the list with the event name as the key
  394. model_summaries[[event]] <- summary_tidy
  395. }
  396. summaries_FFA_scene_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  397. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  398. }))
  399. summaries_FFA_scene_events$region <- "FFA"
  400. summaries_FFA_scene_events <- summaries_FFA_scene_events %>%
  401. mutate(sig = case_when(
  402. p.value < 0.001 ~ "***",
  403. p.value < 0.01 ~ "**",
  404. p.value < 0.05 ~ "*",
  405. TRUE ~ ""
  406. ))
  407. ### Interaction effects are non-sig for all events, repeat models without interaction effects for all events
  408. model_summaries <- list()
  409. for (event in events) {
  410. #### Create the formula dynamically
  411. formula <- as.formula(paste(event, "~ age + hemisphere + meanFD + (1 | ID)"))
  412. #### Fit the model
  413. model <- lmer(formula, data = eventsmagnitude_data_FFA)
  414. #### Summarise the model
  415. summary_tidy <- tidy(model)
  416. #### Store the summary in the list with the event name as the key
  417. model_summaries[[event]] <- summary_tidy
  418. }
  419. summaries_FFA_scene_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  420. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  421. }))
  422. summaries_FFA_scene_events$region <- "FFA"
  423. summaries_FFA_scene_events$p.value <- as.numeric(summaries_FFA_scene_events$p.value)
  424. summaries_FFA_scene_events <- summaries_FFA_scene_events %>%
  425. mutate(sig = case_when(
  426. p.value < 0.004 ~ "significant after correction",
  427. TRUE ~ ""
  428. ))
  429. ```
  430. MMPFC
  431. ```{r}
  432. # Read data
  433. eventsmagnitude_data_MMPFC <- subset(eventsmagnitude_data, region == "MMPFC")
  434. ## Face events
  435. ### Models with interaction effects
  436. events <- paste0("F", sprintf("%02d", 1:12))
  437. model_summaries <- list()
  438. for (event in events) {
  439. #### Create the formula dynamically
  440. formula <- as.formula(paste(event, "~ age*hemisphere + meanFD + (1 | ID)"))
  441. #### Fit the model
  442. model <- lmer(formula, data = eventsmagnitude_data_MMPFC)
  443. #### Summarise the model
  444. summary_tidy <- tidy(model)
  445. #### Store the summary in the list with the event name as the key
  446. model_summaries[[event]] <- summary_tidy
  447. }
  448. summaries_MMPFC_face_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  449. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  450. }))
  451. summaries_MMPFC_face_events$region <- "MMPFC"
  452. summaries_MMPFC_face_events <- summaries_MMPFC_face_events %>%
  453. mutate(sig = case_when(
  454. p.value < 0.001 ~ "***",
  455. p.value < 0.01 ~ "**",
  456. p.value < 0.05 ~ "*",
  457. TRUE ~ ""
  458. ))
  459. ### Interaction effects are non-sig for all events but F05, F07, F10 and F12, repeat models without interaction effects for all events but F05, F06, F07, F10 and F12
  460. events <- c("F01", "F02", "F03", "F04", "F06", "F08", "F09", "F11")
  461. model_summaries <- list()
  462. for (event in events) {
  463. #### Create the formula dynamically
  464. formula <- as.formula(paste(event, "~ age + hemisphere + meanFD + (1 | ID)"))
  465. #### Fit the model
  466. model <- lmer(formula, data = eventsmagnitude_data_MMPFC)
  467. #### Summarise the model
  468. summary_tidy <- tidy(model)
  469. #### Store the summary in the list with the event name as the key
  470. model_summaries[[event]] <- summary_tidy
  471. }
  472. summaries_MMPFC_face_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  473. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  474. }))
  475. model_MMPFC_F05 <- lmer(F05 ~ age*hemisphere + meanFD + (1 | ID), data = eventsmagnitude_data_MMPFC)
  476. summary_MMPFC_F05 <- tidy(model_MMPFC_F05)
  477. summary_MMPFC_F05 <- summary_MMPFC_F05 %>% mutate(event = "F05")
  478. model_MMPFC_F07 <- lmer(F07 ~ age*hemisphere + meanFD + (1 | ID), data = eventsmagnitude_data_MMPFC)
  479. summary_MMPFC_F07 <- tidy(model_MMPFC_F07)
  480. summary_MMPFC_F07 <- summary_MMPFC_F07 %>% mutate(event = "F07")
  481. model_MMPFC_F10 <- lmer(F10 ~ age*hemisphere + meanFD + (1 | ID), data = eventsmagnitude_data_MMPFC)
  482. summary_MMPFC_F10 <- tidy(model_MMPFC_F10)
  483. summary_MMPFC_F10 <- summary_MMPFC_F10 %>% mutate(event = "F10")
  484. model_MMPFC_F12 <- lmer(F12 ~ age*hemisphere + meanFD + (1 | ID), data = eventsmagnitude_data_MMPFC)
  485. summary_MMPFC_F12 <- tidy(model_MMPFC_F12)
  486. summary_MMPFC_F12 <- summary_MMPFC_F12 %>% mutate(event = "F12")
  487. summaries_MMPFC_face_events <- bind_rows(summaries_MMPFC_face_events, summary_MMPFC_F05, summary_MMPFC_F07, summary_MMPFC_F10, summary_MMPFC_F12)
  488. summaries_MMPFC_face_events$region <- "MMPFC"
  489. summaries_MMPFC_face_events$p.value <- as.numeric(summaries_MMPFC_face_events$p.value)
  490. summaries_MMPFC_face_events <- summaries_MMPFC_face_events %>%
  491. mutate(sig = case_when(
  492. p.value < 0.004 ~ "significant after correction",
  493. TRUE ~ ""
  494. ))
  495. ## Scene events
  496. ### Models with interaction effects
  497. events <- paste0("S", sprintf("%02d", 1:12)) #generates S01, S02, ..., S12
  498. model_summaries <- list()
  499. for (event in events) {
  500. #### Create the formula dynamically
  501. formula <- as.formula(paste(event, "~ age*hemisphere + meanFD + (1 | ID)"))
  502. #### Fit the model
  503. model <- lmer(formula, data = eventsmagnitude_data_MMPFC)
  504. #### Summarise the model
  505. summary_tidy <- tidy(model)
  506. #### Store the summary in the list with the event name as the key
  507. model_summaries[[event]] <- summary_tidy
  508. }
  509. summaries_MMPFC_scene_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  510. model_summaries[[event]] %>% mutate(event = event) # Add the event name to each summary
  511. }))
  512. summaries_MMPFC_scene_events$region <- "MMPFC"
  513. summaries_MMPFC_scene_events <- summaries_MMPFC_scene_events %>%
  514. mutate(sig = case_when(
  515. p.value < 0.001 ~ "***",
  516. p.value < 0.01 ~ "**",
  517. p.value < 0.05 ~ "*",
  518. TRUE ~ ""
  519. ))
  520. ### Interaction effects are non-sig for all events, repeat models without interaction effects for all events
  521. model_summaries <- list()
  522. for (event in events) {
  523. #### Create the formula dynamically
  524. formula <- as.formula(paste(event, "~ age + hemisphere + meanFD + (1 | ID)"))
  525. #### Fit the model
  526. model <- lmer(formula, data = eventsmagnitude_data_MMPFC)
  527. #### Summarise the model
  528. summary_tidy <- tidy(model)
  529. #### Store the summary in the list with the event name as the key
  530. model_summaries[[event]] <- summary_tidy
  531. }
  532. summaries_MMPFC_scene_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  533. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  534. }))
  535. summaries_MMPFC_scene_events$region <- "MMPFC"
  536. summaries_MMPFC_scene_events$p.value <- as.numeric(summaries_MMPFC_scene_events$p.value)
  537. summaries_MMPFC_scene_events <- summaries_MMPFC_scene_events %>%
  538. mutate(sig = case_when(
  539. p.value < 0.004 ~ "significant after correction",
  540. TRUE ~ ""
  541. ))
  542. ```
  543. Amygdala
  544. ```{r}
  545. # Read data
  546. eventsmagnitude_data_amygdala <- subset(eventsmagnitude_data, region == "amygdala")
  547. ## Face events
  548. ### Models with interaction effects
  549. events <- paste0("F", sprintf("%02d", 1:12))
  550. model_summaries <- list()
  551. for (event in events) {
  552. #### Create the formula dynamically
  553. formula <- as.formula(paste(event, "~ age*hemisphere + meanFD + (1 | ID)"))
  554. #### Fit the model
  555. model <- lmer(formula, data = eventsmagnitude_data_amygdala)
  556. #### Summarise the model
  557. summary_tidy <- tidy(model)
  558. #### Store the summary in the list with the event name as the key
  559. model_summaries[[event]] <- summary_tidy
  560. }
  561. summaries_amygdala_face_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  562. model_summaries[[event]] %>% mutate(event = event) # Add the event name to each summary
  563. }))
  564. summaries_amygdala_face_events$region <- "amygdala"
  565. summaries_amygdala_face_events <- summaries_amygdala_face_events %>%
  566. mutate(sig = case_when(
  567. p.value < 0.001 ~ "***",
  568. p.value < 0.01 ~ "**",
  569. p.value < 0.05 ~ "*",
  570. TRUE ~ ""
  571. ))
  572. ### Interaction effects are non-sig for all events but F04, repeat models without interaction effects for all events but F04
  573. events <- c("F01", "F02", "F03", "F05", "F06", "F07", "F08", "F09", "F10", "F11", "F12")
  574. model_summaries <- list()
  575. for (event in events) {
  576. #### Create the formula dynamically
  577. formula <- as.formula(paste(event, "~ age + hemisphere + meanFD + (1 | ID)"))
  578. #### Fit the model
  579. model <- lmer(formula, data = eventsmagnitude_data_amygdala)
  580. #### Summarise the model
  581. summary_tidy <- tidy(model)
  582. #### Store the summary in the list with the event name as the key
  583. model_summaries[[event]] <- summary_tidy
  584. }
  585. summaries_amygdala_face_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  586. model_summaries[[event]] %>% mutate(event = event) # Add the event name to each summary
  587. }))
  588. model_amygdala_F04 <- lmer(F04 ~ age*hemisphere + meanFD + (1 | ID), data = eventsmagnitude_data_amygdala)
  589. summary_amygdala_F04 <- tidy(model_amygdala_F04)
  590. summary_amygdala_F04 <- summary_amygdala_F04 %>% mutate(event = "F04")
  591. summaries_amygdala_face_events <- bind_rows(summaries_amygdala_face_events, summary_amygdala_F04)
  592. summaries_amygdala_face_events$region <- "amygdala"
  593. summaries_amygdala_face_events$p.value <- as.numeric(summaries_amygdala_face_events$p.value)
  594. summaries_amygdala_face_events <- summaries_amygdala_face_events %>%
  595. mutate(sig = case_when(
  596. p.value < 0.004 ~ "significant after correction",
  597. TRUE ~ ""
  598. ))
  599. ## Scene events
  600. ### Models with interaction effects
  601. events <- paste0("S", sprintf("%02d", 1:12)) #generates S01, S02, ..., S12
  602. model_summaries <- list()
  603. for (event in events) {
  604. #### Create the formula dynamically
  605. formula <- as.formula(paste(event, "~ age*hemisphere + meanFD + (1 | ID)"))
  606. #### Fit the model
  607. model <- lmer(formula, data = eventsmagnitude_data_amygdala)
  608. #### Summarise the model
  609. summary_tidy <- tidy(model)
  610. #### Store the summary in the list with the event name as the key
  611. model_summaries[[event]] <- summary_tidy
  612. }
  613. summaries_amygdala_scene_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  614. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  615. }))
  616. summaries_amygdala_scene_events$region <- "amygdala"
  617. summaries_amygdala_scene_events <- summaries_amygdala_scene_events %>%
  618. mutate(sig = case_when(
  619. p.value < 0.001 ~ "***",
  620. p.value < 0.01 ~ "**",
  621. p.value < 0.05 ~ "*",
  622. TRUE ~ ""
  623. ))
  624. ### Interaction effects are non-sig for all events, repeat models without interaction effects for all events
  625. model_summaries <- list()
  626. for (event in events) {
  627. # Create the formula dynamically
  628. formula <- as.formula(paste(event, "~ age + hemisphere + meanFD + (1 | ID)"))
  629. # Fit the model
  630. model <- lmer(formula, data = eventsmagnitude_data_amygdala)
  631. # Summarise the model
  632. summary_tidy <- tidy(model)
  633. # Store the summary in the list with the event name as the key
  634. model_summaries[[event]] <- summary_tidy
  635. }
  636. summaries_amygdala_scene_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  637. model_summaries[[event]] %>% mutate(event = event) # Add the event name to each summary
  638. }))
  639. summaries_amygdala_scene_events$region <- "amygdala"
  640. summaries_amygdala_scene_events$p.value <- as.numeric(summaries_amygdala_scene_events$p.value)
  641. summaries_amygdala_scene_events <- summaries_amygdala_scene_events %>%
  642. mutate(sig = case_when(
  643. p.value < 0.004 ~ "significant after correction",
  644. TRUE ~ ""
  645. ))
  646. ```
  647. STS
  648. ```{r}
  649. # Read data
  650. eventsmagnitude_data_STS <- subset(eventsmagnitude_data, region == "STS")
  651. ## Face events
  652. ### Models with interaction effects
  653. events <- paste0("F", sprintf("%02d", 1:12))
  654. model_summaries <- list()
  655. for (event in events) {
  656. #### Create the formula dynamically
  657. formula <- as.formula(paste(event, "~ age*hemisphere + meanFD + (1 | ID)"))
  658. #### Fit the model
  659. model <- lmer(formula, data = eventsmagnitude_data_STS)
  660. #### Summarise the model
  661. summary_tidy <- tidy(model)
  662. #### Store the summary in the list with the event name as the key
  663. model_summaries[[event]] <- summary_tidy
  664. }
  665. summaries_STS_face_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  666. model_summaries[[event]] %>% mutate(event = event) # Add the event name to each summary
  667. }))
  668. summaries_STS_face_events$region <- "STS"
  669. summaries_STS_face_events <- summaries_STS_face_events %>%
  670. mutate(sig = case_when(
  671. p.value < 0.001 ~ "***",
  672. p.value < 0.01 ~ "**",
  673. p.value < 0.05 ~ "*",
  674. TRUE ~ ""
  675. ))
  676. ### Interaction effects are non-sig for all events but F09, repeat models without interaction effects for all events but F09
  677. events <- c("F01", "F02", "F03", "F04", "F05", "F06", "F07", "F08", "F10", "F11", "F12")
  678. model_summaries <- list()
  679. for (event in events) {
  680. #### Create the formula dynamically
  681. formula <- as.formula(paste(event, "~ age + hemisphere + meanFD + (1 | ID)"))
  682. #### Fit the model
  683. model <- lmer(formula, data = eventsmagnitude_data_STS)
  684. #### Summarise the model
  685. summary_tidy <- tidy(model)
  686. #### Store the summary in the list with the event name as the key
  687. model_summaries[[event]] <- summary_tidy
  688. }
  689. summaries_STS_face_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  690. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  691. }))
  692. model_STS_F09 <- lmer(F09 ~ age*hemisphere + meanFD + (1 | ID), data = eventsmagnitude_data_STS)
  693. summary_STS_F09 <- tidy(model_STS_F09)
  694. summary_STS_F09 <- summary_STS_F09 %>% mutate(event = "F09")
  695. summaries_STS_face_events <- bind_rows(summaries_STS_face_events, summary_STS_F09)
  696. summaries_STS_face_events$region <- "STS"
  697. summaries_STS_face_events$p.value <- as.numeric(summaries_STS_face_events$p.value)
  698. summaries_STS_face_events <- summaries_STS_face_events %>%
  699. mutate(sig = case_when(
  700. p.value < 0.004 ~ "significant after correction",
  701. TRUE ~ ""
  702. ))
  703. ## Scene events
  704. ### Models with interaction effects
  705. events <- paste0("S", sprintf("%02d", 1:12)) #generates S01, S02, ..., S12
  706. model_summaries <- list()
  707. for (event in events) {
  708. #### Create the formula dynamically
  709. formula <- as.formula(paste(event, "~ age*hemisphere + meanFD + (1 | ID)"))
  710. #### Fit the model
  711. model <- lmer(formula, data = eventsmagnitude_data_STS)
  712. #### Summarise the model
  713. summary_tidy <- tidy(model)
  714. #### Store the summary in the list with the event name as the key
  715. model_summaries[[event]] <- summary_tidy
  716. }
  717. summaries_STS_scene_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  718. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  719. }))
  720. summaries_STS_scene_events$region <- "STS"
  721. summaries_STS_scene_events <- summaries_STS_scene_events %>%
  722. mutate(sig = case_when(
  723. p.value < 0.001 ~ "***",
  724. p.value < 0.01 ~ "**",
  725. p.value < 0.05 ~ "*",
  726. TRUE ~ ""
  727. ))
  728. ### Interaction effects are non-sig for all events but S02 and S08, repeat models without interaction effects for all events except S02 and S08
  729. events <- c("S01", "S03", "S04", "S05", "S06", "S07", "S09", "S10", "S11", "S12")
  730. model_summaries <- list()
  731. for (event in events) {
  732. #### Create the formula dynamically
  733. formula <- as.formula(paste(event, "~ age + hemisphere + meanFD + (1 | ID)"))
  734. #### Fit the model
  735. model <- lmer(formula, data = eventsmagnitude_data_STS)
  736. #### Summarise the model
  737. summary_tidy <- tidy(model)
  738. #### Store the summary in the list with the event name as the key
  739. model_summaries[[event]] <- summary_tidy
  740. }
  741. summaries_STS_scene_events <- do.call(rbind, lapply(names(model_summaries), function(event) {
  742. model_summaries[[event]] %>% mutate(event = event) #add the event name to each summary
  743. }))
  744. model_STS_S02 <- lmer(S02 ~ age*hemisphere + meanFD + (1 | ID), data = eventsmagnitude_data_STS)
  745. summary_STS_S02 <- tidy(model_STS_S02)
  746. summary_STS_S02 <- summary_STS_S02 %>% mutate(event = "S02")
  747. model_STS_S08 <- lmer(S08 ~ age*hemisphere + meanFD + (1 | ID), data = eventsmagnitude_data_STS)
  748. summary_STS_S08 <- tidy(model_STS_S08)
  749. summary_STS_S08 <- summary_STS_S08 %>% mutate(event = "S08")
  750. summaries_STS_scene_events <- bind_rows(summaries_STS_scene_events, summary_STS_S02, summary_STS_S08)
  751. summaries_STS_scene_events$region <- "STS"
  752. summaries_STS_scene_events$p.value <- as.numeric(summaries_STS_scene_events$p.value)
  753. summaries_STS_scene_events <- summaries_STS_scene_events %>%
  754. mutate(sig = case_when(
  755. p.value < 0.004 ~ "significant after correction",
  756. TRUE ~ ""
  757. ))
  758. ```
  759. Save results
  760. ```{r}
  761. all_summaries <- bind_rows(
  762. summaries_FFA_face_events,
  763. summaries_FFA_scene_events,
  764. summaries_MMPFC_face_events,
  765. summaries_MMPFC_scene_events,
  766. summaries_amygdala_face_events,
  767. summaries_amygdala_scene_events,
  768. summaries_STS_face_events,
  769. summaries_STS_scene_events
  770. )
  771. all_summaries$p.value <- as.numeric(all_summaries$p.value)
  772. event_order <- c(
  773. "F01","F09","F11","F04","F10","F12","F07","F08","F06","F05","F02","F03",
  774. "S07","S02","S03","S04","S12","S11","S10","S06","S08","S05","S09","S01"
  775. )
  776. all_summaries <- all_summaries %>%
  777. mutate(
  778. across(c(estimate, std.error, statistic, df), ~ round(.x, 2)),
  779. p.value = if_else(p.value < 0.001,
  780. format(p.value, scientific = TRUE, digits = 4),
  781. sprintf("%.3f", p.value)),
  782. event = factor(event, levels = event_order)
  783. ) %>%
  784. arrange(region, event)
  785. write.csv(all_summaries, here("results/RQ1", "Eventmagnitude_age.csv"), row.names = FALSE)
  786. ```
  787. ### Model plots
  788. Timecourses and events
  789. ```{r}
  790. # Data wrangling
  791. TCs_total <- read_csv(here("processed_data", "TCs_total.csv"), col_names=TRUE)
  792. data_plot <- merge(TCs_total, participants[,c(1,3,4)], by="ID")
  793. data_plot <- data_plot %>%
  794. filter(Child_Adult == "child")
  795. data_plot$Child_Adult <- NULL
  796. data_plot <- data_plot %>%
  797. group_by(time, AgeGroup) %>%
  798. summarise(
  799. avg_lamygdala = mean(lamygdala, na.rm = TRUE),
  800. avg_ramygdala = mean(ramygdala, na.rm = TRUE),
  801. avg_lFFA = mean(lFFA, na.rm = TRUE),
  802. avg_rFFA = mean(rFFA, na.rm = TRUE),
  803. avg_lMMPFC = mean(lMMPFC, na.rm = TRUE),
  804. avg_rMMPFC = mean(rMMPFC, na.rm = TRUE),
  805. avg_lSTS = mean(lSTS, na.rm = TRUE),
  806. avg_rSTS = mean(rSTS, na.rm = TRUE), .groups = "drop"
  807. )
  808. avg_adult_TCs_trad_fROIs <- read.csv(here("raw_data/timecourses", "TCs_trad_adult_fROIs.csv"), header = TRUE, dec = ".", sep = ",")
  809. avg_adult_TCs_trad_fROIs <- avg_adult_TCs_trad_fROIs %>%
  810. mutate(time = time - 1)
  811. TCs_amygdala <- read.csv(here("raw_data/timecourses", "TCs_amygdala.csv"), header = TRUE, dec = ".", sep = ",")
  812. TCs_amygdala$sub <- gsub("[pixar]", "", TCs_amygdala$sub)
  813. TCs_amygdala$sub <- trimws(TCs_amygdala$sub)
  814. TCs_amygdala$run <- NULL
  815. colnames(TCs_amygdala) <- c("time", "ID", "lamygdala", "ramygdala")
  816. adult_list <- row.names(participants)[which(participants[["Child_Adult"]] == "adult")]
  817. adult_TCs_amygdala <- TCs_amygdala %>%
  818. filter(ID %in% adult_list)
  819. avg_adult_TCs_amygdala <- adult_TCs_amygdala %>%
  820. group_by(time) %>%
  821. summarise(
  822. avg_lamygdala = mean(lamygdala, na.rm = TRUE),
  823. avg_ramygdala = mean(ramygdala, na.rm = TRUE)
  824. )
  825. avg_adult_TCs <- merge(avg_adult_TCs_amygdala, avg_adult_TCs_trad_fROIs, by="time")
  826. avg_adult_TCs$AgeGroup <- "adult"
  827. avg_adult_TCs <- avg_adult_TCs[, c("time", "AgeGroup", "avg_lamygdala", "avg_ramygdala", "avg_lFFA", "avg_rFFA", "avg_lMMPFC", "avg_rMMPFC", "avg_lSTS", "avg_rSTS")]
  828. data_plot <- rbind(data_plot, avg_adult_TCs)
  829. colnames(data_plot) <- c("time", "AgeGroup", "lamygdala", "ramygdala", "lFFA", "rFFA", "lMMPFC", "rMMPFC", "lSTS", "rSTS")
  830. data_plot <- data_plot %>%
  831. pivot_longer(
  832. cols = c(lamygdala, ramygdala, lFFA, rFFA, lMMPFC, rMMPFC, lSTS, rSTS),
  833. names_to = "Temp",
  834. values_to = "value"
  835. ) %>%
  836. ## Separate the "Temp" column into "Hemisphere" and "Region"
  837. mutate(
  838. hemisphere = if_else(str_starts(Temp, "l"), "left", "right"),
  839. region = case_when(
  840. str_detect(Temp, "amygdala") ~ "amygdala",
  841. str_detect(Temp, "FFA") ~ "FFA",
  842. str_detect(Temp, "MMPFC") ~ "MMPFC",
  843. str_detect(Temp, "STS") ~ "STS",
  844. TRUE ~ NA_character_
  845. )
  846. ) %>%
  847. select(-Temp) #drop the temporary column
  848. data_plot_FFA <- subset(data_plot, region == "FFA")
  849. data_plot_MMPFC <- subset(data_plot, region == "MMPFC")
  850. data_plot_amygdala <- subset(data_plot, region == "amygdala")
  851. data_plot_STS <- subset(data_plot, region == "STS")
  852. events_TRs <- read_csv(here("processed_data", "events_TRs.csv"), col_names=TRUE)
  853. events_TRs <- events_TRs[,c(1:3)]
  854. custom_colours <- c("#D34835", "#DDA86A", "#83AC66", "#87BACF", "#6989F3", "#652975")
  855. # Figure in main text
  856. Plot_rFFA <- ggplot(
  857. subset(data_plot_FFA, hemisphere == "right")
  858. ) +
  859. geom_line(aes(x = time, y = value, color = AgeGroup, group = AgeGroup),
  860. size = 1.2, alpha = 1) +
  861. geom_hline(yintercept = 0, color = "grey", linewidth = 0.8) +
  862. geom_text(
  863. data = subset(
  864. data_plot_FFA,
  865. hemisphere == "right" &
  866. time %in% c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55,
  867. 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110,
  868. 115, 120, 125, 130, 135, 140, 145, 150, 155,
  869. 160, 165, 170)
  870. ),
  871. aes(x = time, y = -0.07, label = time),
  872. vjust = 0,
  873. size = 3,
  874. color = "black",
  875. angle = 90
  876. ) +
  877. theme_classic() +
  878. scale_color_manual(values = custom_colours) +
  879. scale_y_continuous(limits = c(-1.1, 1.1)) +
  880. scale_x_continuous(
  881. breaks = c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55,
  882. 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110,
  883. 115, 120, 125, 130, 135, 140, 145, 150, 155,
  884. 160, 165, 170),
  885. labels = NULL
  886. ) +
  887. labs(
  888. title = "Right FFA",
  889. y = "Response magnitude",
  890. color = "Age group"
  891. ) +
  892. theme(
  893. legend.position = "bottom",
  894. axis.title.x = element_text(size = 16),
  895. axis.title.y = element_text(size = 16),
  896. axis.text.x = element_text(size = 16),
  897. axis.text.y = element_text(size = 16),
  898. axis.ticks.x = element_blank(),
  899. axis.line.x = element_blank(),
  900. plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
  901. axis.line.y = element_line(linewidth = 0.8, color = "black"),
  902. legend.text = element_text(size = 16),
  903. legend.title = element_text(size = 16, face = "bold")
  904. )
  905. ## Add the shaded rectangles for events
  906. Plot_rFFA <- Plot_rFFA +
  907. geom_rect(data = events_TRs,
  908. aes(xmin = onset, xmax = end, ymin = -1.1, ymax = 1.1, fill = trial_type),
  909. alpha = 0.3, color = NA) +
  910. scale_fill_manual(values = c("faces" = "#F9B6C7", "scenes" = "#3B7D23")) +
  911. guides(fill = guide_legend(title = "Event type"), title.theme = element_text(face = "bold"))
  912. print(Plot_rFFA)
  913. ## Save plot
  914. ggsave(here("results/figures", "TCs_rFFA.png"), plot = Plot_rFFA, width = 15, height = 5, units = "in", dpi = 300)
  915. # Supplementary figure
  916. ## FFA Timecourses, including events
  917. ### Plot
  918. Plot_FFA <- ggplot(data_plot_FFA) +
  919. geom_line(aes(x = time, y = value, color = AgeGroup, group = AgeGroup), size = 1.2, alpha = 1) +
  920. geom_hline(yintercept = 0, color = "grey", linewidth = 0.8) +
  921. geom_text(data = data_plot_FFA[data_plot_FFA$time %in% c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170),],
  922. aes(x = time, y = -0.07, label = time), vjust = 0, size = 3, color = "black", angle = 90) +
  923. facet_wrap(~ hemisphere, scales = "free_y") +
  924. theme_classic() +
  925. scale_color_manual(values = custom_colours) +
  926. scale_y_continuous(limits = c(-1.1, 1.1)) +
  927. scale_x_continuous(
  928. breaks = c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170),
  929. labels = NULL
  930. ) +
  931. labs(
  932. title = "FFA",
  933. y = "Response magnitude",
  934. color = "Age group"
  935. ) +
  936. theme(
  937. strip.text = element_text(size = 16),
  938. strip.background = element_blank(),
  939. legend.position = "bottom",
  940. axis.title.x = element_text(size = 16),
  941. axis.title.y = element_text(size = 16),
  942. axis.text.x = element_text(size = 16),
  943. axis.text.y = element_text(size = 16),
  944. axis.ticks.x = element_blank(),
  945. axis.line.x = element_blank(),
  946. plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
  947. axis.line.y = element_line(linewidth = 0.8, color = "black"),
  948. legend.text = element_text(size = 16),
  949. legend.title = element_text(size = 16, face = "bold")
  950. )
  951. ### Add the shaded rectangles for events
  952. Plot_FFA <- Plot_FFA +
  953. geom_rect(data = events_TRs,
  954. aes(xmin = onset, xmax = end, ymin = -1.1, ymax = 1.1, fill = trial_type),
  955. alpha = 0.3, color = NA) +
  956. scale_fill_manual(values = c("faces" = "#F9B6C7", "scenes" = "#3B7D23")) +
  957. guides(fill = guide_legend(title = "Event type"), title.theme = element_text(face = "bold"))
  958. print(Plot_FFA)
  959. ## MMPFC Timecourses, including events
  960. Plot_MMPFC <- ggplot(data_plot_MMPFC) +
  961. geom_line(aes(x = time, y = value, color = AgeGroup, group = AgeGroup), size = 1.2, alpha = 1) +
  962. geom_hline(yintercept = 0, color = "grey", linewidth = 0.8) +
  963. geom_text(data = data_plot_MMPFC[data_plot_MMPFC$time %in% c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170),],
  964. aes(x = time, y = -0.07, label = time), vjust = 0, size = 3, color = "black", angle = 90) +
  965. facet_wrap(~ hemisphere, scales = "free_y") +
  966. theme_classic() +
  967. scale_color_manual(values = custom_colours) +
  968. scale_y_continuous(limits = c(-1.1, 1.1)) +
  969. scale_x_continuous(
  970. breaks = c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170),
  971. labels = NULL
  972. ) +
  973. labs(
  974. title = "MMPFC",
  975. y = "Response magnitude",
  976. color = "Age group"
  977. ) +
  978. theme(
  979. strip.text = element_text(size = 16),
  980. strip.background = element_blank(),
  981. legend.position = "bottom",
  982. axis.title.x = element_text(size = 16),
  983. axis.title.y = element_text(size = 16),
  984. axis.text.x = element_text(size = 16),
  985. axis.text.y = element_text(size = 16),
  986. axis.ticks.x = element_blank(),
  987. axis.line.x = element_blank(),
  988. plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
  989. axis.line.y = element_line(linewidth = 0.8, color = "black"),
  990. legend.text = element_text(size = 16),
  991. legend.title = element_text(size = 16, face = "bold")
  992. )
  993. ### Add the shaded rectangles for events
  994. Plot_MMPFC <- Plot_MMPFC +
  995. geom_rect(data = events_TRs,
  996. aes(xmin = onset, xmax = end, ymin = -1.1, ymax = 1.1, fill = trial_type),
  997. alpha = 0.3, color = NA) +
  998. scale_fill_manual(values = c("faces" = "#F9B6C7", "scenes" = "#3B7D23")) +
  999. guides(fill = guide_legend(title = "Event type"), title.theme = element_text(face = "bold"))
  1000. ### Get the 10th "faces" event (pink box)
  1001. pink_boxes <- events_TRs[events_TRs$trial_type == "faces",]
  1002. tenth_pink_box <- pink_boxes[10,] #get the 10th pink box
  1003. ### Adjust the y position for the asterisk to stay within the valid range
  1004. Plot_MMPFC <- Plot_MMPFC +
  1005. geom_text(data = tenth_pink_box,
  1006. aes(x = (onset + end) / 2 - 1, y = 1.05, label = "*"), # Adjust y to fit within plot limits
  1007. size = 8, color = "black", fontface = "bold", vjust = 0)
  1008. print(Plot_MMPFC)
  1009. ## Amygdala Timecourses, including events
  1010. Plot_amygdala <- ggplot(data_plot_amygdala) +
  1011. geom_line(aes(x = time, y = value, color = AgeGroup, group = AgeGroup), size = 1.2, alpha = 1) +
  1012. geom_hline(yintercept = 0, color = "grey", linewidth = 0.8) +
  1013. geom_text(data = data_plot_amygdala[data_plot_amygdala$time %in% c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170),],
  1014. aes(x = time, y = -0.07, label = time), vjust = 0, size = 3, color = "black", angle = 90) +
  1015. facet_wrap(~ hemisphere, scales = "free_y") +
  1016. theme_classic() +
  1017. scale_color_manual(values = custom_colours) +
  1018. scale_y_continuous(limits = c(-1.1, 1.1)) +
  1019. scale_x_continuous(
  1020. breaks = c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170),
  1021. labels = NULL
  1022. ) +
  1023. labs(
  1024. title = "Amygdala",
  1025. y = "Response magnitude",
  1026. color = "Age group"
  1027. ) +
  1028. theme(
  1029. strip.text = element_text(size = 16),
  1030. strip.background = element_blank(),
  1031. legend.position = "bottom",
  1032. axis.title.x = element_text(size = 16),
  1033. axis.title.y = element_text(size = 16),
  1034. axis.text.x = element_text(size = 16),
  1035. axis.text.y = element_text(size = 16),
  1036. axis.ticks.x = element_blank(),
  1037. axis.line.x = element_blank(),
  1038. plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
  1039. axis.line.y = element_line(linewidth = 0.8, color = "black"),
  1040. legend.text = element_text(size = 16),
  1041. legend.title = element_text(size = 16, face = "bold")
  1042. )
  1043. ### Add the shaded rectangles for events
  1044. Plot_amygdala <- Plot_amygdala +
  1045. geom_rect(data = events_TRs,
  1046. aes(xmin = onset, xmax = end, ymin = -1.1, ymax = 1.1, fill = trial_type),
  1047. alpha = 0.3, color = NA) +
  1048. scale_fill_manual(values = c("faces" = "#F9B6C7", "scenes" = "#3B7D23")) +
  1049. guides(fill = guide_legend(title = "Event type"), title.theme = element_text(face = "bold"))
  1050. print(Plot_amygdala)
  1051. ## STS Timecourses, including events
  1052. Plot_STS <- ggplot(data_plot_STS) +
  1053. geom_line(aes(x = time, y = value, color = AgeGroup, group = AgeGroup), size = 1.2, alpha = 1) +
  1054. geom_hline(yintercept = 0, color = "grey", linewidth = 0.8) +
  1055. geom_text(data = data_plot_STS[data_plot_STS$time %in% c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170),],
  1056. aes(x = time, y = -0.07, label = time), vjust = 0, size = 3, color = "black", angle = 90) +
  1057. facet_wrap(~ hemisphere, scales = "free_y") +
  1058. theme_classic() +
  1059. scale_color_manual(values = custom_colours) +
  1060. scale_y_continuous(limits = c(-1.1, 1.1)) +
  1061. scale_x_continuous(
  1062. breaks = c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170),
  1063. labels = NULL
  1064. ) +
  1065. labs(
  1066. title = "STS",
  1067. y = "Response magnitude",
  1068. color = "Age group"
  1069. ) +
  1070. theme(
  1071. strip.text = element_text(size = 16),
  1072. strip.background = element_blank(),
  1073. legend.position = "bottom",
  1074. axis.title.x = element_text(size = 16),
  1075. axis.title.y = element_text(size = 16),
  1076. axis.text.x = element_text(size = 16),
  1077. axis.text.y = element_text(size = 16),
  1078. axis.ticks.x = element_blank(),
  1079. axis.line.x = element_blank(),
  1080. plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
  1081. axis.line.y = element_line(linewidth = 0.8, color = "black"),
  1082. legend.text = element_text(size = 16),
  1083. legend.title = element_text(size = 16, face = "bold")
  1084. )
  1085. ### Add the shaded rectangles for events
  1086. Plot_STS <- Plot_STS +
  1087. geom_rect(data = events_TRs,
  1088. aes(xmin = onset, xmax = end, ymin = -1.1, ymax = 1.1, fill = trial_type),
  1089. alpha = 0.3, color = NA) +
  1090. scale_fill_manual(values = c("faces" = "#F9B6C7", "scenes" = "#3B7D23")) +
  1091. guides(fill = guide_legend(title = "Event type"), title.theme = element_text(face = "bold"))
  1092. ### Get the 10th "faces" event (pink box)
  1093. pink_boxes <- events_TRs[events_TRs$trial_type == "faces",]
  1094. fifth_pink_box <- pink_boxes[5,] #get the 10th pink box
  1095. ### Adjust the y position for the asterisk to stay within the valid range
  1096. Plot_STS <- Plot_STS +
  1097. geom_text(data = fifth_pink_box,
  1098. aes(x = (onset + end) / 2 - 1, y = 1.05, label = "*"), # Adjust y to fit within plot limits
  1099. size = 8, color = "black", fontface = "bold", vjust = 0)
  1100. print(Plot_STS)
  1101. ## Save plots
  1102. ggsave(here("results/figures", "TCs_FFA.png"), plot = Plot_FFA, width = 15, height = 5, units = "in", dpi = 300)
  1103. ggsave(here("results/figures", "TCs_MMPFC.png"), plot = Plot_MMPFC, width = 15, height = 5, units = "in", dpi = 300)
  1104. ggsave(here("results/figures", "TCs_amygdala.png"), plot = Plot_amygdala, width = 15, height = 5, units = "in", dpi = 300)
  1105. ggsave(here("results/figures", "TCs_STS.png"), plot = Plot_STS, width = 15, height = 5, units = "in", dpi = 300)
  1106. ```
  1107. Age and response magnitude of MMPFC and STS to face events
  1108. Histogram
  1109. ```{r}
  1110. # Data wrangling
  1111. data_plot <- merge(long_data, participants[,c(1,3)], by="ID")
  1112. data_plot$ID <- as.factor(data_plot$ID)
  1113. data_plot$hemisphere <- as.factor(data_plot$hemisphere)
  1114. data_plot$region <- as.factor(data_plot$region)
  1115. data_plot_MMPFC <- subset(data_plot, region == "MMPFC")
  1116. data_plot_STS <- subset(data_plot, region == "STS")
  1117. # Colours
  1118. custom_colours <- c("#D34835", "#DDA86A", "#83AC66", "#87BACF", "#6989F3", "#652975")
  1119. # Make individual plots for MMPFC and STS and events with significant age associations
  1120. ## General plotting function
  1121. generate_plot <- function(data, feature, subtitle, show_legend = TRUE, y_label = NULL, legend_pos = "bottom") {
  1122. summary_data <- data %>%
  1123. filter(!is.na(.data[[feature]]) & !is.na(AgeGroup)) %>%
  1124. group_by(hemisphere, AgeGroup) %>%
  1125. summarise(
  1126. mean_value = mean(.data[[feature]], na.rm = TRUE),
  1127. se_value = sd(.data[[feature]], na.rm = TRUE) / sqrt(n()),
  1128. .groups = "drop"
  1129. )
  1130. ggplot(summary_data, aes(x = as.factor(AgeGroup), y = mean_value, fill = as.factor(AgeGroup))) +
  1131. geom_bar(stat = "identity", position = "stack", width = 1, alpha = 1) +
  1132. geom_errorbar(aes(ymin = mean_value - se_value, ymax = mean_value + se_value), width = 0) +
  1133. geom_hline(yintercept = 0, linetype = "solid", size = 1, color = "black") +
  1134. scale_fill_manual(values = custom_colours) +
  1135. labs(title = subtitle, x = "Age Group", y = y_label, fill = "Age Group") +
  1136. facet_wrap(~hemisphere) +
  1137. coord_cartesian(ylim = c(-1, 1)) +
  1138. theme_classic() +
  1139. theme(
  1140. strip.text = element_text(size = 14, face = "bold"),
  1141. strip.background= element_blank(),
  1142. legend.position = ifelse(show_legend, legend_pos, "right"),
  1143. legend.text = element_text(size = 16),
  1144. legend.title = element_text(size = 14, face = "bold"),
  1145. axis.title = element_text(size = 14, face = "bold"),
  1146. axis.text = element_text(size = 12),
  1147. axis.title.x = element_blank(),
  1148. axis.text.x = element_blank(),
  1149. axis.ticks.x = element_blank(),
  1150. axis.line.x = element_blank(),
  1151. plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
  1152. axis.line.y = element_line(size = 0.8, color = "black")
  1153. )
  1154. }
  1155. ## Generate plots with one reusable function
  1156. plot_F05 <- generate_plot(data_plot_MMPFC, "F05", "MMPFC response to F05", show_legend = TRUE, y_label = "Response magnitude")
  1157. plot_F10 <- generate_plot(data_plot_MMPFC, "F10", "STS response to F10", show_legend = TRUE)
  1158. ## Combine plots with a common centered title and shared legend
  1159. final_plot <- (plot_F05 | plot_F10 ) +
  1160. plot_layout(ncol = 2, guides = "collect") +
  1161. plot_annotation(
  1162. title = ""
  1163. ) &
  1164. theme(
  1165. legend.position = "bottom",
  1166. legend.text = element_text(size = 12),
  1167. legend.title = element_text(size = 12),
  1168. plot.title = element_text(size = 14, hjust = 0.5, face = "bold")
  1169. )
  1170. # Save plot
  1171. ggsave(here("results/figures", "Event_age_histogram.png"), plot = final_plot, width = 7, height = 5, units = "in", dpi = 300)
  1172. ```
  1173. ## Section III: Age effects on functional connectivity (Fx) between regions
  1174. ### Read data
  1175. ```{r}
  1176. # Covariates
  1177. participants <- read_tsv(here("raw_data/covariates", "participants.tsv"), col_names=TRUE)
  1178. colnames(participants)[which(names(participants) == "participant_id")] <- "ID"
  1179. ## Remove the word "pixar" and trim extra whitespace
  1180. participants$ID <- gsub("[sub-pixar-]","", participants$ID)
  1181. participants$ID <- trimws(participants$ID)
  1182. # Motion
  1183. motion <- read_tsv(here("raw_data/motion", "outlier_info.tsv"), col_names=TRUE)
  1184. colnames(motion)[which(names(motion) == "subject")] <- "ID"
  1185. ## Remove the word "pixar" and trim extra whitespace
  1186. motion$ID <- gsub("[sub-pixar-]","", motion$ID)
  1187. motion$ID <- trimws(motion$ID)
  1188. # Child fMRI variables
  1189. child_MRIvariables <- read_csv(here("processed_data", "child_MRIvariables.csv"), col_names=TRUE)
  1190. # Complete dataset
  1191. child_MRIvariables <- child_MRIvariables[,c(1, 202:229)]
  1192. child_MRIvariables <- merge(child_MRIvariables, participants[,c(1:2)], by="ID")
  1193. child_MRIvariables <- merge(child_MRIvariables, motion[,c(4,6)], by="ID")
  1194. child_MRIvariables <- child_MRIvariables %>%
  1195. rename_with(~ sub("^zscored_", "", .x), starts_with("zscored_"))
  1196. child_MRIvariables[, c("Age", "MeanFD")] <- scale(child_MRIvariables[, c("Age", "MeanFD")])
  1197. ```
  1198. ### Plot correlations between regions by age
  1199. ```{r}
  1200. # Read data
  1201. TCs_total <- read_csv(here("processed_data", "TCs_total.csv"), col_names = TRUE)
  1202. TCs_total$time <- as.numeric(as.character(TCs_total$time))
  1203. data_plot <- merge(TCs_total, participants[, c("ID", "AgeGroup")], by = "ID")
  1204. data_plot <- data_plot %>%
  1205. filter(AgeGroup != "Adult")
  1206. # Rename timecourse columns
  1207. data_plot <- data_plot %>%
  1208. rename_with(~ c("left amygdala", "right amygdala", "left FFA", "right FFA", "left MMPFC", "right MMPFC", "left STS", "right STS"),
  1209. .cols = c("lamygdala", "ramygdala", "lFFA", "rFFA", "lMMPFC", "rMMPFC", "lSTS", "rSTS"))
  1210. # Define regions of interest
  1211. regions <- c("left amygdala", "right amygdala", "left FFA", "right FFA", "left MMPFC", "right MMPFC", "left STS", "right STS")
  1212. # Calculate correlation matrices per AgeGroup
  1213. age_groups <- unique(data_plot$AgeGroup)
  1214. correlation_by_agegroup <- map_df(age_groups, function(age) {
  1215. group_data <- data_plot[data_plot$AgeGroup == age, ]
  1216. seg1 <- group_data[group_data$time >= 1 & group_data$time <= 82, ]
  1217. seg2 <- group_data[group_data$time >= 87 & group_data$time <= 168, ]
  1218. mat1 <- cor(seg1[, regions], use = "pairwise.complete.obs")
  1219. mat2 <- cor(seg2[, regions], use = "pairwise.complete.obs")
  1220. mat_avg <- (mat1 + mat2) / 2
  1221. # Replace amygdala correlation with full-time tc version
  1222. amyg_cor <- cor(group_data[, c("left amygdala", "right amygdala")], use = "pairwise.complete.obs")[1, 2]
  1223. mat_avg["left amygdala", "right amygdala"] <- amyg_cor
  1224. mat_avg["right amygdala", "left amygdala"] <- amyg_cor
  1225. tibble(AgeGroup = age, cor_matrix = list(mat_avg))
  1226. })
  1227. # Convert matrices to long format for ggplot2
  1228. data_plot <- correlation_by_agegroup %>%
  1229. mutate(cor_matrix = map(cor_matrix, ~ melt(.x, varnames = c("Region1", "Region2")))) %>%
  1230. unnest(cor_matrix)
  1231. data_plot <- data_plot %>%
  1232. filter(Region1 != Region2)
  1233. # Order factor levels for consistent axis order
  1234. region_levels <- regions
  1235. data_plot$Region1 <- factor(data_plot$Region1, levels = region_levels)
  1236. data_plot$Region2 <- factor(data_plot$Region2, levels = region_levels)
  1237. # Identify upper triangle entries (for displaying numbers only)
  1238. data_plot <- data_plot %>%
  1239. mutate(is_upper = as.numeric(Region1) < as.numeric(Region2))
  1240. # Plot heatmap
  1241. heatmap_plot <- ggplot(data_plot, aes(x = Region1, y = Region2, fill = value)) +
  1242. geom_tile() +
  1243. geom_tile(data = data_plot %>% filter(is_upper), fill = "white", color = NA) +
  1244. geom_text(data = data_plot %>% filter(is_upper),
  1245. aes(label = sprintf("%.2f", value)), size = 2.8) +
  1246. facet_wrap(~ AgeGroup, nrow=1) +
  1247. scale_fill_gradient2(
  1248. low = "#D34835", high = "#87BACF", mid = "white",
  1249. midpoint = 0, limit = c(-0.80, 0.80),
  1250. name = "Correlation",
  1251. breaks = c(-0.80, -0.40, 0, 0.40, 0.80),
  1252. labels = c("-0.80", "-0.40", "0.00", "0.40", "0.80")
  1253. ) +
  1254. theme_minimal() +
  1255. labs(
  1256. title = "",
  1257. x = "",
  1258. y = ""
  1259. ) +
  1260. theme(
  1261. legend.position = "right",
  1262. axis.text.x = element_text(angle = 45, hjust = 1, face = "bold", color = "darkgrey", size = 13),
  1263. axis.text.y = element_text(face = "bold", color = "darkgrey", size = 13),
  1264. strip.text = element_text(size = 14),
  1265. plot.title = element_text(hjust = 0.5, size = 19),
  1266. plot.title.position = "plot"
  1267. )
  1268. # Display the plot
  1269. print(heatmap_plot)
  1270. # Save plot
  1271. ggsave(here("results/figures", "Fx_AgeGroup.png"), plot = heatmap_plot, width = 13, height = 4, units = "in", dpi = 300)
  1272. ```
  1273. ### Descriptive statistics: Mean FFA Fx values with ROIs (non-preregistered)
  1274. ```{r}
  1275. # Read data
  1276. TCs_total <- read_csv(here("processed_data", "TCs_total.csv"), col_names = TRUE)
  1277. TCs_total$time <- as.numeric(as.character(TCs_total$time))
  1278. data_plot <- merge(TCs_total, participants[, c("ID", "AgeGroup")], by = "ID")
  1279. data_plot <- data_plot %>%
  1280. filter(AgeGroup != "Adult")
  1281. # Rename ROI columns
  1282. data_plot <- data_plot %>%
  1283. rename_with(
  1284. ~ c("left amygdala", "right amygdala",
  1285. "left FFA", "right FFA",
  1286. "left MMPFC", "right MMPFC",
  1287. "left STS", "right STS"),
  1288. .cols = c("lamygdala", "ramygdala",
  1289. "lFFA", "rFFA",
  1290. "lMMPFC", "rMMPFC",
  1291. "lSTS", "rSTS")
  1292. )
  1293. # Define ROIs
  1294. regions <- c("left amygdala", "right amygdala",
  1295. "left FFA", "right FFA",
  1296. "left MMPFC", "right MMPFC",
  1297. "left STS", "right STS")
  1298. # Compute subject-level correlation matrices
  1299. cor_subject <- data_plot %>%
  1300. group_by(ID) %>%
  1301. group_map(~ {
  1302. seg1 <- .x %>% filter(time >= 1 & time <= 82)
  1303. seg2 <- .x %>% filter(time >= 87 & time <= 168)
  1304. mat1 <- cor(seg1[, regions], use = "pairwise.complete.obs")
  1305. mat2 <- cor(seg2[, regions], use = "pairwise.complete.obs")
  1306. mat_avg <- (mat1 + mat2) / 2
  1307. # Replace amygdala correlation with full-time tc version
  1308. amyg_cor <- cor(.x[, c("left amygdala", "right amygdala")], use = "pairwise.complete.obs")[1,2]
  1309. mat_avg["left amygdala", "right amygdala"] <- amyg_cor
  1310. mat_avg["right amygdala", "left amygdala"] <- amyg_cor
  1311. mat_avg
  1312. })
  1313. # Convert all matrices to long format
  1314. cor_long <- map_dfr(cor_subject, ~ as.data.frame(as.table(.x)) %>%
  1315. setNames(c("Region1", "Region2", "r")))
  1316. # Exclude diagonal (self-correlations)
  1317. cor_long <- cor_long %>%
  1318. filter(Region1 != Region2)
  1319. # Filter for left/right FFA and compute summary + t-tests
  1320. cor_summary_with_t <- cor_long %>%
  1321. filter(Region1 %in% c("right FFA", "left FFA")) %>%
  1322. group_by(Region1, Region2) %>%
  1323. summarise(
  1324. mean_r = mean(r, na.rm = TRUE),
  1325. sd_r = sd(r, na.rm = TRUE),
  1326. n = n(),
  1327. t_test = list(t.test(r, mu = 0, na.rm = TRUE)),
  1328. .groups = "drop"
  1329. ) %>%
  1330. mutate(
  1331. mean_sd = sprintf("%.2f (%.2f)", mean_r, sd_r),
  1332. t_test = map(t_test, broom::tidy)
  1333. ) %>%
  1334. unnest(t_test) %>%
  1335. mutate(
  1336. p.value = ifelse(p.value < 0.001,
  1337. format(p.value, scientific = TRUE, digits = 4),
  1338. sprintf("%.3f", p.value)),
  1339. t_stat = sprintf("%.2f", statistic),
  1340. df = sprintf("%.0f", parameter)
  1341. ) %>%
  1342. select(Region1, Region2, mean_sd, t_stat, df, p.value) %>%
  1343. arrange(Region1, Region2)
  1344. # View final table
  1345. cor_summary_with_t
  1346. # Save the table
  1347. write_csv(cor_summary_with_t, here("results/RQ1", "FFA_fx_descriptive_table.csv"))
  1348. ```
  1349. ### Models (models including STS were not pre-registered)
  1350. #### Primary models
  1351. Linear models
  1352. ```{r}
  1353. model_rlFFA <- lm(fx_rFFA_lFFA ~ Age + MeanFD, data = child_MRIvariables)
  1354. summary(model_rlFFA)
  1355. ```
  1356. Linear-mixed-effect models
  1357. ```{r}
  1358. # Make long data for rest of the models:
  1359. long_data <- child_MRIvariables %>%
  1360. pivot_longer(
  1361. cols = starts_with("fx_"),
  1362. names_to = c("Region1", "Region2"),
  1363. names_pattern = "fx_(.+)_(.+)",
  1364. values_to = "zscored_value"
  1365. )
  1366. long_data <- long_data %>%
  1367. mutate(
  1368. Hemisphere = ifelse(grepl("^r", Region2), "right", "left"),
  1369. Region2 = sub("^[rl]", "", Region2)
  1370. )
  1371. # Models with interaction effects
  1372. ## rFFA and r/lMMPFC
  1373. model_data <- long_data %>%
  1374. filter(Region1 == "rFFA" & Region2 == "MMPFC")
  1375. model_rFFA_rlMMPFC <- lmer(
  1376. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1377. data = model_data
  1378. )
  1379. summary(model_rFFA_rlMMPFC)
  1380. ## rFFA and r/lamygdala
  1381. model_data <- long_data %>%
  1382. filter(Region1 == "rFFA" & Region2 == "amygdala")
  1383. model_rFFA_rlamygdala <- lmer(
  1384. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1385. data = model_data
  1386. )
  1387. summary(model_rFFA_rlamygdala)
  1388. ## rFFA and r/lSTS
  1389. model_data <- long_data %>%
  1390. filter(Region1 == "rFFA" & Region2 == "STS")
  1391. model_rFFA_rlSTS <- lmer(
  1392. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1393. data = model_data
  1394. )
  1395. summary(model_rFFA_rlSTS)
  1396. # No interaction effects are significant except for rFFA and r/lMMPFC, so models are run without interaction effects except for rFFA and r/lMMPFC
  1397. ## rFFA and r/lMMPFC
  1398. model_data <- long_data %>%
  1399. filter(Region1 == "rFFA" & Region2 == "MMPFC")
  1400. model_rFFA_rlMMPFC <- lmer(
  1401. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1402. data = model_data
  1403. )
  1404. summary(model_rFFA_rlMMPFC)
  1405. ## rFFA and r/lamygdala
  1406. model_data <- long_data %>%
  1407. filter(Region1 == "rFFA" & Region2 == "amygdala")
  1408. model_rFFA_rlamygdala <- lmer(
  1409. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1410. data = model_data
  1411. )
  1412. summary(model_rFFA_rlamygdala)
  1413. ## rFFA and r/lSTS
  1414. model_data <- long_data %>%
  1415. filter(Region1 == "rFFA" & Region2 == "STS")
  1416. model_rFFA_rlSTS <- lmer(
  1417. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1418. data = model_data
  1419. )
  1420. summary(model_rFFA_rlSTS)
  1421. ```
  1422. #### Secondary models
  1423. Linear models
  1424. ```{r}
  1425. # rlMMPFC
  1426. model_rlMMPFC <- lm(fx_lMMPFC_rMMPFC ~ Age + MeanFD, data = child_MRIvariables)
  1427. summary(model_rlMMPFC)
  1428. # rlamygdala
  1429. model_rlamygdala <- lm(fx_lamygdala_ramygdala ~ Age + MeanFD, data = child_MRIvariables)
  1430. summary(model_rlamygdala)
  1431. # rlSTS
  1432. model_rlSTS <- lm(fx_lSTS_rSTS ~ Age + MeanFD, data = child_MRIvariables)
  1433. summary(model_rlSTS)
  1434. ```
  1435. Linear-mixed-effect models
  1436. ```{r}
  1437. # Models with interaction effects
  1438. ## lFFA and r/lMMPFC
  1439. model_data <- long_data %>%
  1440. filter(Region1 == "lFFA" & Region2 == "MMPFC")
  1441. model_lFFA_rlMMPFC <- lmer(
  1442. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1443. data = model_data
  1444. )
  1445. summary(model_lFFA_rlMMPFC)
  1446. ## lFFA and r/lamygdala
  1447. model_data <- long_data %>%
  1448. filter(Region1 == "lFFA" & Region2 == "amygdala")
  1449. model_lFFA_rlamygdala <- lmer(
  1450. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1451. data = model_data
  1452. )
  1453. summary(model_lFFA_rlamygdala)
  1454. ## lFFA and r/lSTS
  1455. model_data <- long_data %>%
  1456. filter(Region1 == "lFFA" & Region2 == "STS")
  1457. model_lFFA_rlSTS <- lmer(
  1458. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1459. data = model_data
  1460. )
  1461. summary(model_lFFA_rlSTS)
  1462. ## rMMPFC and r/lamygdala
  1463. model_data <- long_data %>%
  1464. filter(Region1 == "rMMPFC" & Region2 == "amygdala")
  1465. model_rMMPFC_rlamygdala <- lmer(
  1466. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1467. data = model_data
  1468. )
  1469. summary(model_rMMPFC_rlamygdala)
  1470. ## rMMPFC and r/lSTS
  1471. model_data <- long_data %>%
  1472. filter(Region1 == "rMMPFC" & Region2 == "STS")
  1473. model_rMMPFC_rlSTS <- lmer(
  1474. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1475. data = model_data
  1476. )
  1477. summary(model_rMMPFC_rlSTS)
  1478. ## lMMPFC and r/lamygdala
  1479. model_data <- long_data %>%
  1480. filter(Region1 == "lMMPFC" & Region2 == "amygdala")
  1481. model_lMMPFC_rlamygdala <- lmer(
  1482. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1483. data = model_data
  1484. )
  1485. summary(model_lMMPFC_rlamygdala)
  1486. ## lMMPFC and r/lSTS
  1487. model_data <- long_data %>%
  1488. filter(Region1 == "lMMPFC" & Region2 == "STS")
  1489. model_lMMPFC_rlSTS <- lmer(
  1490. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1491. data = model_data
  1492. )
  1493. summary(model_lMMPFC_rlSTS)
  1494. ## ramygdala and r/lSTS
  1495. model_data <- long_data %>%
  1496. filter(Region1 == "ramygdala" & Region2 == "STS")
  1497. model_ramygdala_rlSTS <- lmer(
  1498. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1499. data = model_data
  1500. )
  1501. summary(model_ramygdala_rlSTS)
  1502. ## lamygdala and r/lSTS
  1503. model_data <- long_data %>%
  1504. filter(Region1 == "lamygdala" & Region2 == "STS")
  1505. model_lamygdala_rlSTS <- lmer(
  1506. zscored_value ~ Age*Hemisphere + MeanFD + (1 | ID),
  1507. data = model_data
  1508. )
  1509. summary(model_lamygdala_rlSTS)
  1510. # No interaction effects are significant, so models are run without interaction effects
  1511. ## lFFA and r/lMMPFC
  1512. model_data <- long_data %>%
  1513. filter(Region1 == "lFFA" & Region2 == "MMPFC")
  1514. model_lFFA_rlMMPFC <- lmer(
  1515. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1516. data = model_data
  1517. )
  1518. summary(model_lFFA_rlMMPFC)
  1519. ## lFFA and r/lamygdala
  1520. model_data <- long_data %>%
  1521. filter(Region1 == "lFFA" & Region2 == "amygdala")
  1522. model_lFFA_rlamygdala <- lmer(
  1523. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1524. data = model_data
  1525. )
  1526. summary(model_lFFA_rlamygdala)
  1527. ## lFFA and r/lSTS
  1528. model_data <- long_data %>%
  1529. filter(Region1 == "lFFA" & Region2 == "STS")
  1530. model_lFFA_rlSTS <- lmer(
  1531. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1532. data = model_data
  1533. )
  1534. summary(model_lFFA_rlSTS)
  1535. ## rMMPFC and r/lamygdala
  1536. model_data <- long_data %>%
  1537. filter(Region1 == "rMMPFC" & Region2 == "amygdala")
  1538. model_rMMPFC_rlamygdala <- lmer(
  1539. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1540. data = model_data
  1541. )
  1542. summary(model_rMMPFC_rlamygdala)
  1543. ## rMMPFC and r/lSTS
  1544. model_data <- long_data %>%
  1545. filter(Region1 == "rMMPFC" & Region2 == "STS")
  1546. model_rMMPFC_rlSTS <- lmer(
  1547. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1548. data = model_data
  1549. )
  1550. summary(model_rMMPFC_rlSTS)
  1551. ## lMMPFC and r/lamygdala
  1552. model_data <- long_data %>%
  1553. filter(Region1 == "lMMPFC" & Region2 == "amygdala")
  1554. model_lMMPFC_rlamygdala <- lmer(
  1555. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1556. data = model_data
  1557. )
  1558. summary(model_lMMPFC_rlamygdala)
  1559. ## lMMPFC and r/lSTS
  1560. model_data <- long_data %>%
  1561. filter(Region1 == "lMMPFC" & Region2 == "STS")
  1562. model_lMMPFC_rlSTS <- lmer(
  1563. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1564. data = model_data
  1565. )
  1566. summary(model_lMMPFC_rlSTS)
  1567. ## ramygdala and r/lSTS
  1568. model_data <- long_data %>%
  1569. filter(Region1 == "ramygdala" & Region2 == "STS")
  1570. model_ramygdala_rlSTS <- lmer(
  1571. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1572. data = model_data
  1573. )
  1574. summary(model_ramygdala_rlSTS)
  1575. ## lamygdala and r/lSTS
  1576. model_data <- long_data %>%
  1577. filter(Region1 == "lamygdala" & Region2 == "STS")
  1578. model_lamygdala_rlSTS <- lmer(
  1579. zscored_value ~ Age + Hemisphere + MeanFD + (1 | ID),
  1580. data = model_data
  1581. )
  1582. summary(model_lamygdala_rlSTS)
  1583. ```
  1584. Save models
  1585. ```{r}
  1586. # Create tidy summaries of all models
  1587. ## Primary models
  1588. model_rlFFA_summary <- tidy(model_rlFFA) %>%
  1589. mutate(Model = "rlFFA", Label = "Primary")
  1590. model_rFFA_rlMMPFC_summary <- tidy(model_rFFA_rlMMPFC, effects = "fixed") %>%
  1591. mutate(Model = "rFFA_rlMMPFC", Label = "Primary")
  1592. model_rFFA_rlamygdala_summary <- tidy(model_rFFA_rlamygdala, effects = "fixed") %>%
  1593. mutate(Model = "rFFA_rlamygdala", Label = "Primary")
  1594. model_rFFA_rlSTS_summary <- tidy(model_rFFA_rlSTS, effects = "fixed") %>%
  1595. mutate(Model = "rFFA_rlSTS", Label = "Primary")
  1596. ## Secondary models
  1597. model_lFFA_rlMMPFC_summary <- tidy(model_lFFA_rlMMPFC, effects = "fixed") %>%
  1598. mutate(Model = "lFFA_rlMMPFC", Label = "Secondary")
  1599. model_lFFA_rlamygdala_summary <- tidy(model_lFFA_rlamygdala, effects = "fixed") %>%
  1600. mutate(Model = "lFFA_rlamygdala", Label = "Secondary")
  1601. model_lFFA_rlSTS_summary <- tidy(model_lFFA_rlSTS, effects = "fixed") %>%
  1602. mutate(Model = "lFFA_rlSTS", Label = "Secondary")
  1603. model_rlMMPFC_summary <- tidy(model_rlMMPFC) %>%
  1604. mutate(Model = "rlMMPFC", Label = "Secondary")
  1605. model_rMMPFC_rlamygdala_summary <- tidy(model_rMMPFC_rlamygdala, effects = "fixed") %>%
  1606. mutate(Model = "rMMPFC_rlamygdala", Label = "Secondary")
  1607. model_rMMPFC_rlSTS_summary <- tidy(model_rMMPFC_rlSTS, effects = "fixed") %>%
  1608. mutate(Model = "rMMPFC_rlSTS", Label = "Secondary")
  1609. model_lMMPFC_rlamygdala_summary <- tidy(model_lMMPFC_rlamygdala, effects = "fixed") %>%
  1610. mutate(Model = "lMMPFC_rlamygdala", Label = "Secondary")
  1611. model_lMMPFC_rlSTS_summary <- tidy(model_lMMPFC_rlSTS, effects = "fixed") %>%
  1612. mutate(Model = "lMMPFC_rlSTS", Label = "Secondary")
  1613. model_rlamygdala_summary <- tidy(model_rlamygdala, effects = "fixed") %>%
  1614. mutate(Model = "rlamygdala", Label = "Secondary")
  1615. model_ramygdala_rlSTS_summary <- tidy(model_ramygdala_rlSTS, effects = "fixed") %>%
  1616. mutate(Model = "ramygdala_rlSTS", Label = "Secondary")
  1617. model_lamygdala_rlSTS_summary <- tidy(model_lamygdala_rlSTS, effects = "fixed") %>%
  1618. mutate(Model = "lamygdala_rlSTS", Label = "Secondary")
  1619. model_rlSTS_summary <- tidy(model_rlSTS, effects = "fixed") %>%
  1620. mutate(Model = "rlSTS", Label = "Secondary")
  1621. # Combine all summaries into one data frame
  1622. all_model_summaries <- bind_rows(
  1623. model_rlFFA_summary,
  1624. model_rFFA_rlMMPFC_summary,
  1625. model_rFFA_rlamygdala_summary,
  1626. model_rFFA_rlSTS_summary,
  1627. model_lFFA_rlMMPFC_summary,
  1628. model_lFFA_rlamygdala_summary,
  1629. model_lFFA_rlSTS_summary,
  1630. model_rlMMPFC_summary,
  1631. model_lMMPFC_rlamygdala_summary,
  1632. model_lMMPFC_rlSTS_summary,
  1633. model_rMMPFC_rlamygdala_summary,
  1634. model_rMMPFC_rlSTS_summary,
  1635. model_rlamygdala_summary,
  1636. model_lamygdala_rlSTS_summary,
  1637. model_ramygdala_rlSTS_summary,
  1638. model_rlSTS_summary
  1639. ) %>%
  1640. mutate(
  1641. estimate = round(estimate, 2),
  1642. std.error = round(std.error, 2),
  1643. statistic = round(statistic, 2),
  1644. df = round(df, 0),
  1645. p.value = if_else(
  1646. p.value < 0.001,
  1647. format(p.value, scientific = TRUE, digits = 4),
  1648. sprintf("%.3f", p.value)
  1649. ))
  1650. # Save the combined summaries to a CSV file
  1651. write.csv(all_model_summaries, here("results/RQ1", "Fx_age.csv"), row.names = FALSE)
  1652. ```
  1653. ### Model plots
  1654. ```{r}
  1655. # Data wrangling
  1656. data_plot <- read_csv(here("processed_data", "child_MRIvariables.csv"), col_names=TRUE)
  1657. data_plot <- data_plot[,c(1, 202:229)]
  1658. data_plot <- merge(data_plot, participants[,c(1:2)], by="ID")
  1659. data_plot <- merge(data_plot, motion[,c(4,6)], by="ID")
  1660. data_plot <- data_plot %>%
  1661. rename_with(~ sub("^zscored_", "", .x), starts_with("zscored_"))
  1662. data_plot <- data_plot %>%
  1663. pivot_longer(
  1664. cols = starts_with("fx_"),
  1665. names_to = c("Region1", "Region2"),
  1666. names_pattern = "fx_(.+)_(.+)",
  1667. values_to = "zscored_value"
  1668. )
  1669. data_plot <- data_plot %>%
  1670. mutate(
  1671. Hemisphere = ifelse(grepl("^r", Region2), "right", "left"),
  1672. Region2 = sub("^[rl]", "", Region2)
  1673. )
  1674. data_plot_rFFA_MMPFC <- data_plot %>%
  1675. filter(Region2 == "MMPFC", Region1 == "rFFA")
  1676. data_plot_lFFA_MMPFC <- data_plot %>%
  1677. filter(Region2 == "MMPFC", Region1 == "lFFA")
  1678. data_plot_rFFA_amygdala <- data_plot %>%
  1679. filter(Region2 == "amygdala", Region1 == "rFFA")
  1680. data_plot_lFFA_amygdala <- data_plot %>%
  1681. filter(Region2 == "amygdala", Region1 == "lFFA")
  1682. data_plot_rFFA_STS <- data_plot %>%
  1683. filter(Region2 == "STS", Region1 == "rFFA")
  1684. data_plot_lFFA_STS <- data_plot %>%
  1685. filter(Region2 == "STS", Region1 == "lFFA")
  1686. # Colours
  1687. custom_colours_amygdala <- c("#E1AA9F", "#D34835")
  1688. custom_colours_MMPFC <- c("#87BACF", "#578CAD")
  1689. custom_colours_FFA <- c("#C39EDA", "#652975")
  1690. custom_colours_STS <- c("#FFF865", "#FFCE1B")
  1691. # Define a base theme to avoid repetition
  1692. base_theme <- theme_minimal() +
  1693. theme(
  1694. strip.text = element_text(size = 14, face = "bold"),
  1695. strip.background = element_blank(),
  1696. axis.line.x = element_blank(),
  1697. axis.text = element_text(size = 12),
  1698. axis.title = element_text(size = 14, face = "bold"),
  1699. plot.title = element_text(size = 14, hjust = 0.5, face = "bold"),
  1700. legend.text = element_text(size = 12),
  1701. legend.title = element_text(size = 12, face = "bold")
  1702. )
  1703. # Function to create ggplots efficiently
  1704. create_plot <- function(data, color_values, legend_title, x_label = "", y_label = "", legend_position = "none", remove_y_axis_labels = FALSE) {
  1705. plot <- ggplot(data, aes(x = Age, y = zscored_value, color = Hemisphere)) +
  1706. geom_point(alpha = 0.8, size = 2.8, shape = 19) +
  1707. geom_smooth(method = "lm", se = FALSE, linetype = "solid", linewidth = 2) +
  1708. scale_color_manual(values = color_values) +
  1709. labs(x = x_label, y = y_label, color = legend_title) +
  1710. base_theme +
  1711. theme(legend.position = legend_position) +
  1712. coord_cartesian(ylim = c(-4, 3)) +
  1713. scale_x_continuous(breaks = c(3, 6, 9, 12))
  1714. if (remove_y_axis_labels) {
  1715. plot <- plot + theme(
  1716. axis.title.y = element_blank(),
  1717. axis.text.y = element_blank()
  1718. )
  1719. }
  1720. return(plot)
  1721. }
  1722. # Create all plots using the function
  1723. plot_lFFA_MMPFC <- create_plot(data_plot_lFFA_MMPFC, custom_colours_MMPFC, "MMPFC Hemisphere", "", "Functional connectivity between\n left FFA and region")
  1724. plot_lFFA_amygdala <- create_plot(data_plot_lFFA_amygdala, custom_colours_amygdala, "Amygdala Hemisphere", "Age", legend_position = "bottom", remove_y_axis_labels = TRUE)
  1725. plot_lFFA_STS <- create_plot(data_plot_lFFA_STS, custom_colours_STS, "STS Hemisphere", "", legend_position = "bottom", remove_y_axis_labels = TRUE)
  1726. plot_rFFA_MMPFC <- create_plot(data_plot_rFFA_MMPFC, custom_colours_MMPFC, "MMPFC Hemisphere", "", "Functional connectivity between\n right FFA and region", legend_position = "bottom") +
  1727. annotate(
  1728. "text",
  1729. x = Inf, y = -Inf, # bottom right corner
  1730. label = "Age x Hemisphere (right) effect: \nβ(SE)=0.16(0.07), p=0.016", # your annotation
  1731. hjust = 1.05, # slightly >1 for proper right alignment
  1732. vjust = -0.5, # nudges text up from bottom
  1733. size = 4,
  1734. color = "black"
  1735. )
  1736. plot_rFFA_amygdala <- create_plot(data_plot_rFFA_amygdala, custom_colours_amygdala, "Amygdala Hemisphere", "Age", remove_y_axis_labels = TRUE)
  1737. plot_rFFA_STS <- create_plot(data_plot_rFFA_STS, custom_colours_STS, "STS Hemisphere", "", remove_y_axis_labels = TRUE)
  1738. # Combine plots into a single row
  1739. fx_rFFA_age <- (
  1740. plot_rFFA_MMPFC | plot_rFFA_amygdala | plot_rFFA_STS
  1741. ) +
  1742. plot_layout(
  1743. ncol = 3,
  1744. widths = rep(1, 3),
  1745. guides = "collect"
  1746. ) +
  1747. plot_annotation(
  1748. title = "Functional Connectivity Across Age",
  1749. theme = theme(
  1750. plot.title = element_text(size = 14, hjust = 0.5, face = "bold")
  1751. )
  1752. ) &
  1753. theme(legend.position = "bottom")
  1754. fx_rFFA_age
  1755. fx_lFFA_age <- (
  1756. plot_lFFA_MMPFC | plot_lFFA_amygdala | plot_lFFA_STS ) +
  1757. plot_layout(
  1758. ncol = 3,
  1759. widths = rep(1, 3),
  1760. guides = "collect"
  1761. ) +
  1762. plot_annotation(
  1763. title = "Functional Connectivity Across Age",
  1764. theme = theme(
  1765. plot.title = element_text(size = 14, hjust = 0.5, face = "bold")
  1766. )
  1767. ) &
  1768. theme(legend.position = "bottom")
  1769. fx_lFFA_age
  1770. # Save plots
  1771. ggsave(here("results/figures", "Fx_rFFA_age.png"), plot = fx_rFFA_age, width = 10, height = 5, units = "in", dpi = 300)
  1772. ggsave(here("results/figures", "Fx_lFFA_age.png"), plot = fx_lFFA_age, width = 10, height = 5, units = "in", dpi = 300)
  1773. ```
  1774. ## Section IV: Age effects on lateralisation of face response (non-preregistered)
  1775. ### Descriptive statistics: Mean LI values
  1776. ```{r}
  1777. # Read data
  1778. child_LI <- read_csv(here("processed_data", "child_MRIvariables.csv"), col_names=TRUE)
  1779. child_LI <- child_LI[,c(1, 230:231)]
  1780. # Define variables to analyse
  1781. li_columns <- c("lateralisation_index_0.05", "lateralisation_index_0.1")
  1782. # Compute sunnary + t-tests
  1783. li_summary <- li_columns %>%
  1784. lapply(function(col_name) {
  1785. values <- child_LI[[col_name]]
  1786. # Compute mean and SD
  1787. mean_val <- mean(values, na.rm = TRUE)
  1788. sd_val <- sd(values, na.rm = TRUE)
  1789. # One-sample t-test against 0
  1790. t_test <- t.test(values, mu = 0, na.rm = TRUE)
  1791. # Format p-value
  1792. p_val <- ifelse(t_test$p.value < 0.001,
  1793. format(t_test$p.value, scientific = TRUE, digits = 3),
  1794. sprintf("%.3f", t_test$p.value))
  1795. # Create output row
  1796. tibble(
  1797. Measure = col_name,
  1798. mean_sd = sprintf("%.2f (%.2f)", mean_val, sd_val),
  1799. t_stat = sprintf("%.2f", t_test$statistic),
  1800. df = sprintf("%.0f", t_test$parameter),
  1801. p_value = p_val
  1802. )
  1803. }) %>%
  1804. bind_rows() # Combine into a single table
  1805. # View final table
  1806. li_summary
  1807. # Save the table
  1808. write_csv(li_summary, here("results/RQ1", "LI_descriptive_table.csv"))
  1809. ```
  1810. ### Age effects on lateralisation of face response
  1811. #### Read data
  1812. ```{r}
  1813. # Covariates
  1814. participants <- read_tsv(here("raw_data/covariates", "participants.tsv"), col_names=TRUE)
  1815. colnames(participants)[which(names(participants) == "participant_id")] <- "ID"
  1816. ## Remove the word "pixar" and trim extra whitespace
  1817. participants$ID <- gsub("[sub-pixar-]","", participants$ID)
  1818. participants$ID <- trimws(participants$ID)
  1819. # Motion
  1820. motion <- read_tsv(here("raw_data/motion", "outlier_info.tsv"), col_names=TRUE)
  1821. colnames(motion)[which(names(motion) == "subject")] <- "ID"
  1822. ## Remove the word "pixar" and trim extra whitespace
  1823. motion$ID <- gsub("[sub-pixar-]","", motion$ID)
  1824. motion$ID <- trimws(motion$ID)
  1825. # Child fMRI variables
  1826. child_LI <- read_csv(here("processed_data", "child_MRIvariables.csv"), col_names=TRUE)
  1827. # Complete dataset
  1828. child_LI <- child_LI[,c(1, 230:231)]
  1829. child_LI <- merge(child_LI, participants[,c(1:2)], by="ID")
  1830. child_LI <- merge(child_LI, motion[,c(4,6)], by="ID")
  1831. child_LI <- child_LI %>%
  1832. rename_with(~ sub("^zscored_", "", .x), starts_with("zscored_"))
  1833. child_LI[, c("lateralisation_index_0.05","lateralisation_index_0.1", "Age", "MeanFD")] <- scale(child_LI[, c("lateralisation_index_0.05","lateralisation_index_0.1", "Age", "MeanFD")])
  1834. ```
  1835. #### Models
  1836. ##### Threshold 0.05
  1837. ```{r}
  1838. model_LI_0.05_age <- lm(lateralisation_index_0.05 ~ Age + MeanFD, data = child_LI)
  1839. summary_LI_0.05_age <- tidy(model_LI_0.05_age)
  1840. ```
  1841. ##### Threshold 0.1
  1842. ```{r}
  1843. model_LI_0.1_age <- lm(lateralisation_index_0.1 ~ Age + MeanFD, data = child_LI)
  1844. summary_LI_0.1_age <- tidy(model_LI_0.1_age)
  1845. ```
  1846. ##### Save models
  1847. ```{r}
  1848. all_summaries <- bind_rows(
  1849. summary_LI_0.05_age %>% mutate(threshold = "0.05"),
  1850. summary_LI_0.1_age %>% mutate(threshold = "0.1")
  1851. ) %>%
  1852. mutate(
  1853. p.value = as.numeric(p.value), # ensure numeric
  1854. estimate = round(estimate, 2),
  1855. std.error = round(std.error, 2),
  1856. statistic = round(statistic, 2),
  1857. p.value = if_else(
  1858. p.value < 0.001,
  1859. format(p.value, scientific = TRUE, digits = 4),
  1860. sprintf("%.3f", p.value)
  1861. )
  1862. )
  1863. write.csv(all_summaries, here("results/RQ1", "LI_age.csv"), row.names = FALSE)
  1864. ```
  1865. #### Model plots
  1866. LI and age
  1867. ```{r}
  1868. # Read data
  1869. child_LI <- read_csv(here("processed_data", "child_MRIvariables.csv"), col_names=TRUE)
  1870. child_LI <- child_LI[,c(1, 230:231)]
  1871. child_LI <- merge(child_LI, participants[,c(1:2)], by="ID")
  1872. child_LI <- merge(child_LI, motion[,c(4,6)], by="ID")
  1873. child_LI <- child_LI %>%
  1874. rename_with(~ sub("^zscored_", "", .x), starts_with("zscored_"))
  1875. # LI of children at threshold p<0.05
  1876. ## Plot
  1877. LI_age_children <- ggplot(
  1878. child_LI,
  1879. aes(x = Age, y = lateralisation_index_0.05)
  1880. ) +
  1881. geom_point(color = "#9463A7", alpha = 0.8, size = 2.9, shape = 19) +
  1882. geom_smooth(method = "lm", se = FALSE, color = "#9463A7", linewidth = 2) +
  1883. labs(
  1884. title = "Lateralisation Index of FFA Across Age",
  1885. x = "Age",
  1886. y = "Lateralisation Index (p<0.05)"
  1887. ) +
  1888. theme_minimal() +
  1889. theme(
  1890. plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
  1891. axis.title.x = element_text(size = 16, face = "bold"),
  1892. axis.title.y = element_text(size = 16, face = "bold"),
  1893. axis.text = element_text(size = 12)
  1894. )
  1895. ggsave(here("results/figures", "LI_0.05_age.png"), plot = LI_age_children, width = 5, height = 5, units = "in", dpi = 300)
  1896. # LI of children at threshold p<0.01
  1897. LI_age_children <- ggplot(
  1898. child_LI,
  1899. aes(x = Age, y = lateralisation_index_0.1)
  1900. ) +
  1901. geom_point(color = "#9463A7", alpha = 0.8, size = 2.9, shape = 19) +
  1902. geom_smooth(method = "lm", se = FALSE, color = "#9463A7", linewidth = 2) +
  1903. labs(
  1904. title = "Lateralisation Index of FFA Across Age",
  1905. x = "Age",
  1906. y = "Lateralisation Index (p<0.10)"
  1907. ) +
  1908. theme_minimal() +
  1909. theme(
  1910. plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
  1911. axis.title.x = element_text(size = 16, face = "bold"),
  1912. axis.title.y = element_text(size = 16, face = "bold"),
  1913. axis.text = element_text(size = 12)
  1914. )
  1915. ggsave(here("results/figures", "LI_0.1_age.png"), plot = LI_age_children, width = 5, height = 5, units = "in", dpi = 300)
  1916. ```
  1917. ```{r}
  1918. markdownToHTML("FRIC_DevelopmentalChange.Rmd",output ="FRIC_DevelopmentalChange.html")
  1919. ```

FRIC_DevelopmentalChange.Rmd at commit 9982f13, no license · at the source

Overview

Authors: Lorena Jiménez-Sánchez1, Melissa Thye1, Hilary Richardson1
  1. Department of Psychology, School of Philosophy, Psychology and Language Sciences, University of Edinburgh, Edinburgh, UK
Institutions: University of Edinburgh (United Kingdom)
Journal: Developmental cognitive neuroscience, volume 80, article 101765
Dates: received 12 March 2026; accepted 16 June 2026; published online 18 June 2026; in print August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.dcn.2026.101765 · PMID 42322775 · PMCID PMC13312595 · OpenAlex W7165176499
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: fMRI (modality), human (organism)
Methods: Statistics, Preprocessing, Connectivity, fMRI & imaging
Keywords: Face selectivity, development, social cognition, fusiform face area, paediatric fMRI
MeSH: Amygdala*, Facial Recognition*, Prefrontal Cortex*, Social Perception*, Temporal Lobe*, Adult, Brain Mapping, Child, Child, Preschool, Cross-Sectional Studies, Female, Humans, Magnetic Resonance Imaging, Male, Young Adult (* major topic)
Topic: Face Recognition and Perception (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: University of Edinburgh (108890/Z/15/Z); Wellcome Trust (108890/Z/15/Z)
Citations: not cited yet (Europe PMC); 82 references in the paper

Abstract

The fusiform face area (FFA) preferentially responds to faces within the first months of life. One hypothesis is that higher-order social responses in middle medial prefrontal cortex (MMPFC) or face responses in superior temporal sulcus (STS) drive the development of face-selective responses in FFA, with right-hemisphere dominance in FFA eventually arising from lateralised connections to these regions. Another hypothesis proposes an innate face template in the amygdala guides attention to face-like shapes. This study opportunistically examined the development of the FFA, MMPFC, STS, and amygdala in childhood using an open cross-sectional movie-viewing fMRI dataset with 3–12-year-olds (N = 117, M = 6.77 years) and adults (N = 33, M = 24.77 years). We tested for correlations between FFA development and development in MMPFC, STS, and amygdala on the premise that associations between these regions may be observable even in children, and such associations could constrain hypotheses and analytic approaches in future studies with infants. First, we measured functional maturity- i.e., how similar each child’s response to the movie was to an adult average response timecourse. In all regions, older children’s responses were more adult-like. Next, we tested whether FFA maturity correlated with functional connectivity with, or functional maturity of, MMPFC, STS, or amygdala. Children with more mature right FFA responses showed stronger right FFA-right MMPFC connectivity. Children with more mature FFA responses also had more mature STS responses, bilaterally. This study provides preliminary evidence that FFA co-develops with higher-order social brain regions.

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

Repositories

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

hrichardsonlab/fmri-analysis

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 15a4d88696a4daa910c428d636a1d5ad90c581d5, 10 September 2026
Languages: Python (24), Shell (18), R (1)
Size: 310 files, 43 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, 1 notebook
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (24 files), pandas (24 files), Nipype (14 files), Nilearn (11 files), NiBabel (9 files), SciPy (8 files), FSL (6 files), PyBIDS (5 files), fMRIPrep (4 files), Dcm2Bids (2 files), FreeSurfer (2 files), ANTs (1 file), Matplotlib (1 file), MRIQC (1 file), scikit-learn (1 file), tedana (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
44 files

LorenaJS/FRIC

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 9982f13ac4e7399ce88801d7264b89bb0c45139c, 22 May 2026
Languages: R (5)
Size: 18 files, 5 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, 5 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (5 files), ggpubr (4 files), patchwork (4 files), broom (3 files), car (3 files), lme4 (3 files), lmerTest (3 files), reshape2 (3 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
6 files

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 48 scripts, each with its path and the digest of its content;
  • 10 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

Datasets cited

Data availability

The movie fMRI data analysed during the current study (originally collected by Richardson et al., 2018) are publicly available on OpenNeuro (https://openneuro.org/datasets/ds000228). Further details about the fMRI analysis pipeline, and publicly available code, can be found at https://github.com/hrichardsonlab/fmri-analysis. Average timecourses from the referent adult populations, used as regressors in this study, are available through OSF (https://osf.io/7a8w5/ for FFA and LOC, and https://osf.io/mxkag/ for MMPFC and S2). Search spaces are also available through OSF (https://osf.io/mxkag/ for MMPFC and S2 split by hemisphere, and left FFA defined by flipping the right FFA search space). Data and scripts used for this analysis can be accessed at https://github.com/LorenaJS/FRIC.

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, 3 authors, 5 keywords, 15 MeSH terms, 2 funders, 78 references.

Cite

This paper

Jiménez-Sánchez, L., Thye, M., & Richardson, H. (2026). Fusiform face area development correlates with development in higher-order social brain regions. Developmental cognitive neuroscience, 80, 101765. https://doi.org/10.1016/j.dcn.2026.101765

BibTeX

@article{jimenezsanchez2026fusiform,
author = {Jiménez-Sánchez, Lorena and Thye, Melissa and Richardson, Hilary},
title = {{Fusiform face area development correlates with development in higher-order social brain regions}},
journal = {Developmental cognitive neuroscience},
year = {2026},
month = jun,
volume = {80},
pages = {101765},
publisher = {Elsevier},
issn = {1878-9293},
doi = {10.1016/j.dcn.2026.101765},
url = {https://doi.org/10.1016/j.dcn.2026.101765},
pmid = {42322775},
pmcid = {PMC13312595}
}

RIS

TY - JOUR
AU - Jiménez-Sánchez, Lorena
AU - Thye, Melissa
AU - Richardson, Hilary
TI - Fusiform face area development correlates with development in higher-order social brain regions
T2 - Developmental cognitive neuroscience
J2 - Dev Cogn Neurosci
PY - 2026
DA - 2026/06/18
VL - 80
SP - 101765
SN - 1878-9293
PB - Elsevier
DO - 10.1016/j.dcn.2026.101765
UR - https://doi.org/10.1016/j.dcn.2026.101765
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.dcn.2026.101765",
"type": "article-journal",
"title": "Fusiform face area development correlates with development in higher-order social brain regions",
"container-title": "Developmental cognitive neuroscience",
"author": [
{
"family": "Jiménez-Sánchez",
"given": "Lorena"
},
{
"family": "Thye",
"given": "Melissa"
},
{
"family": "Richardson",
"given": "Hilary"
}
],
"container-title-short": "Dev Cogn Neurosci",
"volume": "80",
"page": "101765",
"DOI": "10.1016/j.dcn.2026.101765",
"PMID": "42322775",
"PMCID": "PMC13312595",
"ISSN": "1878-9293",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.dcn.2026.101765",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
18
]
]
}
}

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

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1162/imag.a.1347 [code]
Neural and behavioural correlates of theory of mind reasoning in five-year-old children born preterm.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: MRIQC, Dcm2Bids, fMRIPrep, 18 other tools, OpenNeuro ds000228, fMRI, 5 references, author Lorena Jiménez-Sánchez
[2] doi:10.1162/imag.a.1252 [code]
Does the brain's E:I balance really shape long-range temporal correlations? Lessons learned from 3T MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: MRIQC, Dcm2Bids, fMRIPrep, 12 other tools, fMRI, 1 reference
[3] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: fMRIPrep, Nipype, ANTs, 9 other tools, 5 references
[4] doi:10.1162/netn.a.547 [code]
An evaluation of the efficacy of single-echo and multi-echo fMRI denoising strategies.
Journal: Network neuroscience (Cambridge, Mass.)
In common: fMRIPrep, tedana, Nipype, 10 other tools, fMRI, 2 references
[5] doi:10.1162/imag.a.1198 [code]
MEPrep: A robust pipeline for multi-echo fMRI denoising and preprocessing.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: fMRIPrep, tedana, PyBIDS, 9 other tools, fMRI, 2 references
[6] doi:10.1162/imag.a.1245 [code]
Towards precision EEG connectomics: Evaluating the benefits of dense sampling.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nipype, ANTs, FreeSurfer, 13 other tools
[7] doi:10.1038/s42003-026-10843-3 [code]
Self-supervised learning yields representational signatures of category-selective cortex.
Journal: Communications biology
In common: 13 references
[8] doi:10.1093/nc/niag029 [code]
A data-driven approach to identifying and evaluating connectivity-based neural correlates of conscious visual perception.
Journal: Neuroscience of consciousness
In common: car, broom, FreeSurfer, 12 other tools, 1 reference
[9] doi:10.1038/s41467-026-71830-0 [code]
Predicting individual differences of fear and cognitive learning and extinction.
Journal: Nature communications
In common: Nipype, ANTs, car, 11 other tools
[10] doi:10.1016/j.dcn.2026.101703 [code]
Visual Word Form Area demonstrates individual and task-agnostic consistency but inter-individual variability.
Journal: Developmental cognitive neuroscience
In common: FreeSurfer, lmerTest, Nilearn, 7 other tools, fMRI, 5 references

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.