#load required packages
library(bootnet) 
library(networktools)
library(NetworkComparisonTest)
library(qgraph)

#Load data
Irinadata <- read.table("C:/Users/lapab/Dropbox/EAT Lab/AED Webinar/clinicaldata.csv", header=TRUE, sep=",", na = "NA")
#check data 
summary(Irinadata)

#Assign names to nodes
mynames <-c("restrict",	"fast",	"avoidfood",	"foodrules", "concentrate_cal",
          "losecontrol", "binge",	"eatsecret",	"flatstomach",	"emptystom",
          "concentrate_ws", "feargain",	"feelfat",	"desirelose",	"guilty",
          "judgeweight",	"judgeshape",	"weigh",	"dissatisweight",
          "dissatisfshape",	"seeeat",	"seebody",	"othersee")

#Estimate network using default methods
mynetwork <- estimateNetwork(Irinadata, default="EBICglasso") #use polychoric correlation

#Use Spearman correlation
mynetwork <- estimateNetwork(data, default="EBICglasso", corMethod = "cor", corArgs =
                              list(method = "spearman", use = "pairwise.complete.obs"))

#Plot network
plot(mynetwork, layout="spring", vsize=6, color="lightblue", border.color="black", nodeNames = names, legend=FALSE)

#Save plot as pdf
setwd("C:/Users/lapab/Dropbox/EAT Lab/AED Webinar/Results") #set directory for where to save file
pdf("MyNetwork.pdf")
myplot <- plot(mynetwork, layout="spring", vsize=6, color="lightblue", border.color="black", nodeNames = mynames, legend=FALSE)
dev.off()

#Create centrality plot (will show strength centrality)
pdf("MyCentrality.pdf",width=4)
c1 <- centralityPlot(myplot)
dev.off()

#Expected influence plot
pdf("MyExpectedInfluence.pdf", width=4)
C2 <- centralityPlot(myplot, include = "ExpectedInfluence")
dev.off()

#Save centrality values
CentralityTable <- centralityTable(mynetwork)
write.csv(CentralityTable, "MyCentralityTable.csv")

#Constructing a partial correlation matrix
N1edges <-getWmat(mynetwork)
write.csv(N1edges, "NetworkEdges.csv")

#Estimating Network Stability
b1 <- bootnet(mynetwork, boots=1000,nCores=4, statistics=c("strength", "expectedInfluence", "edge"))
b2 <- bootnet(mynetwork, boots=1000,nCores=4, type="case", statistics=c("strength", "expectedInfluence", "edge"))

#Save bootstrapped files
save(b1, file = "b1.Rdata")
save(b2, file = "b2.Rdata")

#load bootstrapped files
setwd("C:/Users/lapab/Dropbox/EAT Lab/AED Webinar")
load("b1.Rdata")
load("b2.Rdata")

#Get centrality stability coefficient
corStability(b2) 

#Save edge stability graph
pdf("EdgeStability.pdf")
plot(b1, labels = FALSE, order = "sample") 
dev.off()

#Save centrality stability graph
pdf("CentrStability.pdf") 
plot(b2)
dev.off()

# Strength Centrality diff test
pdf("CentraityDifference.pdf")
plot(b1, "strength", order="sample", labels=TRUE) 
dev.off()

# EI diff test
pdf("EIDifference.pdf")
plot(b1, "expectedInfluence", order="sample", labels=TRUE) 
dev.off()

#Edge weights diff test
pdf("Difftest.pdf")
plot(b1, "edge", plot = "difference", onlyNonZero = TRUE, order = "sample")
dev.off()

###END OF STABILITY ANALYSIS###

#Identifying bridge symptoms

#Load data
Irinadata <- read.table("C:/Users/lapab/Dropbox/EAT Lab/AED Webinar/bridgedata.csv", header=TRUE, sep=",", na = "NA")
#check data 
summary(Irinadata)

