############################################################################################################################
################################# Correlation and adjacency matrices #######################################################
#################################             Tables S1-S4           #######################################################
############################################################################################################################

# Correlation matrix for the between-person network:
write(print(xtable(corBetween), type = "html"), "corBetween.html")

# Adjacency matrix for the between-person network:
write(print(xtable(networkBetween), type = "html"), "adjBetween.html")

# Correlation matrix for the cross-sectional network:
write(print(xtable(DBQBGCor), type = "html"), "corDBQ_BG_CrossSect.html")

# Adjacency matrix for the cross-sectional network:
write(print(xtable(DBQBGNet$graph), type = "html"), "adjDBQ_BG_CrossSect.html")


############################################################################################
############################                      ##########################################
############################                      ##########################################
############################ Between-person model ##########################################
############################                      ##########################################
############################                      ##########################################
############################################################################################
############################  Bootstrap analyses  ##########################################
############################################################################################

# Prepare bootnet analysis for the between-person network:
BetweenNetForBootnet <- estimateNetwork(data = as.data.frame(select(DBQBetween, n1.m:n31.m)), prepFun = "cor_auto", 
                                        prepArgs = list(ordinalLevelMax = 7), estFun = "EBICglasso",
                                        estArgs = list(n = nrow(DBQLongit.m)))

# Edge weight accuracy, between-person network:
set.seed(2018)
BootnetBetween$BetweenEdgeWeight <- bootnet(BetweenNetForBootnet, nBoots = 2500)

# Centrality stability, between-person network:
set.seed(2018)
BootnetBetween$BetweenStability <- bootnet(BetweenNetForBootnet, nBoots = 2500, type = "case")
set.seed(2018)
corStability(BootnetBetween$BetweenStability)

################### Accuracy & stability plot (Supplementary Figure S1): #########################################################
pdf("Between_Stability.pdf", width = 12)
grid.arrange(
  plot(BootnetBetween$BetweenEdgeWeight, order = "sample", labels = FALSE, samplelwd = .5,panels = FALSE),
  plot(BootnetBetween$BetweenStability),
  ncol = 2)
dev.off()
###################################################################################################################################

# Accuracy plot in a longer format so it is easier to examine individual edge weight estimates:
pdf("Between_long.pdf", height = 45)
plot(BootnetBetween$BetweenEdgeWeight, order = "sample", labels = TRUE, samplelwd = .5,panels = TRUE)
dev.off()


###### Between-subject network in summary: Closeness and strength indices can be interpreted
############################################################################################


############################################################################################
############################                      ##########################################
###########################                        #########################################
########################### Cross-sectional model  #########################################
###########################                        #########################################
############################                      ##########################################
############################################################################################
###################  Choosing the background variables  ####################################
############################################################################################
################################ DBQ and Attitudes #########################################
############################################################################################

# Select the attitude variables:
ltdqAttVars <- as.character(expression(V00095,V00096,V00097,V00098,V00099,V00100,V00101,V00102))
LTDQAttAll <- LTDQ[,ltdqAttVars]

# The variables are counterintuitively coded, let's change:
LTDQAttAll <- as.data.frame(lapply(LTDQAttAll, function(x) 6-x))

DBQ1temp2 <- DEQ1[,c(1,52:87, 40, 88, 107)] # 52 = too high a gear; 40 = slippery road, 88 = drugs, 107 = swerve to avoid accident
DBQ1temp2 <- DBQ1temp2[,-DBQ_Drop]
DBQ1temp2$ID <- NULL
NetworkDataDBQAtt <- cbind(DBQ1temp2,LTDQAttAll)
NetworkDataDBQAtt <- na.omit(NetworkDataDBQAtt)
dim(NetworkDataDBQAtt) # 8966 39

# Choosing the cases this way means that the numbers of cases will slightly differ 
# for the different background variable analyses. On the other hand, carrying out listwise deletion
# for all background variables combined would be quite wasteful.

