#Network Analysis Script
#This script performs the network analysis for the manuscript Bower et al. (under review).,
#"A Network Analysis of Children’s Emerging Place Value Concepts"
#For any questions or comments on this script, please contact Corinne Bower, cbower@umd.edu

#Load Packages
library("haven")
library("ggplot2")
library("bootnet")
library("qgraph")
library("igraph")
library("NetworkComparisonTest")

#read data
df <- read_sav("~/Box Sync/PVL_Study 1 and 3/T1-T2-T3 Study1/Data/Combined PVL Longitudinal_rev11Sept2019PW_converted to SPSS 16Oct2019_FixedNLEerror_AddedB10vars_12Mar2021_CB.sav")
View(df)

#Create Time 1 (kindergarten) and Time 2 (first grade) data frames.
df.PVt1.noMR <- df[c("B10_count_total_t1", "tc_total_t1", "en_total_t1", "mc_total_t1", "wn_a_total_t1", "nle_avg_t1r")]
colnames(df.PVt1.noMR) <- c("B10", "RW", "EN", "MC", "CU", "NLE")
df.PVt2.noMR <- df[c("B10_total_t2", "tc_total_t2", "en_total_t2", "mc_total_t2", "wn_a_total_t2", "nle_avg_t2r")]
colnames(df.PVt2.noMR) <- c("B10", "RW", "EN", "MC", "CU", "NLE")

######TIME 1 Network##########
#Network Estimation
Network.t1.noMR <- estimateNetwork(df.PVt1.noMR, default = "EBICglasso")
#quick plotting of network (not for publication)
plot(Network.t1.noMR, layout = "spring", labels = colnames(df.PVt1.noMR))
qgraph(Network.t1.noMR$graph, layout="spring")
#For varying thresholds, alter threshold level (manuscript levels= 0, .2, .4)
qgraph(Network.t1.noMR$graph, layout="circle", threshold=.2)

#Centrality Indices
centralityPlot(Network.t1.noMR, include =c("Betweenness","Closeness","Strength"))
#Plot centrality indices as bar graphs instead of default graph
library("magrittr")
library("dplyr")
T1centralityZ <- centralityTable(Network.t1.noMR$graph, standardized = TRUE,  relative = FALSE, weighted =
                                   TRUE, signed = TRUE)
T1centrality.Zdf <- as.data.frame(T1centralityZ)
T1centrality.Zdf <- T1centrality.Zdf[,c("node", "measure", "value")]
T1centrality.Zdf <- T1centrality.Zdf %>% mutate(node = factor(node, 
                                                              levels = c("B10", "EN", "CU", "RW", "MC", "NLE"),
                                                              ordered = TRUE))
#Strength Plot
ggplot(data=T1centrality.Zdf[T1centrality.Zdf$measure=="Strength",], aes(x=node, y=value)) +
  geom_bar(stat="identity") + theme_classic() + labs(title="Strength", x="Node", y = "Z-scores") + geom_hline(yintercept=0) +
  theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold")) + 
  scale_y_continuous(breaks=seq(-2,2,.5),limits = c(-2,2))  
#Closeness Plot
ggplot(data=T1centrality.Zdf[T1centrality.Zdf$measure=="Closeness",], aes(x=node, y=value)) +
  geom_bar(stat="identity") + theme_classic() + labs(title="Closeness", x="Node", y = "Z-scores") + geom_hline(yintercept=0) +
  theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold")) + 
  scale_y_continuous(breaks=seq(-2,2,.5),limits = c(-2,2))  
#Betweenness Plot
ggplot(data=T1centrality.Zdf[T1centrality.Zdf$measure=="Betweenness",], aes(x=node, y=value)) +
  geom_bar(stat="identity") + theme_classic() + labs(title="Betweenness", x="Node", y = "Z-scores") + geom_hline(yintercept=0) +
  theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold")) + 
  scale_y_continuous(breaks=seq(-2,2,.5),limits = c(-2,2)) 

