PACt-MD Multishell GBSS Notebook

Motivation

Two alternative pathways (frontostriatal and medial-temporal) are thought to reflect the pathophysiology of major depressive disorder (MDD). The frontostriatal pathway, connecting the frontal lobes and basal ganglia, governs executive functions, including, selection and perception of important information, manipulation of information in working memory, planning and organization, behavioral control, adaptation to changes, and decision making, and is thought to be particularly vulnerable to cardiovascular risk factors. In contrast, the medial-temporal lobe circuit, which is responsible for binding information and memory, is primarily affected by HPA-axis glucocorticoids. Elevated levels of cortisol, for example, have been shown to shrink hippocampal tissues, and lead to lower memory performance. Importantly, depression has been linked to both …

Methods Overview

This set of analyses was inspired by work by Arash Nazeri who first proposed adapting the TBSS methodology for use with grey matter (https://github.com/arash-n/GBSS). The data were preprocessed using the following scripts: https://github.com/johnaeanderson/dwi_preprocessing. NODDI parameters were estimated using the MDT toolbox https://github.com/cbclab/MDT.

Figure 2. Methods overview. A) (adapted from ___) describes the NODDI model and derived parameter estimates (NDI, ODI, ISO). B) provides a brief visual overview of the preprocessing steps and GBSS (Gray-matter Based Spatial Statistics, NODDI-GBSS). Panel i) shows the raw DWI data, panel ii) shows denoised, EDDY/motion/susceptibility corrected data, panel iii) shows the NODDI model fISO estimate, and panel iv) shows the pseudo-T1 image estimated via the GBSS algorithm. C) shows the group-derived grey matter skeleton which was used to constrain all analyses.

Descriptive Statistics:

• include consort diagram here

#load in the PLS results 
#wm_CSF_Age_Mat <-("/Volumes/Seagate_Backup/stats/subset_with_controls/analysis_age_HcDepMciMciDep_CSF_wmskel_HC_STRUCTresult.mat")
#gm_CSF_Age_Mat <-("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/02_ControllingSexEducationNoise_Age_Covariate/Grey_Matter/CSF/analysis_CSF_age_ci83_STRUCTresult.mat")

#wm_CSF_Age_BS <- Beh_PLS_BrainScore(wm_CSF_Age_Mat,"CSF",1)
#gm_CSF_Age_BS <- Beh_PLS_BrainScore(gm_CSF_Age_Mat,"CSF",1)
descriptives <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/descriptives.xls")
Demo_Data_with_Reserve_Score <- read_csv("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/Demo_Data_with_Reserve_Score.csv", 
    col_types = cols(X1 = col_skip()))
age <- subset(Demo_Data_with_Reserve_Score, select = c("demo_id","age"))
colnames(age) <- c("subjects","age")

gm_CSF_Age_BS <- read_csv("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/gm_CSF_Age_BS.csv")
gm_CSF_Age_BS$subjects <- noquote(gm_CSF_Age_BS$subjects)

Extract demographic scores

Combined Parameter Analysis

source("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/functions.R")
descriptives <- read_excel("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/Centered_Scores_PLS_Order_6_groups.xlsx")
descriptives$subjects <- NULL

#read in "tissue types"
Tissue_Types <- read_csv("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/Tissue_Types.csv")
#Pull out the brain scores for the combined data
gm_CSF_ODI_NDI_Age_no_groups_Mat <- "/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/02_ControllingSexEducationNoise_Age_Covariate/Grey_Matter/Combined/Age_3_tissues_no_groups_STRUCTresult.mat"

gm_CSF_ODI_NDI_Age_no_groups_BS <- Beh_PLS_BrainScore(gm_CSF_ODI_NDI_Age_no_groups_Mat,"All",1)
gm_CSF_ODI_NDI_Age_no_groups_BS$Tissue <- Tissue_Types$Tissue
#invert LV1
gm_CSF_ODI_NDI_Age_no_groups_BS$LV1 <- gm_CSF_ODI_NDI_Age_no_groups_BS$LV1*-1

#Plot
p7 <- ggscatter(gm_CSF_ODI_NDI_Age_no_groups_BS, x = "Age", y = "LV1", color = "Tissue",palette = c("#00AFBB", "#E7B800", "#FC4E07") )+facet_wrap(~Tissue)+theme_cowplot()+geom_smooth(method ="rlm", se=FALSE, colour = "black")+ theme(legend.position = "none")#+stat_cor(label.x =60, label.y = 170)


#now plot the raincloudplots
a <- subset(gm_CSF_ODI_NDI_Age_no_groups_BS, select = c("Subject","Tissue","Age","LV1","LV2","LV3"))

## Permutation test similar to ANOVA
#permuted <- a %>%
#dplyr::mutate(Tissue = factor(Tissue)) %>%
#specify(LV1 ~ Tissue) %>%
#hypothesize(null = "independence") %>%
#generate(reps = 2000, type = "permute") %>%
#calculate(stat = "F")

b <- a %>%
  group_by(Tissue)%>%
  do(data.frame(val=bootstraps(., times = 2000)))
#checking to see what the data looks like for one resample...
#analysis(bt_data$splits[[1]]) %>% as_tibble()
#define the correlation function
get_cor <- function(val.splits) {
  # access to the sample data 
  split_data <- analysis(val.splits)
  # calculate the correlation value 
  split_cor<- cor(split_data$Age, split_data$LV1)
  return(split_cor)
}

b$bt_cors <- map_dbl(b$val.splits, get_cor)

#Confidence Intervals
confidence_intervals <- b %>%
  group_by(Tissue) %>%
  group_modify(~ {
     quantile(.x$bt_cors, probs = c(0.025, 0.5, 0.975)) %>%
     tibble::enframe(name = "prob", value = "quantile")
  })
confidence_intervals_wide <- dcast(data = confidence_intervals, formula = Tissue ~ prob, value.var = "quantile")

CI_Table <- ggtexttable(confidence_intervals_wide, rows = NULL, theme = ttheme("blank"))

#Rainclouds with boxplots
p8 <- ggplot(b,aes(x=Tissue,y=bt_cors, fill = Tissue, colour = Tissue))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2)+
  geom_point(position = position_jitter(width = .15), size = .25)+
#note that here we need to set the x-variable to a numeric variable and bump it to get the boxplots to line up with the rainclouds. 
  geom_boxplot(aes(x = as.numeric(as.factor(Tissue))+0.25, y = bt_cors),outlier.shape = NA, alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('Score')+xlab('Tissue')+theme_cowplot()+guides(fill = FALSE, colour = FALSE) +
  geom_errorbar(data = confidence_intervals_wide, aes(x = Tissue, y = `2.5%`, ymin = `2.5%`, ymax = `97.5%`), position = position_nudge(-0.25), colour = "BLACK", width = 0.2, size = 0.8)+ylab("Correlation")+scale_color_manual(values=c("#00AFBB", "#E7B800", "#FC4E07"))+ scale_fill_manual(values=c("#00AFBB", "#E7B800", "#FC4E07"))+geom_hline(yintercept =0, linetype=2)

p9 <- ggdraw() + draw_image("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/04_Results/01_Figures/GM_fISO_ODI_NDI_Age_combined_LV1.png")

ggarrange(p9, ggarrange(p7, p8, labels=c("b","c"), ncol = 2), nrow = 2, labels ="a")

Given that the NDI parameter neurite density, does not correlate with age within grey matter, we do not consider this parameter further.

fISO Analysis

# load in the data here...
source("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/functions.R")

gm_CSF_Age_Mat <- ("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/02_ControllingSexEducationNoise_Age_Covariate/Grey_Matter/CSF/HC_MDD_aMCI_NaMCI_aMCI_MDD_NaMCI_MDD_fISO_age_STRUCTresult.mat")

#then white matter and age data



descriptives <- read_excel("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/Centered_Scores_PLS_Order_6_groups.xlsx")
descriptives$subjects <- NULL

#Now creating the gm_age datamats
#Beh_PLS_Result <- function(x, LV, label, Invert)
gm_CSF_Age_Dat <-  Beh_PLS_Result(gm_CSF_Age_Mat,1,"fISO",1)[[4]]

bar_age_gm <- ggplot(gm_CSF_Age_Dat, aes(x = group, y = Correlation, fill= group))+geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=Correlation-lower_adj, ymax=Correlation + upper_adj), width=0,
                 position=position_dodge(.9), colour="black", size = 1)  + geom_hline(yintercept = 0, colour="black")+
theme(text = element_text(size=20)) + scale_color_manual(values=c("#56C1FF", "#00A2FF", "#0076BA","#004D7F"))+ scale_fill_manual(values=c("#56C1FF", "#00A2FF", "#0076BA","#004D7F","#0076BA","#004D7F")) + theme(legend.position = "none")#+ facet_wrap(~Label)
#plot(bar_age_gm)
#ggarrange(bar_age_gm, bar_age_wm, ncol=1, nrow=2, align="v", common.legend = TRUE, legend ="right", labels = c("grey matter", "white matter"))


################################################
# Pull out the brain scores and use vs the behavioural data
################################################

#Now creating the gm_age datamats

gm_CSF_Age_BS <- Beh_PLS_BrainScore(gm_CSF_Age_Mat,"CSF",1)
gm_CSF_Age_BS$Group <- factor(gm_CSF_Age_BS$Group, levels = c("HC","naMCI","MDD_naMCI", "MDD", "aMCI", "MDD_aMCI" ))


