## 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)

IL10 analysis

library(readxl)
library(tidyr)

IL_10 <- read_excel("IL_10.xlsx") 
IL_10[IL_10 == 0.39] <- NA


IL_10_long <- pivot_longer(IL_10, 
                       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_10") %>%
  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_10 ~ time_since_inj_h*Condition + (1 | ID), data = IL_10_long)
plot_model(m0, type = "pred",terms = c("time_since_inj_h","Condition"),show.data = TRUE)

#linear mixed model - 27 subs
m1<- lmer(IL_10 ~ time_since_inj_h*Condition + (1 | ID), data = IL_10_long %>% filter(!ID %in% c("008","022","023")))
plot_model(m1, type = "pred",terms = c("time_since_inj_h","Condition"),show.data = TRUE)

#linear mixed model - 25 subs
m2 <- lmer(IL_10 ~ time_since_inj_h*Condition + (1 | ID), data = IL_10_long %>% filter(!ID %in% c("008","022","023","025","002")))
plot_model(m2, type = "pred",terms = c("time_since_inj_h","Condition"),show.data = TRUE)

tab_model(m0,m1,m2)
  IL_10 IL_10 IL_10
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 1.56 0.76 – 2.36 <0.001 1.56 0.76 – 2.36 <0.001 1.61 0.80 – 2.43 <0.001
time since inj h 0.03 -0.12 – 0.19 0.640 0.03 -0.12 – 0.19 0.640 0.02 -0.13 – 0.17 0.782
Condition [Plac] -1.28 -3.87 – 1.31 0.322 -1.28 -3.87 – 1.31 0.322 -1.33 -3.91 – 1.24 0.298
time since inj h ×
Condition [Plac] 		0.33 -0.16 – 0.82 0.179 0.33 -0.16 – 0.82 0.179 0.34 -0.14 – 0.83 0.160
Random Effects
σ2 0.87 0.87 0.84
Ï„00 0.16 ID 0.16 ID 0.21 ID
ICC 0.16 0.16 0.20
N 17 ID 17 ID 16 ID
Observations 35 35 33
Marginal R2 / Conditional R2 0.102 / 0.245 0.102 / 0.245 0.102 / 0.280

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/", 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  21
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/", 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/", 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/sub-0122025_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0122026_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0122027_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0222028_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0222029_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0322030_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0322031_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0322032_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0421002_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0422033_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0422034_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0422035_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0521003_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0522036_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0721005_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0821006_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0821007_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0921009_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0921010_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0921012_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0921014_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0921015_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0921016_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-0921017_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-1021019_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/sub-1121020_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                     0 
## conmats_filtered_tedana_bp/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.5123e+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 = 2.2357e+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.5573e+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/", 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)`:TNFd_Z, 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)

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) %>% 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) %>% 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) %>% 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) %>% 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) %>% 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) -23.36 -30.83 – -15.88 <0.001 -27.04 -37.37 – -16.71 <0.001 -24.05 -34.43 – -13.68 <0.001
condition [Placebo] -2.31 -3.80 – -0.82 0.002 -2.31 -3.80 – -0.82 0.002 -0.39 -2.46 – 1.68 0.711
quartile 38.63 38.25 – 39.02 <0.001 38.63 38.25 – 39.02 <0.001 37.81 37.28 – 38.35 <0.001
condition [Placebo] ×
quartile 		3.65 3.10 – 4.19 <0.001 3.65 3.10 – 4.19 <0.001 2.12 1.37 – 2.88 <0.001
Age cat [Y] 7.64 -7.18 – 22.46 0.312 1.44 -13.51 – 16.40 0.850
condition [Placebo] × Age
cat [Y] 		-3.99 -6.97 – -1.02 0.009
quartile × Age cat [Y] 1.70 0.93 – 2.47 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		3.17 2.08 – 4.26 <0.001
Random Effects
σ2 1604.35 1604.35 1599.31
Ï„00 384.85 Subj_ID 384.53 Subj_ID 384.53 Subj_ID
ICC 0.19 0.19 0.19
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 66420 66420 66420
Marginal R2 / Conditional R2 0.509 / 0.604 0.511 / 0.605 0.512 / 0.607 
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) -27.46 -30.69 – -24.22 <0.001 -29.00 -33.43 – -24.57 <0.001 -28.39 -32.88 – -23.90 <0.001
condition [Placebo] -2.55 -3.69 – -1.41 <0.001 -2.55 -3.69 – -1.41 <0.001 -2.13 -3.71 – -0.54 0.008
quartile 26.48 26.19 – 26.78 <0.001 26.48 26.19 – 26.78 <0.001 26.37 25.96 – 26.78 <0.001
condition [Placebo] ×
quartile 		2.34 1.92 – 2.76 <0.001 2.34 1.92 – 2.76 <0.001 1.92 1.34 – 2.50 <0.001
Age cat [Y] 3.20 -3.09 – 9.49 0.318 1.94 -4.54 – 8.41 0.557
condition [Placebo] × Age
cat [Y] 		-0.88 -3.16 – 1.40 0.449
quartile × Age cat [Y] 0.24 -0.35 – 0.83 0.420
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.88 0.04 – 1.71 0.039
Random Effects
σ2 313.50 313.50 313.22
Ï„00 68.94 Subj_ID 68.95 Subj_ID 68.95 Subj_ID
ICC 0.18 0.18 0.18
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.715 / 0.767 0.716 / 0.767 0.716 / 0.767 
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) -26.51 -34.54 – -18.48 <0.001 -30.41 -41.50 – -19.31 <0.001 -27.20 -38.35 – -16.05 <0.001
condition [Placebo] -2.64 -4.38 – -0.91 0.003 -2.64 -4.38 – -0.91 0.003 -0.28 -2.68 – 2.11 0.816
quartile 42.39 41.94 – 42.84 <0.001 42.39 41.94 – 42.84 <0.001 41.51 40.89 – 42.13 <0.001
condition [Placebo] ×
quartile 		3.98 3.35 – 4.62 <0.001 3.98 3.35 – 4.62 <0.001 2.24 1.36 – 3.11 <0.001
Age cat [Y] 8.10 -7.80 – 23.99 0.318 1.43 -14.64 – 17.50 0.861
condition [Placebo] × Age
cat [Y] 		-4.90 -8.36 – -1.44 0.006
quartile × Age cat [Y] 1.84 0.94 – 2.73 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		3.63 2.37 – 4.90 <0.001
Random Effects
σ2 725.31 725.31 719.12
Ï„00 442.48 Subj_ID 442.54 Subj_ID 442.55 Subj_ID
ICC 0.38 0.38 0.38
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.680 / 0.801 0.681 / 0.802 0.683 / 0.804 
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) -16.10 -27.34 – -4.87 0.005 -21.70 -37.23 – -6.17 0.006 -16.56 -32.15 – -0.98 0.037
condition [Placebo] -1.74 -3.78 – 0.29 0.092 -1.74 -3.78 – 0.29 0.092 1.24 -1.56 – 4.04 0.385
quartile 47.02 46.49 – 47.54 <0.001 47.02 46.49 – 47.54 <0.001 45.56 44.84 – 46.29 <0.001
condition [Placebo] ×
quartile 		4.62 3.88 – 5.37 <0.001 4.62 3.88 – 5.37 <0.001 2.22 1.20 – 3.24 <0.001
Age cat [Y] 11.62 -10.68 – 33.92 0.307 0.96 -21.51 – 23.43 0.933
condition [Placebo] × Age
cat [Y] 		-6.20 -10.24 – -2.16 0.003
quartile × Age cat [Y] 3.02 1.98 – 4.06 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		4.99 3.52 – 6.47 <0.001
Random Effects
σ2 993.79 993.79 980.02
Ï„00 872.76 Subj_ID 871.29 Subj_ID 871.30 Subj_ID
ICC 0.47 0.47 0.47
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.622 / 0.799 0.625 / 0.800 0.628 / 0.803 
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 - Pearsons