#Stability of Centrality Indices
CentralStability.t1.noMR <- bootnet(Network.t1.noMR, nBoots = 1000, type = "case", statistics = c("edge", "strength", "closeness", "betweenness"))
#Produce CS coefficients for each index
corStability(CentralStability.t1.noMR)
summary(CentralStability.t1.noMR)
#Plot stabilities
plot(CentralStability.t1.noMR, statistics = c("strength", "closeness", "betweenness"))

#Edge-weight accuracy
EdgeWgt.t1.noMR<- bootnet(Network.t1.noMR, nBoots = 2500, statistics = c("edge", "strength", "closeness", "betweenness"))
df.EdgeWgt.t1.noMR <- as.data.frame(summary(EdgeWgt.t1.noMR))
print(EdgeWgt.t1.noMR)
#Number of non-zero edges in sample: 15/15; Mean weight of sample: 0.1548166 
summary(EdgeWgt.t1.noMR)
plot(EdgeWgt.t1.noMR, labels = TRUE, order = "sample")

#Testing for significance differences in edge weights
#STRENGTH difference between nodes of B10 and RW.
differenceTest(EdgeWgt.t1.noMR, 1, 2, "strength") #(CI: -.388, .206)
#matrix plot of strength diff tests of nodes:
plot(EdgeWgt.t1.noMR, "strength")
plot(EdgeWgt.t1.noMR, "closeness")
#matrix plot of strength diff tests of non-zero edges:
plot(EdgeWgt.t1.noMR, "edge", plot = "difference", onlyNonZero = TRUE, order = "sample")
#difference between edge B10-RW and B10-EN
differenceTest(EdgeWgt.t1.noMR, "B10--RW", "B10--EN", "edge")
#NODE CLOSENESS.
differenceTest(EdgeWgt.t1.noMR, "B10", "RW", measure = "closeness")

#Community Detection using spinglass algorithm
graph1<-qgraph(Network.t1.noMR$graph, layout="spring", vsize=7, cut=0, maximum=.45, border.width=1.5)
g1 = as.igraph(graph1, attributes=TRUE)
sgc1 <- spinglass.community(g1)
sgc1$membership
sgc1$vcount
#Result: 2 clusters (RW, MC, NLE; B10, EN, CU)
#Color graph based on these findings:
group.spinglass1<- list(c(1,3,5), c(2,4,6))
qgraph(Network.t1.noMR$graph, layout="spring", vsize=9, cut=0, maximum=.45, border.width=1.5, groups=group.spinglass1, palette = 'ggplot2')
#Visualize using network with MDS qgraph (layout was created below first)
qgraph(Network.t1.noMR$graph, layout=Time1.noMR_MDS_interval$conf, groups=group.spinglass1, palette = 'ggplot2')
#At different thresholds with changing border and fill colors based on hypotheses/results
qgraph(Network.t1.noMR$graph, layout=Time1.noMR_MDS_interval$conf, vsize=7,
       threshold=0, color=c("lightpink","turquoise"),
       border.color=c("red4","turquoise4"), border.width=6)
#Run several simulations with spinglass
matrix_spinglass1 <- matrix(NA, nrow=1,ncol=1000)
for (i in 1:1000) {
  set.seed(i)
  spinglass1 <- spinglass.community(g1)
  matrix_spinglass1[1,i] <- max(spinglass1$membership)
}
mean(as.vector(matrix_spinglass1)) #2
max(as.vector(matrix_spinglass1)) #3
min(as.vector(matrix_spinglass1)) #2
median(as.vector(matrix_spinglass1)) #2

#Visualizing network using Jones' (2018) multidimensional scaling function
#plot edges using partial-correlations from graphical LASSO network in a correlational space.
library("networktools")
#reconfigure similarity matrix (correlations) to dissimilarity matrix
library("smacof")
Time1corr.noMR <- cor(df.PVt1.noMR, use = "complete.obs")
qgraph(Time1corr.noMR, layout="spring")
dissimilarity_Time1.noMR <- sim2diss(Time1corr.noMR)
Time1.noMR_MDS <- mds(dissimilarity_Time1.noMR)
plot(Time1.noMR_MDS, plot.type = "Shepard")
#Test to see which transformation is best:
#ordinal
Time1.noMR_MDS_ordinal <- mds(dissimilarity_Time1.noMR,
                              type="ordinal")