#gm_NDI_Age_BS <- Beh_PLS_BrainScore(gm_NDI_Age_Mat,"NDI",1)
#
#
#wm_age_BS<- rbind(wm_ODI_Age_BS,wm_NDI_Age_BS,wm_CSF_Age_BS)
#wm_age_BS$modality <- rep("wm", 828)
#gm_age_BS<- rbind(gm_ODI_Age_BS,gm_NDI_Age_BS,gm_CSF_Age_BS)
#gm_age_BS$modality <- rep("gm", 828)
#age_BS <- rbind(wm_age_BS, gm_age_BS)
#age_BS_Long <- melt(age_BS, measure.var =c("LV1","LV2","LV3"))

#Scatterplot of CSF LV1
CSF_Age_Scatter <- ggscatter(gm_CSF_Age_BS, x = "Age", y = "LV1", color = "Group",palette = c("#56C1FF", "#00A2FF", "#0076BA","#004D7F","#0076BA","#004D7F") )+facet_wrap(~Group)+theme_cowplot()+geom_smooth(method ="rlm", se=FALSE, aes(colour = Group))+ theme(legend.position = "none")


#p2 <- ggarrange(CSF_Age_Scatter, bar_age_gm,common.legend = TRUE)

p3 <- ggdraw() + draw_image("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/04_Results/01_Figures/GM_fISO_Age_6_groups_LV1.png")
#plot_grid(p3, CSF_Age_Scatter, bar_age_gm, labels = "AUTO", ncol=3, nrow = 1)

a <- subset(gm_CSF_Age_BS, select = c("Subject","Group","Age","LV1","LV2","LV3"))

library(dplyr)
b <- a %>%
  group_by(Group)%>%
  do(data.frame(val=bootstraps(., times = 2000)))
#checking to see what the data looks like for one resample...
#analysis(bt_data$splits[[1]]) %>% as_tibble()
#define the correlation function
get_cor <- function(val.splits) {
  # access to the sample data 
  split_data <- analysis(val.splits)
  # calculate the correlation value 
  split_cor<- cor(split_data$Age, split_data$LV1)
  return(split_cor)
}

b$bt_cors <- map_dbl(b$val.splits, get_cor)

#now plot the results

p4 <- ggscatter(gm_CSF_Age_BS, x = "Age", y = "LV1", color = "Group",add = "reg.line",  # Add regression line
   add.params = list(color = "black"), # Customize regression line
   conf.int = FALSE)+facet_wrap(~Group)+theme_cowplot()+ theme(legend.position = "none") + scale_color_manual(values = c("black","HC" = "#3F7F4F", "MDD" = "#775577", "naMCI" = "#A82A00", "MDD_naMCI" = "#116688","aMCI" = "#BFBF1F", "MDD_aMCI" = "#4C413D"))   


c <- summarySE(data=b, measurevar = "bt_cors", groupvars = "Group")

confidence_intervals <- b %>%
  group_by(Group) %>%
  group_modify(~ {
     quantile(.x$bt_cors, probs = c(0.025, 0.5, 0.975)) %>%
     tibble::enframe(name = "prob", value = "quantile")
  })
confidence_intervals_wide <- dcast(data = confidence_intervals, formula = Group ~ prob, value.var = "quantile")

#Rainclouds with boxplots
library(forcats)

p5 <- b %>%
  mutate(Group = fct_reorder(Group, bt_cors, .fun='mean')) %>%
 ggplot( aes(x=Group,y=bt_cors, fill = Group, colour = Group))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2)+
  geom_point(position = position_jitter(width = .15), size = .25)+
#note that here we need to set the x-variable to a numeric variable and bump it to get the boxplots to line up with the rainclouds. 
  geom_boxplot(aes(x = as.numeric(Group)+0.25, y = bt_cors),outlier.shape = NA, alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('Score')+xlab('Group')+theme_cowplot()+guides(fill = FALSE, colour = FALSE) +
  geom_errorbar(data = confidence_intervals_wide, aes(x = Group, y = `2.5%`, ymin = `2.5%`, ymax = `97.5%`), position = position_nudge(-0.25), colour = "BLACK", width = 0.2, size = 0.8)+ylab("Correlation")+ scale_color_manual(values = c("HC" = "#3F7F4F", "MDD" = "#775577", "naMCI" = "#A82A00", "MDD_naMCI" = "#116688","aMCI" = "#BFBF1F", "MDD_aMCI" = "#4C413D"))+ scale_fill_manual(values = c("HC" = "#3F7F4F", "MDD" = "#775577", "naMCI" = "#A82A00", "MDD_naMCI" = "#116688","aMCI" = "#BFBF1F", "MDD_aMCI" = "#4C413D")) + coord_flip()

#plot(p5)
#plot_grid(p3, CSF_Age_Scatter, p5, labels = "AUTO", ncol=3, nrow = 1)


#group differences (from bootstrap):

#      [,1]       [,2]       
# [1,] "HC"       "MDD"      
# [2,] "HC"       "aMCI"     
# [3,] "HC"       "naMCI"    
# [4,] "HC"       "MDD_aMCI" 
# [5,] "HC"       "MDD_naMCI"
# [6,] "MDD"      "aMCI"     
# [7,] "MDD"      "naMCI"    
# [8,] "MDD"      "MDD_aMCI" 
# [9,] "MDD"      "MDD_naMCI"
#[10,] "aMCI"     "naMCI"    
#[11,] "aMCI"     "MDD_aMCI" 
#[12,] "aMCI"     "MDD_naMCI"
#[13,] "naMCI"    "MDD_aMCI" 
#[14,] "naMCI"    "MDD_naMCI"
#[15,] "MDD_aMCI" "MDD_naMCI"

HC_sub_MDD <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="MDD")$bt_cors
HC_sub_aMCI <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="aMCI")$bt_cors
HC_sub_naMCI <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="naMCI")$bt_cors
HC_sub_MDD_aMCI <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="MDD_aMCI")$bt_cors
HC_sub_MDD_naMCI <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors
MDD_sub_aMCI <- subset(b, Group =="MDD")$bt_cors - subset(b, Group =="aMCI")$bt_cors
MDD_sub_naMCI <- subset(b, Group =="MDD")$bt_cors - subset(b, Group =="naMCI")$bt_cors
MDD_sub_MDD_aMCI <- subset(b, Group =="MDD")$bt_cors - subset(b, Group =="MDD_aMCI")$bt_cors
MDD_sub_MDD_naMCI <- subset(b, Group =="MDD")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors
aMCI_sub_naMCI <- subset(b, Group =="aMCI")$bt_cors - subset(b, Group =="naMCI")$bt_cors
aMCI_sub_MDD_aMCI <- subset(b, Group =="aMCI")$bt_cors - subset(b, Group =="MDD_aMCI")$bt_cors
aMCI_sub_MDD_naMCI <- subset(b, Group =="aMCI")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors
naMCI_sub_MDD_aMCI <- subset(b, Group =="naMCI")$bt_cors - subset(b, Group =="MDD_aMCI")$bt_cors
naMCI_sub_MDD_naMCI <- subset(b, Group =="naMCI")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors
MDD_aMCI_sub_MDD_naMCI <- subset(b, Group =="MDD_aMCI")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors

all_comparisons <- cbind(HC_sub_MDD ,HC_sub_aMCI ,HC_sub_naMCI ,HC_sub_MDD_aMCI ,HC_sub_MDD_naMCI ,MDD_sub_aMCI ,MDD_sub_naMCI ,MDD_sub_MDD_aMCI ,MDD_sub_MDD_naMCI ,aMCI_sub_naMCI, aMCI_sub_MDD_aMCI ,aMCI_sub_MDD_naMCI ,naMCI_sub_MDD_aMCI ,naMCI_sub_MDD_naMCI, MDD_aMCI_sub_MDD_naMCI)
#quantile(all_comparisons, probs = c(0.025, 0.5, 0.975))
library(matrixStats)
probs <- c(0.025, 0.5, 0.975)

#robust cohen's d
elimna <-
function(m){
#
# remove any rows of data having missing values
#
if(is.null(dim(m)))m<-as.matrix(m)
ikeep<-c(1:nrow(m))
for(i in 1:nrow(m))if(sum(is.na(m[i,])>=1))ikeep[i]<-0
elimna<-m[ikeep[ikeep>=1],]
elimna
}

pbvar <-
function(x,beta=.2){
#   Compute the percentage bend midvariance
#
#   beta is the bending constant for omega sub N.
#
pbvar=0
x=elimna(x)
w<-abs(x-median(x))
w<-sort(w)
m<-floor((1-beta)*length(x)+.5)
omega<-w[m]
if(omega>0){
y<-(x-median(x))/omega
z<-ifelse(y>1,1,y)
z<-ifelse(z<(-1),-1,z)
pbvar<-length(x)*omega^2*sum(z^2)/(length(x[abs(y)<1]))^2
}
pbvar
}


cohensd <- 
  function(x){
    d <- round(mean(x)/sqrt(pbvar(x, 0.2)),digits = 2)
d
}
# Row quantiles
q <- rowQuantiles(t(all_comparisons), probs = probs)

r <- apply(t(all_comparisons), 1, cohensd)
r <- as.data.frame(r)
colnames(r) <- "d"
combo <- cbind(q, r)

