graph measures

knitr::opts_chunk$set(echo = TRUE, cache = TRUE, fig.fullwidth=TRUE, out.width = "100%")
library(R.matlab)

## R.matlab v3.7.0 (2022-08-25 21:52:34 UTC) successfully loaded. See ?R.matlab for help.

## 
## Attaching package: 'R.matlab'

## The following objects are masked from 'package:base':
## 
##     getOption, isOpen

library(RCurl)
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2
## ──

## ✔ ggplot2 3.4.0      ✔ purrr   1.0.1 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.3.0      ✔ stringr 1.5.0 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::complete() masks RCurl::complete()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ dplyr::lag()      masks stats::lag()

library(ggstatsplot)

## You can cite this package as:
##      Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
##      Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167

library(afex)

## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: emmeans() for ANOVA models now uses model = 'multivariate' as default.
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
## 
## Attaching package: 'afex'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer

library(broom)
library(ComplexHeatmap)

## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.14.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##     genomic data. Bioinformatics 2016.
## 
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================

library(circlize)

## ========================================
## circlize version 0.4.15
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
##   in R. Bioinformatics 2014.
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(circlize))
## ========================================

library(gridExtra)

## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

library(report)
library(emmeans)
library(sjPlot)
library(DT)
library(igraph)

## 
## Attaching package: 'igraph'
## 
## The following object is masked from 'package:circlize':
## 
##     degree
## 
## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union
## 
## The following objects are masked from 'package:purrr':
## 
##     compose, simplify
## 
## The following object is masked from 'package:tidyr':
## 
##     crossing
## 
## The following object is masked from 'package:tibble':
## 
##     as_data_frame
## 
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## 
## The following object is masked from 'package:base':
## 
##     union

library(scales)

## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor

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_2/", 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)
# absolute transformed
result_matrix <- abs(result_matrix)

dim(result_matrix)

## [1] 394 394  56

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_2/", 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)
result_interferon_matrix <- abs(result_interferon_matrix)
# List to store all 2D matrices
placebo_matrices <- list()
# Loop through all .tsv files in the folder
files <- list.files(path = "conmats_filtered_2/", 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)
result_placebo_matrix <- abs(result_placebo_matrix)
# nnodes <- length(brainregions$X2)
# tri_pos <- which(upper.tri(matrix(nrow = nnodes, ncol = nnodes)), arr.ind = T)

mean_placebo_funcconn <- apply(result_placebo_matrix, c(1, 2), mean)
diag(mean_placebo_funcconn) <- 0
z_mean_placebo_funcconn <- atanh(mean_placebo_funcconn)

mean_interferon_funcconn <- apply(result_interferon_matrix, c(1, 2), mean)
diag(mean_interferon_funcconn) <- 0
z_mean_interferon_funcconn <- atanh(mean_interferon_funcconn)

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

## The automatically generated colors map from the 1^st and 99^th of the
## values in the matrix. There are outliers in the matrix whose patterns
## might be hidden by this color mapping. You can manually set the color
## to `col` argument.
## 
## Use `suppressMessages()` to turn off this message.

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

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

## The automatically generated colors map from the 1^st and 99^th of the
## values in the matrix. There are outliers in the matrix whose patterns
## might be hidden by this color mapping. You can manually set the color
## to `col` argument.
## 
## Use `suppressMessages()` to turn off this message.

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

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

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

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

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

## Warning: The heatmap has not been initialized. You might have different results
## if you repeatedly execute this function, e.g. when row_km/column_km was
## set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.

h1+h2+h3_v2_adjusted

## Warning: Heatmap/annotation names are duplicated: Fisher Z

## Warning: Heatmap/annotation names are duplicated: Fisher Z, Fisher Z

h1+h2

## Warning: Heatmap/annotation names are duplicated: Fisher Z

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)

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)  

DT::datatable(df_diff_pMinusI  %>% filter(abs_weight_decile >= 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()

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)

## Joining, by = "fileIndex"

Grouped by scale, condition, subject

Unthresholded

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

## `summarise()` has grouped output by 'Subj_ID'. You can override using the
## `.groups` argument.

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

knitr::kable(anova_summary)

Effect	df	MSE	F	ges	p.value
condition	1, 27	0.00	16.11 ***	.096	<.001

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

## Adding missing grouping variables: `Subj_ID`

Absolute weight

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

## `summarise()` has grouped output by 'Subj_ID'. You can override using the
## `.groups` argument.

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

knitr::kable(anova_summary)

Effect	df	MSE	F	ges	p.value
condition	1, 27	0.00	16.11 ***	.096	<.001

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

## Adding missing grouping variables: `Subj_ID`

Proportional thresholding based on abs weight

# total edges per session: 77421
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(77421*x[i]),weight) %>% dplyr::summarise(meanFuncCon = mean(weight)),
    x                = condition,
    y                = meanFuncCon,
    type             = "p",
    title =x[i] ,
    bf.message = FALSE)
  
    plots_list[[i]] <- p

}

## `summarise()` has grouped output by 'Subj_ID'. You can override using the
## `.groups` argument.
## Adding missing grouping variables: `Subj_ID`
## `summarise()` has grouped output by 'Subj_ID'. You can override using the
## `.groups` argument.
## Adding missing grouping variables: `Subj_ID`
## `summarise()` has grouped output by 'Subj_ID'. You can override using the
## `.groups` argument.
## Adding missing grouping variables: `Subj_ID`
## `summarise()` has grouped output by 'Subj_ID'. You can override using the
## `.groups` argument.
## Adding missing grouping variables: `Subj_ID`
## `summarise()` has grouped output by 'Subj_ID'. You can override using the
## `.groups` argument.
## Adding missing grouping variables: `Subj_ID`
## `summarise()` has grouped output by 'Subj_ID'. You can override using the
## `.groups` argument.
## Adding missing grouping variables: `Subj_ID`

# 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

## `summarise()` has grouped output by 'Subj_ID', 'condition'. You can override
## using the `.groups` argument.

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

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: Age_cat

knitr::kable(anova_summary)

Effect	df	MSE	F	ges	p.value
Age_cat	1, 26	0.00	4.10 +	.114	.053
condition	1, 26	0.00	15.85 ***	.103	<.001
Age_cat:condition	1, 26	0.00	2.11	.015	.158

Grouped by scale, condition, subject and network:yeo

• Unthresholded

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)
knitr::kable(df_counts)

	Interferon	Placebo