varnamesAttAll <- c(
  "Decreasing motorway speed limit is a good idea",
  "Even at night-time on quiet roads important to keep within speed limit",
  "Drivers who cause accidents by reckless driving should be banned for life",
  "People should drive slower than speed limits when it is raining",
  "Cars should never overtake on inside lane",
  "In towns with lot of pedestrians the speed limit should be 20mph",
  "Penalties for speeding should be more severe",
  "Increasing motorway speed limit is a good idea")

NetworkDataNamesDBQAtt <- c(dbqnames,varnamesAttAll)
labelsDBQAtt <- c(namesNGroups$names, paste("a", seq(1:length(varnamesAttAll)), sep = ""))
NetworkGroupsDBQAtt <- c(namesNGroups$groups,rep("Attitudes", 8))
                  
NetworkCorDBQAtt <- cor_auto(NetworkDataDBQAtt, ordinalLevelMax = 7)

CrossSectNetworkDBQAtt <- qgraph(NetworkCorDBQAtt, graph = "glasso", sampleSize = nrow(NetworkDataDBQAtt), layout = "spring",
                                 vsize = 3, nodeNames = NetworkDataNamesDBQAtt, groups = NetworkGroupsDBQAtt, legend.cex = .3,
                                 labels = labelsDBQAtt, color = c("#00ff00", "#00B2FFFF", "#FF0000FF"))

############### Create Supplementary Figure S3: ###############

pdf("CrossSectNetworkDBQAtt.pdf", width = 11.5)
makeBW(CrossSectNetworkDBQAtt)
dev.off()
###############################################################
# The attention variables to include:
includesAtt <- c("V00095", "V00099")
includesAttNames <- varnamesAttAll[which(names(LTDQAttAll) %in% includesAtt)]
LTDQAttInclude <- LTDQ[,includesAtt]

############################################################################################
############################ DBQ and Improvement needs######################################
############################################################################################

# Select the improvement needs variables:

ImpVars <- as.character(expression(V00245,V00246,V00247,V00248,V00249,V00250,V00251,V00252,V00253,V00254,V00255,V00256,V00257,V00258,V00259,V00260,
                                   V00261,V00262,V00263,V00264))

IMPS <- DEQ1[,ImpVars]

NetworkDataDBQImpr <- cbind(DBQ1temp2,IMPS)
NetworkDataDBQImpr <- na.omit(NetworkDataDBQImpr)
dim(NetworkDataDBQImpr) #8605 51

varnamesImprovement <- c(
  "improve car controls",
  "improve on pulling out of junctions",
  "improve on reversing",
  "improve parking",
  "improve judging speed of others",
  "improve judging what others going to do",
  "improve on spotting hazards",
  "improve on driving in heavy traffic",
  "improve on driving in the dark",
  "improve on overtaking",
  "improve on using roundabouts",
  "improve on joining with moving traffic",
  "improve on changing lanes",
  "improve on driving on high speed roads",
  "improve on driving on country roads",
  "improve on driving in heavy rain",
  "improve on driving in thick fog",
  "improve on driving on snow or ice",
  "improve on turning right",
  "improve on knowing what speed is safe")

NetworkDataNamesDBQImpr <- c(dbqnames,varnamesImprovement)
labelsDBQImpr <- c(namesNGroups$names, paste("i", seq(1:length(varnamesImprovement)), sep = ""))
NetworkGroupsDBQImpr <- c(namesNGroups$groups,rep("Improvement needs", length(varnamesImprovement)))

NetworkCorDBQImpr <- cor_auto(NetworkDataDBQImpr, ordinalLevelMax = 7)

CrossSectNetworkDBQImpr <- qgraph(NetworkCorDBQImpr, graph = "glasso", sampleSize = nrow(NetworkDataDBQImpr), 
                                  layout = "spring",vsize = 3, nodeNames = NetworkDataNamesDBQImpr, 
                                  groups = NetworkGroupsDBQImpr, labels = labelsDBQImpr, legend.cex = .3,
                                  color = c("#00B2FFFF", "#00ff00", "#FF0000FF"))

############### Create Supplementary Figure S4: ###############                                  
pdf("CrossSectNetworkDBQImpr.pdf", width = 10)
makeBW(CrossSectNetworkDBQImpr)
dev.off()
###############################################################

