Distinct Representations of Irrelevant Emotional Information in the Visual Network Are Associated With Psychopathology in Youth.

Code ↔ Paper

The paper beside its authors' code: matches between them have not been computed for this paper yet.

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 · 1,900 lines · 89 KB · no license

---
title: "visual_emotion_processing_psychopathology"
output: html_document
date: "2025-11-25"
---
# Main steps
- part 1 - general setup
- part 2 - dissimilarity and total problem score
- part 3 - dissimilarity and performance
- part 4 - effect sizes
- part 5 - performance and total problem score
- part 6 - demographic table
- part 7 - mediation
note: Pearson correlations was calculated by python RSAtoolbox package
last update on: 2026-04-24
author: Yen-Chu Lin
# part 1 - general setup
- (1) load packages
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(lme4)
library(lmerTest)
library(sjPlot)
require(MuMIn)
library(effsize) # for cohen.d
library(tableone)
library(vtable)
library(summarytools)
library(MatchIt)
```
- (2) functions
```{r}
select_condition_and_combine_dataframes_v1 <- function(
r_visual, r_limbic, r_dorsal, r_ventral, r_frontal, r_default, r_somatic, r_v1, col_name){
# the function subsets dataframe based on condition selected & combine all dataframes including v1
# also see select_conditions & combine_networks
# subset data based on col_name
rvis <- r_visual[c("src_subject_id",col_name)]
rlim <- r_limbic[c("src_subject_id",col_name)]
rdor <- r_dorsal[c("src_subject_id",col_name)]
rven <- r_ventral[c("src_subject_id",col_name)]
rfro <- r_frontal[c("src_subject_id",col_name)]
rdef <- r_default[c("src_subject_id",col_name)]
rsom <- r_somatic[c("src_subject_id",col_name)]
rv1 <- r_v1[c("src_subject_id",col_name)]
colnames(rvis)[2] <- "vis"
colnames(rlim)[2] <- "lim"
colnames(rdor)[2] <- "dor"
colnames(rven)[2] <- "ven"
colnames(rfro)[2] <- "fro"
colnames(rdef)[2] <- "def"
colnames(rsom)[2] <- "som"
colnames(rv1)[2] <- "v1"
# combine all 7 networks & v1 into one dataframe
data <- merge(rvis, rlim, by="src_subject_id")
data <- merge(data, rdor, by="src_subject_id")
data <- merge(data, rven, by="src_subject_id")
data <- merge(data, rfro, by="src_subject_id")
data <- merge(data, rdef, by="src_subject_id")
data <- merge(data, rsom, by="src_subject_id")
data <- merge(data, rv1, by="src_subject_id")
# return data
return(data)
}
func_merge_df <- function(data_r, fam_site_id, nback) {
# merge with family ad site ids
data_full <- merge(data_r, fam_site_id, by="src_subject_id")
# merge with nback
data_full <- merge(data_full, nback, by="src_subject_id")
return(data_full)
}
vector_to_factor <- function(data_long) {
data_long$src_subject_id <- as.factor(data_long$src_subject_id)
data_long$condition <- as.factor(data_long$condition)
data_long$rel_family_id <- as.factor(data_long$rel_family_id)
data_long$site_id_l <- as.factor(data_long$site_id_l)
data_long$network <- as.factor(data_long$network)
return(data_long)
}
```
- (3) inputs
```{r}
# file path
r_file_path <- "/Volumes/FABlab drive 2/Processed data/RSA_Pearson_correlations/"
abcd_file_path <- "/Volumes/FABlab drive 2/Data/abcd-data-release-4.0/"
subj_list_path <- "/Volumes/FABlab drive 2/Processed data/"
# dissimilarity (Pearson correlation files)
file1 <- "visual_pearson_corr.csv"
file2 <- "limbic_pearson_corr_0back.csv"
file3 <- "limbic_pearson_corr_2back.csv"
file4 <- "doratt_pearson_corr_0back.csv"
file5 <- "doratt_pearson_corr_2back.csv"
file6 <- "venatt_pearson_corr_0back.csv"
file7 <- "venatt_pearson_corr_2back.csv"
file8 <- "sommot_pearson_corr_0back.csv"
file9 <- "sommot_pearson_corr_2back.csv"
file10 <- "control_pearson_corr_0back.csv"
file11 <- "control_pearson_corr_2back.csv"
file12 <- "default_pearson_corr_0back.csv"
file13 <- "default_pearson_corr_2back.csv"
file14 <- "v1_pearson_corr.csv"
# CBCL scores
cbcl_filename <- "abcd_cbcls01.txt"
cbcl_col_select <-
c("subjectkey",
"cbcl_scr_dsm5_adhd_t","cbcl_scr_dsm5_adhd_r",
"cbcl_scr_syn_attention_t","cbcl_scr_syn_attention_r",
"cbcl_scr_syn_external_t","cbcl_scr_syn_external_r",
"cbcl_scr_syn_internal_t","cbcl_scr_syn_internal_r",
"cbcl_scr_syn_totprob_t","cbcl_scr_syn_totprob_r",
"cbcl_scr_syn_anxdep_r","cbcl_scr_syn_withdep_r",
"cbcl_scr_syn_somatic_r","cbcl_scr_syn_social_r",
"cbcl_scr_syn_thought_r","cbcl_scr_syn_rulebreak_r",
"cbcl_scr_syn_aggressive_r"
)
# N-back behavioral data (accuracy)
naback_behav_filename <- "abcd_mrinback02.txt"
nback_col_select <-
c("subjectkey",
"tfmri_nb_all_beh_c0bpf_rate","tfmri_nb_all_beh_c0bnf_rate","tfmri_nb_all_beh_c0bngf_rate",
"tfmri_nb_all_beh_c2bpf_rate","tfmri_nb_all_beh_c2bnf_rate","tfmri_nb_all_beh_c2bngf_rate",
"tfmri_nb_all_beh_c0b_rate","tfmri_nb_all_beh_c2b_rate"
)
# list of clean subjects
subj_filename <- "cleanedsubjectkeys.csv"
# family income, caregiver education, and child race/ethnicity
pdemo_filename <- "pdem02.txt"
pdemo_col_select <-
c("subjectkey","demo_prnt_ed_v2","demo_prtnr_ed_v2",
"demo_prnt_income_v2", "demo_prtnr_income_v2","demo_comb_income_v2",
"demo_race_a_p___10", "demo_race_a_p___11", "demo_race_a_p___12", "demo_race_a_p___13",
"demo_race_a_p___14", "demo_race_a_p___15", "demo_race_a_p___16", "demo_race_a_p___17",
"demo_race_a_p___18", "demo_race_a_p___19", "demo_race_a_p___20", "demo_race_a_p___21",
"demo_race_a_p___22", "demo_race_a_p___23", "demo_race_a_p___24", "demo_race_a_p___25",
"demo_race_a_p___77", "demo_race_a_p___99","demo_ethn_v2")
```
- (4) load files
```{r}
# load the dissimilarity files
r_vis <- read.csv(paste(r_file_path, file1, sep=""), header = TRUE)
r_lim_0b <- read.csv(paste(r_file_path, file2, sep=""), header = TRUE)
r_lim_2b <- read.csv(paste(r_file_path, file3, sep=""), header = TRUE)
r_dor_0b <- read.csv(paste(r_file_path, file4, sep=""), header = TRUE)
r_dor_2b <- read.csv(paste(r_file_path, file5, sep=""), header = TRUE)
r_ven_0b <- read.csv(paste(r_file_path, file6, sep=""), header = TRUE)
r_ven_2b <- read.csv(paste(r_file_path, file7, sep=""), header = TRUE)
r_som_0b <- read.csv(paste(r_file_path, file8, sep=""), header = TRUE)
r_som_2b <- read.csv(paste(r_file_path, file9, sep=""), header = TRUE)
r_fro_0b <- read.csv(paste(r_file_path, file10, sep=""), header = TRUE)
r_fro_2b <- read.csv(paste(r_file_path, file11, sep=""), header = TRUE)
r_def_0b <- read.csv(paste(r_file_path, file12, sep=""), header = TRUE)
r_def_2b <- read.csv(paste(r_file_path, file13, sep=""), header = TRUE)
r_v1 <- read.csv(paste(r_file_path, file14, sep=""), header = TRUE)
r_vis_0b <- r_vis[c("src_subject_id","X0happy0fear","X0happy0neu","X0neu0fear")]
r_vis_2b <- r_vis[c("src_subject_id","X2happy2fear","X2happy2neu","X2neu2fear")]
r_v1_0b <- r_v1[c("src_subject_id","X0happy0fear","X0happy0neu","X0neu0fear")]
r_v1_2b <- r_v1[c("src_subject_id","X2happy2fear","X2happy2neu","X2neu2fear")]
# load cbcl
abcd_cbcl <- read.delim(paste(abcd_file_path, cbcl_filename, sep=""), header = TRUE) %>% filter(eventname == "2_year_follow_up_y_arm_1")
abcd_cbcl <- abcd_cbcl[cbcl_col_select]
abcd_cbcl$subjectkey <- gsub("NDAR_", "", abcd_cbcl$subjectkey)
names(abcd_cbcl)[names(abcd_cbcl) == 'subjectkey'] <- 'src_subject_id'
# load the file with family_id and subset the columns we need
acspsw <- read.delim(paste(abcd_file_path, "acspsw03.txt", sep=""), header = TRUE) %>%
slice(-1) %>%
filter(eventname == "baseline_year_1_arm_1") %>% # only has data from baseline & year 1
subset(select = c("subjectkey", "rel_family_id"))
# load the file with site_id and subset the columns we need
abcd_site <- read.delim(paste(abcd_file_path, "abcd_lt01.txt", sep=""), header = TRUE) %>%
slice(-1) %>%
filter(eventname == "2_year_follow_up_y_arm_1") %>%
subset(select = c("subjectkey", "site_id_l","interview_age","sex"))
# combine family id and site id
fam_site_id <- merge(acspsw, abcd_site, by = "subjectkey")
fam_site_id$subjectkey <- gsub("NDAR_", "", fam_site_id$subjectkey)
names(fam_site_id)[names(fam_site_id) == 'subjectkey'] <- 'src_subject_id'
# load a clean list of subject to use in the neuroimaging analysis
subj_list <- read.csv(paste(subj_list_path, subj_filename, sep=""), header = FALSE) %>%
slice(-1)
# load N-back behavioral data
nback <- read.delim(paste(abcd_file_path, naback_behav_filename, sep=""), header = TRUE) %>% filter(eventname == "2_year_follow_up_y_arm_1")
nback <- nback[nback_col_select] %>% na.omit()
nback$subjectkey <- gsub("NDAR_", "", nback$subjectkey)
names(nback)[names(nback) == 'subjectkey'] <- 'src_subject_id'
# load family income and caregiver education data
abcd_pdem <- read.delim(paste(abcd_file_path, pdemo_filename, sep=""), header = TRUE) %>%
slice(-1) %>%
select(all_of(pdemo_col_select))
abcd_pdem <- abcd_pdem[pdemo_col_select] %>% na.omit()
abcd_pdem$subjectkey <- gsub("NDAR_", "", abcd_pdem$subjectkey)
names(abcd_pdem)[names(abcd_pdem) == 'subjectkey'] <- 'src_subject_id'
# time 3 (year 3 follow-up) CBCL scores
abcd_cbcl_3year <- read.delim(paste(abcd_file_path, cbcl_filename, sep=""), header = TRUE) %>%
filter(eventname == "3_year_follow_up_y_arm_1") %>%
select(all_of(c("src_subject_id","cbcl_scr_syn_totprob_r")))
abcd_cbcl_3year$src_subject_id <- gsub("NDAR_", "", abcd_cbcl_3year$src_subject_id)
abcd_cbcl_3year$cbcl_scr_syn_totprob_r <- as.numeric(abcd_cbcl_3year$cbcl_scr_syn_totprob_r)
# rename the columns
names(abcd_cbcl_3year)[names(abcd_cbcl_3year) == 'cbcl_scr_syn_totprob_r'] <- 'cbcl_scr_syn_totprob_r_3year'
# check for completeness
abcd_cbcl_3year_totprob_r <- abcd_cbcl_3year[complete.cases(abcd_cbcl_3year$cbcl_scr_syn_totprob_r_3year), ]
dim(abcd_cbcl_3year_totprob_r) # 6133 2
```
- (5) preprocess dataframes
```{r}
# dissimilarity
## SUBSETTING DATA BASED ON CONDITION SELECTED
### 0-back ###
col_name1 <- "X0happy0fear"
col_name2 <- "X0happy0neu"
col_name3 <- "X0neu0fear"
# subset and combine dataframes
data_hf_0b <- select_condition_and_combine_dataframes_v1(
r_vis_0b, r_lim_0b, r_dor_0b, r_ven_0b, r_fro_0b, r_def_0b, r_som_0b, r_v1_0b, col_name1)
data_hn_0b <- select_condition_and_combine_dataframes_v1(
r_vis_0b, r_lim_0b, r_dor_0b, r_ven_0b, r_fro_0b, r_def_0b, r_som_0b, r_v1_0b, col_name2)
data_fn_0b <- select_condition_and_combine_dataframes_v1(
r_vis_0b, r_lim_0b, r_dor_0b, r_ven_0b, r_fro_0b, r_def_0b, r_som_0b, r_v1_0b, col_name3)
### 2-back ###
col_name4 <- "X2happy2fear"
col_name5 <- "X2happy2neu"
col_name6 <- "X2neu2fear"
# subset and combine dataframes
data_hf_2b <- select_condition_and_combine_dataframes_v1(
r_vis_2b, r_lim_2b, r_dor_2b, r_ven_2b, r_fro_2b, r_def_2b, r_som_2b, r_v1_2b, col_name4)
data_hn_2b <- select_condition_and_combine_dataframes_v1(
r_vis_2b, r_lim_2b, r_dor_2b, r_ven_2b, r_fro_2b, r_def_2b, r_som_2b, r_v1_2b, col_name5)
data_fn_2b <- select_condition_and_combine_dataframes_v1(
r_vis_2b, r_lim_2b, r_dor_2b, r_ven_2b, r_fro_2b, r_def_2b, r_som_2b, r_v1_2b, col_name6)
# sanity check
# exclude subjects with NA for dissimilarity
# dim(data_hf_0b) # 4952 9
# dim(data_hf_2b) # 4952 9
# dim(data_hn_0b) # 4952 9
# dim(data_hn_2b) # 4952 9
# dim(data_fn_0b) # 4952 9
# dim(data_fn_2b) # 4952 9
# check completeness
# dim(data_hf_0b[complete.cases(data_hf_0b), ]) # 4952 9
# dim(data_hf_2b[complete.cases(data_hf_2b), ]) # 4952 9
# dim(data_hn_0b[complete.cases(data_hn_0b), ]) # 4952 9
# dim(data_hn_2b[complete.cases(data_hn_2b), ]) # 4952 9
# dim(data_fn_0b[complete.cases(data_fn_0b), ]) # 4952 9
# dim(data_fn_2b[complete.cases(data_fn_2b), ]) # 4952 9
# head(data_hf_0b)
```
# part 2 - dissimilarity and CBCL score
- (1) process datasets for dissimilarity and total problems score
```{r}
# merge dataframes
### 0-back ###
# (1) add a column to specify the condition
data_hf_0b$condition <- "hf"
data_hn_0b$condition <- "hn"
data_fn_0b$condition <- "fn"
# (2) merge distance with family id and site id
data_hf0b <- func_merge_df(data_hf_0b, fam_site_id, nback)
data_hn0b <- func_merge_df(data_hn_0b, fam_site_id, nback)
data_fn0b <- func_merge_df(data_fn_0b, fam_site_id, nback)
# (3) merge with cbcl
data_hf0b <- merge(data_hf0b, abcd_cbcl, by="src_subject_id")
data_hn0b <- merge(data_hn0b, abcd_cbcl, by="src_subject_id")
data_fn0b <- merge(data_fn0b, abcd_cbcl, by="src_subject_id")
# (4) convert the data into long format
colname_list <- c("vis","lim","dor","ven","fro","def","som","v1")
data_hf0b_long <-
pivot_longer(data_hf0b, cols = all_of(colname_list), names_to = "network", values_to = "dis")
data_hn0b_long <-
pivot_longer(data_hn0b, cols = all_of(colname_list), names_to = "network", values_to = "dis")
data_fn0b_long <-
pivot_longer(data_fn0b, cols = all_of(colname_list), names_to = "network", values_to = "dis")
# (5) COMBINE ALL THREE CONDITIONS
data_0back <- rbind(data_hf0b_long, data_hn0b_long, data_fn0b_long)
# (6) make age and CBCL scores numeric
cbcl_col_numeric <-
c("interview_age",
"cbcl_scr_dsm5_adhd_t","cbcl_scr_dsm5_adhd_r",
"cbcl_scr_syn_attention_t","cbcl_scr_syn_attention_r",
"cbcl_scr_syn_external_t","cbcl_scr_syn_external_r",
"cbcl_scr_syn_internal_t","cbcl_scr_syn_internal_r",
"cbcl_scr_syn_totprob_t","cbcl_scr_syn_totprob_r",
"cbcl_scr_syn_anxdep_r","cbcl_scr_syn_withdep_r",
"cbcl_scr_syn_somatic_r","cbcl_scr_syn_social_r",
"cbcl_scr_syn_thought_r","cbcl_scr_syn_rulebreak_r",
"cbcl_scr_syn_aggressive_r"
)
data_hf0b_long <- data_hf0b_long %>%
mutate_at(cbcl_col_numeric, as.numeric) %>%
mutate_at(nback_col_select[-1], as.numeric)
data_hn0b_long <- data_hn0b_long %>%
mutate_at(cbcl_col_numeric, as.numeric) %>%
mutate_at(nback_col_select[-1], as.numeric)
data_fn0b_long <- data_fn0b_long %>%
mutate_at(cbcl_col_numeric, as.numeric) %>%
mutate_at(nback_col_select[-1], as.numeric)
data_0back <- data_0back %>%
mutate_at(cbcl_col_numeric, as.numeric) %>%
mutate_at(nback_col_select[-1], as.numeric)
### 2-back ###
# (1) add a column to specify the condition
data_hf_2b$condition <- "hf"
data_hn_2b$condition <- "hn"
data_fn_2b$condition <- "fn"
# (2) merge distance with family id and site id
data_hf2b <- func_merge_df(data_hf_2b, fam_site_id, nback)
data_hn2b <- func_merge_df(data_hn_2b, fam_site_id, nback)
data_fn2b <- func_merge_df(data_fn_2b, fam_site_id, nback)
# (3) merge with cbcl
data_hf2b <- merge(data_hf2b, abcd_cbcl, by="src_subject_id")
data_hn2b <- merge(data_hn2b, abcd_cbcl, by="src_subject_id")
data_fn2b <- merge(data_fn2b, abcd_cbcl, by="src_subject_id")
# (4) convert the data into long format
colname_list <- c("vis","lim","dor","ven","fro","def","som","v1")
data_hf2b_long <-
pivot_longer(data_hf2b, cols = all_of(colname_list), names_to = "network", values_to = "dis")
data_hn2b_long <-
pivot_longer(data_hn2b, cols = all_of(colname_list), names_to = "network", values_to = "dis")
data_fn2b_long <-
pivot_longer(data_fn2b, cols = all_of(colname_list), names_to = "network", values_to = "dis")
# (5) COMBINE ALL THREE CONDITIONS
data_2back <- rbind(data_hf2b_long, data_hn2b_long, data_fn2b_long)
# (6) make age and CBCL scores numeric
data_hf2b_long <- data_hf2b_long %>%
mutate_at(cbcl_col_numeric, as.numeric) %>%
mutate_at(nback_col_select[-1], as.numeric)
data_hn2b_long <- data_hn2b_long %>%
mutate_at(cbcl_col_numeric, as.numeric) %>%
mutate_at(nback_col_select[-1], as.numeric)
data_fn2b_long <- data_fn2b_long %>%
mutate_at(cbcl_col_numeric, as.numeric) %>%
mutate_at(nback_col_select[-1], as.numeric)
data_2back <- data_2back %>%
mutate_at(cbcl_col_numeric, as.numeric) %>%
mutate_at(nback_col_select[-1], as.numeric)
```
- (2) mixed model - dissimilarity vs. cbcl
one model per network
```{r}
vector_to_factor <- function(data_long) {
data_long$src_subject_id <- as.factor(data_long$src_subject_id)
data_long$rel_family_id <- as.factor(data_long$rel_family_id)
data_long$site_id_l <- as.factor(data_long$site_id_l)
return(data_long)
}
data_0back <- vector_to_factor(data_0back)
data_2back <- vector_to_factor(data_2back)
### total problem score ###
## 0-back ##
# vis
lme2_0back_vis_tol <- lmer(cbcl_scr_syn_totprob_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_tol) # 0.940 ; no interactions
# lim
lme2_0back_lim_tol <- lmer(cbcl_scr_syn_totprob_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_tol) # 0.487 ; no interactions
# v1
lme2_0back_v1_tol <- lmer(cbcl_scr_syn_totprob_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_tol) # 0.0359 * ; no interactions
## 2-back ##
# vis
lme2_2back_vis_tol <- lmer(cbcl_scr_syn_totprob_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_tol) # 0.00016 *** ; no interactions
# lim
lme2_2back_lim_tol <- lmer(cbcl_scr_syn_totprob_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_tol) # 0.477 ; no interactions
# v1
lme2_2back_v1_tol <- lmer(cbcl_scr_syn_totprob_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_tol) # 6.73e-05 *** ; no interactions
### internalization score ###
## 0-back ##
# vis
lme2_0back_vis_int <- lmer(cbcl_scr_syn_internal_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_int) # 0.0145 * ; no interactions
# lim
lme2_0back_lim_int <- lmer(cbcl_scr_syn_internal_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_int) # 0.834 ; no interactions
# v1
lme2_0back_v1_int <- lmer(cbcl_scr_syn_internal_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_int) # 0.473 ; no interactions
## 2-back ##
# vis
lme2_2back_vis_int <- lmer(cbcl_scr_syn_internal_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_int) # 0.325 ; no interactions
# lim
lme2_2back_lim_int <- lmer(cbcl_scr_syn_internal_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_int) # 0.970 ; no interactions
# v1
lme2_2back_v1_int <- lmer(cbcl_scr_syn_internal_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_int) # 0.624 ; no interactions
### externalization score ###
## 0-back ##
# vis
lme2_0back_vis_ext <- lmer(cbcl_scr_syn_external_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_ext) # 0.168 ; no interactions
# lim
lme2_0back_lim_ext <- lmer(cbcl_scr_syn_external_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_ext) # 0.265 ; no interactions
# v1
lme2_0back_v1_ext <- lmer(cbcl_scr_syn_external_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_ext) # 0.0217 * ; no interactions
## 2-back ##
# vis
lme2_2back_vis_ext <- lmer(cbcl_scr_syn_external_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_ext) # 0.001170 ** ; no interactions
# lim
lme2_2back_lim_ext <- lmer(cbcl_scr_syn_external_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_ext) # 0.554 ; no interactions
# v1
lme2_2back_v1_ext <- lmer(cbcl_scr_syn_external_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_ext) # 1.01e-05 *** ; no interactions
```
additional models for CBCL sub-scales
### anxiety/depression score ###
```{r}
## 0-back ##
# lim
lme2_0back_lim_anxdep <- lmer(cbcl_scr_syn_anxdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_anxdep) # 0.9227 ; no interactions
# vis
lme2_0back_vis_anxdep <- lmer(cbcl_scr_syn_anxdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_anxdep) # 0.1643 ; no interactions
# v1
lme2_0back_v1_anxdep <- lmer(cbcl_scr_syn_anxdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_anxdep) # 0.8706 ; no interactions
## 2-back ##
# lim
lme2_2back_lim_anxdep <- lmer(cbcl_scr_syn_anxdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_anxdep) # 0.5283 ; no interactions
# vis
lme2_2back_vis_anxdep <- lmer(cbcl_scr_syn_anxdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_anxdep) # 0.6005 ; no interactions
# v1
lme2_2back_v1_anxdep <- lmer(cbcl_scr_syn_anxdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_anxdep) # 0.7857 ; no interactions
```
### withdrawn/depression score ###
```{r}
## 0-back ##
# lim
lme2_0back_lim_withdep <- lmer(cbcl_scr_syn_withdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_withdep) # 0.37028 ; no interactions
# vis
lme2_0back_vis_withdep <- lmer(cbcl_scr_syn_withdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_withdep) # 0.00141 ** ; no interactions
# v1
lme2_0back_v1_withdep <- lmer(cbcl_scr_syn_withdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_withdep) # 0.02769 *; no interactions
## 2-back ##
# lim
lme2_2back_lim_withdep <- lmer(cbcl_scr_syn_withdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_withdep) # 0.39684 ; no interactions
# vis
lme2_2back_vis_withdep <- lmer(cbcl_scr_syn_withdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_withdep) # 0.64523 ; no interactions
# v1
lme2_2back_v1_withdep <- lmer(cbcl_scr_syn_withdep_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_withdep) # 0.63588 ; no interactions
```
### somatic complaints score ###
```{r}
## 0-back ##
# lim
lme2_0back_lim_somatic <- lmer(cbcl_scr_syn_somatic_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_somatic) # 0.899 ; no interactions
# vis
lme2_0back_vis_somatic <- lmer(cbcl_scr_syn_somatic_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_somatic) # 0.206 ; no interactions
# v1
lme2_0back_v1_somatic <- lmer(cbcl_scr_syn_somatic_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_somatic) # 0.985 ; no interactions
## 2-back ##
# lim
lme2_2back_lim_somatic <- lmer(cbcl_scr_syn_somatic_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_somatic) # 0.977 ; no interactions
# vis
lme2_2back_vis_somatic <- lmer(cbcl_scr_syn_somatic_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_somatic) # 0.157 ; no interactions
# v1
lme2_2back_v1_somatic <- lmer(cbcl_scr_syn_somatic_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_somatic) # 0.722 ; no interactions
```
### rule-breaking behavior score ###
```{r}
## 0-back ##
# lim
lme2_0back_lim_rulebreak <- lmer(cbcl_scr_syn_rulebreak_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_rulebreak) # 0.646 ; no interactions
# vis
lme2_0back_vis_rulebreak <- lmer(cbcl_scr_syn_rulebreak_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_rulebreak) # 0.160 ; no interactions
# v1
lme2_0back_v1_rulebreak <- lmer(cbcl_scr_syn_rulebreak_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_rulebreak) # 0.523 ; no interactions
## 2-back ##
# lim
lme2_2back_lim_rulebreak <- lmer(cbcl_scr_syn_rulebreak_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_rulebreak) # 0.374 ; no interactions
# vis
lme2_2back_vis_rulebreak <- lmer(cbcl_scr_syn_rulebreak_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_rulebreak) # 0.0691 ; no interactions
# v1
lme2_2back_v1_rulebreak <- lmer(cbcl_scr_syn_rulebreak_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_rulebreak) # 0.0189 * ; no interactions
```
### aggressive behavior score ###
```{r}
## 0-back ##
# lim
lme2_0back_lim_aggressive <- lmer(cbcl_scr_syn_aggressive_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_aggressive) # 0.20935 ; no interactions
# vis
lme2_0back_vis_aggressive <- lmer(cbcl_scr_syn_aggressive_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_aggressive) # 0.01586 * ; no interactions
# v1
lme2_0back_v1_aggressive <- lmer(cbcl_scr_syn_aggressive_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_aggressive) # 0.00101 ** ; no interactions
## 2-back ##
# lim
lme2_2back_lim_aggressive <- lmer(cbcl_scr_syn_aggressive_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_aggressive) # 0.68314 ; no interactions
# vis
lme2_2back_vis_aggressive <- lmer(cbcl_scr_syn_aggressive_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_aggressive) # 0.000485 *** ; no interactions
# v1
lme2_2back_v1_aggressive <- lmer(cbcl_scr_syn_aggressive_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_aggressive) # 1.73e-06 *** ; no interactions
```
### social problems score ###
```{r}
## 0-back ##
# lim
lme2_0back_lim_social <- lmer(cbcl_scr_syn_social_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_social) # 0.168 ; no interactions
# vis
lme2_0back_vis_social <- lmer(cbcl_scr_syn_social_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_social) # 0.0688 ; no interactions
# v1
lme2_0back_v1_social <- lmer(cbcl_scr_syn_social_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_social) # 0.00103 ** ; no interactions
## 2-back ##
# lim
lme2_2back_lim_social <- lmer(cbcl_scr_syn_social_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_social) # 0.288 ; no interactions
# vis
lme2_2back_vis_social <- lmer(cbcl_scr_syn_social_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_social) # 3.32e-05 *** ; no interactions
# v1
lme2_2back_v1_social <- lmer(cbcl_scr_syn_social_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_social) # 4.57e-05 *** ; no interactions
```
### thought problems score ###
```{r}
## 0-back ##
# lim
lme2_0back_lim_thought <- lmer(cbcl_scr_syn_thought_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_thought) # 0.70670 ; no interactions
# vis
lme2_0back_vis_thought <- lmer(cbcl_scr_syn_thought_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_thought) # 0.576680 ; no interactions
# v1
lme2_0back_v1_thought <- lmer(cbcl_scr_syn_thought_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_thought) # 0.043311 * ; no interactions
## 2-back ##
# lim
lme2_2back_lim_thought <- lmer(cbcl_scr_syn_thought_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_thought) # 0.524190 ; no interactions
# vis
lme2_2back_vis_thought <- lmer(cbcl_scr_syn_thought_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_thought) # 0.00191 ** ; no interactions
# v1
lme2_2back_v1_thought <- lmer(cbcl_scr_syn_thought_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_thought) # 0.00732 ** ; no interactions
```
### attention problems score ###
```{r}
## 0-back ##
# lim
lme2_0back_lim_attention <- lmer(cbcl_scr_syn_attention_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lme2_0back_lim_attention) # 0.87337 ; no interactions
# vis
lme2_0back_vis_attention <- lmer(cbcl_scr_syn_attention_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lme2_0back_vis_attention) # 0.1514 ; no interactions
# v1
lme2_0back_v1_attention <- lmer(cbcl_scr_syn_attention_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lme2_0back_v1_attention) # 0.000538 *** ; no interactions
## 2-back ##
# lim
lme2_2back_lim_attention <- lmer(cbcl_scr_syn_attention_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lme2_2back_lim_attention) # 0.2165 ; no interactions
# vis
lme2_2back_vis_attention <- lmer(cbcl_scr_syn_attention_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lme2_2back_vis_attention) # 4.2e-05 *** ; no interactions
# v1
lme2_2back_v1_attention <- lmer(cbcl_scr_syn_attention_r ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lme2_2back_v1_attention) # 0.000116 *** ; no interactions
```
- (3) plot - dissimilarity vs. cbcl
(3.1) cbcl score distributions
```{r}
# preprocess dataset
cbcl_col_list <-
c(
"cbcl_scr_syn_attention_t", "cbcl_scr_syn_attention_r",
"cbcl_scr_syn_external_t", "cbcl_scr_syn_external_r",
"cbcl_scr_syn_internal_t", "cbcl_scr_syn_internal_r",
"cbcl_scr_syn_totprob_t", "cbcl_scr_syn_totprob_r",
"cbcl_scr_syn_anxdep_r","cbcl_scr_syn_withdep_r",
"cbcl_scr_syn_somatic_r","cbcl_scr_syn_social_r",
"cbcl_scr_syn_thought_r","cbcl_scr_syn_rulebreak_r",
"cbcl_scr_syn_aggressive_r"
)
data_cbcl <- data_hf0b %>%
mutate_at(cbcl_col_list, as.numeric)
f1_totptob_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_totprob_r)) +
geom_histogram(binwidth=2, color = "black", fill = "white") +
scale_y_continuous(expand = c(0,0), limits = c(0,600)) +
labs(x = "total problems raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_ext_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_external_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
scale_y_continuous(expand = c(0,0), limits = c(0,1700)) +
labs(x = "externalizing problems raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_int_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_internal_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
scale_y_continuous(expand = c(0,0), limits = c(0,1000)) +
labs(x = "internalizing problems raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_totptob_t <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_totprob_t)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
scale_x_continuous(expand = c(0,0), limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,350)) +
labs(x = "total problems t-score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_ext_t <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_external_t)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
scale_x_continuous(expand = c(0,0), limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,1300)) +
labs(x = "externalizing problems t-score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_int_t <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_internal_t)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
scale_x_continuous(expand = c(0,0),limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,600)) +
labs(x = "internalizing problems t-score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_totptob_r
f1_ext_r
f1_int_r
f1_totptob_t
f1_ext_t
f1_int_t
# sanity check
# min(data_0back["cbcl_scr_syn_totprob_r"]) # 0
# max(data_0back["cbcl_scr_syn_totprob_r"]) # 161
# min(data_0back["cbcl_scr_syn_internal_r"]) # 0
# max(data_0back["cbcl_scr_syn_internal_r"]) # 50
# min(data_0back["cbcl_scr_syn_external_r"]) # 0
# max(data_0back["cbcl_scr_syn_external_r"]) # 50
# min(data_2back["cbcl_scr_syn_totprob_r"]) # 0
# max(data_2back["cbcl_scr_syn_totprob_r"]) # 161
# min(data_2back["cbcl_scr_syn_internal_r"]) # 0
# max(data_2back["cbcl_scr_syn_internal_r"]) # 50
# min(data_2back["cbcl_scr_syn_external_r"]) # 0
# max(data_2back["cbcl_scr_syn_external_r"]) # 50
# min(data_2back["cbcl_scr_syn_totprob_t"]) # 24
# max(data_2back["cbcl_scr_syn_totprob_t"]) # 88
# min(data_2back["cbcl_scr_syn_internal_t"]) # 33
# max(data_2back["cbcl_scr_syn_internal_t"]) # 90
# min(data_2back["cbcl_scr_syn_external_t"]) # 33
# max(data_2back["cbcl_scr_syn_external_t"]) # 83
# min(data_0back["dis"]) # 0.00586022
# max(data_0back["dis"]) # 1.757153
# min(data_2back["dis"]) # 0.006938755
# max(data_2back["dis"]) # 1.900827
# head(data_0back)
```
```{r}
f1_anxdep_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_anxdep_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
# scale_x_continuous(expand = c(0,0), limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,2000)) +
labs(x = "anxious/depressed raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_anxdep_r
```
```{r}
f1_withdep_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_withdep_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
# scale_x_continuous(expand = c(0,0), limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,3000)) +
labs(x = "withdrawn/depressed raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_withdep_r
```
```{r}
f1_somatic_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_somatic_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
# scale_x_continuous(expand = c(0,0), limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,2500)) +
labs(x = "somatic complaints raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_somatic_r
```
```{r}
f1_social_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_social_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
# scale_x_continuous(expand = c(0,0), limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,3000)) +
labs(x = "social problems raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_social_r
```
```{r}
f1_thought_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_thought_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
# scale_x_continuous(expand = c(0,0), limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,2500)) +
labs(x = "thought problems raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_thought_r
```
```{r}
f1_rulebreak_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_rulebreak_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
# scale_x_continuous(expand = c(0,0), limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,3000)) +
labs(x = "rule-breaking behavior raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_rulebreak_r
```
```{r}
f1_aggressive_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_aggressive_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
# scale_x_continuous(expand = c(0,0), limits = c(20,95)) +
scale_y_continuous(expand = c(0,0), limits = c(0,2000)) +
labs(x = "aggressive behavior raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_aggressive_r
```
```{r}
f1_attention_r <- ggplot(data_cbcl, aes(x=cbcl_scr_syn_attention_r)) +
geom_histogram(binwidth=1, color = "black", fill = "white") +
#scale_x_continuous(expand = c(0,0), limits = c(-1,21)) +
scale_y_continuous(expand = c(0,0), limits = c(0,2000)) +
labs(x = "attention problems raw score",
y = "number of subjects") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18), aspect.ratio = 1)
f1_attention_r
```
(3.2) subjects with both ext and int problems (check overlaps)
```{r}
# head(data_cbcl)
sum((data_cbcl$cbcl_scr_syn_external_t>=60) & (data_cbcl$cbcl_scr_syn_internal_t>=60)) # 204
length(data_cbcl$src_subject_id[(data_cbcl$cbcl_scr_syn_external_t>=60) & (data_cbcl$cbcl_scr_syn_internal_t>=60)]) # 204
# data_cbcl$src_subject_id[(data_cbcl$cbcl_scr_syn_external_t>=60) & (data_cbcl$cbcl_scr_syn_internal_t>=60)]
```
(3.3) total problems score plot
```{r}
# cbcl_scr_syn_totprob_r
### 0-back ###
# vis
f2_0back_vis_tol <-
ggplot(data_0back[data_0back$network == "vis",], aes(cbcl_scr_syn_totprob_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="total problems score", y="dissimilarities", title="Vis", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-10,180)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# lim
f2_0back_lim_tol <-
ggplot(data_0back[data_0back$network == "lim",], aes(cbcl_scr_syn_totprob_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="total problems score", y="dissimilarities", title="Lim", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-10,180)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# v1
f2_0back_v1_tol <-
ggplot(data_0back[data_0back$network == "v1",], aes(cbcl_scr_syn_totprob_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="total problems score", y="dissimilarities", title="V1", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-10,180)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
### 2-back ###
# vis
f2_2back_vis_tol <-
ggplot(data_2back[data_2back$network == "vis",], aes(cbcl_scr_syn_totprob_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="total problems score", y="dissimilarities", title="Vis", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-10,180)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# lim
f2_2back_lim_tol <-
ggplot(data_2back[data_2back$network == "lim",], aes(cbcl_scr_syn_totprob_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="total problems score", y="dissimilarities", title="Lim", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-10,180)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# v1
f2_2back_v1_tol <-
ggplot(data_2back[data_2back$network == "v1",], aes(cbcl_scr_syn_totprob_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="total problems score", y="dissimilarities", title="V1", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-10,180)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
f2_0back_vis_tol
f2_0back_lim_tol
f2_0back_v1_tol
f2_2back_vis_tol
f2_2back_lim_tol
f2_2back_v1_tol
```
(3.4) internalizating score plot
```{r}
# cbcl_scr_syn_internal_r
### 0-back ###
# vis
f2_0back_vis_int <-
ggplot(data_0back[data_0back$network == "vis",], aes(cbcl_scr_syn_internal_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="internalizing problems score", y="dissimilarities", title="Vis", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,60)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# lim
f2_0back_lim_int <-
ggplot(data_0back[data_0back$network == "lim",], aes(cbcl_scr_syn_internal_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="internalizing problems score", y="dissimilarities", title="Lim", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,60)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# v1
f2_0back_v1_int <-
ggplot(data_0back[data_0back$network == "v1",], aes(cbcl_scr_syn_internal_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="internalizing problems score", y="dissimilarities", title="V1", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
### 2-back ###
# vis
f2_2back_vis_int <-
ggplot(data_2back[data_2back$network == "vis",], aes(cbcl_scr_syn_internal_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="internalizing problems score", y="dissimilarities", title="Vis", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# lim
f2_2back_lim_int <-
ggplot(data_2back[data_2back$network == "lim",], aes(cbcl_scr_syn_internal_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="internalizing problems score", y="dissimilarities", title="Lim", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# v1
f2_2back_v1_int <-
ggplot(data_2back[data_2back$network == "v1",], aes(cbcl_scr_syn_internal_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="internalizing problems score", y="dissimilarities", title="V1", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
f2_0back_vis_int
f2_0back_lim_int
f2_0back_v1_int
f2_2back_vis_int
f2_2back_lim_int
f2_2back_v1_int
```
(3.5) externalizating score plot
```{r}
# cbcl_scr_syn_external_r
### 0-back ###
# vis
f2_0back_vis_ext <-
ggplot(data_0back[data_0back$network == "vis",], aes(cbcl_scr_syn_external_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="externalizing problems score", y="dissimilarities", title="Vis", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# lim
f2_0back_lim_ext <-
ggplot(data_0back[data_0back$network == "lim",], aes(cbcl_scr_syn_external_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="externalizing problems score", y="dissimilarities", title="Lim", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# v1
f2_0back_v1_ext <-
ggplot(data_0back[data_0back$network == "v1",], aes(cbcl_scr_syn_external_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="externalizing problems score", y="dissimilarities", title="V1", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
### 2-back ###
# vis
f2_2back_vis_ext <-
ggplot(data_2back[data_2back$network == "vis",], aes(cbcl_scr_syn_external_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="externalizing problems score", y="dissimilarities", title="Vis", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# lim
f2_2back_lim_ext <-
ggplot(data_2back[data_2back$network == "lim",], aes(cbcl_scr_syn_external_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="externalizing problems score", y="dissimilarities", title="Lim", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
# v1
f2_2back_v1_ext <-
ggplot(data_2back[data_2back$network == "v1",], aes(cbcl_scr_syn_external_r, dis)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(x="externalizing problems score", y="dissimilarities", title="V1", fill="count") +
scale_x_continuous(expand = c(0,0), limits = c(-2,55)) +
scale_y_continuous(expand = c(0,0), limits = c(-0.1,2.1)) +
theme_classic() +
theme(text = element_text(size = 20),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 18))
f2_0back_vis_ext
f2_0back_lim_ext
f2_0back_v1_ext
f2_2back_vis_ext
f2_2back_lim_ext
f2_2back_v1_ext
```
- (4) predicting CBCL scores
process dataframes
```{r}
# create a data frame for CBCL prediction
data_0back_3y_cbcl <- merge(data_0back, abcd_cbcl_3year_totprob_r, by="src_subject_id")
data_2back_3y_cbcl <- merge(data_2back, abcd_cbcl_3year_totprob_r, by="src_subject_id")
# keep complete data
data_0back_3y_cbcl <- data_0back_3y_cbcl[!is.na(data_0back_3y_cbcl$cbcl_scr_syn_totprob_r_3year),]
data_2back_3y_cbcl <- data_2back_3y_cbcl[!is.na(data_2back_3y_cbcl$cbcl_scr_syn_totprob_r_3year),]
# sanity check
dim(data_0back_3y_cbcl) # 94032 27
dim(data_2back_3y_cbcl) # 94032 27
head(data_0back_3y_cbcl)
head(data_2back_3y_cbcl)
```
mixed models - dissimilarity vs. future cbcl
```{r}
### future total problem score ###
## 0-back ##
# vis
lme2_0back_vis_tol_3yr <- lmer(cbcl_scr_syn_totprob_r_3year ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back_3y_cbcl[data_0back_3y_cbcl$network == "vis",])
summary(lme2_0back_vis_tol_3yr) # 0.4407 ; no interactions
# lim
lme2_0back_lim_tol_3yr <- lmer(cbcl_scr_syn_totprob_r_3year ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back_3y_cbcl[data_0back_3y_cbcl$network == "lim",])
summary(lme2_0back_lim_tol_3yr) # 0.9499 ; no interactions
# v1
lme2_0back_v1_tol_3yr <- lmer(cbcl_scr_syn_totprob_r_3year ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back_3y_cbcl[data_0back_3y_cbcl$network == "v1",])
summary(lme2_0back_v1_tol_3yr) # 0.000457 *** ; no interactions
## 2-back ##
# vis
lme2_2back_vis_tol_3yr <- lmer(cbcl_scr_syn_totprob_r_3year ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back_3y_cbcl[data_2back_3y_cbcl$network == "vis",])
summary(lme2_2back_vis_tol_3yr) # 0.000525 *** ; no interactions
# lim
lme2_2back_lim_tol_3yr <- lmer(cbcl_scr_syn_totprob_r_3year ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back_3y_cbcl[data_2back_3y_cbcl$network == "lim",])
summary(lme2_2back_lim_tol_3yr) # 0.8226 ; no interactions
# v1
lme2_2back_v1_tol_3yr <- lmer(cbcl_scr_syn_totprob_r_3year ~ dis * condition + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back_3y_cbcl[data_2back_3y_cbcl$network == "v1",])
summary(lme2_2back_v1_tol_3yr) # 1.04e-06 *** ; no interactions
```
- (5) group analysis (using > 63 as the clinical threshold)
process dataset
```{r}
# group match with age, sex, & site id
low_score = 50
cut_off_score = 63
### for 0-back ###
# subsetting data
data_0back_cbclgroup <- data_0back[(data_0back$network == "vis" & data_0back$condition == "hf"),]
# check number of subjects
sum(data_0back_cbclgroup$cbcl_scr_syn_totprob_t <= low_score) # 3573
sum(data_0back_cbclgroup$cbcl_scr_syn_totprob_t > cut_off_score) # 264
# add a column to label the cbcl group (0 - low, 1 - high)
data_0back_cbclgroup$cbcl_group_clinical <-
ifelse(data_0back_cbclgroup$cbcl_scr_syn_totprob_t <= low_score, 0, # is less than or equals to low_score, set as 0
ifelse(data_0back_cbclgroup$cbcl_scr_syn_totprob_t > cut_off_score, 1, NA)) # is equals to or more than cut_off_score, set as 1; NA otherwise
# data cleaning
# save the subset data of only low and high groups
data_0back_cbclgroup <- data_0back_cbclgroup[!is.na(data_0back_cbclgroup$cbcl_group_clinical),]
data_0back_cbclgroup$site_id_l <- as.character(data_0back_cbclgroup$site_id_l)
# data matching
matchit_model <- matchit(
cbcl_group_clinical ~ interview_age + sex + site_id_l,
data = data_0back_cbclgroup,
method = "nearest", # nearest neighbor matching
distance = "logit", # propensity score via logistic regression
ratio = 1, # 1:1 matching
replace = FALSE # without replacement
)
summary(matchit_model)
# get matched dataframe
df_match_0back <- match_data(matchit_model,
data = data_0back_cbclgroup,
distance = "prop.score")
# get the list of the subject ids
matched_subject_id <- df_match_0back$src_subject_id
# get the matched data
data_0back_matched <- data_0back[data_0back$src_subject_id %in% matched_subject_id,]
data_2back_matched <- data_2back[data_2back$src_subject_id %in% matched_subject_id,]
# add a column to label the cbcl group (0 - low, 1 - high)
data_0back_matched$cbcl_group_clinical <-
ifelse(data_0back_matched$cbcl_scr_syn_totprob_t <= low_score, 0,
ifelse(data_0back_matched$cbcl_scr_syn_totprob_t > cut_off_score, 1, NA))
data_2back_matched$cbcl_group_clinical <-
ifelse(data_2back_matched$cbcl_scr_syn_totprob_t <= low_score, 0,
ifelse(data_2back_matched$cbcl_scr_syn_totprob_t > cut_off_score, 1, NA))
# sanity check
dim(data_0back_matched) # 12672 27 (264*2*8*3
dim(data_2back_matched) # 12672 27
head(data_0back_matched)
head(data_2back_matched)
```
compare dissimilarities between groups
```{r}
# t-tests
# t-test of all dissimilarities (all emotion comparisons)
### 0-back ###
ttest_vis_0b <- t.test(dis ~ cbcl_group_clinical, data = data_0back_matched[data_0back_matched$network == "vis",])
ttest_lim_0b <- t.test(dis ~ cbcl_group_clinical, data = data_0back_matched[data_0back_matched$network == "lim",])
ttest_v1_0b <- t.test(dis ~ cbcl_group_clinical, data = data_0back_matched[data_0back_matched$network == "v1",])
ttest_vis_0b # 0.001511
ttest_lim_0b # 0.05925
ttest_v1_0b # 0.004075
### 2-back ###
ttest_vis_2b <- t.test(dis ~ cbcl_group_clinical, data = data_2back_matched[data_2back_matched$network == "vis",])
ttest_lim_2b <- t.test(dis ~ cbcl_group_clinical, data = data_2back_matched[data_2back_matched$network == "lim",])
ttest_v1_2b <- t.test(dis ~ cbcl_group_clinical, data = data_2back_matched[data_2back_matched$network == "v1",])
ttest_vis_2b # 3.222e-06
ttest_lim_2b # 0.0005018
ttest_v1_2b # 0.004273
# ttest_vis_2b$statistic
# ttest_vis_2b$p.value
# effect sizes
### 0-back ###
# a function to get a sub-dataframe for low group
func_low_group_0b <- function(network_name) {
group1 <- data_0back_matched[(data_0back_matched$network == network_name & data_0back_matched$cbcl_group_clinical == 0),]
return(group1$dis)
}
# a function to get a sub-dataframe for high group
func_high_group_0b <- function(network_name) {
group2 <- data_0back_matched[(data_0back_matched$network == network_name & data_0back_matched$cbcl_group_clinical == 1),]
return(group2$dis)
}
cohen.d(func_low_group_0b("vis"), func_high_group_0b("vis")) # -0.1597133 (negligible)
cohen.d(func_low_group_0b("lim"), func_high_group_0b("lim")) # -0.09485952 (negligible)
cohen.d(func_low_group_0b("v1"), func_high_group_0b("v1")) # -0.1445535 (negligible)
### 2-back ###
# a function to get a sub-dataframe for low group
func_low_group_2b <- function(network_name) {
group1 <- data_2back_matched[(data_2back_matched$network == network_name & data_2back_matched$cbcl_group_clinical == 0),]
return(group1$dis)
}
# a function to get a sub-dataframe for high group
func_high_group_2b <- function(network_name) {
group2 <- data_2back_matched[(data_2back_matched$network == network_name & data_2back_matched$cbcl_group_clinical == 1),]
return(group2$dis)
}
cohen.d(func_low_group_2b("vis"), func_high_group_2b("vis")) # -0.2348411 (small)
cohen.d(func_low_group_2b("lim"), func_high_group_2b("lim")) # -0.1752295 (negligible)
cohen.d(func_low_group_2b("v1"), func_high_group_2b("v1")) # -0.1438009 (negligible)
```
# part 3 - dissimilarity and performance
- (1) process datasets
```{r}
# Calculate the mean accuracy of all faces
# 0-back
data_0back <- data_0back %>%
mutate(data_0back, allface_0b_rate = rowMeans(select(data_0back, c("tfmri_nb_all_beh_c0bpf_rate","tfmri_nb_all_beh_c0bnf_rate", "tfmri_nb_all_beh_c0bngf_rate"))))
# 2-back
data_2back <- data_2back %>%
mutate(data_2back, allface_2b_rate = rowMeans(select(data_2back, c("tfmri_nb_all_beh_c2bpf_rate","tfmri_nb_all_beh_c2bnf_rate","tfmri_nb_all_beh_c2bngf_rate"))))
# head(data_0back)
# head(data_2back)
```
- (2) mixed models
```{r}
###### controlled for overall performance ######
## 0-back ##
# vis #
lmer3_0back_vis <- lmer(allface_0b_rate ~ dis + tfmri_nb_all_beh_c0b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lmer3_0back_vis) # 0.8330
# lim #
lmer3_0back_lim <- lmer(allface_0b_rate ~ dis + tfmri_nb_all_beh_c0b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lmer3_0back_lim) # 0.1511
# v1
lmer3_0back_v1 <- lmer(allface_0b_rate ~ dis + tfmri_nb_all_beh_c0b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lmer3_0back_v1) # 0.391
## 2-back ##
# vis #
lmer3_2back_vis <- lmer(allface_2b_rate ~ dis + tfmri_nb_all_beh_c2b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lmer3_2back_vis) # 0.000977 ***
# lim
lmer3_2back_lim <- lmer(allface_2b_rate ~ dis + tfmri_nb_all_beh_c2b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lmer3_2back_lim) # 0.00492 **
# v1
lmer3_2back_v1 <- lmer(allface_2b_rate ~ dis + tfmri_nb_all_beh_c2b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lmer3_2back_v1) # 0.0215 *
```
- (3) plots
```{r}
# scatter plot of all three dissimilarity and cbcl scores
## 0-back
# vis
f3_0back_vis <-
ggplot(data_0back[data_0back$network == "vis",], aes(x=dis, y=allface_0b_rate)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(title = "Vis (0-back)",
x = "Dissimilarity",
y = "Accuracy") +
scale_x_continuous(expand = c(0,0), limits = c(0,2)) +
scale_y_continuous(expand = c(0,0), limits = c(0.5,1.1)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18))
# lim
f3_0back_lim <-
ggplot(data_0back[data_0back$network == "lim",], aes(x=dis, y=allface_0b_rate)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(title = "Lim (0-back)",
x = "Dissimilarity",
y = "Accuracy") +
scale_x_continuous(expand = c(0,0), limits = c(0,2)) +
scale_y_continuous(expand = c(0,0), limits = c(0.5,1.1)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18))
# v1
f3_0back_v1 <-
ggplot(data_0back[data_0back$network == "v1",], aes(x=dis, y=allface_0b_rate)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(title = "V1 (0-back)",
x = "Dissimilarity",
y = "Accuracy") +
scale_x_continuous(expand = c(0,0), limits = c(0,2)) +
scale_y_continuous(expand = c(0,0), limits = c(0.5,1.1)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18))
## 2-back
# vis
f3_2back_vis <-
ggplot(data_2back[data_2back$network == "vis",], aes(x=dis, y=allface_2b_rate)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(title = "Vis (2-back)",
x = "Dissimilarity",
y = "Accuracy") +
scale_x_continuous(expand = c(0,0), limits = c(0,2)) +
scale_y_continuous(expand = c(0,0), limits = c(0.5,1.1)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18))
# lim
f3_2back_lim <-
ggplot(data_2back[data_2back$network == "lim",], aes(x=dis, y=allface_2b_rate)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(title = "Lim (2-back)",
x = "Dissimilarity",
y = "Accuracy") +
scale_x_continuous(expand = c(0,0), limits = c(0,2)) +
scale_y_continuous(expand = c(0,0), limits = c(0.5,1.1)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18))
# v1
f3_2back_v1 <-
ggplot(data_2back[data_2back$network == "v1",], aes(x=dis, y=allface_2b_rate)) +
geom_hex() + stat_smooth(color="red", fill="red", method=lm, size=0.5, level=0.999) +
labs(title = "V1 (2-back)",
x = "Dissimilarity",
y = "Accuracy") +
scale_x_continuous(expand = c(0,0), limits = c(0,2)) +
scale_y_continuous(expand = c(0,0), limits = c(0.5,1.1)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 18))
f3_0back_vis
f3_0back_lim
f3_0back_v1
f3_2back_vis
f3_2back_lim
f3_2back_v1
```
- (4) mixed models with three levels of emotion comparisons
```{r}
###### controlled for overall performance ######
## 0-back ##
# vis #
lmer3_0back_vis <- lmer(allface_0b_rate ~ dis * condition + tfmri_nb_all_beh_c0b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "vis",])
summary(lmer3_0back_vis) # 0.7076
# lim #
lmer3_0back_lim <- lmer(allface_0b_rate ~ dis * condition + tfmri_nb_all_beh_c0b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "lim",])
summary(lmer3_0back_lim) # 0.1559
# v1
lmer3_0back_v1 <- lmer(allface_0b_rate ~ dis * condition + tfmri_nb_all_beh_c0b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_0back[data_0back$network == "v1",])
summary(lmer3_0back_v1) # 0.4662
## 2-back ##
# vis #
lmer3_2back_vis <- lmer(allface_2b_rate ~ dis * condition + tfmri_nb_all_beh_c2b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "vis",])
summary(lmer3_2back_vis) # 0.019417 *
# lim
lmer3_2back_lim <- lmer(allface_2b_rate ~ dis * condition + tfmri_nb_all_beh_c2b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "lim",])
summary(lmer3_2back_lim) # 0.0878 .
# v1
lmer3_2back_v1 <- lmer(allface_2b_rate ~ dis * condition + tfmri_nb_all_beh_c2b_rate + sex + interview_age + (1|site_id_l/rel_family_id), data = data_2back[data_2back$network == "v1",])
summary(lmer3_2back_v1) # 0.157
```
# part 4 - effect sizes
- (1) process dataset for effect sizes
```{r}
# merge the dataframes
data_hf_0b_eff <- merge(data_hf_0b, fam_site_id, by="src_subject_id")
data_hn_0b_eff <- merge(data_hn_0b, fam_site_id, by="src_subject_id")
data_fn_0b_eff <- merge(data_fn_0b, fam_site_id, by="src_subject_id")
data_hf_2b_eff <- merge(data_hf_2b, fam_site_id, by="src_subject_id")
data_hn_2b_eff <- merge(data_hn_2b, fam_site_id, by="src_subject_id")
data_fn_2b_eff <- merge(data_fn_2b, fam_site_id, by="src_subject_id")
data_hf_0b_eff$nback <- "0_back"
data_hn_0b_eff$nback <- "0_back"
data_fn_0b_eff$nback <- "0_back"
data_hf_2b_eff$nback <- "2_back"
data_hn_2b_eff$nback <- "2_back"
data_fn_2b_eff$nback <- "2_back"
data_hf_0b_eff$condition <- "hf"
data_hn_0b_eff$condition <- "hn"
data_fn_0b_eff$condition <- "fn"
data_hf_2b_eff$condition <- "hf"
data_hn_2b_eff$condition <- "hn"
data_fn_2b_eff$condition <- "fn"
data_hf_0b_eff$v1 <- NULL
data_hn_0b_eff$v1 <- NULL
data_fn_0b_eff$v1 <- NULL
data_hf_2b_eff$v1 <- NULL
data_hn_2b_eff$v1 <- NULL
data_fn_2b_eff$v1 <- NULL
# (6.1) sanity check
dim(data_hf_0b_eff) # 4952 14
dim(data_hf_2b_eff) # 4952 14
dim(data_hn_0b_eff) # 4952 14
dim(data_hn_2b_eff) # 4952 14
dim(data_fn_0b_eff) # 4952 14
dim(data_fn_2b_eff) # 4952 14
# head(data_hf_0b_eff)
```
- (2) Process dataframe for mixed effect models
one model for each network - association between 1-r and each network (for example, vis vs. non-vis)
```{r}
list_eff <- c("vis","lim","dor","ven","fro","def","som")
### 0-back ###
# convert the data into long format
data_hf_0b_long <- pivot_longer(data_hf_0b_eff, cols = all_of(list_eff), names_to = "network", values_to = "dis")
data_hn_0b_long <- pivot_longer(data_hn_0b_eff, cols = all_of(list_eff), names_to = "network", values_to = "dis")
data_fn_0b_long <- pivot_longer(data_fn_0b_eff, cols = all_of(list_eff), names_to = "network", values_to = "dis")
# COMBINE ALL THREE CONDITIONS
data_0b <- rbind(data_hf_0b_long, data_hn_0b_long, data_fn_0b_long)
# convert vector objects to factors
data_0b <- vector_to_factor(data_0b)
# setup the data for mixed effect model
data_0b$vis <- (data_0b$network == "vis") * 1
data_0b$lim <- (data_0b$network == "lim") * 1
data_0b$dor <- (data_0b$network == "dor") * 1
data_0b$ven <- (data_0b$network == "ven") * 1
data_0b$fro <- (data_0b$network == "fro") * 1
data_0b$def <- (data_0b$network == "def") * 1
data_0b$som <- (data_0b$network == "som") * 1
data_0b$vis <- as.factor(data_0b$vis)
data_0b$lim <- as.factor(data_0b$lim)
data_0b$dor <- as.factor(data_0b$dor)
data_0b$ven <- as.factor(data_0b$ven)
data_0b$fro <- as.factor(data_0b$fro)
data_0b$def <- as.factor(data_0b$def)
data_0b$som <- as.factor(data_0b$som)
### 2-back ###
# convert the data into long format
data_hf_2b_long <- pivot_longer(data_hf_2b_eff, cols = all_of(list_eff), names_to = "network", values_to = "dis")
data_hn_2b_long <- pivot_longer(data_hn_2b_eff, cols = all_of(list_eff), names_to = "network", values_to = "dis")
data_fn_2b_long <- pivot_longer(data_fn_2b_eff, cols = all_of(list_eff), names_to = "network", values_to = "dis")
# COMBINE ALL THREE CONDITIONS
data_2b <- rbind(data_hf_2b_long, data_hn_2b_long, data_fn_2b_long)
# convert vector objects to factors
data_2b <- vector_to_factor(data_2b)
# setup the data for mixed effect model
data_2b$vis <- (data_2b$network == "vis") * 1
data_2b$lim <- (data_2b$network == "lim") * 1
data_2b$dor <- (data_2b$network == "dor") * 1
data_2b$ven <- (data_2b$network == "ven") * 1
data_2b$fro <- (data_2b$network == "fro") * 1
data_2b$def <- (data_2b$network == "def") * 1
data_2b$som <- (data_2b$network == "som") * 1
data_2b$vis <- as.factor(data_2b$vis)
data_2b$lim <- as.factor(data_2b$lim)
data_2b$dor <- as.factor(data_2b$dor)
data_2b$ven <- as.factor(data_2b$ven)
data_2b$fro <- as.factor(data_2b$fro)
data_2b$def <- as.factor(data_2b$def)
data_2b$som <- as.factor(data_2b$som)
```
- (3) run mixed models
use R marginal squared as effect sizes of each network
```{r}
### 0-back ###
# run mixed effect model for seven networks separately
lme1_vis_0b <- lmer(dis ~ vis * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_0b)
lme1_lim_0b <- lmer(dis ~ lim * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_0b)
lme1_dor_0b <- lmer(dis ~ dor * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_0b)
lme1_ven_0b <- lmer(dis ~ ven * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_0b)
lme1_fro_0b <- lmer(dis ~ fro * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_0b)
lme1_def_0b <- lmer(dis ~ def * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_0b)
lme1_som_0b <- lmer(dis ~ som * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_0b)
# get the R squared marginal values
m <- matrix(ncol = 0, nrow = 0) # create an empty matrix
r_squared_0back <- data.frame(m) # empty matrix to dataframe
r_squared_0back <- r.squaredGLMM(lme1_vis_0b) %>%
rbind(r_squared_0back, r.squaredGLMM(lme1_lim_0b)) %>%
rbind(r_squared_0back, r.squaredGLMM(lme1_dor_0b)) %>%
rbind(r_squared_0back, r.squaredGLMM(lme1_ven_0b)) %>%
rbind(r_squared_0back, r.squaredGLMM(lme1_fro_0b)) %>%
rbind(r_squared_0back, r.squaredGLMM(lme1_def_0b)) %>%
rbind(r_squared_0back, r.squaredGLMM(lme1_som_0b))
row.names(r_squared_0back) <- list_eff
### 2-back ###
# run mixed effect model for seven networks separately
lme1_vis_2b <- lmer(dis ~ vis * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_2b)
lme1_lim_2b <- lmer(dis ~ lim * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_2b)
lme1_dor_2b <- lmer(dis ~ dor * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_2b)
lme1_ven_2b <- lmer(dis ~ ven * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_2b)
lme1_fro_2b <- lmer(dis ~ fro * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_2b)
lme1_def_2b <- lmer(dis ~ def * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_2b)
lme1_som_2b <- lmer(dis ~ som * condition + (1|src_subject_id) + (1|site_id_l/rel_family_id), data = data_2b)
# get the R squared marginal values
m <- matrix(ncol = 0, nrow = 0) # create an empty matrix
r_squared_2back <- data.frame(m) # empty matrix to dataframe
r_squared_2back <- r.squaredGLMM(lme1_vis_2b) %>%
rbind(r_squared_2back, r.squaredGLMM(lme1_lim_2b)) %>%
rbind(r_squared_2back, r.squaredGLMM(lme1_dor_2b)) %>%
rbind(r_squared_2back, r.squaredGLMM(lme1_ven_2b)) %>%
rbind(r_squared_2back, r.squaredGLMM(lme1_fro_2b)) %>%
rbind(r_squared_2back, r.squaredGLMM(lme1_def_2b)) %>%
rbind(r_squared_2back, r.squaredGLMM(lme1_som_2b))
row.names(r_squared_2back) <- list_eff
r_squared_0back
r_squared_2back
```
- (4) plot - effect sizes bar plots
```{r}
df_plot = data.frame(network=rownames(r_squared_0back),
r2m_0back=as.numeric(r_squared_0back[-1][,1]),
r2m_2back=as.numeric(r_squared_2back[-1][,1]))
# list_eff <- c("vis","lim","dor","ven","fro","def","som") # defined before
plot_order_of_networks <- c("dor","ven","fro","lim","def","vis","som")
# 0-back
f1_0back <- ggplot(data=df_plot, aes(x = network, y = r2m_0back)) +
geom_bar(stat = "identity", position = "dodge", color = "black", fill = "white") +
aes(x = factor(network, plot_order_of_networks), y = r2m_0back) +
labs(x = "Functional networks",
y = "Effect sizes") +
scale_y_continuous(expand = c(0,0), limits = c(0,0.6)) +
labs(title = "0-back") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 20), axis.text.x = element_text(angle = 90, hjust = 1), aspect.ratio = 1)
f1_0back
# 2-back
f1_2back <- ggplot(data=df_plot, aes(x = network, y = r2m_2back)) +
geom_bar(stat = "identity", position = "dodge", color = "black", fill = "white") +
aes(x = factor(network, plot_order_of_networks), y = r2m_2back) +
labs(x = "Functional networks",
y = "Effect sizes") +
scale_y_continuous(expand = c(0,0), limits = c(0,0.6)) +
labs(title = "2-back") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), text = element_text(size = 20), axis.text.x = element_text(angle = 90, hjust = 1), aspect.ratio = 1)
f1_2back
```
- (5) cohen.d
```{r}
# 0-back
cohen.d(data_0back$dis[data_0back$network == "vis"], data_0back$dis[!(data_0back$network == "vis")]) # -0.7819427 (medium)
cohen.d(data_0back$dis[data_0back$network == "lim"], data_0back$dis[!(data_0back$network == "lim")]) # 1.161041 (large)
cohen.d(data_0back$dis[data_0back$network == "fro"], data_0back$dis[!(data_0back$network == "fro")]) # 0.1922416 (negligible)
cohen.d(data_0back$dis[data_0back$network == "dor"], data_0back$dis[!(data_0back$network == "dor")]) # -0.1743875 (negligible)
cohen.d(data_0back$dis[data_0back$network == "ven"], data_0back$dis[!(data_0back$network == "ven")]) # 0.2764237 (small)
cohen.d(data_0back$dis[data_0back$network == "def"], data_0back$dis[!(data_0back$network == "def")]) # 0.4370847 (small)
cohen.d(data_0back$dis[data_0back$network == "som"], data_0back$dis[!(data_0back$network == "som")]) # 0.4547063 (small)
# 2-back
cohen.d(data_2back$dis[data_2back$network == "vis"], data_2back$dis[!(data_2back$network == "vis")]) # -0.5052931 (medium)
cohen.d(data_2back$dis[data_2back$network == "lim"], data_2back$dis[!(data_2back$network == "lim")]) # 1.417457 (large)
cohen.d(data_2back$dis[data_2back$network == "fro"], data_2back$dis[!(data_2back$network == "fro")]) # -0.1085989 (negligible)
cohen.d(data_2back$dis[data_2back$network == "dor"], data_2back$dis[!(data_2back$network == "dor")]) # -0.3534341 (small)
cohen.d(data_2back$dis[data_2back$network == "ven"], data_2back$dis[!(data_2back$network == "ven")]) # 0.01807474 (negligible)
cohen.d(data_2back$dis[data_2back$network == "def"], data_2back$dis[!(data_2back$network == "def")]) # 0.2594864 (small)
cohen.d(data_2back$dis[data_2back$network == "som"], data_2back$dis[!(data_2back$network == "som")]) # 0.527744 (medium)
```
# part 5 - performance and total problem score
- (1) data processing
```{r}
# for three way interactions
# preprocess subj_list
names(subj_list)[names(subj_list) == 'V1'] <- 'src_subject_id'
subj_list$src_subject_id <- gsub("NDAR_", "", subj_list$src_subject_id)
# merge dataframes
nback_sub <- nback[,c("src_subject_id","tfmri_nb_all_beh_c0bpf_rate","tfmri_nb_all_beh_c0bnf_rate","tfmri_nb_all_beh_c0bngf_rate","tfmri_nb_all_beh_c2bpf_rate","tfmri_nb_all_beh_c2bnf_rate","tfmri_nb_all_beh_c2bngf_rate")]
perf_psy <- merge(fam_site_id, nback_sub, by="src_subject_id")
perf_psy <- merge(perf_psy, abcd_cbcl, by="src_subject_id")
perf_psy <- merge(perf_psy, subj_list, by="src_subject_id")
perf_psy <- perf_psy %>%
mutate_at(c("interview_age","cbcl_scr_syn_totprob_t","cbcl_scr_syn_totprob_r","tfmri_nb_all_beh_c0bpf_rate","tfmri_nb_all_beh_c0bnf_rate","tfmri_nb_all_beh_c0bngf_rate","tfmri_nb_all_beh_c2bpf_rate","tfmri_nb_all_beh_c2bnf_rate","tfmri_nb_all_beh_c2bngf_rate"), as.numeric)
# head(perf_psy)
dim(perf_psy) # 5003 21
# select subjects with clean imaging data
perf_psy <- merge(data_hf_0b["src_subject_id"], perf_psy, by="src_subject_id")
# head(perf_psy)
dim(perf_psy) # 4952 21
# pivot longer
perf_psy_long <-
pivot_longer(perf_psy, cols = c("tfmri_nb_all_beh_c0bpf_rate","tfmri_nb_all_beh_c0bnf_rate","tfmri_nb_all_beh_c0bngf_rate","tfmri_nb_all_beh_c2bpf_rate","tfmri_nb_all_beh_c2bnf_rate","tfmri_nb_all_beh_c2bngf_rate"), names_to = "nback_condition", values_to = "nback_perf")
# extract 0-back and 2-back & emotions
perf_psy_long$nback_condition <- gsub("tfmri_nb_all_beh_","",as.character(perf_psy_long$nback_condition))
perf_psy_long$nback_condition <- gsub("_rate","",as.character(perf_psy_long$nback_condition))
perf_psy_long$nback_condition <- gsub("b","_b",as.character(perf_psy_long$nback_condition))
perf_psy_long <- perf_psy_long %>% separate(nback_condition,c("nback","emotion"))
perf_psy_long <- perf_psy_long %>%
mutate_at(c("nback","emotion"), as.factor)
# head(perf_psy_long)
dim(perf_psy_long) # 30018 18
```
- (2) full mix model
```{r}
lmer6_acc_psych_threelevel <- lmer(nback_perf ~ cbcl_scr_syn_totprob_r * nback * emotion + sex + interview_age + (1|site_id_l/rel_family_id), data = perf_psy_long)
summary(lmer6_acc_psych_threelevel)
tab_model(lmer6_acc_psych_threelevel, digits = 4)
```
- (3) three-way anova
```{r}
anova_model1 <- aov(nback_perf ~ cbcl_scr_syn_totprob_r * nback * emotion, data = perf_psy_long)
summary(anova_model1)
tab_model(anova_model1)
```
- (4) plot
```{r}
# Create the scatter plot and add the regression line
p <- ggplot(data = perf_psy_long, aes(x=cbcl_scr_syn_totprob_r, y=nback_perf, color=emotion, shape=emotion)) +
geom_point(aes(color=emotion),alpha = 0.5) +
facet_wrap(~nback) +
geom_smooth(method = "lm", se = TRUE, alpha = 0.1) +
scale_color_manual(values=c('#999999','#56B4E9','#E69F00')) +
labs(title = "", x = "CBCL total problems score", y = "Accuracy") +
theme_classic() +
theme(text = element_text(size = 20))
p
```
# part 6 - demographics table
- (1) preprocessing
```{r}
# demographic table includes
# (1) age
# (2) sex assigned at birth; female (%)
# (3) race/ethnicity (%)
# (4) family income (%)
# (5) caregiver education (%)
# select subjects with clean imaging data
abcd_pdem <- merge(perf_psy["src_subject_id"], abcd_pdem, by="src_subject_id")
abcd_pdem <- merge(abcd_pdem, fam_site_id, by="src_subject_id")
abcd_pdem <- abcd_pdem %>%
mutate_at("interview_age", as.numeric)
# sanity check
dim(abcd_pdem) # 4952 47
# head(abcd_pdem)
# re-code household income into three levels
household.income = abcd_pdem$demo_comb_income_v2
household.income[abcd_pdem$demo_comb_income_v2 == "1"] = 1 # "[<50K]"
household.income[abcd_pdem$demo_comb_income_v2 == "2"] = 1 # "[<50K]"
household.income[abcd_pdem$demo_comb_income_v2 == "3"] = 1 # "[<50K]"
household.income[abcd_pdem$demo_comb_income_v2 == "4"] = 1 # "[<50K]"
household.income[abcd_pdem$demo_comb_income_v2 == "5"] = 1 # "[<50K]"
household.income[abcd_pdem$demo_comb_income_v2 == "6"] = 1 # "[<50K]"
household.income[abcd_pdem$demo_comb_income_v2 == "7"] = 2 # "[>=50K & <100K]"
household.income[abcd_pdem$demo_comb_income_v2 == "8"] = 2 # "[>=50K & <100K]"
household.income[abcd_pdem$demo_comb_income_v2 == "9"] = 3 # "[>=100K]"
household.income[abcd_pdem$demo_comb_income_v2 == "10"] = 3 # "[>=100K]"
household.income[abcd_pdem$demo_comb_income_v2 == "777"] = NA
household.income[abcd_pdem$demo_comb_income_v2 == "999"] = NA
household.income[household.income %in% c(NA, "999", "777")] = NA
abcd_pdem$household.income = factor(household.income, levels= 1:3, labels = c("[<50K]", "[>=50K & <100K]", "[>=100K]") )
# re-code caregiver education into five levels
high.educ1 = abcd_pdem$demo_prnt_ed_v2
high.educ2 = abcd_pdem$demo_prtnr_ed_v2
high.educ1[which(high.educ1 == "999")] = NA
high.educ2[which(high.educ2 == "999")] = NA
high.educ1[which(high.educ1 == "777")] = NA
high.educ2[which(high.educ2 == "777")] = NA
high.educ = pmax(as.numeric(as.character(high.educ1)), as.numeric(as.character(high.educ2)), na.rm=T)
idx <- which(high.educ %in% 0:12, arr.ind = TRUE)
high.educ[idx] = 1 # "< HS Diploma"
idx <- which(high.educ %in% 13:14, arr.ind = TRUE)
high.educ[idx] = 2 # "HS Diploma/GED"
idx <- which(high.educ %in% 15:17, arr.ind = TRUE)
high.educ[idx] = 3 # "Some College"
idx <- which(high.educ == 18, arr.ind = TRUE)
high.educ[idx] = 4 # "Bachelor"
idx <- which(high.educ %in% 19:21, arr.ind = TRUE)
high.educ[idx] = 5 # "Post Graduate Degree"
high.educ[which(high.educ == "999")]=NA
high.educ[which(high.educ == "777")]=NA
abcd_pdem$high.educ = factor(high.educ, levels= 1:5, labels = c("< HS Diploma","HS Diploma/GED","Some College","Bachelor","Post Graduate Degree") )
# re-code racial categories with more easily identifiable variable names (into 8 levels)
abcd_pdem$white= (abcd_pdem$demo_race_a_p___10 == 1)*1
abcd_pdem$black= (abcd_pdem$demo_race_a_p___11 == 1)*1
abcd_pdem$asian = 0
abcd_pdem$asian[abcd_pdem$demo_race_a_p___18 == 1 | abcd_pdem$demo_race_a_p___19 == 1 |
abcd_pdem$demo_race_a_p___20 == 1 | abcd_pdem$demo_race_a_p___21 == 1 |
abcd_pdem$demo_race_a_p___22 == 1 | abcd_pdem$demo_race_a_p___23 == 1 |
abcd_pdem$demo_race_a_p___24==1] = 1
abcd_pdem$aian = 0
abcd_pdem$aian[abcd_pdem$demo_race_a_p___12 == 1 | abcd_pdem$demo_race_a_p___13 == 1] = 1
abcd_pdem$nhpi = 0
abcd_pdem$nhpi[abcd_pdem$demo_race_a_p___14 == 1 | abcd_pdem$demo_race_a_p___15 == 1 |
abcd_pdem$demo_race_a_p___16 == 1 | abcd_pdem$demo_race_a_p___17 == 1] = 1
abcd_pdem$other = 0
abcd_pdem$other[abcd_pdem$demo_race_a_p___25 == 1] = 1
# generate mixed category
abcd_pdem$mixed =
abcd_pdem$white + abcd_pdem$black + abcd_pdem$asian + abcd_pdem$aian + abcd_pdem$nhpi + abcd_pdem$other
abcd_pdem$mixed[abcd_pdem$mixed <= 1] = 0
abcd_pdem$mixed[abcd_pdem$mixed > 1] = 1
# generate single race.ethnicity variable
# note: If you want fewer categories, change code here so that different categories (e.g., other, mixed) have same value (e.g., both = 7)
abcd_pdem$race.eth = NA
abcd_pdem$race.eth[abcd_pdem$demo_ethn_v2 == 1] = 1
abcd_pdem$race.eth[abcd_pdem$white == 1] = 2
abcd_pdem$race.eth[abcd_pdem$black == 1] = 3
abcd_pdem$race.eth[abcd_pdem$asian == 1] = 4
abcd_pdem$race.eth[abcd_pdem$aian == 1] = 5
abcd_pdem$race.eth[abcd_pdem$nhpi == 1] = 6
abcd_pdem$race.eth[abcd_pdem$other == 1] = 7
abcd_pdem$race.eth[abcd_pdem$mixed == 1] = 8
abcd_pdem$race.eth <- factor(abcd_pdem$race.eth, levels = 1:8,
labels = c("Hispanic", "White", "Black", "Asian", "AIAN", "NHPI", "Other", "Mixed"))
# head(abcd_pdem)
```
- (2) generate the table
```{r}
# create a summary/demographics table
vars1 <- c("interview_age", "sex", "race.eth", "high.educ", "household.income")
tab1 <- CreateTableOne(vars = vars1, data = abcd_pdem)
tabAsStringMatrix <- print(tab1, printToggle = FALSE, noSpaces = TRUE)
tab1 = knitr::kable(tabAsStringMatrix)
tab1
```
# part 7 - mediation
- (1) mediation analysis - all three emotional comparisons
```{r}
library(mediation)
detach("package:lmerTest", unload = TRUE) # must unload lmerTest to run mediation
# detach("package:psych", unload = TRUE) # also be sure psych is detached
# using average accuracy of all faces
# IV - acc, DV - psy, med - sqED
### 0-back ###
# vis
lme4_0back_vis.0 <-
lmer(cbcl_scr_syn_totprob_r ~ allface_0b_rate + sex + interview_age + (1|site_id_l), data = data_0back[data_0back$network == "vis",])
lme4_0back_vis.M <-
lmer(dis ~ allface_0b_rate + sex + interview_age + (1|site_id_l), data = data_0back[data_0back$network == "vis",])
lme4_0back_vis.Y <-
lmer(cbcl_scr_syn_totprob_r ~ allface_0b_rate + dis + sex + interview_age + (1|site_id_l), data = data_0back[data_0back$network == "vis",])
lme4_med_0back_vis <-
mediate(model.m = lme4_0back_vis.M, model.y = lme4_0back_vis.Y, treat='allface_0b_rate', mediator = 'dis', data=data_0back[data_0back$network == "vis",], sims=1000, dropobs=TRUE)
summary(lme4_0back_vis.0)
summary(lme4_0back_vis.M)
summary(lme4_0back_vis.Y)
summary(lme4_med_0back_vis) # ACME -0.65252 -1.11414 -0.17 0.006 **
# lim
lme4_0back_lim.0 <-
lmer(cbcl_scr_syn_totprob_r ~ allface_0b_rate + sex + interview_age + (1|site_id_l), data = data_0back[data_0back$network == "lim",])
lme4_0back_lim.M <-
lmer(dis ~ allface_0b_rate + sex + interview_age + (1|site_id_l), data = data_0back[data_0back$network == "lim",])
lme4_0back_lim.Y <-
lmer(cbcl_scr_syn_totprob_r ~ allface_0b_rate + dis + sex + interview_age + (1|site_id_l), data = data_0back[data_0back$network == "lim",])
lme4_med_0back_lim <-
mediate(model.m = lme4_0back_lim.M, model.y = lme4_0back_lim.Y, treat='allface_0b_rate', mediator = 'dis', data=data_0back[data_0back$network == "lim",], sims=1000, dropobs=TRUE)
summary(lme4_0back_lim.0)
summary(lme4_0back_lim.M)
summary(lme4_0back_lim.Y)
summary(lme4_med_0back_lim) # ACME -0.08545 -0.26251 0.06 0.26
# v1
lme4_0back_v1.0 <-
lmer(cbcl_scr_syn_totprob_r ~ allface_0b_rate + sex + interview_age + (1|site_id_l), data = data_0back[data_0back$network == "v1",])
lme4_0back_v1.M <-
lmer(dis ~ allface_0b_rate + sex + interview_age + (1|site_id_l), data = data_0back[data_0back$network == "v1",])
lme4_0back_v1.Y <-
lmer(cbcl_scr_syn_totprob_r ~ allface_0b_rate + dis + sex + interview_age + (1|site_id_l), data = data_0back[data_0back$network == "v1",])
lme4_med_0back_v1 <-
mediate(model.m = lme4_0back_v1.M, model.y = lme4_0back_v1.Y, treat='allface_0b_rate', mediator = 'dis', data=data_0back[data_0back$network == "v1",], sims=1000, dropobs=TRUE)
summary(lme4_0back_v1.0)
summary(lme4_0back_v1.M)
summary(lme4_0back_v1.Y)
summary(lme4_med_0back_v1) # ACME -0.6463 -1.0603 -0.27 <2e-16 ***
## 2-back ###
# vis
lme4_2back_vis.0 <-
lmer(cbcl_scr_syn_totprob_r ~ allface_2b_rate + sex + interview_age + (1|site_id_l), data = data_2back[data_2back$network == "vis",])
lme4_2back_vis.M <-
lmer(dis ~ allface_2b_rate + sex + interview_age + (1|site_id_l), data = data_2back[data_2back$network == "vis",])
lme4_2back_vis.Y <-
lmer(cbcl_scr_syn_totprob_r ~ allface_2b_rate + dis + sex + interview_age + (1|site_id_l), data = data_2back[data_2back$network == "vis",])
lme4_med_2back_vis <-
mediate(model.m = lme4_2back_vis.M, model.y = lme4_2back_vis.Y, treat='allface_2b_rate', mediator = 'dis', data=data_2back[data_2back$network == "vis",], sims=1000, dropobs=TRUE)
summary(lme4_2back_vis.0)
summary(lme4_2back_vis.M)
summary(lme4_2back_vis.Y)
summary(lme4_med_2back_vis) # ACME -1.1415 -1.6326 -0.66 <2e-16 ***
# lim
lme4_2back_lim.0 <-
lmer(cbcl_scr_syn_totprob_r ~ allface_2b_rate + sex + interview_age + (1|site_id_l), data = data_2back[data_2back$network == "lim",])
lme4_2back_lim.M <-
lmer(dis ~ allface_2b_rate + sex + interview_age + (1|site_id_l), data = data_2back[data_2back$network == "lim",])
lme4_2back_lim.Y <-
lmer(cbcl_scr_syn_totprob_r ~ allface_2b_rate + dis + sex + interview_age + (1|site_id_l), data = data_2back[data_2back$network == "lim",])
lme4_med_2back_lim <-
mediate(model.m = lme4_2back_lim.M, model.y = lme4_2back_lim.Y, treat='allface_2b_rate', mediator = 'dis', data=data_2back[data_2back$network == "lim",], sims=1000, dropobs=TRUE)
summary(lme4_2back_lim.0)
summary(lme4_2back_lim.M)
summary(lme4_2back_lim.Y)
summary(lme4_med_2back_lim) # ACME -0.5729 -0.9724 -0.20 <2e-16 ***
# v1
lme4_2back_v1.0 <-
lmer(cbcl_scr_syn_totprob_r ~ allface_2b_rate + sex + interview_age + (1|site_id_l), data = data_2back[data_2back$network == "v1",])
lme4_2back_v1.M <-
lmer(dis ~ allface_2b_rate + sex + interview_age + (1|site_id_l), data = data_2back[data_2back$network == "v1",])
lme4_2back_v1.Y <-
lmer(cbcl_scr_syn_totprob_r ~ allface_2b_rate + dis + sex + interview_age + (1|site_id_l), data = data_2back[data_2back$network == "v1",])
lme4_med_2back_v1 <-
mediate(model.m = lme4_2back_v1.M, model.y = lme4_2back_v1.Y, treat='allface_2b_rate', mediator = 'dis', data=data_2back[data_2back$network == "v1",], sims=1000, dropobs=TRUE)
summary(lme4_2back_v1.0)
summary(lme4_2back_v1.M)
summary(lme4_2back_v1.Y)
summary(lme4_med_2back_v1) # ACME -0.47589 -0.77228 -0.19 0.002 **
```

ABCD_emo_psy_vis.Rmd at commit 1099ea8, no license · at the source

Overview

Authors: Yen-Chu Lin¹, Qingyang Meng¹, Angelica Lopez-Tucker¹, Janya Utkarsh¹, Fin Sterner¹, May I. Conley², Lena J. Skalaban³, Richard Watts⁴, Dylan G. Gee⁵, Arielle Baskin-Sommers⁵, B.J. Casey¹

Institutions: Columbia University (United States); Barnard College (United States); NYU Langone Health (United States); New York University (United States); University of Oregon (United States); University of Canterbury (New Zealand); Yale University (United States)

Journal: Biological psychiatry global open science, volume 6, issue 6, article 100799

Dates: received 25 March 2026; accepted 15 July 2026; published online 29 July 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1016/j.bpsgos.2026.100799 · PMID 42741195 · PMCID PMC13572311 · OpenAlex W7171673176

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: human (organism)

Methods: Connectivity, Statistics, fMRI & imaging, Single-unit activity, calcium imaging, Machine learning

Keywords: Adolescence, Cognition, Emotion, Neuroimaging, Psychopathology, Visual cortex

Journal subjects: Archival Report

Topic: Neural and Behavioral Psychology Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: NIH (U01DA041174)

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

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above.

yenchu-lin/abcd_vis_psychopathology

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 1099ea8d8a8558c240376c456e542d096a4cebb2, 23 June 2026

Languages: R (1)

Size: 3 files, 1 script

Software Heritage: not archived

Found in: the text, “Associations of Neural Dissimilarity With Cognit”

Holds: README, 1 notebook

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

Tools: lme4 (1 file), lmerTest (1 file), tidyverse (1 file)

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

27 September 2026: the link answers

2 files

ABCD_emo_psy_vis.Rmd, R, 1,900 lines
README.md, Text, 1 line

At the source: github.com/yenchu-lin/abcd_vis_psychopathology

Tracing map

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

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
1 script, each with its path and the digest of its content;
no match between paragraphs and code yet;
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.

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 11 authors, 6 keywords, 1 funder, 89 references.

Cite

This paper

Lin, Y.-C., Meng, Q., Lopez-Tucker, A., Utkarsh, J., Sterner, F., Conley, M. I., Skalaban, L. J., Watts, R., Gee, D. G., Baskin-Sommers, A., & Casey, B. (2026). Distinct Representations of Irrelevant Emotional Information in the Visual Network Are Associated With Psychopathology in Youth. Biological psychiatry global open science, 6(6), 100799. https://doi.org/10.1016/j.bpsgos.2026.100799

BibTeX

@article{lin2026distinct, author = {Lin, Yen-Chu and Meng, Qingyang and Lopez-Tucker, Angelica and Utkarsh, Janya and Sterner, Fin and Conley, May I. and Skalaban, Lena J. and Watts, Richard and Gee, Dylan G. and Baskin-Sommers, Arielle and Casey, B.J.}, title = {{Distinct Representations of Irrelevant Emotional Information in the Visual Network Are Associated With Psychopathology in Youth}}, journal = {Biological psychiatry global open science}, year = {2026}, month = jul, volume = {6}, number = {6}, pages = {100799}, publisher = {Elsevier}, issn = {2667-1743}, doi = {10.1016/j.bpsgos.2026.100799}, url = {https://doi.org/10.1016/j.bpsgos.2026.100799}, pmid = {42741195}, pmcid = {PMC13572311} }

RIS

TY - JOUR AU - Lin, Yen-Chu AU - Meng, Qingyang AU - Lopez-Tucker, Angelica AU - Utkarsh, Janya AU - Sterner, Fin AU - Conley, May I. AU - Skalaban, Lena J. AU - Watts, Richard AU - Gee, Dylan G. AU - Baskin-Sommers, Arielle AU - Casey, B.J. TI - Distinct Representations of Irrelevant Emotional Information in the Visual Network Are Associated With Psychopathology in Youth T2 - Biological psychiatry global open science J2 - Biol Psychiatry Glob Open Sci PY - 2026 DA - 2026/07/29 VL - 6 IS - 6 SP - 100799 SN - 2667-1743 PB - Elsevier DO - 10.1016/j.bpsgos.2026.100799 UR - https://doi.org/10.1016/j.bpsgos.2026.100799 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.bpsgos.2026.100799", "type": "article-journal", "title": "Distinct Representations of Irrelevant Emotional Information in the Visual Network Are Associated With Psychopathology in Youth", "container-title": "Biological psychiatry global open science", "author": [ { "family": "Lin", "given": "Yen-Chu" }, { "family": "Meng", "given": "Qingyang" }, { "family": "Lopez-Tucker", "given": "Angelica" }, { "family": "Utkarsh", "given": "Janya" }, { "family": "Sterner", "given": "Fin" }, { "family": "Conley", "given": "May I." }, { "family": "Skalaban", "given": "Lena J." }, { "family": "Watts", "given": "Richard" }, { "family": "Gee", "given": "Dylan G." }, { "family": "Baskin-Sommers", "given": "Arielle" }, { "family": "Casey", "given": "B.J." } ], "container-title-short": "Biol Psychiatry Glob Open Sci", "volume": "6", "issue": "6", "page": "100799", "DOI": "10.1016/j.bpsgos.2026.100799", "PMID": "42741195", "PMCID": "PMC13572311", "ISSN": "2667-1743", "publisher": "Elsevier", "URL": "https://doi.org/10.1016/j.bpsgos.2026.100799", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 29 ] ] } }

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.