Default2Default	92988	92988
Default2Dorsal Attention	101024	101024
Default2Frontoparietal	103320	103320
Default2Limbic	66584	66584
Default2Somatomotor	121688	121688
Default2Subcortical	78064	78064
Default2Ventral Attention	110208	110208
Default2Visual	135464	135464
Dorsal Attention2Dorsal Attention	26488	26488
Dorsal Attention2Frontoparietal	55440	55440
Dorsal Attention2Limbic	35728	35728
Dorsal Attention2Somatomotor	65296	65296
Dorsal Attention2Subcortical	41888	41888
Dorsal Attention2Ventral Attention	59136	59136
Dorsal Attention2Visual	72688	72688
Frontoparietal2Frontoparietal	27720	27720
Frontoparietal2Limbic	36540	36540
Frontoparietal2Somatomotor	66780	66780
Frontoparietal2Subcortical	42840	42840
Frontoparietal2Ventral Attention	60480	60480
Frontoparietal2Visual	74340	74340
Limbic2Limbic	11368	11368
Limbic2Somatomotor	43036	43036
Limbic2Subcortical	27608	27608
Limbic2Ventral Attention	38976	38976
Limbic2Visual	47908	47908
Somatomotor2Somatomotor	38584	38584
Somatomotor2Subcortical	50456	50456
Somatomotor2Ventral Attention	71232	71232
Somatomotor2Visual	87556	87556
Subcortical2Subcortical	15708	15708
Subcortical2Ventral Attention	45696	45696
Subcortical2Visual	56168	56168
Ventral Attention2Ventral Attention	31584	31584
Ventral Attention2Visual	79296	79296
Visual2Visual	47908	47908

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

## `summarise()` has grouped output by 'Subj_ID', 'condition'. You can override
## using the `.groups` argument.

# 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"))
knitr::kable(nice(main_anova_summary))

Effect	df	MSE	F	ges	p.value
condition	1, 27	0.00	18.57 ***	.023	<.001
from2to	35, 945	0.00	142.45 ***	.745	<.001
condition:from2to	35, 945	0.00	1.81 **	.014	.003

emmip(main_anova_summary, from2to ~ condition)

## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 36. Consider
## specifying shapes manually if you must have them.

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)

# corrected and uncorrected
knitr::kable(left_join(EMM.unadj,EMM.fdr))

## Joining, by = c("contrast", "from2to", "estimate", "SE", "df", "t.ratio")

contrast	from2to	estimate	SE	df	t.ratio	p.value.unadj	p.value.fdr
Interferon - Placebo	Default2Default	-0.0190538	0.0099067	27	-1.9233255	0.0650437	0.1204645
Interferon - Placebo	Default2Dorsal.Attention	-0.0197879	0.0087745	27	-2.2551543	0.0324366	0.1204645
Interferon - Placebo	Default2Frontoparietal	-0.0111147	0.0052667	27	-2.1103712	0.0442322	0.1204645
Interferon - Placebo	Default2Limbic	-0.0125644	0.0057663	27	-2.1789163	0.0382382	0.1204645
Interferon - Placebo	Default2Somatomotor	-0.0118285	0.0057666	27	-2.0512036	0.0500652	0.1204645
Interferon - Placebo	Default2Subcortical	-0.0063649	0.0046541	27	-1.3676128	0.1827103	0.2529835
Interferon - Placebo	Default2Ventral.Attention	-0.0206605	0.0081451	27	-2.5365481	0.0172871	0.1037228
Interferon - Placebo	Default2Visual	-0.0182318	0.0059788	27	-3.0494058	0.0050885	0.0457966
Interferon - Placebo	Dorsal.Attention2Dorsal.Attention	-0.0161942	0.0089819	27	-1.8029714	0.0825655	0.1415408
Interferon - Placebo	Dorsal.Attention2Frontoparietal	-0.0049945	0.0042842	27	-1.1658072	0.2538880	0.3264274
Interferon - Placebo	Dorsal.Attention2Limbic	-0.0023448	0.0038210	27	-0.6136473	0.5445842	0.6126572
Interferon - Placebo	Dorsal.Attention2Somatomotor	-0.0138993	0.0061182	27	-2.2717887	0.0312817	0.1204645
Interferon - Placebo	Dorsal.Attention2Subcortical	-0.0073025	0.0037568	27	-1.9437780	0.0624123	0.1204645
Interferon - Placebo	Dorsal.Attention2Ventral.Attention	-0.0129097	0.0067620	27	-1.9091434	0.0669247	0.1204645
Interferon - Placebo	Dorsal.Attention2Visual	-0.0161409	0.0073717	27	-2.1895855	0.0373738	0.1204645
Interferon - Placebo	Frontoparietal2Frontoparietal	-0.0046194	0.0066871	27	-0.6907855	0.4955967	0.5755317
Interferon - Placebo	Frontoparietal2Limbic	-0.0034150	0.0042632	27	-0.8010602	0.4300901	0.5339050
Interferon - Placebo	Frontoparietal2Somatomotor	-0.0117760	0.0058918	27	-1.9987264	0.0557981	0.1204645
Interferon - Placebo	Frontoparietal2Subcortical	-0.0076420	0.0038586	27	-1.9805126	0.0579189	0.1204645
Interferon - Placebo	Frontoparietal2Ventral.Attention	-0.0160793	0.0051668	27	-3.1120378	0.0043580	0.0457966
Interferon - Placebo	Frontoparietal2Visual	-0.0106668	0.0053892	27	-1.9792827	0.0580646	0.1204645
Interferon - Placebo	Limbic2Limbic	-0.0009078	0.0187898	27	-0.0483149	0.9618209	0.9618209
Interferon - Placebo	Limbic2Somatomotor	-0.0022179	0.0050063	27	-0.4430188	0.6612815	0.7213980
Interferon - Placebo	Limbic2Subcortical	0.0025821	0.0080916	27	0.3191046	0.7521029	0.7963443
Interferon - Placebo	Limbic2Ventral.Attention	-0.0106189	0.0049624	27	-2.1398783	0.0415560	0.1204645
Interferon - Placebo	Limbic2Visual	-0.0005289	0.0053970	27	-0.0980026	0.9226542	0.9490157
Interferon - Placebo	Somatomotor2Somatomotor	-0.0453309	0.0139491	27	-3.2497355	0.0030886	0.0457966
Interferon - Placebo	Somatomotor2Subcortical	-0.0042145	0.0056045	27	-0.7519776	0.4585702	0.5502842
Interferon - Placebo	Somatomotor2Ventral.Attention	-0.0231803	0.0080127	27	-2.8929332	0.0074575	0.0536943
Interferon - Placebo	Somatomotor2Visual	-0.0098405	0.0071267	27	-1.3807841	0.1786691	0.2529835
Interferon - Placebo	Subcortical2Subcortical	-0.0129955	0.0087041	27	-1.4930403	0.1470212	0.2301201
Interferon - Placebo	Subcortical2Ventral.Attention	-0.0102303	0.0059058	27	-1.7322523	0.0946397	0.1548649
Interferon - Placebo	Subcortical2Visual	-0.0077813	0.0053641	27	-1.4506368	0.1584003	0.2376005
Interferon - Placebo	Ventral.Attention2Ventral.Attention	-0.0402414	0.0124304	27	-3.2373289	0.0031866	0.0457966
Interferon - Placebo	Ventral.Attention2Visual	-0.0114958	0.0059146	27	-1.9436455	0.0624290	0.1204645
Interferon - Placebo	Visual2Visual	-0.0153088	0.0115526	27	-1.3251387	0.1962305	0.2616407