#my_comparisons=list(c("HC","MDD"), c("HC","aMCI"), c("HC","naMCI"), c("HC","MDD_aMCI"), c("HC","MDD_naMCI"), #c("naMCI","MDD_naMCI"),c("naMCI","MDD"),c("naMCI","aMCI"), ,c("naMCI","aMCI"), ,c("naMCI","MDD_aMCI"), ) 
#stat.test <- compare_means(bt_cors ~ Group, data = b,method = "t.test")%>%
#  mutate(y.position = c(1.6, 1.8, 2,1.3,1.5,1.1))%>%
#  mutate(c_d = r)
#stat.test
##  
#
#
#p6 <- p5+
#  stat_pvalue_manual(
#    data = stat.test, label = "c_d",
#    xmin = "group1", xmax = "group2",
#    y.position = "y.position",
#    inherit.aes = FALSE
#    )
#
##    y.position = "y.position",
#table_comparisons <- ggtexttable(round(combo,digits = 2),rows = rownames(combo))
# 
plot_grid(p3, plot_grid(p4, p5, labels = c("b", "c"), ncol=2), nrow = 2, labels = "a")

kable(combo) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")

#use this to fortify the ggplot object

These results suggest a difference of kind for the fISO/age relationship between healthier individuals (rMDD, and HC) and those with more severe diagnoses. This fits with the idea that there is an observable break in the normal process of aging.

ODI Analysis

source("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/functions.R")


# load in the data here...

# Then grey matter and age data
gm_ODI_Age_Mat <-("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/02_ControllingSexEducationNoise_Age_Covariate/Grey_Matter/ODI/analysis_Grey_ODI_age_ci83_Amnestic_STRUCTresult.mat")

#then white matter and age data



descriptives <- read_excel("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/Centered_Scores_PLS_Order_6_groups.xlsx")
descriptives$subjects <- NULL



#Now creating the gm_age datamats
#Beh_PLS_Result <- function(x, LV, label, Invert)
gm_ODI_Age_Dat <- Beh_PLS_Result(gm_ODI_Age_Mat,1,"ODI",1)[[4]]

bar_age_gm <- ggplot(gm_ODI_Age_Dat, aes(x = group, y = Correlation, fill= group))+geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=Correlation-lower_adj, ymax=Correlation + upper_adj), width=0,
                 position=position_dodge(.9), colour="black", size = 1)  + geom_hline(yintercept = 0, colour="black")+
theme(text = element_text(size=20)) + scale_color_manual(values=c("#56C1FF", "#00A2FF", "#0076BA","#004D7F"))+ scale_fill_manual(values=c("#56C1FF", "#00A2FF", "#0076BA","#004D7F","#0076BA","#004D7F")) + theme(legend.position = "none")#+ facet_wrap(~Label)
#plot(bar_age_gm)
#ggarrange(bar_age_gm, bar_age_wm, ncol=1, nrow=2, align="v", common.legend = TRUE, legend ="right", labels = c("grey matter", "white matter"))


################################################
# Pull out the brain scores and use vs the behavioural data
################################################


#Now creating the gm_age datamats

gm_ODI_Age_BS <- Beh_PLS_BrainScore(gm_ODI_Age_Mat,"ODI",1)
gm_ODI_Age_BS$Group <- factor(gm_ODI_Age_BS$Group, levels = c("HC","naMCI","MDD_naMCI", "MDD", "aMCI", "MDD_aMCI" ))



#gm_NDI_Age_BS <- Beh_PLS_BrainScore(gm_NDI_Age_Mat,"NDI",1)
#
#
#wm_age_BS<- rbind(wm_ODI_Age_BS,wm_NDI_Age_BS,wm_CSF_Age_BS)
#wm_age_BS$modality <- rep("wm", 828)
#gm_age_BS<- rbind(gm_ODI_Age_BS,gm_NDI_Age_BS,gm_CSF_Age_BS)
#gm_age_BS$modality <- rep("gm", 828)
#age_BS <- rbind(wm_age_BS, gm_age_BS)
#age_BS_Long <- melt(age_BS, measure.var =c("LV1","LV2","LV3"))

#Scatterplot of CSF LV1
#CSF_Age_Scatter <- ggscatter(gm_ODI_Age_BS, x = "Age", y = "LV1", color = "Group",palette = c("#56C1FF", "#00A2FF", "#0076BA","#004D7F","#0076BA","#004D7F") )+facet_wrap(~Group)+theme_cowplot()+geom_smooth(method ="rlm", se=FALSE, aes(colour = Group))+ theme(legend.position = "none")


#p2 <- ggarrange(CSF_Age_Scatter, bar_age_gm,common.legend = TRUE)

p3 <- ggdraw() + draw_image("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/04_Results/01_Figures/GM_ODI_Age_6_groups_LV1.png")
#plot_grid(p3, CSF_Age_Scatter, bar_age_gm, labels = "AUTO", ncol=3, nrow = 1)

a <- subset(gm_ODI_Age_BS, select = c("Subject","Group","Age","LV1","LV2","LV3"))

library(dplyr)
b <- a %>%
  group_by(Group)%>%
  do(data.frame(val=bootstraps(., times = 2000)))
#checking to see what the data looks like for one resample...
#analysis(bt_data$splits[[1]]) %>% as_tibble()
#define the correlation function
get_cor <- function(val.splits) {
  # access to the sample data 
  split_data <- analysis(val.splits)
  # calculate the correlation value 
  split_cor<- cor(split_data$Age, split_data$LV1)
  return(split_cor)
}

b$bt_cors <- map_dbl(b$val.splits, get_cor)

b <- b %>%
  mutate(GroupNew = fct_reorder(Group, bt_cors, .fun='mean'))
#now plot the results

gm_ODI_Age_BS$Group <- factor(gm_ODI_Age_BS$Group, ordered = TRUE, 
                                levels = levels(b$GroupNew))


p4 <- ggscatter(gm_ODI_Age_BS, x = "Age", y = "LV1", color = "Group",add = "reg.line",  # Add regression line
   add.params = list(color = "black"), # Customize regression line
   conf.int = FALSE)+facet_wrap(~Group)+theme_cowplot()+ theme(legend.position = "none") + scale_color_manual(values = c("black","HC" = "#3F7F4F", "MDD" = "#775577", "naMCI" = "#A82A00", "MDD_naMCI" = "#116688","aMCI" = "#BFBF1F", "MDD_aMCI" = "#4C413D"))  

c <- summarySE(data=b, measurevar = "bt_cors", groupvars = "GroupNew")

confidence_intervals <- b %>%
  group_by(Group) %>%
  group_modify(~ {
     quantile(.x$bt_cors, probs = c(0.025, 0.5, 0.975)) %>%
     tibble::enframe(name = "prob", value = "quantile")
  })
confidence_intervals_wide <- dcast(data = confidence_intervals, formula = Group ~ prob, value.var = "quantile")

#Rainclouds with boxplots
library(forcats)

p5 <- b %>%
  ungroup(Group) %>%
  mutate(Group = fct_reorder(Group, bt_cors, .fun='mean')) %>%
 ggplot( aes(x=Group,y=bt_cors, fill = Group, colour = Group))+
  geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2)+
  geom_point(position = position_jitter(width = .15), size = .25)+
#note that here we need to set the x-variable to a numeric variable and bump it to get the boxplots to line up with the rainclouds. 
  geom_boxplot(aes(x = as.numeric(Group)+0.25, y = bt_cors),outlier.shape = NA, alpha = 0.3, width = .1, colour = "BLACK") +
  ylab('Score')+xlab('Group')+theme_cowplot()+guides(fill = FALSE, colour = FALSE) +
  geom_errorbar(data = confidence_intervals_wide, aes(x = Group, y = `2.5%`, ymin = `2.5%`, ymax = `97.5%`), position = position_nudge(-0.25), colour = "BLACK", width = 0.2, size = 0.8)+ylab("Correlation")+ scale_color_manual(values = c("HC" = "#3F7F4F", "MDD" = "#775577", "naMCI" = "#A82A00", "MDD_naMCI" = "#116688","aMCI" = "#BFBF1F", "MDD_aMCI" = "#4C413D"))+ scale_fill_manual(values = c("HC" = "#3F7F4F", "MDD" = "#775577", "naMCI" = "#A82A00", "MDD_naMCI" = "#116688","aMCI" = "#BFBF1F", "MDD_aMCI" = "#4C413D")) + coord_flip() 

#plot(p5)
#plot_grid(p3, CSF_Age_Scatter, p5, labels = "AUTO", ncol=3, nrow = 1)


#group differences (from bootstrap):
#      [,1]       [,2]       
# [1,] "HC"       "MDD"      
# [2,] "HC"       "aMCI"     
# [3,] "HC"       "naMCI"    
# [4,] "HC"       "MDD_aMCI" 
# [5,] "HC"       "MDD_naMCI"
# [6,] "MDD"      "aMCI"     
# [7,] "MDD"      "naMCI"    
# [8,] "MDD"      "MDD_aMCI" 
# [9,] "MDD"      "MDD_naMCI"
#[10,] "aMCI"     "naMCI"    
#[11,] "aMCI"     "MDD_aMCI" 
#[12,] "aMCI"     "MDD_naMCI"
#[13,] "naMCI"    "MDD_aMCI" 
#[14,] "naMCI"    "MDD_naMCI"
#[15,] "MDD_aMCI" "MDD_naMCI"