Overall

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,variables_ext) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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))

data.aggr.withvars_diff %>% filter(threshold==0.25,measure=="IL-6_Delta(time2-baseline)") %>% ggscatterstats(
  ## arguments relevant for ggscatterstats
  data = .,
  x = meanStrength_diff,
  y = value_diff)

# IL-6 NS GEff 25%



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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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))

data.aggr.withvars_diff %>% filter(threshold==0.25,measure=="IL-6_Delta(time2-baseline)") %>% grouped_ggscatterstats(
  ## arguments relevant for ggscatterstats
  data = .,
  x = meanStrength_diff,
  y = value_diff,
  grouping.var = name)

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) %>% 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) %>% 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) %>% 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) %>% 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) %>% 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) 62.76 31.78 – 93.73 <0.001 43.85 1.35 – 86.36 0.043 44.90 2.36 – 87.44 0.039
condition [Placebo] 17.68 14.99 – 20.37 <0.001 17.68 14.99 – 20.37 <0.001 11.44 7.74 – 15.15 <0.001
quartile 41.04 40.35 – 41.74 <0.001 41.04 40.35 – 41.74 <0.001 42.53 41.57 – 43.49 <0.001
condition [Placebo] ×
quartile 		1.65 0.66 – 2.63 0.001 1.65 0.66 – 2.63 0.001 0.34 -1.02 – 1.69 0.625
Age cat [Y] 39.26 -21.94 – 100.46 0.209 37.09 -24.22 – 98.40 0.236
condition [Placebo] × Age
cat [Y] 		12.95 7.60 – 18.29 <0.001
quartile × Age cat [Y] -3.08 -4.47 – -1.70 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		2.72 0.76 – 4.67 0.006
Random Effects
σ2 1742.87 1742.87 1717.28
Ï„00 6716.48 Subj_ID 6569.67 Subj_ID 6569.70 Subj_ID
ICC 0.79 0.79 0.79
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.215 / 0.838 0.245 / 0.842 0.247 / 0.844 
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) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanStrength = mean(strength)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% 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) %>% 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) %>% 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) %>% 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) -386.48 -400.64 – -372.32 <0.001 -384.53 -401.82 – -367.23 <0.001 -378.93 -398.76 – -359.10 <0.001
condition [Placebo] -35.00 -49.87 – -20.13 <0.001 -35.00 -49.87 – -20.13 <0.001 -42.13 -62.77 – -21.48 <0.001
quartile 317.81 313.96 – 321.65 <0.001 317.81 313.96 – 321.65 <0.001 313.61 308.27 – 318.95 <0.001
condition [Placebo] ×
quartile 		17.69 12.25 – 23.13 <0.001 17.69 12.25 – 23.13 <0.001 24.45 16.89 – 32.00 <0.001
Age cat [Y] -4.05 -24.32 – 16.22 0.695 -15.68 -44.25 – 12.89 0.282
condition [Placebo] × Age
cat [Y] 		14.81 -14.95 – 44.56 0.329
quartile × Age cat [Y] 8.71 1.02 – 16.41 0.027
(condition [Placebo] ×
quartile) × Age cat [Y] 		-14.04 -24.92 – -3.15 0.011
Random Effects
σ2 159760.05 159760.05 159725.49
Ï„00 632.33 Subj_ID 655.80 Subj_ID 655.81 Subj_ID
ICC 0.00 0.00 0.00
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 66420 66420 66420
Marginal R2 / Conditional R2 0.454 / 0.456 0.454 / 0.456 0.454 / 0.457 
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) -674.49 -704.99 – -644.00 <0.001 -675.73 -711.65 – -639.81 <0.001 -669.34 -712.00 – -626.68 <0.001
condition [Placebo] -56.01 -91.32 – -20.71 0.002 -56.01 -91.32 – -20.71 0.002 -67.28 -116.30 – -18.26 0.007
quartile 493.56 484.43 – 502.69 <0.001 493.56 484.43 – 502.69 <0.001 486.99 474.32 – 499.67 <0.001
condition [Placebo] ×
quartile 		31.45 18.54 – 44.36 <0.001 31.45 18.54 – 44.36 <0.001 43.98 26.05 – 61.90 <0.001
Age cat [Y] 2.57 -36.07 – 41.21 0.896 -10.71 -72.19 – 50.77 0.733
condition [Placebo] × Age
cat [Y] 		23.40 -47.25 – 94.04 0.516
quartile × Age cat [Y] 13.64 -4.63 – 31.91 0.143
(condition [Placebo] ×
quartile) × Age cat [Y] 		-26.01 -51.85 – -0.18 0.048
Random Effects
σ2 300262.95 300262.95 300142.98
Ï„00 2154.77 Subj_ID 2253.84 Subj_ID 2253.97 Subj_ID
ICC 0.01 0.01 0.01
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.518 / 0.521 0.518 / 0.521 0.518 / 0.521 
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) -297.97 -310.40 – -285.54 <0.001 -293.46 -307.96 – -278.96 <0.001 -289.48 -306.69 – -272.26 <0.001
condition [Placebo] -20.52 -34.75 – -6.29 0.005 -20.52 -34.75 – -6.29 0.005 -25.28 -45.04 – -5.52 0.012
quartile 264.61 260.93 – 268.29 <0.001 264.61 260.93 – 268.29 <0.001 261.89 256.78 – 267.00 <0.001
condition [Placebo] ×
quartile 		9.72 4.52 – 14.93 <0.001 9.72 4.52 – 14.93 <0.001 13.87 6.64 – 21.10 <0.001
Age cat [Y] -9.38 -24.99 – 6.24 0.239 -17.64 -42.45 – 7.17 0.163
condition [Placebo] × Age
cat [Y] 		9.88 -18.59 – 38.36 0.496
quartile × Age cat [Y] 5.64 -1.73 – 13.00 0.133
(condition [Placebo] ×
quartile) × Age cat [Y] 		-8.62 -19.03 – 1.80 0.105
Random Effects
σ2 48784.03 48784.04 48775.86
Ï„00 374.69 Subj_ID 368.35 Subj_ID 368.36 Subj_ID
ICC 0.01 0.01 0.01
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.649 / 0.651 0.649 / 0.651 0.649 / 0.652 
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) -186.97 -195.70 – -178.24 <0.001 -184.40 -194.55 – -174.24 <0.001 -177.97 -190.09 – -165.85 <0.001
condition [Placebo] -28.46 -38.61 – -18.31 <0.001 -28.46 -38.61 – -18.31 <0.001 -33.83 -47.92 – -19.73 <0.001
quartile 195.25 192.63 – 197.88 <0.001 195.25 192.63 – 197.88 <0.001 191.95 188.31 – 195.60 <0.001
condition [Placebo] ×
quartile 		11.89 8.18 – 15.60 <0.001 11.89 8.18 – 15.60 <0.001 15.49 10.33 – 20.64 <0.001
Age cat [Y] -5.34 -16.12 – 5.43 0.331 -18.69 -36.16 – -1.22 0.036
condition [Placebo] × Age
cat [Y] 		11.14 -9.17 – 31.45 0.282
quartile × Age cat [Y] 6.85 1.60 – 12.10 0.011
(condition [Placebo] ×
quartile) × Age cat [Y] 		-7.47 -14.90 – -0.04 0.049
Random Effects
σ2 24826.21 24826.21 24818.66
Ï„00 173.07 Subj_ID 173.50 Subj_ID 173.51 Subj_ID
ICC 0.01 0.01 0.01
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.670 / 0.672 0.670 / 0.672 0.670 / 0.672 
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) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% 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) %>% 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) %>% 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) %>% 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) -2989.59 -3196.26 – -2782.92 <0.001 -2875.63 -3103.51 – -2647.75 <0.001 -3131.56 -3412.92 – -2850.20 <0.001
condition [Placebo] 106.39 -147.02 – 359.79 0.411 106.39 -147.02 – 359.79 0.411 -16.73 -368.37 – 334.92 0.926
quartile 1897.97 1832.44 – 1963.51 <0.001 1897.97 1832.44 – 1963.51 <0.001 1984.91 1893.96 – 2075.85 <0.001
condition [Placebo] ×
quartile 		-68.98 -161.66 – 23.70 0.145 -68.98 -161.66 – 23.70 0.145 11.64 -116.97 – 140.25 0.859
Age cat [Y] -236.69 -452.90 – -20.48 0.032 294.86 -110.63 – 700.35 0.154
condition [Placebo] × Age
cat [Y] 		255.70 -251.07 – 762.47 0.323
quartile × Age cat [Y] -180.55 -311.61 – -49.49 0.007
(condition [Placebo] ×
quartile) × Age cat [Y] 		-167.45 -352.80 – 17.90 0.077
Random Effects
σ2 15468924.61 15468924.63 15445362.69
Ï„00 74522.52 Subj_ID 63152.48 Subj_ID 63181.21 Subj_ID
ICC 0.00 0.00 0.00
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.218 / 0.222 0.219 / 0.222 0.220 / 0.224 
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) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanbwc = mean(bwc)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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 Pearsons