# 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
# 
# knitr::kable(anova_summary)

• 10% threshold

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

## `summarise()` has grouped output by 'Subj_ID', 'condition'. You can override
## using the `.groups` argument.

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

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

	Interferon	Placebo
Default2Default	92988	92988
Default2Dorsal Attention	101024	101024
Default2Frontoparietal	103320	103320
Default2Limbic	66584	66584
Default2Somatomotor	121688	121688
Default2Subcortical	78064	78064
Default2Ventral Attention	110208	110208
Default2Visual	135464	135464
Dorsal Attention2Dorsal Attention	26488	26488
Dorsal Attention2Frontoparietal	55440	55440
Dorsal Attention2Limbic	35728	35728
Dorsal Attention2Somatomotor	65296	65296
Dorsal Attention2Subcortical	41888	41888
Dorsal Attention2Ventral Attention	59136	59136
Dorsal Attention2Visual	72688	72688
Frontoparietal2Frontoparietal	27720	27720
Frontoparietal2Limbic	36540	36540
Frontoparietal2Somatomotor	66780	66780
Frontoparietal2Subcortical	42840	42840
Frontoparietal2Ventral Attention	60480	60480
Frontoparietal2Visual	74340	74340
Limbic2Limbic	11368	11368
Limbic2Somatomotor	43036	43036
Limbic2Subcortical	27608	27608
Limbic2Ventral Attention	38976	38976
Limbic2Visual	47908	47908
Somatomotor2Somatomotor	38584	38584
Somatomotor2Subcortical	50456	50456
Somatomotor2Ventral Attention	71232	71232
Somatomotor2Visual	87556	87556
Subcortical2Subcortical	15708	15708
Subcortical2Ventral Attention	45696	45696
Subcortical2Visual	56168	56168
Ventral Attention2Ventral Attention	31584	31584
Ventral Attention2Visual	79296	79296
Visual2Visual	47908	47908

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","from2to"))

## Warning: Missing values for following ID(s):
## 322030, 1021019
## Removing those cases from the analysis.

knitr::kable(nice(main_anova_summary))

Effect	df	MSE	F	ges	p.value
condition	1, 25	0.03	18.22 ***	.094	<.001
from2to	35, 875	0.00	131.20 ***	.507	<.001
condition:from2to	35, 875	0.00	1.47 *	.005	.041

emmip(main_anova_summary, from2to ~ condition)

## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 36. Consider
## specifying shapes manually if you must have them.

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)

# corrected and uncorrected
knitr::kable(left_join(EMM.unadj,EMM.fdr))

## Joining, by = c("contrast", "from2to", "estimate", "SE", "df", "t.ratio")

