# Create new temporary DBQ datasets for data manipulation, replacing the previous DBQ1temp: DBQ1temp <- DEQ1[,c(1, 52:87, 40, 88, 107)] # 52 = too high a gear; 40 = slippery road, 88 = drugs, 107 = swerve to avoid accident DBQ2temp <- DEQ2[,c(1, 140:177, 184)] # 140 = too high a gear; 176 = slippery road, 177 = drugs, 184 = swerve to avoid accident DBQ3temp <- DEQ3[,c(1, 181:218, 140)]# 181 = too high a gear; 217 = slippery road, 218 = drugs, 140 = swerve to avoid accident DBQ4temp <- DEQ4[,c(1, 181:218, 140)] # Include variables of DBQ27, an item related to mobile phone use (19), items related to alcohol (22) and drugs (38) and # item related to near misses (39). In other words, drop the following: DBQ_Drop <- c(4,9,11,17,24,29,30,34) DBQ1temp <- DBQ1temp[,-DBQ_Drop] DBQ2temp <- DBQ2temp[,-DBQ_Drop] DBQ3temp <- DBQ3temp[,-DBQ_Drop] DBQ4temp <- DBQ4temp[,-DBQ_Drop] nomissings <- which(!is.na(rowSums( cbind(DBQ1temp, DBQ2temp, DBQ3temp, DBQ4temp)))) # indexes without any missing values DBQ1temp <- DBQ1temp[nomissings,] DBQ2temp <- DBQ2temp[nomissings,] DBQ3temp <- DBQ3temp[nomissings,] DBQ4temp <- DBQ4temp[nomissings,] DBQ1temp$time <- 1 DBQ2temp$time <- 2 DBQ3temp$time <- 3 DBQ4temp$time <- 4 # Move variable "time" into the beginning of the datasets DBQ1temp <- DBQ1temp[,c(1,33,2:32)] DBQ2temp <- DBQ2temp[,c(1,33,2:32)] DBQ3temp <- DBQ3temp[,c(1,33,2:32)] DBQ4temp <- DBQ4temp[,c(1,33,2:32)] # Which of the variables are errors? errorindex <- c(3, 5, 6, 8, 9, 10, 12, 14, 16, 19, 22, 23, 25, 26, 27, 28, 30, 31, 33) errors <- names(DBQ1temp)[errorindex] violations <- names(DBQ1temp)[!(names(DBQ1temp) %in% c("ID","time",errors))] seq(1:33)[!seq(1:33)%in%errorindex] # In fact, let's name the variables differently (V = violation, E = error) # Too tedious to do by hand, so let's write a function: name <- function(inputNames) { outputNames <- character() NetworkGroups <- character() j <- k <- 1 for (i in 1:length(inputNames)) { #for (j in 1:length(errors)) { if (inputNames[i] %in% errors) { outputNames[i] <- paste("e", j , sep= "") NetworkGroups[i] <- "Error" j <- j + 1 } else if (inputNames[i] %in% violations) { outputNames[i] <- paste("v", k , sep= "") NetworkGroups[i] <- "Violation" k <- k + 1 } } outputNames[1] <- "ID" outputNames[2] <- "time" return(list(names = outputNames, groups = NetworkGroups)) } namesNGroups <- name(names(DBQ1temp)) namesNGroups$groups <- namesNGroups$groups[3:33] namesNGroups$names <- namesNGroups$names[3:33] # Rename the variables in the dataframes to enable using rbind to combine the datasets: names(DBQ1temp) <- c("ID", "time", paste("n",seq(1:31),sep="")) names(DBQ2temp) <- names(DBQ3temp) <- names(DBQ4temp) <- names(DBQ1temp) # Combine the datasets: DBQtemp <- rbind(DBQ1temp,DBQ2temp,DBQ3temp,DBQ4temp) table(DBQtemp$time) # 1173 observations # Create between-subject variables (means across the four data points): DBQBetween <- group_by(DBQtemp, ID) %>% summarize( n1.m = mean(n1), n2.m = mean(n2), n3.m = mean(n3), n4.m = mean(n4), n5.m = mean(n5), n6.m = mean(n6), n7.m = mean(n7), n8.m = mean(n8), n9.m = mean(n9), n10.m = mean(n10), n11.m = mean(n11), n12.m = mean(n12), n13.m = mean(n13), n14.m = mean(n14), n15.m = mean(n15), n16.m = mean(n16), n17.m = mean(n17), n18.m = mean(n18), n19.m = mean(n19), n20.m = mean(n20), n21.m = mean(n21), n22.m = mean(n22), n23.m = mean(n23), n24.m = mean(n24), n25.m = mean(n25), n26.m = mean(n26), n27.m = mean(n27), n28.m = mean(n28), n29.m = mean(n29), n30.m = mean(n30), n31.m = mean(n31)) # Calculate polychoric correlations for the between-subjects variables: corBetween <- select(DBQBetween, n1.m:n31.m) %>% cor_auto() # Create between-subjects network model: networkBetween <- EBICglasso(S = corBetween, n = nrow(DBQBetween)) # Create a vector for the DBQ variable names: dbqnames <- c( "drive away from traffic lights at too high a gear", "overtake a slow driver on inside", "attempt to overtake and hadn't noticed signalling right", "forget where left car in carpark", "sound horn to indicate annoyance", "switch on one thing when meant to switch on other", "pull out of junction so far that driver has to let you out", "realised have no recollection of road been travelling", "crossed junction knowing lights have turned against you", "failed to notice people crossing when turned into sidestreet", "become angered by driver and given chase", "misread signs and taken wrong turning off roundabout", "disregarded speed limit on residential road", "turning left, nearly hit cyclist on inside", "used mobile phone without hands free kit", "stay in motorway lane you know will be closed", "when queuing to turn left nearly hit car in front", "drive when suspect over legal alcohol limit", "become angered by driver and indicated hostility", "misjudged speed of oncoming vehicle when overtaking", "hit something when reversing that hadn't seen", "raced away from traffic lights to beat other driver", "noticed ending up on a different road than intended", "get into wrong lane when approaching roundabout/junction", "drive so close to car that would not be able to stop", "missed giveway signs and avoided colliding with traffic", "disregarded speed limit on motorway", "failed to check rear-view mirror before manoeuvring", "brake too quickly on slippery road or steer wrong in skid", "drove after taking drugs which affected you", "had to brake or swerve to avoid accident" ) # Check that things match: dim(networkBetween) length(dbqnames) length(namesNGroups$groups) # Have a quick look at the resulting network: DBQBetweenNet <- qgraph(networkBetween, layout = "spring", nodeNames = dbqnames, groups = namesNGroups$groups, labels = namesNGroups$names, vsize = 3, filename = "DBQBetween", filetype = "pdf", legend.cex = .23, legend = TRUE) # Have a look at the between-subjects centrality plot: pdf("centrality_Between.pdf") centralityPlot(DBQBetweenNet) dev.off() # Figure 3 is created after estimating the cross-sectional network, since the two networks are scaled to # the same maximum value ################################ Cross-sectional network model, Fig 5 ################################ NetworkDataDBQ_BGVars <- cbind(DBQ1temp2,LTDQAttInclude,SIsInclude,ImprInclude) NetworkDataDBQ_BGVars <- na.omit(NetworkDataDBQ_BGVars) dim(NetworkDataDBQ_BGVars) #8865 40 NetworkGroupsDBQ_BGVars <- c(namesNGroups$groups, rep("Attitude",2), rep("Self-image", 2), rep("Improvement needs",5)) NetworkDataNamesAll <- c(dbqnames,includesAttNames,includesSINames, includesImprNames) # It is perhaps more intuitive to interpret the irritable / placid self image variable coded the other way around: NetworkDataDBQ_BGVars$V00174 <- 7 - NetworkDataDBQ_BGVars$V00174 # And these need to be recoded: NetworkDataDBQ_BGVars$V00095 <- 6 - NetworkDataDBQ_BGVars$V00095 NetworkDataDBQ_BGVars$V00099 <- 6 - NetworkDataDBQ_BGVars$V00099 labelsDBQBG <- c(namesNGroups$names, paste("a", seq(1:length(includesAtt)), sep = ""), paste("si", seq(1:length(includesSI)), sep = ""), paste("in", seq(1:length(includesImpr)), sep = "")) DBQBGCor <- cor_auto(NetworkDataDBQ_BGVars, ordinalLevelMax = 7) DBQBGNet <- estimateNetwork(data = NetworkDataDBQ_BGVars, prepFun = "cor_auto", prepArgs = list(ordinalLevelMax = 7), estFun = "EBICglasso", estArgs = list(n = nrow(NetworkDataDBQ_BGVars))) # Find the maximum value of the two networks combined maxTwoNetworks <- max(abs(c(DBQBGNet$graph,networkBetween))) maxTwoNetworks # [1] 0.4771615 Fig5 <- qgraph(DBQBGCor, graph = "glasso", sampleSize = nrow(NetworkDataDBQ_BGVars), layout = "spring", vsize = 3, nodeNames = NetworkDataNamesAll, groups = NetworkGroupsDBQ_BGVars, legend.cex = .3, labels = labelsDBQBG, maximum = maxTwoNetworks, color = c("#ffcc00", "#00B2FFFF", "#00FF66FF", "#CC00FFFF", "#FF0000FF")) pdf("Fig5.pdf", width = 10) makeBW(Fig5) dev.off() ########################################################################################### ##################################### ##################################### ##################################### Create Figure 3 ##################################### ##################################### ##################################### ########################################################################################### Fig3 <- qgraph(networkBetween, layout = "spring", nodeNames = dbqnames, groups = namesNGroups$groups, labels = namesNGroups$names, legend.cex = .3, legend = TRUE, maximum = maxTwoNetworks, vsize = 3, color = c("#00B2FFFF", "#FF0000FF")) pdf("Fig3.pdf", width = 9.5) makeBW(Fig3) dev.off() ########################################################################################### ########################################################################################### ##################################### ##################################### ##################################### Create Figure 4 ##################################### ##################################### ##################################### ########################################################################################### # Compute centralities and the Zhang signed clustering coefficient centrBetween <- centralityTable(Fig3) #strength closeness centrBetween <- centrBetween[c(centrBetween$measure == "Strength" | centrBetween$measure == "Closeness"),] centrBetween$type <- "Between-person" clusterBetween <- clusteringTable(Fig3, signed = TRUE) %>% filter(measure == "Zhang") clusterBetween$type <- "Between-person" clusterBetween$measure <- "Clust. coef." # Arrange the nodes correctly centrClustBetween <- rbind(centrBetween, clusterBetween) %>% arrange(node, measure) centrClustBetweenList <- list( "Closeness" = centrClustBetween %>% filter(measure == "Closeness") %>% arrange(desc(value)), "Strength" = centrClustBetween %>% filter(measure == "Strength") %>% arrange(desc(value)), "Clust" = centrClustBetween %>% filter(measure == "Clust. coef.") %>% arrange(desc(value))) centrClustBetweenList$plotClose <- ggplot(centrClustBetweenList$Closeness, aes(x = value, y = reorder(node,value), group = 1)) + geom_path(alpha = 1, size = .5) + geom_point(size = 1.5) + xlab("") + ylab("") + ggtitle("Closeness") + theme_bw() + theme(plot.title = element_text(hjust = .5 )) centrClustBetweenList$plotStrength <- ggplot(centrClustBetweenList$Strength, aes(x = value, y = reorder(node,value), group = 1)) + geom_path(alpha = 1, size = .5) + geom_point(size = 1.5) + xlab("") + ylab("") + ggtitle("Strength") + theme_bw() + theme(plot.title = element_text(hjust = .5 )) centrClustBetweenList$plotClust <- ggplot(centrClustBetweenList$Clust, aes(x = value, y = reorder(node,value), group = 1)) + geom_path(alpha = 1, size = .5) + geom_point(size = 1.5) + xlab("") + ylab("") + ggtitle("Clustering coefficient") + theme_bw() + theme(plot.title = element_text(hjust = .5 )) grid.arrange(centrClustBetweenList$plotClose,centrClustBetweenList$plotStrength,centrClustBetweenList$plotClust, ncol=3) # Save the figure as a pdf pdf("Fig4.pdf", width = 8, height = 8) grid.arrange(centrClustBetweenList$plotClose,centrClustBetweenList$plotStrength,centrClustBetweenList$plotClust, ncol=3) dev.off() ########################################################################################### ##################################### ##################################### ##################################### Create Figure 6 ##################################### ##################################### ##################################### ########################################################################################### # Compute centralities and the Zhang signed clustering coefficient centrCrossSect <- centralityTable(Fig5) centrCrossSect$type <- "Cross-sectional" clusterCrossSect <- clusteringTable(Fig5, signed = TRUE) %>% filter(measure == "Zhang") clusterCrossSect$type <- "Cross-sectional" clusterCrossSect$measure <- "Clust. coef." # Arrange the nodes correctly centrClustCrossSect <- rbind(centrCrossSect, clusterCrossSect) %>% arrange(node, measure) centrClustCrossSectList <- list( "Betweenness" = centrClustCrossSect %>% filter(measure == "Betweenness") %>% arrange(desc(value)), "Closeness" = centrClustCrossSect %>% filter(measure == "Closeness") %>% arrange(desc(value)), "Strength" = centrClustCrossSect %>% filter(measure == "Strength") %>% arrange(desc(value)), "Clust" = centrClustCrossSect %>% filter(measure == "Clust. coef.") %>% arrange(desc(value))) centrClustCrossSectList$plotBetween <- ggplot(centrClustCrossSectList$Betweenness, aes(x = value, y = reorder(node,value), group = 1)) + geom_line(alpha = 1, size = .5) + geom_point(size = 1.5) + xlab("") + ylab("") + ggtitle("Betweenness") + theme_bw() + theme(plot.title = element_text(hjust = .5 )) centrClustCrossSectList$plotClose <- ggplot(centrClustCrossSectList$Closeness, aes(x = value, y = reorder(node,value), group = 1)) + geom_path(alpha = 1, size = .5) + geom_point(size = 1.5) + xlab("") + ylab("") + ggtitle("Closeness") + theme_bw() + theme(plot.title = element_text(hjust = .5 )) centrClustCrossSectList$plotStrength <- ggplot(centrClustCrossSectList$Strength, aes(x = value, y = reorder(node,value), group = 1)) + geom_path(alpha = 1, size = .5) + geom_point(size = 1.5) + xlab("") + ylab("") + ggtitle("Strength") + theme_bw() + theme(plot.title = element_text(hjust = .5 )) centrClustCrossSectList$plotClust <- ggplot(centrClustCrossSectList$Clust, aes(x = value, y = reorder(node,value), group = 1)) + geom_path(alpha = 1, size = .5) + geom_point(size = 1.5) + xlab("") + ylab("") + ggtitle("Clustering coefficient") + theme_bw() + theme(plot.title = element_text(hjust = .5 )) # Save the figure as a pdf pdf("Fig6.pdf", width = 12, height = 8) grid.arrange(centrClustCrossSectList$plotBetween, centrClustCrossSectList$plotClose,centrClustCrossSectList$plotStrength, centrClustCrossSectList$plotClust, ncol=4) dev.off()