Overall

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

data.aggr %>% left_join(.,variables_ext) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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))

data.aggr.withvars %>%
  dplyr::group_by(Subj_ID,measure,threshold) %>%
  summarise(meangeff_diff = diff(meangeff), value_diff=diff(value)) %>% arrange(Subj_ID,measure,threshold)
## # A tibble: 1,620 × 5
## # Groups:   Subj_ID, measure [540]
##    Subj_ID measure                    threshold meangeff_diff value_diff
##      <dbl> <chr>                          <dbl>         <dbl>      <dbl>
##  1  122025 Avoid_Loss                      0.1        -0.0287      8.40 
##  2  122025 Avoid_Loss                      0.25       -0.0717      8.40 
##  3  122025 Avoid_Loss                      0.35       -0.0785      8.40 
##  4  122025 Gain                            0.1        -0.0287      8.3  
##  5  122025 Gain                            0.25       -0.0717      8.3  
##  6  122025 Gain                            0.35       -0.0785      8.3  
##  7  122025 IL-6_Delta(time2-baseline)      0.1        -0.0287     -3.65 
##  8  122025 IL-6_Delta(time2-baseline)      0.25       -0.0717     -3.65 
##  9  122025 IL-6_Delta(time2-baseline)      0.35       -0.0785     -3.65 
## 10  122025 IL6_Z                           0.1        -0.0287     -0.410
## # ℹ 1,610 more rows
data.aggr.withvars_diff %>% filter(threshold==0.25,measure=="IL-6_Delta(time2-baseline)") %>% ggscatterstats(
  ## arguments relevant for ggscatterstats
  data = .,
  x = meangeff_diff,
  y = value_diff)

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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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)