contrast	from2to	estimate	SE	df	t.ratio	p.value.unadj	p.value.fdr
Interferon - Placebo	Default2Default	-0.0370918	0.0106146	25	-3.494420	0.0017904	0.0033107
Interferon - Placebo	Default2Dorsal.Attention	-0.0343658	0.0091502	25	-3.755758	0.0009251	0.0022136
Interferon - Placebo	Default2Frontoparietal	-0.0336530	0.0093327	25	-3.605918	0.0013525	0.0027049
Interferon - Placebo	Default2Limbic	-0.0388473	0.0092591	25	-4.195559	0.0002994	0.0016804
Interferon - Placebo	Default2Somatomotor	-0.0345566	0.0101964	25	-3.389110	0.0023291	0.0039928
Interferon - Placebo	Default2Subcortical	-0.0237063	0.0087104	25	-2.721619	0.0116590	0.0139908
Interferon - Placebo	Default2Ventral.Attention	-0.0367185	0.0088229	25	-4.161717	0.0003267	0.0016804
Interferon - Placebo	Default2Visual	-0.0302748	0.0086905	25	-3.483662	0.0018393	0.0033107
Interferon - Placebo	Dorsal.Attention2Dorsal.Attention	-0.0360476	0.0099240	25	-3.632351	0.0012651	0.0026790
Interferon - Placebo	Dorsal.Attention2Frontoparietal	-0.0304166	0.0081121	25	-3.749518	0.0009399	0.0022136
Interferon - Placebo	Dorsal.Attention2Limbic	-0.0305318	0.0105145	25	-2.903763	0.0075995	0.0101327
Interferon - Placebo	Dorsal.Attention2Somatomotor	-0.0525249	0.0105019	25	-5.001474	0.0000371	0.0007591
Interferon - Placebo	Dorsal.Attention2Subcortical	-0.0406767	0.0107092	25	-3.798300	0.0008302	0.0022136
Interferon - Placebo	Dorsal.Attention2Ventral.Attention	-0.0301320	0.0106305	25	-2.834481	0.0089536	0.0111148
Interferon - Placebo	Dorsal.Attention2Visual	-0.0339087	0.0110951	25	-3.056179	0.0052736	0.0073020
Interferon - Placebo	Frontoparietal2Frontoparietal	-0.0334120	0.0089539	25	-3.731562	0.0009838	0.0022136
Interferon - Placebo	Frontoparietal2Limbic	-0.0322978	0.0127136	25	-2.540417	0.0176599	0.0198674
Interferon - Placebo	Frontoparietal2Somatomotor	-0.0337240	0.0107344	25	-3.141678	0.0042853	0.0064279
Interferon - Placebo	Frontoparietal2Subcortical	-0.0231354	0.0086892	25	-2.662549	0.0133652	0.0155209
Interferon - Placebo	Frontoparietal2Ventral.Attention	-0.0336492	0.0101485	25	-3.315675	0.0027947	0.0045732
Interferon - Placebo	Frontoparietal2Visual	-0.0322497	0.0082662	25	-3.901402	0.0006379	0.0021342
Interferon - Placebo	Limbic2Limbic	-0.0276898	0.0144947	25	-1.910331	0.0676259	0.0695581
Interferon - Placebo	Limbic2Somatomotor	-0.0418527	0.0104916	25	-3.989171	0.0005094	0.0020377
Interferon - Placebo	Limbic2Subcortical	-0.0322310	0.0130521	25	-2.469414	0.0207145	0.0225976
Interferon - Placebo	Limbic2Ventral.Attention	-0.0391735	0.0088808	25	-4.411016	0.0001715	0.0014291
Interferon - Placebo	Limbic2Visual	-0.0195918	0.0120306	25	-1.628501	0.1159571	0.1159571
Interferon - Placebo	Somatomotor2Somatomotor	-0.0585847	0.0118298	25	-4.952284	0.0000422	0.0007591
Interferon - Placebo	Somatomotor2Subcortical	-0.0308219	0.0100648	25	-3.062352	0.0051955	0.0073020
Interferon - Placebo	Somatomotor2Ventral.Attention	-0.0430604	0.0092785	25	-4.640881	0.0000945	0.0011342
Interferon - Placebo	Somatomotor2Visual	-0.0315085	0.0099448	25	-3.168341	0.0040153	0.0062848
Interferon - Placebo	Subcortical2Subcortical	-0.0380115	0.0101531	25	-3.743838	0.0009536	0.0022136
Interferon - Placebo	Subcortical2Ventral.Attention	-0.0349065	0.0089669	25	-3.892813	0.0006521	0.0021342
Interferon - Placebo	Subcortical2Visual	-0.0251376	0.0112288	25	-2.238665	0.0343201	0.0363389
Interferon - Placebo	Ventral.Attention2Ventral.Attention	-0.0508286	0.0116724	25	-4.354585	0.0001985	0.0014291
Interferon - Placebo	Ventral.Attention2Visual	-0.0405376	0.0099942	25	-4.056129	0.0004289	0.0019300
Interferon - Placebo	Visual2Visual	-0.0400252	0.0139114	25	-2.877159	0.0080947	0.0104075

• 35% threshold

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

## `summarise()` has grouped output by 'Subj_ID', 'condition'. You can override
## using the `.groups` argument.

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

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

	Interferon	Placebo
Default2Default	92988	92988
Default2Dorsal Attention	101024	101024
Default2Frontoparietal	103320	103320
Default2Limbic	66584	66584
Default2Somatomotor	121688	121688
Default2Subcortical	78064	78064
Default2Ventral Attention	110208	110208
Default2Visual	135464	135464
Dorsal Attention2Dorsal Attention	26488	26488
Dorsal Attention2Frontoparietal	55440	55440
Dorsal Attention2Limbic	35728	35728
Dorsal Attention2Somatomotor	65296	65296
Dorsal Attention2Subcortical	41888	41888
Dorsal Attention2Ventral Attention	59136	59136
Dorsal Attention2Visual	72688	72688
Frontoparietal2Frontoparietal	27720	27720
Frontoparietal2Limbic	36540	36540
Frontoparietal2Somatomotor	66780	66780
Frontoparietal2Subcortical	42840	42840
Frontoparietal2Ventral Attention	60480	60480
Frontoparietal2Visual	74340	74340
Limbic2Limbic	11368	11368
Limbic2Somatomotor	43036	43036
Limbic2Subcortical	27608	27608
Limbic2Ventral Attention	38976	38976
Limbic2Visual	47908	47908
Somatomotor2Somatomotor	38584	38584
Somatomotor2Subcortical	50456	50456
Somatomotor2Ventral Attention	71232	71232
Somatomotor2Visual	87556	87556
Subcortical2Subcortical	15708	15708
Subcortical2Ventral Attention	45696	45696
Subcortical2Visual	56168	56168
Ventral Attention2Ventral Attention	31584	31584
Ventral Attention2Visual	79296	79296
Visual2Visual	47908	47908

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","from2to"))
knitr::kable(nice(main_anova_summary))

Effect	df	MSE	F	ges	p.value
condition	1, 27	0.02	15.81 ***	.061	<.001
from2to	35, 945	0.00	178.96 ***	.669	<.001
condition:from2to	35, 945	0.00	2.21 ***	.010	<.001

emmip(main_anova_summary, from2to ~ condition)

## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 36. Consider
## specifying shapes manually if you must have them.

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)

# corrected and uncorrected
knitr::kable(left_join(EMM.unadj,EMM.fdr))

## Joining, by = c("contrast", "from2to", "estimate", "SE", "df", "t.ratio")

