#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