HC_sub_MDD <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="MDD")$bt_cors
HC_sub_aMCI <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="aMCI")$bt_cors
HC_sub_naMCI <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="naMCI")$bt_cors
HC_sub_MDD_aMCI <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="MDD_aMCI")$bt_cors
HC_sub_MDD_naMCI <- subset(b, Group =="HC")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors
MDD_sub_aMCI <- subset(b, Group =="MDD")$bt_cors - subset(b, Group =="aMCI")$bt_cors
MDD_sub_naMCI <- subset(b, Group =="MDD")$bt_cors - subset(b, Group =="naMCI")$bt_cors
MDD_sub_MDD_aMCI <- subset(b, Group =="MDD")$bt_cors - subset(b, Group =="MDD_aMCI")$bt_cors
MDD_sub_MDD_naMCI <- subset(b, Group =="MDD")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors
aMCI_sub_naMCI <- subset(b, Group =="aMCI")$bt_cors - subset(b, Group =="naMCI")$bt_cors
aMCI_sub_MDD_aMCI <- subset(b, Group =="aMCI")$bt_cors - subset(b, Group =="MDD_aMCI")$bt_cors
aMCI_sub_MDD_naMCI <- subset(b, Group =="aMCI")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors
naMCI_sub_MDD_aMCI <- subset(b, Group =="naMCI")$bt_cors - subset(b, Group =="MDD_aMCI")$bt_cors
naMCI_sub_MDD_naMCI <- subset(b, Group =="naMCI")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors
MDD_aMCI_sub_MDD_naMCI <- subset(b, Group =="MDD_aMCI")$bt_cors - subset(b, Group =="MDD_naMCI")$bt_cors

all_comparisons <- cbind(HC_sub_MDD ,HC_sub_aMCI ,HC_sub_naMCI ,HC_sub_MDD_aMCI ,HC_sub_MDD_naMCI ,MDD_sub_aMCI ,MDD_sub_naMCI ,MDD_sub_MDD_aMCI ,MDD_sub_MDD_naMCI ,aMCI_sub_naMCI, aMCI_sub_MDD_aMCI ,aMCI_sub_MDD_naMCI ,naMCI_sub_MDD_aMCI ,naMCI_sub_MDD_naMCI, MDD_aMCI_sub_MDD_naMCI)
#quantile(all_comparisons, probs = c(0.025, 0.5, 0.975))
library(matrixStats)
probs <- c(0.025, 0.5, 0.975)

#robust cohen's d
elimna <-
function(m){
#
# remove any rows of data having missing values
#
if(is.null(dim(m)))m<-as.matrix(m)
ikeep<-c(1:nrow(m))
for(i in 1:nrow(m))if(sum(is.na(m[i,])>=1))ikeep[i]<-0
elimna<-m[ikeep[ikeep>=1],]
elimna
}

pbvar <-
function(x,beta=.2){
#   Compute the percentage bend midvariance
#
#   beta is the bending constant for omega sub N.
#
pbvar=0
x=elimna(x)
w<-abs(x-median(x))
w<-sort(w)
m<-floor((1-beta)*length(x)+.5)
omega<-w[m]
if(omega>0){
y<-(x-median(x))/omega
z<-ifelse(y>1,1,y)
z<-ifelse(z<(-1),-1,z)
pbvar<-length(x)*omega^2*sum(z^2)/(length(x[abs(y)<1]))^2
}
pbvar
}


cohensd <- 
  function(x){
    d <- round(mean(x)/sqrt(pbvar(x, 0.2)),digits = 2)
d
}
# Row quantiles
q <- rowQuantiles(t(all_comparisons), probs = probs)
#print(q)

r <- apply(t(all_comparisons), 1, cohensd)
#print(r)
r <- as.data.frame(r)
colnames(r) <- "d"
combo <- cbind(q, r)
#use this to fortify the ggplot object
#use this to fortify the ggplot object

#my_comparisons=list(c("HC","MDD"), c("HC","aMCI"), c("HC","naMCI"), c("HC","MDD_aMCI"), c("HC","MDD_naMCI"), #c("naMCI","MDD_naMCI"),c("naMCI","MDD"),c("naMCI","aMCI"), ,c("naMCI","aMCI"), ,c("naMCI","MDD_aMCI"), ) 
#stat.test <- compare_means(bt_cors ~ Group, data = b,method = "t.test")%>%
#  mutate(y.position = c(1.6, 1.8, 2,1.3,1.5,1.1))%>%
#  mutate(c_d = r)
#stat.test
##  
#
#
#p6 <- p5+
#  stat_pvalue_manual(
#    data = stat.test, label = "c_d",
#    xmin = "group1", xmax = "group2",
#    y.position = "y.position",
#    inherit.aes = FALSE
#    )
#
##    y.position = "y.position",
#table_comparisons <- ggtexttable(round(combo,digits = 2),rows = rownames(combo))
# 
plot_grid(p3, plot_grid(p4, p5, labels = c("b", "c"), ncol=2), nrow = 2, labels = "a")

kable(combo) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")

##In White Matter, now to look at the tissue types by age collapsing across groups:

# #function to pull data from Mean Behavioural PLS outputs...
# Beh_PLS_BrainScore <- function(x, label, Invert){
# input <- readMat(x)
# res <- input$result
# age <- as.data.frame(res[[10]])
# brain_score <- as.data.frame(res[[8]][,1:3])
# subjects <- as.data.frame(unlist(input$subj.name.lst))
# colnames(subjects) <- "subjects"
# Label <- rep(label, length(brain_score$V1))
# #sub_pre <- rep("PAC01_CMH_", length(brain_score))
# #sub_post <- rep("01",length(brain_score))
# #cols <- c("sub_pre", "subjects", "sub_post")
# #data <- cbind(sub_pre, subjects, sub_post)
# 
# #Sub_ID <- do.call(paste, c(data[cols], sep="")
# 
# bs <- cbind(subjects, Label,descriptives, brain_score)
# colnames(bs) <- c("subjects","label","Subject","Group","Sexc","Agec","Educationc","PC1c","Age","Sex","Education","PC1", "LV1","LV2","LV3")
# return(bs)
# }
# 
# #read in "tissue types"
# Tissue_Types <- read_csv("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/Tissue_Types.csv")
# #Pull out the brain scores for the combined data
# wm_CSF_ODI_NDI_Age_no_groups_BS <- Beh_PLS_BrainScore(wm_CSF_ODI_NDI_Age_no_groups_Mat,"All",1)
# wm_CSF_ODI_NDI_Age_no_groups_BS$Tissue <- Tissue_Types$Tissue
# #invert LV1
# #wm_CSF_ODI_NDI_Age_no_groups_Mat$LV1 <- wm_CSF_ODI_NDI_Age_no_groups_Mat$LV1*-1
# 
# a <- subset(wm_CSF_ODI_NDI_Age_no_groups_BS, select = c("Subject","Tissue","Age","LV1","LV2","LV3"))
# 
# a_long <- melt(data= a, id.vars = c("Subject","Tissue","Age"))
# 
# #Plot
# p7 <- ggscatter(a_long, x = "Age", y = "value", color = "Tissue",palette = c("#00AFBB", "#E7B800", "#FC4E07") )+facet_grid(variable~Tissue, scales="free_y")+theme_cowplot()+geom_smooth(method ="rlm", se=FALSE, colour = "black")+ theme(legend.position = "none")#+stat_cor(label.x =60, label.y = 170)
# 
# #now plot the raincloudplots
# 
# ## Permutation test similar to ANOVA
# #permuted <- a %>%
# #dplyr::mutate(Tissue = factor(Tissue)) %>%
# #specify(LV1 ~ Tissue) %>%
# #hypothesize(null = "independence") %>%
# #generate(reps = 2000, type = "permute") %>%
# #calculate(stat = "F")
# 
# b <- a_long %>%
#   group_by(Tissue, variable)%>%
#   do(data.frame(val=bootstraps(., times = 2000)))
# #checking to see what the data looks like for one resample...
# #analysis(bt_data$splits[[1]]) %>% as_tibble()
# #define the correlation function
# get_cor <- function(val.splits) {
#   # access to the sample data 
#   split_data <- analysis(val.splits)
#   # calculate the correlation value 
#   split_cor<- cor(split_data$Age, split_data$value)
#   return(split_cor)
# }
# 
# b$bt_cors <- map_dbl(b$val.splits, get_cor)
# 
# #Confidence Intervals
# confidence_intervals <- b %>%
#   group_by(Tissue, variable) %>%
#   group_modify(~ {
#      quantile(.x$bt_cors, probs = c(0.025, 0.5, 0.975)) %>%
#      tibble::enframe(name = "prob", value = "quantile")
#   })
# confidence_intervals_wide <- dcast(data = confidence_intervals, formula = Tissue + variable ~ prob, value.var = "quantile")
# 
# CI_Table <- ggtexttable(confidence_intervals_wide, rows = NULL, theme = ttheme("blank"))
# 
# #Rainclouds with boxplots
# p8 <- ggplot(b,aes(x=Tissue,y=bt_cors, fill = Tissue, colour = Tissue))+
#   geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2)+
#   geom_point(position = position_jitter(width = .15), size = .25)+
# #note that here we need to set the x-variable to a numeric variable and bump it to get the boxplots to line up with the rainclouds. 
#   geom_boxplot(aes(x = as.numeric(as.factor(Tissue))+0.25, y = bt_cors),outlier.shape = NA, alpha = 0.3, width = .1, colour = "BLACK") +
#   ylab('Score')+xlab('Tissue')+theme_cowplot()+guides(fill = FALSE, colour = FALSE) +
#   geom_errorbar(data = confidence_intervals_wide, aes(x = Tissue, y = `2.5%`, ymin = `2.5%`, ymax = `97.5%`), position = position_nudge(-0.25), colour = "BLACK", width = 0.2, size = 0.8)+ylab("Correlation")+scale_color_manual(values=c("#00AFBB", "#E7B800", "#FC4E07"))+ scale_fill_manual(values=c("#00AFBB", "#E7B800", "#FC4E07"))+geom_hline(yintercept =0, linetype=2)+facet_grid(rows = vars(variable)) 
# 
# plot(p8)
# 
# ggarrange(p7, p8, labels=c("b","c"))
# ```
# ##Task PLS looking at difference in CSF NDI and ODI across groups
# #First in Grey matter
# 
# ```{r}
# gm_All_Mat <-("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/01_ControllingAgeSexEducationNoise/Grey_Matter/Combined/All_Grey_no_age_all_tissues_group_comparison_STRUCTresult.mat")
# 
# #function to pull data from Mean Behavioural PLS outputs...
# MCPLS_BrainScore <- function(x, label, Invert){
# input <- readMat(x)
# res <- input$result
# brain_score <- as.data.frame(res[[6]][,1])
# subjects <- as.data.frame(unlist(input$subj.name.lst))
# colnames(subjects) <- "subjects"
# n <- 3
# 
# HC <- subset(descriptives, Group =="HC")
# HC2 <- do.call("rbind", replicate(n, HC, simplify = FALSE))
# MDD <- subset(descriptives, Group =="MDD")
# MDD2 <- do.call("rbind", replicate(n, MDD, simplify = FALSE))
# MCI <- subset(descriptives, Group =="MCI")
# MCI2 <- do.call("rbind", replicate(n, MCI, simplify = FALSE))
# MDD_MCI <- subset(descriptives, Group =="MDD_MCI")
# MDD_MCI2 <- do.call("rbind", replicate(n, MDD_MCI, simplify = FALSE))
# 
# Descriptives2 <- rbind(HC2, MDD2, MCI2, MDD_MCI2)
# 
# 
# #sub_pre <- rep("PAC01_CMH_", length(brain_score))
# #sub_post <- rep("01",length(brain_score))
# #cols <- c("sub_pre", "subjects", "sub_post")
# #data <- cbind(sub_pre, subjects, sub_post)
# 
# #Sub_ID <- do.call(paste, c(data[cols], sep="")
# 
# bs <- cbind(brain_score, Descriptives2)
# colnames(bs) <- c("LV1", "subjects","Group","Sexc","Agec","Educationc","PC1c","Age","Sex","Education","PC1")
# return(bs)
# }
# 
# Grey_MCPLS_All_Tissues <- MCPLS_BrainScore(gm_All_Mat,"all",1)
# #read in "tissue types"
# Tissue_Types_by_Group <- read_csv("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/Tissue_Types_by_Group.csv")
# #Pull out the brain scores for the combined data
# Grey_MCPLS_All_Tissues$Tissue <- Tissue_Types_by_Group$x
# 
# 
# a <- subset(Grey_MCPLS_All_Tissues, select = c("subjects","Tissue", "Group","Age","LV1"))
# 
# Grand_Mean <- mean(a$LV1)
# a$LV1 <- a$LV1 - Grand_Mean
# 
# b <- a %>%
#   group_by(Tissue, Group)%>%
#   do(data.frame(val=bootstraps(., times = 2000)))
# #checking to see what the data looks like for one resample...
# #analysis(bt_data$splits[[1]]) %>% as_tibble()
# #define the correlation function
# get_mean <- function(val.splits) {
#   # access to the sample data 
#   split_data <- analysis(val.splits)
#   # calculate the correlation value 
#   split_mean<- mean(split_data$LV1)
#   return(split_mean)
# }
# 
# b$split_mean <- map_dbl(b$val.splits, get_mean)
# 
# #Confidence Intervals
# confidence_intervals <- b %>%
#   group_by(Tissue, Group) %>%
#   group_modify(~ {
#      quantile(.x$split_mean, probs = c(0.025, 0.5, 0.975)) %>%
#      tibble::enframe(name = "prob", value = "quantile")
#   })
# confidence_intervals_wide <- dcast(data = confidence_intervals, formula = Tissue + Group ~ prob, value.var = "quantile")
# 
# CI_Table <- ggtexttable(confidence_intervals_wide, rows = NULL, theme = ttheme("blank"))
# 
# #Rainclouds with boxplots
# p8 <- ggplot(b,aes(x=Group,y=split_mean, fill = Group, colour = Group))+
#   geom_flat_violin(position = position_nudge(x = .25, y = 0),adjust =2)+
#   geom_point(position = position_jitter(width = .15), size = .25)+
# #note that here we need to set the x-variable to a numeric variable and bump it to get the boxplots to line up with the rainclouds. 
#   geom_boxplot(aes(x = as.numeric(as.factor(Group))+0.25, y = split_mean),outlier.shape = NA, alpha = 0.3, width = .1, colour = "BLACK") +
#   ylab('Score')+xlab('Group')+theme_cowplot()+guides(fill = FALSE, colour = FALSE) +
#   geom_errorbar(data = confidence_intervals_wide, aes(x = Group, y = `2.5%`, ymin = `2.5%`, ymax = `97.5%`), position = position_nudge(-0.25), colour = "BLACK", width = 0.2, size = 0.8)+ylab("Brain Score")+scale_color_manual(values=c("#56C1FF", "#00A2FF", "#0076BA","#004D7F"))+ scale_fill_manual(values=c("#56C1FF", "#00A2FF", "#0076BA","#004D7F"))+geom_hline(yintercept =0, linetype=2)+facet_grid(~Tissue, scales="free") 
# plot(p8)
# 
# #group differences (from bootstrap):
# HC_sub_MDD_CSF <- subset(b, Group =="HC" & Tissue =="CSF")$split_mean - subset(b, Group =="MDD" & Tissue =="CSF")$split_mean
# HC_sub_MCI_CSF <- subset(b, Group =="HC" & Tissue =="CSF")$split_mean - subset(b, Group =="MCI" & Tissue =="CSF")$split_mean
# HC_sub_MDD_MCI_CSF <- subset(b, Group =="HC" & Tissue =="CSF")$split_mean - subset(b, Group =="MDD_MCI" & Tissue =="CSF")$split_mean
# MDD_sub_MCI_CSF <- subset(b, Group =="MDD" & Tissue =="CSF")$split_mean - subset(b, Group =="MCI" & Tissue =="CSF")$split_mean
# MDD_sub_MDD_MCI_CSF <- subset(b, Group =="MDD" & Tissue =="CSF")$split_mean - subset(b, Group =="MDD_MCI" & Tissue =="CSF")$split_mean
# MCI_sub_MDD_MCI_CSF <- subset(b, Group =="MCI" & Tissue =="CSF")$split_mean - subset(b, Group =="MDD_MCI" & Tissue =="CSF")$split_mean
# 
# HC_sub_MDD_NDI <- subset(b, Group =="HC" & Tissue =="NDI")$split_mean - subset(b, Group =="MDD" & Tissue =="NDI")$split_mean
# HC_sub_MCI_NDI <- subset(b, Group =="HC" & Tissue =="NDI")$split_mean - subset(b, Group =="MCI" & Tissue =="NDI")$split_mean
# HC_sub_MDD_MCI_NDI <- subset(b, Group =="HC" & Tissue =="NDI")$split_mean - subset(b, Group =="MDD_MCI" & Tissue =="NDI")$split_mean
# MDD_sub_MCI_NDI <- subset(b, Group =="MDD" & Tissue =="NDI")$split_mean - subset(b, Group =="MCI" & Tissue =="NDI")$split_mean
# MDD_sub_MDD_MCI_NDI <- subset(b, Group =="MDD" & Tissue =="NDI")$split_mean - subset(b, Group =="MDD_MCI" & Tissue =="NDI")$split_mean
# MCI_sub_MDD_MCI_NDI <- subset(b, Group =="MCI" & Tissue =="NDI")$split_mean - subset(b, Group =="MDD_MCI" & Tissue =="NDI")$split_mean
# 
# HC_sub_MDD_ODI <- subset(b, Group =="HC" & Tissue =="ODI")$split_mean - subset(b, Group =="MDD" & Tissue =="ODI")$split_mean
# HC_sub_MCI_ODI <- subset(b, Group =="HC" & Tissue =="ODI")$split_mean - subset(b, Group =="MCI" & Tissue =="ODI")$split_mean
# HC_sub_MDD_MCI_ODI <- subset(b, Group =="HC" & Tissue =="ODI")$split_mean - subset(b, Group =="MDD_MCI" & Tissue =="ODI")$split_mean
# MDD_sub_MCI_ODI <- subset(b, Group =="MDD" & Tissue =="ODI")$split_mean - subset(b, Group =="MCI" & Tissue =="ODI")$split_mean
# MDD_sub_MDD_MCI_ODI <- subset(b, Group =="MDD" & Tissue =="ODI")$split_mean - subset(b, Group =="MDD_MCI" & Tissue =="ODI")$split_mean
# MCI_sub_MDD_MCI_ODI <- subset(b, Group =="MCI" & Tissue =="ODI")$split_mean - subset(b, Group =="MDD_MCI" & Tissue =="ODI")$split_mean
# 
# all_comparisons <- cbind(HC_sub_MDD_CSF,HC_sub_MCI_CSF,HC_sub_MDD_MCI_CSF,MDD_sub_MCI_CSF,MDD_sub_MDD_MCI_CSF,MCI_sub_MDD_MCI_CSF,HC_sub_MDD_ODI,HC_sub_MCI_ODI,HC_sub_MDD_MCI_ODI,MDD_sub_MCI_ODI,MDD_sub_MDD_MCI_ODI,MCI_sub_MDD_MCI_ODI,HC_sub_MDD_NDI,HC_sub_MCI_NDI,HC_sub_MDD_MCI_NDI,MDD_sub_MCI_NDI,MDD_sub_MDD_MCI_NDI,MCI_sub_MDD_MCI_NDI)
# 
# quantile(all_comparisons, probs = c(0.025, 0.5, 0.975))
# library(matrixStats)
# probs <- c(0.025, 0.5, 0.975)
# 
# q <- rowQuantiles(t(all_comparisons), probs = probs)
# print(q)
# 
# cohensd <- 
#   function(x){
#     d <- round(mean(x)/sqrt(pbvar(x, 0.2)),digits = 2)
# d
# }
# # Row quantiles
# 
# r <- apply(t(all_comparisons), 1, cohensd)
# print(r)
# 
# #use this to fortify the ggplot object
# 
# my_comparisons=list(c("HC","MDD"), c("HC","MCI"), c("HC","MDD_MCI"), c("MDD","MCI"), c("MDD","MDD_MCI"), c("MCI","MDD_MCI")) 
# stat.test <- compare_means(split_mean ~ Group, data = b,method = "t.test")%>%
#   mutate(y.position = c(1.6, 1.8, 2,1.3,1.5,1.1))%>%
#   mutate(c_d = r)
# stat.test
# #  

