## Warning: package 'reactable' was built under R version 4.2.3

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/", 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  54
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/", 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/", 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_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/sub-0122025_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0122026_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0122027_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0222028_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0222029_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0322030_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0322031_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0322032_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0421002_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0422033_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0422034_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0422035_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0521003_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0522036_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0721005_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0821006_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0821007_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0921009_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0921010_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0921012_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0921014_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0921015_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0921016_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-0921017_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-1021019_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/sub-1121020_ses-I_task-rest_space-MNI152NLin2009cAsym_atlas-Combined_measure-pearsoncorrelation_conmat.tsv 
##                                                                                                                                  0 
## conmats_filtered_tedana/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.5494e+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.2814e+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.6141e+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/", pattern = "ses-P",full.names = TRUE)
sub_numbers <- gsub(".*sub-(\\d+)_ses.*", "\\1", files)
sub_numbers_df <- data.frame(sub = 1:length(sub_numbers), Subj_ID = as.numeric(sub_numbers))

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

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

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

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

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

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

List of most Interferon affected connections in graph

Where is the most condition related change happening?

library(igraph)

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

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

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


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

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

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

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

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

df_list <- list()

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

data.funccon <- bind_rows(df_list) %>% left_join(.,variables_ext)

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

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=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=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 + (1 | Subj_ID), data.aggr )

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

