Brain Shuffling and Unshuffling (Nested 20 Times): Power and Type 1 Error

library(reshape2)
library(ggplot2)
library(colorspace)
nshuffH1<-c(rep(0,6),rep(50,5), rep(150,4), rep(200,3), rep(250,2), rep(350,1))
nMC <- 5
files <- 20

total <- nMC * files
#loading 20 files with 5nMC each (total of 100 nMC) and each of the 20 files is bootstrapped 200 times (4000 bootstraps in total)

load("Results/PPP2aSeparate2.RData")
PPP2a_Sep2 <- PPP2a_Sep1
load("Results/PPP2aSeparate3.RData")
PPP2a_Sep3 <- PPP2a_Sep1
load("Results/PPP2aSeparate4.RData")
PPP2a_Sep4 <- PPP2a_Sep1
load("Results/PPP2aSeparate5.RData")
PPP2a_Sep5 <- PPP2a_Sep1
load("Results/PPP2aSeparate6.RData")
PPP2a_Sep6 <- PPP2a_Sep1
load("Results/PPP2aSeparate7.RData")
PPP2a_Sep7 <- PPP2a_Sep1
load("Results/PPP2aSeparate8.RData")
PPP2a_Sep8 <- PPP2a_Sep1
load("Results/PPP2aSeparate9.RData")
PPP2a_Sep9 <- PPP2a_Sep1
load("Results/PPP2aSeparate10.RData")
PPP2a_Sep10 <- PPP2a_Sep1
load("Results/PPP2aSeparate11.RData")
PPP2a_Sep11 <- PPP2a_Sep1
load("Results/PPP2aSeparate12.RData")
PPP2a_Sep12 <- PPP2a_Sep1
load("Results/PPP2aSeparate13.RData")
PPP2a_Sep13 <- PPP2a_Sep1
load("Results/PPP2aSeparate14.RData")
PPP2a_Sep14 <- PPP2a_Sep1
load("Results/PPP2aSeparate15.RData")
PPP2a_Sep15 <- PPP2a_Sep1
load("Results/PPP2aSeparate16.RData")
PPP2a_Sep16 <- PPP2a_Sep1
load("Results/PPP2aSeparate17.RData")
PPP2a_Sep17 <- PPP2a_Sep1
load("Results/PPP2aSeparate18.RData")
PPP2a_Sep18 <- PPP2a_Sep1
load("Results/PPP2aSeparate19.RData")
PPP2a_Sep19 <- PPP2a_Sep1
load("Results/PPP2aSeparate20.RData")
PPP2a_Sep20 <- PPP2a_Sep1
load("Results/PPP2aSeparate1.RData")

PPP2a <- (PPP2a_Sep1 + PPP2a_Sep2 + PPP2a_Sep3 + PPP2a_Sep4 + PPP2a_Sep5 + PPP2a_Sep6 + PPP2a_Sep7 + PPP2a_Sep8 + PPP2a_Sep9 + PPP2a_Sep10 + PPP2a_Sep11 + PPP2a_Sep12 + PPP2a_Sep13 + PPP2a_Sep14 + PPP2a_Sep15 + PPP2a_Sep16 + PPP2a_Sep17 + PPP2a_Sep18 + + PPP2a_Sep19 + PPP2a_Sep20)/total