contrast	from2to	estimate	SE	df	t.ratio	p.value.unadj	p.value.fdr
Interferon - Placebo	Default2Default	-0.0253441	0.0097983	27	-2.5865862	0.0154048	0.0241119
Interferon - Placebo	Default2Dorsal.Attention	-0.0230500	0.0084613	27	-2.7241696	0.0111662	0.0194059
Interferon - Placebo	Default2Frontoparietal	-0.0224312	0.0069379	27	-3.2331636	0.0032201	0.0096403
Interferon - Placebo	Default2Limbic	-0.0211105	0.0068733	27	-3.0713750	0.0048199	0.0096403
Interferon - Placebo	Default2Somatomotor	-0.0226373	0.0071263	27	-3.1765936	0.0037106	0.0096403
Interferon - Placebo	Default2Subcortical	-0.0136979	0.0063012	27	-2.1738477	0.0386551	0.0434870
Interferon - Placebo	Default2Ventral.Attention	-0.0279336	0.0084437	27	-3.3081957	0.0026649	0.0096403
Interferon - Placebo	Default2Visual	-0.0228276	0.0069819	27	-3.2695529	0.0029382	0.0096403
Interferon - Placebo	Dorsal.Attention2Dorsal.Attention	-0.0227624	0.0091610	27	-2.4846961	0.0194607	0.0280233
Interferon - Placebo	Dorsal.Attention2Frontoparietal	-0.0186098	0.0059301	27	-3.1382205	0.0040834	0.0096403
Interferon - Placebo	Dorsal.Attention2Limbic	-0.0174400	0.0064156	27	-2.7183770	0.0113201	0.0194059
Interferon - Placebo	Dorsal.Attention2Somatomotor	-0.0289719	0.0089641	27	-3.2319977	0.0032295	0.0096403
Interferon - Placebo	Dorsal.Attention2Subcortical	-0.0172189	0.0055447	27	-3.1054797	0.0044296	0.0096403
Interferon - Placebo	Dorsal.Attention2Ventral.Attention	-0.0216781	0.0085433	27	-2.5374466	0.0172515	0.0258773
Interferon - Placebo	Dorsal.Attention2Visual	-0.0268588	0.0080606	27	-3.3320940	0.0025083	0.0096403
Interferon - Placebo	Frontoparietal2Frontoparietal	-0.0163363	0.0074496	27	-2.1928967	0.0371091	0.0430944
Interferon - Placebo	Frontoparietal2Limbic	-0.0150786	0.0070983	27	-2.1242587	0.0429541	0.0468590
Interferon - Placebo	Frontoparietal2Somatomotor	-0.0223249	0.0075321	27	-2.9639615	0.0062752	0.0118899
Interferon - Placebo	Frontoparietal2Subcortical	-0.0180557	0.0058787	27	-3.0713549	0.0048201	0.0096403
Interferon - Placebo	Frontoparietal2Ventral.Attention	-0.0255342	0.0073728	27	-3.4633012	0.0017949	0.0096403
Interferon - Placebo	Frontoparietal2Visual	-0.0230964	0.0057036	27	-4.0494687	0.0003882	0.0046588
Interferon - Placebo	Limbic2Limbic	-0.0098725	0.0140553	27	-0.7024046	0.4884388	0.4884388
Interferon - Placebo	Limbic2Somatomotor	-0.0159618	0.0068728	27	-2.3224769	0.0279886	0.0359853
Interferon - Placebo	Limbic2Subcortical	-0.0130999	0.0087638	27	-1.4947847	0.1465675	0.1507552
Interferon - Placebo	Limbic2Ventral.Attention	-0.0236318	0.0060683	27	-3.8943228	0.0005849	0.0052645
Interferon - Placebo	Limbic2Visual	-0.0112084	0.0057258	27	-1.9575273	0.0606960	0.0642663
Interferon - Placebo	Somatomotor2Somatomotor	-0.0556599	0.0112351	27	-4.9541198	0.0000345	0.0012403
Interferon - Placebo	Somatomotor2Subcortical	-0.0161122	0.0071445	27	-2.2551959	0.0324337	0.0402625
Interferon - Placebo	Somatomotor2Ventral.Attention	-0.0335587	0.0078437	27	-4.2784157	0.0002111	0.0038003
Interferon - Placebo	Somatomotor2Visual	-0.0210954	0.0086904	27	-2.4274411	0.0221516	0.0306715
Interferon - Placebo	Subcortical2Subcortical	-0.0320967	0.0103876	27	-3.0898895	0.0046041	0.0096403
Interferon - Placebo	Subcortical2Ventral.Attention	-0.0225049	0.0070271	27	-3.2025849	0.0034770	0.0096403
Interferon - Placebo	Subcortical2Visual	-0.0125308	0.0053156	27	-2.3573474	0.0259097	0.0345463
Interferon - Placebo	Ventral.Attention2Ventral.Attention	-0.0424107	0.0119656	27	-3.5443964	0.0014570	0.0096403
Interferon - Placebo	Ventral.Attention2Visual	-0.0168930	0.0075464	27	-2.2385479	0.0336280	0.0403535
Interferon - Placebo	Visual2Visual	-0.0291114	0.0110810	27	-2.6271532	0.0140205	0.0229426

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

## `summarise()` has grouped output by 'Subj_ID', 'condition', 'from2to'. You can
## override using the `.groups` argument.

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

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: Age_cat

knitr::kable(nice(main_anova_summary))

Effect	df	MSE	F	ges	p.value
Age_cat	1, 26	0.02	3.38 +	.023	.077
condition	1, 26	0.00	17.77 ***	.024	<.001
Age_cat:condition	1, 26	0.00	0.73	.001	.402
from2to	35, 910	0.00	143.69 ***	.756	<.001
Age_cat:from2to	35, 910	0.00	1.57 *	.033	.020
condition:from2to	35, 910	0.00	1.76 **	.015	.005
Age_cat:condition:from2to	35, 910	0.00	1.50 *	.013	.032

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)

# corrected and uncorrected
knitr::kable(left_join(EMM.unadj,EMM.fdr))

## Joining, by = c("contrast", "estimate", "SE", "df", "t.ratio")

contrast	estimate	SE	df	t.ratio	p.value.unadj	p.value.fdr
Interferon - Placebo	-0.0120455	0.0028576	26	-4.215204	0.0002662	0.0002662

# 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
# 
# knitr::kable(anova_summary)

Summary

Node degree

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

## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID'. You can override
## using the `.groups` argument.

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

knitr::kable(nice(anova_summary))

Effect	df	MSE	F	ges	p.value
condition	1, 27	2.13	23.23 ***	.287	<.001
threshold	1, 27	0.34	796630.79 ***	>.999	<.001
condition:threshold	1, 27	0.50	48.77 ***	.165	<.001

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

knitr::kable(anova_summary)

Effect	df	MSE	F	ges	p.value
Age_cat	1, 26	1.62	0.41	.005	.530
condition	1, 26	2.21	22.17 ***	.287	<.001
Age_cat:condition	1, 26	2.21	0.03	<.001	.873
threshold	1, 26	0.34	784321.09 ***	>.999	<.001
Age_cat:threshold	1, 26	0.34	0.72	.002	.405
condition:threshold	1, 26	0.51	48.21 ***	.167	<.001
Age_cat:condition:threshold	1, 26	0.51	0.46	.002	.505

# 
# 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
# 
# knitr::kable(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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID', 'condition'. You
## can override using the `.groups` argument.