tab_model(fit1,fit2,fit3)
  strength strength strength
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -22.05 -29.70 – -14.39 <0.001 -16.25 -34.72 – 2.22 0.085 -19.52 -38.26 – -0.78 0.041
condition [Placebo] -2.88 -5.00 – -0.77 0.008 -2.88 -5.00 – -0.77 0.008 -5.00 -10.13 – 0.12 0.056
quartile 37.62 37.07 – 38.17 <0.001 37.62 37.07 – 38.17 <0.001 38.45 37.12 – 39.78 <0.001
condition [Placebo] ×
quartile 		4.27 3.50 – 5.04 <0.001 4.27 3.50 – 5.04 <0.001 6.08 4.20 – 7.95 <0.001
Age -0.12 -0.49 – 0.24 0.499 -0.05 -0.42 – 0.31 0.772
condition [Placebo] × Age 0.05 -0.05 – 0.15 0.373
quartile × Age -0.02 -0.04 – 0.01 0.177
(condition [Placebo] ×
quartile) × Age 		-0.04 -0.08 – -0.00 0.038
Random Effects
σ2 2156.86 2156.86 2155.73
Ï„00 396.03 Subj_ID 404.49 Subj_ID 404.49 Subj_ID
ICC 0.16 0.16 0.16
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 44280 44280 44280
Marginal R2 / Conditional R2 0.439 / 0.526 0.439 / 0.527 0.439 / 0.528 
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 + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit3 <- lmer(strength ~ condition*quartile*Age + (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) -28.08 -31.48 – -24.68 <0.001 -25.61 -33.74 – -17.47 <0.001 -25.09 -33.42 – -16.76 <0.001
condition [Placebo] -3.32 -4.51 – -2.13 <0.001 -3.32 -4.51 – -2.13 <0.001 -3.13 -6.01 – -0.25 0.033
quartile 27.15 26.84 – 27.46 <0.001 27.15 26.84 – 27.46 <0.001 26.81 26.06 – 27.55 <0.001
condition [Placebo] ×
quartile 		2.88 2.44 – 3.31 <0.001 2.88 2.44 – 3.31 <0.001 3.07 2.02 – 4.13 <0.001
Age -0.05 -0.21 – 0.11 0.511 -0.06 -0.23 – 0.10 0.440
condition [Placebo] × Age -0.00 -0.06 – 0.05 0.888
quartile × Age 0.01 -0.01 – 0.02 0.323
(condition [Placebo] ×
quartile) × Age 		-0.00 -0.02 – 0.02 0.688
Random Effects
σ2 340.47 340.47 340.48
Ï„00 76.31 Subj_ID 78.02 Subj_ID 78.02 Subj_ID
ICC 0.18 0.19 0.19
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.712 / 0.764 0.711 / 0.765 0.711 / 0.765 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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

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

fit3 <- lmer(strength ~ condition*quartile*Age + (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.01 -27.85 – -4.17 0.008 -6.90 -35.71 – 21.91 0.639 -13.95 -42.94 – 15.03 0.345
condition [Placebo] -2.45 -4.56 – -0.33 0.024 -2.45 -4.56 – -0.33 0.024 -6.88 -12.00 – -1.76 0.008
quartile 48.09 47.54 – 48.64 <0.001 48.09 47.54 – 48.64 <0.001 50.09 48.77 – 51.42 <0.001
condition [Placebo] ×
quartile 		5.67 4.89 – 6.44 <0.001 5.67 4.89 – 6.44 <0.001 9.09 7.21 – 10.96 <0.001
Age -0.20 -0.76 – 0.37 0.496 -0.04 -0.61 – 0.52 0.879
condition [Placebo] × Age 0.10 -0.01 – 0.20 0.063
quartile × Age -0.04 -0.07 – -0.02 0.001
(condition [Placebo] ×
quartile) × Age 		-0.07 -0.11 – -0.04 <0.001
Random Effects
σ2 1080.26 1080.26 1075.25
Ï„00 969.00 Subj_ID 989.42 Subj_ID 989.43 Subj_ID
ICC 0.47 0.48 0.48
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.616 / 0.798 0.615 / 0.799 0.616 / 0.800 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

Behavior - paired difference correlation analysis

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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% 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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

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 + (1 | Subj_ID), data.aggr )

fit3 <- lmer(strength ~ condition*quartile*Age + (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.70 29.96 – 95.44 <0.001 96.85 17.39 – 176.31 0.017 97.86 18.29 – 177.43 0.016
condition [Placebo] 20.76 17.96 – 23.57 <0.001 20.76 17.96 – 23.57 <0.001 34.09 27.32 – 40.86 <0.001
quartile 42.36 41.63 – 43.09 <0.001 42.36 41.63 – 43.09 <0.001 38.72 36.97 – 40.47 <0.001
condition [Placebo] ×
quartile 		2.19 1.16 – 3.21 <0.001 2.19 1.16 – 3.21 <0.001 3.31 0.83 – 5.78 0.009
Age -0.73 -2.29 – 0.82 0.355 -0.76 -2.31 – 0.80 0.342
condition [Placebo] × Age -0.29 -0.42 – -0.15 <0.001
quartile × Age 0.08 0.04 – 0.11 <0.001
(condition [Placebo] ×
quartile) × Age 		-0.02 -0.07 – 0.02 0.331
Random Effects
σ2 1893.62 1893.62 1878.04
Ï„00 7505.03 Subj_ID 7546.96 Subj_ID 7546.96 Subj_ID
ICC 0.80 0.80 0.80
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.212 / 0.841 0.227 / 0.845 0.228 / 0.846 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

Grouped by network

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "strength") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% 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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

Betweenness centrality - using BCT Toolbox’s betweenness_wei.m

Thresholded

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% 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=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=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 + (1 | Subj_ID), data.aggr )

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

tab_model(fit1,fit2,fit3)
  bwc bwc bwc
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -430.91 -448.95 – -412.87 <0.001 -436.56 -467.82 – -405.29 <0.001 -448.86 -492.85 – -404.87 <0.001
condition [Placebo] -42.09 -62.80 – -21.38 <0.001 -42.09 -62.80 – -21.38 <0.001 -8.50 -58.70 – 41.69 0.740
quartile 346.16 340.81 – 351.52 <0.001 346.16 340.81 – 351.52 <0.001 357.35 344.37 – 370.33 <0.001
condition [Placebo] ×
quartile 		21.54 13.97 – 29.12 <0.001 21.54 13.97 – 29.12 <0.001 -4.44 -22.79 – 13.92 0.636
Age 0.12 -0.43 – 0.67 0.664 0.39 -0.48 – 1.25 0.380
condition [Placebo] × Age -0.72 -1.70 – 0.26 0.150
quartile × Age -0.24 -0.49 – 0.01 0.064
(condition [Placebo] ×
quartile) × Age 		0.56 0.20 – 0.92 0.002
Random Effects
σ2 206703.10 206703.10 206622.92
Ï„00 779.09 Subj_ID 808.23 Subj_ID 808.27 Subj_ID
ICC 0.00 0.00 0.00
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 44280 44280 44280
Marginal R2 / Conditional R2 0.435 / 0.437 0.434 / 0.437 0.435 / 0.437 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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

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

fit3 <- lmer(bwc ~ condition*quartile*Age + (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) -670.52 -700.35 – -640.68 <0.001 -674.27 -725.32 – -623.22 <0.001 -690.11 -762.94 – -617.28 <0.001
condition [Placebo] -60.77 -95.54 – -26.00 0.001 -60.77 -95.54 – -26.00 0.001 -10.48 -94.72 – 73.75 0.807
quartile 494.35 485.36 – 503.34 <0.001 494.35 485.36 – 503.34 <0.001 511.22 489.43 – 533.00 <0.001
condition [Placebo] ×
quartile 		33.41 20.69 – 46.12 <0.001 33.41 20.69 – 46.12 <0.001 -7.78 -38.59 – 23.03 0.620
Age 0.08 -0.81 – 0.97 0.858 0.42 -1.00 – 1.85 0.563
condition [Placebo] × Age -1.08 -2.73 – 0.57 0.199
quartile × Age -0.36 -0.79 – 0.06 0.096
(condition [Placebo] ×
quartile) × Age 		0.89 0.28 – 1.49 0.004
Random Effects
σ2 291189.41 291189.41 290974.94
Ï„00 2008.99 Subj_ID 2100.43 Subj_ID 2100.69 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.527 / 0.530 0.527 / 0.531 0.528 / 0.531 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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

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

fit3 <- lmer(bwc ~ condition*quartile*Age + (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) -191.30 -200.21 – -182.39 <0.001 -198.85 -213.98 – -183.71 <0.001 -207.61 -229.14 – -186.08 <0.001
condition [Placebo] -23.41 -33.66 – -13.15 <0.001 -23.41 -33.66 – -13.15 <0.001 -6.53 -31.37 – 18.31 0.607
quartile 197.97 195.32 – 200.62 <0.001 197.97 195.32 – 200.62 <0.001 203.48 197.06 – 209.91 <0.001
condition [Placebo] ×
quartile 		9.68 5.93 – 13.42 <0.001 9.68 5.93 – 13.42 <0.001 -1.09 -10.18 – 7.99 0.814
Age 0.16 -0.10 – 0.43 0.227 0.35 -0.07 – 0.77 0.103
condition [Placebo] × Age -0.36 -0.85 – 0.12 0.144
quartile × Age -0.12 -0.24 – 0.01 0.065
(condition [Placebo] ×
quartile) × Age 		0.23 0.05 – 0.41 0.011
Random Effects
σ2 25311.62 25311.62 25302.49
Ï„00 189.26 Subj_ID 185.45 Subj_ID 185.46 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.669 / 0.671 0.669 / 0.671 0.669 / 0.671 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

Grouped by network

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% 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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

Unthresholded

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,variables_ext) %>% 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 + (1 | Subj_ID), data.aggr )

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

tab_model(fit1,fit2,fit3)
  bwc bwc bwc
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -2903.28 -3106.19 – -2700.37 <0.001 -3099.69 -3408.03 – -2791.36 <0.001 -2578.36 -3065.60 – -2091.11 <0.001
condition [Placebo] -1.52 -254.08 – 251.03 0.991 -1.52 -254.08 – 251.03 0.991 56.66 -555.25 – 668.57 0.856
quartile 1853.48 1788.17 – 1918.80 <0.001 1853.48 1788.17 – 1918.80 <0.001 1655.79 1497.54 – 1814.04 <0.001
condition [Placebo] ×
quartile 		-15.45 -107.82 – 76.93 0.743 -15.45 -107.82 – 76.93 0.743 -61.26 -285.07 – 162.54 0.592
Age 4.22 -0.80 – 9.25 0.100 -6.98 -16.53 – 2.56 0.151
condition [Placebo] × Age -1.25 -13.23 – 10.73 0.838
quartile × Age 4.25 1.15 – 7.35 0.007
(condition [Placebo] ×
quartile) × Age 		0.98 -3.40 – 5.37 0.660
Random Effects
σ2 15365542.06 15365542.05 15354832.13
Ï„00 65224.49 Subj_ID 60035.09 Subj_ID 60048.16 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.216 / 0.220 0.217 / 0.220 0.217 / 0.220 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

Grouped by network

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "bwc") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% 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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

Global efficiency - using BCT Toolbox’s efficiency_wei.m

Thresholded

Behavior - paired difference correlation analysis

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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% 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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% 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=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=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 + (1 | Subj_ID), data.aggr )

