References are linked throughout the text
This script reproduces key analyses from “Networks of Adversity in Childhood and Adolescence and their Relationship to Adult Mental Health”. Please cite the paper if you are re-using some of the code below.
Here, we clean and visualize data from the Avon Longitudinal Study of Parents and Children (ALSPAC) and investigate adverse experiences in childhood and adolescence and its effects on mental health in early adulthood using network analysis. Data to reproduce the results can be obtained upon request from the ALSPAC.
We conducted two network analyses to compare the interrelations of childhood and adolescent adversities as well as their relationships to early adulthood mental health and well-being.
Packages used in this code include: knitr, tidyverse, qgraph, bootnet, lavaan, glasso, dplyr, EGAnet and hfgolino/EGA, NetworkComparisonTest, networktools
Please refer to this tutorial by E. Fried on clustering methods in network analysis and this paper by Golino & Epskamp (2017) specifically on the EGA algorithm.
#clustering
EGAClaACE<-EGA(FinalData[1:10], plot.EGA = TRUE, corr = c("spearman"), model = c("glasso"), algorithm = c("walktrap"), plot.type = c("qgraph"))
EGAClaACE$dim.variables
## items dimension
## 6 6 1
## 7 7 1
## 8 8 1
## 9 9 1
## 10 10 1
## 1 1 2
## 2 2 2
## 3 3 2
## 4 4 2
## 5 5 2
Bootnet allows to estimate and bootstrap networks based on the Gaussian graphical model and visualized in qgraph. Please refer to Epskamp et al. (2018). To visualize the network, we are using the colorblind-friendly colormap “viridis”.
#adaption of ACE data
ACEData <- dplyr::select(FinalData, 1:10,33:39)
colnames(ACEData) <- c("1","2","3","4","5","6","7","8","9", "10","11","12","13","14","15","16","17")
Groups <- c(rep("1. Cluster - Direkt Abuse",5),rep("2. Cluster - Family Factors",5), rep("Mental Health",5), rep("Wellbeing Factors",2))
nodeNamesACEB <- c(nodeNames <- c("Physical Abuse inside of Family","Physical Abuse outside of Family","Mental Abuse inside of Family", "Mental Abuse outside of Family", "Sexual Abuse", "Substance Abuse by Mother","Substance Abuse by Mother's Partner","Abuse of Mother", "Criminality of Carer", "Psychopathology of Carer", "Bipolar Disorder", "Depression", "Chronic Fatigue Syndrom","Alcohol Use Disorder","Drug Use", "Happiness","Life Satisfaction"))
#stepwise model estimation
NetworkACEMH<- estimateNetwork(ACEData, default = "ggmModSelect", corMethod = "spearman")
plot(NetworkACEMH, groups = Groups, nodeNames = nodeNamesACEB, legend.cex = 0.3, layoutScale = c(0.9,0.9), posCol = "gray25", negCol = "grey", negDashed = TRUE, fade = FALSE, esize = 10, color = c("#440154FF", "#21908CFF","#3B528BFF", "#5DC863FF", "#9FDA3AFF"), border.color = "white",border.width = 2, label.color = "white", vsize = 7,curve = 0.1, curveAll = T)
centralityPlot(NetworkACEMH, include = c("Strength", "Closeness","Betweenness","ExpectedInfluence"))
We calculate bridge nodes/centrality using bridge strength. Bridge strength refers to the sum of the absolute value of all edges between a node of a cluster to all the nodes of the opposing cluster. Please refer to networktools, Jones et al. (2019) and Vanzulah (2017).
ACEedges <- getWmat(NetworkACEMH)
write.csv(ACEedges, "ACEedges.csv")
ACEedgeplot <- plot(NetworkACEMH, groups = Groups, nodeNames = nodeNamesACEB, legend.cex = 0.3, layoutScale = c(0.8,0.8), posCol = "gray25", negCol = "grey", negDashed = TRUE, fade = FALSE, esize = 10, color = c("#440154FF", "#3B528BFF", "#21908CFF", "#5DC863FF", "#9FDA3AFF"), border.color = "white",border.width = 2, label.color = "white", vsize = 7,curve = 0.1, curveAll = T)
ACEbridge <- bridge(ACEedgeplot, communities = c('1','1','1','1','1','1','1','1','1','1','2','2','2','2','2','2','2'), useCommunities = "all", directed = NULL, nodes = NULL)
pdf("bridgecentralityACE.pdf", width = 5)
plot(ACEbridge, include = "Bridge Strength")
dev.off()
#Bridge highlighted
Groups5 <- c(rep("ACEs Cluster",2),
rep("1. Bridge Nodes of ACEs Cluster:
3: Mental Abuse inside of Family,
4: Mental Abuse outside of Family,
7: Substance Abuse by Mother's Partner", 2),
rep("ACEs Cluster",2),
rep("1. Bridge Nodes of ACEs Cluster:
3: Mental Abuse inside of Family,
4: Mental Abuse outside of Family,
7: Substance Abuse by Mother's Partner",1),
rep("ACEs Cluster",3),
rep("Mental Health and Wellbeing Cluster",1),
rep("2. Bridge Nodes of Mental Health Cluster:
12: Depression,
15: Drug Use,
17: Life Satisfaction",1),
rep("Mental Health and Wellbeing Cluster",2),
rep("2. Bridge Nodes of Mental Health Cluster:
12: Depression,
15: Drug Use,
17: Life Satisfaction",1),
rep("Mental Health and Wellbeing Cluster",1),
rep("2. Bridge Nodes of Mental Health Cluster:
12: Depression,
15: Drug Use,
17: Life Satisfaction",1))
NetworkACEMH6<- estimateNetwork(ACEData, default = "ggmModSelect", corMethod = "spearman")
plot(NetworkACEMH6, legend.cex = 0.25, groups = Groups5, layoutScale = c(0.9,0.9), posCol = "gray25", negCol = "grey", negDashed = TRUE, fade = FALSE, esize = 10, color = c("#21908CFF","#5DC863FF","snow4","snow3"),border.color = "white",border.width = 2, label.color = "white", vsize = 7,curve = 0.1, curveAll = T)
Bootstrapping & structure stability:
Bootstrapping: Repeatedly constructing CIs, show accuracy of edge-weight estimates & compare edges to one-another.
Correlation stability coefficient (CS-coefficient): Quantify the stability of centrality indices using subset bootstraps.
For more information, please refer to Epskamp et al. (2018).
#Please remove "#" to run this analysis
#BNetworkACEs<- bootnet(NetworkACEMH, default = "ggmModSelect", nBoots = 2000, nCores = 8)
#save(BNetworkACEs, file = "BNetworkACEs.Rdata")
load("BNetworkACEs.Rdata")
plot(BNetworkACEs, order = "sample", labels = FALSE)
#Please remove "#" to run this analysis
#B2NetworkACEs<- bootnet(NetworkACEMH, default = "ggmModSelect", type = ("case"), nBoots = 2000, nCores = 8, statistics = c("edge", "strength", "betweenness", "closeness","expectedInfluence"))
#save(B2NetworkACEs, file = "B2NetworkACEs.Rdata")
load("B2NetworkACEs.Rdata")
plot(B2NetworkACEs, statistics = c("strength", "closeness","betweenness", "expectedInfluence"))
corStabNetworkACEs<- corStability(B2NetworkACEs, cor = 0.7)
## === Correlation Stability Analysis ===
##
## Sampling levels tested:
## nPerson Drop% n
## 1 3911 75.0 197
## 2 5127 67.2 195
## 3 6344 59.4 194
## 4 7561 51.7 186
## 5 8777 43.9 185
## 6 9994 36.1 208
## 7 11211 28.3 210
## 8 12427 20.6 195
## 9 13644 12.8 202
## 10 14861 5.0 228
##
## Maximum drop proportions to retain correlation of 0.7 in at least 95% of the samples:
##
## betweenness: 0.594
## - For more accuracy, run bootnet(..., caseMin = 0.517, caseMax = 0.672)
##
## closeness: 0.206
## - For more accuracy, run bootnet(..., caseMin = 0.128, caseMax = 0.283)
##
## edge: 0.75 (CS-coefficient is highest level tested)
## - For more accuracy, run bootnet(..., caseMin = 0.672, caseMax = 1)
##
## expectedInfluence: 0.75 (CS-coefficient is highest level tested)
## - For more accuracy, run bootnet(..., caseMin = 0.672, caseMax = 1)
##
## strength: 0.75 (CS-coefficient is highest level tested)
## - For more accuracy, run bootnet(..., caseMin = 0.672, caseMax = 1)
##
## Accuracy can also be increased by increasing both 'nBoots' and 'caseN'.
ClaAAEDataA <- dplyr::select(FinalData, 10:32)
colnames(ClaAAEDataA) <- c("9","1","2","3","4","5","6","7","8","10","11","12","13","14","15","16","17","18","19","20","21","22","23")
order2 <- c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","23")
ClaAAEData <- ClaAAEDataA[, order2]
EGAClaAAE<-EGA(ClaAAEData, plot.EGA = TRUE, corr = c("spearman"), model = c("glasso"), algorithm = c("walktrap"), plot.type = c("qgraph"))
EGAClaAAE$dim.variables
## items dimension
## 7 7 1
## 8 8 1
## 9 9 1
## 11 11 1
## 13 13 1
## 19 19 1
## 20 20 1
## 22 22 1
## 1 1 2
## 2 2 2
## 3 3 2
## 4 4 2
## 5 5 2
## 10 10 2
## 12 12 2
## 21 21 2
## 23 23 2
## 14 14 3
## 16 16 3
## 17 17 3
## 18 18 3
## 6 6 NA
## 15 15 NA
AAEDataC <- dplyr::select(FinalData,10:39)
colnames(AAEDataC) <- c("9","1","2","3","4","5","6","7","8","10","11","12","13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30")
#order according to clusters
order3 <- c("1","2","3","4","5","10","12","21","23","7","8","9","11","13","19","20","22","14","16","17","18","6","15","24","25","26","27","28","29","30")
AAEDataMH <- AAEDataC[, order3]
nodeNamesAAEs <- c("Physical Abuse inside of Family","Physical Abuse outside of Family","Mental Abuse inside of Family","Mental Abuse outside of Family","Sexual Abuse","Conflicts with Parents","Trouble with Police","Abuse by Romantic Partner","Occupational Issues", "Abuse of Mother","Criminality of Carer","Carer Psychopathology","Divorce of Parents","Health Issues","Experience of Life/Death Situation","Housing Issues","Death to Close Contact","Educational Issues","Being Bullied","Feeling Lonely","Number of Friends", "Carer Substance Abuse","Teenager became Parent","Bipolar Disorder", "Depression", "Chronic Fatigue Syndrom","Alcohol Use Disorder","Drug Use", "Happiness","Life Satisfaction")
Groups2 <- c(rep("1. Cluster - Direkt Abuse",9),rep("2. Cluster - Family Factors",8),rep("3. Cluster - Educational and Social Factors",4),rep("4. No cluster",1), rep("5. No cluster",1), rep("Mental Health",5), rep("Wellbeing Factors",2))
NetworkAAEMH<- estimateNetwork(AAEDataMH, default = "ggmModSelect", corMethod = "spearman")
plot(NetworkAAEMH, groups = Groups2, nodeNames = nodeNamesAAEs, legend.cex = 0.25, posCol = "gray25", negCol = "grey", negDashed = TRUE, fade = FALSE, esize = 8, color = c("#440154FF","#31688EFF", "#26828EFF", "snow3", "snow3", "#6DCD59FF", "#B4DE2CFF"), border.color = "white",border.width = 2, label.color = "white")
centralityPlot(NetworkAAEMH, include = c("Strength", "Closeness","Betweenness","ExpectedInfluence"))
AAEedges <- getWmat(NetworkAAEMH)
write.csv(AAEedges, "AAEedges.csv")
AAEedgeplot <- plot(NetworkAAEMH,layoutScale = c(0.9,0.9), posCol = "gray25", negCol = "grey", negDashed = TRUE, fade = FALSE, esize = 10, color = c("#440154FF", "#3B528BFF", "#21908CFF", "#5DC863FF", "#9FDA3AFF"),vsize = 7)
AAEbridge <- bridge(AAEedgeplot, 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','2','2','2','2','2','2','2'), useCommunities = "all", directed = NULL, nodes = NULL)
pdf("bridgecentralityAAE.pdf", width = 5)
plot(AAEbridge, include = "Bridge Strength")
dev.off()
Groups8 <- c(rep("2. AAEs",2),
rep("1. Bridge Nodes AAEs Cluster:
3: Mental Abuse inside of Family,
12: Trouble with Police,
14: Educational Issues,
21: Abuse by Romantic Partner", 1),
rep("2. AAEs",3),
rep("1. Bridge Nodes AAEs Cluster:
3: Mental Abuse inside of Family,
12: Trouble with Police,
14: Educational Issues,
21: Abuse by Romantic Partner",2), #8
rep("2. AAEs",9),
rep("1. Bridge Nodes AAEs Cluster:
3: Mental Abuse inside of Family,
12: Trouble with Police,
14: Educational Issues,
21: Abuse by Romantic Partner",1), #18
rep("2. AAEs",5), #23
rep("4. Mental Health & Wellbeing",1),
rep("3. Bridge Nodes Mental Health Cluster:
25: Depression,
27: Alcohol Use Disorder,
28: Drug Use,
30: Life Satisfaction",1),
rep("4. Mental Health & Wellbeing",1),
rep("3. Bridge Nodes Mental Health Cluster:
25: Depression,
27: Alcohol Use Disorder,
28: Drug Use,
30: Life Satisfaction",2),
rep("4. Mental Health & Wellbeing",1),
rep("3. Bridge Nodes Mental Health Cluster:
25: Depression,
27: Alcohol Use Disorder,
28: Drug Use,
30: Life Satisfaction",1))
plot(NetworkAAEMH, legend.cex = 0.3, groups = Groups8, posCol = "gray25", negCol = "grey", negDashed = TRUE, fade = FALSE, esize = 10, color = c("#440154FF","snow4","#21908CFF","snow3"),border.color = "white",border.width = 2, label.color = "white", vsize = 7,curve = 0.1, curveAll = T)
#Remove "#" to run analysis
#BNetworkAAEs<- bootnet(NetworkAAEMH, default = "ggmModSelect", nBoots = 2000, nCores = 8)
#save(BNetworkAAEs , file = "BNetworkAAEs.Rdata")
load("BNetworkAAEs.Rdata")
plot(BNetworkAAEs, order = "sample", labels = FALSE)
#B2NetworkAAE<- bootnet(NetworkAAEMH, default = "ggmModSelect", type = "case", nBoots = 2000, nCores = 8, statistics = c("edge", "strength", "betweenness", "closeness","expectedInfluence"))
#save(B2NetworkAAE , file = "B2NetworkAAE.Rdata")
load("B2NetworkAAE.Rdata")
plot(B2NetworkAAE, statistics = c("strength", "closeness","betweenness", "expectedInfluence"))
corStabNetworkAAE<- corStability(B2NetworkAAE, cor = 0.7)
## === Correlation Stability Analysis ===
##
## Sampling levels tested:
## nPerson Drop% n
## 1 3911 75.0 24
## 2 5127 67.2 69
## 3 6344 59.4 97
## 4 7561 51.7 156
## 5 8777 43.9 209
## 6 9994 36.1 234
## 7 11211 28.3 261
## 8 12427 20.6 280
## 9 13644 12.8 343
## 10 14861 5.0 327
##
## Maximum drop proportions to retain correlation of 0.7 in at least 95% of the samples:
##
## betweenness: 0.361
## - For more accuracy, run bootnet(..., caseMin = 0.283, caseMax = 0.439)
##
## closeness: 0.439
## - For more accuracy, run bootnet(..., caseMin = 0.361, caseMax = 0.517)
##
## edge: 0.672
## - For more accuracy, run bootnet(..., caseMin = 0.594, caseMax = 0.75)
##
## expectedInfluence: 0.75 (CS-coefficient is highest level tested)
## - For more accuracy, run bootnet(..., caseMin = 0.672, caseMax = 1)
##
## strength: 0.594
## - For more accuracy, run bootnet(..., caseMin = 0.517, caseMax = 0.672)
##
## Accuracy can also be increased by increasing both 'nBoots' and 'caseN'.
Network Comparison Test (NCT) is a permutation test. The NCT estimates the network structure and calculates a metric that functions as the observed test statistic. Group membership is multiple times rearranged via permutation, followed by a recalculation of the network structure and test statistic. This is resulting in a reference distribution. The NCT then compares the first observed test statistic with this reference distribution, indicating whether the observed test statistic is significantly different. For network comparisons, please refer to Borkulo et al. (2017).
#Adapt ACEs and AAEs to the same number of nodes and names
ACEData2 <- data.frame(ACEData,SUBACE=rowMeans(ACEData[6:7])) #this adds x in front of the column names
ACEData3 <- ACEData2[, -c(6:7)]
orderACEData <- c("X1","X2","X3","X4","X5","SUBACE","X8","X9","X10","X11","X12","X13","X14","X15","X16","X17")
ACEData4 <- ACEData3[, orderACEData]
colnames(ACEData4) <- c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16")
nodeNamesCom <- c("Physical Abuse inside of Family","Physical Abuse outside of Family","Mental Abuse inside of Family", "Mental Abuse outside of Family", "Sexual Abuse", "Carer Substance Abuse", "Abuse of Mother", "Criminality of Carer", "Psychopathology of Carer", "Bipolar Disorder", "Depression", "Chronic Fatigue Syndrom","Alcohol Use Disorder","Drug Use", "Happiness","Life Satisfaction")
Groups4 <- c(rep("ACEs",9),rep("Mental Health & Wellbeing",7))
BNetworkACEModSelSmall<- estimateNetwork(ACEData4, default = "ggmModSelect", corMethod = "spearman")
plot(BNetworkACEModSelSmall, groups = Groups4, layout = "circle", edge.labels = FALSE, edge.label.cex = 0.5, edge.label.color = "black", GLratio = 2.5, nodeNames = nodeNamesCom, legend.cex = 0.4, layoutOffset = -0.05,vsize = 4.5, label.cex = 2, curve = 0.3, curveAll = T, edge.label.position = 0.5, posCol = "gray25", negCol = "grey", negDashed = TRUE, fade = FALSE, esize = 8, color = c("#5DC863FF", "#9FDA3AFF"), border.color = "white",border.width = 2, label.color = "white")
centralityPlot(BNetworkACEModSelSmall, include = c("Strength", "Closeness","Betweenness","ExpectedInfluence"))
#Adapt AAEs
AAEDataSmallB <- dplyr::select(FinalData, 10:18,33:39)
colnames(AAEDataSmallB) <- c("9","1","2","3","4","5","6","7","8","10","11","12","13","14","15","16")
order2 <- c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16")
AAEDataSmall <- AAEDataSmallB[, order2]
Groups3 <- c(rep("AAEs",9),rep("Mental Health & Wellbeing",7))
BNetworkAAEModSelSmall<- estimateNetwork(AAEDataSmall, default = "ggmModSelect", corMethod = "spearman")
plot(BNetworkAAEModSelSmall, groups = Groups3, layout = "circle",edge.labels = FALSE, edge.label.cex = 0.5, edge.label.color = "black", GLratio = 2.5, nodeNames = nodeNamesCom, legend.cex = 0.4, layoutOffset = -0.05,vsize = 4.5, label.cex = 2, curve = 0.3, curveAll = T, edge.label.position = 0.5, posCol = "gray25", negCol = "grey", negDashed = TRUE, fade = FALSE, esize = 8, color = c("#440154FF", "#21908CFF"), border.color = "white",border.width = 2, label.color = "white")
centralityPlot(BNetworkAAEModSelSmall, include = c("Strength", "Closeness","Betweenness","ExpectedInfluence"))
#Please remove "#" to run analysis
#NetworkCom <- NCT(BNetworkACEModSelSmall, BNetworkAAEModSelSmall, it=2000, binary.data=FALSE, weighted=TRUE, test.edges=FALSE, edges="ALL",paired=TRUE)
#save(NetworkCom, file = "NetworkCom.Rdata")
load("NetworkCom.Rdata")
summary(NetworkCom)
##
## NETWORK INVARIANCE TEST
## Test statistic M: 0.0919
## p-value 0.248
##
## GLOBAL STRENGTH INVARIANCE TEST
## Global strength per group: 3.53 3.78
## Test statistic S: 0.243
## p-value 0.142
##
## EDGE INVARIANCE TEST
## Edges tested:
## Test statistic E:
## p-value
##
## CENTRALITY INVARIANCE TEST
## Nodes tested:
## Centralities tested:
## Test statistic C:
## p-value
#Network structure invariance test
plot(NetworkCom, what="network")
#Global strength invariance test
plot(NetworkCom, what="strength")
#Edge invariance test
#plot(NetworkCom, what="edge")
We computed path diagrams in childhood and adolescence based on the shortest pathway analysis. For more information, please refer to Fritz et al. (2019) and Isvoranu et al. (2020).
ACEData5 <- dplyr::select(ACEData4,1:9,14)
ACEData5B <- dplyr::select(ACEData4,1:9,11)
AAEDataSmallC <- dplyr::select(AAEDataSmall,1:9,11)
AAEDataSmallCB <- dplyr::select(AAEDataSmall,1:9,14)
ACECorMat <- cor(ACEData5, use = "pairwise.complete.obs")
AAEcorMat <- cor(AAEDataSmallC, use = "pairwise.complete.obs")
ACECorMat2 <- cor(ACEData5B, use = "pairwise.complete.obs")
AAEcorMat2 <- cor(AAEDataSmallCB, use = "pairwise.complete.obs")
ACEGraph <- qgraph(ACECorMat, graph = "glasso", layout = "spring", tuning = 0.25, sampleSize = nrow(ACEData4), minimum = 0,cut = 0.15, maximum = 1, details = TRUE, esize = 20, posCol = "gray25", negCol = "grey", fade = FALSE, esize = 8, color = c("#5DC863FF"), border.color = "white",border.width = 2, label.color = "white")
ACEGraph2<- qgraph(ACECorMat2, graph = "glasso", layout = "spring", tuning = 0.25, sampleSize = nrow(ACEData4), minimum = 0,cut = 0.15, maximum = 1, details = TRUE, esize = 20, posCol = "gray25", negCol = "grey", fade = FALSE, esize = 8, color = c("#5DC863FF"), border.color = "white",border.width = 2, label.color = "white")
AAEGraph <- qgraph(AAEcorMat, graph = "glasso", layout = "spring", tuning = 0.25, sampleSize = nrow(AAEDataSmall), minimum = 0, cut = 0.15, maximum = 1, details = TRUE, esize = 20,posCol = "gray25", negCol = "grey", fade = FALSE, esize = 8, color = c("#21908CFF"), border.color = "white",border.width = 2, label.color = "white")
AAEGraph2 <- qgraph(AAEcorMat2, graph = "glasso", layout = "spring", tuning = 0.25, sampleSize = nrow(AAEDataSmall), minimum = 0, cut = 0.15, maximum = 1, details = TRUE, esize = 20,posCol = "gray25", negCol = "grey", fade = FALSE, esize = 8, color = c("#21908CFF"), border.color = "white",border.width = 2, label.color = "white")
pathways(ACEGraph, from = c(1,2,3,4,5,6,7,8,9), to = c(10))
pathways(AAEGraph, from = c(1,2,3,4,5,6,7,8,9), to = c(10))
flow(ACEGraph, from = c(10), horizontal = TRUE, equalize = TRUE, minCurve = 1, maxCurve = 4, unfadeFirst = FALSE, fade =TRUE, sampleSize=nrow(ACEData5))
flow(ACEGraph2, from = c(10), horizontal = TRUE, equalize = TRUE, minCurve = 1, maxCurve = 4, unfadeFirst = FALSE, fade =TRUE, sampleSize=nrow(ACEData5))
flow(AAEGraph, from = c(10), horizontal = TRUE, equalize = TRUE, minCurve = 1, maxCurve = 4, unfadeFirst = FALSE, fade =TRUE, sampleSize=nrow(AAEDataSmallC))
flow(AAEGraph2, from = c(10), horizontal = TRUE, equalize = TRUE, minCurve = 1, maxCurve = 4, unfadeFirst = FALSE, fade =TRUE, sampleSize=nrow(AAEDataSmallC))
If you have any questions, please don’t hesitate to get in touch.