# load the required packages
library("qgraph")
library("dplyr")
library("glmnet")
library("corpcor")
library("bootnet")
library("glasso")
library("mgm")
library("igraph")
library("ppls")
library("Parker")
library("ggplot2")
library("NetworkComparisonTest")


# load the data
Data1 <- read.csv(file = ".csv")
str(Data1)

group1 <- c(
  rep("Maternal coparenting",4),
  rep("Maternal involvement",4),
  rep("Maternal parenting styles",3),
  rep("Paternal coparenting",4),
  rep("Paternal involvement",4),
  rep("Paternal parenting styles",3))

abbreviation1 <- list("Maternal integrity","Maternal reprimand","Maternal conflict","Maternal disparagement","Maternal living care","Maternal intellectual support","Maternal emotional interaction","Maternal discipline instruction","Maternal rejection","Maternal warmth","Maternal overprotection","Paternal integrity","Paternal reprimand","Paternal conflict","Paternal disparagement","Paternal living care","Paternal intellectual support","Paternal emotional interaction","Paternal discipline instruction","Paternal rejection","Paternal warmth","Paternal overprotection")

# EBIC network
Network <- estimateNetwork(Data1, default = "EBICglasso")
plot(Network, layout = "spring", groups = group1, nodeNames = abbreviation1, legend.cex =.45, vsize=7, minimum=0, maximum = 0.6, cut = 0, edge.width=1, palette = "pastel")
summary(Network)
Network$graph

label1 <- c("MC1:Maternal integrity","MC2:Maternal reprimand","MC3:Maternal conflict","MC4:Maternal disparagement","MI1:Maternal living care","MI2:Maternal intellectual support","MI3:Maternal emotional interaction","MI4:Maternal discipline instruction","MS1:Maternal rejection","MS2:Maternal warmth","MS3:Maternal overprotection","PC1:Paternal integrity","PC2:Paternal reprimand","PC3:Paternal conflict","PC4:Paternal disparagement","PI1:Paternal living care","PI2:Paternal intellectual support","PI3:Paternal emotional interaction","PI4:Paternal discipline instruction","PS1:Paternal rejection","PS2:Paternal warmth","PS3:Paternal overprotection")

# Centrality 
centralityPlot(Network,labels = label1,include = c("Strength","Closeness", "ExpectedInfluence","Betweenness"), decreasing = TRUE)
centrality$node.centrality
nc <- centrality$node.centrality
ebc <- centrality$edge.betweenness.centrality
write.csv(nc, file = ".csv")
write.csv(ebc, file = ".csv")


# accuracy
acc1 <- bootnet(Data, nBoots = 1000,  nCores = 8, default = "EBICglasso",
                 statistics = c("edge","strength","closeness", "betweenness", "expectedInfluence"),
                 verbose = FALSE)
plot(acc1, labels = FALSE, order = "sample") 
plot(acc1, "edge", plot = "difference", onlyNonZero = TRUE, order = "sample")
plot(acc1, "expectedInfluence", plot = "difference")
plot(acc1, "betweenness", plot = "difference")

# Stability
sta1 <- bootnet(Data, nBoots = 1000, default = "EBICglasso", 
                statistics = c ("Betweenness","expectedInfluence"),
                type = "case", nCores = 8)
corStability(sta1)
plot(sta1, c("Betweenness","expectedInfluence"))


##——————————————————————————network comparison——————————————————————————##

# load the data
Data3 <- read.csv(file = ".csv")
str(Data3)

# Split file 
early <- subset(Data3, agg == 1 )
middle <- subset(Data3, agg == 2 )
late <- subset(Data3, agg == 3 )

# Subsetting
All <- subset(Data3, select = c((1:22)))
subset_early <- subset(early, select = c((1:22)))
subset_middle <- subset(middle, select = c((1:22)))
subset_late <- subset(late, select = c((1:22)))

sub_group3 <- c(
  rep("Maternal coparenting",4),
  rep("Maternal involvement",4),
  rep("Maternal parenting styles",3),
  rep("Paternal coparenting",4),
  rep("Paternal involvement",4),
  rep("Paternal parenting styles",3))

abbreviation3 <- list("Maternal integrity","Maternal reprimand","Maternal conflict","Maternal disparagement","Maternal living care","Maternal intellectual support","Maternal emotional interaction","Maternal discipline instruction","Maternal rejection","Maternal warmth","Maternal overprotection","Paternal integrity","Paternal reprimand","Paternal conflict","Paternal disparagement","Paternal living care","Paternal intellectual support","Paternal emotional interaction","Paternal discipline instruction","Paternal rejection","Paternal warmth","Paternal overprotection")