plot(Time1.noMR_MDS_ordinal, plot.type = "Shepard",
     main="Ordinal")
text(1.1,0.3, paste("Stress =",
                    round(Time1.noMR_MDS_ordinal$stress,2)))
Time1.noMR_MDS_ordinal$stress #.00
#ratio
Time1.noMR_MDS_ratio <- mds(dissimilarity_Time1.noMR,
                            type="ratio")
plot(Time1.noMR_MDS_ratio, plot.type = "Shepard",
     main="Ratio")
text(1.1,0.3, paste("Stress =",
                    round(Time1.noMR_MDS_ratio$stress,2)))
Time1.noMR_MDS_ratio$stress #.22
#interval
Time1.noMR_MDS_interval <- mds(dissimilarity_Time1.noMR,
                               type="interval")
plot(Time1.noMR_MDS_interval, plot.type = "Shepard",
     main="Interval")
text(1.1,0.3, paste("Stress =",
                    round(Time1.noMR_MDS_interval$stress,2)))
Time1.noMR_MDS_interval$stress #.09
#mspline
Time1.noMR_MDS_mspline <- mds(dissimilarity_Time1.noMR,
                              type="mspline")
plot(Time1.noMR_MDS_mspline, plot.type = "Shepard",
     main="Spline")
text(1.1,0.3, paste("Stress =",
                    round(Time1.noMR_MDS_mspline$stress,2)))
Time1.noMR_MDS_mspline$stress #.04
#INTERVAL is best b/c lowish stress and distances match the edge weights (not the case in ordinal)
#+T1 qgraph using LASSO network with correlational space
qgraph(Network.t1.noMR$graph, layout=Time1.noMR_MDS_interval$conf, threshold=0)
text(-1,-1, paste("Stress = ",
                  round(Time1.noMR_MDS_interval$stress,2)))




######TIME 2 Network##########
#Network Estimation
Network.t2.noMR <- estimateNetwork(df.PVt2.noMR, default = "EBICglasso")
#quick plotting of network (not for paper)
plot(Network.t2.noMR, layout = "spring", labels = colnames(df.PVt2.noMR))
qgraph(Network.t2.noMR$graph, layout="spring")
#For varying thresholds, alter threshold level (manuscript levels= 0, .2, .4)
qgraph(Network.t2.noMR$graph, layout="circle", threshold=.2)

#Centrality Indices
centralityPlot(Network.t2.noMR, include =c("Betweenness","Closeness","Strength"))
#Plot centrality indices as bar graphs instead of default graph
T2centralityZ <- centralityTable(Network.t2.noMR$graph, standardized = TRUE,  relative = FALSE, weighted =
                                   TRUE, signed = TRUE)
T2centrality.Zdf <- as.data.frame(T2centralityZ)
T2centrality.Zdf <- T2centrality.Zdf[,c("node", "measure", "value")]
T2centrality.Zdf <- T2centrality.Zdf %>% mutate(node = factor(node, 
                                                              levels = c("B10", "EN", "CU", "RW", "MC", "NLE"),
                                                              ordered = TRUE))
#Strength Plot
ggplot(data=T2centrality.Zdf[T2centrality.Zdf$measure=="Strength",], aes(x=node, y=value)) +
  geom_bar(stat="identity") + theme_classic() + labs(title="Strength", x="Node", y = "Z-scores") + geom_hline(yintercept=0) +
  theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold")) + 
  scale_y_continuous(breaks=seq(-2,2,.5),limits = c(-2,2))  
#Closeness Plot
ggplot(data=T2centrality.Zdf[T2centrality.Zdf$measure=="Closeness",], aes(x=node, y=value)) +
  geom_bar(stat="identity") + theme_classic() + labs(title="Closeness", x="Node", y = "Z-scores") + geom_hline(yintercept=0) +
  theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold")) + 
  scale_y_continuous(breaks=seq(-2,2,.5),limits = c(-2,2))  