load("Results/PPP2nSeparate2.RData")
PPP2n_Sep2 <- PPP2n_Sep1
load("Results/PPP2nSeparate3.RData")
PPP2n_Sep3 <- PPP2n_Sep1
load("Results/PPP2nSeparate4.RData")
PPP2n_Sep4 <- PPP2n_Sep1
load("Results/PPP2nSeparate5.RData")
PPP2n_Sep5 <- PPP2n_Sep1
load("Results/PPP2nSeparate6.RData")
PPP2n_Sep6 <- PPP2n_Sep1
load("Results/PPP2nSeparate7.RData")
PPP2n_Sep7 <- PPP2n_Sep1
load("Results/PPP2nSeparate8.RData")
PPP2n_Sep8 <- PPP2n_Sep1
load("Results/PPP2nSeparate9.RData")
PPP2n_Sep9 <- PPP2n_Sep1
load("Results/PPP2nSeparate10.RData")
PPP2n_Sep10 <- PPP2n_Sep1
load("Results/PPP2nSeparate11.RData")
PPP2n_Sep11 <- PPP2n_Sep1
load("Results/PPP2nSeparate12.RData")
PPP2n_Sep12 <- PPP2n_Sep1
load("Results/PPP2nSeparate13.RData")
PPP2n_Sep13 <- PPP2n_Sep1
load("Results/PPP2nSeparate14.RData")
PPP2n_Sep14 <- PPP2n_Sep1
load("Results/PPP2nSeparate15.RData")
PPP2n_Sep15 <- PPP2n_Sep1
load("Results/PPP2nSeparate16.RData")
PPP2n_Sep16 <- PPP2n_Sep1
load("Results/PPP2nSeparate17.RData")
PPP2n_Sep17 <- PPP2n_Sep1
load("Results/PPP2nSeparate18.RData")
PPP2n_Sep18 <- PPP2n_Sep1
load("Results/PPP2nSeparate19.RData")
PPP2n_Sep19 <- PPP2n_Sep1
load("Results/PPP2nSeparate20.RData")
PPP2n_Sep20 <- PPP2n_Sep1
load("Results/PPP2nSeparate1.RData")

PPP2n <- (PPP2n_Sep1 + PPP2n_Sep2 + PPP2n_Sep3 + PPP2n_Sep4 + PPP2n_Sep5 + PPP2n_Sep6 + PPP2n_Sep7 + PPP2n_Sep8 + PPP2n_Sep9 + PPP2n_Sep10 + PPP2n_Sep11 + PPP2n_Sep12 + PPP2n_Sep13 + PPP2n_Sep14 + PPP2n_Sep15 + PPP2n_Sep16 + PPP2n_Sep17 + PPP2n_Sep18 + + PPP2n_Sep19 + PPP2n_Sep20)/total


load("Results/PPP2a_unshuffSeparate2.RData")
PPP2a_unshuff_Sep2 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate3.RData")
PPP2a_unshuff_Sep3 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate4.RData")
PPP2a_unshuff_Sep4 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate5.RData")
PPP2a_unshuff_Sep5 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate6.RData")
PPP2a_unshuff_Sep6 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate7.RData")
PPP2a_unshuff_Sep7 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate8.RData")
PPP2a_unshuff_Sep8 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate9.RData")
PPP2a_unshuff_Sep9 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate10.RData")
PPP2a_unshuff_Sep10 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate11.RData")
PPP2a_unshuff_Sep11 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate12.RData")
PPP2a_unshuff_Sep12 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate13.RData")
PPP2a_unshuff_Sep13 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate14.RData")
PPP2a_unshuff_Sep14 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate15.RData")
PPP2a_unshuff_Sep15 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate16.RData")
PPP2a_unshuff_Sep16 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate17.RData")
PPP2a_unshuff_Sep17 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate18.RData")
PPP2a_unshuff_Sep18 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate19.RData")
PPP2a_unshuff_Sep19 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate20.RData")
PPP2a_unshuff_Sep20 <- PPP2a_unshuff_Sep1
load("Results/PPP2a_unshuffSeparate1.RData")

PPP2a_unshuff <- (PPP2a_unshuff_Sep1 + PPP2a_unshuff_Sep2 + PPP2a_unshuff_Sep3 + PPP2a_unshuff_Sep4 + PPP2a_unshuff_Sep5 + PPP2a_unshuff_Sep6 + PPP2a_unshuff_Sep7 + PPP2a_unshuff_Sep8 + PPP2a_unshuff_Sep9 + PPP2a_unshuff_Sep10 + PPP2a_unshuff_Sep11 + PPP2a_unshuff_Sep12 + PPP2a_unshuff_Sep13 + PPP2a_unshuff_Sep14 + PPP2a_unshuff_Sep15 + PPP2a_unshuff_Sep16 + PPP2a_unshuff_Sep17 + PPP2a_unshuff_Sep18 + + PPP2a_unshuff_Sep19 + PPP2a_unshuff_Sep20)/total