# IL-6 NS GEff 25%

data.aggr.withvars_diff %>% filter(threshold==0.25,measure=="IL-6_Delta(time2-baseline)") %>% grouped_ggscatterstats(
  ## arguments relevant for ggscatterstats
  data = .,
  x = meangeff_diff,
  y = value_diff,
  grouping.var = name)

# 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) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meangeff = mean(geff)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% 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) %>% 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) %>% 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) %>% 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.19 0.13 – 0.24 <0.001 0.16 0.09 – 0.22 <0.001 0.17 0.09 – 0.24 <0.001
condition [Placebo] 0.01 -0.03 – 0.04 0.741 0.01 -0.03 – 0.04 0.741 0.01 -0.04 – 0.06 0.643
quartile 0.15 0.14 – 0.16 <0.001 0.15 0.14 – 0.16 <0.001 0.15 0.14 – 0.16 <0.001
condition [Placebo] ×
quartile 		0.02 0.01 – 0.03 0.002 0.02 0.01 – 0.03 0.002 0.02 0.00 – 0.04 0.050
Age cat [Y] 0.07 -0.02 – 0.16 0.153 0.04 -0.06 – 0.15 0.391
condition [Placebo] × Age
cat [Y] 		-0.01 -0.08 – 0.06 0.744
quartile × Age cat [Y] 0.01 -0.01 – 0.03 0.381
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.01 -0.02 – 0.03 0.705
Random Effects
σ2 0.02 0.02 0.02
Ï„00 0.01 Subj_ID 0.01 Subj_ID 0.01 Subj_ID
ICC 0.45 0.44 0.44
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 1296 1296 1296
Marginal R2 / Conditional R2 0.512 / 0.730 0.524 / 0.732 0.524 / 0.732 
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.05 -0.00 – 0.10 0.055 0.02 -0.04 – 0.08 0.500 0.03 -0.04 – 0.11 0.355
condition [Placebo] 0.00 -0.06 – 0.06 0.986 0.00 -0.06 – 0.06 0.986 0.00 -0.08 – 0.08 0.933
quartile 0.18 0.16 – 0.19 <0.001 0.18 0.16 – 0.19 <0.001 0.17 0.15 – 0.19 <0.001
condition [Placebo] ×
quartile 		0.02 -0.00 – 0.04 0.069 0.02 -0.00 – 0.04 0.069 0.02 -0.01 – 0.05 0.227
Age cat [Y] 0.06 -0.00 – 0.13 0.066 0.04 -0.07 – 0.14 0.485
condition [Placebo] × Age
cat [Y] 		-0.01 -0.12 – 0.11 0.919
quartile × Age cat [Y] 0.01 -0.02 – 0.04 0.495
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.00 -0.04 – 0.05 0.887
Random Effects
σ2 0.02 0.02 0.02
Ï„00 0.01 Subj_ID 0.01 Subj_ID 0.01 Subj_ID
ICC 0.34 0.31 0.31
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 432 432 432
Marginal R2 / Conditional R2 0.656 / 0.771 0.668 / 0.772 0.668 / 0.772 
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.23 0.17 – 0.29 <0.001 0.20 0.12 – 0.27 <0.001 0.20 0.11 – 0.29 <0.001
condition [Placebo] 0.00 -0.05 – 0.06 0.887 0.00 -0.05 – 0.06 0.887 0.02 -0.06 – 0.09 0.688
quartile 0.15 0.13 – 0.16 <0.001 0.15 0.13 – 0.16 <0.001 0.14 0.12 – 0.16 <0.001
condition [Placebo] ×
quartile 		0.02 0.00 – 0.04 0.036 0.02 0.00 – 0.04 0.036 0.02 -0.01 – 0.05 0.225
Age cat [Y] 0.07 -0.03 – 0.16 0.177 0.06 -0.07 – 0.18 0.358
condition [Placebo] × Age
cat [Y] 		-0.02 -0.14 – 0.09 0.666
quartile × Age cat [Y] 0.00 -0.03 – 0.03 0.789
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.01 -0.03 – 0.05 0.674
Random Effects
σ2 0.01 0.01 0.02
Ï„00 0.02 Subj_ID 0.02 Subj_ID 0.02 Subj_ID
ICC 0.51 0.51 0.50
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 432 432 432
Marginal R2 / Conditional R2 0.510 / 0.762 0.522 / 0.764 0.522 / 0.763 
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.28 0.22 – 0.35 <0.001 0.25 0.16 – 0.33 <0.001 0.26 0.17 – 0.35 <0.001
condition [Placebo] 0.01 -0.04 – 0.07 0.623 0.01 -0.04 – 0.07 0.623 0.02 -0.06 – 0.09 0.679
quartile 0.13 0.12 – 0.15 <0.001 0.13 0.12 – 0.15 <0.001 0.13 0.11 – 0.15 <0.001
condition [Placebo] ×
quartile 		0.02 0.00 – 0.04 0.043 0.02 0.00 – 0.04 0.043 0.02 -0.01 – 0.05 0.180
Age cat [Y] 0.07 -0.04 – 0.18 0.219 0.04 -0.09 – 0.17 0.546
condition [Placebo] × Age
cat [Y] 		-0.00 -0.11 – 0.10 0.930
quartile × Age cat [Y] 0.01 -0.02 – 0.04 0.464
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.00 -0.04 – 0.04 0.871
Random Effects
σ2 0.01 0.01 0.01
Ï„00 0.02 Subj_ID 0.02 Subj_ID 0.02 Subj_ID
ICC 0.60 0.59 0.59
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 432 432 432
Marginal R2 / Conditional R2 0.447 / 0.777 0.460 / 0.780 0.461 / 0.780 
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) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meangeff = mean(geff)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% 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) %>% 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) %>% 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) %>% 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.16 -0.29 – -0.04 0.012 -0.23 -0.40 – -0.07 0.005 -0.21 -0.39 – -0.04 0.017
condition [Placebo] 0.01 -0.09 – 0.12 0.801 0.01 -0.09 – 0.12 0.801 -0.04 -0.18 – 0.10 0.585
quartile 0.15 0.13 – 0.18 <0.001 0.15 0.13 – 0.18 <0.001 0.15 0.11 – 0.19 <0.001
condition [Placebo] ×
quartile 		0.03 -0.00 – 0.07 0.073 0.03 -0.00 – 0.07 0.073 0.05 -0.01 – 0.10 0.089
Age cat [Y] 0.14 -0.07 – 0.35 0.180 0.10 -0.15 – 0.36 0.416
condition [Placebo] × Age
cat [Y] 		0.11 -0.10 – 0.32 0.295
quartile × Age cat [Y] 0.01 -0.05 – 0.06 0.852
(condition [Placebo] ×
quartile) × Age cat [Y] 		-0.02 -0.10 – 0.05 0.556
Random Effects
σ2 0.05 0.05 0.05
Ï„00 0.08 Subj_ID 0.07 Subj_ID 0.07 Subj_ID
ICC 0.61 0.60 0.60
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 432 432 432
Marginal R2 / Conditional R2 0.234 / 0.702 0.262 / 0.707 0.262 / 0.706 
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) %>% 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) %>% 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) %>% 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) %>% 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) %>% 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.48 0.42 – 0.55 <0.001 0.44 0.35 – 0.53 <0.001 0.44 0.35 – 0.53 <0.001
condition [Placebo] 0.05 0.04 – 0.06 <0.001 0.05 0.04 – 0.06 <0.001 0.06 0.05 – 0.07 <0.001
quartile 0.09 0.09 – 0.09 <0.001 0.09 0.09 – 0.09 <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.756
Age cat [Y] 0.08 -0.05 – 0.21 0.233 0.09 -0.04 – 0.22 0.186
condition [Placebo] × Age
cat [Y] 		-0.02 -0.03 – -0.00 0.010
quartile × Age cat [Y] -0.01 -0.01 – -0.00 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.02 0.01 – 0.02 <0.001
Random Effects
σ2 0.04 0.04 0.04
Ï„00 0.03 Subj_ID 0.03 Subj_ID 0.03 Subj_ID
ICC 0.45 0.44 0.44
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 66420 66420 66420
Marginal R2 / Conditional R2 0.152 / 0.532 0.169 / 0.538 0.170 / 0.539 
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.38 0.32 – 0.44 <0.001 0.34 0.26 – 0.42 <0.001 0.32 0.24 – 0.39 <0.001
condition [Placebo] 0.05 0.03 – 0.06 <0.001 0.05 0.03 – 0.06 <0.001 0.08 0.05 – 0.10 <0.001
quartile 0.16 0.15 – 0.16 <0.001 0.16 0.15 – 0.16 <0.001 0.17 0.16 – 0.17 <0.001
condition [Placebo] ×
quartile 		0.01 0.00 – 0.01 0.003 0.01 0.00 – 0.01 0.003 -0.00 -0.01 – 0.00 0.381
Age cat [Y] 0.09 -0.02 – 0.20 0.126 0.13 0.02 – 0.25 0.021
condition [Placebo] × Age
cat [Y] 		-0.06 -0.09 – -0.03 <0.001
quartile × Age cat [Y] -0.02 -0.03 – -0.01 <0.001
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.02 0.01 – 0.04 <0.001
Random Effects
σ2 0.06 0.06 0.05
Ï„00 0.02 Subj_ID 0.02 Subj_ID 0.02 Subj_ID
ICC 0.29 0.28 0.28
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.299 / 0.505 0.314 / 0.508 0.315 / 0.508 
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.54 0.47 – 0.61 <0.001 0.51 0.41 – 0.60 <0.001 0.51 0.41 – 0.60 <0.001
condition [Placebo] 0.05 0.04 – 0.06 <0.001 0.05 0.04 – 0.06 <0.001 0.05 0.04 – 0.07 <0.001
quartile 0.06 0.06 – 0.07 <0.001 0.06 0.06 – 0.07 <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.970
Age cat [Y] 0.07 -0.06 – 0.21 0.288 0.07 -0.07 – 0.20 0.344
condition [Placebo] × Age
cat [Y] 		-0.01 -0.03 – 0.01 0.352
quartile × Age cat [Y] -0.00 -0.01 – 0.00 0.441
(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.03 Subj_ID 0.03 Subj_ID 0.03 Subj_ID
ICC 0.59 0.59 0.59
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.107 / 0.637 0.127 / 0.644 0.128 / 0.645 
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.53 0.46 – 0.60 <0.001 0.49 0.39 – 0.59 <0.001 0.50 0.40 – 0.59 <0.001
condition [Placebo] 0.05 0.04 – 0.06 <0.001 0.05 0.04 – 0.06 <0.001 0.04 0.03 – 0.05 <0.001
quartile 0.05 0.05 – 0.05 <0.001 0.05 0.05 – 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.385
Age cat [Y] 0.08 -0.07 – 0.22 0.295 0.06 -0.08 – 0.20 0.385
condition [Placebo] × Age
cat [Y] 		0.01 -0.01 – 0.03 0.237
quartile × Age cat [Y] -0.00 -0.01 – 0.00 0.564
(condition [Placebo] ×
quartile) × Age cat [Y] 		0.01 0.00 – 0.01 0.003
Random Effects
σ2 0.02 0.02 0.02
Ï„00 0.04 Subj_ID 0.03 Subj_ID 0.03 Subj_ID
ICC 0.69 0.69 0.69
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.081 / 0.718 0.104 / 0.725 0.105 / 0.726 
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 - Pearsons

Overall

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>% left_join(.,variables_ext) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% 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) %>% 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) %>% 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) %>% 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) %>% 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) %>% group_by(threshold,Subj_ID,condition) %>% summarise(meanleff = mean(leff)) -> data.aggr

data.aggr %>% left_join(.,variables_ext) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanleff = mean(leff)) -> data.aggr

colnames(data.aggr)

data.aggr %>% left_join(.,variables_ext) %>% pivot_longer(cols = `IL-6_Delta(time2-baseline)`:TNFd_Z, 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)