network1 <-estimateNetwork(subset_early, default = "EBICglasso")
network2 <-estimateNetwork(subset_middle, default = "EBICglasso")
network3 <-estimateNetwork(subset_late, default = "EBICglasso")
network1$graph
network2$graph
network3$graph

centrality1 <- centrality_auto(network1)
nc1 <- centrality1$node.centrality
write.csv(nc1, file = ".csv")

centrality2 <- centrality_auto(network2)
nc2 <- centrality2$node.centrality
write.csv(nc2, file = ".csv")

centrality3 <- centrality_auto(network3)
nc3 <- centrality3$node.centrality
write.csv(nc3, file = ".csv")

L <-averageLayout(network1, network2,network3)

plot(network1, groups = sub_group3, layout = L,
     legend.cex =0.45, nodeNames = abbreviation3,
     border.width=1, border.color="black", vsize=6, minimum=0,maximum =0.5, cut = 0,
     palette = "pastel")
summary (network1)

plot(network2, groups = sub_group3, layout = L, 
     legend.cex =0.45, nodeNames = abbreviation3,
     border.width=1, border.color="black", vsize=6, minimum=0,maximum =0.5, cut = 0,
     palette = "pastel")
summary (network2)

plot(network3, groups = sub_group3, layout = L, 
     legend.cex =0.45, nodeNames = abbreviation3,
     border.width=1, border.color="black", vsize=6, minimum=0,maximum =0.5, cut = 0,
     palette = "pastel")
summary (network3)

centralityPlot(list(early = network1, middle = network2, late=network3),
               decreasing = TRUE, include = c("ExpectedInfluence","Betweenness"),labels = label1)

boot1 <- bootnet(network1, nBoots = 1000, 
                 statistics = c ("Betweenness","expectedInfluence"),
                 type = "case", nCores = 8) 
boot2 <- bootnet(network2, nBoots = 1000, 
                 statistics = c ("Betweenness","expectedInfluence"),
                 type = "case", nCores = 8) 
boot3 <- bootnet(network3, nBoots = 1000, 
                 statistics = c ("Betweenness","expectedInfluence"),
                 type = "case", nCores = 8)

corStability(boot1)
corStability(boot2)
corStability(boot3)

plot(boot1, c("Betweenness","expectedInfluence"))
plot(boot2, c("expectedInfluence","Betweenness"))
plot(boot3,  c("expectedInfluence","Betweenness"))

# accuracy
acc1 <- bootnet(subset_early,nBoots = 1000, default = "EBICglasso", nCores = 8, verbose = FALSE)
acc2 <- bootnet(subset_middle,nBoots = 1000, default = "EBICglasso", nCores = 8, verbose = FALSE)
acc3 <- bootnet(subset_late,nBoots = 1000, default = "EBICglasso", nCores = 8, verbose = FALSE)

plot(acc1, labels = FALSE, order = "sample") 
plot(acc2, labels = FALSE, order = "sample") 
plot(acc3, labels = FALSE, order = "sample") 

# comparison
nct_12 <- NCT(subset_early, subset_middle, gamma = 0.5, it = 1000, 
              binary.data=FALSE, weighted = TRUE, progressbar = FALSE,
              test.edges=TRUE, test.centrality=TRUE, 
              centrality=c("betweenness","expectedInfluence"),nodes="all")
nct_13 <- NCT(subset_early, subset_late, gamma = 0.5, it = 1000, 
              binary.data=FALSE, weighted = TRUE, progressbar = FALSE,
              test.edges=TRUE, test.centrality=TRUE, 
              centrality=c("betweenness","expectedInfluence"),nodes="all")
nct_23 <- NCT(subset_middle, subset_late, gamma = 0.5, it = 1000, 
              binary.data=FALSE, weighted = TRUE, progressbar = FALSE,
              test.edges=TRUE, test.centrality=TRUE, 
              centrality=c("betweenness","expectedInfluence"),nodes="all")

nct_12$glstrinv.sep    
nct_12$glstrinv.pval   
nct_12$nwinv.pval      
nct_12$diffcen.pval   
nct_12$einv.pvals     

nct_13$glstrinv.sep   
nct_13$glstrinv.pval   
nct_13$nwinv.pval      
nct_13$diffcen.pval   
nct_13$einv.pvals     

nct_23$glstrinv.sep 
nct_23$glstrinv.pval  
nct_23$nwinv.pval      
nct_23$diffcen.pval   
nct_23$einv.pvals    

plot(nct_12, what="strength")
plot(nct_12, what="edge")
plot(nct_12, what="network")

plot(nct_13, what="strength")
plot(nct_13, what="edge")
plot(nct_13, what="network")

plot(nct_23, what="strength")
plot(nct_23, what="edge")
plot(nct_23, what="network")