load("Results/PPP2n_unshuffSeparate2.RData")
PPP2n_unshuff_Sep2 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate3.RData")
PPP2n_unshuff_Sep3 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate4.RData")
PPP2n_unshuff_Sep4 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate5.RData")
PPP2n_unshuff_Sep5 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate6.RData")
PPP2n_unshuff_Sep6 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate7.RData")
PPP2n_unshuff_Sep7 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate8.RData")
PPP2n_unshuff_Sep8 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate9.RData")
PPP2n_unshuff_Sep9 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate10.RData")
PPP2n_unshuff_Sep10 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate11.RData")
PPP2n_unshuff_Sep11 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate12.RData")
PPP2n_unshuff_Sep12 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate13.RData")
PPP2n_unshuff_Sep13 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate14.RData")
PPP2n_unshuff_Sep14 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate15.RData")
PPP2n_unshuff_Sep15 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate16.RData")
PPP2n_unshuff_Sep16 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate17.RData")
PPP2n_unshuff_Sep17 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate18.RData")
PPP2n_unshuff_Sep18 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate19.RData")
PPP2n_unshuff_Sep19 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate20.RData")
PPP2n_unshuff_Sep20 <- PPP2n_unshuff_Sep1
load("Results/PPP2n_unshuffSeparate1.RData")

PPP2n_unshuff <- (PPP2n_unshuff_Sep1 + PPP2n_unshuff_Sep2 + PPP2n_unshuff_Sep3 + PPP2n_unshuff_Sep4 + PPP2n_unshuff_Sep5 + PPP2n_unshuff_Sep6 + PPP2n_unshuff_Sep7 + PPP2n_unshuff_Sep8 + PPP2n_unshuff_Sep9 + PPP2n_unshuff_Sep10 + PPP2n_unshuff_Sep11 + PPP2n_unshuff_Sep12 + PPP2n_unshuff_Sep13 + PPP2n_unshuff_Sep14 + PPP2n_unshuff_Sep15 + PPP2n_unshuff_Sep16 + PPP2n_unshuff_Sep17 + PPP2n_unshuff_Sep18 + + PPP2n_unshuff_Sep19 + PPP2n_unshuff_Sep20)/total

############################### Power Testing Brain Graphs Shuffled ##############################
pow_a <- na.omit(melt(PPP2a))$value
  k <- c(0, 50, 150, 200, 250, 350)

X<-cbind(pow_a,nshuffH1)
nshuffH0<-c(0,50,150,200,250,350, 50,150,200,250,350, 150,200,250,350, 200,250,350, 250,350, 350)
nshuffH1<-c(rep(0,6),rep(50,5), rep(150,4), rep(200,3), rep(250,2), rep(350,1))
X<-cbind(X,nshuffH0)
X<-data.frame(X)
  
this <- rev(qualitative_hcl(6, "Dark 3"))
pd <- position_dodge(.65)
X$nshuffH0<-as.factor(X$nshuffH0)

gpa<-ggplot(X, aes(x=nshuffH1, y=pow_a, color=nshuffH0, shape=nshuffH0))+theme_bw()
gpa<-gpa + geom_line(aes(color=nshuffH0),size=1.05) +  
  geom_point(aes(shape=nshuffH0),size=4) +
  geom_errorbar(aes(ymin=pow_a-2*sqrt(pow_a*(1-pow_a)/100),ymax=pow_a+2*sqrt(pow_a*(1-pow_a)/100)),size=1.05,width=10 )+
  labs(family="serif", title="Power Testing with Brain Graphs", x="Number of Vertices Shuffled in H1", y="Power", color="Number shuffled\n in H0",shape="Number shuffled\n in H0") +  theme_classic() + 
    theme(plot.title = element_text(family="serif", face = "bold", size = 17, hjust=0.5), 
        panel.background = element_rect(fill = "#FFFFFF"), 
        legend.title = element_text(family="serif", face="bold", color = "#000000", size = 12),
        legend.background = element_rect(fill = "#FFFFFF"), 
        text = element_text(size=20, family = "serif"),
        axis.title = element_text(family = "serif", size = 15, colour = "#000000"),
        plot.background = element_rect(fill = "#FFFFFF"),
        axis.line = element_line(colour = "#000000", size = 1), 
        legend.text = element_text( family="serif", color = "#000000", size=15)) +  
        scale_color_manual(values=c("0" = "red", "50"="#B675E0","150"="#00A6CA", "200"="#00AA5A", "250"="#AA9000", "350"="#E16A86"))+
        scale_shape_manual(values=c("0" = 1, "50"= 15,"150"=5, "200"= 20, "250"= 8, "350"= 20))