includesImpr <- c("V00245", "V00248", "V00251", "V00257", "V00264")
includesImprNames <- varnamesImprovement[which(names(IMPS) %in% includesImpr)]
ImprInclude <- IMPS[includesImpr]

############################################################################################
############################### DBQ and self-image ########################################
############################################################################################

# Select the self-image variables:

SelfImageVars <- as.character(expression(V00316,V00177,V00178,V00179,V00180,V00185,V00150,V00154,V00155,V00161,V00174,V00186,V00288,V00298,V00302,V00303,V00304,
                                         V00309))
SIs <- DEQ1[,SelfImageVars]

NetworkDataDBQSI <- cbind(DBQ1temp2,SIs)
NetworkDataDBQSI <- na.omit(NetworkDataDBQSI)
dim(NetworkDataDBQSI) # 8872 49

varnamesSelfImage <- c(
  "generally, confident are you in your driving ability",
  "cautious a driver are you",
  "likely to be involved in accident, compared with others",
  "likely are you to avoid risky driving situations",
  "skilled a driver are you",
  "compare driving to others drivers in general",
  "Kind of driver - attentive/inattentive",
  "Kind of driver - careful/careless",
  "Kind of driver - decisive/indecisive",
  "Kind of driver - experienced/inexperienced",
  "Kind of driver - irritable/placid",
  "Kind of driver - nervous/confident",
  "Kind of driver - patient/impatient",
  "Kind of driver - responsible/irresponsible",
  "Kind of driver - safe/risky",
  "Kind of driver - selfish/considerate",
  "Kind of driver - slow/fast",
  "Kind of driver - tolerant/intolerant")

NetworkDataNamesDBQSI <- c(dbqnames,varnamesSelfImage)
labelsDBQSI <- c(namesNGroups$names, paste("si", seq(1:length(varnamesSelfImage)), sep = ""))
NetworkGroupsDBQSI <- c(namesNGroups$groups,rep("Self-Image", length(varnamesSelfImage)))

NetworkCorDBQSI <- cor_auto(NetworkDataDBQSI, ordinalLevelMax = 7)

CrossSectNetworkDBQSI <- qgraph(NetworkCorDBQSI, graph = "glasso", sampleSize = nrow(NetworkDataDBQSI), layout = "spring",
                                vsize = 3, nodeNames = NetworkDataNamesDBQSI, groups = NetworkGroupsDBQSI, legend.cex = .3,
                                labels = labelsDBQSI, color = c("#00ff00", "#00B2FFFF", "#FF0000FF"))

############### Create Supplementary Figure S5: ###############
pdf("CrossSectNetworkDBQSI.pdf", width = 10)
makeBW(CrossSectNetworkDBQSI)
dev.off()
###############################################################
includesSI <- c("V00174", "V00304")
includesSINames <- varnamesSelfImage[which(names(SIs) %in% includesSI)]
SIsInclude <- SIs[includesSI]


############################################################################################
######################### Bootstrap analyses, cross-sectional model#########################
############################################################################################

CrossSectNetForBootnet <- estimateNetwork(NetworkDataDBQ_BGVars, default = "EBICglasso")

set.seed(2018)
CrossSectEdgeWeightAccuracy <- bootnet(CrossSectNetForBootnet, nBoots = 2500)
plot(CrossSectEdgeWeightAccuracy, labels = FALSE, order = "sample")

set.seed(2018)
CrossSectCentralityStability <- bootnet(CrossSectNetForBootnet, nBoots = 2500, type = "case")
plot(CrossSectCentralityStability)

# Calculate CS-index:
corStability(CrossSectCentralityStability)
plot(CentralityStability,perNode = TRUE) 

######################### Create Supplementary Figure S2: ###########################################################################
pdf("CrossSect_Accuracy_Stability.pdf", width = 12)
grid.arrange(
  plot(CrossSectEdgeWeightAccuracy, order = "sample", labels = FALSE, samplelwd = .5, panels = FALSE),
  plot(CrossSectCentralityStability),
  ncol = 2)
dev.off()
#####################################################################################################################################
# Cross-sectional network in summary: All indices can be interpreted

################################################################################################################################################################