Precision fMRI reveals that the language network exhibits adult-like left-hemispheric lateralization by 4 years of age.

Code ↔ Paper

2 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 2 matches

[1] § Methods › Do response magnitude and inter-regional correlations exhibit age-related changes? ↔ source data and code/OzernovPalchik_OBrien_source_data_code_FINAL_100325.zip/code/NCOMM_OzernovPalchik_OBrien_SI_code_091025.Rmd, lines 3219–3247 · score 0.61 · correlation strength, inter regional, pairwise comparisons, network, models, adult
[2] § Methods › Do children (of different ages) show left-hemispheric bias, and is the bias adult-like? ↔ source data and code/OzernovPalchik_OBrien_source_data_code_FINAL_100325.zip/code/NCOMM_OzernovPalchik_OBrien_main_code_091025.Rmd, lines 707–791 · score 0.57 · linear regression model, predicting LI, LH bias, fit, volume, RH

Paper

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

The paper is loaded when this pane is shown.

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

The authors' code

R Markdown · 3,309 lines · 112 KB · no license · 1 match

---
title: "NComm_SI_Alice_power_extracted"
output:
html_document:
toc: true
toc_float: true
theme: flatly
date: "2025-03-10"
editor_options:
chunk_output_type: console
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(
echo = TRUE,
warning = FALSE,
message = FALSE,
fig.width = 8,
fig.height = 6
)
```
###########################################
#SI 0
###########################################
# SI-0: Set Up
## Load Required Libraries
```{r load-libraries}
# Load required libraries
library(lme4)
library(lsmeans)
library(dplyr)
library(kableExtra)
library(tidyr)
library(MuMIn) # for r.squaredGLMM
library(pwr) # for power analysis
library(ggplot2)
library(ggsignif)
library(ggpubr)
library(cowplot)
library(stringr)
library(brms)
# Note: Ensure your working directory is set to the code/ folder
# or that data files are accessible via ../data/ relative path
```
## Load Datasets
```{r load-data}
# Read in datasets from ../data/ folder
ds1_eff <- read.csv("../data/effect_ds1_090125.csv")
ds2_eff <- read.csv("../data/effect_ds2_031025.csv")
ds1_vol <- read.csv("../data/volume_ds1_031025.csv")
ds2_vol <- read.csv("../data/volume_ds2_031025.csv")
ds1_demo <- read.csv("../data/ds1_demo_data_031425.csv")
ds2_demo <- read.csv("../data/ds2_demo_090225.csv")
beh <- read.csv('../data/MNM_behavfMRI_180825.csv')
cat("Data loaded successfully!\n")
cat("DS1 Effect data dimensions:", dim(ds1_eff), "\n")
cat("DS2 Effect data dimensions:", dim(ds2_eff), "\n")
# Report unique N for each dataset
cat("\nUnique subject counts:\n")
cat("DS1 Effect N =", length(unique(ds1_eff$Subject)), "\n")
cat("DS2 Effect N =", length(unique(ds2_eff$Subject)), "\n")
cat("DS1 Volume N =", length(unique(ds1_vol$Subject)), "\n")
cat("DS2 Volume N =", length(unique(ds2_vol$Subject)), "\n")
cat("DS1 Demo N =", length(unique(ds1_demo$Subject)), "\n")
cat("DS2 Demo N =", length(unique(ds2_demo$Subject)), "\n")
cat("DS2 Beh N =", length(unique(beh$Subject)), "\n")
```
## Convert Variables to Factors
```{r factor-conversion}
# Define columns to convert to factors for effect datasets
cols_to_factor <- c("age", "front", "hemi", "ROI")
# Convert specified columns to factors in both datasets
ds1_eff[cols_to_factor] <- lapply(ds1_eff[cols_to_factor], as.factor)
ds2_eff[cols_to_factor] <- lapply(ds2_eff[cols_to_factor], as.factor)
ds1_vol[cols_to_factor] <- lapply(ds1_vol[cols_to_factor], as.factor)
ds2_vol[cols_to_factor] <- lapply(ds2_vol[cols_to_factor], as.factor)
cat("Factor conversion completed.\n")
cat("Age levels in DS1:", levels(ds1_eff$age), "\n")
cat("ROI levels in DS1:", length(levels(ds1_eff$ROI)), "regions\n")
cat("ROI levels in DS2:", length(levels(ds2_eff$ROI)), "regions\n")
```
## ROI Name Mapping for Dataset 2
```{r roi-mapping}
# Create mapping vector for ROI names
roi_map <- c(
LH_AntTemp = "AntTemp",
LH_IFG = "IFG",
LH_IFGorb = "IFGorb",
LH_MFG = "MFG",
LH_PostTemp = "PostTemp",
RH_AntTemp = "R AntTemp",
RH_IFG = "R IFG",
RH_IFGorb = "R IFG orb",
RH_MFG = "R MFG",
RH_PostTemp = "R Post temp"
)
# Function to update DS2 ROI names to match DS1 format
update_ds2_rois <- function(df_ds2, df_ds1) {
df_ds2$ROI <- factor(
roi_map[as.character(df_ds2$ROI)], # map old → new names
levels = levels(df_ds1$ROI) # same order/levels as ds1
)
return(df_ds2)
}
# Apply ROI name updates
ds2_vol <- update_ds2_rois(ds2_vol, ds1_vol)
ds2_eff <- update_ds2_rois(ds2_eff, ds1_eff)
# Verify that ROI levels match between datasets
stopifnot(
all(levels(ds2_vol$ROI) == levels(ds1_vol$ROI)),
all(levels(ds2_eff$ROI) == levels(ds1_eff$ROI))
)
cat("ROI mapping completed successfully.\n")
cat("ROI levels now match between datasets.\n")
```
## Combine Datasets
```{r combine-datasets}
# Combine effect datasets
ds1_eff$ds <- '1'
ds2_eff$ds <- '2'
ds_eff_combined <- bind_rows(ds1_eff, ds2_eff)
ds_eff_combined$age_years <- NULL
# Combine volume datasets
ds1_vol$ds <- '1'
ds2_vol$ds <- '2'
ds_vol_combined <- bind_rows(ds1_vol, ds2_vol)
ds_vol_combined$group <- NULL
cat("Datasets combined successfully.\n")
cat("Combined effect data dimensions:", dim(ds_eff_combined), "\n")
cat("Combined volume data dimensions:", dim(ds_vol_combined), "\n")
# Check for NAs in key variables
cat("\nNA counts:\n")
cat("Effect data NAs:", sum(is.na(ds_eff_combined$effect)), "\n")
cat("Volume data NAs:", sum(is.na(ds_vol_combined$volume)), "\n")
# Report unique N by age for effect data
cat("\nEffect data - Unique N by age:\n")
age_counts_eff <- ds_eff_combined %>%
distinct(Subject, age) %>%
count(age, name = "N")
print(age_counts_eff)
# Report unique N by age for volume data
cat("\nVolume data - Unique N by age:\n")
age_counts_vol <- ds_vol_combined %>%
distinct(Subject, age) %>%
count(age, name = "N")
print(age_counts_vol)
```
## Create Right Hemisphere Subsets
```{r rh-subsets}
# Dataset 1 - Right Hemisphere
ds1_eff_rh <- ds1_eff %>%
filter(hemi == "rh", condition == "language") %>%
select(Subject, ROI, age, front, hemi, contrast) %>%
mutate(dataset = "DS1")
# Create age group subsets for DS1
ds1_early_rh <- filter(ds1_eff_rh, age == "early")
ds1_middle_rh <- filter(ds1_eff_rh, age == "middle")
ds1_late_rh <- filter(ds1_eff_rh, age == "late")
ds1_adult_rh <- filter(ds1_eff_rh, age == "adult")
# Dataset 2 - Right Hemisphere
ds2_eff_rh <- ds2_eff %>%
filter(hemi == "rh", condition == "EffectSize.mentalsocialphysical") %>%
select(Subject, ROI, age, front, hemi, contrast) %>%
mutate(dataset = "DS2")
# Create age group subsets for DS2
ds2_early_rh <- filter(ds2_eff_rh, age == "early")
ds2_middle_rh <- filter(ds2_eff_rh, age == "middle")
ds2_late_rh <- filter(ds2_eff_rh, age == "late")
ds2_adult_rh <- filter(ds2_eff_rh, age == "adult")
```
###########################################
#SI 1
###########################################
# SI-1.B: Behavioral performance of Dataset2
## Load and Prepare Behavioral Data
```{r behavioral-prep}
# Get subjects from neural data
neural_subjects <- unique(ds2_eff$Subject)
# Get subjects from behavioral data
behavioral_subjects <- unique(beh$Subject)
# Find subjects in neural data but missing from behavioral data
setdiff(neural_subjects, behavioral_subjects)
# Merge datasets (only subjects with both neural and behavioral data)
beh_merged <- beh %>%
inner_join(ds2_eff %>% dplyr::select(Subject, age) %>% distinct(), by = "Subject")
colSums(is.na(beh_merged))
# Report final sample
cat("Final sample size:", length(unique(beh_merged$Subject)), "subjects\n")
```
## Set Up Custom Contrasts for Age Groups
```{r age-contrasts}
# # Create custom contrasts matrix for age group comparisons
# custom_contrasts <- matrix(c(
# 1, 0, 0, # early vs middle
# -1, 1, 0, # middle vs late
# 0, -1, 1, # late vs adult
# 0, 0, -1 # adult baseline
# ), ncol = 3, byrow = TRUE)
#
# rownames(custom_contrasts) <- c("early", "middle", "late", "adult")
# colnames(custom_contrasts) <- c("early_vs_middle", "middle_vs_late", "late_vs_adult")
#
# # Apply contrasts to age factor
# beh_merged$age <- as.factor(beh_merged$age)
# contrasts(beh_merged$age) <- custom_contrasts
#
# #check contrasts
# kable(custom_contrasts, caption = "Custom Age Group Contrast Matrix") %>%
# kable_styling(bootstrap_options = c("striped", "hover"))
```
## SI-1B Behavioral Performance Summary
### Summary behavioral data
```{r behavioral-summary}
# Create participant info dataframe
participant_info <- beh_merged %>%
dplyr::select(Subject, age) %>%
distinct()
# Aggregate behavioral performance (excluding Music condition)
agg_summary <- beh_merged %>%
filter(Cond != "Music") %>%
mutate(CondGroup = case_when(
Cond == "Foreign" ~ "Foreign",
Cond %in% c("Mental", "Physical", "Social") ~ "MPS"
)) %>%
filter(!is.na(CondGroup)) %>%
group_by(Subject, CondGroup) %>%
summarise(
Total = n(),
Sum_Correct = sum(Correct, na.rm = TRUE),
Sum_Wrong = sum(Wrong, na.rm = TRUE),
Sum_Miss = sum(Miss, na.rm = TRUE),
.groups = "drop"
) %>%
left_join(participant_info, by = "Subject")
# Calculate percent correct for each subject and condition
agg_per <- agg_summary %>%
mutate(Percent_Correct = (Sum_Correct / Total) * 100)
# Display summary statistics
summary_stats <- agg_per %>%
group_by(CondGroup, age) %>%
summarise(
N = n(),
Mean_PC = round(mean(Percent_Correct, na.rm = TRUE), 2),
SD_PC = round(sd(Percent_Correct, na.rm = TRUE), 2),
SE_PC = round(sd(Percent_Correct, na.rm = TRUE) / sqrt(n()), 2),
.groups = "drop"
)
kable(summary_stats,
caption = "Behavioral Performance Summary by Condition and Age Group",
col.names = c("Condition", "Age", "N", "Mean %", "SD", "SE")) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
```
### Analysis of behavioral data
```{r behavioral-analysis}
# Test age and condition effects using mixed-effects model
m_beh_1 <- lmer(Percent_Correct ~ CondGroup * age + (1 | Subject), data = agg_per)
# Display model summary
cat("Mixed-Effects Model Results:\n")
cat("============================\n")
print(summary(m_beh_1))
# Get estimated marginal means and all pairwise comparisons for age
age_emm <- emmeans(m_beh_1, ~ age)
age_contrasts <- pairs(age_emm, adjust = "tukey")
print(age_contrasts)
# Calculate R-squared
r2_values <- r.squaredGLMM(m_beh_1)
cat("\nR-squared values:\n")
cat("Marginal R²:", round(r2_values[1], 3), "\n")
cat("Conditional R²:", round(r2_values[2], 3), "\n")
# Post-hoc comparisons for condition groups
cat("Post-hoc comparisons for Condition Groups:\n")
cat("==========================================\n")
condition_comparisons <- emmeans(m_beh_1, pairwise ~ CondGroup, adjust = "tukey")
print(condition_comparisons)
# Age group comparisons - formatted as requested table
cat("\n\n**Accuracy on the in-scanner task**\n")
cat("====================================\n")
# Get age group comparisons
age_comparisons <- emmeans(m_beh_1, pairwise ~ age, adjust = "tukey")
age_contrasts <- as.data.frame(age_comparisons$contrasts)
# Create formatted table
age_results <- age_contrasts %>%
mutate(
Comparison = gsub(" - ", " vs. ", contrast),
B = round(estimate, 1),
SE = round(SE, 2),
t = round(t.ratio, 2),
p = round(p.value, 2),
p_formatted = case_when(
p.value < 0.001 ~ "**<0.001**",
p.value < 0.01 ~ paste0("**", sprintf("%.2f", p.value), "**"),
p.value < 0.05 ~ paste0("**", sprintf("%.2f", p.value), "**"),
TRUE ~ sprintf("%.2f", p.value)
)
) %>%
select(Comparison, B, SE, t, p_formatted) %>%
arrange(Comparison)
# Display as kable
kable(age_results,
col.names = c("", "***B***", "***SE***", "***t***", "***p***"),
caption = "**Accuracy on the in-scanner task**",
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
# Analysis within language condition only (excluding Foreign)
agg_per_lang <- agg_per %>% filter(CondGroup != "Foreign")
# Fit linear model for language condition
lm_lang <- lm(Percent_Correct ~ age, data = agg_per_lang)
cat("\n\nLanguage Condition Analysis (MPS only):\n")
cat("=======================================\n")
print(summary(lm_lang))
# Calculate effect sizes
anova_results <- anova(lm_lang)
eta_squared <- anova_results$"Sum Sq"[1] / sum(anova_results$"Sum Sq")
cat("\nEta-squared (effect size):", round(eta_squared, 3), "\n")
```
###########################################
#SI-2: The Language Network’s Topography in Children
###########################################
#SI-2B Within- and between-participant similarity of activation patterns for the Language > Control contrast.
# Complete SI-2B: Within- and between-participant similarity analysis with tables
```{r}
## 1. WITHIN-PARTICIPANT CORRELATIONS
# A. Whole brain within-participant correlations
within_corr <- read.csv("../data/corr_analyses/alice_correlations_within.csv")
# Ensure Subject column exists
if("Subject" %in% names(within_corr)) {
# Already has Subject column
} else if("Participant" %in% names(within_corr)) {
names(within_corr)[names(within_corr) == "Participant"] <- "Subject"
} else {
# Assume first column is Subject ID
names(within_corr)[1] <- "Subject"
}
within_corr$group <- factor(tolower(within_corr$Group), levels = c("early", "middle", "late", "adult"))
# Statistical models - whole brain within
m_within_wb <- lm(Corr ~ group, data = within_corr)
# Continuous age analysis (children only)
cont_within_wb <- merge(within_corr, ds1_demo, by = "Subject") %>% filter(Set != 'ADULT')
wb_within_cont <- lm(Corr ~ Age, data = cont_within_wb)
# B. ROI-based within-participant correlations (LH language network)
all_sp_corr <- read.csv("../data/corr_analyses/sp_corr_data.csv")
all_sp_corr_lh <- all_sp_corr %>%
filter(variable %in% c("IFG", "IFGorb", "MFG", "AntTemp", "PostTemp")) %>%
group_by(Participant, group) %>%
summarise(corr = mean(value, na.rm = TRUE), .groups = "drop") %>%
rename(Subject = Participant)
all_sp_corr_lh$group <- factor(all_sp_corr_lh$group, levels = c("early", "middle", "late", "adult"))
# Statistical models - ROI within
m_within_roi <- lm(corr ~ group, data = all_sp_corr_lh)
# Continuous age analysis (children only)
cont_within_roi <- merge(all_sp_corr_lh, ds1_demo, by = "Subject") %>%
filter(group != 'adult')
roi_within_cont <- lm(corr ~ Age, data = cont_within_roi)
## 2. BETWEEN-PARTICIPANT CORRELATIONS
# Process between-participant data function
process_between_data <- function(df, group_name) {
df$r <- rowMeans(df[, -1], na.rm = TRUE)
df$group <- group_name
# Standardize first column name to "Subject"
if("Subject" %in% names(df)) {
# Already has Subject column
} else if("Participant" %in% names(df)) {
names(df)[names(df) == "Participant"] <- "Subject"
} else if("participant" %in% names(df)) {
names(df)[names(df) == "participant"] <- "Subject"
} else {
# Assume first column is Subject ID
names(df)[1] <- "Subject"
}
if(group_name == "adult") df$Subject <- gsub("_", "", df$Subject)
return(df[, c("Subject", "r", "group")])
}
# C. Whole brain between-participant correlations
df1_bt2_avg <- read.csv("../data/corr_analyses/alice_correlations_between_early_2.csv")
df2_bt2_avg <- read.csv("../data/corr_analyses/alice_correlations_between_middle_2_all.csv")
df3_bt2_avg <- read.csv("../data/corr_analyses/alice_correlations_between_late_2.csv")
df4_bt2_avg <- read.csv("../data/corr_analyses/alice_correlations_between_adult_2_updated_092823.csv")
all_between_wb <- rbind(
process_between_data(df1_bt2_avg, "early"),
process_between_data(df2_bt2_avg, "middle"),
process_between_data(df3_bt2_avg, "late"),
process_between_data(df4_bt2_avg, "adult")
)
all_between_wb$group <- factor(all_between_wb$group, levels = c("early", "middle", "late", "adult"))
# Statistical models - whole brain between
m_between_wb <- lm(r ~ group, data = all_between_wb)
# Continuous age analysis (children only)
cont_between_wb <- merge(all_between_wb, ds1_demo, by = "Subject") %>%
filter(group != 'adult')
wb_between_cont <- lm(r ~ Age, data = cont_between_wb)
# D. ROI-based between-participant correlations
df1_roi_avg <- read.csv("../data/corr_analyses/ROI/avg_corr_matrix_early_lh.csv")
df2_roi_avg <- read.csv("../data/corr_analyses/ROI/avg_corr_matrix_middle_lh.csv")
df3_roi_avg <- read.csv("../data/corr_analyses/ROI/avg_corr_matrix_late_lh.csv")
df4_roi_avg <- read.csv("../data/corr_analyses/ROI/avg_corr_matrix_adult_lh_updated092823.csv")
all_between_roi <- rbind(
process_between_data(df1_roi_avg, "early"),
process_between_data(df2_roi_avg, "middle"),
process_between_data(df3_roi_avg, "late"),
process_between_data(df4_roi_avg, "adult")
)
all_between_roi$group <- factor(all_between_roi$group, levels = c("early", "middle", "late", "adult"))
# Statistical models - ROI between
m_between_roi <- lm(r ~ group, data = all_between_roi)
# Continuous age analysis (children only)
cont_between_roi <- merge(all_between_roi, ds1_demo, by = "Subject") %>%
filter(group != 'adult')
roi_between_cont <- lm(r ~ Age, data = cont_between_roi)
## DEBUGGING SECTION - ADD THIS TO TROUBLESHOOT
# Debug function to check continuous models
debug_continuous_model <- function(model, model_name, data) {
cat("\n=== Debugging", model_name, "===\n")
# Check if Age variable exists and its properties
cat("Age variable summary:\n")
if("Age" %in% names(data)) {
print(summary(data$Age))
cat("Age class:", class(data$Age), "\n")
cat("Age range:", range(data$Age, na.rm = TRUE), "\n")
cat("Any NAs in Age:", sum(is.na(data$Age)), "\n")
} else {
cat("WARNING: No 'Age' variable found in data!\n")
cat("Available variables:", names(data), "\n")
}
# Check model summary
cat("\nModel summary:\n")
model_summary <- summary(model)
print(model_summary)
# Check coefficient names
cat("\nCoefficient names:\n")
print(rownames(model_summary$coefficients))
# Try to extract Age coefficient
cat("\nAge coefficient extraction:\n")
if("Age" %in% rownames(model_summary$coefficients)) {
age_coef <- model_summary$coefficients["Age", ]
print(age_coef)
} else {
cat("WARNING: 'Age' not found in model coefficients!\n")
}
cat("========================\n\n")
}
# Check data merging for continuous models
cat("=== CHECKING DATA MERGING ===\n")
cat("ds1_demo columns:", names(ds1_demo), "\n")
if("Age" %in% names(ds1_demo)) {
cat("ds1_demo Age summary:", summary(ds1_demo$Age), "\n")
} else {
cat("WARNING: No Age column in ds1_demo!\n")
}
# Check each merged dataset
check_merged_data <- function(merged_data, name) {
cat("\n", name, ":\n")
cat("Dimensions:", dim(merged_data), "\n")
cat("Age column exists:", "Age" %in% names(merged_data), "\n")
if("Age" %in% names(merged_data)) {
cat("Age summary:", summary(merged_data$Age), "\n")
}
}
check_merged_data(cont_within_roi, "cont_within_roi")
check_merged_data(cont_within_wb, "cont_within_wb")
check_merged_data(cont_between_roi, "cont_between_roi")
check_merged_data(cont_between_wb, "cont_between_wb")
# Run debugging for all continuous models
debug_continuous_model(roi_within_cont, "ROI Within Continuous", cont_within_roi)
debug_continuous_model(wb_within_cont, "WB Within Continuous", cont_within_wb)
debug_continuous_model(roi_between_cont, "ROI Between Continuous", cont_between_roi)
debug_continuous_model(wb_between_cont, "WB Between Continuous", cont_between_wb)
# Robust coefficient extraction function
extract_age_coefficient <- function(model) {
model_summary <- summary(model)
coef_matrix <- model_summary$coefficients
# Check different possible names for Age variable
age_names <- c("Age", "age", "AGE")
age_row <- NULL
for(name in age_names) {
if(name %in% rownames(coef_matrix)) {
age_row <- name
break
}
}
if(is.null(age_row)) {
cat("WARNING: No Age variable found in model coefficients.\n")
cat("Available coefficients:", rownames(coef_matrix), "\n")
return(c(Estimate = 0, `Std. Error` = 0, `t value` = 0, `Pr(>|t|)` = 1))
}
return(coef_matrix[age_row, ])
}
## END DEBUGGING SECTION
## 3. TABLE GENERATION FUNCTIONS
# Function to format p-values with bold for significance
format_p <- function(p) {
ifelse(p < 0.001, "**<0.001**",
ifelse(p < 0.01, paste0("**", sprintf("%.3f", p), "**"),
ifelse(p < 0.05, paste0("**", sprintf("%.3f", p), "**"),
sprintf("%.3f", p))))
}
# Function to create formatted table
create_results_table <- function(continuous_model, categorical_model, title) {
# Extract continuous age effect using robust method
cont_coef <- extract_age_coefficient(continuous_model)
# Extract pairwise comparisons vs Adult using emmeans
emm_result <- emmeans(categorical_model, ~ group)
contrasts_result <- pairs(emm_result, adjust = "sidak")
contrast_df <- as.data.frame(contrasts_result)
# Find comparisons against adult (adult should be reference)
# Look for patterns like "early - adult", "middle - adult", "late - adult"
early_vs_adult <- contrast_df[grepl("early.*-.*adult", contrast_df$contrast), ]
middle_vs_adult <- contrast_df[grepl("middle.*-.*adult", contrast_df$contrast), ]
late_vs_adult <- contrast_df[grepl("late.*-.*adult", contrast_df$contrast), ]
# If no matches found, try alternative approach with manual contrasts
if(nrow(early_vs_adult) == 0) {
# Create custom contrasts comparing each group to adult
contrast_list <- list(
"early_vs_adult" = c(1, 0, 0, -1), # early - adult
"middle_vs_adult" = c(0, 1, 0, -1), # middle - adult
"late_vs_adult" = c(0, 0, 1, -1) # late - adult
)
custom_contrasts <- contrast(emm_result, contrast_list, adjust = "sidak")
contrast_df <- as.data.frame(custom_contrasts)
early_vs_adult <- contrast_df[1, ]
middle_vs_adult <- contrast_df[2, ]
late_vs_adult <- contrast_df[3, ]
}
# Create results dataframe
results <- data.frame(
Comparison = c("Age = Continuous", "Early vs. Adult", "Middle vs. Adult", "Late vs. Adult"),
B = c(round(cont_coef["Estimate"], 2),
round(early_vs_adult$estimate, 2),
round(middle_vs_adult$estimate, 2),
round(late_vs_adult$estimate, 2)),
SE = c(round(cont_coef["Std. Error"], 2),
round(early_vs_adult$SE, 2),
round(middle_vs_adult$SE, 2),
round(late_vs_adult$SE, 2)),
t = c(round(cont_coef["t value"], 2),
round(early_vs_adult$t.ratio, 2),
round(middle_vs_adult$t.ratio, 2),
round(late_vs_adult$t.ratio, 2)),
p = c(format_p(cont_coef["Pr(>|t|)"]),
format_p(early_vs_adult$p.value),
format_p(middle_vs_adult$p.value),
format_p(late_vs_adult$p.value))
)
# Create kable table
kable(results,
col.names = c("", "***B***", "***SE***", "***t***", "***p***"),
caption = paste0("**", title, "**"),
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
}
# Alternative simple table generation (use if main function fails)
create_simple_table <- function(continuous_model, categorical_model, title) {
# Continuous age effect
cont_summary <- summary(continuous_model)
age_coef <- cont_summary$coefficients["Age", ]
# Get group means and use adult as reference
group_summary <- summary(categorical_model)
coef_matrix <- group_summary$coefficients
# Extract coefficients (adult is reference, so intercept = adult mean)
# Other groups show difference from adult
intercept <- coef_matrix["(Intercept)", ]
early_coef <- if("groupearly" %in% rownames(coef_matrix)) coef_matrix["groupearly", ] else c(0, 0, 0, 1)
middle_coef <- if("groupmiddle" %in% rownames(coef_matrix)) coef_matrix["groupmiddle", ] else c(0, 0, 0, 1)
late_coef <- if("grouplate" %in% rownames(coef_matrix)) coef_matrix["grouplate", ] else c(0, 0, 0, 1)
# Create results
results <- data.frame(
Comparison = c("Age = Continuous", "Early vs. Adult", "Middle vs. Adult", "Late vs. Adult"),
B = c(round(age_coef["Estimate"], 2),
round(early_coef["Estimate"], 2),
round(middle_coef["Estimate"], 2),
round(late_coef["Estimate"], 2)),
SE = c(round(age_coef["Std. Error"], 2),
round(early_coef["Std. Error"], 2),
round(middle_coef["Std. Error"], 2),
round(late_coef["Std. Error"], 2)),
t = c(round(age_coef["t value"], 2),
round(early_coef["t value"], 2),
round(middle_coef["t value"], 2),
round(late_coef["t value"], 2)),
p = c(format_p(age_coef["Pr(>|t|)"]),
format_p(early_coef["Pr(>|t|)"]),
format_p(middle_coef["Pr(>|t|)"]),
format_p(late_coef["Pr(>|t|)"]))
)
kable(results,
col.names = c("", "***B***", "***SE***", "***t***", "***p***"),
caption = paste0("**", title, "**"),
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
}
# Debug function to check contrasts (optional - remove if not needed)
debug_contrasts <- function(categorical_model, model_name) {
cat("\n=== Debug for", model_name, "===\n")
emm_result <- emmeans(categorical_model, ~ group)
cat("Estimated marginal means:\n")
print(emm_result)
contrasts_result <- pairs(emm_result, adjust = "sidak")
cat("All pairwise contrasts:\n")
print(contrasts_result)
cat("========================\n")
}
# Run debug for all models (comment out if not needed)
debug_contrasts(m_within_roi, "ROI Within")
debug_contrasts(m_within_wb, "Whole Brain Within")
debug_contrasts(m_between_roi, "ROI Between")
debug_contrasts(m_between_wb, "Whole Brain Between")
# Table A: Within-participant correlations (LH language network)
table_a <- create_results_table(
roi_within_cont, m_within_roi,
"A. Age differences in within-participant correlations (indexing activation stability across runs) within the LH language network."
)
# Table B: Within-participant correlations (whole brain)
cat("Generating Table B...\n")
table_b <- create_results_table(
wb_within_cont, m_within_wb,
"B. Age differences in within-participant correlations across the brain."
)
# Table C: Between-participant correlations (LH language network)
cat("Generating Table C...\n")
table_c <- create_results_table(
roi_between_cont, m_between_roi,
"C. Age differences in between-participant correlations (indexing inter-individual topographic variability) within the LH language network."
)
# Table D: Between-participant correlations (whole brain)
cat("Generating Table D...\n")
table_d <- create_results_table(
wb_between_cont, m_between_wb,
"D. Age differences in between-participant correlations across the brain."
)
# Display all tables
cat("## SI-2B Statistical Results Tables\n\n")
print(table_a)
cat("\n")
print(table_b)
cat("\n")
print(table_c)
cat("\n")
print(table_d)
## 5. COMBINED FIGURE (OPTIONAL)
# Prepare data for plotting
prepare_plot_data <- function(within_data, between_data, within_col, between_col, within_group_col, between_group_col) {
within_plot <- within_data %>%
select(all_of(c(within_group_col, within_col))) %>%
rename(group = !!within_group_col, Corr = !!within_col) %>%
mutate(comparison = "within")
between_plot <- between_data %>%
select(all_of(c(between_group_col, between_col))) %>%
rename(group = !!between_group_col, Corr = !!between_col) %>%
mutate(comparison = "between")
return(rbind(within_plot, between_plot))
}
# ROI data
roi_plot_data <- prepare_plot_data(all_sp_corr_lh, all_between_roi, "corr", "r", "group", "group")
# Whole brain data
wb_plot_data <- prepare_plot_data(within_corr, all_between_wb, "Corr", "r", "group", "group")
# Create plots
create_corr_plot <- function(data, title) {
data$combined <- paste(data$group, data$comparison, sep = "_")
data$combined <- factor(data$combined,
levels = c("early_within", "early_between", "middle_within", "middle_between",
"late_within", "late_between", "adult_within", "adult_between"))
ggbarplot(data, x = "combined", y = "Corr", fill = "comparison",
add = "mean_se", palette = c("gray33", "gray70"), width = 0.6,
ylab = "Correlations", ylim = c(0, 1), xlab = "", title = title) +
scale_x_discrete(labels = rep(c("w/in", "b/t"), 4)) +
theme_classic() + theme(text = element_text(size = 18))
}
roi_plot <- create_corr_plot(roi_plot_data, "A. LH Language Parcels")
wb_plot <- create_corr_plot(wb_plot_data, "B. Whole Brain")
cat("\n## Plots\n")
print(roi_plot)
print(wb_plot)
```
###########################################
#SI-3: Additional Analyses Related to Language Lateralization
###########################################
## SI-3.A: Lateralization Index and Magnitude by Component (Frontal and Temporal)
### Data Preparation for Lateralization Analysis
```{r lateralization-data-prep}
# Summarize volume data across ROIs by summing effects for each Subject, age, hemi, and front
ds1_vol_sum <- ds1_vol %>%
group_by(Subject, age, hemi, front) %>%
summarise(contrast = sum(effect, na.rm = TRUE), .groups = "drop")
ds2_vol_sum <- ds2_vol %>%
group_by(Subject, age, hemi, front) %>%
summarise(contrast = sum(effect, na.rm = TRUE), .groups = "drop")
# Filter effect data for language conditions
####Wide to Long DS1
# Filter effect data for language conditions
ds1_eff_lang <- ds1_eff %>%
filter(condition == "language")
# Filter for language condition
levels(ds1_eff_lang$ROI) <- gsub("R ", "", levels(ds1_eff_lang$ROI))
levels(ds1_eff_lang$ROI) <- gsub("IFG orb", "IFGorb", levels(ds1_eff_lang$ROI))
levels(ds1_eff_lang$ROI) <- gsub("Post temp", "PostTemp", levels(ds1_eff_lang$ROI))
## Remove num_outliers, which has some NAs
ds1_eff_lang <- ds1_eff_lang %>% select(-num_outliers)
## Make sure there are no single (null) observations
oneobs <- ds1_eff_lang %>%
group_by(Subject) %>%
summarize(n = n()) %>%
filter(n==1) %>%
pull(Subject)
ds1_eff_lang <- ds1_eff_lang %>%
filter(!Subject %in% oneobs)
na_rows <- ds1_eff_lang %>%
filter(if_any(everything(), is.na))
na_rows
## Inspect data structure
table(ds1_eff_lang$ROI, ds1_eff_lang$hemi)
####Wide to Long DS2
#Extract hemisphere info first, then clean ROI names
ds2_eff_lang <- ds2_eff %>%
dplyr::filter(condition == "EffectSize.mentalsocialphysical")
ds2_eff_lang <- ds2_eff_lang %>%
mutate(
# Extract hemisphere from ROI names
hemi_new = case_when(
stringr::str_detect(as.character(ROI), "^LH_") ~ "lh",
stringr::str_detect(as.character(ROI), "^RH_") ~ "rh",
TRUE ~ as.character(hemi)
),
# Clean ROI names by removing prefixes
ROI_new = str_replace(as.character(ROI), "^(LH_|RH_)", "")
) %>%
# Replace columns
select(-hemi, -ROI) %>%
rename(hemi = hemi_new, ROI = ROI_new) %>%
mutate(
hemi = as.factor(hemi),
ROI = as.factor(ROI)
)
# Now apply the same standardization as DS1
levels(ds2_eff_lang$ROI) <- gsub("R ", "", levels(ds2_eff_lang$ROI))
levels(ds2_eff_lang$ROI) <- gsub("IFG orb", "IFGorb", levels(ds2_eff_lang$ROI))
levels(ds2_eff_lang$ROI) <- gsub("Post temp", "PostTemp", levels(ds2_eff_lang$ROI))
# Check the results
print("DS2 ROI levels after cleaning:")
print(levels(ds2_eff_lang$ROI))
print("DS2 hemisphere levels:")
print(levels(ds2_eff_lang$hemi))
## Inspect data structure
table(ds2_eff_lang$ROI, ds2_eff_lang$hemi)
cat("Data preparation completed for lateralization analysis.\n")
cat("DS1 volume summary dimensions:", dim(ds1_vol_sum), "\n")
cat("DS2 volume summary dimensions:", dim(ds2_vol_sum), "\n")
```
### LI Analysis Functions
```{r lateralization-functions}
# Function to analyze Lateralization Index (LI) for a specific region
analyze_LI_component <- function(data, dataset_label = "Dataset", region = "front") {
# Filter for the specified region (front or post)
data_region <- data %>%
filter(front == region)
# Pivot to wide format: separate 'lh' and 'rh' columns, then calculate LI
data_wide <- data_region %>%
pivot_wider(names_from = hemi, values_from = contrast) %>%
mutate(LI = (lh - rh) / (lh + rh))
# Fit linear model: LI predicted by age
model <- lm(LI ~ age, data = data_wide)
# Estimated marginal means and contrasts (each age vs. adult)
emm <- emmeans(model, ~ age)
contrast_res <- contrast(emm, method = "trt.vs.ctrl", ref = "adult", adjust = "sidak")
# Return results
list(
dataset = dataset_label,
region = region,
model = model,
data_wide = data_wide,
contrasts = contrast_res
)
}
# Function to analyze magnitude (contrast values) for a region
analyze_magnitude_component <- function(data, dataset_label = "Dataset", region = "front") {
# Filter for the specified region (front or post)
data_region <- data %>%
filter(front == region)
# Fit linear mixed model: contrast ~ hemi * age + (hemi | Subject) + (hemi | ROI)
m <- lmer(contrast ~ hemi * age + (hemi | Subject) + (hemi | ROI), data = data_region)
# Estimated marginal means and contrasts for interaction
emm <- emmeans(m, ~ hemi:age, lmer.df = "kenward-roger")
contrast_res <- contrast(emm, interaction = c("pairwise", "trt.vs.ctrl"), adjust = "sidak")
# Return results
list(
dataset = dataset_label,
region = region,
model = m,
contrasts = contrast_res
)
}
cat("Lateralization analysis functions defined.\n")
```
### Lateralization Index Analysis - Dataset 1
```{r li-analysis-ds1}
# Dataset 1 - Frontal regions
cat("Dataset 1 - Frontal LI Analysis:\n")
cat("=================================================\n")
res_li_ds1_front <- analyze_LI_component(ds1_vol_sum, "DS1", region = "front")
print(summary(res_li_ds1_front$model))
cat("\nContrasts (vs. Adult):\n")
print(res_li_ds1_front$contrasts)
cat("\n\nDataset 1 - Temporal LI Analysis:\n")
cat("==================================================\n")
res_li_ds1_post <- analyze_LI_component(ds1_vol_sum, "DS1", region = "post")
print(summary(res_li_ds1_post$model))
cat("\nContrasts (vs. Adult):\n")
print(res_li_ds1_post$contrasts)
```
### LI Analysis - Dataset 2
```{r li-analysis-ds2}
# Dataset 2 - Frontal regions
cat("Dataset 2 - Frontal Lateralization Index Analysis:\n")
cat("=================================================\n")
res_li_ds2_front <- analyze_LI_component(ds2_vol_sum, "DS2", region = "front")
print(summary(res_li_ds2_front$model))
cat("\nContrasts (vs. Adult):\n")
print(res_li_ds2_front$contrasts)
cat("\n\nDataset 2 - Temporal Lateralization Index Analysis:\n")
cat("==================================================\n")
res_li_ds2_post <- analyze_LI_component(ds2_vol_sum, "DS2", region = "post")
print(summary(res_li_ds2_post$model))
cat("\nContrasts (vs. Adult):\n")
print(res_li_ds2_post$contrasts)
```
### Magnitude Analysis - Dataset 1
```{r magnitude-analysis-ds1}
# Dataset 1 - Frontal regions magnitude
cat("Dataset 1 - Frontal Magnitude Analysis:\n")
cat("=====================================\n")
res_mag_ds1_front <- analyze_magnitude_component(ds1_eff_lang, "DS1", region = "front")
print(summary(res_mag_ds1_front$model))
cat("\nContrasts:\n")
print(res_mag_ds1_front$contrasts)
cat("\n\nDataset 1 - Temporal Magnitude Analysis:\n")
cat("======================================\n")
res_mag_ds1_post <- analyze_magnitude_component(ds1_eff_lang, "DS1", region = "post")
print(summary(res_mag_ds1_post$model))
cat("\nContrasts:\n")
print(res_mag_ds1_post$contrasts)
```
### Magnitude Analysis - Dataset 2
```{r magnitude-analysis-ds2}
# Dataset 2 - Frontal regions magnitude
cat("Dataset 2 - Frontal Magnitude Analysis:\n")
cat("=====================================\n")
res_mag_ds2_front <- analyze_magnitude_component(ds2_eff_lang, "DS2", region = "front")
print(summary(res_mag_ds2_front$model))
cat("\nContrasts:\n")
print(res_mag_ds2_front$contrasts)
cat("\n\nDataset 2 - Temporal Magnitude Analysis:\n")
cat("======================================\n")
res_mag_ds2_post <- analyze_magnitude_component(ds2_eff_lang, "DS2", region = "post")
print(summary(res_mag_ds2_post$model))
cat("\nContrasts:\n")
print(res_mag_ds2_post$contrasts)
```
# SI-3B: Bayesian Analysis (commented because long runtime)
```{r Bayes-Factors}
ds1_eff_lang
ds2_eff_lang
## Prepare data for monotonic function
ds1_eff_lang <- ds1_eff_lang %>%
mutate(
age_ordered_factor = factor(age, levels = c("early", "middle", "late", "adult"), ordered = TRUE)
)
ds2_eff_lang <- ds2_eff_lang %>%
mutate(
age_ordered_factor = factor(age, levels = c("early", "middle", "late", "adult"), ordered = TRUE)
)
save(
ds1_eff_lang, ds2_eff_lang,
file="data_to_openmind.RData"
)
# dir.create("brms_models", showWarnings = FALSE)
#
# base_args <- list(
# chains = 4,
# cores = 4,
# iter = 20000,
# warmup = 4000,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.999),
# seed = 123
# )
#
# fit_or_load <- function(name, formula, data, priors, models_dir = "brms_models") {
# f <- file.path(models_dir, paste0(name, ".RDS"))
# if (file.exists(f)) return(readRDS(f))
# fit <- do.call(brm, c(list(formula = formula, data = data, prior = priors), base_args))
# saveRDS(fit, f)
# fit
# }
#
#
# ## ---- Model Fitting ----
# ### ---- DS1 ----
# #### ---- Monotonic ----
# ds1_monotonic_null <- fit_or_load(
# formula = contrast ~ hemi + mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_monotonic_null"
# )
#
# ds1_monotonic_liberal <- fit_or_load(
# formula = contrast ~ hemi * mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, .5), coef="moage_ordered_factor:hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_monotonic_liberal"
# )
#
# ds1_monotonic_moderate <- fit_or_load(
# formula = contrast ~ hemi * mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, .2), coef="moage_ordered_factor:hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_monotonic_moderate"
# )
#
# ds1_monotonic_conservative <- fit_or_load(
# formula = contrast ~ hemi * mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, .1), coef="moage_ordered_factor:hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_monotonic_conservative"
# )
#
# ds1_monotonic_extraconservative <- fit_or_load(
# formula = contrast ~ hemi * mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, .05), coef="moage_ordered_factor:hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_monotonic_extraconservative"
# )
#
# ### ---- Categorical ----
# ds1_categorical_null <- fit_or_load(
# formula = contrast ~ hemi + age +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_categorical_null"
# )
#
# ds1_categorical_liberal <- fit_or_load(
# formula = contrast ~ hemi * age +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 0.5), class = "b"),
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_categorical_liberal"
# )
#
# ds1_categorical_moderate <- fit_or_load(
# formula = contrast ~ hemi * age +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 0.2), class = "b"),
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_categorical_moderate"
# )
#
# ds1_categorical_conservative <- fit_or_load(
# formula = contrast ~ hemi * age +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 0.1), class = "b"),
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_categorical_conservative"
# )
#
# ds1_categorical_extraconservative <- fit_or_load(
# formula = contrast ~ hemi * age +
# (hemi | Subject) + (hemi | ROI),
# data=ds1_eff_lang,
# priors=c(
# prior(normal(0, 0.05), class = "b"),
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds1_categorical_extraconservative"
# )
#
# ### ---- DS2 ----
# #### ---- Monotonic ----
# ds2_monotonic_null <- fit_or_load(
# formula = contrast ~ hemi + mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_monotonic_null"
# )
#
# ds2_monotonic_liberal <- fit_or_load(
# formula = contrast ~ hemi * mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, .5), coef="moage_ordered_factor:hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_monotonic_liberal"
# )
#
# ds2_monotonic_moderate <- fit_or_load(
# formula = contrast ~ hemi * mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, .2), coef="moage_ordered_factor:hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_monotonic_moderate"
# )
#
# ds2_monotonic_conservative <- fit_or_load(
# formula = contrast ~ hemi * mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, .1), coef="moage_ordered_factor:hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_monotonic_conservative"
# )
#
# ds2_monotonic_extraconservative <- fit_or_load(
# formula = contrast ~ hemi * mo(age_ordered_factor) +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 1), class="b"),
# prior(normal(0, .05), coef="moage_ordered_factor:hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_monotonic_extraconservative"
# )
#
# ### ---- Categorical ----
# ds2_categorical_null <- fit_or_load(
# formula = contrast ~ hemi + age +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_categorical_null"
# )
#
# ds2_categorical_liberal <- fit_or_load(
# formula = contrast ~ hemi * age +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 0.5), class = "b"),
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_categorical_liberal"
# )
#
# ds2_categorical_moderate <- fit_or_load(
# formula = contrast ~ hemi * age +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 0.2), class = "b"),
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_categorical_moderate"
# )
#
# ds2_categorical_conservative <- fit_or_load(
# formula = contrast ~ hemi * age +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 0.1), class = "b"),
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_categorical_conservative"
# )
#
# ds2_categorical_extraconservative <- fit_or_load(
# formula = contrast ~ hemi * age +
# (hemi | Subject) + (hemi | ROI),
# data=ds2_eff_lang,
# priors=c(
# prior(normal(0, 0.05), class = "b"),
# prior(normal(0, 1), coef="ageearly"),
# prior(normal(0, 1), coef="agemiddle"),
# prior(normal(0, 1), coef="agelate"),
# prior(normal(0, 1), coef="hemirh"),
# prior(normal(0, 2), class="Intercept"),
# prior(normal(0, 2), class="sigma")
# ),
# name = "ds2_categorical_extraconservative"
# )
#
# ## ---- Bayes Factors ----
# ### ---- DS1 ----
# if (file.exists("brms_models/ds1_bayes_factors.RData")) {
# load("brms_models/ds1_bayes_factors.RData")
# } else {
# BF01_monotonic_liberal_ds1 <- bayes_factor(ds1_monotonic_null, ds1_monotonic_liberal)
# BF01_monotonic_moderate_ds1 <- bayes_factor(ds1_monotonic_null, ds1_monotonic_moderate)
# BF01_monotonic_conservative_ds1 <- bayes_factor(ds1_monotonic_null, ds1_monotonic_conservative)
# BF01_monotonic_extraconservative_ds1 <- bayes_factor(ds1_monotonic_null, ds1_monotonic_extraconservative)
#
# BF01_categorical_liberal_ds1 <- bayes_factor(ds1_categorical_null, ds1_categorical_liberal)
# BF01_categorical_moderate_ds1 <- bayes_factor(ds1_categorical_null, ds1_categorical_moderate)
# BF01_categorical_conservative_ds1 <- bayes_factor(ds1_categorical_null, ds1_categorical_conservative)
# BF01_categorical_extraconservative_ds1 <- bayes_factor(ds1_categorical_null, ds1_categorical_extraconservative)
#
# save(
# BF01_monotonic_liberal_ds1, BF01_monotonic_moderate_ds1,
# BF01_monotonic_conservative_ds1, BF01_monotonic_extraconservative_ds1,
# BF01_categorical_liberal_ds1, BF01_categorical_moderate_ds1,
# BF01_categorical_conservative_ds1, BF01_categorical_extraconservative_ds1,
# file="brms_models/ds1_bayes_factors.RData"
# )
#
# }
#
# ### ---- DS2 ----
# if (file.exists("brms_models/ds2_bayes_factors.RData")) {
# load("brms_models/ds2_bayes_factors.RData")
# } else {
# BF01_monotonic_liberal_ds2 <- bayes_factor(ds2_monotonic_null, ds2_monotonic_liberal)
# BF01_monotonic_moderate_ds2 <- bayes_factor(ds2_monotonic_null, ds2_monotonic_moderate)
# BF01_monotonic_conservative_ds2 <- bayes_factor(ds2_monotonic_null, ds2_monotonic_conservative)
# BF01_monotonic_extraconservative_ds2 <- bayes_factor(ds2_monotonic_null, ds2_monotonic_extraconservative)
#
# BF01_categorical_liberal_ds2 <- bayes_factor(ds2_categorical_null, ds2_categorical_liberal)
# BF01_categorical_moderate_ds2 <- bayes_factor(ds2_categorical_null, ds2_categorical_moderate)
# BF01_categorical_conservative_ds2 <- bayes_factor(ds2_categorical_null, ds2_categorical_conservative)
# BF01_categorical_extraconservative_ds2 <- bayes_factor(ds2_categorical_null, ds2_categorical_extraconservative)
#
# save(
# BF01_monotonic_liberal_ds2, BF01_monotonic_moderate_ds2,
# BF01_monotonic_conservative_ds2, BF01_monotonic_extraconservative_ds2,
# BF01_categorical_liberal_ds2, BF01_categorical_moderate_ds2,
# BF01_categorical_conservative_ds2, BF01_categorical_extraconservative_ds2,
# file="brms_models/ds2_bayes_factors.RData"
# )
#
# }
#
# ## Bayes Factors with monotonic age function
# BF01_monotonic_liberal_ds1$bf
# BF01_monotonic_moderate_ds1$bf
# BF01_monotonic_conservative_ds1$bf
# BF01_monotonic_extraconservative_ds1$bf
#
# 1/BF01_monotonic_liberal_ds1$bf
# 1/BF01_monotonic_moderate_ds1$bf
# 1/BF01_monotonic_conservative_ds1$bf
# 1/BF01_monotonic_extraconservative_ds1$bf
#
#
# BF01_monotonic_liberal_ds2$bf
# BF01_monotonic_moderate_ds2$bf
# BF01_monotonic_conservative_ds2$bf
# BF01_monotonic_extraconservative_ds2$bf
#
# 1/BF01_monotonic_liberal_ds2$bf
# 1/BF01_monotonic_moderate_ds2$bf
# 1/BF01_monotonic_conservative_ds2$bf
# 1/BF01_monotonic_extraconservative_ds2$bf
#
#
# ## Bayes factors with categorical age
# BF01_categorical_liberal_ds1$bf
# BF01_categorical_moderate_ds1$bf
# BF01_categorical_conservative_ds1$bf
# BF01_categorical_extraconservative_ds1$bf
#
# 1/BF01_categorical_liberal_ds1$bf
# 1/BF01_categorical_moderate_ds1$bf
# 1/BF01_categorical_conservative_ds1$bf
# 1/BF01_categorical_extraconservative_ds1$bf
#
# BF01_categorical_liberal_ds2$bf
# BF01_categorical_moderate_ds2$bf
# BF01_categorical_conservative_ds2$bf
# BF01_categorical_extraconservative_ds2$bf
#
# 1/BF01_categorical_liberal_ds2$bf
# 1/BF01_categorical_moderate_ds2$bf
# 1/BF01_categorical_conservative_ds2$bf
# 1/BF01_categorical_extraconservative_ds2$bf
#
# ## Select model outputs
# ds1_monotonic_extraconservative
# ds2_monotonic_conservative
# ds1_categorical_moderate
# ds1_categorical_conservative
#
#
# ## REMOVE IF NOT USING :: PP_CHECK
# lmerfit <- lmer(contrast ~ hemi * age +
# (hemi | Subject) + (hemi | ROI), data=ds1_eff_lang, control=lmerControl(optimizer = "bobyqa"))
#
# ysim <- tibble(
# sample = numeric(),
# yhat = numeric()
# )
#
# for(i in 1:100) {
# ysim <- ysim %>%
# add_row(
# yhat=unlist(simulate(lmerfit, 1)),
# sample=i
# )
# }
#
# ggplot() +
# geom_density(data=ysim,
# aes(x=yhat, group=sample),
# color="cadetblue",
# linewidth=.1) +
# geom_density(aes(x=ds1_eff_lang$contrast))
#
# ggsave("posterior_predictive_check.png")
```
# SI-3.C: Lateralization Index Analysis Using LI Toolbox
### Load and Prepare LI Toolbox Data
### Load and Prepare LI Toolbox Data
```{r li-toolbox-data-prep}
# Load LI toolbox output data for both datasets
ds1_li_toolbox <- read.csv("../data/ds1_li_toolbox_output.csv")
ds2_li_toolbox <- read.csv("../data/ds2_li_toolbox_output.csv")
# Set age factor levels for both datasets
ds1_li_toolbox$age <- factor(ds1_li_toolbox$Group, levels = c("early", "middle", "late", "adult"))
ds2_li_toolbox$age <- factor(ds2_li_toolbox$Group, levels = c("early", "middle", "late", "adult"))
cat("LI Toolbox data loaded successfully.\n")
cat("DS1 LI data points:", nrow(ds1_li_toolbox), "\n")
cat("DS2 LI data points:", nrow(ds2_li_toolbox), "\n")
```
```{r li-toolbox-analysis}
# SI-3.C: Lateralization Index Analysis Using LI Toolbox
## Load and Prepare LI Toolbox Data
# Standardize Subject ID columns function
standardize_subject_column <- function(data) {
if("Subject" %in% names(data)) {
return(data)
} else if("Participant" %in% names(data)) {
return(data %>% rename(Subject = Participant))
} else if("participant" %in% names(data)) {
return(data %>% rename(Subject = participant))
} else {
# Assume first column is Subject ID
names(data)[1] <- "Subject"
return(data)
}
}
# Apply standardization to LI datasets
ds1_li_toolbox <- standardize_subject_column(ds1_li_toolbox)
ds2_li_toolbox <- standardize_subject_column(ds2_li_toolbox)
# Ensure consistent factor levels for age groups
age_levels <- c("early", "middle", "late", "adult")
ds1_li_toolbox$age <- factor(ds1_li_toolbox$age, levels = age_levels)
ds2_li_toolbox$age <- factor(ds2_li_toolbox$age, levels = age_levels)
## Visualization Function
create_li_barplot <- function(data, title) {
ggbarplot(data,
x = "age",
y = "LI_overall",
add = c("mean_se", "jitter"),
add.params = list(size = 1, alpha = 0.2),
width = 0.6,
size = 0.5,
fill = "age",
palette = c("gray", "gray", "gray", "black"),
alpha = 0.6,
title = title,
ylab = "Lateralization Index (LI)",
xlab = "Age Group") +
geom_hline(yintercept = 0, linetype = "solid", color = "black") +
theme_classic() +
theme(legend.position = "none",
text = element_text(size = 18),
plot.title = element_text(hjust = 0.5, face = "bold")) +
scale_y_continuous(limits = c(-1, NA))
}
## Statistical Analysis Function
analyze_overall_li <- function(data, dataset_name) {
# Fit model with adult as reference
model_li <- lm(LI_overall ~ age, data = data)
# Get emmeans and contrasts vs adult
emm_result <- emmeans(model_li, ~ age)
contrasts_result <- contrast(emm_result, method = "trt.vs.ctrl", ref = "adult", adjust = "sidak")
# ANOVA for overall effect
anova_results <- anova(model_li)
eta_squared <- anova_results$"Sum Sq"[1] / sum(anova_results$"Sum Sq")
# Output results
cat("\n", dataset_name, " Statistical Analysis:\n")
cat("==============================\n")
cat("ANOVA Results:\n")
print(anova_results)
cat("\nPost-hoc Comparisons (vs. Adult):\n")
print(contrasts_result)
cat("\nEffect size (eta-squared):", round(eta_squared, 4), "\n\n")
return(list(
model = model_li,
contrasts = contrasts_result,
eta_squared = eta_squared,
anova = anova_results
))
}
## Table Generation Function
create_li_contrast_table <- function(contrast_result, title) {
contrast_df <- as.data.frame(contrast_result)
# Format results
results <- contrast_df %>%
mutate(
Comparison = gsub(" - adult", " vs. Adult", contrast),
B = round(estimate, 2),
SE = round(SE, 2),
t = round(t.ratio, 2),
p = case_when(
p.value < 0.001 ~ "**<0.001**",
p.value < 0.01 ~ paste0("**", sprintf("%.3f", p.value), "**"),
p.value < 0.05 ~ paste0("**", sprintf("%.3f", p.value), "**"),
TRUE ~ sprintf("%.3f", p.value)
)
) %>%
select(Comparison, B, SE, t, p)
# Create kable table
kable(results,
col.names = c("", "***B***", "***SE***", "***t***", "***p***"),
caption = paste0("**", title, "**"),
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
}
## Execute Dataset 1 Analysis
cat("## Dataset 1 - Lateralization Index Analysis\n")
# Create visualization
mean_lat_ds1 <- create_li_barplot(ds1_li_toolbox, "Dataset 1: Lateralization Index")
print(mean_lat_ds1)
# Run statistical analysis
ds1_results <- analyze_overall_li(ds1_li_toolbox, "Dataset 1")
# Create summary table
ds1_table <- create_li_contrast_table(ds1_results$contrasts,
"Dataset 1 - Lateralization Index Comparisons (vs. Adult)")
print(ds1_table)
## Execute Dataset 2 Analysis
cat("## Dataset 2 - Lateralization Index Analysis\n")
# Create visualization
mean_lat_ds2 <- create_li_barplot(ds2_li_toolbox, "Dataset 2: Lateralization Index")
print(mean_lat_ds2)
# Run statistical analysis
ds2_results <- analyze_overall_li(ds2_li_toolbox, "Dataset 2")
# Create summary table
ds2_table <- create_li_contrast_table(ds2_results$contrasts,
"Dataset 2 - Lateralization Index Comparisons (vs. Adult)")
print(ds2_table)
## Summary Statistics
cat("## Summary Statistics\n")
cat("Dataset 1 eta-squared:", round(ds1_results$eta_squared, 4), "\n")
cat("Dataset 2 eta-squared:", round(ds2_results$eta_squared, 4), "\n")
```
## SI-3.D: Lateralization Analysis Combining Datasets 1 and 2
### Data Preparation for Combined Analysis
```{r combined-lat-data-prep}
# SI-3.D: Lateralization Analysis Combining Datasets 1 and 2
## Data Preparation for Combined Analysis
# Filter combined effect data for language conditions from both datasets
ds_eff_combined_lang <- ds_eff_combined %>%
filter(condition == "language" | condition == "EffectSize.mentalsocialphysical") %>%
mutate(age = factor(age, levels = c("early", "middle", "late", "adult")))
# Prepare combined volume data for LI calculation
ds_vol_sum_combined <- ds_vol_combined %>%
group_by(Subject, age, hemi, front, ds) %>%
dplyr::summarise(vol = sum(effect, na.rm = TRUE), .groups = "drop") %>%
mutate(age = factor(age, levels = c("early", "middle", "late", "adult")))
cat("Combined dataset prepared for lateralization analysis.\n")
cat("Effect data - Total observations:", nrow(ds_eff_combined_lang), "\n")
cat("Effect data - Number of subjects:", length(unique(ds_eff_combined_lang$Subject)), "\n")
cat("Volume data - Total observations:", nrow(ds_vol_sum_combined), "\n")
cat("Volume data - Number of subjects:", length(unique(ds_vol_sum_combined$Subject)), "\n")
cat("\nDataset distribution:\n")
print(table(ds_eff_combined_lang$ds))
cat("\nAge group distribution:\n")
print(table(ds_eff_combined_lang$age))
cat("\nHemisphere distribution:\n")
print(table(ds_eff_combined_lang$hemi))
## 1. MIXED-EFFECTS MODEL ANALYSIS
# Fit mixed-effects model for activation magnitude
m_combined_effect <- lmer(contrast ~ hemi * age + (hemi | Subject) + (hemi | ROI),
data = ds_eff_combined_lang, REML = FALSE)
# Model summary and R-squared
cat("\n## Combined Dataset Mixed-Effects Model\n")
cat("=====================================\n")
print(summary(m_combined_effect))
r2_values <- r.squaredGLMM(m_combined_effect)
cat("\nModel R-squared Values:\n")
cat("Marginal R² (fixed effects):", round(r2_values[1], 4), "\n")
cat("Conditional R² (fixed + random):", round(r2_values[2], 4), "\n")
# Post-hoc contrasts
cat("\nPost-hoc Contrasts:\n")
emm_combined <- emmeans(m_combined_effect, ~ hemi:age, lmer.df = "kenward-roger")
contrasts_combined <- contrast(emm_combined, interaction = c("pairwise", "trt.vs.ctrl"),
by = NULL, adjust = "sidak")
print(contrasts_combined)
## 2. LATERALIZATION INDEX ANALYSIS
# Calculate LI from combined volume data
data_wide_li <- ds_vol_sum_combined %>%
pivot_wider(names_from = hemi, values_from = vol) %>%
mutate(
LI = ifelse(!is.na(lh) & !is.na(rh) & (lh + rh) != 0,
(lh - rh) / (lh + rh), NA_real_)
) %>%
filter(!is.na(LI) & is.finite(LI))
cat("\n## Lateralization Index Analysis\n")
cat("Lateralization Index calculated for", nrow(data_wide_li), "subjects.\n")
# Statistical analysis
model_li <- lm(LI ~ age, data = data_wide_li)
anova_li <- anova(model_li)
eta_squared_li <- anova_li$"Sum Sq"[1] / sum(anova_li$"Sum Sq")
cat("\nLI Model Results:\n")
print(summary(model_li))
cat("\nANOVA Results:\n")
print(anova_li)
cat("\nEffect size (eta-squared):", round(eta_squared_li, 4), "\n")
# Contrasts vs. adult
emm_li <- emmeans(model_li, ~ age)
contrasts_li <- contrast(emm_li, method = "trt.vs.ctrl", ref = "adult", adjust = "sidak")
cat("\nContrasts vs. Adult:\n")
print(contrasts_li)
## 3. DATA PREPARATION FOR VISUALIZATION
# LI plot data - one point per participant (average across regions/datasets if multiple)
combined_li_plot_data <- data_wide_li %>%
group_by(Subject, age) %>%
summarise(mean_LI = mean(LI, na.rm = TRUE), .groups = "drop") %>%
filter(!is.na(mean_LI) & is.finite(mean_LI))
# Magnitude difference data - one point per participant
combined_mag_data <- ds_eff_combined_lang %>%
group_by(Subject, age, hemi) %>%
summarise(contrast = mean(contrast, na.rm = TRUE), .groups = "drop") %>%
pivot_wider(names_from = hemi, values_from = contrast) %>%
mutate(difference = lh - rh) %>%
filter(!is.na(difference) & is.finite(difference)) %>%
group_by(Subject, age) %>%
summarise(difference = mean(difference, na.rm = TRUE), .groups = "drop")
cat("\nVisualization data prepared:\n")
cat("LI plot data:", nrow(combined_li_plot_data), "subjects\n")
cat("Magnitude plot data:", nrow(combined_mag_data), "subjects\n")
## 4. STATISTICAL ANALYSIS FUNCTIONS
analyze_combined_metric <- function(data, metric_col, title) {
# Fit model
formula_str <- paste(metric_col, "~ age")
model <- lm(as.formula(formula_str), data = data)
# ANOVA and effect size
anova_result <- anova(model)
eta_squared <- anova_result$"Sum Sq"[1] / sum(anova_result$"Sum Sq")
# Emmeans and contrasts
emm_result <- emmeans(model, ~ age)
contrasts_result <- contrast(emm_result, method = "trt.vs.ctrl", ref = "adult", adjust = "sidak")
# Output results
cat("\n", title, " Analysis:\n")
cat("========================\n")
print(summary(model))
cat("\nANOVA Results:\n")
print(anova_result)
cat("\nEffect size (eta-squared):", round(eta_squared, 4), "\n")
cat("\nContrasts vs. Adult:\n")
print(contrasts_result)
return(list(
model = model,
anova = anova_result,
contrasts = contrasts_result,
eta_squared = eta_squared
))
}
## 5. VISUALIZATION FUNCTIONS
create_combined_plot <- function(data, y_var, y_label, y_limits) {
ggplot(data, aes_string(x = "age", y = y_var, fill = "age")) +
geom_hline(yintercept = 0, linetype = "solid", color = "black", linewidth = 1) +
geom_bar(stat = "summary", fun = "mean", width = 0.6, alpha = 0.6,
color = "black", linewidth = 0.5) +
geom_jitter(width = 0.2, alpha = 0.2, size = 1) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2, linewidth = 1) +
scale_fill_manual(values = c("gray", "gray", "gray", "black")) +
scale_y_continuous(limits = y_limits, expand = expansion(mult = c(0.02, 0.02))) +
labs(y = y_label, x = "Age Group") +
theme_classic() +
theme(
legend.position = "none",
text = element_text(size = 18),
axis.text = element_text(size = 14),
axis.title = element_text(size = 16),
plot.margin = margin(10, 10, 20, 10)
)
}
## 6. CREATE FIGURES
# Panel A: Lateralization Index
fig_li <- create_combined_plot(
data = combined_li_plot_data,
y_var = "mean_LI",
y_label = "Lateralization Index\n(LH - RH) / (LH + RH)",
y_limits = c(-1, 1)
)
# Panel B: Magnitude Difference
fig_magnitude <- create_combined_plot(
data = combined_mag_data,
y_var = "difference",
y_label = "Magnitude Difference\nLH - RH",
y_limits = c(-4, 4)
)
# Combine figures
combined_figure <- plot_grid(fig_li, fig_magnitude, ncol = 2, align = "hv",
labels = c("A", "B"), label_size = 18)
cat("\n## Combined Figure\n")
print(combined_figure)
## 7. SUMMARY TABLES FUNCTIONS
create_summary_table <- function(data, metric_col, caption) {
summary_stats <- data %>%
group_by(age) %>%
summarise(
N = n(),
Mean = round(mean(.data[[metric_col]], na.rm = TRUE), 4),
SD = round(sd(.data[[metric_col]], na.rm = TRUE), 4),
SE = round(sd(.data[[metric_col]], na.rm = TRUE) / sqrt(n()), 4),
.groups = "drop"
)
kable(summary_stats,
caption = caption,
col.names = c("Age", "N", "Mean", "SD", "SE")) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
}
create_contrast_table <- function(contrast_result, title) {
contrast_df <- as.data.frame(contrast_result) %>%
mutate(
Comparison = gsub(" - adult", " vs. Adult", contrast),
B = round(estimate, 2),
SE = round(SE, 2),
t = round(t.ratio, 2),
p = case_when(
p.value < 0.001 ~ "**<0.001**",
p.value < 0.01 ~ paste0("**", sprintf("%.3f", p.value), "**"),
p.value < 0.05 ~ paste0("**", sprintf("%.3f", p.value), "**"),
TRUE ~ sprintf("%.3f", p.value)
)
) %>%
select(Comparison, B, SE, t, p)
kable(contrast_df,
col.names = c("", "***B***", "***SE***", "***t***", "***p***"),
caption = paste0("**", title, "**"),
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
}
## 8. GENERATE SUMMARY TABLES
cat("\n## Summary Tables\n")
# LI summary statistics
li_summary_table <- create_summary_table(combined_li_plot_data, "mean_LI",
"Combined Dataset - Lateralization Index Summary by Age Group")
print(li_summary_table)
# LI contrasts table
li_contrast_table <- create_contrast_table(contrasts_li,
"Combined Dataset - Lateralization Index Contrasts (vs. Adult)")
print(li_contrast_table)
# Magnitude summary statistics
mag_summary_table <- create_summary_table(combined_mag_data, "difference",
"Combined Dataset - Magnitude Difference Summary by Age Group")
print(mag_summary_table)
# Analyze magnitude differences
mag_results <- analyze_combined_metric(combined_mag_data, "difference", "Magnitude Difference")
# Magnitude contrasts table
mag_contrast_table <- create_contrast_table(mag_results$contrasts,
"Combined Dataset - Magnitude Difference Contrasts (vs. Adult)")
print(mag_contrast_table)
## 9. ANALYSIS SUMMARY
cat("\n## Analysis Summary\n")
cat("==================\n")
cat("Mixed-effects model R² (marginal):", round(r2_values[1], 4), "\n")
cat("Mixed-effects model R² (conditional):", round(r2_values[2], 4), "\n")
cat("LI analysis effect size (eta²):", round(eta_squared_li, 4), "\n")
cat("Magnitude analysis effect size (eta²):", round(mag_results$eta_squared, 4), "\n")
cat("LI sample size:", nrow(combined_li_plot_data), "subjects\n")
cat("Magnitude sample size:", nrow(combined_mag_data), "subjects\n")
# Check data distribution by dataset
cat("\nData distribution by dataset:\n")
li_by_ds <- data_wide_li %>%
group_by(ds, age) %>%
summarise(n = n(), .groups = "drop")
print(li_by_ds)
# Store results for further use
combined_lateralization_results <- list(
mixed_effects_model = m_combined_effect,
li_model = model_li,
magnitude_model = mag_results$model,
li_contrasts = contrasts_li,
magnitude_contrasts = mag_results$contrasts,
combined_figure = combined_figure,
li_data = combined_li_plot_data,
magnitude_data = combined_mag_data
)
```
```{r}
# Combined Dataset Analysis: Mixed-Effects Models and Lateralization
## 1. MIXED-EFFECTS MODEL ANALYSIS
# Fit mixed-effects model for activation magnitude
m_combined_effect <- lmer(contrast ~ hemi * age + (hemi | Subject) + (hemi | ROI),
data = ds_eff_combined_lang, REML = FALSE)
# Model summary and R-squared
cat("## Combined Dataset Mixed-Effects Model\n")
cat("=====================================\n")
print(summary(m_combined_effect))
r2_values <- r.squaredGLMM(m_combined_effect)
cat("\nModel R-squared Values:\n")
cat("Marginal R² (fixed effects):", round(r2_values[1], 4), "\n")
cat("Conditional R² (fixed + random):", round(r2_values[2], 4), "\n")
# Post-hoc contrasts
cat("\nPost-hoc Contrasts:\n")
emm_combined <- emmeans(m_combined_effect, ~ hemi:age, lmer.df = "kenward-roger")
contrasts_combined <- contrast(emm_combined, interaction = c("pairwise", "trt.vs.ctrl"),
by = NULL, adjust = "sidak")
print(contrasts_combined)
## 2. LATERALIZATION INDEX ANALYSIS
# Calculate LI from volume data
data_wide_li <- ds_vol_sum_combined %>%
pivot_wider(names_from = hemi, values_from = vol) %>%
mutate(
LI = ifelse(!is.na(lh) & !is.na(rh) & (lh + rh) != 0,
(lh - rh) / (lh + rh), NA_real_)
) %>%
filter(!is.na(LI))
# Ensure factor levels
data_wide_li$age <- factor(data_wide_li$age, levels = c("early", "middle", "late", "adult"))
cat("\n## Lateralization Index Analysis\n")
cat("Lateralization Index calculated for", nrow(data_wide_li), "subjects.\n")
# Statistical analysis
model_li <- lm(LI ~ age, data = data_wide_li)
anova_li <- anova(model_li)
eta_squared_li <- anova_li$"Sum Sq"[1] / sum(anova_li$"Sum Sq")
cat("\nLI Model Results:\n")
print(summary(model_li))
cat("\nANOVA Results:\n")
print(anova_li)
cat("\nEffect size (eta-squared):", round(eta_squared_li, 4), "\n")
# Contrasts vs. adult
emm_li <- emmeans(model_li, ~ age)
contrasts_li <- contrast(emm_li, method = "trt.vs.ctrl", ref = "adult", adjust = "sidak")
cat("\nContrasts vs. Adult:\n")
print(contrasts_li)
## 3. DATA PREPARATION FOR VISUALIZATION
# LI plot data - one point per participant
combined_li_plot_data <- data_wide_li %>%
group_by(Subject, age) %>%
summarise(mean_LI = mean(LI, na.rm = TRUE), .groups = "drop") %>%
na.omit()
# Magnitude difference data - one point per participant
combined_mag_data <- ds_eff_combined_lang %>%
group_by(Subject, age, hemi) %>%
summarise(contrast = mean(contrast, na.rm = TRUE), .groups = "drop") %>%
pivot_wider(names_from = hemi, values_from = contrast) %>%
mutate(difference = lh - rh) %>%
na.omit() %>%
group_by(Subject, age) %>%
summarise(difference = mean(difference, na.rm = TRUE), .groups = "drop")
# Ensure consistent factor levels
combined_li_plot_data$age <- factor(combined_li_plot_data$age, levels = c("early", "middle", "late", "adult"))
combined_mag_data$age <- factor(combined_mag_data$age, levels = c("early", "middle", "late", "adult"))
## 4. VISUALIZATION FUNCTION
create_aligned_plot <- function(data, y_var, y_label, y_limits) {
ggplot(data, aes_string(x = "age", y = y_var, fill = "age")) +
geom_hline(yintercept = 0, linetype = "solid", color = "black", linewidth = 1) +
geom_bar(stat = "summary", fun = "mean", width = 0.6, alpha = 0.6,
color = "black", linewidth = 0.5) +
geom_jitter(width = 0.2, alpha = 0.2, size = 1) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2, linewidth = 1) +
scale_fill_manual(values = c("gray", "gray", "gray", "black")) +
scale_y_continuous(limits = y_limits, expand = expansion(mult = c(0.02, 0.02))) +
labs(y = y_label, x = "Age Group") +
theme_classic() +
theme(
legend.position = "none",
text = element_text(size = 18),
axis.text = element_text(size = 14),
axis.title = element_text(size = 16),
plot.margin = margin(10, 10, 20, 10)
)
}
## 5. CREATE FIGURES
# Panel A: Lateralization Index
fig_li <- create_aligned_plot(
data = combined_li_plot_data,
y_var = "mean_LI",
y_label = "Lateralization Index\n(LH - RH) / (LH + RH)",
y_limits = c(-1, 1)
)
# Panel B: Magnitude Difference
fig_magnitude <- create_aligned_plot(
data = combined_mag_data,
y_var = "difference",
y_label = "Magnitude Difference\nLH - RH",
y_limits = c(-4, 4)
)
# Combine figures
combined_figure <- plot_grid(fig_li, fig_magnitude, ncol = 2, align = "hv",
labels = c("A", "B"), label_size = 18)
cat("\n## Combined Figure\n")
print(combined_figure)
## 6. SUMMARY TABLES
# LI summary statistics
li_summary <- data_wide_li %>%
group_by(age) %>%
summarise(
N = n(),
Mean_LI = round(mean(LI, na.rm = TRUE), 4),
SD_LI = round(sd(LI, na.rm = TRUE), 4),
SE_LI = round(sd(LI, na.rm = TRUE) / sqrt(n()), 4),
.groups = "drop"
)
# LI contrasts table
li_contrast_df <- as.data.frame(contrasts_li) %>%
mutate(
Comparison = gsub(" - adult", " vs. Adult", contrast),
B = round(estimate, 2),
SE = round(SE, 2),
t = round(t.ratio, 2),
p = case_when(
p.value < 0.001 ~ "**<0.001**",
p.value < 0.01 ~ paste0("**", sprintf("%.3f", p.value), "**"),
p.value < 0.05 ~ paste0("**", sprintf("%.3f", p.value), "**"),
TRUE ~ sprintf("%.3f", p.value)
)
) %>%
select(Comparison, B, SE, t, p)
# Display tables
cat("\n## Summary Tables\n")
print(kable(li_summary, caption = "Lateralization Index Summary by Age Group",
col.names = c("Age", "N", "Mean LI", "SD", "SE")) %>%
kable_styling(bootstrap_options = c("striped", "hover")))
print(kable(li_contrast_df, caption = "Lateralization Index Contrasts (vs. Adult)",
col.names = c("", "***B***", "***SE***", "***t***", "***p***"),
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")))
```
#SI-3E. Cross-hemispheric functional correlations
```{r}
# Cross-Hemispheric Functional Correlations Analysis
## Load required libraries
library(ggsignif)
## Data Preparation and Standardization
# Load cross-hemispheric correlation data
rs_w <- read.csv("../data/ds1_rs_w.csv")
# Standardize Subject column name
if("Subject" %in% names(rs_w)) {
# Already standardized
} else if("Participant" %in% names(rs_w)) {
rs_w <- rs_w %>% rename(Subject = Participant)
} else {
names(rs_w)[1] <- "Subject"
}
# Compute correlations between each LH-RH pair for each subject
roi_cor_df <- rs_w %>%
group_by(Subject) %>%
summarize(
IFGorb = cor(value.Lang_LH_IFGorb, value.Lang_RH_IFGorb, use = "pairwise.complete.obs"),
IFG = cor(value.Lang_LH_IFG, value.Lang_RH_IFG, use = "pairwise.complete.obs"),
MFG = cor(value.Lang_LH_MFG, value.Lang_RH_MFG, use = "pairwise.complete.obs"),
AntTemp = cor(value.Lang_LH_AntTemp, value.Lang_RH_AntTemp, use = "pairwise.complete.obs"),
PostTemp = cor(value.Lang_LH_PostTemp, value.Lang_RH_PostTemp, use = "pairwise.complete.obs"),
.groups = "drop"
)
# Compute mean correlation across ROIs per participant
roi_cor_summary <- roi_cor_df %>%
rowwise() %>%
mutate(mean_rs = mean(c(IFGorb, IFG, MFG, AntTemp, PostTemp), na.rm = TRUE)) %>%
select(Subject, mean_rs) %>%
ungroup()
# Add age group with consistent factor levels
roi_cor_summary$age <- as.factor(
ifelse(grepl("inside", roi_cor_summary$Subject), "late",
ifelse(grepl("reader", roi_cor_summary$Subject), "middle",
ifelse(grepl("FACT", roi_cor_summary$Subject), "early", "adult")))
)
roi_cor_summary$age <- factor(roi_cor_summary$age, levels = c("early", "middle", "late", "adult"))
cat("Cross-hemispheric correlations calculated for", nrow(roi_cor_summary), "subjects.\n")
## Statistical Analysis Function
analyze_cross_hemi_corr <- function(data, title) {
# Fit model with adult as reference
model_rs <- lm(mean_rs ~ age, data = data)
# ANOVA and effect size
anova_rs <- anova(model_rs)
eta_squared <- anova_rs$"Sum Sq"[1] / sum(anova_rs$"Sum Sq")
# Emmeans and contrasts vs adult
emm_result <- emmeans(model_rs, ~ age)
contrasts_result <- contrast(emm_result, method = "trt.vs.ctrl", ref = "adult", adjust = "sidak")
# Output results
cat("\n", title, " Statistical Analysis:\n")
cat("=====================================\n")
cat("Model Summary:\n")
print(summary(model_rs))
cat("\nANOVA Results:\n")
print(anova_rs)
cat("\nEffect size (eta-squared):", round(eta_squared, 4), "\n")
cat("\nContrasts vs. Adult Group:\n")
print(contrasts_result)
return(list(
model = model_rs,
anova = anova_rs,
contrasts = contrasts_result,
eta_squared = eta_squared
))
}
## Summary Statistics Function
create_rs_summary_table <- function(data) {
rs_summary <- data %>%
group_by(age) %>%
summarise(
N = n(),
Mean_RS = round(mean(mean_rs, na.rm = TRUE), 4),
SD_RS = round(sd(mean_rs, na.rm = TRUE), 4),
SE_RS = round(sd(mean_rs, na.rm = TRUE) / sqrt(n()), 4),
.groups = "drop"
)
kable(rs_summary,
caption = "Cross-Hemispheric Correlation Summary by Age Group",
col.names = c("Age", "N", "Mean RS", "SD", "SE")) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
}
## Contrast Results Table Function
create_rs_contrast_table <- function(contrast_result, title) {
contrast_df <- as.data.frame(contrast_result)
results <- contrast_df %>%
mutate(
Comparison = gsub(" - adult", " vs. Adult", contrast),
B = round(estimate, 2),
SE = round(SE, 2),
t = round(t.ratio, 2),
p = case_when(
p.value < 0.001 ~ "**<0.001**",
p.value < 0.01 ~ paste0("**", sprintf("%.3f", p.value), "**"),
p.value < 0.05 ~ paste0("**", sprintf("%.3f", p.value), "**"),
TRUE ~ sprintf("%.3f", p.value)
)
) %>%
select(Comparison, B, SE, t, p)
kable(results,
col.names = c("", "***B***", "***SE***", "***t***", "***p***"),
caption = paste0("**", title, "**"),
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
}
## Execute Analysis
cat("## Cross-Hemispheric Functional Correlations Analysis\n")
# Summary statistics table
rs_summary_table <- create_rs_summary_table(roi_cor_summary)
print(rs_summary_table)
# Statistical analysis
rs_results <- analyze_cross_hemi_corr(roi_cor_summary, "Cross-Hemispheric Correlations")
# Create figure (keeping original style as requested)
cross_hemi_corr_plot <- ggbarplot(roi_cor_summary,
x = "age",
y = "mean_rs",
add = c("mean_se", "jitter"),
add.params = list(size = 1, alpha = 0.2),
width = 0.6,
size = 0.5,
fill = "age",
palette = c("gray", "gray", "gray", "black"),
alpha = 0.6,
ylab = "Cross-hemispheric\nCorrelation",
xlab = "Age Group") +
geom_hline(yintercept = 0, linetype = "solid", color = "black") +
geom_signif(comparisons = list(c("middle", "adult")),
annotations = "***",
y_position = 0.9,
tip_length = 0.02,
textsize = 8,
family = "Times New Roman") +
theme_classic() +
theme(legend.position = "none",
text = element_text(size = 30, family = "Times New Roman"),
plot.title = element_text(hjust = 0.5, size = 24, family = "Times New Roman")) +
coord_cartesian(ylim = c(-0.5, 1))
print(cross_hemi_corr_plot)
# Create contrast results table
rs_contrast_table <- create_rs_contrast_table(rs_results$contrasts,
"Cross-Hemispheric Correlations Comparisons (vs. Adult)")
print(rs_contrast_table)
## Analysis Summary
cat("\n## Analysis Summary\n")
cat("==================\n")
cat("Effect size (eta-squared):", round(rs_results$eta_squared, 4), "\n")
cat("Sample size:", nrow(roi_cor_summary), "subjects\n")
# Store results for further use
rs_analysis_results <- list(
model = rs_results$model,
contrasts = rs_results$contrasts,
eta_squared = rs_results$eta_squared,
plot = cross_hemi_corr_plot,
data = roi_cor_summary
)
```
##SI-3F. Ordered Age Group Contrasts Analysis
## Compare LH>RH differences between adjacent age groups
```{r}
# Adjacent Age Group Contrasts Analysis - Children Only
# Prepare data for DS1 (exclude adults)
ds1_eff_lang <- ds1_eff %>%
filter(condition == "language", age != "adult") %>%
mutate(
age = factor(age, levels = c("early", "middle", "late")),
hemi = factor(hemi, levels = c("lh", "rh"))
)
# Prepare data for DS2 (exclude adults)
ds2_eff_lang <- ds2_eff %>%
filter(condition == "EffectSize.mentalsocialphysical", age != "adult") %>%
mutate(
age = factor(age, levels = c("early", "middle", "late")),
hemi = factor(hemi, levels = c("lh", "rh"))
)
cat("DS1 observations:", nrow(ds1_eff_lang), "\n")
cat("DS2 observations:", nrow(ds2_eff_lang), "\n")
# Custom contrasts for LH>RH differences
# Order: early:lh, early:rh, middle:lh, middle:rh, late:lh, late:rh
age_contrasts <- list(
"early_vs_middle" = c( 1, -1, -1, 1, 0, 0),
"middle_vs_late" = c( 0, 0, 1, -1, -1, 1),
"early_vs_late" = c( 1, -1, 0, 0, -1, 1)
)
# DS1 Analysis
cat("\n=== DS1 ANALYSIS ===\n")
m1 <- lmer(contrast ~ hemi * age + (hemi | Subject) + (hemi | ROI), data = ds1_eff_lang)
emm1 <- emmeans(m1, ~ hemi:age)
contrast(emm1, method = age_contrasts, adjust = "sidak")
# DS2 Analysis
cat("\n=== DS2 ANALYSIS ===\n")
m2 <- lmer(contrast ~ hemi * age + (hemi | Subject) + (hemi | ROI), data = ds2_eff_lang)
emm2 <- emmeans(m2, ~ hemi:age)
contrast(emm2, method = age_contrasts, adjust = "sidak")
```
# SI-3F Lateralization Analysis - Children Only (Consecutive Age Comparisons)
## LI Analysis Consecutive
```{r}
## Data Preparation and Standardization Functions
standardize_li_data <- function(data, dataset_label) {
# Standardize Subject column
if("Subject" %in% names(data)) {
# Already standardized
} else if("Participant" %in% names(data)) {
data <- data %>% rename(Subject = Participant)
} else {
names(data)[1] <- "Subject"
}
# Ensure consistent factor levels (children only)
data <- data %>%
filter(age != "adult") %>%
mutate(age = factor(age, levels = c("early", "middle", "late")))
return(data)
}
## Statistical Analysis Function
analyze_children_li <- function(data, dataset_label) {
# Standardize data
data <- standardize_li_data(data, dataset_label)
# Pivot to wide format and calculate LI
data_wide <- data %>%
pivot_wider(names_from = hemi, values_from = contrast) %>%
mutate(LI = (lh - rh) / (lh + rh)) %>%
filter(is.finite(LI))
# Store in global environment
assign(paste0("data_wide_li_", dataset_label), data_wide, envir = .GlobalEnv)
# Fit linear model
model_li <- lm(LI ~ age, data = data_wide)
# ANOVA and effect size
anova_li <- anova(model_li)
eta_squared <- anova_li$"Sum Sq"[1] / sum(anova_li$"Sum Sq")
# Emmeans and custom contrasts
emm_li <- emmeans(model_li, ~ age)
# Custom contrasts for consecutive age comparisons
age_contrasts <- list(
"early_vs_middle" = c( 1, -1, 0),
"middle_vs_late" = c( 0, 1, -1),
"early_vs_late" = c( 1, 0, -1)
)
contrasts_li <- contrast(emm_li, method = age_contrasts, adjust = "sidak")
# Output results
cat("\n=== ", toupper(dataset_label), " LATERALIZATION INDEX ANALYSIS ===\n")
cat("================================================\n")
print(summary(model_li))
cat("\nANOVA Results:\n")
print(anova_li)
cat("\nEffect size (eta-squared):", round(eta_squared, 4), "\n")
cat("\nConsecutive Age Group Contrasts:\n")
print(contrasts_li)
return(list(
dataset = dataset_label,
model = model_li,
anova = anova_li,
contrasts = contrasts_li,
eta_squared = eta_squared,
data = data_wide
))
}
## Summary Statistics Function
create_li_children_summary <- function(data, dataset_label) {
# Prepare data
data <- standardize_li_data(data, dataset_label)
data_wide <- data %>%
pivot_wider(names_from = hemi, values_from = contrast) %>%
mutate(LI = (lh - rh) / (lh + rh)) %>%
filter(is.finite(LI))
# Summary statistics
li_summary <- data_wide %>%
group_by(age) %>%
summarise(
N = n(),
Mean_LI = round(mean(LI, na.rm = TRUE), 4),
SD_LI = round(sd(LI, na.rm = TRUE), 4),
SE_LI = round(sd(LI, na.rm = TRUE) / sqrt(n()), 4),
.groups = "drop"
)
kable(li_summary,
caption = paste(dataset_label, "Lateralization Index Summary (Children Only)"),
col.names = c("Age", "N", "Mean LI", "SD", "SE")) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
}
## Contrast Results Table Function
create_li_children_contrast_table <- function(contrast_result, dataset_label) {
contrast_df <- as.data.frame(contrast_result)
results <- contrast_df %>%
mutate(
Comparison = case_when(
contrast == "early_vs_middle" ~ "Early vs. Middle",
contrast == "middle_vs_late" ~ "Middle vs. Late",
contrast == "early_vs_late" ~ "Early vs. Late",
TRUE ~ contrast
),
B = round(estimate, 2),
SE = round(SE, 2),
t = round(t.ratio, 2),
p = case_when(
p.value < 0.001 ~ "**<0.001**",
p.value < 0.01 ~ paste0("**", sprintf("%.3f", p.value), "**"),
p.value < 0.05 ~ paste0("**", sprintf("%.3f", p.value), "**"),
TRUE ~ sprintf("%.3f", p.value)
)
) %>%
select(Comparison, B, SE, t, p)
kable(results,
col.names = c("", "***B***", "***SE***", "***t***", "***p***"),
caption = paste0("**", dataset_label, " - Consecutive Age Group Comparisons**"),
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
}
## Execute Analyses
# Dataset 1 Analysis
cat("## Dataset 1 - Volume-Based Lateralization Index Analysis\n")
# Summary table
ds1_li_summary <- create_li_children_summary(ds1_vol_sum, "Dataset 1")
print(ds1_li_summary)
# Statistical analysis
ds1_li_results <- analyze_children_li(ds1_vol_sum, "ds1")
# Contrast results table
ds1_li_contrast_table <- create_li_children_contrast_table(ds1_li_results$contrasts, "Dataset 1")
print(ds1_li_contrast_table)
# Dataset 2 Analysis
cat("\n## Dataset 2 - Volume-Based Lateralization Index Analysis\n")
# Summary table
ds2_li_summary <- create_li_children_summary(ds2_vol_sum, "Dataset 2")
print(ds2_li_summary)
# Statistical analysis
ds2_li_results <- analyze_children_li(ds2_vol_sum, "ds2")
# Contrast results table
ds2_li_contrast_table <- create_li_children_contrast_table(ds2_li_results$contrasts, "Dataset 2")
print(ds2_li_contrast_table)
## Analysis Summary
cat("\n## Analysis Summary\n")
cat("==================\n")
cat("Dataset 1:\n")
cat(" Effect size (eta-squared):", round(ds1_li_results$eta_squared, 4), "\n")
cat(" Sample size:", nrow(ds1_li_results$data), "subjects\n")
cat("\nDataset 2:\n")
cat(" Effect size (eta-squared):", round(ds2_li_results$eta_squared, 4), "\n")
cat(" Sample size:", nrow(ds2_li_results$data), "subjects\n")
# Store results for further use
li_children_analysis_results <- list(
ds1 = ds1_li_results,
ds2 = ds2_li_results
)
```
## Effect Analysis Consecutive
```{r}
# Custom contrasts for LH>RH differences between consecutive age groups
# Order: early:lh, early:rh, middle:lh, middle:rh, late:lh, late:rh
age_contrasts <- list(
"early_vs_middle" = c( 1, -1, -1, 1, 0, 0), # (early_lh - early_rh) - (middle_lh - middle_rh)
"middle_vs_late" = c( 0, 0, 1, -1, -1, 1), # (middle_lh - middle_rh) - (late_lh - late_rh)
"early_vs_late" = c( 1, -1, 0, 0, -1, 1) # (early_lh - early_rh) - (late_lh - late_rh)
)
# DS1 Analysis
cat("\n=== DS1 ANALYSIS ===\n")
m1 <- lmer(contrast ~ hemi * age + (hemi | Subject) + (hemi | ROI), data = ds1_eff_lang)
emm1 <- emmeans(m1, ~ hemi:age)
contrast(emm1, method = age_contrasts, adjust = "sidak")
# DS2 Analysis
cat("\n=== DS2 ANALYSIS ===\n")
m2 <- lmer(contrast ~ hemi * age + (hemi | Subject) + (hemi | ROI), data = ds2_eff_lang)
emm2 <- emmeans(m2, ~ hemi:age)
contrast(emm2, method = age_contrasts, adjust = "sidak")
```
# SI-4: Additional Analyses Related to Other Aspects of the Data
##SI-4A
```{r}
ds1_eff_rh <- ds1_eff %>% filter(hemi=='rh')
ds2_eff_rh <- ds2_eff %>% filter(hemi=='rh')
run_model <- function(df, dataset_label, age, region_label) {
if (nrow(df) == 0) {
return(data.frame(
Dataset = as.character(dataset_label),
Age = as.character(age),
Region = as.character(region_label),
# pairwise contrast stats
Effect_b = NA_real_, Effect_SE = NA_real_, Effect_t = NA_real_, Effect_p = NA_real_,
# omnibus/main effect stats
Main_F = NA_real_, Main_NumDF = NA_real_, Main_DenDF = NA_real_, Main_p = NA_real_,
stringsAsFactors = FALSE
))
}
# Fit the linear mixed model
mod <- lmer(effect ~ condition + (1|Subject) + (1|ROI), data = df)
# 1) extract the pairwise condition contrast (same as before)
s <- summary(mod)$coefficients
condition_rows <- s[rownames(s) != "(Intercept)", , drop=FALSE]
# Check if condition row exists
if (nrow(condition_rows) == 0) {
return(data.frame(
Dataset = as.character(dataset_label),
Age = as.character(age),
Region = as.character(region_label),
# pairwise contrast stats
Effect_b = NA_real_, Effect_SE = NA_real_, Effect_t = NA_real_, Effect_p = NA_real_,
# omnibus/main effect stats
Main_F = NA_real_, Main_NumDF = NA_real_, Main_DenDF = NA_real_, Main_p = NA_real_,
stringsAsFactors = FALSE
))
}
condition_row <- condition_rows[1, ]
# 2) extract omnibus/main effect from the ANOVA table
anov <- anova(mod)
# Check if "condition" row exists in ANOVA table
if (!"condition" %in% rownames(anov)) {
return(data.frame(
Dataset = as.character(dataset_label),
Age = as.character(age),
Region = as.character(region_label),
# pairwise contrast stats
Effect_b = NA_real_, Effect_SE = NA_real_, Effect_t = NA_real_, Effect_p = NA_real_,
# omnibus/main effect stats
Main_F = NA_real_, Main_NumDF = NA_real_, Main_DenDF = NA_real_, Main_p = NA_real_,
stringsAsFactors = FALSE
))
}
me <- anov["condition", ]
# Extract values safely with explicit conversion
effect_b <- as.numeric(condition_row["Estimate"])
effect_se <- as.numeric(condition_row["Std. Error"])
effect_t <- as.numeric(condition_row["t value"])
effect_p <- as.numeric(condition_row["Pr(>|t|)"])
main_f <- if (!is.null(me$`F value`)) as.numeric(me$`F value`) else NA_real_
main_numdf <- if (!is.null(me$NumDF)) as.numeric(me$NumDF) else NA_real_
main_dendf <- if (!is.null(me$DenDF)) as.numeric(me$DenDF) else NA_real_
main_p <- if (!is.null(me$`Pr(>F)`)) as.numeric(me$`Pr(>F)`) else NA_real_
data.frame(
Dataset = as.character(dataset_label),
Age = as.character(age),
Region = as.character(region_label),
# pairwise contrast
Effect_b = effect_b,
Effect_SE = effect_se,
Effect_t = effect_t,
Effect_p = effect_p,
# omnibus/main effect
Main_F = main_f,
Main_NumDF = main_numdf,
Main_DenDF = main_dendf,
Main_p = main_p,
stringsAsFactors = FALSE
)
}
analyze_dataset <- function(data, dataset_label) {
age_levels <- sort(unique(data$age))
out <- NULL
for (ag in age_levels) {
sub_data <- filter(data, age==ag)
res_full <- run_model(sub_data, dataset_label, ag, "Full")
res_frontal<- run_model(filter(sub_data, front=="front"), dataset_label, ag, "Frontal")
res_temporal<-run_model(filter(sub_data, front=="post"), dataset_label, ag, "Temporal")
out <- bind_rows(out, res_full, res_frontal, res_temporal)
}
out
}
# run on each dataset
res1_rh <- analyze_dataset(ds1_eff_rh, "ds1_rh_effect")
res2_rh <- analyze_dataset(ds2_eff_rh, "ds2_rh_effect")
final_results_rh <- bind_rows(res1_rh, res2_rh)
kable(final_results_rh, digits=3,
caption="SI-2: Pairwise condition estimates *and* omnibus main‐effect stats") %>%
kable_styling(full_width=TRUE)
```
## SI-4B: Analyses that treat age as a continuous variable
### DS1 Continous
```{r}
### Data Preparation for Continuous Age Analysis
# Use existing LI data created in previous sections
# From SI-3.A lateralization analysis - use the data_wide objects that were created
# Check if LI data exists from previous analyses
if(exists("data_wide_li_ds1")) {
ds1_li_wide <- data_wide_li_ds1
} else {
# Fallback: create from ds1_vol_sum if it exists
if(exists("ds1_vol_sum")) {
ds1_li_wide <- ds1_vol_sum %>%
pivot_wider(names_from = hemi, values_from = contrast) %>%
mutate(LI = (lh - rh) / (lh + rh)) %>%
filter(is.finite(LI))
} else {
stop("Cannot find existing DS1 LI data. Please run SI-3.A section first.")
}
}
if(exists("data_wide_li_ds2")) {
ds2_li_wide <- data_wide_li_ds2
} else {
# Fallback: create from ds2_vol_sum if it exists
if(exists("ds2_vol_sum")) {
ds2_li_wide <- ds2_vol_sum %>%
pivot_wider(names_from = hemi, values_from = contrast) %>%
mutate(LI = (lh - rh) / (lh + rh)) %>%
filter(is.finite(LI))
} else {
stop("Cannot find existing DS2 LI data. Please run SI-3.A section first.")
}
}
# Or use the combined LI data if available
if(exists("combined_li_plot_data")) {
# Use the combined data and split by dataset if needed
combined_li_cont_data <- combined_li_plot_data
cat("Using combined LI data from previous analysis:\n")
cat("Total subjects:", nrow(combined_li_cont_data), "\n")
} else if(exists("data_wide_li")) {
# Use the data_wide_li from the combined analysis
combined_li_cont_data <- data_wide_li %>%
group_by(Subject, age) %>%
summarise(mean_LI = mean(LI, na.rm = TRUE), .groups = "drop")
cat("Using data_wide_li from previous analysis:\n")
cat("Total subjects:", nrow(combined_li_cont_data), "\n")
}
# Compute subject-level summaries if using individual datasets
if(exists("ds1_li_wide") && !"mean_LI" %in% names(ds1_li_wide)) {
ds1_li_summary <- ds1_li_wide %>%
group_by(Subject, age) %>%
summarise(mean_LI = mean(LI, na.rm = TRUE), .groups = "drop")
} else if(exists("ds1_li_wide")) {
ds1_li_summary <- ds1_li_wide
}
if(exists("ds2_li_wide") && !"mean_LI" %in% names(ds2_li_wide)) {
ds2_li_summary <- ds2_li_wide %>%
group_by(Subject, age) %>%
summarise(mean_LI = mean(LI, na.rm = TRUE), .groups = "drop")
} else if(exists("ds2_li_wide")) {
ds2_li_summary <- ds2_li_wide
}
cat("LI data preparation completed using existing objects.\n")
if(exists("ds1_li_summary")) cat("DS1:", nrow(ds1_li_summary), "subjects\n")
if(exists("ds2_li_summary")) cat("DS2:", nrow(ds2_li_summary), "subjects\n")
if(exists("combined_li_cont_data")) cat("Combined:", nrow(combined_li_cont_data), "subjects\n")
```
### Dataset 1 Continuous Age Analysis
```{r ds1-continuous}
#A. Volume-based lateralization index (LI)
# Merge DS1 LI data with demographic data (children only)
# Use ds1_li_summary created in previous chunk
if (exists("ds1_li_summary")) {
ds1_li_cont <- ds1_li_summary %>%
left_join(ds1_demo, by = "Subject") %>%
filter(age != "adult") %>%
mutate(age_centered = Age - mean(Age, na.rm = TRUE))
} else if (exists("ds1_li_wide")) {
ds1_li_cont <- ds1_li_wide %>%
group_by(Subject, age) %>%
summarise(mean_LI = mean(LI, na.rm = TRUE), .groups = "drop") %>%
left_join(ds1_demo, by = "Subject") %>%
filter(age != "adult") %>%
mutate(age_centered = Age - mean(Age, na.rm = TRUE))
} else {
stop("DS1 LI data not available. Please run previous sections first.")
}
cat("## Dataset 1 - Continuous Age Analysis\n")
cat("=====================================\n")
cat("Sample size:", nrow(ds1_li_cont), "children\n")
cat("Age range:", round(min(ds1_li_cont$Age, na.rm = TRUE), 1), "-",
round(max(ds1_li_cont$Age, na.rm = TRUE), 1), "years\n")
# Model 1: LI predicted by continuous age
m1_ds1_li <- lm(mean_LI ~ age_centered, data = ds1_li_cont)
cat("\n### DS1 Lateralization Index ~ Age (Continuous)\n")
print(summary(m1_ds1_li))
# Calculate effect size
anova_ds1_li <- anova(m1_ds1_li)
eta_squared_ds1_li <- anova_ds1_li$"Sum Sq"[1] / sum(anova_ds1_li$"Sum Sq")
cat("\nEffect size (eta-squared):", round(eta_squared_ds1_li, 4), "\n")
#B. Magnitude-based lateralization
# Merge DS1 effect data with demographic data (children only)
ds1_eff_cont <- ds1_eff_lang %>%
left_join(ds1_demo, by = "Subject") %>%
filter(age != "adult" & !is.na(Age)) %>%
mutate(age_centered = Age - mean(Age, na.rm = TRUE))
cat("\nEffect data sample size:", length(unique(ds1_eff_cont$Subject)), "children\n")
# Model 2: Effect magnitude predicted by hemisphere and continuous age
m2_ds1_eff <- lmerTest::lmer(contrast ~ hemi * age_centered + (hemi | Subject) + (hemi | ROI),
data = ds1_eff_cont)
cat("\n### DS1 Effect Magnitude ~ Hemisphere * Age (Continuous)\n")
print(summary(m2_ds1_eff))
# Extract R-squared for mixed model
r2_ds1_eff <- MuMIn::r.squaredGLMM(m2_ds1_eff)
cat("\nMixed model R-squared:\n")
cat("Marginal R²:", round(r2_ds1_eff[1], 4), "\n")
cat("Conditional R²:", round(r2_ds1_eff[2], 4), "\n")
#C. Magnitude of the Language > Control contrast in the LH language network
```
### Dataset 2 Continuous Age Analysis
```{r ds2-continuous}
#A. Volume-based lateralization index (LI)
# Merge DS2 LI data with demographic data (children only)
ds2_cont_LI <- read.csv("../data/ds2_cont_age.csv")
# Compute subject-level mean neural effects and center age
ds2_cont_lat <- ds2_cont_LI %>%
group_by(Subject) %>%
summarise(
age_years = first(age_years),
age = first(age),
mean_LI = mean(lat, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(age_centered = age_years - mean(age_years, na.rm = TRUE))
cat("\n## Dataset 2 - Continuous Age Analysis\n")
cat("=====================================\n")
cat("Sample size:", nrow(ds2_cont_lat), "children\n")
cat("Age range:", round(min(ds2_cont_lat$age_years, na.rm = TRUE), 1), "-",
round(max(ds2_cont_lat$age_years, na.rm = TRUE), 1), "years\n")
# Model 3: LI predicted by continuous age (using subject-level data)
m3_ds2_li <- lm(mean_LI ~ age_centered, data = ds2_cont_lat)
cat("\n### DS2 Lateralization Index ~ Age (Continuous)\n")
print(summary(m3_ds2_li))
# Calculate effect size
anova_ds2_li <- anova(m3_ds2_li)
eta_squared_ds2_li <- anova_ds2_li$"Sum Sq"[1] / sum(anova_ds2_li$"Sum Sq")
cat("\nEffect size (eta-squared):", round(eta_squared_ds2_li, 4), "\n")
#B. Magnitude-based lateralization
# Merge DS2 effect data with demographic data (children only)
ds2_eff_cont <- ds2_eff_lang %>%
filter(age != "adult") %>%
mutate(age_centered = age_years - mean(age_years, na.rm = TRUE))
cat("\nEffect data sample size:", length(unique(ds2_eff_cont$Subject)), "children\n")
# Model 4: Effect magnitude predicted by hemisphere and continuous age
m4_ds2_eff <- lmerTest::lmer(contrast ~ hemi * age_centered + (hemi | Subject) + (hemi | ROI),
data = ds2_eff_cont)
cat("\n### DS2 Effect Magnitude ~ Hemisphere * Age (Continuous)\n")
print(summary(m4_ds2_eff))
# Extract R-squared for mixed model
r2_ds2_eff <- MuMIn::r.squaredGLMM(m4_ds2_eff)
cat("\nMixed model R-squared:\n")
cat("Marginal R²:", round(r2_ds2_eff[1], 4), "\n")
cat("Conditional R²:", round(r2_ds2_eff[2], 4), "\n")
```
### Combined Continuous Age Analysis
```{r combined-continuous}
# Combine both datasets for meta-analysis approach
# Check if datasets exist and have data
if (exists("ds1_li_cont") && nrow(ds1_li_cont) > 0) {
ds1_li_cont$dataset <- "DS1"
} else {
warning("ds1_li_cont is empty or doesn't exist")
}
if (exists("ds2_cont_lat") && nrow(ds2_cont_lat) > 0) {
ds2_cont_lat$dataset <- "DS2"
} else {
warning("ds2_cont_lat is empty or doesn't exist")
}
# Combine datasets
combined_li_cont <- bind_rows(ds1_li_cont, ds2_cont_lat)
# Ensure dataset is a factor with at least 2 levels
if (length(unique(combined_li_cont$dataset)) < 2) {
warning("Only one dataset available for combined analysis. Skipping combined models.")
cat("\n## Combined Dataset - Continuous Age Analysis\n")
cat("=============================================\n")
cat("WARNING: Only one dataset available. Combined analysis skipped.\n")
} else {
combined_li_cont$dataset <- as.factor(combined_li_cont$dataset)
cat("\n## Combined Dataset - Continuous Age Analysis\n")
cat("=============================================\n")
cat("Combined sample size:", nrow(combined_li_cont), "children\n")
cat("DS1 contribution:", sum(combined_li_cont$dataset == "DS1"), "subjects\n")
cat("DS2 contribution:", sum(combined_li_cont$dataset == "DS2"), "subjects\n")
# Model 5: Combined LI analysis with dataset as covariate
m5_combined_li <- lm(mean_LI ~ age_centered + dataset, data = combined_li_cont)
cat("\n### Combined Lateralization Index ~ Age + Dataset\n")
print(summary(m5_combined_li))
# Calculate effect size for age effect
anova_combined_li <- anova(m5_combined_li)
# Use correct row name (age_centered, not Age)
if ("age_centered" %in% rownames(anova_combined_li)) {
eta_squared_age <- anova_combined_li["age_centered", "Sum Sq"] / sum(anova_combined_li$"Sum Sq")
cat("\nAge effect size (eta-squared):", round(eta_squared_age, 4), "\n")
}
# Test for dataset interaction with age (use age_centered, not Age)
m5b_combined_li <- lm(mean_LI ~ age_centered * dataset, data = combined_li_cont)
cat("\n### Testing Age x Dataset Interaction\n")
print(anova(m5_combined_li, m5b_combined_li))
}
# Combined effect analysis
ds1_eff_cont$dataset <- "DS1"
ds2_eff_cont$dataset <- "DS2"
combined_eff_cont <- bind_rows(ds1_eff_cont, ds2_eff_cont)
combined_eff_cont$dataset<-as.factor(combined_eff_cont$dataset)
cat("\nCombined effect data sample size:", length(unique(combined_eff_cont$Subject)), "children\n")
# Model 6: Combined effect analysis
m6_combined_eff <- lmerTest::lmer(contrast ~ hemi * age_centered + dataset +
(hemi | Subject) + (hemi | ROI),
data = combined_eff_cont)
cat("\n### Combined Effect Magnitude ~ Hemisphere * Age + Dataset\n")
print(summary(m6_combined_eff))
# Extract R-squared for combined mixed model
r2_combined_eff <- MuMIn::r.squaredGLMM(m6_combined_eff)
cat("\nCombined mixed model R-squared:\n")
cat("Marginal R²:", round(r2_combined_eff[1], 4), "\n")
cat("Conditional R²:", round(r2_combined_eff[2], 4), "\n")
```
### Summary Table for Continuous Age Effects
```{r continuous-summary-table}
# Only create summary if combined models exist
if (exists("m5_combined_li") && exists("m6_combined_eff")) {
# Recompute combined LI eta² for the age term to avoid chunk-order issues
eta_squared_age <- {
a <- anova(m5_combined_li)
as.numeric(a["age_centered", "Sum Sq"] / sum(a[, "Sum Sq"]))
}
# Create summary table of all continuous age effects
continuous_results <- data.frame(
Analysis = c("DS1 - LI ~ Age", "DS1 - Effect ~ Hemi*Age",
"DS2 - LI ~ Age", "DS2 - Effect ~ Hemi*Age",
"Combined - LI ~ Age", "Combined - Effect ~ Hemi*Age"),
Dataset = c("DS1", "DS1", "DS2", "DS2", "Combined", "Combined"),
N_Subjects = c(length(unique(ds1_li_cont$Subject)),
length(unique(ds1_eff_cont$Subject)),
length(unique(ds2_cont_lat$Subject)),
length(unique(ds2_eff_cont$Subject)),
length(unique(combined_li_cont$Subject)),
length(unique(combined_eff_cont$Subject))),
age_centered_coefficient = c(
round(coef(m1_ds1_li)["age_centered"], 4),
round(fixef(m2_ds1_eff)["age_centered"], 4),
round(coef(m3_ds2_li)["age_centered"], 4),
round(fixef(m4_ds2_eff)["age_centered"], 4),
round(coef(m5_combined_li)["age_centered"], 4),
round(fixef(m6_combined_eff)["age_centered"], 4)
),
age_centered_p_value = c(
round(summary(m1_ds1_li)$coefficients["age_centered", "Pr(>|t|)"], 4),
round(summary(m2_ds1_eff)$coefficients["age_centered", "Pr(>|t|)"], 4),
round(summary(m3_ds2_li)$coefficients["age_centered", "Pr(>|t|)"], 4),
round(summary(m4_ds2_eff)$coefficients["age_centered", "Pr(>|t|)"], 4),
round(summary(m5_combined_li)$coefficients["age_centered", "Pr(>|t|)"], 4),
round(summary(m6_combined_eff)$coefficients["age_centered", "Pr(>|t|)"], 4)
),
Effect_Size = c(
round(eta_squared_ds1_li, 4),
round(r2_ds1_eff[1], 4),
round(eta_squared_ds2_li, 4),
round(r2_ds2_eff[1], 4),
round(eta_squared_age, 4),
round(r2_combined_eff[1], 4)
)
)
# Format p-values for significance
continuous_results$age_centered_p_formatted <- case_when(
continuous_results$age_centered_p_value < 0.001 ~ "**<0.001**",
continuous_results$age_centered_p_value < 0.01 ~ paste0("**", sprintf("%.3f", continuous_results$age_centered_p_value), "**"),
continuous_results$age_centered_p_value < 0.05 ~ paste0("**", sprintf("%.3f", continuous_results$age_centered_p_value), "**"),
TRUE ~ sprintf("%.3f", continuous_results$age_centered_p_value)
)
} else {
# Create summary without combined models
continuous_results <- data.frame(
Analysis = c("DS1 - LI ~ Age", "DS1 - Effect ~ Hemi*Age",
"DS2 - LI ~ Age", "DS2 - Effect ~ Hemi*Age"),
Dataset = c("DS1", "DS1", "DS2", "DS2"),
N_Subjects = c(length(unique(ds1_li_cont$Subject)),
length(unique(ds1_eff_cont$Subject)),
length(unique(ds2_cont_lat$Subject)),
length(unique(ds2_eff_cont$Subject))),
age_centered_coefficient = c(
round(coef(m1_ds1_li)["age_centered"], 4),
round(fixef(m2_ds1_eff)["age_centered"], 4),
round(coef(m3_ds2_li)["age_centered"], 4),
round(fixef(m4_ds2_eff)["age_centered"], 4)
),
age_centered_p_value = c(
round(summary(m1_ds1_li)$coefficients["age_centered", "Pr(>|t|)"], 4),
round(summary(m2_ds1_eff)$coefficients["age_centered", "Pr(>|t|)"], 4),
round(summary(m3_ds2_li)$coefficients["age_centered", "Pr(>|t|)"], 4),
round(summary(m4_ds2_eff)$coefficients["age_centered", "Pr(>|t|)"], 4)
),
Effect_Size = c(
round(eta_squared_ds1_li, 4),
round(r2_ds1_eff[1], 4),
round(eta_squared_ds2_li, 4),
round(r2_ds2_eff[1], 4)
)
)
# Format p-values for significance
continuous_results$age_centered_p_formatted <- case_when(
continuous_results$age_centered_p_value < 0.001 ~ "**<0.001**",
continuous_results$age_centered_p_value < 0.01 ~ paste0("**", sprintf("%.3f", continuous_results$age_centered_p_value), "**"),
continuous_results$age_centered_p_value < 0.05 ~ paste0("**", sprintf("%.3f", continuous_results$age_centered_p_value), "**"),
TRUE ~ sprintf("%.3f", continuous_results$age_centered_p_value)
)
}
# Create formatted table
kable(continuous_results %>% select(-age_centered_p_value),
col.names = c("Analysis", "Dataset", "N", "Age β", "Age p", "Effect Size"),
caption = "**Summary of Continuous Age Effects Across All Analyses**",
escape = FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
add_header_above(c(" " = 3, "Age Effect" = 2, " " = 1))
cat("\n## Key Findings from Continuous Age Analysis\n")
cat("==========================================\n")
cat("- All models used children only (adults excluded)\n")
cat("- Age effects tested as continuous predictor variable\n")
cat("- LI = Lateralization Index, Effect = Activation magnitude\n")
cat("- Effect sizes: eta² for LI models, marginal R² for mixed models\n")
# Store results for further analysis
continuous_age_results <- list(
ds1_li_model = m1_ds1_li,
ds1_effect_model = m2_ds1_eff,
ds2_li_model = m3_ds2_li,
ds2_effect_model = m4_ds2_eff,
combined_li_model = m5_combined_li,
combined_effect_model = m6_combined_eff,
summary_table = continuous_results
)
```
#SI-4C: Age-related changes in the magnitude of the Language > Control contrast in the LH language network controlling for motion. [parallels main Results section 3]
## 2) Do different properties of the language network change between child ages (controlling for motion)?
### Effect Size Analysis
#### Dataset 1
```{r effect-size-ds1}
# Filter for language condition, left hemisphere, children only (exclude adults)
data_all_l_eff_lh_ds1 <- ds1_eff %>%
filter(condition == 'language',
hemi == 'lh') %>%
# age != 'adult') %>%
dplyr::select(Subject, ROI, age, front, hemi, contrast, num_outliers)
# Merge with demographic data and create mean_cont_df_ds1
# ds1_demo already has Subject column from CSV, no need to create from ID
cont_df_ds1 <- merge(data_all_l_eff_lh_ds1, ds1_demo, by = "Subject")
# Calculate subject-level means
mean_cont_df_ds1 <- cont_df_ds1 %>%
group_by(Subject) %>%
summarize(mean_num_outliers = mean(num_outliers),
mean_con = mean(contrast),
mean_age = mean(Age),
age_group = first(Set))
# Now run models that depend on mean_cont_df_ds1
m3_ds1 <- lm(mean_num_outliers ~ mean_age, data = mean_cont_df_ds1)
summary(m3_ds1)
motion_group_ds1 <- lm(mean_num_outliers ~ age_group, data = mean_cont_df_ds1)
summary(motion_group_ds1)
lsmeans(motion_group_ds1, list(pairwise ~ age_group), adjust = "tukey")
# Mixed effects model: effect size ~ motion + age
m1_ds1 <- lmer(contrast ~ num_outliers + age + (1|Subject) + (1|ROI),
data = data_all_l_eff_lh_ds1)
emm1 <- emmeans(m1_ds1, ~ age)
contrast(emm1, method = "trt.vs.ctrl", ref = "adult", adjust = "sidak")
```
#### Dataset 2
```{r effect-size-ds2}
ds2_outliers<-read.csv("../data/ds2_outliers.csv")
ds2_outliers$Subject<-ds2_outliers$SubID
# Run motion comparison by age
ds2_eff_v2<-merge(ds2_eff, ds2_outliers,by="Subject")%>% filter(condition=="EffectSize.mentalsocialphysical")
ds2_participant <- ds2_eff_v2 %>%
group_by(Subject) %>%
summarise(
# Average the outcome variable
num_outliers_avg = mean(num_outliers, na.rm = TRUE),
# Keep essential variables only
age = first(age),
# Count how many observations we averaged (useful for checking)
n_observations = n(),
.groups = 'drop'
)
# Check the result
head(ds2_participant)
cat("Original data rows:", nrow(ds2_eff_v2), "\n")
cat("Participant-level data rows:", nrow(ds2_participant), "\n")
cat("Observations per participant (should be consistent):", unique(ds2_participant$n_observations), "\n")
# Now run a simple linear model
m_ds2_simple <- lm(num_outliers_avg ~ age, data = ds2_participant)
anova(m_ds2_simple)
data_all_l_eff_lh_ds2 <- ds2_eff_v2 %>%
dplyr::filter(hemi == 'lh') %>%
dplyr::select(Subject, ROI, age, front, hemi, contrast, num_outliers)
# Mixed effects model: effect size ~ motion + age
m1_ds2 <- lmer(contrast ~ num_outliers + age + (1|Subject) + (1|ROI),
data = data_all_l_eff_lh_ds2)
summary(m1_ds2)
emm2 <- emmeans(m1_ds2, ~ age)
contrast(emm2, method = "trt.vs.ctrl", ref = "adult", adjust = "sidak")
```
#### Dataset 1
```{r continuous-age-ds1}
# SI-4B Part C: Continuous age analysis for magnitude (parallel to Table SI-4B part C)
# Linear regression: magnitude ~ age (children only, using mean_cont_df_ds1)
# Filter mean_cont_df_ds1 to children only (exclude adults)
# mean_cont_df_ds1 already created in effect-size-ds1 chunk with mean_con column
mean_cont_df_ds1_children <- mean_cont_df_ds1 %>%
filter(age_group != "adult")
# Fit linear regression model
m_cont_age_ds1 <- lm(mean_con ~ mean_age, data = mean_cont_df_ds1_children)
cat("\n## SI-4B Part C: Age as continuous predictor of LH magnitude (DS1)\n")
cat("Linear regression: magnitude ~ Age (children only)\n")
summary(m_cont_age_ds1)
```
### Resting State (Dataset 1 only)
```{r resting-state-ds1}
# SI-4C Part B: IRFC analysis controlling for motion
# Table SI-4C: Inter-regional functional correlation strength in LH network, controlling for motion
# Load pre-summarized resting state data with age
# This file has: mean_rs, sub (Subject), age, hemi
ds1_rs_sum <- read.csv("../data/ds1_rs_sum_031025.csv")
# Filter to LH only and rename columns
ds1_rs_lh <- ds1_rs_sum %>%
filter(hemi == "lh") %>%
rename(Subject = sub, mean_r = mean_rs)
cat("\n## SI-4C Part B: IRFC controlling for motion (DS1)\n")
cat("Note: Analysis uses age groups; motion covariate integration pending\n")
# Model with age groups (parallel to Table SI-4C part B)
m_rs_age <- lm(mean_r ~ age, data = ds1_rs_lh)
summary(m_rs_age)
# Pairwise comparisons
emm_rs <- emmeans(m_rs_age, ~ age)
contrast(emm_rs, method = "trt.vs.ctrl", ref = "adult", adjust = "sidak")
# Store for use in SI-4D
ds1_rs_summary <- ds1_rs_lh
```
#SI-4D. Ordered age comparisons for the magnitude of the Language > Control contrast and the strength of functional correlations in the LH language network.
## A. Ordered age comparisons for magnitude
```{r si-4d-magnitude}
# Table SI-4D Part A: Ordered age comparisons for LH magnitude
# Linear mixed-effects: magnitude ~ age, pairwise comparisons among child groups
# Dataset 1: Use data_all_l_eff_lh_ds1 (created in effect-size-ds1 chunk)
# Filter to children only
ds1_mag_children <- data_all_l_eff_lh_ds1 %>%
filter(age != "adult")
# Fit linear mixed-effects model
m_mag_ord_ds1 <- lmer(contrast ~ age + (1|Subject) + (1|ROI), data = ds1_mag_children)
summary(m_mag_ord_ds1)
# Pairwise comparisons among child age groups
emm_mag_ds1 <- emmeans(m_mag_ord_ds1, ~ age)
pairs_mag_ds1 <- pairs(emm_mag_ds1, adjust = "sidak")
cat("\n## SI-4D Part A: Ordered age comparisons for magnitude (DS1)\n")
print(pairs_mag_ds1)
# Dataset 2: Similar analysis
ds2_mag_children <- data_all_l_eff_lh_ds2 %>%
filter(age != "adult")
m_mag_ord_ds2 <- lmer(contrast ~ age + (1|Subject) + (1|ROI), data = ds2_mag_children)
summary(m_mag_ord_ds2)
emm_mag_ds2 <- emmeans(m_mag_ord_ds2, ~ age)
pairs_mag_ds2 <- pairs(emm_mag_ds2, adjust = "sidak")
cat("\n## SI-4D Part A: Ordered age comparisons for magnitude (DS2)\n")
print(pairs_mag_ds2)
```
## B. Ordered age comparisons for IRFC
```{r si-4d-irfc}
# Table SI-4D Part B: Ordered age comparisons for IRFC strength
# Linear regression: IRFC ~ age, pairwise comparisons among child groups (DS1 only)
# Use ds1_rs_summary created in resting-state-ds1 chunk
# Filter to children only
ds1_rs_children <- ds1_rs_summary %>%
filter(age != "adult")
# Fit linear regression model
m_irfc_ord_ds1 <- lm(mean_r ~ age, data = ds1_rs_children)
summary(m_irfc_ord_ds1)
# Pairwise comparisons among child age groups
emm_irfc_ds1 <- emmeans(m_irfc_ord_ds1, ~ age)
pairs_irfc_ds1 <- pairs(emm_irfc_ds1, adjust = "sidak")
cat("\n## SI-4D Part B: Ordered age comparisons for IRFC (DS1 only)\n")
print(pairs_irfc_ds1)
```

NCOMM_OzernovPalchik_OBrien_SI_code_091025.Rmd, no license · at the source

Overview

Authors: Ola Ozernov-Palchik^1,2, Amanda M. O’Brien^1,3, Elizabeth Jiachen Lee^1,4, Hilary Richardson⁵, Rachel Romeo⁶, Moshe Poliak⁴, Benjamin Lipkin^1,4, Hannah Small⁷, Jimmy Capella⁸, Alfonso Nieto-Castañón¹, Rebecca Saxe^1,4, John D. E. Gabrieli^1,4, Evelina Fedorenko^1,3,4

ORCID iDs: Ola Ozernov-Palchik, Hilary Richardson, Rachel Romeo, Rebecca Saxe, Evelina Fedorenko

Institutions: Massachusetts Institute of Technology (United States); Wheelock College (United States); Harvard University (United States); University of Edinburgh (United Kingdom); University of Maryland, College Park (United States); Johns Hopkins University (United States); University of North Carolina at Chapel Hill (United States)

Journal: Nature communications, volume 17, issue 1, article 6505

Dates: received 16 July 2024; accepted 23 April 2026; published online 16 May 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41467-026-72916-5 · PMID 42143030 · PMCID PMC13376885 · OpenAlex W4396973866

Open access: gold, a free copy (OpenAlex)

Preprint: doi:10.1101/2024.05.15.594172

Status: code verified

Categories: fMRI (modality), human (organism)

Methods: Spectral & time-frequency, Connectivity, Statistics, Machine learning, fMRI & imaging

Keywords: Language, Human behaviour

MeSH: Brain*, Functional Laterality*, Language*, Language Development*, Magnetic Resonance Imaging*, Adolescent, Brain Mapping, Child, Child, Preschool, Female, Humans, Male (* major topic)

Topic: Hemispheric Asymmetry in Neuroscience (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: U.S. Department of Health & Human Services | NIH | Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD) (HD100064)

Citations: cited by 11 papers (Europe PMC); 124 references in the paper

Abstract

In adults, left hemisphere (LH) damage often leads to aphasia, but many cases of early damage leave linguistic processing intact, with a functional language system developing in the right hemisphere. To explain this early apparent equipotentiality of the two hemispheres for language, some have proposed that the language system is more bilateral during early development and becomes increasingly left-lateralized with age. We examined language lateralization using fMRI in two large developmental cohorts (total n = 273 children aged 4-16 years; n = 107 adults). Strong, adult-like LH lateralization (in response magnitude and activation volume) was evident by age 4, although other features of the LH language network showed protracted development, including the magnitude of language response and the strength of functional connectivity. Thus, although the RH can take over language function in some cases of early brain damage, this plasticity occurs in spite of adult-level LH bias present by age 4 years.

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 2 matches between paragraphs and lines of code.

OSF 3mvpx

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Size: 39 files

Software Heritage: not checked

Found in: “Code availability”

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

Tools: ggplot2 (5 files), ggpubr (5 files), lme4 (5 files), lmerTest (5 files), tidyverse (5 files), cowplot (3 files), brms (2 files), emmeans (2 files), SPM (2 files)

Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)

28 September 2026: the link answers (HTTP 200)

8 files

At the source: osf.io/3mvpx/

web.conn-toolbox.org/resources/conn-extensions

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: the link answers

Software Heritage: not checked

Found in: “Code availability”

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

Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)

28 September 2026: the link answers (HTTP 200)

At the source: web.conn-toolbox.org/resources/conn-extensions/evlab

Code availability

Custom scripts used for fMRI preprocessing, functional localization, lateralization analyses, and Bayesian modeling are available at https://web.conn-toolbox.org/resources/conn-extensions/evlab and https://osf.io/3mvpx/. Additional analysis pipelines rely on publicly available software, including SPM12 and Wilke’s LI Toolbox.

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

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;
7 scripts, each with its path and the digest of its content;
2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

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

Data

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

Data availability

The source data underlying all main figures and tables are provided in a source data file accompanying this paper. Individual whole-brain activation maps and extracted data frames for all reported analyses have been deposited in the Open Science Framework (https://osf.io/3mvpx/). Consistent with historical data-sharing practices and participant consent constraints, raw fMRI data are not publicly released; for most legacy datasets, only activation maps are available, whereas raw data sharing for more recent participants depends on individual consent and must be evaluated case by case. To support reproducibility of the present findings, we provide group-level statistical maps for the language > control contrast. Data from the Olulade et al. (2020) replication analyses (Supplementary Information Sections 1–3) remain under the stewardship of the original author team and may be requested from Elissa Newport, with responses typically provided within 30 days. Source data are provided with this paper.

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, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 13 authors, 2 keywords, 12 MeSH terms, 1 funder, 112 references.

Cite

This paper

Ozernov-Palchik, O., O’Brien, A. M., Lee, E. J., Richardson, H., Romeo, R., Poliak, M., Lipkin, B., Small, H., Capella, J., Nieto-Castañón, A., Saxe, R., Gabrieli, J. D. E., & Fedorenko, E. (2026). Precision fMRI reveals that the language network exhibits adult-like left-hemispheric lateralization by 4 years of age. Nature communications, 17(1), 6505. https://doi.org/10.1038/s41467-026-72916-5

BibTeX

@article{ozernovpalchik2026precision, author = {Ozernov-Palchik, Ola and O’Brien, Amanda M. and Lee, Elizabeth Jiachen and Richardson, Hilary and Romeo, Rachel and Poliak, Moshe and Lipkin, Benjamin and Small, Hannah and Capella, Jimmy and Nieto-Castañón, Alfonso and Saxe, Rebecca and Gabrieli, John D. E. and Fedorenko, Evelina}, title = {{Precision fMRI reveals that the language network exhibits adult-like left-hemispheric lateralization by 4 years of age}}, journal = {Nature communications}, year = {2026}, month = may, volume = {17}, number = {1}, pages = {6505}, publisher = {Nature Publishing Group}, issn = {2041-1723}, doi = {10.1038/s41467-026-72916-5}, url = {https://doi.org/10.1038/s41467-026-72916-5}, pmid = {42143030}, pmcid = {PMC13376885} }

RIS

TY - JOUR AU - Ozernov-Palchik, Ola AU - O’Brien, Amanda M. AU - Lee, Elizabeth Jiachen AU - Richardson, Hilary AU - Romeo, Rachel AU - Poliak, Moshe AU - Lipkin, Benjamin AU - Small, Hannah AU - Capella, Jimmy AU - Nieto-Castañón, Alfonso AU - Saxe, Rebecca AU - Gabrieli, John D. E. AU - Fedorenko, Evelina TI - Precision fMRI reveals that the language network exhibits adult-like left-hemispheric lateralization by 4 years of age T2 - Nature communications J2 - Nat Commun PY - 2026 DA - 2026/05/16 VL - 17 IS - 1 SP - 6505 SN - 2041-1723 PB - Nature Publishing Group DO - 10.1038/s41467-026-72916-5 UR - https://doi.org/10.1038/s41467-026-72916-5 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41467-026-72916-5", "type": "article-journal", "title": "Precision fMRI reveals that the language network exhibits adult-like left-hemispheric lateralization by 4 years of age", "container-title": "Nature communications", "author": [ { "family": "Ozernov-Palchik", "given": "Ola" }, { "family": "O’Brien", "given": "Amanda M." }, { "family": "Lee", "given": "Elizabeth Jiachen" }, { "family": "Richardson", "given": "Hilary" }, { "family": "Romeo", "given": "Rachel" }, { "family": "Poliak", "given": "Moshe" }, { "family": "Lipkin", "given": "Benjamin" }, { "family": "Small", "given": "Hannah" }, { "family": "Capella", "given": "Jimmy" }, { "family": "Nieto-Castañón", "given": "Alfonso" }, { "family": "Saxe", "given": "Rebecca" }, { "family": "Gabrieli", "given": "John D. E." }, { "family": "Fedorenko", "given": "Evelina" } ], "container-title-short": "Nat Commun", "volume": "17", "issue": "1", "page": "6505", "DOI": "10.1038/s41467-026-72916-5", "PMID": "42143030", "PMCID": "PMC13376885", "ISSN": "2041-1723", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41467-026-72916-5", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 16 ] ] } }

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

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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

Precision fMRI reveals that the language network exhibits adult-like left-hemispheric lateralization by 4 years of age.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

OSF 3mvpx

web.conn-toolbox.org/resources/conn-extensions

Code availability

Tracing map

Data

Data availability

Versions

Version 1, 28 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity