## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## Warning: package 'broom' was built under R version 4.2.3
## Warning: package 'emmeans' was built under R version 4.2.3
## Warning: package 'sjPlot' was built under R version 4.2.3
## Warning: package 'reactable' was built under R version 4.2.3

IL6 analysis

library(readxl)
library(tidyr)

IL_6 <- read_excel("IL_6.xlsx") 


IL_6_long <- pivot_longer(IL_6, 
                       cols = ends_with(c("Plac", "IFN")), 
                       names_to = c("Time_point", ".value"), 
                       names_pattern = "(.*)_(.*)") %>% select(-starts_with("Delta")) %>%  pivot_longer(cols = c("Plac", "IFN"), names_to = "Condition", values_to="IL_6") %>%
  mutate(time_since_inj_h = case_when(
    Time_point == "Baseline" ~ "0",
    Time_point == "Time 1" ~ "4",
    Time_point == "Time 2" ~ "6.5",
    TRUE ~ "Unknown"
  )) %>% mutate(time_since_inj_h=as.numeric(time_since_inj_h))


library(lme4)
library(sjPlot)

#linear mixed model - all sub
m0 <- lmer(IL_6 ~ time_since_inj_h*Condition + (1 | ID), data = IL_6_long)
tab_model(m0)
  IL_6
Predictors Estimates CI p
(Intercept) 2.12 -0.06 – 4.29 0.056
time since inj h 1.84 1.37 – 2.31 <0.001
Condition [Plac] 0.09 -2.84 – 3.02 0.951
time since inj h ×
Condition [Plac] 		-1.00 -1.66 – -0.33 0.004
Random Effects
σ2 36.67
Ï„00 ID 3.29
ICC 0.08
N ID 30
Observations 180
Marginal R2 / Conditional R2 0.306 / 0.363 
plot_model(m0, type = "pred",terms = c("time_since_inj_h","Condition"),show.data = TRUE)

#linear mixed model - 27 subs
m0 <- lmer(IL_6 ~ time_since_inj_h*Condition + (1 | ID), data = IL_6_long %>% filter(!ID %in% c("008","022","023")))
tab_model(m0)
  IL_6
Predictors Estimates CI p
(Intercept) 1.95 -0.20 – 4.09 0.075
time since inj h 1.73 1.27 – 2.19 <0.001
Condition [Plac] 0.09 -2.80 – 2.97 0.953
time since inj h ×
Condition [Plac] 		-0.90 -1.55 – -0.24 0.008
Random Effects
σ2 31.92
Ï„00 ID 3.03
ICC 0.09
N ID 27
Observations 162
Marginal R2 / Conditional R2 0.309 / 0.369 
plot_model(m0, type = "pred",terms = c("time_since_inj_h","Condition"),show.data = TRUE)

#linear mixed model - 25 subs
m0 <- lmer(IL_6 ~ time_since_inj_h*Condition + (1 | ID), data = IL_6_long %>% filter(!ID %in% c("008","022","023","025","002")))
tab_model(m0)
  IL_6
Predictors Estimates CI p
(Intercept) 1.40 -0.82 – 3.62 0.213
time since inj h 1.80 1.33 – 2.27 <0.001
Condition [Plac] 0.68 -2.27 – 3.63 0.648
time since inj h ×
Condition [Plac] 		-0.94 -1.61 – -0.27 0.006
Random Effects
σ2 30.84
Ï„00 ID 3.65
ICC 0.11
N ID 25
Observations 150
Marginal R2 / Conditional R2 0.318 / 0.390 
plot_model(m0, type = "pred",terms = c("time_since_inj_h","Condition"),show.data = TRUE)

Functional connectivity

Mean correlation maps

library(abind)

# List to store all 2D matrices
matrices <- list()

# Loop through all .tsv files in the folder
files <- list.files(path = "conmats_filtered_tedana_bp_25/", pattern = "*.tsv",full.names = TRUE)
for (file in files) {
  mat <- read.table(file, sep = "\t", header = FALSE)
  matrices[[file]] <- mat
}

# Combine all matrices into a single 3D matrix
result_matrix <- abind(matrices, along = 3)
dim(result_matrix)
## [1] 410 410  50
results_brain_labs <- brainregions$X2

# List to store all 2D matrices
interferon_matrices <- list()
# Loop through all .tsv files in the folder
files <- list.files(path = "conmats_filtered_tedana_bp_25/", pattern = "ses-I",full.names = TRUE)
for (file in files) {
  mat <- read.table(file, sep = "\t", header = FALSE)
  interferon_matrices[[file]] <- mat
}
# Combine all matrices into a single 3D matrix
result_interferon_matrix <- abind(interferon_matrices, along = 3)

# List to store all 2D matrices
placebo_matrices <- list()
# Loop through all .tsv files in the folder
files <- list.files(path = "conmats_filtered_tedana_bp_25/", pattern = "ses-P",full.names = TRUE)
for (file in files) {
  mat <- read.table(file, sep = "\t", header = FALSE)
  placebo_matrices[[file]] <- mat
}
# Combine all matrices into a single 3D matrix
result_placebo_matrix <- abind(placebo_matrices, along = 3)

# nnodes <- length(brainregions$X2)
# tri_pos <- which(upper.tri(matrix(nrow = nnodes, ncol = nnodes)), arr.ind = T)
# Load the R.matlab package
library(R.matlab)

# Specify the file path and name
file <- "C:\\Users\\saptaf1\\Downloads\\archives\\NBS_benchmarking-master\\interferon_files\\tedananogsr_bp_funccon_54.mat"

# Write the 3D matrix to the .mat file
writeMat(funccon=atanh(result_matrix), file)

file <- "C:\\Users\\saptaf1\\Downloads\\archives\\NBS_benchmarking-master\\interferon_files\\yeo_410.mat"

writeMat(yeo410=adjacency_matrix, file)
n_missing <- sum(is.na(result_placebo_matrix))

mean_placebo_funcconn <- apply(result_placebo_matrix, c(1, 2), function(x) mean(x, na.rm = TRUE))
diag(mean_placebo_funcconn) <- 0
z_mean_placebo_funcconn <- atanh(mean_placebo_funcconn)
z_mean_placebo_funcconn[z_mean_placebo_funcconn == 0] <- .Machine$double.xmin
  

# problems with 1121023_ses-I,0821008_ses-I,1121022_ses-I
n_missing <- sum(is.na(result_interferon_matrix))
n_missing <- apply(result_interferon_matrix, c(3), function(x) sum(is.na(x)))
print(n_missing)
## conmats_filtered_tedana_bp_25/sub-0122026_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0122027_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0222028_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0222029_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0322030_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0322031_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0322032_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0422033_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0422034_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0422035_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0521003_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0522036_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0721005_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0821006_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0821007_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0921009_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0921010_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0921012_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0921014_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0921015_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0921016_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-0921017_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-1021019_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-1121020_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0 
## conmats_filtered_tedana_bp_25/sub-1121021_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                        0
mean_interferon_funcconn <- apply(result_interferon_matrix, c(1, 2), function(x) mean(x, na.rm = TRUE))
diag(mean_interferon_funcconn) <- 0
z_mean_interferon_funcconn <- atanh(mean_interferon_funcconn)
z_mean_interferon_funcconn[z_mean_interferon_funcconn == 0] <- .Machine$double.xmin


mean_diff_funcconn <- apply(result_placebo_matrix-result_interferon_matrix, c(1, 2), mean)
diag(mean_diff_funcconn) <- 0
#z_mean_diff_funcconn <- atanh(mean_diff_funcconn)
z_mean_diff_funcconn <- z_mean_placebo_funcconn-z_mean_interferon_funcconn

h1 <- Heatmap(z_mean_placebo_funcconn, name = "Fisher Z", row_split = network.combined$name, column_split = network.combined$name, row_title_rot = 0, show_column_dend = FALSE, show_row_dend = FALSE, show_column_names = FALSE, column_title = "Placebo", column_title_gp = gpar(fontsize = 20, fontface = "bold"), clustering_distance_rows="pearson", clustering_distance_columns ="pearson")

h2 <- Heatmap(z_mean_interferon_funcconn, name = "Fisher Z", row_split= factor(network.combined$name, levels = (names(row_order(h1)))), column_split= factor(network.combined$name, levels = (names(row_order(h1)))), row_title_rot = 0, show_column_dend = FALSE, show_row_dend = FALSE, show_column_names = FALSE, column_title = "Interferon", column_title_gp = gpar(fontsize = 20, fontface = "bold"), cluster_rows = FALSE, cluster_columns = FALSE, cluster_row_slices = FALSE, cluster_column_slices = FALSE, row_order = unlist(row_order(h1)), column_order = unlist(row_order(h1)) )

h3_v2 <- Heatmap(z_mean_diff_funcconn, name = "Fisher Z", row_split= factor(network.combined$name, levels = (names(row_order(h1)))), column_split= factor(network.combined$name, levels = (names(row_order(h1)))),  row_title_rot = 0, show_column_dend = FALSE, show_row_dend = FALSE, show_column_names = FALSE, column_title = "Placebo-Interferon", column_title_gp = gpar(fontsize = 20, fontface = "bold"), cluster_rows = FALSE, cluster_columns = FALSE, cluster_row_slices = FALSE, cluster_column_slices = FALSE, row_order = unlist(row_order(h1)), column_order = unlist(row_order(h1)))

h3_v2_adjusted <- Heatmap(z_mean_diff_funcconn, name = "Fisher Z", row_split= factor(network.combined$name, levels = (names(row_order(h1)))), column_split= factor(network.combined$name, levels = (names(row_order(h1)))),  row_title_rot = 0, show_column_dend = FALSE, show_row_dend = FALSE, show_column_names = FALSE, column_title = "Placebo-Interferon", column_title_gp = gpar(fontsize = 20, fontface = "bold"), cluster_rows = FALSE, cluster_columns = FALSE, cluster_row_slices = FALSE, cluster_column_slices = FALSE, row_order = unlist(row_order(h1)), column_order = unlist(row_order(h1)), col=colorRamp2(c(-1, 0, 1), c("blue", "white", "red")))

h1+h2+h3_v2_adjusted

h1+h2

h3_v2

hist(z_mean_diff_funcconn[upper.tri(z_mean_diff_funcconn, diag = FALSE)], main = "Histogram of functional connectivity difference values")

g1  <- graph.adjacency(z_mean_placebo_funcconn,weighted=TRUE,mode = "upper", diag = FALSE)
g2  <- graph.adjacency(z_mean_interferon_funcconn,weighted=TRUE,mode = "upper", diag = FALSE)

df <- bind_rows(get.data.frame(g1) %>% mutate(condition="Placebo"), get.data.frame(g2) %>% mutate(condition="Interferon"))

ggplot(df, aes(x=weight, fill=condition)) +
  geom_density(alpha=.25)

df %>% mutate(from = gsub("V", "ROI_", from), to = gsub("V", "ROI_", to)) %>% rowwise() %>% mutate(from_network = ifelse(is.na(match(from, network.combined$ROI)), from, network.combined$name[match(from, network.combined$ROI)]), to_network = ifelse(is.na(match(to, network.combined$ROI)), to, network.combined$name[match(to, network.combined$ROI)]), from2to = paste0(sort(c(from_network,to_network))[1],"2",sort(c(from_network,to_network))[2])) -> df_ext

ggplot(df_ext, aes(x=weight, fill=condition)) +
  geom_density(alpha=.25) + facet_wrap(~from2to)

df_ext %>% mutate(fromNum=parse_number(from), toNum = parse_number(to)) -> df_ext

df_ext$from2toIndex <- as.integer(factor(df_ext$from2to))

nodes <- unique(c(df_ext$fromNum, df_ext$toNum))
adjacency_matrix <- matrix(0, nrow = length(nodes), ncol = length(nodes))
colnames(adjacency_matrix) <- rownames(adjacency_matrix) <- nodes
adjacency_matrix[cbind(match(df_ext$fromNum, nodes), match(df_ext$toNum, nodes))] <- df_ext$from2toIndex

write.csv(adjacency_matrix, file = "yeo_edge_group.csv", row.names = FALSE)

Wilcoxon

tmp = data.frame()
# run np test within-subject differences in edges
for (i in seq(1,dim(result_placebo_matrix)[3])) {
  # print(i)
  A <- atanh(result_placebo_matrix[,,i])
  A[A == 0] <- .Machine$double.xmin
  
  B <- atanh(result_interferon_matrix[,,i])
  B[B == 0] <- .Machine$double.xmin
  
  g1  <- graph.adjacency(A, weighted=TRUE, mode = "upper", diag = FALSE)
  g2  <- graph.adjacency(B, weighted=TRUE, mode = "upper", diag = FALSE)
  
  if (nrow(get.data.frame(g1)) != nrow(get.data.frame(g2))) {
    break
  }
  
  df <- bind_rows(get.data.frame(g1) %>% mutate(condition="Placebo"), get.data.frame(g2) %>% mutate(condition="Interferon")) %>% mutate(sub = i)
  
  tmp = bind_rows(tmp,df)
}
# Perform Wilcoxon signed-rank test
#DT::datatable(tmp %>% wilcox_test(weight ~ condition, paired = TRUE, detailed = TRUE, exact=TRUE))

# Perform Wilcoxon signed-rank test - unthresholded
print(wilcox.test(tmp %>% filter(condition=="Placebo")%>%.$weight, tmp %>% filter(condition=="Interferon")%>%.$weight, paired = TRUE))
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  tmp %>% filter(condition == "Placebo") %>% .$weight and tmp %>% filter(condition == "Interferon") %>% .$weight
## V = 1.3057e+12, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
EDGES <- 83845
dataframe1 <- tmp %>% filter(condition=="Placebo")%>% group_by(sub) %>% top_n(round(EDGES*0.1),weight)
dataframe2 <- tmp %>% filter(condition=="Interferon")

columns_to_subset <- c("from", "to", "sub") # replace with your desired subset of columns

filtered_dataframe2 <- dataframe2 %>%
  left_join(dataframe1, by = columns_to_subset) %>% na.omit() # filter rows to only include exact matches in dataframe1 and retain all columns from dataframe2

print(wilcox.test(dataframe1%>%.$weight, filtered_dataframe2%>%.$weight.x, paired = TRUE))
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  dataframe1 %>% .$weight and filtered_dataframe2 %>% .$weight.x
## V = 1.9185e+10, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
ggplot(filtered_dataframe2 %>% 
  pivot_longer(cols = starts_with("weight"), 
               names_to = c(".value", "group"), 
               names_pattern = "(weight).(.)") %>% mutate(group=ifelse(group=="x","Interferon","Placebo")), aes(x=weight, fill=group)) +
  geom_density(alpha=.25)

EDGES <- 83845
dataframe1 <- tmp %>% filter(condition=="Placebo")%>% group_by(sub) %>% top_n(round(EDGES*0.35),weight)
dataframe2 <- tmp %>% filter(condition=="Interferon")

columns_to_subset <- c("from", "to", "sub") # replace with your desired subset of columns

filtered_dataframe2 <- dataframe2 %>%
  left_join(dataframe1, by = columns_to_subset) %>% na.omit() # filter rows to only include exact matches in dataframe1 and retain all columns from dataframe2


ggplot(filtered_dataframe2 %>% 
  pivot_longer(cols = starts_with("weight"), 
               names_to = c(".value", "group"), 
               names_pattern = "(weight).(.)") %>% mutate(group=ifelse(group=="x","Interferon","Placebo")), aes(x=weight, fill=group)) +
  geom_density(alpha=.25)

print(wilcox.test(dataframe1%>%.$weight, filtered_dataframe2%>%.$weight.x, paired = TRUE))
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  dataframe1 %>% .$weight and filtered_dataframe2 %>% .$weight.x
## V = 2.1934e+11, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
rm(tmp)

Density plots

tmp = data.frame()
# run np test within-subject differences in edges
for (i in seq(1,dim(result_placebo_matrix)[3])) {
  # print(i)
  A <- atanh(result_placebo_matrix[,,i])
  A[A == 0] <- .Machine$double.xmin
  
  B <- atanh(result_interferon_matrix[,,i])
  B[B == 0] <- .Machine$double.xmin
  
  g1  <- graph.adjacency(A, weighted=TRUE, mode = "upper", diag = FALSE)
  g2  <- graph.adjacency(B, weighted=TRUE, mode = "upper", diag = FALSE)
  
  if (nrow(get.data.frame(g1)) != nrow(get.data.frame(g2))) {
    break
  }
  
  df <- bind_rows(get.data.frame(g1) %>% mutate(condition="Placebo"), get.data.frame(g2) %>% mutate(condition="Interferon")) %>% mutate(sub = i)
  
  tmp = bind_rows(tmp,df)
}

tmp.propthr <- tmp %>% group_by(sub) %>% top_n(round(EDGES*0.1),weight)

ggplot(tmp.propthr, aes(x=weight, fill=condition)) +
  geom_density(alpha=.25)

tmp.propthr <- tmp %>% group_by(sub) %>% top_n(round(EDGES*0.35),weight)

ggplot(tmp.propthr, aes(x=weight, fill=condition)) +
  geom_density(alpha=.25)

### Functional connectivity: correlation with measures
files <- list.files(path = "conmats_filtered_tedana_bp_25/", pattern = "ses-P",full.names = TRUE)
sub_numbers <- gsub(".*sub-(\\d+)_ses.*", "\\1", files)
sub_numbers_df <- data.frame(sub = 1:length(sub_numbers), Subj_ID = as.numeric(sub_numbers))

dataframe1 <- tmp %>% filter(condition=="Placebo")%>% group_by(sub) %>% top_n(round(EDGES*0.1),weight)
dataframe2 <- tmp %>% filter(condition=="Interferon")
columns_to_subset <- c("from", "to", "sub") # replace with your desired subset of columns
filtered_dataframe2 <- dataframe2 %>%
  left_join(dataframe1, by = columns_to_subset) %>% na.omit() 

connectivity_diff <- filtered_dataframe2 %>% group_by(sub,from,to) %>% mutate(weight_diff =weight.y-weight.x)
connectivity_diff<- connectivity_diff %>% left_join(sub_numbers_df)

variables_ext %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") %>% group_by(Subj_ID, measure) %>% summarise(value_diff = diff(value)) -> var_diff

# for (n in unique(var_diff$measure)) {
#   curr_diff <- var_diff %>% filter(measure==n)
#   connectivity_diff <- connectivity_diff %>% left_join(curr_diff) %>% mutate(from2to=paste(from,to,sep = "_"))
#   m1<-lmer(weight_diff~value_diff+from2to+(1|Subj_ID),data = connectivity_diff)
# }

for (n in unique(var_diff$measure)) {
  print(n)
  start_time <- Sys.time()
  curr_diff <- var_diff %>% filter(measure==n)
  corr_df <- connectivity_diff %>% left_join(curr_diff, by=c("Subj_ID"="Subj_ID"), multiple = "all")  %>% mutate(from2to=factor(paste(from,to,sep = "_")), measure=factor(measure)) %>% group_by(measure,from2to) %>% rstatix::cor_test(weight_diff, value_diff,use = "pairwise.complete.obs")
  end_time <- Sys.time()
  print(start_time-end_time)
}

# correlation
corr_df <- connectivity_diff %>% left_join(var_diff, by=c("Subj_ID"="Subj_ID"),multiple = "all")  %>% mutate(from2to=paste(from,to,sep = "_")) %>% group_by(measure,from2to) %>% rstatix::cor_test(weight_diff, value_diff,use = "pairwise.complete.obs")

List of most Interferon affected connections in graph

Where is the most condition related change happening?

library(igraph)

g_diff_pMinusI  <- graph.adjacency(z_mean_diff_funcconn,weighted=TRUE,mode = "upper", diag = FALSE)

df_diff_pMinusI <- get.data.frame(g_diff_pMinusI) %>% mutate(from = gsub("V", "ROI_", from), to = gsub("V", "ROI_", to)) %>% rowwise() %>% mutate(from_network = ifelse(is.na(match(from, network.combined$ROI)), from, network.combined$name[match(from, network.combined$ROI)]), to_network = ifelse(is.na(match(to, network.combined$ROI)), to, network.combined$name[match(to, network.combined$ROI)]), from2to = paste0(sort(c(from_network,to_network))[1],"2",sort(c(from_network,to_network))[2])) %>% mutate(from_roi = ifelse(is.na(match(from, brainregions$ROI_X)), from, brainregions$X2[match(from, brainregions$ROI_X)]), to_roi = ifelse(is.na(match(to, brainregions$ROI_X)), from, brainregions$X2[match(to, brainregions$ROI_X)]), abs_weight = abs(weight))

df_diff_pMinusI$abs_weight_decile <- ntile(df_diff_pMinusI$abs_weight, 10)  


# top 20% of affected edges
df_diff_pMinusI %>% filter(abs_weight_decile >= 9) %>% ggplot(., aes(abs_weight, fill = from2to)) +
  geom_histogram(binwidth = 0.07)

df_diff_pMinusI %>% filter(abs_weight_decile >= 9)  %>% ggplot(., aes(from2to)) +
  geom_bar() + coord_flip()

# top 10% of affected edges
df_diff_pMinusI %>% filter(abs_weight_decile >= 10) %>% ggplot(., aes(abs_weight, fill = from2to)) +
  geom_histogram(binwidth = 0.1)

df_diff_pMinusI %>% filter(abs_weight_decile >= 10)  %>% ggplot(., aes(from2to)) +
  geom_bar() + coord_flip()

DT::datatable(df_diff_pMinusI  %>% filter(abs_weight_decile >= 10))
library(igraph)

df_list <- list()

for (i in 1:dim(result_matrix)[3]) {
  #print(i)
  g  <- graph.adjacency(result_matrix[,,i],weighted=TRUE,mode = "upper", diag = FALSE)
  df_list[[i]] <- get.data.frame(g) %>% mutate(fileIndex=i)
}

data.funccon <- bind_rows(df_list) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002))

Grouped by scale, condition, subject

Unthresholded

data.funccon %>% mutate(weight = atanh(weight)) %>% group_by(Subj_ID,condition) %>% summarise(meanFuncCon = mean(weight)) -> data.aggr


anova_summary <- nice(aov_ez("Subj_ID", "meanFuncCon", data.aggr, within ="condition"))

ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanFuncCon,
  type             = "p",
  bf.message = FALSE)

DT::datatable(anova_summary)

Only positive weight > 0

data.funccon %>% mutate(weight = atanh(weight)) %>% group_by(Subj_ID,condition) %>% filter(weight>0) %>% summarise(meanFuncCon = mean(weight))-> data.aggr


anova_summary <- nice(aov_ez("Subj_ID", "meanFuncCon", data.aggr, within ="condition"))

ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanFuncCon,
  type             = "p",
  bf.message = FALSE)

DT::datatable(anova_summary)

Proportional thresholding based on on weight (not abs_weight)

# total edges per session: EDGES
plots_list <- list()
x <- c(0.1,0.15,0.2,0.25,0.3,0.35)
for (i in 1:6) {
  
  p <- ggwithinstats(
      data             = data.funccon %>% mutate(weight = atanh(weight), abs_weight = abs(weight)) %>% group_by(Subj_ID,condition) %>% top_n(round(EDGES*x[i]),weight) %>% dplyr::summarise(meanFuncCon = mean(weight)),
    x                = condition,
    y                = meanFuncCon,
    type             = "p",
    title =x[i] ,
    bf.message = FALSE)
  
    plots_list[[i]] <- p

}

# arrange the plots in a grid using grid.arrange from gridExtra
gridExtra::grid.arrange(grobs = plots_list)

Grouped by scale, condition, subject and age_cat

data.funccon %>% mutate(weight = atanh(weight)) %>% group_by(Subj_ID,condition,Age_cat) %>% summarise(meanFuncCon = mean(weight)) -> data.aggr

anova_summary <- nice(aov_ez("Subj_ID", "meanFuncCon", data.aggr, within ="condition", between="Age_cat"))
DT::datatable(anova_summary)

Grouped by scale, condition, subject and network:yeo

  • Unthresholded (long time to run)
data.funccon %>% mutate(weight = atanh(weight), from = gsub("V", "ROI_", from), to = gsub("V", "ROI_", to)) %>% rowwise() %>% mutate(from_network = ifelse(is.na(match(from, network.combined$ROI)), from, network.combined$name[match(from, network.combined$ROI)]), to_network = ifelse(is.na(match(to, network.combined$ROI)), to, network.combined$name[match(to, network.combined$ROI)]), from2to = paste0(sort(c(from_network,to_network))[1],"2",sort(c(from_network,to_network))[2])) -> data.funccon.ext


data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition,from2to)  -> tmp

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to, tmp$condition)
df_counts <- table(tmp$from2to)

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to)
# Convert the counts to a data frame and add column names
df_counts <- as.data.frame(df_counts)
names(df_counts) <- c("network", "count")

# Create pie chart faceted by condition
ggplot(df_counts, aes(x = "", y = count, fill = network)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  theme_void()

data.funccon.ext %>% group_by(Subj_ID,condition,from2to) %>% summarise(meanFuncConn = mean(weight)) -> data.funccon.aggr

# grouped_ggwithinstats(
#     data             = data.aggr,
#   x                = condition,
#   y                = meanFuncConn,
#   grouping.var = from2to,
#   type             = "p",
#   bf.message = FALSE)

# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"))

emmip(main_anova_summary, from2to ~ condition)

EMM <- emmeans(main_anova_summary, ~  condition | from2to)   # where treat has 2 levels
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# aov_test <- aov(meanFuncConn ~ condition*from2to + Error(Subj_ID/(condition*from2to)), data=data.funccon.aggr)
# report(aov_test)

# uncorrected p values
# data.funccon.aggr %>% ungroup() %>%
#   nest_by(from2to) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanFuncConn", data, within ="condition")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
DT::datatable(df_counts)
DT::datatable(nice(main_anova_summary))
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))

10% threshold

data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition) %>% top_n(round(EDGES*0.1),weight) %>% group_by(Subj_ID,condition, from2to) %>% dplyr::summarise(meanFuncConn = mean(weight)) -> data.funccon.aggr

data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition,from2to) %>% top_n(round(EDGES*0.1),weight) -> tmp

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to, tmp$condition)
df_counts <- table(tmp$from2to)

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to)
# Convert the counts to a data frame and add column names
df_counts <- as.data.frame(df_counts)
names(df_counts) <- c("network", "count")

# Create pie chart faceted by condition
ggplot(df_counts, aes(x = "", y = count, fill = network)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  theme_void()

# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"))

main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"))

emmip(main_anova_summary, from2to ~ condition)

EMM <- emmeans(main_anova_summary, ~  condition | from2to)   # where treat has 2 levels
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(df_counts)
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"))
DT::datatable(nice(main_anova_summary))
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"))
DT::datatable(nice(main_anova_summary))
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
data.funccon.aggr <- data.funccon.aggr %>% left_join(variables)

data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition,from2to) %>% top_n(round(EDGES*0.1),weight) -> tmp

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to, tmp$condition)
df_counts <- table(tmp$from2to)

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to)
# Convert the counts to a data frame and add column names
df_counts <- as.data.frame(df_counts)
names(df_counts) <- c("network", "count")

# Create pie chart faceted by condition
ggplot(df_counts, aes(x = "", y = count, fill = network)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  theme_void()

# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"), between="Age_cat")

main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"), between="Age_cat")

emmip(main_anova_summary, from2to ~ condition | Age_cat)

EMM <- emmeans(main_anova_summary, ~condition:Age_cat | from2to ) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(df_counts)
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"), between="Age_cat")
DT::datatable(nice(main_anova_summary))
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"), between="Age_cat")
DT::datatable(nice(main_anova_summary))
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))

25% threshold

data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition) %>% top_n(round(EDGES*0.25),weight) %>% group_by(Subj_ID,condition, from2to) %>% dplyr::summarise(meanFuncConn = mean(weight)) -> data.funccon.aggr

data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition,from2to) %>% top_n(round(EDGES*0.25),weight) -> tmp

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to, tmp$condition)
df_counts <- table(tmp$from2to)

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to)
# Convert the counts to a data frame and add column names
df_counts <- as.data.frame(df_counts)
names(df_counts) <- c("network", "count")

# Create pie chart faceted by condition
ggplot(df_counts, aes(x = "", y = count, fill = network)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  theme_void()

# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"))

main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"))

emmip(main_anova_summary, from2to ~ condition)

EMM <- emmeans(main_anova_summary, ~  condition | from2to)   # where treat has 2 levels
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(df_counts)
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"))
DT::datatable(nice(main_anova_summary))
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"))
DT::datatable(nice(main_anova_summary))
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
data.funccon.aggr <- data.funccon.aggr %>% left_join(variables)

data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition,from2to) %>% top_n(round(EDGES*0.25),weight) -> tmp

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to, tmp$condition)
df_counts <- table(tmp$from2to)

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to)
# Convert the counts to a data frame and add column names
df_counts <- as.data.frame(df_counts)
names(df_counts) <- c("network", "count")

# Create pie chart faceted by condition
ggplot(df_counts, aes(x = "", y = count, fill = network)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  theme_void()

# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"), between="Age_cat")

main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"), between="Age_cat")

emmip(main_anova_summary, from2to ~ condition | Age_cat)

EMM <- emmeans(main_anova_summary, ~condition:Age_cat | from2to ) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(df_counts)
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"), between="Age_cat")
DT::datatable(nice(main_anova_summary))
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"), between="Age_cat")
DT::datatable(nice(main_anova_summary))
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))

35% threshold

data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition) %>% top_n(round(EDGES*0.35),weight) %>% group_by(Subj_ID,condition, from2to) %>% dplyr::summarise(meanFuncConn = mean(weight)) -> data.funccon.aggr

data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition,from2to) %>% top_n(round(EDGES*0.35),weight) -> tmp

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to, tmp$condition)
df_counts <- table(tmp$from2to)

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to)
# Convert the counts to a data frame and add column names
df_counts <- as.data.frame(df_counts)
names(df_counts) <- c("network", "count")

# Create pie chart faceted by condition
ggplot(df_counts, aes(x = "", y = count, fill = network)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  theme_void()

# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"))

# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"))

emmip(main_anova_summary, from2to ~ condition)

EMM <- emmeans(main_anova_summary, ~  condition | from2to)   # where treat has 2 levels
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(df_counts)
# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"))
DT::datatable(nice(main_anova_summary))
# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"))
DT::datatable(nice(main_anova_summary))
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
data.funccon.aggr <- data.funccon.aggr %>% left_join(variables)

data.funccon.ext %>% mutate(abs_weight = abs(weight)) %>% group_by(Subj_ID,condition,from2to) %>% top_n(round(EDGES*0.1),weight) -> tmp

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to, tmp$condition)
df_counts <- table(tmp$from2to)

# Calculate the number of occurrences of each network within each condition
df_counts <- table(tmp$from2to)
# Convert the counts to a data frame and add column names
df_counts <- as.data.frame(df_counts)
names(df_counts) <- c("network", "count")

# Create pie chart faceted by condition
ggplot(df_counts, aes(x = "", y = count, fill = network)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  theme_void()

# main anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"), between="Age_cat")

main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"), between="Age_cat")

emmip(main_anova_summary, from2to ~ condition | Age_cat)

EMM <- emmeans(main_anova_summary, ~condition:Age_cat | from2to ) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(df_counts)
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition"), between="Age_cat")
DT::datatable(nice(main_anova_summary))
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.aggr,within =c("condition","from2to"), between="Age_cat")
DT::datatable(nice(main_anova_summary))
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))

Grouped by scale, condition, subject and network:yeo and age

#slow to calc
data.funccon %>% mutate(weight = atanh(weight), from = gsub("V", "ROI_", from), to = gsub("V", "ROI_", to)) %>% rowwise() %>% mutate(from_network = ifelse(is.na(match(from, network.combined$ROI)), from, network.combined$name[match(from, network.combined$ROI)]), to_network = ifelse(is.na(match(to, network.combined$ROI)), to, network.combined$name[match(to, network.combined$ROI)]), from2to = paste0(sort(c(from_network,to_network))[1],"2",sort(c(from_network,to_network))[2])) %>% group_by(Subj_ID,condition,from2to, Age_cat) %>% summarise(meanFuncConn = mean(weight)) -> data.funccon.network.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanFuncConn", data.funccon.network.aggr,within =c("condition","from2to"), between = "Age_cat")

EMM <- emmeans(main_anova_summary, ~  condition)   
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)



# data.funccon.network.aggr %>% ungroup() %>%
#   nest_by(from2to) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanFuncConn", data, within ="condition", between = "Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
DT::datatable(nice(main_anova_summary))
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))

Summary

NBS

Node strength - using BCT Toolbox’s strengths_und.m

Thresholded

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanStrength = mean(strength)) -> data.aggr

anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr, within =c("condition","threshold"))


emmip(anova_summary, condition~threshold )

grouped_ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanStrength,
  grouping.var     = threshold,
  type             = "p",
  bf.message = FALSE,
  exact = FALSE
  )

DT::datatable(nice(anova_summary))

Grouped by scale, condition, subject and age_cat

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition, Age_cat) %>% summarise(meanStrength = mean(strength)) -> data.aggr

anova_summary <-aov_ez("Subj_ID", "meanStrength", data.aggr, within =c("condition","threshold"), between = "Age_cat")


# 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
# 
# # aw <- aov_ez("Subj_ID", "meanDegree", data.aggr, within ="condition", between="Age_cat")
# # p_an <- afex_plot(aw, x = "condition", trace = "Age_cat")  + ggtitle("test") +   theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm"))
# 
# # plot anovas
# 
# d_nested <- data.aggr %>% ungroup() %>%
#   nest_by(threshold)
# 
# d_plots <- 
#   d_nested %>% 
#   mutate(aov = list(map(data,~aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>%
#   mutate(plot = list(map2(aov,threshold,~afex_plot(., x = "condition", trace = "Age_cat") +ggtitle(paste("threshold:",threshold)) + theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm")))))
# 
# d_plots$plot_one <- lapply(d_plots$plot, "[[", 1)
# 
# library(gridExtra)
# do.call(grid.arrange, c(d_plots$plot_one))
DT::datatable(nice(anova_summary))

Grouped by scale, condition, subject and network:yeo

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanStrength = mean(strength)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
#emmip(main_anova_summary, condition~threshold)

# average across threshold
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name"))

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=1
DT::datatable(nice(aov_ez("Subj_ID", "meanStrength", data.aggr %>% filter(threshold=="0.1"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.1"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=25
DT::datatable(nice(aov_ez("Subj_ID", "meanStrength", data.aggr %>% filter(threshold=="0.25"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.25"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=3.5
DT::datatable(nice(aov_ez("Subj_ID", "meanStrength", data.aggr %>% filter(threshold=="0.35"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.35"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Grouped by scale, condition, subject, age_cat and network:yeo

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name,Age_cat) %>% summarise(meanStrength = mean(strength)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name","threshold"), between = "Age_cat")

DT::datatable(nice(main_anova_summary)) 
emmip(main_anova_summary, name~condition | Age_cat)

emmip(main_anova_summary, name~Age_cat)

main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name"),between = "Age_cat")

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name","threshold"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
# focus on threshold=1
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr %>% filter(threshold==0.1),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=1
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr %>% filter(threshold==0.25),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=3.5
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr %>% filter(threshold==0.35),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition",between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Hubness?

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% mutate(quartile = ntile(strength, 4)) -> data.aggr

fit1 <- lmer(strength ~ condition*quartile + (1 | Subj_ID), data.aggr )

fit2 <- lmer(strength ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr )

fit3 <- lmer(strength ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr )

tab_model(fit1,fit2,fit3)
  strength strength strength
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -20.07 -23.95 – -16.19 <0.001 -18.98 -24.18 – -13.79 <0.001 -24.05 -29.30 – -18.80 <0.001
condition [Placebo] -2.40 -3.74 – -1.06 <0.001 -2.40 -3.74 – -1.06 <0.001 -0.39 -2.17 – 1.39 0.667
quartile 35.32 34.97 – 35.66 <0.001 35.32 34.97 – 35.66 <0.001 37.81 37.35 – 38.27 <0.001
condition [Placebo] ×
quartile 		3.86 3.37 – 4.34 <0.001 3.86 3.37 – 4.34 <0.001 2.12 1.47 – 2.77 <0.001
Age cat [Y] -2.47 -10.18 – 5.24 0.530 9.05 1.13 – 16.97 0.025
condition [Placebo] × Age
cat [Y] 		-4.57 -7.25 – -1.88 0.001
quartile × Age cat [Y] -5.67 -6.36 – -4.98 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		3.94 2.95 – 4.92 <0.001
Random Effects
σ2 1194.77 1194.77 1187.71
Ï„00 92.35 Subj_ID 94.76 Subj_ID 94.76 Subj_ID
ICC 0.07 0.07 0.07
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 61500 61500 61500
Marginal R2 / Conditional R2 0.576 / 0.607 0.576 / 0.607 0.579 / 0.610 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

fit1 <- lmer(strength ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit2 <- lmer(strength ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit3 <- lmer(strength ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

tab_model(fit1,fit2,fit3)
  strength strength strength
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -24.97 -26.73 – -23.20 <0.001 -24.50 -26.81 – -22.19 <0.001 -28.39 -30.78 – -26.00 <0.001
condition [Placebo] -3.04 -4.07 – -2.01 <0.001 -3.04 -4.07 – -2.01 <0.001 -2.13 -3.49 – -0.76 0.002
quartile 24.64 24.37 – 24.90 <0.001 24.64 24.37 – 24.90 <0.001 26.37 26.01 – 26.72 <0.001
condition [Placebo] ×
quartile 		2.62 2.24 – 2.99 <0.001 2.62 2.24 – 2.99 <0.001 1.92 1.42 – 2.42 <0.001
Age cat [Y] -1.06 -4.38 – 2.26 0.532 7.79 4.19 – 11.39 <0.001
condition [Placebo] × Age
cat [Y] 		-2.06 -4.12 – -0.00 0.050
quartile × Age cat [Y] -3.93 -4.46 – -3.40 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		1.59 0.84 – 2.35 <0.001
Random Effects
σ2 236.57 236.57 233.16
Ï„00 16.91 Subj_ID 17.36 Subj_ID 17.37 Subj_ID
ICC 0.07 0.07 0.07
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.770 / 0.785 0.769 / 0.785 0.772 / 0.788 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

fit1 <- lmer(strength ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

fit2 <- lmer(strength ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

fit3 <- lmer(strength ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

tab_model(fit1,fit2,fit3)
  strength strength strength
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -22.68 -26.86 – -18.50 <0.001 -21.47 -27.05 – -15.89 <0.001 -27.20 -32.85 – -21.55 <0.001
condition [Placebo] -2.67 -4.14 – -1.19 <0.001 -2.67 -4.14 – -1.19 <0.001 -0.28 -2.24 – 1.67 0.775
quartile 38.72 38.34 – 39.10 <0.001 38.72 38.34 – 39.10 <0.001 41.51 41.00 – 42.02 <0.001
condition [Placebo] ×
quartile 		4.18 3.64 – 4.72 <0.001 4.18 3.64 – 4.72 <0.001 2.24 1.52 – 2.95 <0.001
Age cat [Y] -2.75 -11.03 – 5.53 0.515 10.28 1.76 – 18.79 0.018
condition [Placebo] × Age
cat [Y] 		-5.42 -8.37 – -2.47 <0.001
quartile × Age cat [Y] -6.35 -7.11 – -5.58 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		4.42 3.35 – 5.50 <0.001
Random Effects
σ2 486.47 486.47 477.81
Ï„00 106.64 Subj_ID 109.28 Subj_ID 109.29 Subj_ID
ICC 0.18 0.18 0.19
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.780 / 0.820 0.779 / 0.820 0.783 / 0.823 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

fit1 <- lmer(strength ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

fit2 <- lmer(strength ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

fit3 <- lmer(strength ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

tab_model(fit1,fit2,fit3)
  strength strength strength
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -12.57 -18.37 – -6.77 <0.001 -10.98 -18.76 – -3.20 0.006 -16.56 -24.41 – -8.72 <0.001
condition [Placebo] -1.50 -3.21 – 0.22 0.087 -1.50 -3.21 – 0.22 0.087 1.24 -1.03 – 3.51 0.283
quartile 42.60 42.16 – 43.04 <0.001 42.60 42.16 – 43.04 <0.001 45.56 44.98 – 46.15 <0.001
condition [Placebo] ×
quartile 		4.77 4.14 – 5.39 <0.001 4.77 4.14 – 5.39 <0.001 2.22 1.39 – 3.05 <0.001
Age cat [Y] -3.60 -15.20 – 7.99 0.542 9.08 -2.75 – 20.91 0.132
condition [Placebo] × Age
cat [Y] 		-6.22 -9.64 – -2.80 <0.001
quartile × Age cat [Y] -6.73 -7.62 – -5.85 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		5.79 4.54 – 7.04 <0.001
Random Effects
σ2 654.03 654.03 642.83
Ï„00 209.25 Subj_ID 214.90 Subj_ID 214.91 Subj_ID
ICC 0.24 0.25 0.25
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.748 / 0.809 0.747 / 0.810 0.750 / 0.813 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

Behavior - paired difference correlation analysis - Pearson

Overall

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanStrength = mean(strength)) -> data.aggr

colnames(data.aggr)
## [1] "threshold"    "Subj_ID"      "condition"    "meanStrength"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meanStrength_diff = diff(meanStrength), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meanStrength_diff, value_diff)

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanStrength = mean(strength)) -> data.aggr

colnames(data.aggr)
## [1] "threshold"    "Subj_ID"      "condition"    "name"         "meanStrength"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meanStrength_diff = diff(meanStrength), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meanStrength_diff, value_diff)

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Behavior - paired difference correlation analysis - Spearman

Overall

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanStrength = mean(strength)) -> data.aggr

colnames(data.aggr)
## [1] "threshold"    "Subj_ID"      "condition"    "meanStrength"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meanStrength_diff = diff(meanStrength), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meanStrength_diff, value_diff, method = "spearman")

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanStrength = mean(strength)) -> data.aggr

colnames(data.aggr)
## [1] "threshold"    "Subj_ID"      "condition"    "name"         "meanStrength"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meanStrength_diff = diff(meanStrength), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meanStrength_diff, value_diff, method = "spearman")

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Unthresholded

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanStrength = mean(strength)) -> data.aggr

anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr, within =c("condition"))

DT::datatable(nice(anova_summary))
emmip(anova_summary, ~condition )

ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanStrength,
  type             = "p",
  bf.message = FALSE,
  exact = FALSE
  )

Grouped by scale, condition, subject and age_cat

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition, Age_cat) %>% summarise(meanStrength = mean(strength)) -> data.aggr

anova_summary <- nice(aov_ez("Subj_ID", "meanStrength", data.aggr, within =c("condition"), between = "Age_cat"))

DT::datatable(anova_summary)
# 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
# 
# # aw <- aov_ez("Subj_ID", "meanDegree", data.aggr, within ="condition", between="Age_cat")
# # p_an <- afex_plot(aw, x = "condition", trace = "Age_cat")  + ggtitle("test") +   theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm"))
# 
# # plot anovas
# 
# d_nested <- data.aggr %>% ungroup() %>%
#   nest_by(threshold)
# 
# d_plots <- 
#   d_nested %>% 
#   mutate(aov = list(map(data,~aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>%
#   mutate(plot = list(map2(aov,threshold,~afex_plot(., x = "condition", trace = "Age_cat") +ggtitle(paste("threshold:",threshold)) + theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm")))))
# 
# d_plots$plot_one <- lapply(d_plots$plot, "[[", 1)
# 
# library(gridExtra)
# do.call(grid.arrange, c(d_plots$plot_one))

Grouped by scale, condition, subject and network:yeo

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition,name) %>% summarise(meanStrength = mean(strength)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name"))
DT::datatable(nice(main_anova_summary))
#emmip(main_anova_summary, condition~threshold)

# average across threshold
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name"))

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Grouped by scale, condition, subject, age_cat and network:yeo

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition,name,Age_cat) %>% summarise(meanStrength = mean(strength)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name"), between = "Age_cat")

DT::datatable(nice(main_anova_summary)) 
emmip(main_anova_summary, name~condition | Age_cat)

emmip(main_anova_summary, name~Age_cat)

main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name"),between = "Age_cat")

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
main_anova_summary <- aov_ez("Subj_ID", "meanStrength", data.aggr,within =c("condition","name"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))

Hubness?

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition) %>% mutate(quartile = ntile(strength, 4)) -> data.aggr

fit1 <- lmer(strength ~ condition*quartile + (1 | Subj_ID), data.aggr )

fit2 <- lmer(strength ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr )

fit3 <- lmer(strength ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr )

tab_model(fit1,fit2,fit3)
  strength strength strength
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 45.06 28.87 – 61.24 <0.001 46.00 23.98 – 68.02 <0.001 44.90 22.82 – 66.97 <0.001
condition [Placebo] 19.21 16.51 – 21.91 <0.001 19.21 16.51 – 21.91 <0.001 11.44 7.88 – 15.01 <0.001
quartile 39.91 39.22 – 40.61 <0.001 39.91 39.22 – 40.61 <0.001 42.53 41.61 – 43.45 <0.001
condition [Placebo] ×
quartile 		1.57 0.58 – 2.56 0.002 1.57 0.58 – 2.56 0.002 0.34 -0.97 – 1.64 0.611
Age cat [Y] -2.14 -35.22 – 30.94 0.899 0.36 -32.92 – 33.64 0.983
condition [Placebo] × Age
cat [Y] 		17.65 12.28 – 23.02 <0.001
quartile × Age cat [Y] -5.94 -7.33 – -4.55 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		2.80 0.84 – 4.77 0.005
Random Effects
σ2 1629.90 1629.90 1585.72
Ï„00 1680.65 Subj_ID 1752.58 Subj_ID 1752.64 Subj_ID
ICC 0.51 0.52 0.52
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.400 / 0.704 0.395 / 0.708 0.403 / 0.716 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

Behavior - paired difference correlation analysis

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanStrength = mean(strength)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meanStrength_diff = diff(meanStrength), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meanStrength_diff, value_diff)

# all correlations between *measures* and node_strength
datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanStrength = mean(strength)) -> data.aggr

colnames(data.aggr)
## [1] "threshold"    "Subj_ID"      "condition"    "name"         "meanStrength"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meanStrength_diff = diff(meanStrength), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meanStrength_diff, value_diff)

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Betweenness centrality - using BCT Toolbox’s betweenness_wei.m

Thresholded

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr, within =c("condition","threshold"))

DT::datatable(nice(anova_summary))
emmip(anova_summary, condition~threshold )

grouped_ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanbwc,
  grouping.var     = threshold,
  type             = "p",
  bf.message = FALSE,
  exact = FALSE
  )

Grouped by scale, condition, subject and age_cat

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition, Age_cat) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

anova_summary <-aov_ez("Subj_ID", "meanbwc", data.aggr, within =c("condition","threshold"), between = "Age_cat")

DT::datatable(nice(anova_summary))
# 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
# 
# # aw <- aov_ez("Subj_ID", "meanDegree", data.aggr, within ="condition", between="Age_cat")
# # p_an <- afex_plot(aw, x = "condition", trace = "Age_cat")  + ggtitle("test") +   theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm"))
# 
# # plot anovas
# 
# d_nested <- data.aggr %>% ungroup() %>%
#   nest_by(threshold)
# 
# d_plots <- 
#   d_nested %>% 
#   mutate(aov = list(map(data,~aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>%
#   mutate(plot = list(map2(aov,threshold,~afex_plot(., x = "condition", trace = "Age_cat") +ggtitle(paste("threshold:",threshold)) + theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm")))))
# 
# d_plots$plot_one <- lapply(d_plots$plot, "[[", 1)
# 
# library(gridExtra)
# do.call(grid.arrange, c(d_plots$plot_one))

Grouped by scale, condition, subject and network:yeo

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
#emmip(main_anova_summary, condition~threshold)

# average across threshold
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name"))

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=1
DT::datatable(nice(aov_ez("Subj_ID", "meanbwc", data.aggr %>% filter(threshold=="0.1"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.1"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=25
DT::datatable(nice(aov_ez("Subj_ID", "meanbwc", data.aggr %>% filter(threshold=="0.25"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.25"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=3.5
DT::datatable(nice(aov_ez("Subj_ID", "meanbwc", data.aggr %>% filter(threshold=="0.35"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.35"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Grouped by scale, condition, subject, age_cat and network:yeo

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name,Age_cat) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name","threshold"), between = "Age_cat")

DT::datatable(nice(main_anova_summary)) 
emmip(main_anova_summary, name~condition | Age_cat)

emmip(main_anova_summary, name~Age_cat)

main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name"),between = "Age_cat")

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name","threshold"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
# focus on threshold=1
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr %>% filter(threshold==0.1),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=25
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr %>% filter(threshold==0.25),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=3.5
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr %>% filter(threshold==0.35),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition",between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Hubness?

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% mutate(quartile = ntile(bwc, 4)) -> data.aggr

fit1 <- lmer(bwc ~ condition*quartile + (1 | Subj_ID), data.aggr )

fit2 <- lmer(bwc ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr )

fit3 <- lmer(bwc ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr )

tab_model(fit1,fit2,fit3)
  bwc bwc bwc
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -381.72 -396.33 – -367.10 <0.001 -382.24 -399.62 – -364.86 <0.001 -378.93 -398.66 – -359.20 <0.001
condition [Placebo] -38.91 -54.46 – -23.36 <0.001 -38.91 -54.46 – -23.36 <0.001 -42.13 -62.90 – -21.35 <0.001
quartile 316.73 312.71 – 320.75 <0.001 316.73 312.71 – 320.75 <0.001 313.61 308.24 – 318.99 <0.001
condition [Placebo] ×
quartile 		19.57 13.88 – 25.26 <0.001 19.57 13.88 – 25.26 <0.001 24.45 16.85 – 32.04 <0.001
Age cat [Y] 1.19 -19.67 – 22.05 0.911 -6.34 -36.08 – 23.40 0.676
condition [Placebo] × Age
cat [Y] 		7.31 -24.01 – 38.63 0.647
quartile × Age cat [Y] 7.09 -1.01 – 15.19 0.086
(condition [Placebo] ×
quartile) × Age cat [Y] 		-11.08 -22.54 – 0.37 0.058
Random Effects
σ2 161764.81 161764.81 161737.02
Ï„00 603.27 Subj_ID 631.98 Subj_ID 631.99 Subj_ID
ICC 0.00 0.00 0.00
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 61500 61500 61500
Marginal R2 / Conditional R2 0.451 / 0.453 0.451 / 0.453 0.451 / 0.453 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

fit1 <- lmer(bwc ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit2 <- lmer(bwc ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit3 <- lmer(bwc ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

tab_model(fit1,fit2,fit3)
  bwc bwc bwc
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -672.01 -703.68 – -640.34 <0.001 -677.59 -713.89 – -641.28 <0.001 -669.34 -711.86 – -626.82 <0.001
condition [Placebo] -60.06 -96.98 – -23.15 0.001 -60.06 -96.98 – -23.15 0.001 -67.28 -116.59 – -17.96 0.008
quartile 494.12 484.58 – 503.67 <0.001 494.12 484.58 – 503.67 <0.001 486.99 474.24 – 499.75 <0.001
condition [Placebo] ×
quartile 		33.43 19.93 – 46.93 <0.001 33.43 19.93 – 46.93 <0.001 43.98 25.94 – 62.01 <0.001
Age cat [Y] 12.67 -27.03 – 52.37 0.532 -6.08 -70.18 – 58.02 0.853
condition [Placebo] × Age
cat [Y] 		16.39 -57.95 – 90.74 0.666
quartile × Age cat [Y] 16.21 -3.02 – 35.44 0.098
(condition [Placebo] ×
quartile) × Age cat [Y] 		-23.96 -51.15 – 3.23 0.084
Random Effects
σ2 303896.44 303896.44 303775.06
Ï„00 2092.59 Subj_ID 2156.68 Subj_ID 2156.83 Subj_ID
ICC 0.01 0.01 0.01
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.516 / 0.520 0.516 / 0.520 0.517 / 0.520 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

fit1 <- lmer(bwc ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

fit2 <- lmer(bwc ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

fit3 <- lmer(bwc ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

tab_model(fit1,fit2,fit3)
  bwc bwc bwc
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -291.48 -304.47 – -278.49 <0.001 -288.34 -303.30 – -273.37 <0.001 -289.48 -306.89 – -272.07 <0.001
condition [Placebo] -24.15 -38.98 – -9.33 0.001 -24.15 -38.98 – -9.33 0.001 -25.28 -45.09 – -5.48 0.012
quartile 262.34 258.50 – 266.17 <0.001 262.34 258.50 – 266.17 <0.001 261.89 256.77 – 267.01 <0.001
condition [Placebo] ×
quartile 		11.62 6.20 – 17.04 <0.001 11.62 6.20 – 17.04 <0.001 13.87 6.63 – 21.12 <0.001
Age cat [Y] -7.14 -23.88 – 9.60 0.403 -4.55 -30.79 – 21.70 0.734
condition [Placebo] × Age
cat [Y] 		2.56 -27.30 – 32.42 0.867
quartile × Age cat [Y] 1.01 -6.71 – 8.73 0.797
(condition [Placebo] ×
quartile) × Age cat [Y] 		-5.13 -16.05 – 5.79 0.358
Random Effects
σ2 48995.08 48995.08 48993.03
Ï„00 383.96 Subj_ID 389.61 Subj_ID 389.61 Subj_ID
ICC 0.01 0.01 0.01
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.646 / 0.648 0.646 / 0.648 0.646 / 0.648 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

fit1 <- lmer(bwc ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

fit2 <- lmer(bwc ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

fit3 <- lmer(bwc ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

tab_model(fit1,fit2,fit3)
  bwc bwc bwc
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -181.67 -190.55 – -172.78 <0.001 -180.80 -190.92 – -170.67 <0.001 -177.97 -189.92 – -166.02 <0.001
condition [Placebo] -32.51 -43.10 – -21.92 <0.001 -32.51 -43.10 – -21.92 <0.001 -33.83 -47.98 – -19.67 <0.001
quartile 193.74 191.00 – 196.48 <0.001 193.74 191.00 – 196.48 <0.001 191.95 188.29 – 195.61 <0.001
condition [Placebo] ×
quartile 		13.66 9.79 – 17.53 <0.001 13.66 9.79 – 17.53 <0.001 15.49 10.31 – 20.66 <0.001
Age cat [Y] -1.97 -12.75 – 8.80 0.720 -8.39 -26.41 – 9.63 0.361
condition [Placebo] × Age
cat [Y] 		2.98 -18.35 – 24.32 0.784
quartile × Age cat [Y] 4.05 -1.47 – 9.57 0.150
(condition [Placebo] ×
quartile) × Age cat [Y] 		-4.16 -11.96 – 3.65 0.297
Random Effects
σ2 25022.28 25022.28 25020.04
Ï„00 148.85 Subj_ID 155.61 Subj_ID 155.61 Subj_ID
ICC 0.01 0.01 0.01
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.667 / 0.669 0.667 / 0.669 0.667 / 0.669 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

Behavior - paired difference correlation analysis

Overall

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meanbwc_diff = diff(meanbwc), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meanbwc_diff, value_diff)

# all correlations between *measures* and node_strength
DT::datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

colnames(data.aggr)
## [1] "threshold" "Subj_ID"   "condition" "name"      "meanbwc"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meanbwc_diff = diff(meanbwc), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meanbwc_diff, value_diff)

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Unthresholded

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr, within =c("condition"))

DT::datatable(nice(anova_summary))
emmip(anova_summary, ~condition )

ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanbwc,
  type             = "p",
  bf.message = FALSE,
  exact = FALSE
  )

Grouped by scale, condition, subject and age_cat

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition, Age_cat) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

anova_summary <- nice(aov_ez("Subj_ID", "meanbwc", data.aggr, within =c("condition"), between = "Age_cat"))

DT::datatable(anova_summary)
# 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
# 
# # aw <- aov_ez("Subj_ID", "meanDegree", data.aggr, within ="condition", between="Age_cat")
# # p_an <- afex_plot(aw, x = "condition", trace = "Age_cat")  + ggtitle("test") +   theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm"))
# 
# # plot anovas
# 
# d_nested <- data.aggr %>% ungroup() %>%
#   nest_by(threshold)
# 
# d_plots <- 
#   d_nested %>% 
#   mutate(aov = list(map(data,~aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>%
#   mutate(plot = list(map2(aov,threshold,~afex_plot(., x = "condition", trace = "Age_cat") +ggtitle(paste("threshold:",threshold)) + theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm")))))
# 
# d_plots$plot_one <- lapply(d_plots$plot, "[[", 1)
# 
# library(gridExtra)
# do.call(grid.arrange, c(d_plots$plot_one))

Grouped by scale, condition, subject and network:yeo

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition,name) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name"))
DT::datatable(nice(main_anova_summary))
#emmip(main_anova_summary, condition~threshold)

# average across threshold
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name"))

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Grouped by scale, condition, subject, age_cat and network:yeo

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition,name,Age_cat) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name"), between = "Age_cat")

DT::datatable(nice(main_anova_summary)) 
emmip(main_anova_summary, name~condition | Age_cat)

emmip(main_anova_summary, name~Age_cat)

main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name"),between = "Age_cat")

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
main_anova_summary <- aov_ez("Subj_ID", "meanbwc", data.aggr,within =c("condition","name"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))

Hubness?

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition) %>% mutate(quartile = ntile(bwc, 4)) -> data.aggr

fit1 <- lmer(bwc ~ condition*quartile + (1 | Subj_ID), data.aggr )

fit2 <- lmer(bwc ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr )

fit3 <- lmer(bwc ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr )

tab_model(fit1,fit2,fit3)
  bwc bwc bwc
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -3106.77 -3276.30 – -2937.23 <0.001 -3068.36 -3243.49 – -2893.24 <0.001 -3131.56 -3357.30 – -2905.82 <0.001
condition [Placebo] 117.39 -120.19 – 354.96 0.333 117.39 -120.19 – 354.96 0.333 -16.73 -334.10 – 300.65 0.918
quartile 1976.52 1915.08 – 2037.96 <0.001 1976.52 1915.08 – 2037.96 <0.001 1984.91 1902.83 – 2066.99 <0.001
condition [Placebo] ×
quartile 		-76.00 -162.89 – 10.90 0.086 -76.00 -162.89 – 10.90 0.086 11.64 -104.44 – 127.72 0.844
Age cat [Y] -87.28 -191.81 – 17.25 0.102 56.35 -283.97 – 396.67 0.746
condition [Placebo] × Age
cat [Y] 		304.80 -173.66 – 783.26 0.212
quartile × Age cat [Y] -19.06 -142.80 – 104.68 0.763
(condition [Placebo] ×
quartile) × Age cat [Y] 		-199.17 -374.16 – -24.18 0.026
Random Effects
σ2 12589317.97 12589317.98 12581484.04
Ï„00 3391.39 Subj_ID 2166.28 Subj_ID 2175.83 Subj_ID
ICC 0.00 0.00 0.00
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.272 / 0.272 0.272 / 0.272 0.272 / 0.273 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

Behavior - paired difference correlation analysis

Overall

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meanbwc_diff = diff(meanbwc), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meanbwc_diff, value_diff)

# all correlations between *measures* and node_strength
DT::datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

colnames(data.aggr)
## [1] "threshold" "Subj_ID"   "condition" "name"      "meanbwc"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meanbwc_diff = diff(meanbwc), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meanbwc_diff, value_diff)

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Global efficiency - using BCT Toolbox’s efficiency_wei.m

Thresholded

Behavior - paired difference correlation analysis - Pearson

Overall

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meangeff = mean(geff)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meangeff_diff = diff(meangeff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meangeff_diff, value_diff)

# all correlations between *measures* and node_strength
DT::datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,Network) %>% summarise(meangeff = mean(geff)) %>% mutate(name=Network) -> data.aggr

colnames(data.aggr)
## [1] "threshold" "Subj_ID"   "condition" "Network"   "meangeff"  "name"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meangeff_diff = diff(meangeff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meangeff_diff, value_diff)

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Behavior - paired difference correlation analysis - Spearman

Overall

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meangeff = mean(geff)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meangeff_diff = diff(meangeff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meangeff_diff, value_diff, method = "spearman")

# all correlations between *measures* and node_strength
DT::datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,Network) %>% summarise(meangeff = mean(geff)) %>% mutate(name=Network) -> data.aggr

colnames(data.aggr)
## [1] "threshold" "Subj_ID"   "condition" "Network"   "meangeff"  "name"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meangeff_diff = diff(meangeff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meangeff_diff, value_diff, method = "spearman")

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meangeff = mean(geff)) -> data.aggr

anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr, within =c("condition","threshold"))

DT::datatable(nice(anova_summary))
emmip(anova_summary, condition~threshold )

grouped_ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meangeff,
  grouping.var     = threshold,
  type             = "p",
  bf.message = FALSE,
  exact = FALSE
  )

Grouped by scale, condition, subject and age_cat

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition, Age_cat) %>% summarise(meangeff = mean(geff)) -> data.aggr

anova_summary <-aov_ez("Subj_ID", "meangeff", data.aggr, within =c("condition","threshold"), between = "Age_cat")

DT::datatable(nice(anova_summary))
# 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
# 
# # aw <- aov_ez("Subj_ID", "meanDegree", data.aggr, within ="condition", between="Age_cat")
# # p_an <- afex_plot(aw, x = "condition", trace = "Age_cat")  + ggtitle("test") +   theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm"))
# 
# # plot anovas
# 
# d_nested <- data.aggr %>% ungroup() %>%
#   nest_by(threshold)
# 
# d_plots <- 
#   d_nested %>% 
#   mutate(aov = list(map(data,~aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>%
#   mutate(plot = list(map2(aov,threshold,~afex_plot(., x = "condition", trace = "Age_cat") +ggtitle(paste("threshold:",threshold)) + theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm")))))
# 
# d_plots$plot_one <- lapply(d_plots$plot, "[[", 1)
# 
# library(gridExtra)
# do.call(grid.arrange, c(d_plots$plot_one))

Grouped by scale, condition, subject and network:yeo

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,Network) %>% summarise(meangeff = mean(geff)) %>% mutate(name=Network) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
#emmip(main_anova_summary, condition~threshold)

# average across threshold
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name"))

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=1
DT::datatable(nice(aov_ez("Subj_ID", "meangeff", data.aggr %>% filter(threshold=="0.1"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.1"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=25
DT::datatable(nice(aov_ez("Subj_ID", "meangeff", data.aggr %>% filter(threshold=="0.25"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.25"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=3.5
DT::datatable(nice(aov_ez("Subj_ID", "meangeff", data.aggr %>% filter(threshold=="0.35"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.35"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Grouped by scale, condition, subject, age_cat and network:yeo

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,Network,Age_cat) %>% summarise(meangeff = mean(geff)) %>% mutate(name=Network) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name","threshold"), between = "Age_cat")

DT::datatable(nice(main_anova_summary)) 
emmip(main_anova_summary, name~condition | Age_cat)

emmip(main_anova_summary, name~Age_cat)

main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name"),between = "Age_cat")

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name","threshold"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
# focus on threshold=1
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr %>% filter(threshold==0.1),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=25
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr %>% filter(threshold==0.25),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=3.5
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr %>% filter(threshold==0.35),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition",between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Hubness?

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% mutate(quartile = ntile(geff, 4)) -> data.aggr

fit1 <- lmer(geff ~ condition*quartile + (1 | Subj_ID), data.aggr )

fit2 <- lmer(geff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr )

fit3 <- lmer(geff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr )

tab_model(fit1,fit2,fit3)
  geff geff geff
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.18 0.15 – 0.22 <0.001 0.18 0.14 – 0.22 <0.001 0.17 0.12 – 0.21 <0.001
condition [Placebo] 0.00 -0.03 – 0.03 0.955 0.00 -0.03 – 0.03 0.955 0.01 -0.03 – 0.06 0.600
quartile 0.14 0.13 – 0.15 <0.001 0.14 0.13 – 0.15 <0.001 0.15 0.14 – 0.16 <0.001
condition [Placebo] ×
quartile 		0.02 0.01 – 0.04 <0.001 0.02 0.01 – 0.04 <0.001 0.02 0.00 – 0.03 0.027
Age cat [Y] 0.00 -0.04 – 0.05 0.868 0.04 -0.02 – 0.11 0.198
condition [Placebo] × Age
cat [Y] 		-0.02 -0.09 – 0.04 0.468
quartile × Age cat [Y] -0.02 -0.03 – -0.00 0.045
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.01 -0.01 – 0.04 0.262
Random Effects
σ2 0.01 0.01 0.01
Ï„00 0.00 Subj_ID 0.00 Subj_ID 0.00 Subj_ID
ICC 0.18 0.19 0.19
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 1200 1200 1200
Marginal R2 / Conditional R2 0.638 / 0.703 0.636 / 0.704 0.637 / 0.705 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

fit1 <- lmer(geff ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit2 <- lmer(geff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit3 <- lmer(geff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

tab_model(fit1,fit2,fit3)
  geff geff geff
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.06 0.02 – 0.10 0.004 0.05 0.01 – 0.09 0.020 0.03 -0.02 – 0.09 0.218
condition [Placebo] -0.01 -0.06 – 0.04 0.716 -0.01 -0.06 – 0.04 0.716 0.00 -0.07 – 0.08 0.926
quartile 0.17 0.15 – 0.18 <0.001 0.17 0.15 – 0.18 <0.001 0.17 0.15 – 0.19 <0.001
condition [Placebo] ×
quartile 		0.02 0.00 – 0.04 0.015 0.02 0.00 – 0.04 0.015 0.02 -0.01 – 0.04 0.180
Age cat [Y] 0.02 -0.02 – 0.05 0.341 0.06 -0.03 – 0.14 0.174
condition [Placebo] × Age
cat [Y] 		-0.03 -0.14 – 0.08 0.583
quartile × Age cat [Y] -0.02 -0.05 – 0.01 0.235
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.01 -0.03 – 0.05 0.474
Random Effects
σ2 0.01 0.01 0.01
Ï„00 0.00 Subj_ID 0.00 Subj_ID 0.00 Subj_ID
ICC 0.08 0.08 0.08
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 400 400 400
Marginal R2 / Conditional R2 0.745 / 0.766 0.745 / 0.767 0.745 / 0.766 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

fit1 <- lmer(geff ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

fit2 <- lmer(geff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

fit3 <- lmer(geff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

tab_model(fit1,fit2,fit3)
  geff geff geff
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.22 0.18 – 0.27 <0.001 0.22 0.17 – 0.27 <0.001 0.20 0.14 – 0.26 <0.001
condition [Placebo] 0.00 -0.05 – 0.05 0.943 0.00 -0.05 – 0.05 0.943 0.02 -0.05 – 0.09 0.647
quartile 0.13 0.12 – 0.15 <0.001 0.13 0.12 – 0.15 <0.001 0.14 0.13 – 0.16 <0.001
condition [Placebo] ×
quartile 		0.03 0.01 – 0.04 0.009 0.03 0.01 – 0.04 0.009 0.02 -0.01 – 0.04 0.166
Age cat [Y] 0.00 -0.05 – 0.05 0.978 0.05 -0.04 – 0.14 0.258
condition [Placebo] × Age
cat [Y] 		-0.03 -0.14 – 0.07 0.542
quartile × Age cat [Y] -0.02 -0.05 – 0.01 0.116
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.02 -0.02 – 0.05 0.395
Random Effects
σ2 0.01 0.01 0.01
Ï„00 0.00 Subj_ID 0.00 Subj_ID 0.00 Subj_ID
ICC 0.22 0.23 0.23
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 400 400 400
Marginal R2 / Conditional R2 0.657 / 0.733 0.654 / 0.734 0.655 / 0.734 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

fit1 <- lmer(geff ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

fit2 <- lmer(geff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

fit3 <- lmer(geff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

tab_model(fit1,fit2,fit3)
  geff geff geff
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.27 0.23 – 0.31 <0.001 0.27 0.22 – 0.32 <0.001 0.26 0.20 – 0.32 <0.001
condition [Placebo] 0.01 -0.04 – 0.06 0.656 0.01 -0.04 – 0.06 0.656 0.02 -0.05 – 0.08 0.634
quartile 0.12 0.11 – 0.14 <0.001 0.12 0.11 – 0.14 <0.001 0.13 0.11 – 0.15 <0.001
condition [Placebo] ×
quartile 		0.02 0.01 – 0.04 0.010 0.02 0.01 – 0.04 0.010 0.02 -0.01 – 0.04 0.123
Age cat [Y] -0.01 -0.06 – 0.05 0.847 0.02 -0.07 – 0.11 0.632
condition [Placebo] × Age
cat [Y] 		-0.01 -0.11 – 0.09 0.828
quartile × Age cat [Y] -0.01 -0.04 – 0.01 0.280
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.01 -0.03 – 0.05 0.563
Random Effects
σ2 0.01 0.01 0.01
Ï„00 0.00 Subj_ID 0.00 Subj_ID 0.00 Subj_ID
ICC 0.30 0.31 0.31
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 400 400 400
Marginal R2 / Conditional R2 0.621 / 0.736 0.617 / 0.737 0.617 / 0.737 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

Unthresholded

Behavior - paired difference correlation analysis

OVerall

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meangeff = mean(geff)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meangeff_diff = diff(meangeff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meangeff_diff, value_diff)

# all correlations between *measures* and node_strength
DT::datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,Network) %>% summarise(meangeff = mean(geff)) %>% mutate(name=Network) -> data.aggr

colnames(data.aggr)
## [1] "threshold" "Subj_ID"   "condition" "Network"   "meangeff"  "name"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanStrength =as.numeric(meanStrength ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meangeff_diff = diff(meangeff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meangeff_diff, value_diff)

# all correlations between *measures* and node_strength
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meangeff = mean(geff)) -> data.aggr

anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr, within =c("condition"))

DT::datatable(nice(anova_summary))
emmip(anova_summary, ~condition )

ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meangeff,
  type             = "p",
  bf.message = FALSE,
  exact = FALSE
  )

Grouped by scale, condition, subject and age_cat

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition, Age_cat) %>% summarise(meangeff = mean(geff)) -> data.aggr

anova_summary <- nice(aov_ez("Subj_ID", "meangeff", data.aggr, within =c("condition"), between = "Age_cat"))

DT::datatable(anova_summary)
# 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
# 
# # aw <- aov_ez("Subj_ID", "meanDegree", data.aggr, within ="condition", between="Age_cat")
# # p_an <- afex_plot(aw, x = "condition", trace = "Age_cat")  + ggtitle("test") +   theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm"))
# 
# # plot anovas
# 
# d_nested <- data.aggr %>% ungroup() %>%
#   nest_by(threshold)
# 
# d_plots <- 
#   d_nested %>% 
#   mutate(aov = list(map(data,~aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>%
#   mutate(plot = list(map2(aov,threshold,~afex_plot(., x = "condition", trace = "Age_cat") +ggtitle(paste("threshold:",threshold)) + theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm")))))
# 
# d_plots$plot_one <- lapply(d_plots$plot, "[[", 1)
# 
# library(gridExtra)
# do.call(grid.arrange, c(d_plots$plot_one))

Grouped by scale, condition, subject and network:yeo

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition,Network) %>% summarise(meangeff = mean(geff)) %>% mutate(name=Network) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name"))
DT::datatable(nice(main_anova_summary))
#emmip(main_anova_summary, condition~threshold)

# average across threshold
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name"))

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Grouped by scale, condition, subject, age_cat and network:yeo

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition,Network,Age_cat) %>% summarise(meangeff = mean(geff)) %>% mutate(name=Network) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name"), between = "Age_cat")

DT::datatable(nice(main_anova_summary)) 
emmip(main_anova_summary, name~condition | Age_cat)

emmip(main_anova_summary, name~Age_cat)

main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name"),between = "Age_cat")

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","name"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))

Hubness?

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition) %>% mutate(quartile = ntile(geff, 4)) -> data.aggr

fit1 <- lmer(geff ~ condition*quartile + (1 | Subj_ID), data.aggr )

fit2 <- lmer(geff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr )

fit3 <- lmer(geff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr )

tab_model(fit1,fit2,fit3)
  geff geff geff
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.19 -0.28 – -0.10 <0.001 -0.19 -0.29 – -0.09 <0.001 -0.21 -0.33 – -0.10 <0.001
condition [Placebo] -0.00 -0.10 – 0.09 0.932 -0.00 -0.10 – 0.09 0.932 -0.04 -0.17 – 0.09 0.546
quartile 0.13 0.11 – 0.16 <0.001 0.13 0.11 – 0.16 <0.001 0.15 0.12 – 0.18 <0.001
condition [Placebo] ×
quartile 		0.04 0.01 – 0.08 0.019 0.04 0.01 – 0.08 0.019 0.05 -0.00 – 0.09 0.060
Age cat [Y] 0.00 -0.11 – 0.12 0.961 0.06 -0.12 – 0.23 0.513
condition [Placebo] × Age
cat [Y] 		0.08 -0.11 – 0.28 0.417
quartile × Age cat [Y] -0.04 -0.09 – 0.02 0.170
(condition [Placebo] ×
quartile) × Age cat [Y] 		-0.01 -0.08 – 0.06 0.861
Random Effects
σ2 0.04 0.04 0.04
Ï„00 0.02 Subj_ID 0.02 Subj_ID 0.02 Subj_ID
ICC 0.31 0.32 0.32
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 400 400 400
Marginal R2 / Conditional R2 0.362 / 0.558 0.358 / 0.562 0.365 / 0.568 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

Local efficiency - using BCT Toolbox’s efficiency_wei.m

Thresholded

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanleff = mean(leff)) -> data.aggr

anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr, within =c("condition","threshold"))


emmip(anova_summary, condition~threshold )

grouped_ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanleff,
  grouping.var     = threshold,
  type             = "p",
  bf.message = FALSE,
  exact = FALSE
  )

DT::datatable(nice(anova_summary))

Grouped by scale, condition, subject and age_cat

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition, Age_cat) %>% summarise(meanleff = mean(leff)) -> data.aggr

anova_summary <-aov_ez("Subj_ID", "meanleff", data.aggr, within =c("condition","threshold"), between = "Age_cat")


# 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
# 
# # aw <- aov_ez("Subj_ID", "meanDegree", data.aggr, within ="condition", between="Age_cat")
# # p_an <- afex_plot(aw, x = "condition", trace = "Age_cat")  + ggtitle("test") +   theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm"))
# 
# # plot anovas
# 
# d_nested <- data.aggr %>% ungroup() %>%
#   nest_by(threshold)
# 
# d_plots <- 
#   d_nested %>% 
#   mutate(aov = list(map(data,~aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>%
#   mutate(plot = list(map2(aov,threshold,~afex_plot(., x = "condition", trace = "Age_cat") +ggtitle(paste("threshold:",threshold)) + theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm")))))
# 
# d_plots$plot_one <- lapply(d_plots$plot, "[[", 1)
# 
# library(gridExtra)
# do.call(grid.arrange, c(d_plots$plot_one))
DT::datatable(nice(anova_summary))

Grouped by scale, condition, subject and network:yeo

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>%   mutate(ROI = gsub("Var", "ROI_", ROI)) %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanleff = mean(leff)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
#emmip(main_anova_summary, condition~threshold)

# average across threshold
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name"))

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=1
DT::datatable(nice(aov_ez("Subj_ID", "meanleff", data.aggr %>% filter(threshold=="0.1"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.1"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=25
DT::datatable(nice(aov_ez("Subj_ID", "meanleff", data.aggr %>% filter(threshold=="0.25"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.25"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=3.5
DT::datatable(nice(aov_ez("Subj_ID", "meanleff", data.aggr %>% filter(threshold=="0.35"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X0.35"))
EMM.unadj <- summary(pairs(EMM), by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(pairs(EMM), by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Grouped by scale, condition, subject, age_cat and network:yeo

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>%   mutate(ROI = gsub("Var", "ROI_", ROI))%>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name,Age_cat) %>% summarise(meanleff = mean(leff)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name","threshold"), between = "Age_cat")

DT::datatable(nice(main_anova_summary)) 
emmip(main_anova_summary, name~condition | Age_cat)

emmip(main_anova_summary, name~Age_cat)

main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name"),between = "Age_cat")

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name","threshold"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
# focus on threshold=1
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr %>% filter(threshold==0.1),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=25
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr %>% filter(threshold==0.25),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# focus on threshold=3.5
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr %>% filter(threshold==0.35),within =c("name","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name) %>% contrast(interaction = c( "consec", "consec"))
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)
# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))
# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition",between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Hubness?

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% mutate(quartile = ntile(leff, 4)) -> data.aggr

fit1 <- lmer(leff ~ condition*quartile + (1 | Subj_ID), data.aggr )

fit2 <- lmer(leff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr )

fit3 <- lmer(leff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr )

tab_model(fit1,fit2,fit3)
  leff leff leff
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.45 0.42 – 0.49 <0.001 0.46 0.41 – 0.50 <0.001 0.44 0.40 – 0.48 <0.001
condition [Placebo] 0.05 0.05 – 0.06 <0.001 0.05 0.05 – 0.06 <0.001 0.06 0.05 – 0.07 <0.001
quartile 0.08 0.08 – 0.08 <0.001 0.08 0.08 – 0.08 <0.001 0.09 0.09 – 0.10 <0.001
condition [Placebo] ×
quartile 		0.01 0.00 – 0.01 <0.001 0.01 0.00 – 0.01 <0.001 -0.00 -0.00 – 0.00 0.730
Age cat [Y] -0.01 -0.08 – 0.06 0.760 0.03 -0.03 – 0.10 0.337
condition [Placebo] × Age
cat [Y] 		-0.01 -0.03 – 0.00 0.090
quartile × Age cat [Y] -0.02 -0.03 – -0.02 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.02 0.01 – 0.02 <0.001
Random Effects
σ2 0.03 0.03 0.03
Ï„00 0.01 Subj_ID 0.01 Subj_ID 0.01 Subj_ID
ICC 0.18 0.19 0.19
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 61500 61500 61500
Marginal R2 / Conditional R2 0.223 / 0.365 0.223 / 0.369 0.226 / 0.372 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

fit1 <- lmer(leff ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit2 <- lmer(leff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit3 <- lmer(leff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

tab_model(fit1,fit2,fit3)
  leff leff leff
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.37 0.34 – 0.40 <0.001 0.37 0.33 – 0.40 <0.001 0.32 0.28 – 0.35 <0.001
condition [Placebo] 0.05 0.03 – 0.06 <0.001 0.05 0.03 – 0.06 <0.001 0.08 0.06 – 0.09 <0.001
quartile 0.14 0.14 – 0.15 <0.001 0.14 0.14 – 0.15 <0.001 0.17 0.16 – 0.17 <0.001
condition [Placebo] ×
quartile 		0.01 0.00 – 0.02 <0.001 0.01 0.00 – 0.02 <0.001 -0.00 -0.01 – 0.00 0.331
Age cat [Y] 0.01 -0.05 – 0.06 0.738 0.12 0.07 – 0.18 <0.001
condition [Placebo] × Age
cat [Y] 		-0.06 -0.09 – -0.03 <0.001
quartile × Age cat [Y] -0.05 -0.06 – -0.04 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.03 0.02 – 0.04 <0.001
Random Effects
σ2 0.05 0.05 0.04
Ï„00 0.00 Subj_ID 0.00 Subj_ID 0.00 Subj_ID
ICC 0.09 0.09 0.10
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.370 / 0.427 0.369 / 0.428 0.374 / 0.434 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

fit1 <- lmer(leff ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

fit2 <- lmer(leff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

fit3 <- lmer(leff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.25))

tab_model(fit1,fit2,fit3)
  leff leff leff
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.51 0.47 – 0.54 <0.001 0.52 0.47 – 0.56 <0.001 0.51 0.46 – 0.56 <0.001
condition [Placebo] 0.05 0.05 – 0.06 <0.001 0.05 0.05 – 0.06 <0.001 0.05 0.04 – 0.06 <0.001
quartile 0.06 0.06 – 0.06 <0.001 0.06 0.06 – 0.06 <0.001 0.06 0.06 – 0.07 <0.001
condition [Placebo] ×
quartile 		0.01 0.00 – 0.01 <0.001 0.01 0.00 – 0.01 <0.001 -0.00 -0.00 – 0.00 0.965
Age cat [Y] -0.02 -0.09 – 0.05 0.595 -0.00 -0.08 – 0.07 0.900
condition [Placebo] × Age
cat [Y] 		0.00 -0.01 – 0.02 0.704
quartile × Age cat [Y] -0.01 -0.02 – -0.01 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.01 0.01 – 0.02 <0.001
Random Effects
σ2 0.02 0.02 0.02
Ï„00 0.01 Subj_ID 0.01 Subj_ID 0.01 Subj_ID
ICC 0.32 0.33 0.33
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.194 / 0.451 0.195 / 0.457 0.199 / 0.461 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

fit1 <- lmer(leff ~ condition*quartile + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

fit2 <- lmer(leff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

fit3 <- lmer(leff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr %>% filter(threshold==0.35))

tab_model(fit1,fit2,fit3)
  leff leff leff
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.49 0.45 – 0.52 <0.001 0.50 0.45 – 0.54 <0.001 0.50 0.45 – 0.54 <0.001
condition [Placebo] 0.05 0.05 – 0.06 <0.001 0.05 0.05 – 0.06 <0.001 0.04 0.03 – 0.05 <0.001
quartile 0.05 0.04 – 0.05 <0.001 0.05 0.04 – 0.05 <0.001 0.05 0.05 – 0.05 <0.001
condition [Placebo] ×
quartile 		0.01 0.00 – 0.01 <0.001 0.01 0.00 – 0.01 <0.001 0.00 -0.00 – 0.01 0.311
Age cat [Y] -0.02 -0.09 – 0.05 0.576 -0.02 -0.10 – 0.05 0.575
condition [Placebo] × Age
cat [Y] 		0.02 0.01 – 0.04 0.001
quartile × Age cat [Y] -0.01 -0.01 – -0.01 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.01 0.00 – 0.01 0.002
Random Effects
σ2 0.01 0.01 0.01
Ï„00 0.01 Subj_ID 0.01 Subj_ID 0.01 Subj_ID
ICC 0.42 0.43 0.43
N 25 Subj_ID 25 Subj_ID 25 Subj_ID
Observations 20500 20500 20500
Marginal R2 / Conditional R2 0.172 / 0.521 0.174 / 0.528 0.180 / 0.534 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

Behavior - paired difference correlation analysis - Pearson

Overall

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanleff = mean(leff)) -> data.aggr

colnames(data.aggr)
## [1] "threshold" "Subj_ID"   "condition" "meanleff"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanleff =as.numeric(meanleff ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meanleff_diff = diff(meanleff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meanleff_diff, value_diff)

# all correlations between *measures* and node_leff
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>%   mutate(ROI = gsub("Var", "ROI_", ROI))%>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanleff = mean(leff)) -> data.aggr

colnames(data.aggr)
## [1] "threshold" "Subj_ID"   "condition" "name"      "meanleff"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanleff =as.numeric(meanleff ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meanleff_diff = diff(meanleff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meanleff_diff, value_diff)

# all correlations between *measures* and node_leff
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Behavior - paired difference correlation analysis - Spearman

Overall

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanleff = mean(leff)) -> data.aggr

colnames(data.aggr)
## [1] "threshold" "Subj_ID"   "condition" "meanleff"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanleff =as.numeric(meanleff ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meanleff_diff = diff(meanleff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meanleff_diff, value_diff, method = "spearman")

# all correlations between *measures* and node_leff
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>%   mutate(ROI = gsub("Var", "ROI_", ROI))%>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanleff = mean(leff)) -> data.aggr

colnames(data.aggr)
## [1] "threshold" "Subj_ID"   "condition" "name"      "meanleff"
data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanleff =as.numeric(meanleff ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meanleff_diff = diff(meanleff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meanleff_diff, value_diff, method = "spearman")

# all correlations between *measures* and node_leff
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)

Unthresholded

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanleff = mean(leff)) -> data.aggr

anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr, within =c("condition"))

DT::datatable(nice(anova_summary))

emmip(anova_summary, ~condition )

ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanleff,
  type             = "p",
  bf.message = FALSE,
  exact = FALSE
  )

Grouped by scale, condition, subject and age_cat

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition, Age_cat) %>% summarise(meanleff = mean(leff)) -> data.aggr

anova_summary <- nice(aov_ez("Subj_ID", "meanleff", data.aggr, within =c("condition"), between = "Age_cat"))

DT::datatable(anova_summary)

# 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)
# 
# # aw <- aov_ez("Subj_ID", "meanDegree", data.aggr, within ="condition", between="Age_cat")
# # p_an <- afex_plot(aw, x = "condition", trace = "Age_cat")  + ggtitle("test") +   theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm"))
# 
# # plot anovas
# 
# d_nested <- data.aggr %>% ungroup() %>%
#   nest_by(threshold)
# 
# d_plots <- 
#   d_nested %>% 
#   mutate(aov = list(map(data,~aov_ez("Subj_ID", "meanDegree", data, within ="condition", between="Age_cat")))) %>%
#   mutate(plot = list(map2(aov,threshold,~afex_plot(., x = "condition", trace = "Age_cat") +ggtitle(paste("threshold:",threshold)) + theme(legend.position = c(0.87, 0.15),legend.key.size = unit(0.3, "cm")))))
# 
# d_plots$plot_one <- lapply(d_plots$plot, "[[", 1)
# 
# library(gridExtra)
# do.call(grid.arrange, c(d_plots$plot_one))

Grouped by scale, condition, subject and network:yeo

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>%   mutate(ROI = gsub("Var", "ROI_", ROI))%>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition,name) %>% summarise(meanleff = mean(leff)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name"))
DT::datatable(nice(main_anova_summary))

#emmip(main_anova_summary, condition~threshold)

# average across threshold
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name"))

DT::datatable(nice(main_anova_summary))

EMM <- emmeans(main_anova_summary, ~  condition | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))

# data.aggr %>% ungroup() %>%
#   nest_by(threshold,name) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanDegree", data, within ="condition")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Grouped by scale, condition, subject, age_cat and network:yeo

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff")%>%   mutate(ROI = gsub("Var", "ROI_", ROI)) %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition,name,Age_cat) %>% summarise(meanleff = mean(leff)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name"), between = "Age_cat")

DT::datatable(nice(main_anova_summary)) 

emmip(main_anova_summary, name~condition | Age_cat)

emmip(main_anova_summary, name~Age_cat)

main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name"),between = "Age_cat")

DT::datatable(nice(main_anova_summary))

EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | name)  %>% contrast(interaction = c( "consec", "consec"))
  
EMM.unadj <- summary(EMM, by = NULL, adjust = "none") %>% as.data.frame() %>% dplyr::rename(.,p.value.unadj = p.value)
EMM.fdr <- summary(EMM, by = NULL, adjust = "fdr") %>% as.data.frame() %>% dplyr::rename(.,p.value.fdr = p.value)

# corrected and uncorrected
DT::datatable(left_join(EMM.unadj,EMM.fdr))

main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))

Hubness?

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(Subj_ID,condition) %>% mutate(quartile = ntile(leff, 4)) -> data.aggr

fit1 <- lmer(leff ~ condition*quartile + (1 | Subj_ID), data.aggr )

fit2 <- lmer(leff ~ condition*quartile + Age_cat + (1 | Subj_ID), data.aggr )

fit3 <- lmer(leff ~ condition*quartile*Age_cat + (1 | Subj_ID), data.aggr )

tab_model(fit1,fit2,fit3)

plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

plot_model(fit3, type = "pred", terms = c("condition","Age_cat"))

plot_model(fit3, type = "pred", terms = c("quartile[all]","condition","Age_cat"))

Behavior - paired difference correlation analysis

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanleff = mean(leff)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanleff =as.numeric(meanleff ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meanleff_diff = diff(meanleff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold) %>% rstatix::cor_test(meanleff_diff, value_diff)

# all correlations between *measures* and node_leff
datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>%   mutate(ROI = gsub("Var", "ROI_", ROI)) %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanleff = mean(leff)) -> data.aggr

colnames(data.aggr)

data.aggr %>% left_join(.,variables_ext) %>% filter(!Subj_ID %in% c(122025,421002)) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

# #temporary
# # Fill NAs in each column with random numbers
# data.aggr.withvars_filled <- apply(data.aggr.withvars, 2, function(x) {
#   ifelse(is.na(x), sample(na.omit(x), sum(is.na(x)), replace = TRUE), x)
# }) %>% as.data.frame() %>% mutate(meanleff =as.numeric(meanleff ), value = as.numeric(value) )


# Group the data by the cyl column and calculate the paired differences of mpg
data.aggr.withvars_diff <- data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold,name) %>%
  summarise(meanleff_diff = diff(meanleff), value_diff=diff(value))

cor_results <- data.aggr.withvars_diff %>% group_by(measure,threshold,name) %>% rstatix::cor_test(meanleff_diff, value_diff)

# all correlations between *measures* and node_leff
#DT::datatable(cor_results)
#knitr::kable(cor_results)

datatable(cor_results)