# regular anova
main_anova_summary <- aov_ez("Subj_ID", "meanDegree", data.aggr,within =c("condition","name","threshold"))
knitr::kable(nice(main_anova_summary))

Effect	df	MSE	F	ges	p.value
condition	1, 27	19.11	22.85 ***	.030	<.001
name	4.23, 114.25	57.95	3.43 **	.056	.010
threshold	1, 27	3.84	563946.36 ***	.993	<.001
condition:name	4.43, 119.69	31.14	0.95	.009	.441
condition:threshold	1, 27	4.41	55.10 ***	.017	<.001
name:threshold	5.27, 142.38	12.32	0.64	.003	.681
condition:name:threshold	4.85, 130.87	7.59	1.32	.003	.260

#emmip(main_anova_summary, condition~threshold)

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

## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.

knitr::kable(nice(main_anova_summary))

Effect	df	MSE	F	ges	p.value
condition	1, 27	9.56	22.85 ***	.037	<.001
name	4.23, 114.25	28.97	3.43 **	.069	.010
condition:name	4.43, 119.69	15.57	0.95	.012	.441

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
knitr::kable(left_join(EMM.unadj,EMM.fdr))

## Joining, by = c("condition_consec", "name", "estimate", "SE", "df", "t.ratio")

condition_consec	name	estimate	SE	df	t.ratio	p.value.unadj	p.value.fdr
Placebo - Interferon	Default	0.8456010	0.6432030	27	1.3146722	0.1996784	0.2282038
Placebo - Interferon	Dorsal.Attention	0.9155844	0.6940650	27	1.3191623	0.1981935	0.2282038
Placebo - Interferon	Frontoparietal	1.4936508	0.8938850	27	1.6709652	0.1062816	0.2070502
Placebo - Interferon	Limbic	0.5055419	1.1513801	27	0.4390747	0.6641013	0.6641013
Placebo - Interferon	Somatomotor	1.7661725	0.5739225	27	3.0773710	0.0047490	0.0334381
Placebo - Interferon	Subcortical	3.1759454	1.1160925	27	2.8455934	0.0083595	0.0334381
Placebo - Interferon	Ventral.Attention	1.3188244	0.8431058	27	1.5642455	0.1294063	0.2070502
Placebo - Interferon	Visual	1.1501211	0.5541657	27	2.0754100	0.0476011	0.1269362

# focus on threshold=3.5
knitr::kable(nice(aov_ez("Subj_ID", "meanDegree", data.aggr %>% filter(threshold=="3.5"),within =c("condition","name"))))

Effect	df	MSE	F	ges	p.value
condition	1, 27	16.99	39.17 ***	.060	<.001
name	4.49, 121.34	51.48	1.96 +	.042	.097
condition:name	4.41, 119.08	28.46	1.10	.013	.360

main_anova_summary <- aov_ez("Subj_ID", "meanDegree", data.aggr,within =c("condition","name","threshold"))
knitr::kable(nice(main_anova_summary))

Effect	df	MSE	F	ges	p.value
condition	1, 27	19.11	22.85 ***	.030	<.001
name	4.23, 114.25	57.95	3.43 **	.056	.010
threshold	1, 27	3.84	563946.36 ***	.993	<.001
condition:name	4.43, 119.69	31.14	0.95	.009	.441
condition:threshold	1, 27	4.41	55.10 ***	.017	<.001
name:threshold	5.27, 142.38	12.32	0.64	.003	.681
condition:name:threshold	4.85, 130.87	7.59	1.32	.003	.260

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
knitr::kable(left_join(EMM.unadj,EMM.fdr))

## Joining, by = c("contrast", "name", "estimate", "SE", "df", "t.ratio")

contrast	name	estimate	SE	df	t.ratio	p.value.unadj	p.value.fdr
Interferon - Placebo	Default	-1.338850	0.8325453	27	-1.608141	0.1194368	0.1364992
Interferon - Placebo	Dorsal.Attention	-1.918831	0.9408427	27	-2.039481	0.0512986	0.0683981
Interferon - Placebo	Frontoparietal	-2.447619	1.1422038	27	-2.142892	0.0412910	0.0683981
Interferon - Placebo	Limbic	-1.873153	1.6363108	27	-1.144741	0.2623630	0.2623630
Interferon - Placebo	Somatomotor	-2.337601	0.8380524	27	-2.789326	0.0095655	0.0255081
Interferon - Placebo	Subcortical	-5.243697	1.4819580	27	-3.538358	0.0014799	0.0118389
Interferon - Placebo	Ventral.Attention	-2.363839	1.1455058	27	-2.063577	0.0487919	0.0683981
Interferon - Placebo	Visual	-1.977603	0.6384955	27	-3.097286	0.0045205	0.0180819

# 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
# 
# knitr::kable(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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID', 'condition',
## 'name'. You can override using the `.groups` argument.

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

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: Age_cat

knitr::kable(nice(main_anova_summary)) 

Effect	df	MSE	F	ges	p.value
Age_cat	1, 26	16.43	0.67	<.001	.419
condition	1, 26	19.84	21.84 ***	.034	<.001
Age_cat:condition	1, 26	19.84	0.01	<.001	.929
name	5.26, 136.79	38.90	3.96 **	.061	.002
Age_cat:name	5.26, 136.79	38.90	6.35 ***	.095	<.001
threshold	1, 26	3.85	559124.19 ***	.994	<.001
Age_cat:threshold	1, 26	3.85	0.90	<.001	.351
condition:name	4.83, 125.51	26.76	0.85	.009	.514
Age_cat:condition:name	4.83, 125.51	26.76	2.85 *	.029	.019
condition:threshold	1, 26	4.54	53.61 ***	.019	<.001
Age_cat:condition:threshold	1, 26	4.54	0.18	<.001	.671
name:threshold	5.19, 134.85	11.77	0.73	.004	.607
Age_cat:name:threshold	5.19, 134.85	11.77	2.74 *	.013	.020
condition:name:threshold	4.93, 128.08	7.38	1.27	.004	.281
Age_cat:condition:name:threshold	4.93, 128.08	7.38	1.35	.004	.248

emmip(main_anova_summary, name~condition | Age_cat)

## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 8. Consider
## specifying shapes manually if you must have them.

emmip(main_anova_summary, name~Age_cat)

## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 8. Consider
## specifying shapes manually if you must have them.

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

## Converting to factor: Age_cat

## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.

## Contrasts set to contr.sum for the following variables: Age_cat