fit3 <- lmer(geff ~ condition*quartile*Age + (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.17 0.11 – 0.22 <0.001 0.22 0.11 – 0.33 <0.001 0.22 0.08 – 0.35 0.001
condition [Placebo] 0.01 -0.04 – 0.05 0.764 0.01 -0.04 – 0.05 0.764 -0.01 -0.12 – 0.10 0.886
quartile 0.16 0.14 – 0.17 <0.001 0.16 0.14 – 0.17 <0.001 0.16 0.13 – 0.19 <0.001
condition [Placebo] ×
quartile 		0.02 0.00 – 0.04 0.020 0.02 0.00 – 0.04 0.020 0.02 -0.02 – 0.07 0.227
Age -0.00 -0.00 – 0.00 0.264 -0.00 -0.00 – 0.00 0.426
condition [Placebo] × Age 0.00 -0.00 – 0.00 0.770
quartile × Age -0.00 -0.00 – 0.00 0.796
(condition [Placebo] ×
quartile) × Age 		-0.00 -0.00 – 0.00 0.788
Random Effects
σ2 0.02 0.02 0.02
Ï„00 0.01 Subj_ID 0.01 Subj_ID 0.01 Subj_ID
ICC 0.41 0.41 0.40
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 864 864 864
Marginal R2 / Conditional R2 0.517 / 0.714 0.522 / 0.716 0.522 / 0.715 
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 + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit3 <- lmer(geff ~ condition*quartile*Age + (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.11 0.02 – 0.20 0.021 0.10 -0.02 – 0.23 0.111
condition [Placebo] 0.00 -0.06 – 0.06 0.986 0.00 -0.06 – 0.06 0.986 -0.01 -0.16 – 0.13 0.836
quartile 0.18 0.16 – 0.19 <0.001 0.18 0.16 – 0.19 <0.001 0.18 0.14 – 0.22 <0.001
condition [Placebo] ×
quartile 		0.02 -0.00 – 0.04 0.069 0.02 -0.00 – 0.04 0.069 0.02 -0.03 – 0.08 0.347
Age -0.00 -0.00 – 0.00 0.142 -0.00 -0.00 – 0.00 0.384
condition [Placebo] × Age 0.00 -0.00 – 0.00 0.813
quartile × Age -0.00 -0.00 – 0.00 0.859
(condition [Placebo] ×
quartile) × Age 		-0.00 -0.00 – 0.00 0.833
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.32 0.32
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 432 432 432
Marginal R2 / Conditional R2 0.656 / 0.771 0.663 / 0.772 0.662 / 0.771 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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

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

fit3 <- lmer(geff ~ condition*quartile*Age + (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.34 0.20 – 0.47 <0.001 0.33 0.17 – 0.49 <0.001
condition [Placebo] 0.01 -0.04 – 0.07 0.623 0.01 -0.04 – 0.07 0.623 -0.00 -0.13 – 0.13 0.984
quartile 0.13 0.12 – 0.15 <0.001 0.13 0.12 – 0.15 <0.001 0.14 0.10 – 0.17 <0.001
condition [Placebo] ×
quartile 		0.02 0.00 – 0.04 0.043 0.02 0.00 – 0.04 0.043 0.03 -0.02 – 0.07 0.298
Age -0.00 -0.00 – 0.00 0.376 -0.00 -0.00 – 0.00 0.536
condition [Placebo] × Age 0.00 -0.00 – 0.00 0.807
quartile × Age -0.00 -0.00 – 0.00 0.803
(condition [Placebo] ×
quartile) × Age 		-0.00 -0.00 – 0.00 0.819
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.60 0.60
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 432 432 432
Marginal R2 / Conditional R2 0.447 / 0.777 0.451 / 0.780 0.451 / 0.779 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% 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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

Grouped by scale, condition, subject

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "geff") %>% left_join(.,variables_ext) %>% 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 + (1 | Subj_ID), data.aggr )

fit3 <- lmer(geff ~ condition*quartile*Age + (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.04 -0.31 – 0.23 0.757 -0.06 -0.37 – 0.25 0.718
condition [Placebo] 0.01 -0.09 – 0.12 0.801 0.01 -0.09 – 0.12 0.801 0.14 -0.11 – 0.39 0.258
quartile 0.15 0.13 – 0.18 <0.001 0.15 0.13 – 0.18 <0.001 0.15 0.08 – 0.21 <0.001
condition [Placebo] ×
quartile 		0.03 -0.00 – 0.07 0.073 0.03 -0.00 – 0.07 0.073 0.00 -0.09 – 0.10 0.922
Age -0.00 -0.01 – 0.00 0.310 -0.00 -0.01 – 0.00 0.454
condition [Placebo] × Age -0.00 -0.01 – 0.00 0.260
quartile × Age 0.00 -0.00 – 0.00 0.860
(condition [Placebo] ×
quartile) × Age 		0.00 -0.00 – 0.00 0.480
Random Effects
σ2 0.05 0.05 0.05
Ï„00 0.08 Subj_ID 0.08 Subj_ID 0.08 Subj_ID
ICC 0.61 0.61 0.61
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 432 432 432
Marginal R2 / Conditional R2 0.234 / 0.702 0.249 / 0.707 0.250 / 0.707 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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=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=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 + (1 | Subj_ID), data.aggr )

fit3 <- lmer(leff ~ condition*quartile*Age + (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.47 0.40 – 0.53 <0.001 0.54 0.37 – 0.70 <0.001 0.57 0.41 – 0.74 <0.001
condition [Placebo] 0.06 0.05 – 0.07 <0.001 0.06 0.05 – 0.07 <0.001 0.02 0.00 – 0.05 0.041
quartile 0.11 0.10 – 0.11 <0.001 0.11 0.10 – 0.11 <0.001 0.09 0.08 – 0.09 <0.001
condition [Placebo] ×
quartile 		0.01 0.00 – 0.01 <0.001 0.01 0.00 – 0.01 <0.001 0.02 0.02 – 0.03 <0.001
Age -0.00 -0.00 – 0.00 0.356 -0.00 -0.01 – 0.00 0.154
condition [Placebo] × Age 0.00 0.00 – 0.00 0.002
quartile × Age 0.00 0.00 – 0.00 <0.001
(condition [Placebo] ×
quartile) × Age 		-0.00 -0.00 – -0.00 <0.001
Random Effects
σ2 0.05 0.05 0.05
Ï„00 0.03 Subj_ID 0.03 Subj_ID 0.03 Subj_ID
ICC 0.41 0.41 0.41
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 44280 44280 44280
Marginal R2 / Conditional R2 0.176 / 0.513 0.184 / 0.519 0.184 / 0.519 
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 + (1 | Subj_ID), data.aggr %>% filter(threshold==0.1))

fit3 <- lmer(leff ~ condition*quartile*Age + (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.39 0.33 – 0.45 <0.001 0.47 0.32 – 0.61 <0.001 0.54 0.40 – 0.69 <0.001
condition [Placebo] 0.06 0.05 – 0.08 <0.001 0.06 0.05 – 0.08 <0.001 0.00 -0.04 – 0.04 0.996
quartile 0.16 0.16 – 0.16 <0.001 0.16 0.16 – 0.16 <0.001 0.13 0.12 – 0.14 <0.001
condition [Placebo] ×
quartile 		0.01 0.00 – 0.01 0.004 0.01 0.00 – 0.01 0.004 0.03 0.02 – 0.04 <0.001
Age -0.00 -0.00 – 0.00 0.233 -0.00 -0.01 – -0.00 0.022
condition [Placebo] × Age 0.00 0.00 – 0.00 <0.001
quartile × Age 0.00 0.00 – 0.00 <0.001
(condition [Placebo] ×
quartile) × Age 		-0.00 -0.00 – -0.00 0.001
Random Effects
σ2 0.06 0.06 0.06
Ï„00 0.03 Subj_ID 0.02 Subj_ID 0.02 Subj_ID
ICC 0.30 0.30 0.30
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.298 / 0.510 0.306 / 0.513 0.307 / 0.514 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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

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

fit3 <- lmer(leff ~ condition*quartile*Age + (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.54 0.47 – 0.62 <0.001 0.60 0.42 – 0.79 <0.001 0.60 0.42 – 0.79 <0.001
condition [Placebo] 0.06 0.05 – 0.06 <0.001 0.06 0.05 – 0.06 <0.001 0.05 0.03 – 0.07 <0.001
quartile 0.05 0.05 – 0.05 <0.001 0.05 0.05 – 0.05 <0.001 0.05 0.04 – 0.05 <0.001
condition [Placebo] ×
quartile 		0.01 0.00 – 0.01 <0.001 0.01 0.00 – 0.01 <0.001 0.02 0.01 – 0.02 <0.001
Age -0.00 -0.00 – 0.00 0.480 -0.00 -0.00 – 0.00 0.472
condition [Placebo] × Age 0.00 -0.00 – 0.00 0.516
quartile × Age 0.00 -0.00 – 0.00 0.102
(condition [Placebo] ×
quartile) × Age 		-0.00 -0.00 – -0.00 0.006
Random Effects
σ2 0.02 0.02 0.02
Ï„00 0.04 Subj_ID 0.04 Subj_ID 0.04 Subj_ID
ICC 0.71 0.71 0.71
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.083 / 0.731 0.093 / 0.737 0.093 / 0.738 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

Behavior - paired difference correlation analysis

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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>%   mutate(ROI = gsub("Var", "ROI_", ROI))%>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% 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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

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 + (1 | Subj_ID), data.aggr )

fit3 <- lmer(leff ~ condition*quartile*Age + (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.28 0.21 – 0.35 <0.001 0.35 0.18 – 0.53 <0.001 0.35 0.17 – 0.52 <0.001
condition [Placebo] 0.05 0.05 – 0.06 <0.001 0.05 0.05 – 0.06 <0.001 0.08 0.07 – 0.10 <0.001
quartile 0.07 0.06 – 0.07 <0.001 0.07 0.06 – 0.07 <0.001 0.06 0.06 – 0.07 <0.001
condition [Placebo] ×
quartile 		0.00 0.00 – 0.01 <0.001 0.00 0.00 – 0.01 <0.001 0.00 -0.00 – 0.01 0.107
Age -0.00 -0.00 – 0.00 0.358 -0.00 -0.00 – 0.00 0.399
condition [Placebo] × Age -0.00 -0.00 – -0.00 <0.001
quartile × Age 0.00 0.00 – 0.00 0.014
(condition [Placebo] ×
quartile) × Age 		0.00 -0.00 – 0.00 0.719
Random Effects
σ2 0.01 0.01 0.01
Ï„00 0.04 Subj_ID 0.04 Subj_ID 0.04 Subj_ID
ICC 0.85 0.85 0.85
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 22140 22140 22140
Marginal R2 / Conditional R2 0.139 / 0.872 0.158 / 0.875 0.159 / 0.876 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

Grouped by network

data %>% pivot_longer(
    cols = starts_with("Var"),
    names_to = "ROI",
    values_to = "leff") %>%   mutate(ROI = gsub("Var", "ROI_", ROI)) %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% 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)`:Avoid_Loss, names_to = "measure", values_to = "value") -> data.aggr.withvars 

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


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

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

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

datatable(cor_results)

Node degree - Graphvar

Grouped by scale, condition, subject

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

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

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

grouped_ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanDegree,
  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 = "degree") %>% left_join(.,variables_ext) %>% group_by(threshold,Subj_ID,condition, Age_cat) %>% summarise(meanDegree = mean(degree)) -> data.aggr

anova_summary <- nice(aov_ez("Subj_ID", "meanDegree", data.aggr, within =c("condition","threshold"), 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 = "degree") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% group_by(threshold,Subj_ID,condition,name) %>% summarise(meanDegree = mean(degree)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanDegree", 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", "meanDegree", 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", "meanDegree", data.aggr %>% filter(threshold=="1"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanDegree", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X1"))
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", "meanDegree", data.aggr %>% filter(threshold=="3.5"),within =c("condition","name"))))
main_anova_summary <- aov_ez("Subj_ID", "meanDegree", data.aggr,within =c("condition","name","threshold"))
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition | name, at=list(threshold="X3.5"))
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 = "degree") %>% left_join(.,network.combined) %>% left_join(.,variables_ext) %>% group_by(threshold,Subj_ID,condition,name,Age_cat) %>% summarise(meanDegree = mean(degree)) -> data.aggr

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanDegree", 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", "meanDegree", 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", "meanDegree", 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", "meanDegree", data.aggr %>% filter(threshold==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=3.5
main_anova_summary <- aov_ez("Subj_ID", "meanDegree", data.aggr %>% filter(threshold==3.5),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 = "degree") %>% left_join(.,variables_ext) %>% group_by(threshold,Subj_ID,condition) %>% mutate(quartile = ntile(degree, 4)) -> data.aggr

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

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

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

tab_model(fit1,fit2,fit3)
  degree degree degree
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -14.27 -16.10 – -12.43 <0.001 -13.64 -15.82 – -11.46 <0.001 -13.19 -17.65 – -8.74 <0.001
condition [Placebo] -4.50 -7.10 – -1.90 0.001 -4.50 -7.10 – -1.90 0.001 -6.43 -12.73 – -0.13 0.045
quartile 42.85 42.18 – 43.52 <0.001 42.85 42.18 – 43.52 <0.001 42.79 41.16 – 44.42 <0.001
condition [Placebo] ×
quartile 		1.30 0.35 – 2.25 0.007 1.30 0.35 – 2.25 0.007 1.84 -0.46 – 4.15 0.117
Age -0.01 -0.04 – 0.01 0.295 -0.02 -0.11 – 0.06 0.605
condition [Placebo] × Age 0.04 -0.08 – 0.16 0.511
quartile × Age 0.00 -0.03 – 0.03 0.934
(condition [Placebo] ×
quartile) × Age 		-0.01 -0.06 – 0.03 0.614
Random Effects
σ2 3257.49 3257.48 3257.65
Ï„00 0.00 Subj_ID 0.00 Subj_ID 0.00 Subj_ID
N 27 Subj_ID 27 Subj_ID 27 Subj_ID
Observations 44280 44280 44280
Marginal R2 / Conditional R2 0.421 / NA 0.421 / NA 0.421 / NA 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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

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

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

tab_model(fit1,fit2,fit3)
  degree degree degree
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -26.12 -26.75 – -25.49 <0.001 -25.86 -26.80 – -24.92 <0.001 -25.67 -27.21 – -24.13 <0.001
condition [Placebo] -1.88 -2.70 – -1.07 <0.001 -1.88 -2.70 – -1.07 <0.001 -2.63 -4.60 – -0.65 0.009
quartile 26.98 26.77 – 27.19 <0.001 26.98 26.77 – 27.19 <0.001 26.96 26.45 – 27.47 <0.001
condition [Placebo] ×
quartile 		0.52 0.22 – 0.82 0.001 0.52 0.22 – 0.82 0.001 0.71 -0.02 – 1.43 0.056
Age -0.01 -0.02 – 0.01 0.462 -0.01 -0.04 – 0.02 0.531
condition [Placebo] × Age 0.02 -0.02 – 0.05 0.420
quartile × Age 0.00 -0.01 – 0.01 0.929
(condition [Placebo] ×
quartile) × Age 		-0.00 -0.02 – 0.01 0.574
Random Effects
σ2 160.18 160.18 160.19
Ï„00 0.49 Subj_ID 0.50 Subj_ID 0.50 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.852 / 0.853 0.852 / 0.853 0.852 / 0.853 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

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

fit2 <- lmer(degree ~ condition*quartile + Age + (1 | Subj_ID), data.aggr %>% filter(threshold==3.5))

fit3 <- lmer(degree ~ condition*quartile*Age + (1 | Subj_ID), data.aggr %>% filter(threshold==3.5))

tab_model(fit1,fit2,fit3)
  degree degree degree
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -2.41 -3.64 – -1.18 <0.001 -1.42 -3.57 – 0.73 0.196 -0.72 -3.70 – 2.27 0.638
condition [Placebo] -7.12 -8.51 – -5.74 <0.001 -7.12 -8.51 – -5.74 <0.001 -10.24 -13.59 – -6.89 <0.001
quartile 58.72 58.36 – 59.08 <0.001 58.72 58.36 – 59.08 <0.001 58.62 57.75 – 59.49 <0.001
condition [Placebo] ×
quartile 		2.09 1.58 – 2.59 <0.001 2.09 1.58 – 2.59 <0.001 2.98 1.76 – 4.21 <0.001
Age -0.02 -0.06 – 0.02 0.270 -0.04 -0.09 – 0.02 0.221
condition [Placebo] × Age 0.07 0.00 – 0.13 0.046
quartile × Age 0.00 -0.01 – 0.02 0.797
(condition [Placebo] ×
quartile) × Age 		-0.02 -0.04 – 0.00 0.117
Random Effects
σ2 460.77 460.77 460.71
Ï„00 3.96 Subj_ID 3.92 Subj_ID 3.92 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.906 / 0.907 0.906 / 0.907 0.906 / 0.907 
plot_model(fit1, type = "pred", terms = c("quartile[all]","condition") )

Summary

Betweenness centrality - Graphvar

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

main_anova_summary <- aov_ez("Subj_ID", "meanBWC", data.aggr,within =c("condition","threshold"))
DT::datatable(nice(main_anova_summary))
emmip(main_anova_summary, condition~threshold)

grouped_ggwithinstats(
    data             = data.aggr,
  x                = condition,
  y                = meanBWC,
  grouping.var     = threshold,
  type             = "p",
  bf.message = 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

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

# looking at threshold = 1
main_anova_summary <- aov_ez("Subj_ID", "meanBWC", data.aggr,within =c("condition","threshold"), between = "Age_cat")
EMM <- emmeans(main_anova_summary, ~condition:Age_cat, at=list(threshold="X1")) %>% contrast(interaction = c( "consec", "consec"))
DT::datatable(EMM %>% as.data.frame()) 
# looking at threshold = 3.5
main_anova_summary <- aov_ez("Subj_ID", "meanBWC", data.aggr,within =c("condition","threshold"), between = "Age_cat")
EMM <- emmeans(main_anova_summary, ~condition:Age_cat, at=list(threshold="X3.5")) %>% contrast(interaction = c( "consec", "consec"))
DT::datatable(EMM %>% as.data.frame()) 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanBWC", 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", "meanBWC", 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=="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="X1"))
  
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))
# focus on threshold=1
DT::datatable(nice(aov_ez("Subj_ID", "meanBWC", data.aggr %>% filter(threshold=="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="X1"))
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))
# focus on threshold=3.5
DT::datatable(nice(aov_ez("Subj_ID", "meanBWC", data.aggr %>% filter(threshold=="3.5"),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="X3.5"))
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", "meanBWC", 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==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=3.5
main_anova_summary <- aov_ez("Subj_ID", "meanBWC", data.aggr %>% filter(threshold==3.5),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))

Subgroup analysis: Placebo condition only

# testing

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,Age_cat) %>% summarise(meanBWC = mean(bwc)) %>% filter(condition=="Placebo") -> data.aggr.age

# regular anova
main_anova_summary_a <- aov_ez("Subj_ID", "meanBWC", data.aggr.age, between = "Age_cat")
nice(main_anova_summary_a)
## Anova Table (Type 3 tests)
## 
## Response: meanBWC
##    Effect    df    MSE    F  ges p.value
## 1 Age_cat 1, 25 686.01 0.08 .003    .777
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
afex_plot(main_anova_summary_a, "Age_cat")

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,Age_cat) %>% summarise(meanBWC = mean(bwc)) %>% filter(condition=="Placebo") -> data.aggr.age

ggbetweenstats(
    data             = data.aggr.age,
  x                = Age_cat,
  y                = meanBWC,
  type             = "p",
  bf.message = FALSE,
  )

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,Age) %>% summarise(meanBWC = mean(bwc)) %>% filter(condition=="Placebo") -> data.aggr.age

ggscatterstats(
  data = data.aggr.age, ## data frame from which variables are taken
  x = Age, ## predictor/independent variable
  y = meanBWC ## dependent variabl
)

grouped_ggbetweenstats(
    data             = data.aggr,
  x                = Age_cat,
  y                = meanBWC,
  grouping.var     = threshold,
  type             = "p",
  bf.message = FALSE,
  )

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanBWC", data.aggr,within =c("threshold"), between = "Age_cat")
DT::datatable(nice(main_anova_summary)) 
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)) %>% filter(condition=="Placebo") -> data.aggr

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

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

emmip(main_anova_summary, ~ Age_cat)

# average across threshold
main_anova_summary <- aov_ez("Subj_ID", "meanBWC", data.aggr,within =c("name"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  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", "meanBWC", data.aggr,within =c("name","threshold"), between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  Age_cat | name, at=list(threshold="X1")) %>% 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,within =c("name","threshold"), between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  Age_cat | name, at=list(threshold="X3.5")) %>% 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))

Summary

Global efficiency - Graphvar

Grouped by scale, condition, subject

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

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

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

Grouped by scale, condition, subject and age_cat

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

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","threshold"), between = "Age_cat")
DT::datatable(nice(main_anova_summary))
# anova averaging across threshold
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition"), between = "Age_cat")
DT::datatable(nice(main_anova_summary))
# looking at threshold = 1
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","threshold"), between = "Age_cat")
EMM <- emmeans(main_anova_summary, ~condition:Age_cat, at=list(threshold="X1")) %>% contrast(interaction = c( "consec", "consec"))
DT::datatable(EMM %>% as.data.frame()) 
# looking at threshold = 3.5
main_anova_summary <- aov_ez("Subj_ID", "meangeff", data.aggr,within =c("condition","threshold"), between = "Age_cat")
EMM <- emmeans(main_anova_summary, ~condition:Age_cat, at=list(threshold="X3.5")) %>% contrast(interaction = c( "consec", "consec"))
DT::datatable(EMM %>% as.data.frame()) 
# data.aggr %>% ungroup() %>%
#   nest_by(threshold) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meangeff", 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", "meangeff", 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 = "n_geff") %>% left_join(.,variables_ext) %>% group_by(threshold,Subj_ID,condition) -> data.aggr

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

emmip(main_anova_summary, condition~threshold)

emmip(main_anova_summary, condition~threshold   | network)

# anova averaging across threshold
main_anova_summary <- aov_ez("Subj_ID", "n_geff", data.aggr,within =c("network","condition"))
DT::datatable(nice(main_anova_summary))
main_anova_summary <- aov_ez("Subj_ID", "n_geff", data.aggr,within =c("network","condition"))
EMM <- emmeans(main_anova_summary, ~  condition | network)
  
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))
# looking at threshold = 1
main_anova_summary <- aov_ez("Subj_ID", "n_geff", data.aggr,within =c("network","condition","threshold"))
EMM <- emmeans(main_anova_summary, ~  condition | network, at=list(threshold="X1"))
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))
# looking at threshold = 3.5
main_anova_summary <- aov_ez("Subj_ID", "n_geff", data.aggr,within =c("network","condition","threshold"))
EMM <- emmeans(main_anova_summary, ~  condition | network, at=list(threshold="X3.5"))
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))

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

data %>% pivot_longer(
    cols = starts_with("ROI"),
    names_to = "ROI",
    values_to = "n_geff") %>% left_join(.,variables_ext) %>% group_by(threshold,Subj_ID,condition) -> data.aggr

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

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

emmip(main_anova_summary, network~Age_cat)

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

DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | network)  %>% 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", "n_geff", data.aggr,within =c("condition","network","threshold"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
# focus on threshold=1
main_anova_summary <- aov_ez("Subj_ID", "n_geff", data.aggr %>% filter(threshold==1),within =c("network","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | network) %>% 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", "n_geff", data.aggr %>% filter(threshold==3.5),within =c("network","condition"),between = "Age_cat")
DT::datatable(nice(main_anova_summary))
EMM <- emmeans(main_anova_summary, ~  condition:Age_cat | network) %>% 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))

Summary

Local efficiency - Graphvar

Grouped by scale, condition, subject

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

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

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

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

Grouped by scale, condition, subject and age_cat

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

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","threshold"), between = "Age_cat")
DT::datatable(nice(main_anova_summary))
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","threshold"), between = "Age_cat")
EMM <- emmeans(main_anova_summary, ~condition:Age_cat, at=list(threshold="X1")) %>% 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) %>%
#   mutate(Model = list(nice(aov_ez("Subj_ID", "meanleff", 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", "meanleff", 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 = "leff") %>% 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 <- afex::aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name","threshold"))
# DT::datatable(nice(main_anova_summary) )
# 
# # anova averaging across threshold
# main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name"))
# DT::datatable(nice(main_anova_summary))
# 
# # looking at threshold = 1
# main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr,within =c("condition","name","threshold"))
# EMM <- emmeans(main_anova_summary, ~condition | name, at=list(threshold="X1")) %>% 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", "meanleff", 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 = "leff") %>% 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==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=3.5
main_anova_summary <- aov_ez("Subj_ID", "meanleff", data.aggr %>% filter(threshold==3.5),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", "meanleff", data, within ="condition",between="Age_cat")))) %>% select(-data) %>% unnest(Model) -> anova_summary
# 
# DT::datatable(anova_summary)

Summary