Whereas the fIS/age relationship demarcated a clear break, or difference of kind between healthier and sicker individuals, the relationship between ODI and age presents a more graduated change, thus suggesting a difference of degree, across diagnoses which roughly separates individuals with amnestic MCI from those who are non-amnestic.

PLS-PM brain-behavior analysis

This analysis is designed to allow depression, measured with the MADRS, and health factors, measured with BMI, the Framingham Index and the CIRS-G, to impact both Multishell measures and an emergent cognition score. The hope is that by including the multishell measures in the model, the direct path between health|depression and cognition is weakened (i.e., health’s effect on the brain manifests in measures of fISO and ODI).

Figure 3. PLS-PM analysis using a single “g” factor of cognition

Figure 4. PLS-PM analysis using emergent memory and EF factors

source("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/functions.R")

# laod health variables
#Framingham 
framingham_HC <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 Clinical Control.xlsx", sheet = "Framingham")
framingham_HC$subjects <- noquote(framingham_HC$fr_id)
framingham_Clinical <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 Clinical Case.xlsx", sheet = "Framingham")
  framingham_Clinical$subjects <- noquote(framingham_Clinical$fr_id)
  framingham <- rbind(framingham_Clinical,framingham_HC)
  framingham <- subset(framingham, select=c(fr_id, fr_smoker, fr_total_chol, fr_hdl_chol, fr_diabetes))
  colnames(framingham) <- c("subjects","smoker","total_chol","hdl","diabetes")
  framingham$smoker[framingham$smoker == 2 ] <- 0
  framingham$diabetes[framingham$diabetes == 2] <- 0
  framingham$fr_tot <- framingham$smoker + framingham$total_chol + framingham$hdl + framingham$diabetes 
  framingham <- subset(framingham, select = c("subjects","fr_tot"))
  colnames(framingham) <- c("subjects","Framingham")
  framingham$Framingham <- framingham$Framingham
  framingham$Framingham[framingham$Framingham > 1000] <- NA

#Cumulative Illness Rating Scale-Geriatric (CIRS-G) 
cirsg_HC <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 Clinical Control.xlsx", sheet = "CIRS-G")
cirsg_HC$subjects <- noquote(cirsg_HC$cirsg_id)
cirsg_Clinical <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 Clinical Case.xlsx", sheet = "CIRS-G")
  cirsg_Clinical$subjects <- noquote(cirsg_Clinical$cirsg_id)
  cirsg <- rbind(cirsg_Clinical,cirsg_HC)
  cirsg <- subset(cirsg, select=c(cirsg_id, cirsg_severity_index))
  colnames(cirsg) <- c("subjects","cirsg")
  cirsg$cirsg <- cirsg$cirsg * -1

#BMI
bmi_HC <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 Clinical Control.xlsx", sheet = "Vitals")
bmi_HC$subjects <- noquote(bmi_HC$vs_id)
bmi_Clinical <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 Clinical Case.xlsx", sheet = "Vitals")
  bmi_Clinical$subjects <- noquote(bmi_Clinical$vs_id)
  bmi <- rbind(bmi_Clinical,bmi_HC)
  bmi <- subset(bmi, select=c(vs_id, vs_weight_kg, vs_height_cm))
  colnames(bmi) <- c("subjects","vs_weight_kg","vs_height_cm")
  bmi$BMI <- bmi$vs_weight_kg / ((bmi$vs_height_cm/100)^2) 
  bmi <- subset(bmi, select = c("subjects","BMI"))
  colnames(bmi) <- c("subjects","bmi")
  bmi$bmi[bmi$bmi > 100] <- NA
  
#load depression MADRS score
madrs_HC <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 Clinical Control.xlsx", sheet = "MADRS")
  madrs_HC$subjects <- noquote(madrs_HC$madrs_id)
madrs_Clinical <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 Clinical Case.xlsx", sheet = "MADRS")
  madrs_Clinical$subjects <- noquote(madrs_Clinical$madrs_id)
  madrs <- rbind(madrs_Clinical,madrs_HC)
  madrs <- subset(madrs, select=c(madrs_id, madrs_madrs_total_calc))
  colnames(madrs) <- c("subjects","MADRS")  
  
  
#load fISO & ODI & NDI
gm_CSF_Age_Mat <- ("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/02_ControllingSexEducationNoise_Age_Covariate/Grey_Matter/CSF/HC_MDD_aMCI_NaMCI_aMCI_MDD_NaMCI_MDD_fISO_age_STRUCTresult.mat")
gm_ODI_Age_Mat <-("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/02_ControllingSexEducationNoise_Age_Covariate/Grey_Matter/ODI/analysis_Grey_ODI_age_ci83_Amnestic_STRUCTresult.mat")
descriptives <- read_excel("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/Centered_Scores_PLS_Order_6_groups.xlsx")
age <- subset(descriptives, select = c("subjects","Age"))
descriptives$subjects <- NULL
Mean_NODDI <- read_csv("/Volumes/Seagate_Backup/stats/Analyses_Excluding_Noise_w_Covariates/PLS/Mean_Skeletonised_NODDI_Parameters.csv")

   
#fISO <- gm_CSF_Age_BS <- Beh_PLS_BrainScore(gm_CSF_Age_Mat,"CSF",1) %>%
#  dplyr::select(Subject, LV1)%>%
#  `colnames<-`(c("subjects","fISO"))
Group <- gm_CSF_Age_BS <- Beh_PLS_BrainScore(gm_CSF_Age_Mat,"CSF",1) %>%
  dplyr::select(Subject, Group)%>%
  `colnames<-`(c("subjects","Group"))

#now rename to fit the original data standard
Group$subjects <- str_replace_all(Group$subjects, "sub-CMHH", "PAC01_CMH_")
Group$subjects <- str_replace_all(Group$subjects, "sub-CMH", "PAC01_CMH_")
Group$subjects <- paste(Group$subjects, "_01",sep = "")
Group$subjects <- noquote(Group$subjects)

Mean_NODDI <- merge(Mean_NODDI, Group, by = "subjects")


#ODI <- Beh_PLS_BrainScore(gm_ODI_Age_Mat,"CSF",1)%>%
#  dplyr::select(Subject, LV1)%>%
#  `colnames<-`(c("subjects","ODI"))
#NODDI <- merge(fISO, ODI, by = "subjects")
#NODDI <- merge(NODDI, Group, by = "subjects")
#NODDI$fISO <- NODDI$fISO *-1
#load Cognition Scores
#stroop 
stroop_HC <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 NP Control.xlsx", sheet = "Stroop")
  stroop_HC$subjects <- noquote(stroop_HC$strp_id)
stroop_Clinical <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 NP Case.xlsx", sheet = "Stroop")
  stroop_Clinical$subjects <- noquote(stroop_Clinical$strp_id)
  stroop <- rbind(stroop_Clinical,stroop_HC)
  stroop <- subset(stroop, select=c(strp_id, STRP_CW_Zscore))
  colnames(stroop) <- c("subjects","stroop")
  stroop$stroop <- stroop$stroop 
  
#moca 
moca_HC <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 CLINICAL Control.xlsx", sheet = "MoCA")
  moca_HC$subjects <- noquote(moca_HC$moca_id)
moca_Clinical <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 CLINICAL Case.xlsx", sheet = "MoCA")
  moca_Clinical$subjects <- noquote(moca_Clinical$moca_id)
  moca <- rbind(moca_Clinical,moca_HC)
  moca <- subset(moca, select=c(moca_id, moca_total_score))
  colnames(moca) <- c("subjects","moca")
  moca$moca <- moca$moca
  moca <-replace_with_na(data = moca, replace = list(moca =0))
  