#Betweenness Plot
ggplot(data=T2centrality.Zdf[T2centrality.Zdf$measure=="Betweenness",], aes(x=node, y=value)) +
  geom_bar(stat="identity") + theme_classic() + labs(title="Betweenness", x="Node", y = "Z-scores") + geom_hline(yintercept=0) +
  theme(axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold")) + 
  scale_y_continuous(breaks=seq(-2,2,.5),limits = c(-2,2)) 

#Stability of Centrality Indices
CentralStability.t2.noMR <- bootnet(Network.t2.noMR, nBoots = 1000, type = "case", statistics = c("edge", "strength", "closeness", "betweenness"))
#Produce CS coefficients for each index
corStability(CentralStability.t2.noMR)
summary(CentralStability.t2.noMR)
#Plot stabilities
plot(CentralStability.t2.noMR, statistics = c("strength", "closeness", "betweenness"))

#Edge-weight accuracy
EdgeWgt.t2.noMR<- bootnet(Network.t2.noMR, nBoots = 2500, statistics = c("edge", "strength", "closeness", "betweenness"))
df.EdgeWgt.t2.noMR <- as.data.frame(summary(EdgeWgt.t2.noMR))
print(EdgeWgt.t2.noMR)
summary(EdgeWgt.t2.noMR)
plot(EdgeWgt.t2.noMR, labels = TRUE, order = "sample")

#Testing for significance differences in edge weights
#STRENGTH difference between nodes of B10 and RW.
differenceTest(EdgeWgt.t2.noMR, 1, 2, "strength") #(CI: -.17, .40)
#Strength difference between nodes RW and MC.
differenceTest(EdgeWgt.t2, 2, 4, "strength")
#Strength difference between nodes B10 and MC.
differenceTest(EdgeWgt.t2, 1, 4, "strength")
#Strength difference between nodes RW and EN.
differenceTest(EdgeWgt.t2, 2, 3, "strength")
#matrix plot of strength diff tests of nodes:
plot(EdgeWgt.t2.noMR, "strength")
plot(EdgeWgt.t2.noMR, "closeness")
#matrix plot of strength diff tests of non-zero edges:
plot(EdgeWgt.t2.noMR, "edge", plot = "difference", onlyNonZero = TRUE, order = "sample")
#difference between edge B10-RW and B10-EN
differenceTest(EdgeWgt.t2.noMR, "B10--RW", "B10--EN", "edge") #CI: -.30, .10
#NODE STRENGTH.
differenceTest(EdgeWgt.t2, "RW", "EN", "strength")
differenceTest(EdgeWgt.t2, "RW", "NLE", "strength")
differenceTest(EdgeWgt.t2, "RW", "CU", "strength")
differenceTest(EdgeWgt.t2, "RW", "MC", "strength")

#Community Detection using spinglass algorithm
graph2<-qgraph(Network.t2.noMR$graph, layout="spring", vsize=7, cut=0, maximum=.45, border.width=1.5)
g2 = as.igraph(graph2, attributes=TRUE)
sgc2 <- spinglass.community(g2)
sgc2$membership
sgc2$vcount
#Result: 2 clusters (B10, RW, EN, CU; MC & NLE)
#Color graph based on these findings:
group.spinglass2<- list(c(1,2,3,5), c(4,6))
qgraph(Network.t2.noMR$graph, layout="spring", vsize=9, cut=0, maximum=.45, border.width=1.5, groups=group.spinglass1, palette = 'ggplot2')
#Visualize using network with MDS qgraph (create layout first, see below)
qgraph(Network.t2.noMR$graph, layout=Time2.noMR_MDS_ordinal$conf, groups=group.spinglass2, palette = 'ggplot2')
#At different thresholds with changing border and fill colors based on hypotheses/results
qgraph(Network.t2.noMR$graph, layout=Time2.noMR_MDS_ordinal$conf, vsize=7,
       threshold=0, groups=group.spinglass2, color=c("lightpink","turquoise"),
       border.color=c("black"), border.width=6)
#Run several simulations with spinglass
matrix_spinglass2 <- matrix(NA, nrow=1,ncol=1000)
for (i in 1:1000) {
  set.seed(i)
  spinglass2 <- spinglass.community(g2)
  matrix_spinglass2[1,i] <- max(spinglass2$membership)
}
mean(as.vector(matrix_spinglass2)) #2
max(as.vector(matrix_spinglass2)) #3
min(as.vector(matrix_spinglass2)) #2
median(as.vector(matrix_spinglass2)) #2


#Visualizing network using Jones' (2018) multidimensional scaling function
#plot edges using partial-correlations from graphical LASSO network in a correlational space.
#reconfigure similarity matrix (correlations) to dissimilarity matrix
Time2corr.noMR <- cor(df.PVt2.noMR, use = "complete.obs")
qgraph(Time2corr.noMR, layout="spring")
dissimilarity_Time2.noMR <- sim2diss(Time2corr.noMR)
Time2.noMR_MDS <- mds(dissimilarity_Time2.noMR)
#Test to see which transformation is best:
#ordinal transformation
Time2.noMR_MDS_ordinal <- mds(dissimilarity_Time2.noMR,
                              type="ordinal")
plot(Time2.noMR_MDS_ordinal, plot.type = "Shepard",
     main="Ordinal")
text(1.1,0.3, paste("Stress =",
                    round(Time2.noMR_MDS_ordinal$stress,2)))
Time2.noMR_MDS_ordinal$stress #.02
#ratio transformation
Time2.noMR_MDS_ratio <- mds(dissimilarity_Time2.noMR,
                            type="ratio")
plot(Time2.noMR_MDS_ratio, plot.type = "Shepard",
     main="Ratio")
text(1.1,0.3, paste("Stress =",
                    round(Time2.noMR_MDS_ratio$stress,2)))
Time2.noMR_MDS_ratio$stress #.23
#interval transformation
Time2.noMR_MDS_interval <- mds(dissimilarity_Time2.noMR,
                               type="interval")
plot(Time2.noMR_MDS_interval, plot.type = "Shepard",
     main="Interval")
text(1.1,0.3, paste("Stress =",
                    round(Time2.noMR_MDS_interval$stress,2)))
Time2.noMR_MDS_interval$stress #.08
#mspline transformation
Time2.noMR_MDS_mspline <- mds(dissimilarity_Time2.noMR,
                              type="mspline")
plot(Time2.noMR_MDS_mspline, plot.type = "Shepard",
     main="Spline")
text(1.1,0.3, paste("Stress =",
                    round(Time2.noMR_MDS_mspline$stress,2)))
Time2.noMR_MDS_mspline$stress #.07
#ordinal transformation is best b/c of lowest stress value
#+T2 qgraph using LASSO network with correlational space, with varying threshold levels
qgraph(Network.t2.noMR$graph, layout=Time2.noMR_MDS_ordinal$conf, threshold=0)
text(-1,-1, paste("Stress = ",
                  round(Time2.noMR_MDS_ordinal$stress,2)))


###Network Comparison Test####
df.NCT <- df[c("B10_count_total_t1", "tc_total_t1", "en_total_t1", "mc_total_t1", "wn_a_total_t1", "nle_avg_t1r", 
                  "B10_total_t2", "tc_total_t2", "en_total_t2", "mc_total_t2", "wn_a_total_t2", "nle_avg_t2r")]
colnames(df.NCT) <- c("B101", "RW1", "EN1", "MC1", "CU1", "NLE1", "B102", "RW2", "EN2", "MC2", "CU2", "NLE2")
df.NCT <- na.omit(df.NCT) #n=224
data1 <- df.NCT[c("B101", "RW1", "EN1", "MC1", "CU1", "NLE1")]
data2 <- df.NCT[c("B102", "RW2", "EN2", "MC2", "CU2", "NLE2")]
#rename var names to be the same between data1 and data2
colnames(data1) <- c("B10", "RW", "EN", "MC", "CU", "NLE")
colnames(data2) <- c("B10", "RW", "EN", "MC", "CU", "NLE")
Res_2 <- NCT(data1, data2, it=1000, make.positive.definite = TRUE, binary.data=FALSE, 
             paired=TRUE, weighted=TRUE, test.edges=TRUE, edges="all", test.centrality=TRUE, 
             centrality=c("strength","betweenness", "closeness"),
             nodes="all", progressbar=TRUE)
Res_2$glstrinv.pval
Res_2$nwinv.pval