knitr::kable(nice(main_anova_summary))

Effect	df	MSE	F	ges	p.value
Age_cat	1, 26	8.21	0.67	.001	.419
condition	1, 26	9.92	21.84 ***	.043	<.001
Age_cat:condition	1, 26	9.92	0.01	<.001	.929
name	5.26, 136.79	19.45	3.96 **	.078	.002
Age_cat:name	5.26, 136.79	19.45	6.35 ***	.119	<.001
condition:name	4.83, 125.51	13.38	0.85	.011	.514
Age_cat:condition:name	4.83, 125.51	13.38	2.85 *	.037	.019

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
knitr::kable(left_join(EMM.unadj,EMM.fdr))

## Joining, by = c("condition_consec", "Age_cat_consec", "name", "estimate", "SE",
## "df", "t.ratio")

condition_consec	Age_cat_consec	name	estimate	SE	df	t.ratio	p.value.unadj	p.value.fdr
Placebo - Interferon	Y - O	Default	-1.1646341	1.294269	26	-0.8998393	0.3764650	0.5122366
Placebo - Interferon	Y - O	Dorsal.Attention	1.2571096	1.396602	26	0.9001203	0.3763182	0.5122366
Placebo - Interferon	Y - O	Frontoparietal	-1.5929345	1.799577	26	-0.8851717	0.3841774	0.5122366
Placebo - Interferon	Y - O	Limbic	-5.1687887	2.123050	26	-2.4346054	0.0220788	0.1766300
Placebo - Interferon	Y - O	Somatomotor	0.4846638	1.168848	26	0.4146509	0.6817989	0.6817989
Placebo - Interferon	Y - O	Subcortical	4.1842383	2.127780	26	1.9664811	0.0600034	0.2400137
Placebo - Interferon	Y - O	Ventral.Attention	1.2622329	1.704855	26	0.7403754	0.4657061	0.5322356
Placebo - Interferon	Y - O	Visual	1.1697523	1.108855	26	1.0549194	0.3011669	0.5122366

# focus on threshold=3.5
main_anova_summary <- aov_ez("Subj_ID", "meanDegree", data.aggr,within =c("condition","name","threshold"),between = "Age_cat")

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: Age_cat

knitr::kable(nice(main_anova_summary))

Effect	df	MSE	F	ges	p.value
Age_cat	1, 26	16.43	0.67	<.001	.419
condition	1, 26	19.84	21.84 ***	.034	<.001
Age_cat:condition	1, 26	19.84	0.01	<.001	.929
name	5.26, 136.79	38.90	3.96 **	.061	.002
Age_cat:name	5.26, 136.79	38.90	6.35 ***	.095	<.001
threshold	1, 26	3.85	559124.19 ***	.994	<.001
Age_cat:threshold	1, 26	3.85	0.90	<.001	.351
condition:name	4.83, 125.51	26.76	0.85	.009	.514
Age_cat:condition:name	4.83, 125.51	26.76	2.85 *	.029	.019
condition:threshold	1, 26	4.54	53.61 ***	.019	<.001
Age_cat:condition:threshold	1, 26	4.54	0.18	<.001	.671
name:threshold	5.19, 134.85	11.77	0.73	.004	.607
Age_cat:name:threshold	5.19, 134.85	11.77	2.74 *	.013	.020
condition:name:threshold	4.93, 128.08	7.38	1.27	.004	.281
Age_cat:condition:name:threshold	4.93, 128.08	7.38	1.35	.004	.248

EMM <- emmeans(main_anova_summary, ~  condition: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
knitr::kable(left_join(EMM.unadj,EMM.fdr))

## Joining, by = c("condition_consec", "Age_cat_consec", "name", "estimate", "SE",
## "df", "t.ratio")

condition_consec	Age_cat_consec	name	estimate	SE	df	t.ratio	p.value.unadj	p.value.fdr
Placebo - Interferon	Y - O	Default	-2.1447154	1.648336	26	-1.3011397	0.2046237	0.4092474
Placebo - Interferon	Y - O	Dorsal.Attention	2.5701632	1.855185	26	1.3853946	0.1777031	0.4092474
Placebo - Interferon	Y - O	Frontoparietal	-1.3150997	2.319591	26	-0.5669533	0.5756085	0.6578382
Placebo - Interferon	Y - O	Limbic	-6.7478338	3.070461	26	-2.1976611	0.0370787	0.2966294
Placebo - Interferon	Y - O	Somatomotor	-0.6377358	1.707834	26	-0.3734180	0.7118649	0.7118649
Placebo - Interferon	Y - O	Subcortical	4.9273002	2.869787	26	1.7169568	0.0978788	0.3915153
Placebo - Interferon	Y - O	Ventral.Attention	1.5646368	2.320432	26	0.6742869	0.5060805	0.6578382
Placebo - Interferon	Y - O	Visual	1.2334637	1.282027	26	0.9621200	0.3448526	0.5517642

# 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
# 
# knitr::kable(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

## Joining, by = "fileIndex"

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

## boundary (singular) fit: see help('isSingular')

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

## boundary (singular) fit: see help('isSingular')

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

## boundary (singular) fit: see help('isSingular')

tab_model(fit1,fit2,fit3)

	degree			degree			degree
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	34.64	33.02 – 36.26	<0.001	34.55	32.66 – 36.45	<0.001	34.90	31.15 – 38.65	<0.001
condition [Placebo]	1.26	-1.03 – 3.55	0.282	1.26	-1.03 – 3.55	0.282	0.14	-5.16 – 5.45	0.958
quartile	21.29	20.70 – 21.88	<0.001	21.29	20.70 – 21.88	<0.001	21.16	19.78 – 22.53	<0.001
condition [Placebo] × quartile	0.03	-0.81 – 0.87	0.947	0.03	-0.81 – 0.87	0.947	0.46	-1.48 – 2.41	0.639
Age				0.00	-0.02 – 0.02	0.860	-0.01	-0.08 – 0.07	0.881
condition [Placebo] × Age							0.03	-0.08 – 0.13	0.648
quartile × Age							0.00	-0.02 – 0.03	0.833
(condition [Placebo] × quartile) × Age							-0.01	-0.05 – 0.03	0.626
Random Effects
σ2	2524.01			2524.06			2524.22
τ00	0.00 Subj_ID			0.00 Subj_ID			0.00 Subj_ID
N	28 Subj_ID			28 Subj_ID			28 Subj_ID
Observations	44128			44128			44128
Marginal R2 / Conditional R2	0.184 / NA			0.184 / NA			0.184 / NA

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)	71.53	70.99 – 72.07	<0.001	71.32	70.54 – 72.10	<0.001	71.80	70.54 – 73.07	<0.001
condition [Placebo]	3.22	2.52 – 3.92	<0.001	3.22	2.52 – 3.92	<0.001	1.92	0.29 – 3.54	0.021
quartile	26.01	25.83 – 26.19	<0.001	26.01	25.83 – 26.19	<0.001	25.87	25.45 – 26.29	<0.001
condition [Placebo] × quartile	-0.38	-0.64 – -0.13	0.003	-0.38	-0.64 – -0.13	0.003	0.04	-0.56 – 0.63	0.906
Age				0.00	-0.01 – 0.02	0.463	-0.01	-0.03 – 0.02	0.639
condition [Placebo] × Age							0.03	-0.00 – 0.06	0.081
quartile × Age							0.00	-0.01 – 0.01	0.461
(condition [Placebo] × quartile) × Age							-0.01	-0.02 – 0.00	0.124
Random Effects
σ2	118.07			118.07			118.06
τ00	0.37 Subj_ID			0.38 Subj_ID			0.38 Subj_ID
ICC	0.00			0.00			0.00
N	28 Subj_ID			28 Subj_ID			28 Subj_ID
Observations	22064			22064			22064
Marginal R2 / Conditional R2	0.876 / 0.876			0.876 / 0.876			0.876 / 0.876

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