#CDRS 
cdr_HC <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 CLINICAL Control.xlsx", sheet = "CDR")
  cdr_HC$subjects <- noquote(cdr_HC$cdr_id)
cdr_Clinical <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 CLINICAL Case.xlsx", sheet = "CDR")
  cdr_Clinical$subjects <- noquote(cdr_Clinical$cdr_id)
  cdr <- rbind(cdr_Clinical,cdr_HC)
  cdr <- subset(cdr, select=c(cdr_id, cdr_rating))
  colnames(cdr) <- c("subjects","cdr")
  cdr$cdr <- cdr$cdr * -1

#TMT - trail making task Time-per-connection ratio (B/A)
tmt_HC <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 NP Control.xlsx", sheet = "TMT")
  tmt_HC$subjects <- noquote(tmt_HC$trail_id)
tmt_Clinical <- read_excel("/Volumes/Seagate_Backup/stats/subset_with_controls/Correlation Neuropsych/PACtMD T0 NP Case.xlsx", sheet = "TMT")
  tmt_Clinical$subjects <- noquote(tmt_Clinical$trail_id)
  tmt <- rbind(tmt_Clinical,tmt_HC)
  tmt <- subset(tmt, select=c(trail_id, trail_tpc_ratio))
  colnames(tmt) <- c("subjects","tmt")
  tmt$tmt <- tmt$tmt * -1

  
#bind the data together
pm_data <- Reduce(function(x, y) merge(x, y, all=TRUE, by ="subjects"), list(age, framingham,cirsg, madrs, stroop, moca, cdr,tmt, Mean_NODDI))
pm_data <- pm_data[!is.na(pm_data["Group"]),]

#use MICE package to impute the missing values of framingham and stroop...
tempData <- mice(pm_data,m=5,maxit=50,meth='pmm',seed=500)
completedData <- complete(tempData,1)
#reverese code ODI if in same analysis as fISO
completedData$ODI_GM <- completedData$ODI_GM * -1

#pm_data <- pm_data[complete.cases(pm_data), ]
# take a look at the data
head(completedData)

##               subjects   Age Framingham cirsg MADRS     stroop moca  cdr
## 1 PAC01_CMH_0100008_01 79.05       6.98 -1.13     1 -1.5283447   25  0.0
## 2 PAC01_CMH_0100009_01 70.72       5.00 -1.67     0  0.6737996   25 -0.5
## 3 PAC01_CMH_0100010_01 71.16       6.81 -1.50     0 -0.1678005   25 -0.5
## 4 PAC01_CMH_0100017_01 72.20       6.06 -1.88     0 -0.3378685   23  0.0
## 5 PAC01_CMH_0100027_01 70.92       4.53 -1.60     5 -1.6417234   23  0.0
## 6 PAC01_CMH_0100029_01 74.63       6.60 -1.50     3  0.7959184   25  0.0
##         tmt       ODI_GM      NDI_GM     fISO_GM      ODI_WM     NDI_WM
## 1 -1.459459 -0.008500895 0.008022746 0.002497506 0.008888425 0.01281626
## 2 -2.357143 -0.008998484 0.007633086 0.001779048 0.009320335 0.01291227
## 3 -3.391304 -0.008750399 0.007952635 0.002215957 0.009110087 0.01253527
## 4 -1.912281 -0.008916273 0.007770537 0.001658812 0.009412472 0.01284101
## 5 -2.939394 -0.008896544 0.007731707 0.001645582 0.009353363 0.01258100
## 6 -2.611111 -0.008748828 0.007664412 0.001970879 0.009277749 0.01235271
##       fISO_WM       FA_WM     Group
## 1 0.009899803 0.008594833        HC
## 2 0.008722102 0.008814995      aMCI
## 3 0.009424105 0.008355225      aMCI
## 4 0.008563772 0.008853002  MDD_aMCI
## 5 0.008494837 0.008523343      aMCI
## 6 0.009628969 0.008175443 MDD_naMCI

# path matrix (inner model realtionships)

AGE   =c(0, 0, 0)
#HEALTH
NODDI = c(1,0,0)
#DEPRESSION
COG = c(1, 1,0)

rus_path = rbind(AGE, NODDI, COG)
# add optional column names
colnames(rus_path) = rownames(rus_path)
# how does it look like?
rus_path

##       AGE NODDI COG
## AGE     0     0   0
## NODDI   1     0   0
## COG     1     1   0

# plot the path matrix
innerplot(rus_path)

# list indicating what variables are associated with what latent variables
rus_blocks = list(
c("Age"),  
#c(""),
c("fISO_GM"),
c("stroop", "moca", "cdr", "tmt"))

# all latent variables are measured in a reflective way
rus_modes = rep("A", 3)
# run plspm analysis
rus_pls = plspm(completedData, rus_path, rus_blocks, modes = rus_modes, scaled = TRUE, boot.val = TRUE, br = 1000)
# what's in foot_pls?
rus_pls

## Partial Least Squares Path Modeling (PLS-PM) 
## ---------------------------------------------
##    NAME             DESCRIPTION
## 1  $outer_model     outer model
## 2  $inner_model     inner model
## 3  $path_coefs      path coefficients matrix
## 4  $scores          latent variable scores
## 5  $crossloadings   cross-loadings
## 6  $inner_summary   summary inner model
## 7  $effects         total effects
## 8  $unidim          unidimensionality
## 9  $gof             goodness-of-fit
## 10 $boot            bootstrap results
## 11 $data            data matrix
## ---------------------------------------------
## You can also use the function 'summary'

# path coefficients
rus_pls$path_coefs

##               AGE      NODDI COG
## AGE    0.00000000  0.0000000   0
## NODDI  0.52953573  0.0000000   0
## COG   -0.08405664 -0.3234528   0

rus_pls$inner_model

## $NODDI
##               Estimate Std. Error      t value     Pr(>|t|)
## Intercept 1.843561e-17 0.05441768 3.387798e-16 1.000000e+00
## AGE       5.295357e-01 0.05441768 9.730951e+00 4.133716e-19
## 
## $COG
##                Estimate Std. Error       t value     Pr(>|t|)
## Intercept -1.864672e-17 0.05959638 -3.128834e-16 1.000000e+00
## AGE       -8.405664e-02 0.07025492 -1.196452e+00 2.326911e-01
## NODDI     -3.234528e-01 0.07025492 -4.603987e+00 6.697920e-06

# summarized results
summary(rus_pls)

## PARTIAL LEAST SQUARES PATH MODELING (PLS-PM) 
## 
## ---------------------------------------------------------- 
## MODEL SPECIFICATION 
## 1   Number of Cases      245 
## 2   Latent Variables     3 
## 3   Manifest Variables   6 
## 4   Scale of Data        Standardized Data 
## 5   Non-Metric PLS       FALSE 
## 6   Weighting Scheme     centroid 
## 7   Tolerance Crit       1e-06 
## 8   Max Num Iters        100 
## 9   Convergence Iters    3 
## 10  Bootstrapping        TRUE 
## 11  Bootstrap samples    1000 
## 
## ---------------------------------------------------------- 
## BLOCKS DEFINITION 
##     Block         Type   Size   Mode
## 1     AGE    Exogenous      1      A
## 2   NODDI   Endogenous      1      A
## 3     COG   Endogenous      4      A
## 
## ---------------------------------------------------------- 
## BLOCKS UNIDIMENSIONALITY 
##        Mode  MVs  C.alpha  DG.rho  eig.1st  eig.2nd
## AGE       A    1    1.000   1.000      1.0    0.000
## NODDI     A    1    1.000   1.000      1.0    0.000
## COG       A    4    0.618   0.778      1.9    0.883
## 
## ---------------------------------------------------------- 
## OUTER MODEL 
##              weight  loading  communality  redundancy
## AGE                                                  
##   1 Age       1.000    1.000        1.000      0.0000
## NODDI                                                
##   2 fISO_GM   1.000    1.000        1.000      0.2804
## COG                                                  
##   3 stroop    0.212    0.545        0.297      0.0417
##   3 moca      0.501    0.869        0.755      0.1060
##   3 cdr       0.463    0.800        0.639      0.0898
##   3 tmt       0.187    0.424        0.180      0.0253
## 
## ---------------------------------------------------------- 
## CROSSLOADINGS 
##                  AGE    NODDI     COG
## AGE                                  
##   1 Age       1.0000   0.5295  -0.255
## NODDI                                
##   2 fISO_GM   0.5295   1.0000  -0.368
## COG                                  
##   3 stroop   -0.0657  -0.1769   0.545
##   3 moca     -0.2411  -0.3320   0.869
##   3 cdr      -0.2122  -0.3171   0.800
##   3 tmt      -0.1203  -0.0931   0.424
## 
## ---------------------------------------------------------- 
## INNER MODEL 
## $NODDI
##             Estimate   Std. Error    t value   Pr(>|t|)
## Intercept   1.84e-17       0.0544   3.39e-16   1.00e+00
## AGE         5.30e-01       0.0544   9.73e+00   4.13e-19
## 
## $COG
##              Estimate   Std. Error     t value   Pr(>|t|)
## Intercept   -1.86e-17       0.0596   -3.13e-16   1.00e+00
## AGE         -8.41e-02       0.0703   -1.20e+00   2.33e-01
## NODDI       -3.23e-01       0.0703   -4.60e+00   6.70e-06
## 
## ---------------------------------------------------------- 
## CORRELATIONS BETWEEN LVs 
##           AGE   NODDI     COG
## AGE     1.000   0.529  -0.255
## NODDI   0.529   1.000  -0.368
## COG    -0.255  -0.368   1.000
## 
## ---------------------------------------------------------- 
## SUMMARY INNER MODEL 
##              Type    R2  Block_Communality  Mean_Redundancy    AVE
## AGE     Exogenous  0.00              1.000           0.0000  1.000
## NODDI  Endogenous  0.28              1.000           0.2804  1.000
## COG    Endogenous  0.14              0.468           0.0657  0.468
## 
## ---------------------------------------------------------- 
## GOODNESS-OF-FIT 
## [1]  0.3137
## 
## ---------------------------------------------------------- 
## TOTAL EFFECTS 
##    relationships   direct  indirect   total
## 1   AGE -> NODDI   0.5295     0.000   0.530
## 2     AGE -> COG  -0.0841    -0.171  -0.255
## 3   NODDI -> COG  -0.3235     0.000  -0.323
## 
## --------------------------------------------------------- 
## BOOTSTRAP VALIDATION 
## weights 
##                Original  Mean.Boot  Std.Error  perc.025  perc.975
## AGE-Age           1.000      1.000   1.15e-16    1.0000     1.000
## NODDI-fISO_GM     1.000      1.000   1.02e-16    1.0000     1.000
## COG-stroop        0.212      0.207   7.72e-02    0.0506     0.353
## COG-moca          0.501      0.494   5.08e-02    0.3997     0.592
## COG-cdr           0.463      0.456   5.87e-02    0.3409     0.576
## COG-tmt           0.187      0.197   8.00e-02    0.0466     0.369
## 
## loadings 
##                Original  Mean.Boot  Std.Error  perc.025  perc.975
## AGE-Age           1.000      1.000   1.11e-16     1.000     1.000
## NODDI-fISO_GM     1.000      1.000   1.11e-16     1.000     1.000
## COG-stroop        0.545      0.541   8.14e-02     0.363     0.686
## COG-moca          0.869      0.860   2.90e-02     0.796     0.908
## COG-cdr           0.800      0.792   4.55e-02     0.693     0.871
## COG-tmt           0.424      0.438   8.61e-02     0.266     0.624
## 
## paths 
##               Original  Mean.Boot  Std.Error  perc.025  perc.975
## AGE -> NODDI    0.5295      0.525     0.0528     0.418    0.6231
## AGE -> COG     -0.0841     -0.089     0.0696    -0.228    0.0446
## NODDI -> COG   -0.3235     -0.329     0.0589    -0.441   -0.2174
## 
## rsq 
##        Original  Mean.Boot  Std.Error  perc.025  perc.975
## NODDI      0.28      0.278     0.0550    0.1747     0.388
## COG        0.14      0.153     0.0376    0.0872     0.228
## 
## total.efs 
##               Original  Mean.Boot  Std.Error  perc.025  perc.975
## AGE -> NODDI     0.530      0.525     0.0528     0.418     0.623
## AGE -> COG      -0.255     -0.262     0.0603    -0.379    -0.146
## NODDI -> COG    -0.323     -0.329     0.0589    -0.441    -0.217