############################### Type 1 Error Brain Graphs Shuffled ##############################
pow_n <- na.omit(melt(PPP2n))$value
#nshuffH1<-c(rep(2,6),rep(25,5), rep(50,4), rep(100,3), rep(150,2), rep(250,1))
Y<-cbind(pow_n,nshuffH1)
nshuffH0<-c(0,50,150,200,250,350, 50,150,200,250,350, 150,200,250,350, 200,250,350, 250,350, 350)
nshuffH1<-c(rep(0,6),rep(50,5), rep(150,4), rep(200,3), rep(250,2), rep(350,1))
Y<-cbind(Y,nshuffH0)
Y<-data.frame(Y)
Y$nshuffH0<-as.factor(Y$nshuffH0)

gpa2<-ggplot(Y, aes(x=nshuffH1, y=pow_n, color=nshuffH0, shape=nshuffH0))+theme_bw()
gpa2<-gpa2 + geom_line(aes(color=nshuffH0),size=1.05) +  
  geom_point(aes(shape=nshuffH0),size=4) + ylim(-0.01,0.2) +
  geom_errorbar(aes(ymin=pow_n-2*sqrt(pow_n*(1-pow_n)/100),ymax=pow_n+2*sqrt(pow_n*(1-pow_n)/100)),size=1.05,width=10 )+
  labs(family="serif", title="Type I error: Testing with Brain Graphs", x="Number of Vertices Shuffled in H1", y="Type-I error", color="Number shuffled\n in H0",shape="Number shuffled\n in H0") +  theme_classic() + 
    theme(plot.title = element_text(family="serif", face = "bold", size = 17, hjust=0.5), 
        panel.background = element_rect(fill = "#FFFFFF"), 
        legend.title = element_text(family="serif", face="bold", color = "#000000", size = 12),
        legend.background = element_rect(fill = "#FFFFFF"), 
        text = element_text(size=20, family = "serif"),
        axis.title = element_text(family = "serif", size = 15, colour = "#000000"),
        plot.background = element_rect(fill = "#FFFFFF"),
        axis.line = element_line(colour = "#000000", size = 1), 
        legend.text = element_text( family="serif", color = "#000000", size=15)) +  
        scale_color_manual(values=c("0" = "red", "50"="#B675E0","150"="#00A6CA", "200"="#00AA5A", "250"="#AA9000", "350"="#E16A86"))+
        scale_shape_manual(values=c("0" = 1, "50"= 15,"150"=5, "200"= 20, "250"= 8, "350"= 20)) + 
        geom_hline(yintercept=0.05, linetype="dashed", color = "black", size=1)

############################### Power Testing Brain Graphs Unshuffled ##############################
pow_a_unshuff <- na.omit(melt(PPP2a_unshuff))$value
#nshuffH1<-c(rep(2,6),rep(25,5), rep(50,4), rep(100,3), rep(150,2), rep(250,1))
X<-cbind(pow_a_unshuff,nshuffH1)
nshuffH0<-c(0,50,150,200,250,350, 50,150,200,250,350, 150,200,250,350, 200,250,350, 250,350, 350)
nshuffH1<-c(rep(0,6),rep(50,5), rep(150,4), rep(200,3), rep(250,2), rep(350,1))
X<-cbind(X,nshuffH0)
X<-data.frame(X)
  
library(ggplot2)
library(colorspace)
this <- rev(qualitative_hcl(6, "Dark 3"))
pd <- position_dodge(.65)
X$nshuffH0<-as.factor(X$nshuffH0)