#assign variable names
mynames  <- c("Restrict",	"NotEat",	"AvoidEat",	"EatRules",	"FoodConc",	"StopFear",	
            "NoControl",	"Binge",	"SecretEat",	"FlatStom",	"EmpStom",	"WtShpConc",	
            "FearGain",	"FeelFat",	"WantLose",	"WtShpGuilt",	"WtThink", "ShpThink",	
            "WeighDist",	"WtDiss",	"ShpDiss",	"ThinDesire",	"SeeEat",	"SeeBody",	
            "Othersee",	"IntThought",	"Dreams", "ReExp",	"UpsetRem",	"PhysRem", 
            "AvoidTS",	"AvoidAS",	"MemProb", "LossInt",	"FeelDist",	"EmotNum", 
            "Future", "Sleep", "Irritable", "ConcProb", "HyperV", "Startle")

#create groups
mygroups=list("ED"=c(1:25),"PTSD"=c(26:42)) 

#estimate network
mynetwork <- estimateNetwork(Irinadata, default="EBICglasso")

#Plot network
myplot <-plot(mynetwork, layout="spring", vsize=6, border.color="black", 
                 groups=mygroups, labels=mynames, color=c('#a8e6cf', '#dcedc1'), legend.cex=.4)

#Constructing a partial correlation matrix
myedges <-getWmat(mynetwork)
write.csv(myedges, "MyNetworkEdges.csv")

#Estimate bridge values for each node
bridge(myplot, communities=c('1','1','1','1','1','1','1','1','1','1',
                              '1','1','1','1','1','1','1','1','1','1',
                              '1','1','1','1','1','2','2','2','2','2','2','2',
                              '2','2','2','2','2','2','2','2','2','2'), useCommunities = "all", directed = NULL, nodes = NULL)
#Name our bridge object
mybridge <- bridge(myplot, communities=c('1','1','1','1','1','1','1','1','1','1',
                                        '1','1','1','1','1','1','1','1','1','1',
                                        '1','1','1','1','1','2','2','2','2','2','2','2',
                                        '2','2','2','2','2','2','2','2','2','2'), useCommunities = "all", directed = NULL, nodes = NULL)

#Create bridge graph
pdf("bridgecentrality.pdf", width=4)
plot(mybridge, include = "Bridge Strength")
dev.off()

#Create bridge expected influence graph
pdf("bridgeEI.pdf", width=4)
plot(mybridge, include = "Bridge Expected Influence (1-step)")
dev.off()

#Bridge stability part 1
caseDroppingBoot <- bootnet(mynetwork, boots=1000, type="case", 
                            statistics=c("bridgeStrength", "bridgeExpectedInfluence"), 
                            communities=groupsint)
#get stability coefficients
corStability(caseDroppingBoot)

#Plot centrality stability
plot(caseDroppingBoot, statistics="bridgeExpectedInfluence")

#Bridge stability part 2; centraity difference
nonParametricBoot <- bootnet(mynetwork, boots=1000, type="nonparametric", statistics=c("bridgeStrength", "bridgeExpectedInfluence"), 
                             communities=groupsint)
#Plot centrality difference
plot(nonParametricBoot, statistics="bridgeExpectedInfluence", plot="difference")
plot(nonParametricBoot, statistics="bridgeStrength", plot="difference")

###END OF BRIDGE ANALYSIS###

#Network Comparison Test

#Load data
Irinadata1 <- read.table("C:/Users/lapab/Dropbox/EAT Lab/AED Webinar/clinicaldata.csv", header=TRUE, sep=",", na = "NA")
Irinadata2 <- read.table("C:/Users/lapab/Dropbox/EAT Lab/AED Webinar/nonclinicaldata.csv", header=TRUE, sep=",", na = "NA")

#Omit missing data
newdata1 <- na.omit(Irinadata1)
newdata2 <- na.omit(Irinadata2)

#Estimate networks
mynetwork1 <- estimateNetwork(newdata1, default="EBICglasso")
mynetwork2 <- estimateNetwork(newdata2, default="EBICglasso")

MyNCT <- NCT(mynetwork1, mynetwork2, it=1000, weighted = TRUE, test.edges = FALSE,  edges='ALL')

summary(NCT1)