# plot the results (inner model)
plot(rus_pls)

# plot the loadings of the outer model
plot(rus_pls, what = "loadings", arr.width = 0.1)

# plot the weights of the outer model
plot(rus_pls, what = "weights", arr.width = 0.1)

# load ggplot2 and reshape
library(ggplot2)
library(reshape)
library(ggthemr)
ggthemr("fresh")
# reshape crossloadings data.frame for ggplot
xloads = melt(rus_pls$crossloadings, id.vars = c("name", "block"),
variable_name = "LV")
# bar-charts of crossloadings by block
ggplot(data = xloads,
aes(x = name, y = value, fill = block)) +
geom_hline(yintercept = 0, color = "gray75") +
geom_hline(yintercept = c(-0.5, 0.5), color = "gray70", linetype = 2) +
geom_bar(stat = 'identity', position = 'dodge') +
facet_wrap(block ~ LV) +
theme(axis.text.x = element_text(angle = 90),
line = element_blank()) +
ggtitle("Crossloadings")

 # permutation test with 100 permutations
completedData$MCI <- str_detect(completedData$Group, "MCI")
completedData$MCI <- as.factor(completedData$MCI)
completedData$MCI <- dplyr::recode(completedData$MCI, "TRUE" = "MCI", "FALSE" = "Non-MCI")
MCI_perm = plspm.groups(rus_pls, completedData$MCI,method="permutation", reps=100)
MCI_perm

## GROUP COMPARISON IN PLS-PM FOR PATH COEFFICIENTS 
## 
## Scale of Data:       TRUE 
## Weighting Scheme:    centroid 
## Selected method:     permutation 
## Num of replicates:   100 
## 
## $test 
##              global  group.Non-MCI  group.MCI  diff.abs  p.value  sig.05
## AGE->NODDI   0.5295         0.4700     0.5293    0.0594   0.6436      no
## AGE->COG    -0.0841        -0.1629    -0.0578    0.1051   0.5545      no
## NODDI->COG  -0.3235        -0.2253    -0.3045    0.0792   0.6436      no
## 
## Inner models in the following objects: 
## $global  
## $group1  
## $group2

#
barplot(t(as.matrix(MCI_perm$test[,2:3])), border = NA, beside = TRUE,
col = c("#FEB24C","#74A9CF"), las = 2, ylim = c(-0.1, 1),cex.names = 0.8, col.axis = "gray30", cex.axis = 0.8)
# add horizontal line
abline(h = 0, col = "gray50")
# add itle
title("Path coefficients of MCI and Non-MCI Individuals",
cex.main = 0.95, col.main = "gray30")
# add legend
legend("top", legend = c("Non-MCI", "MCI"), pt.bg = c("#FEB24C", "#A6BDDB"),
ncol = 2, pch = 22, col = c("#FEB24C", "#74A9CF"), bty = "n",
text.col = "gray40")

completedData$MDD <- str_detect(completedData$Group, "MDD")
completedData$MDD <- as.factor(completedData$MDD)
completedData$MDD <- dplyr::recode(completedData$MDD, "TRUE" = "MDD", "FALSE" = "Non-MDD")
MDD_perm = plspm.groups(rus_pls, completedData$MDD,method="permutation", reps=100)
MDD_perm

## GROUP COMPARISON IN PLS-PM FOR PATH COEFFICIENTS 
## 
## Scale of Data:       TRUE 
## Weighting Scheme:    centroid 
## Selected method:     permutation 
## Num of replicates:   100 
## 
## $test 
##              global  group.Non-MDD  group.MDD  diff.abs  p.value  sig.05
## AGE->NODDI   0.5295         0.5887     0.3697    0.2191   0.0297     yes
## AGE->COG    -0.0841        -0.1112    -0.1279    0.0167   0.9109      no
## NODDI->COG  -0.3235        -0.2446    -0.3349    0.0903   0.4752      no
## 
## Inner models in the following objects: 
## $global  
## $group1  
## $group2

barplot(t(as.matrix(MDD_perm$test[,2:3])), border = NA, beside = TRUE,
col = c("#FEB24C","#74A9CF"), las = 2, ylim = c(-0.1, 1),cex.names = 0.8, col.axis = "gray30", cex.axis = 0.8)
# add horizontal line
abline(h = 0, col = "gray50")
# add itle
title("Path coefficients of MDD and Non-MDD Individuals",
cex.main = 0.95, col.main = "gray30")
# add legend
legend("top", legend = c("Non-MDD", "MDD"), pt.bg = c("#FEB24C", "#A6BDDB"),
ncol = 2, pch = 22, col = c("#FEB24C", "#74A9CF"), bty = "n",
text.col = "gray40")

completedData$aMCI <- str_detect(completedData$Group, "aMCI")
completedData$aMCI <- as.factor(completedData$aMCI)
completedData$aMCI <- dplyr::recode(completedData$aMCI, "TRUE" = "aMCI", "FALSE" = "Other")
aMCI_perm = plspm.groups(rus_pls, completedData$aMCI,method="permutation", reps=100)
aMCI_perm

## GROUP COMPARISON IN PLS-PM FOR PATH COEFFICIENTS 
## 
## Scale of Data:       TRUE 
## Weighting Scheme:    centroid 
## Selected method:     permutation 
## Num of replicates:   100 
## 
## $test 
##              global  group.Other  group.aMCI  diff.abs  p.value  sig.05
## AGE->NODDI   0.5295       0.4700      0.5293    0.0594   0.5941      no
## AGE->COG    -0.0841      -0.1629     -0.0578    0.1051   0.5050      no
## NODDI->COG  -0.3235      -0.2253     -0.3045    0.0792   0.5644      no
## 
## Inner models in the following objects: 
## $global  
## $group1  
## $group2

barplot(t(as.matrix(aMCI_perm$test[,2:3])), border = NA, beside = TRUE,
col = c("#FEB24C","#74A9CF"), las = 2, ylim = c(-0.1, 1),cex.names = 0.8, col.axis = "gray30", cex.axis = 0.8)
# add horizontal line
abline(h = 0, col = "gray50")
# add itle
title("Path coefficients of everyone vs aMCI Individuals",
cex.main = 0.95, col.main = "gray30")
# add legend
legend("top", legend = c("aMCI", "Other"), pt.bg = c("#FEB24C", "#A6BDDB"),
ncol = 2, pch = 22, col = c("#FEB24C", "#74A9CF"), bty = "n",
text.col = "gray40")