Summary

Betweenness centrality

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

## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID'. You can override
## using the `.groups` argument.

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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID', 'condition'. You
## can override using the `.groups` argument.

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

## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.

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

## Joining, by = c("condition_consec", "name", "estimate", "SE", "df", "t.ratio")

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

## Joining, by = c("contrast", "name", "estimate", "SE", "df", "t.ratio")

# 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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID', 'condition',
## 'name'. You can override using the `.groups` argument.

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

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: Age_cat

DT::datatable(nice(main_anova_summary)) 

emmip(main_anova_summary, name~Age_cat | threshold)

## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 8. Consider
## specifying shapes manually if you must have them.

emmip(main_anova_summary, ~condition | Age_cat)

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

## Converting to factor: Age_cat

## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.

## Contrasts set to contr.sum for the following variables: 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))

## Joining, by = c("condition_consec", "Age_cat_consec", "name", "estimate", "SE",
## "df", "t.ratio")

# focus on threshold=3.5
main_anova_summary <- aov_ez("Subj_ID", "meanBWC", data.aggr,within =c("condition","name","threshold"), between = "Age_cat")

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: Age_cat

DT::datatable(nice(main_anova_summary))

EMM <- emmeans(main_anova_summary, ~  condition: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))

## Joining, by = c("condition_consec", "Age_cat_consec", "name", "estimate", "SE",
## "df", "t.ratio")

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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID', 'condition'. You
## can override using the `.groups` argument.

# regular anova
main_anova_summary_a <- aov_ez("Subj_ID", "meanBWC", data.aggr.age, between = "Age_cat")

## Converting to factor: Age_cat

## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.

## Contrasts set to contr.sum for the following variables: 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, 26 659.11 1.56 .057    .223
## ---
## 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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'Subj_ID', 'condition'. You can override
## using the `.groups` argument.

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

## Adding missing grouping variables: `Subj_ID`, `condition`

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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'Subj_ID', 'condition'. You can override
## using the `.groups` argument.

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

## Registered S3 method overwritten by 'ggside':
##   method from   
##   +.gg   ggplot2
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

## Converting to factor: Age_cat

## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.

## Contrasts set to contr.sum for the following variables: 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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID', 'condition',
## 'name'. You can override using the `.groups` argument.

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

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: Age_cat

DT::datatable(nice(main_anova_summary)) 

emmip(main_anova_summary, name~Age_cat | threshold)

## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 8. Consider
## specifying shapes manually if you must have them.

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

## Converting to factor: Age_cat

## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.

## Contrasts set to contr.sum for the following variables: 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))

## Joining, by = c("Age_cat_consec", "name", "estimate", "SE", "df", "t.ratio")

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

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: 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))

## Joining, by = c("Age_cat_consec", "name", "estimate", "SE", "df", "t.ratio")

Summary

Global efficiency

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

## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID'. You can override
## using the `.groups` argument.

# 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 = 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 = 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("network","condition","threshold"), between = "Age_cat")
DT::datatable(nice(main_anova_summary))

emmip(main_anova_summary, condition~Age_cat | network)

# anova averaging across threshold
main_anova_summary <- aov_ez("Subj_ID", "n_geff", data.aggr,within =c("network","condition"), between = "Age_cat")
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 = 3.5
main_anova_summary <- aov_ez("Subj_ID", "n_geff", data.aggr,within =c("condition","threshold","network"), between = "Age_cat")
EMM <- emmeans(main_anova_summary, ~condition:Age_cat | network, at=list(threshold="X3.5")) %>% contrast(interaction = c( "consec", "consec"))
DT::datatable(EMM %>% as.data.frame()) 

Summary

Local efficiency

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

## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID'. You can override
## using the `.groups` argument.

# 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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID', 'condition'. You
## can override using the `.groups` argument.

# regular anova
main_anova_summary <- 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"))

## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.

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

## Joining, by = c("condition_consec", "name", "estimate", "SE", "df", "t.ratio")

# 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

## Joining, by = "ROI"
## Joining, by = "fileIndex"
## `summarise()` has grouped output by 'threshold', 'Subj_ID', 'condition',
## 'name'. You can override using the `.groups` argument.

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

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: Age_cat

emmip(main_anova_summary, condition~Age_cat | name)

DT::datatable(nice(main_anova_summary))

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

## Converting to factor: Age_cat

## Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
## To turn off this warning, pass `fun_aggregate = mean` explicitly.

## Contrasts set to contr.sum for the following variables: Age_cat

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"),between = "Age_cat")

## Converting to factor: Age_cat
## Contrasts set to contr.sum for the following variables: Age_cat

EMM <- emmeans(main_anova_summary, ~condition: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))

## Joining, by = c("condition_consec", "Age_cat_consec", "name", "estimate", "SE",
## "df", "t.ratio")

# 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