gpa3<-ggplot(X, aes(x=nshuffH1, y=pow_a_unshuff, color=nshuffH0, shape=nshuffH0))+theme_bw()
gpa3<-gpa3 + geom_line(aes(color=nshuffH0),size=1.05) +  
  geom_errorbar(aes(ymin=pow_a_unshuff-2*sqrt(pow_a_unshuff*(1-pow_a_unshuff)/100),ymax=pow_a_unshuff+2*sqrt(pow_a_unshuff*(1-pow_a_unshuff)/100)),size=1.05,width=10 )+
  geom_point(aes(shape=nshuffH0),size=4) + ylim(0,1) +
  labs(family="serif", title="Power Testing with Brain\n Graphs Unshuffled", x="Number of Vertices Unshuffled in H1", y="Power", color="Number shuffled\n in H0",shape="Number shuffled\n in H0") +  theme_classic() + 
    theme(plot.title = element_text(family="serif", face = "bold", size = 17, hjust=0.5), 
        panel.background = element_rect(fill = "#FFFFFF"), 
        legend.title = element_text(family="serif", face="bold", color = "#000000", size = 12),
        legend.background = element_rect(fill = "#FFFFFF"), 
        text = element_text(size=20, family = "serif"),
        axis.title = element_text(family = "serif", size = 15, colour = "#000000"),
        plot.background = element_rect(fill = "#FFFFFF"),
        axis.line = element_line(colour = "#000000", size = 1), 
        legend.text = element_text( family="serif", color = "#000000", size=15)) +  
        scale_color_manual(values=c("0" = "red", "50"="#B675E0","150"="#00A6CA", "200"="#00AA5A", "250"="#AA9000", "350"="#E16A86"))+
        scale_shape_manual(values=c("0" = 1, "50"= 15,"150"=5, "200"= 20, "250"= 8, "350"= 20))

############################### Type 1 Error Brain Graphs Unshuffled ##############################
pow_n_unshuff <- na.omit(melt(PPP2n_unshuff))$value
#nshuffH1<-c(rep(2,6),rep(25,5), rep(50,4), rep(100,3), rep(150,2), rep(250,1))
Y<-cbind(pow_n_unshuff,nshuffH1)
nshuffH0<-c(0,50,150,200,250,350, 50,150,200,250,350, 150,200,250,350, 200,250,350, 250,350, 350)
nshuffH1<-c(rep(0,6),rep(50,5), rep(150,4), rep(200,3), rep(250,2), rep(350,1))
Y<-cbind(Y,nshuffH0)
Y<-data.frame(Y)
Y$nshuffH0<-as.factor(Y$nshuffH0)

gpa4<-ggplot(Y, aes(x=nshuffH1, y=pow_n_unshuff, color=nshuffH0, shape=nshuffH0))+theme_bw()
gpa4<-gpa4 + geom_line(aes(color=nshuffH0),size=1.05) + 
  geom_errorbar(aes(ymin=pow_n_unshuff-2*sqrt(pow_n_unshuff*(1-pow_n_unshuff)/100),ymax=pow_n_unshuff+2*sqrt(pow_n_unshuff*(1-pow_n_unshuff)/100)),size=1.05,width=10 )+
  geom_point(aes(shape=nshuffH0),size=4) + ylim(0,1) +
  labs(family="serif", title="Type I error: Testing with Brain\n Graphs Unshuffled", x="Number of Vertices Unshuffled in H1", y="Type-I error", color="Number shuffled\n in H0",shape="Number shuffled\n in H0") +  theme_classic() + 
    theme(plot.title = element_text(family="serif", face = "bold", size = 17, hjust=0.5), 
        panel.background = element_rect(fill = "#FFFFFF"), 
        legend.title = element_text(family="serif", face="bold", color = "#000000", size = 12),
        legend.background = element_rect(fill = "#FFFFFF"), 
        text = element_text(size=20, family = "serif"),
        axis.title = element_text(family = "serif", size = 15, colour = "#000000"),
        plot.background = element_rect(fill = "#FFFFFF"),
        axis.line = element_line(colour = "#000000", size = 1), 
        legend.text = element_text( family="serif", color = "#000000", size=15)) +  
        scale_color_manual(values=c("0" = "red", "50"="#B675E0","150"="#00A6CA", "200"="#00AA5A", "250"="#AA9000", "350"="#E16A86"))+
        scale_shape_manual(values=c("0" = 1, "50"= 15,"150"=5, "200"= 20, "250"= 8, "350"= 20)) + 
       geom_hline(yintercept=0.05, linetype="dashed", color = "black", size=1)


library(patchwork)
gpa+ gpa2 + plot_layout(ncol = 2, guides ="collect")

gpa3+ gpa4 + plot_layout(ncol = 2, guides ="collect")