library(“readxl”) library(statnet) library(tidyverse) library(qgraph) library(huge) library(NetworkComparisonTest) library(bootnet) library(matrixcalc) library(glasso) library(“psychTools”) library(dplyr) library(EGAnet) library(qgraph) library(MissMech) library(VIM) library(Amelia) library(mice) library(sjmisc) library(tidyverse) library(data.table) library(reshape) library(ggplot2)
#Raw data BS <- read_excel(“D:/Research/REFERENCES/Statistic/Network analysis/PROJECT/DATA-RM.xlsx”,sheet = “257BS raw”) TM <- read_excel(“D:/Research/REFERENCES/Statistic/Network analysis/PROJECT/DATA-RM.xlsx”,sheet = “257TM raw”)
#257-257 BS <- read_excel(“D:/Research/REFERENCES/Statistic/Network analysis/PROJECT/DATA-RM.xlsx”,sheet = “257BS imputed”) TM <- read_excel(“D:/Research/REFERENCES/Statistic/Network analysis/PROJECT/DATA-RM.xlsx”,sheet = “257TM imputed”)
THR=subset(TM,RM==“1”) TLR=subset(TM,RM==“0”) BHR=subset(BS,RM==“1”) BLR=subset(BS,RM==“0”) TM=subset(TM,RM==“1”|RM==“0”) BS=subset(BS,RM==“1”|RM==“0”) names(BS)
setnames(BHR, old = c(“PANSS-P-TO”,“PANSS-N-TO”, “PANSS-G-TO”, “CDSS TO”, “BCSS-NS”,“BCSS-PS”,“BCSS-NO”,“BCSS-PO”,“BS-EM”, “BS-CO”), new = c(“P”,“N”,“G”, “CS”,“NS”,“PS”,“NO”,“PO”,“Em”, “Co”)) setnames(BLR, old = c(“PANSS-P-TO”, “PANSS-N-TO”, “PANSS-G-TO”, “CDSS TO”, “BCSS-NS”,“BCSS-PS”,“BCSS-NO”,“BCSS-PO”,“BS-EM”, “BS-CO”), new = c(“P”,“N”,“G”, “CS”,“NS”,“PS”,“NO”,“PO”,“Em”, “Co”)) setnames(THR, old = c(“PANSS-P-TO”, “PANSS-N-TO”, “PANSS-G-TO”, “CDSS TO”, “BCSS-NS”,“BCSS-PS”,“BCSS-NO”,“BCSS-PO”,“BS-EM”, “BS-CO”), new = c(“P”,“N”,“G”, “CS”,“NS”,“PS”,“NO”,“PO”,“Em”, “Co”)) setnames(TLR, old = c(“PANSS-P-TO”, “PANSS-N-TO”, “PANSS-G-TO”, “CDSS TO”, “BCSS-NS”,“BCSS-PS”,“BCSS-NO”,“BCSS-PO”,“BS-EM”, “BS-CO”), new = c(“P”,“N”,“G”, “CS”,“NS”,“PS”,“NO”,“PO”,“Em”, “Co”))
#Network analysis
BS= BS[, c(2:10,11)] TM= TM[, c(2:10,11)] BHR= BHR[, c(2:10,11)] THR= THR[, c(2:10,11)] BLR= BLR[, c(2:10,11)] TLR= TLR[, c(2:10,11)] names(TM) names(BS)
#MV SUMMARY
BS=sapply(BS,as.numeric) TM=sapply(TM,as.numeric) BLR=sapply(BLR,as.numeric) BHR=sapply(BHR,as.numeric) THR=sapply(THR,as.numeric) TLR=sapply(TLR,as.numeric)
#Delete MV
BS <- na.omit(BS) TM <- na.omit(TM) BHR <- na.omit(BHR) BLR <- na.omit(BLR) THR <- na.omit(THR) TLR <- na.omit(TLR)
#Transformation BS1 <- huge.npn(BS ) TM1<- huge.npn(TM) BLR1 <- huge.npn(BLR ) BHR1<- huge.npn(BHR) THR1 <- huge.npn(THR ) TLR1<- huge.npn(TLR)
BS2 <- cor_auto(BS1) TM2 <- cor_auto(TM1) BLR2 <- cor_auto(BLR1) BHR2 <- cor_auto(BHR1) TLR2 <- cor_auto(TLR1) THR2 <- cor_auto(THR1)
#Estimate
network5 <- estimateNetwork(BHR1,default = “EBICglasso”) network6<- estimateNetwork(THR1,default = “EBICglasso”)
network7 <- estimateNetwork(BLR1,default = “EBICglasso”) network8<- estimateNetwork(TLR1,default = “EBICglasso”)
p1=plot(network5, layout = “spring”,title = “a)”) p2=plot(network6, layout = “spring”,title = “b)”)
p3=plot(network7, layout = “spring”,title = “c)”) p4=plot(network8, layout = “spring”,title = “d)”)
#strength
Cent2 <- centralityTable(network5, standardized = TRUE) Cent3 <- centralityTable(network6, standardized = TRUE) table.data1 <- data.frame(“Item” = Cent2[Cent2$measure == “Strength”,3], “Strength” = Cent2[Cent2$measure == “Strength”,5]) table.data2 <- data.frame(“Item” = Cent3[Cent3$measure == “Strength”,3], “Strength” = Cent3[Cent3$measure == “Strength”,5]) data3 <- table.data1[,2] data4 <- table.data2[,2] data5=data4-data3 data5
p5=plot(network6, layout = “spring”,color=c(“olivedrab2”, “darkseagreen1”, “wheat3”, “goldenrod2”, “darkorange1”, “darkslategray2”, “mediumpurple”),vsize =data5*10,title = “a)”)
network9 <- qgraph(TLR2, layout = “spring”, graph = “glasso”, tuning = 0.5, sampleSize = nrow(TLR), legend.cex = 0.3,groups = group.item, nodeNames=Names1a,labels = Labels1a,title = “6M follow up”, color=c(“olivedrab2”, “darkseagreen1”, “wheat3”, “goldenrod2”, “darkorange1”, “darkslategray2”, “mediumpurple”),cut = 0.03, maximum = 1,details = TRUE,
borders=FALSE, theme=“colorblind”, usePCH=TRUE,edge.labels=TRUE,esize=14,vsize =data5*10)
table.data1\(Item <- factor(table.data1\)Item,levels = c(“P”,“N”,“G”, “CS”,“NS”,“PS”,“NO”,“PO”,“Em”, “Co”)) centralityPlot(list(‘R-B’ = network5,‘R-6’ = network6),include=c(“Strength”,“Closeness”, “Betweenness”))+ theme(legend.title = element_blank())
centralityPlot(list(‘NR-B’ = network7,‘NR-6’ = network8),include=c(“Strength”,“Closeness”, “Betweenness”))+ theme(legend.title = element_blank())
centralityPlot(…, labels, scale = c(“z-scores”, “raw”, “raw0”,“relative”), include = c(“Degree”,“Strength”,“OutDegree”,“InDegree”,“OutStrength”, “InStrength”, “Closeness”,“Betweenness”), theme_bw = TRUE, print = TRUE, verbose = TRUE, standardized, relative, weighted = TRUE, signed = TRUE, orderBy = “default”, decreasing = FALSE)
p1=ggplot(CEN21, aes(x=Node, y=value, group=variable, color= variable)) + theme_classic()+theme(legend.position = “top”) + theme(legend.title = element_blank())+ geom_line() + geom_point()+xlab("“)+scale_x_discrete(name=”“)+ylab(”Strength“)+ theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) p2=ggplot(CEN22, aes(x=Node, y=value, group=variable, color= variable)) + theme_classic()+theme(legend.position =”none“)+ geom_line() + geom_point()+xlab(”“)+scale_x_discrete(name=”“)+ylab(”Closeness“)+ theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) p3=ggplot(CEN21, aes(x=Node, y=value, group=variable, color= variable)) + theme_classic()+theme(legend.position =”none“)+ geom_line() + geom_point()+xlab(”“)+scale_x_discrete(name=” “)+ylab(”Betweennes")
gridExtra::grid.arrange(p1, p2,p3, nrow =3)
library(“igraph”) igraph_object <- as.igraph(network5)
wc <- walktrap.community(igraph_object) # This is looking for structure in your graph modularity(wc) # The modularity score of that structure
nct.res_long <- NCT(network5, network6, it=1000,paired = TRUE,test.edges = TRUE,edges=“all”, test.centrality=TRUE,centrality=c(“strength”,“closeness”,“betweenness”))
#strength
Cent2 <- centralityTable(network5, standardized = TRUE) Cent3 <- centralityTable(network6, standardized = TRUE) table.data1 <- data.frame(“Item” = Cent2[Cent2$measure == “Strength”,3], “Strength” = Cent2[Cent2$measure == “Strength”,5]) table.data2 <- data.frame(“Item” = Cent3[Cent3$measure == “Strength”,3], “Strength” = Cent3[Cent3$measure == “Strength”,5]) data3 <- table.data1[,2] data4 <- table.data2[,2] data5=data4-data3
Cent2 <- centralityTable(network5, standardized = TRUE) Cent3 <- centralityTable(network6, standardized = TRUE) table.data1 <- data.frame(“Item” = Cent2[Cent2$measure == “Closeness”,3], “Closeness” = Cent2[Cent2$measure == “Closeness”,5]) table.data2 <- data.frame(“Item” = Cent3[Cent3$measure == “Closeness”,3], “Closeness” = Cent3[Cent3$measure == “Closeness”,5]) data3 <- table.data1[,2] data4 <- table.data2[,2] data5=data4-data3
Cent2 <- centralityTable(network5, standardized = TRUE) Cent3 <- centralityTable(network6, standardized = TRUE) table.data1 <- data.frame(“Item” = Cent2[Cent2$measure == “Betweenness”,3], “Betweenness” = Cent2[Cent2$measure == “Betweenness”,5]) table.data2 <- data.frame(“Item” = Cent3[Cent3$measure == “Betweenness”,3], “Betweenness” = Cent3[Cent3$measure == “Betweenness”,5]) data3 <- table.data1[,2] data4 <- table.data2[,2] data5=data4-data3
##Function estimator <- function(Data){ library(“qgraph”) library(“huge”) Network_cor <- huge.npn(Data) Network_cor<- cor_auto(Network_cor) Network <- qgraph(Network_cor, graph=‘glasso’, sampleSize=nrow(Data), layout=‘spring’) return(getWmat(Network))}
set.seed(1) net_boot <- estimateNetwork(BS, fun = estimator) net_boot2 <- estimateNetwork(TM, fun = estimator) net_boot3 <- estimateNetwork(BHR, fun = estimator) net_boot4 <- estimateNetwork(BLR, fun = estimator) net_boot5 <- estimateNetwork(THR, fun = estimator) net_boot6 <- estimateNetwork(TLR, fun = estimator)
stability_men5 <- bootnet(net_boot3, nBoots = 1000, nCores = 8,nodeNames=Names1a,labels = Labels1a) stability_men7 <- bootnet(net_boot4, nBoots = 1000, nCores = 8,nodeNames=Names1a,labels = Labels1a) stability_men9 <- bootnet(net_boot5, nBoots = 1000, nCores = 8,nodeNames=Names1a,labels = Labels1a) stability_men11 <- bootnet(net_boot6, nBoots = 1000, nCores = 8,nodeNames=Names1a,labels = Labels1a)
plot(stability_men5, order=“sample”, plot=“area”, prop0=TRUE) plot(stability_men7, order=“sample”, plot=“area”, prop0=TRUE) plot(stability_men9, order=“sample”, plot=“area”, prop0=TRUE) plot(stability_men11, order=“sample”, plot=“area”, prop0=TRUE)
#different test plot(stability_men5, “edge”, plot = “difference”, onlyNonZero = TRUE, order = “sample”,layout=layout\(A) # difference of edges plot(stability_men7, "edge", plot = "difference", onlyNonZero = TRUE, order = "sample",layout=layout\)A) # difference of edges plot(stability_men9, “edge”, plot = “difference”, onlyNonZero = TRUE, order = “sample”,layout=layout\(A) # difference of edges plot(stability_men11, "edge", plot = "difference", onlyNonZero = TRUE, order = "sample",layout=layout\)A) # difference of edges
plot(stability_men5, “strength”, plot = “difference”) # difference of edges plot(stability_men7, “strength”, plot = “difference”) # difference of edges plot(stability_men9, “strength”, plot = “difference”) # difference of edges plot(stability_men11, “strength”, plot = “difference”) # difference of edges
stability_men6 <- bootnet(net_boot3,statistics=c(‘strength’,‘closeness’,‘betweenness’), nBoots = 1000,type = “case”, nCores = 8) # 3mins stability_men8 <- bootnet(net_boot4,statistics=c(‘strength’,‘closeness’,‘betweenness’), nBoots = 1000,type = “case”, nCores = 8) # 3mins stability_men10 <- bootnet(net_boot5,nBoots = 1000,type = “case”, nCores = 8, statistics=c(‘strength’,‘closeness’,‘betweenness’,‘edge’)) # 3mins stability_men12 <- bootnet(net_boot6,statistics=c(‘strength’,‘closeness’,‘betweenness’,‘edge’), nBoots = 1000,type = “case”, nCores = 8) # 3mins
corStability(stability_men6) corStability(stability_men8) corStability(stability_men10) corStability(stability_men12)
plot(stability_men6,statistics=c(‘strength’,‘closeness’,‘betweenness’))+ theme(legend.title = element_blank()) plot(stability_men8,statistics=c(‘strength’,‘closeness’,‘betweenness’))+ theme(legend.title = element_blank()) plot(stability_men10,statistics=c(‘strength’,‘closeness’,‘betweenness’))+ theme(legend.title = element_blank()) plot(stability_men12,statistics=c(‘strength’,‘closeness’,‘betweenness’))+ theme(legend.title = element_blank())