################################################################################################### # Creating simulated networks ################################################################################################### require(igraph) # GLOBAL VARIABLES GLOBAL_ADJ_BEFORE <- NULL # Adjacency matrix before adding noise GLOBAL_ADJ_AFTER <- NULL # Adjacency matrix after adding noise GLOBAL_ADJ_BEFORE_D <- NULL # Directed versions of adjacency matrix GLOBAL_ADJ_AFTER_D <- NULL # Directed versions of adjacency matrix ################################################################################################### # Main Simulations # mainSimU: # Main simulation function to run the undirected version a certain time # # Inputs: # n: number of times to run the simulation mainSimU <- function(n=1) { all_data <- vector("list", n) # Vector with all simulation data sim_data <- vector("list", 4) # Vector for each graph for(i in 1:n) { t <- proc.time() # Timer myGraph <- SimGraphU(FALSE) # Don't display plots every time X <- myGraph[[1]] bridge_lst <- myGraph[[2]] myComm <- comm_struct(rownames(as_adjacency_matrix(X, type="both"))) sim_data[[1]] <- GLOBAL_ADJ_BEFORE sim_data[[2]] <- GLOBAL_ADJ_AFTER sim_data[[3]] <- myComm sim_data[[4]] <- bridge_lst all_data[[i]] <- sim_data print(proc.time() - t) } return(all_data) } # mainSimD: # Main simulation function to run the directed version a certain time # # Inputs: # n: number of times to run the simulation mainSimD <- function(n=1) { all_data <- vector("list", n) # Vector with all simulation data sim_data <- vector("list", 4) # Vector for each graph for(i in 1:n) { t <- proc.time() # Timer myGraph <- SimGraphD(FALSE) # Don't display plots every time X <- myGraph[[1]] bridge_lst <- myGraph[[2]] myComm <- comm_struct(rownames(as_adjacency_matrix(X, type="both"))) sim_data[[1]] <- GLOBAL_ADJ_BEFORE_D sim_data[[2]] <- GLOBAL_ADJ_AFTER_D sim_data[[3]] <- myComm sim_data[[4]] <- bridge_lst all_data[[i]] <- sim_data print(proc.time() - t) } return(all_data) } ################################################################################################### # Helper Functions # get_comm: # Given a node name, return the community as a scalar (i.e. a=1, b=2, etc.) # # Inputs: # v_name: name of node get_comm <- function(v_name) { com_name <- substr(v_name, 1, 1) # community name of node scalar <- strtoi(charToRaw(com_name)) - strtoi(charToRaw("a")) + 1 return(scalar) } # comm_struct: # Given a list of node names, return the community structure composed of # scalars representing which community node is part of # # Inputs: # v_names: list of node names comm_struct <- function(v_names) { scalar_lst <- lapply(v_names, get_comm) return(scalar_lst) } # Randomize number of vertices in each community # Community_list returns a list of sizes (5-20) of complete graphs # lst[[i]] returns the i^th item in the list community_list <- function(graph_num){ size_list <- vector("list", graph_num) for (i in 1:graph_num) { ComGraph <- sample(6:21, 1) #clusters would be of size 5-20, as Node 1 is the default bridge size_list[i]<-ComGraph } size_list } # todirected: # Returns a directed graph from an undirected graph as input # # Inputs: # X: An undirected graph todirected <- function(X, bridge_lst) { A <- as_adjacency_matrix(X, type="both") Y <- graph_from_adjacency_matrix(A, mode="directed") V(Y)[bridge_lst]$color <- "lightgreen" return(Y) } # graphstructure: # Returns an undirected graph of communities and bridges given size dimensions as inputs # # Inputs: # Com_size number of communities in the graph # size_list list of nodes in each community # bridge_size number of bridges that connect to communities graphstructure <- function(Com_size, size_list, bridge_size) { # Define variables X <- make_empty_graph(n=0, directed = FALSE) # Return graph bridge_lst <- NULL # List of bridges in graph Com_graph_lst <- lapply(size_list, make_full_graph) # Complete graphs of communities sameCom <- 0.75 # Same community ratio diffCom <- 0.25 # Different community ratio # Label vertices in the communities, and add to graph for (i in 1:Com_size) { Com_graph_lst[[i]] = set_vertex_attr(Com_graph_lst[[i]], "name", value=paste(letters[i], 1:size_list[[i]], sep="")) X <- X + Com_graph_lst[[i]] } # Add bridges for (i in 1:bridge_size) { my_com <- sample(1:Com_size, 1) # Pick community randomly my_vertex <- paste(letters[my_com], size_list[[my_com]] + 1, sep="") print(paste("Bridge ", my_vertex, " belongs to community ", letters[my_com], sep="")) # Add vertex X <- add_vertices(X, 1, name=my_vertex) bridge_lst <- c(bridge_lst, my_vertex) # Add edges to all communities for (j in 1:Com_size) { my_comSize <- size_list[[j]] # Size of community bridge_vertices <- NULL # Corresponding vertices that connects with bridge bridge_edges <- NULL # Corresponding edges # Community bridge belongs to if (j == my_com) { bridge_vertices <- sample(1:my_comSize, ceiling(sameCom * my_comSize)) # Not community bridge belongs to } else { bridge_vertices <- sample(1:my_comSize, floor(diffCom * my_comSize)) } for (v in bridge_vertices) { c_vertex <- paste(letters[j], v, sep="") bridge_edges <- c(bridge_edges, my_vertex, c_vertex) } # Add Edges X <- add_edges(X, bridge_edges) print(paste("-- ", my_vertex, ": Added ", length(bridge_vertices), " edges to community ", letters[j], "(", round(length(bridge_vertices)/size_list[[j]] * 100, digits=2), "%)", sep="")) } # Make new vertex available size_list[[my_com]] <- size_list[[my_com]] + 1 } # Print Community sizes print("New community sizes:") for (i in 1:Com_size) { print(paste("Community ", letters[i], ": ", size_list[[i]], sep="")) } V(X)[bridge_lst]$color <- "lightgreen" X <- simplify(X) return(list(X, bridge_lst)) } # isSameCom: # Returns true if v and w belong to the same community, otherwise false # # Inputs: # v: name/label of first node (string type) # w: name/label of second node (string type) isSameCom <- function(v, w) { comV <- substr(v, 1, 1) # community name of node v comW <- substr(w, 1, 1) # community name of node w result <- FALSE if (comV == comW) { result <- TRUE } return(result) } # WeightedGraphU: # Adds weight to undirected input graph X, and returns new graph with weights. # # Inputs: # X: Graph without weights added # bridge_lst: List of bridges WeightedGraphU <- function(X, bridge_lst) { A <- as_adjacency_matrix(X, type="both") m <- nrow(A) n <- ncol(A) rnames <- rownames(A) cnames <- colnames(A) # Assign weight to each edge for (i in 1:m) { for (j in 1:i) { # Skip if self-loop or no edge currently exists if ((i == j) || (A[i,j] == 0)) { next } r <- rnames[i] c <- cnames[j] # Case 1: Edge belongs to bridge and community it belongs to if ((r %in% bridge_lst) && isSameCom(r, c)) { A[i,j] <- runif(1, 0.5, 1) # Random number from 0.5 to 1 } # Case 2: Edge belongs to bridge and community it doesn't belong to else if ((r %in% bridge_lst) && !isSameCom(r, c)) { A[i,j] <- runif(1, 0, 0.5) # Random number from 0 to 0.5 # Less than 0.1, make 0 if (A[i,j] < 0.1) { A[i,j] <- 0 } } # Case 3: Edge between non-bridge vertices else { A[i,j] <- runif(1, 0, 1) # Random number from 0 and 1 # Less than 0.1, make 0 if (A[i,j] < 0.1) { A[i,j] <- 0 } } # "Reflect" edge weights across adjacency matrix A[j,i] <- A[i,j] } } assign("GLOBAL_ADJ_BEFORE", A, envir=.GlobalEnv) # Add noise for (i in 1:m) { for (j in 1:i) { if (i != j) { A[i,j] <- abs(A[i,j] + rnorm(1, mean = 0, sd = 0.05)) # Add noise A[j,i] <- A[i,j] # Reflect edge weights } } } assign("GLOBAL_ADJ_AFTER", A, envir=.GlobalEnv) X <- graph_from_adjacency_matrix(A, mode="undirected", weighted=TRUE) V(X)[bridge_lst]$color <- "lightgreen" return(X) } # WeightedGraphD: # Adds weight to directed input graph X, and returns new graph with weights. # # Inputs: # X: Graph without weights added # bridge_lst: List of bridges WeightedGraphD <- function(X, bridge_lst) { A <- as_adjacency_matrix(X, type="both") m <- nrow(A) n <- ncol(A) rnames <- rownames(A) cnames <- colnames(A) # Assign weight to each edge for (i in 1:m) { for (j in 1:n) { # Skip if self-loop or no edge currently exists if ((i == j) || (A[i,j] == 0)) { next } r <- rnames[i] c <- cnames[j] # Case 1: Edge belongs to bridge and community it belongs to if (((r %in% bridge_lst) && isSameCom(r, c)) || ((c %in% bridge_lst) && isSameCom(c, r))) { A[i,j] <- runif(1, 0.5, 1) # Random number from 0.5 to 1 } # Case 2: Edge belongs to bridge and community it doesn't belong to else if (((r %in% bridge_lst) && !isSameCom(r, c)) || ((c %in% bridge_lst) && !isSameCom(c, r))){ A[i,j] <- runif(1, 0, 0.5) # Random number from 0 to 0.5 # Less than 0.1, make 0 if (A[i,j] < 0.1) { A[i,j] <- 0 } } # Case 3: Edge between non-bridge vertices else { A[i,j] <- runif(1, 0, 1) # Random number from 0 and 1 # Less than 0.1, make 0 if (A[i,j] < 0.1) { A[i,j] <- 0 } } } } assign("GLOBAL_ADJ_BEFORE_D", A, envir=.GlobalEnv) for (i in 1:m) { for (j in 1:n) { if (i != j) { A[i,j] <- abs(A[i,j] + rnorm(1, mean = 0, sd = 0.05)) # Add noise } } } assign("GLOBAL_ADJ_AFTER_D", A, envir=.GlobalEnv) X <- graph_from_adjacency_matrix(A, mode="directed", weighted=TRUE) V(X)[bridge_lst]$color <- "lightgreen" return(X) } # SimGraphU: # Returns an undirected graph of communities and bridges # # Inputs: # showPlot: Displays plot to screen if default or true SimGraphU <- function(showPlot=TRUE) { # Determine graph attributes Com_size <- sample(2:5, 1) # Number of communites in network size_list <- community_list(Com_size) # Size of each community bridge_size <- sample(1:5, 1) # Number of bridges # Print graph attributes print(paste("Graphing with ", Com_size, " communities and ", bridge_size, " bridges", sep="")) for (i in 1:Com_size) { print(paste("Community ", letters[i], ": ", size_list[[i]], sep="")) } # Graph and plot (# WARNING: size_list has changed due to the addition of bridges) myGraph <- graphstructure(Com_size, size_list, bridge_size) X <- myGraph[[1]] bridge_lst <- myGraph[[2]] # Change edge weights X <- WeightedGraphU(X, bridge_lst) if (showPlot == TRUE) { l <- layout_with_fr(X) plot(X, layout=l) } return(list(X, bridge_lst)) } # SimGraphD: # Returns a directed graph of communities and bridges # # Inputs: # None SimGraphD <- function(showPlot=TRUE) { # Determine graph attributes Com_size <- sample(2:5, 1) # Number of communites in network size_list <- community_list(Com_size) # Size of each community bridge_size <- sample(1:5, 1) # Number of bridges # Print graph attributes print(paste("Graphing with ", Com_size, " communities and ", bridge_size, " bridges", sep="")) for (i in 1:Com_size) { print(paste("Community ", letters[i], ": ", size_list[[i]], sep="")) } # Graph and plot (# WARNING: size_list has changed due to the addition of bridges) myGraph <- graphstructure(Com_size, size_list, bridge_size) X <- todirected(myGraph[[1]], myGraph[[2]]) # Convert X to a directed graph bridge_lst <- myGraph[[2]] # Change edge weights X <- WeightedGraphD(X, bridge_lst) if (showPlot == TRUE) { l <- layout_with_fr(X) plot(X, layout=l, edge.arrow.size=0.1) } return(list(X, bridge_lst)) } ## Run the functions to generate simulated networks (takes a long time): #undir_simulations <- mainSimU(500) #dir_simulations <- mainSimD(500) ################################################################################################### # Checking the simulations for sensitivity & specificity ################################################################################################### require(networktools) require(igraph) ################################# ## What is contained in each output: out <- mainSimU(1) as.matrix(out[[1]][[1]]) ## adjacency matrix (no noise) as.matrix(out[[1]][[2]]) ## adjacency matrix (noise) unlist(out[[1]][[3]]) ## community structure out[[1]][[4]] ## true bridges ## Create a representative plot: net <- as.matrix(undir_simulations[[8]][[1]]) bridges <- undir_simulations[[8]][[4]] bridges dim(net) colnames(net) %in% bridges net2 <- as.matrix(undir_simulations[[8]][[2]]) op <- par(mfrow=c(1,2)) qgraph(net, layout="spring", groups=list("Nodes"=c(1:68), "Bridge Nodes"=c(69:73)), colors=c("lightblue", "lightsalmon"), legend=FALSE) qgraph(net2, layout="spring", groups=list("Nodes"=c(1:68), "Bridge Nodes"=c(69:73)), colors=c("lightblue", "lightsalmon"), legend=FALSE) par(op) ################################# ## Necessary functions # "not in" function '%!in%' <- function(x,y)!('%in%'(x,y)) # returns sensitivity and specificity, given character vectors as input sens_spec <- function(all, predicted, actual) { predicted_negatives <- all[all %!in% predicted] actual_negatives <- all[all %!in% actual] true_pos <- sum(predicted %in% actual) true_neg <- sum(predicted_negatives %in% actual_negatives) # sensitivity = true positives over total positives sensitivity <- true_pos/length(actual) # specificity = true negatives over total negatives specificity <- true_neg/length(actual_negatives) return(c(sensitivity, specificity)) } # given centrality values, determines predicted bridges (combined criteria, & separate) predict_bridge <- function(strength, between, close) { strength_c <- names(strength[strength>quantile(strength, probs=0.75, na.rm=TRUE)]) between_c <- names(between[between>quantile(between, probs=0.75, na.rm=TRUE)]) close_c <- names(close[close>quantile(close, probs=0.75, na.rm=TRUE)]) pred_combined <- intersect(intersect(strength_c,between_c),close_c) pred_strength <- names(strength[strength>quantile(strength, probs=0.80, na.rm=TRUE)]) pred_between <- names(between[between>quantile(between, probs=0.80, na.rm=TRUE)]) pred_close <- names(close[close>quantile(close, probs=0.80, na.rm=TRUE)]) return(list(pred_combined=pred_combined,pred_strength=pred_strength,pred_between=pred_between,pred_close=pred_close)) } # returns sensitivity and specificity for each of several listed cases check_bridge <- function(x, directed=FALSE) { out_table <- matrix(rep(NA,32),16,2) try({ rownames(out_table) <- c("1.1","1.2","1.3","1.4","2.1","2.2","2.3","2.4","3.1","3.2","3.3","3.4","4.1","4.2","4.3","4.4") colnames(out_table) <- c("Sensitivity", "Specificity") all <- colnames(as.matrix(x[[1]])) actual <- x[[4]] ## Case 1 = communities prespecified, no noise b1 <- bridge(as.matrix(x[[1]]), communities= unlist(x[[3]]), directed=directed) # w/communities strength <- b1$'Bridge Strength'; between <- b1$'Bridge Betweenness'; close <- b1$'Bridge Closeness' pred_list <- predict_bridge(strength, between, close) ### 1.1 = combined criteria out_table[1, 1:2] <- sens_spec(all, pred_list[[1]], actual) ### 1.2 = strength only out_table[2, 1:2] <- sens_spec(all, pred_list[[2]], actual) ### 1.3 = betweenness only out_table[3, 1:2] <- sens_spec(all, pred_list[[3]], actual) ### 1.4 = closeness only out_table[4, 1:2] <- sens_spec(all, pred_list[[4]], actual) ## Case 2 = communities not specified, no noise b2 <- bridge(as.matrix(x[[1]]), directed=directed) strength <- b2$'Bridge Strength'; between <- b2$'Bridge Betweenness'; close <- b2$'Bridge Closeness' pred_list <- predict_bridge(strength, between, close) out_table[5, 1:2] <- sens_spec(all, pred_list[[1]], actual) out_table[6, 1:2] <- sens_spec(all, pred_list[[2]], actual) out_table[7, 1:2] <- sens_spec(all, pred_list[[3]], actual) out_table[8, 1:2] <- sens_spec(all, pred_list[[4]], actual) ## Case 3 = communities prespecified, noise b3 <- bridge(as.matrix(x[[2]]), communities= unlist(x[[3]]), directed=directed) strength <- b3$'Bridge Strength'; between <- b3$'Bridge Betweenness'; close <- b3$'Bridge Closeness' pred_list <- predict_bridge(strength, between, close) out_table[9, 1:2] <- sens_spec(all, pred_list[[1]], actual) out_table[10, 1:2] <- sens_spec(all, pred_list[[2]], actual) out_table[11, 1:2] <- sens_spec(all, pred_list[[3]], actual) out_table[12, 1:2] <- sens_spec(all, pred_list[[4]], actual) ## Case 4 = communities not specified, noise b4 <- bridge(as.matrix(x[[2]]), directed=directed) strength <- b4$'Bridge Strength'; between <- b4$'Bridge Betweenness'; close <- b4$'Bridge Closeness' pred_list <- predict_bridge(strength, between, close) out_table[13, 1:2] <- sens_spec(all, pred_list[[1]], actual) out_table[14, 1:2] <- sens_spec(all, pred_list[[2]], actual) out_table[15, 1:2] <- sens_spec(all, pred_list[[3]], actual) out_table[16, 1:2] <- sens_spec(all, pred_list[[4]], actual) }) return(out_table) } ################################# ## Testing simulations (takes a long time) undir_tests <- lapply(undir_simulations, FUN=check_bridge) dir_tests <- lapply(undir_simulations, FUN=check_bridge, directed=TRUE) ################################# ## Averages, full table av_mat_undir <- matrix(rep(NA, 32), 16,2) av_mat_dir <- matrix(rep(NA, 32), 16,2) for(i in 1:16) { for(j in 1:2){ vec1 <- sapply(undir_tests, function(l) l[i,j]) av_mat_undir[i,j] <- mean(na.omit(vec1)) vec2 <- sapply(dir_tests, function(l) l[i,j]) av_mat_dir[i,j] <- mean(na.omit(vec2)) } } av_mat_undir av_mat_dir fulltable <- as.data.frame(matrix(rep(NA, 32*6), 32, 6)) colnames(fulltable) <- c("Directed?", "Comm. Prespecified?", "Noise?", "Criteria", "Sensitivity", "Specificity") fulltable[1:16,5:6] <- av_mat_undir fulltable[17:32, 5:6] <- av_mat_dir fulltable$'Directed?' <- c(rep("Undirected", 16), rep("Directed", 16)) fulltable$'Criteria' <- rep(c("Combined", "Strength", "Betweenness", "Closeness"),8) fulltable$'Comm. Prespecified?' <- rep(c(rep(c("Yes"), 4), rep(c("No"), 4)), 4) fulltable$'Noise?' <- rep(c(rep("No", 8), rep("Yes", 8)), 2) fulltable <- fulltable[fulltable$Criteria %!in% c("Combined"),] write.csv(fulltable, file="Sensitivity Specificity, Positive Networks.csv") ## Averages of averages fulltable2 <- matrix(rep(NA, 20),10,2) fulltable2[1,1]<-mean(fulltable$Sensitivity) fulltable2[1,2]<-mean(fulltable$Specificity) fulltable2[2,1]<-mean(fulltable$Sensitivity[fulltable$`Directed?`=="Undirected"]) fulltable2[2,2]<-mean(fulltable$Specificity[fulltable$`Directed?`=="Undirected"]) fulltable2[3,1]<-mean(fulltable$Sensitivity[fulltable$`Directed?`=="Directed"]) fulltable2[3,2]<-mean(fulltable$Specificity[fulltable$`Directed?`=="Directed"]) fulltable2[4,1]<-mean(fulltable$Sensitivity[fulltable$`Comm. Prespecified?`=="Yes"]) fulltable2[4,2]<-mean(fulltable$Specificity[fulltable$`Comm. Prespecified?`=="Yes"]) fulltable2[5,1]<-mean(fulltable$Sensitivity[fulltable$`Comm. Prespecified?`=="No"]) fulltable2[5,2]<-mean(fulltable$Specificity[fulltable$`Comm. Prespecified?`=="No"]) fulltable2[6,1]<-mean(fulltable$Sensitivity[fulltable$`Noise?`=="Yes"]) fulltable2[6,2]<-mean(fulltable$Specificity[fulltable$`Noise?`=="Yes"]) fulltable2[7,1]<-mean(fulltable$Sensitivity[fulltable$`Noise?`=="No"]) fulltable2[7,2]<-mean(fulltable$Specificity[fulltable$`Noise?`=="No"]) fulltable2[8,1]<-mean(fulltable$Sensitivity[fulltable$Criteria=="Strength"]) fulltable2[8,2]<-mean(fulltable$Specificity[fulltable$Criteria=="Strength"]) fulltable2[9,1]<-mean(fulltable$Sensitivity[fulltable$Criteria=="Betweenness"]) fulltable2[9,2]<-mean(fulltable$Specificity[fulltable$Criteria=="Betweenness"]) fulltable2[10,1]<-mean(fulltable$Sensitivity[fulltable$Criteria=="Closeness"]) fulltable2[10,2]<-mean(fulltable$Specificity[fulltable$Criteria=="Closeness"]) fulltable2 ############################################################## ## Simulations with negative edges (bridge expected influence) ei_undirected <- mainSimU(500) ei_directed <- mainSimD(500) ## Turn some of those edges into negatives makesomeneg <- function(x, percent=80){ noiseless <- as.matrix(x[[1]]) noisy <- as.matrix(x[[1]]) for(i in 1:dim(noisy)[1]){ for(j in 1:dim(noisy)[2]){ noisy[i,j] <- makesomeneg2(noisy[i,j], percent) noiseless[i,j] <- makesomeneg2(noiseless[i,j], percent) } } res <- x res[[1]] <- noiseless res[[2]] <- noisy return(res) } makesomeneg2 <- function(x, percent=80){ if(sample(c(0:100),1)>percent){ x <- -1*x }else{} return(x) } undir_sim_neg <- lapply(ei_undirected, FUN=makesomeneg) dir_sim_neg <- lapply(ei_undirected, FUN=makesomeneg) ## Change prediction algorithm so "strength" now means EI1 ## and "between" now means EI2 check_bridge2 <- function(x, directed=FALSE) { out_table <- matrix(rep(NA,32),16,2) try({ rownames(out_table) <- c("1.1","1.2","1.3","1.4","2.1","2.2","2.3","2.4","3.1","3.2","3.3","3.4","4.1","4.2","4.3","4.4") colnames(out_table) <- c("Sensitivity", "Specificity") all <- colnames(as.matrix(x[[1]])) actual <- x[[4]] ## Case 1 = communities prespecified, no noise b1 <- bridge(as.matrix(x[[1]]), communities= unlist(x[[3]]), directed=directed) # w/communities strength <- b1$`Bridge Expected Influence (1-step)`; between <- b1$`Bridge Expected Influence (2-step)`; close <- b1$'Bridge Closeness' pred_list <- predict_bridge(strength, between, close) ### 1.1 = combined criteria out_table[1, 1:2] <- sens_spec(all, pred_list[[1]], actual) ### 1.2 = strength only out_table[2, 1:2] <- sens_spec(all, pred_list[[2]], actual) ### 1.3 = betweenness only out_table[3, 1:2] <- sens_spec(all, pred_list[[3]], actual) ### 1.4 = closeness only out_table[4, 1:2] <- sens_spec(all, pred_list[[4]], actual) ## Case 2 = communities not specified, no noise b2 <- bridge(as.matrix(x[[1]]), directed=directed) strength <- b2$`Bridge Expected Influence (1-step)`; between <- b2$`Bridge Expected Influence (2-step)`; close <- b2$'Bridge Closeness' pred_list <- predict_bridge(strength, between, close) out_table[5, 1:2] <- sens_spec(all, pred_list[[1]], actual) out_table[6, 1:2] <- sens_spec(all, pred_list[[2]], actual) out_table[7, 1:2] <- sens_spec(all, pred_list[[3]], actual) out_table[8, 1:2] <- sens_spec(all, pred_list[[4]], actual) ## Case 3 = communities prespecified, noise b3 <- bridge(as.matrix(x[[2]]), communities= unlist(x[[3]]), directed=directed) strength <- b3$`Bridge Expected Influence (1-step)`; between <- b3$`Bridge Expected Influence (2-step)`; close <- b3$'Bridge Closeness' pred_list <- predict_bridge(strength, between, close) out_table[9, 1:2] <- sens_spec(all, pred_list[[1]], actual) out_table[10, 1:2] <- sens_spec(all, pred_list[[2]], actual) out_table[11, 1:2] <- sens_spec(all, pred_list[[3]], actual) out_table[12, 1:2] <- sens_spec(all, pred_list[[4]], actual) ## Case 4 = communities not specified, noise b4 <- bridge(as.matrix(x[[2]]), directed=directed) strength <- b4$`Bridge Expected Influence (1-step)`; between <- b4$`Bridge Expected Influence (2-step)`; close <- b4$'Bridge Closeness' pred_list <- predict_bridge(strength, between, close) out_table[13, 1:2] <- sens_spec(all, pred_list[[1]], actual) out_table[14, 1:2] <- sens_spec(all, pred_list[[2]], actual) out_table[15, 1:2] <- sens_spec(all, pred_list[[3]], actual) out_table[16, 1:2] <- sens_spec(all, pred_list[[4]], actual) }) return(out_table) } undir_tests <- lapply(undir_sim_neg, FUN=check_bridge2) dir_tests <- lapply(dir_sim_neg, FUN=check_bridge2, directed=TRUE) av_mat_undir <- matrix(rep(NA, 32), 16,2) av_mat_dir <- matrix(rep(NA, 32), 16,2) for(i in 1:16) { for(j in 1:2){ vec1 <- sapply(undir_tests, function(l) l[i,j]) av_mat_undir[i,j] <- mean(na.omit(vec1)) vec2 <- sapply(dir_tests, function(l) l[i,j]) av_mat_dir[i,j] <- mean(na.omit(vec2)) } } av_mat_undir av_mat_dir fulltable_EI <- as.data.frame(matrix(rep(NA, 32*6), 32, 6)) colnames(fulltable_EI) <- c("Directed?", "Comm. Prespecified?", "Noise?", "Criteria", "Sensitivity", "Specificity") fulltable_EI[1:16,5:6] <- av_mat_undir fulltable_EI[17:32, 5:6] <- av_mat_dir fulltable_EI$'Directed?' <- c(rep("Undirected", 16), rep("Directed", 16)) fulltable_EI$'Criteria' <- rep(c("Combined", "EI1", "EI2", "Closeness"),8) fulltable_EI$'Comm. Prespecified?' <- rep(c(rep(c("Yes"), 4), rep(c("No"), 4)), 4) fulltable_EI$'Noise?' <- rep(c(rep("No", 8), rep("Yes", 8)), 2) fulltable_EI <- fulltable_EI[fulltable_EI$Criteria %!in% c("Combined", "Closeness"),] ## Simulation results for bridge strength, betweenness, closeness fulltable2 ################################################################################################### # Study 2: simulation of contagion ################################################################################################### # Simulation 2 predict_bridge2 <- function(strength) { pred_strength <- names(head(sort(strength, decreasing=TRUE), 5)) return(pred_strength) } tensteps_bridge <- function(netdata, directed=F, noiseparam=2) { noisemat <- as.matrix(netdata[[noiseparam]]) b1 <- bridge(as.matrix(netdata[[2]]), communities= unlist(netdata[[3]]), directed=directed) strength <- b1$'Bridge Strength' b2 <- predict_bridge2(strength) values <- vector() for(i in 1:dim(noisemat)[1]) { if(netdata[[3]][[i]][1]==1) { values[i] <- 0.5 } else {values[i] <- 0} if(colnames(netdata[[2]])[i] %in% b2) {values[i] <-0} } step <- function(x) { out <- vector() for(i in 1:dim(noisemat)[1]) { out[i] <- x[i]+(sum(noisemat[-i,i] * x[-i])*0.01) if(colnames(netdata[[2]])[i] %in% b2) {out[i] <-0} } return(out) } values2 <- step(values) values3 <- step(values2) values4 <- step(values3) values5 <- step(values4) values6 <- step(values5) values7 <- step(values6) values8 <- step(values7) values9 <- step(values8) values10 <- step(values9) res1 <- data.frame(values, values2, values3, values4, values5, values6, values7, values8, values9, values10) outnodes <- colSums(res1[unlist(netdata[[3]]) != 1,]) return(outnodes) } tensteps_str <- function(netdata, directed=F, noiseparam=2) { noisemat <- as.matrix(netdata[[noiseparam]]) b1 <- centralityTable(noisemat) b1 <- b1[b1$measure=="Strength",] b2 <- b1$value names(b2) <- b1$node b2 <- names(head(sort(b2, decreasing=TRUE), 5)) values <- vector() for(i in 1:dim(noisemat)[1]) { if(netdata[[3]][[i]][1]==1) { values[i] <- 0.5 } else {values[i] <- 0} if(colnames(netdata[[2]])[i] %in% b2) {values[i] <-0} } step <- function(x) { out <- vector() for(i in 1:dim(noisemat)[1]) { out[i] <- x[i]+(sum(noisemat[-i,i] * x[-i])*0.01) if(colnames(netdata[[2]])[i] %in% b2) {out[i] <-0} } return(out) } values2 <- step(values) values3 <- step(values2) values4 <- step(values3) values5 <- step(values4) values6 <- step(values5) values7 <- step(values6) values8 <- step(values7) values9 <- step(values8) values10 <- step(values9) res1 <- data.frame(values, values2, values3, values4, values5, values6, values7, values8, values9, values10) outnodes <- colSums(res1[unlist(netdata[[3]]) != 1,]) return(outnodes) } tensteps_bet <- function(netdata, directed=F, noiseparam=2) { noisemat <- as.matrix(netdata[[noiseparam]]) b1 <- centralityTable(noisemat) b1 <- b1[b1$measure=="Betweenness",] b2 <- b1$value names(b2) <- b1$node b2 <- names(head(sort(b2, decreasing=TRUE), 5)) values <- vector() for(i in 1:dim(noisemat)[1]) { if(netdata[[3]][[i]][1]==1) { values[i] <- 0.5 } else {values[i] <- 0} if(colnames(netdata[[2]])[i] %in% b2) {values[i] <-0} } step <- function(x) { out <- vector() for(i in 1:dim(noisemat)[1]) { out[i] <- x[i]+(sum(noisemat[-i,i] * x[-i])*0.01) if(colnames(netdata[[2]])[i] %in% b2) {out[i] <-0} } return(out) } values2 <- step(values) values3 <- step(values2) values4 <- step(values3) values5 <- step(values4) values6 <- step(values5) values7 <- step(values6) values8 <- step(values7) values9 <- step(values8) values10 <- step(values9) res1 <- data.frame(values, values2, values3, values4, values5, values6, values7, values8, values9, values10) outnodes <- colSums(res1[unlist(netdata[[3]]) != 1,]) return(outnodes) } tensteps_free <- function(netdata, directed=F, noiseparam=2) { noisemat <- as.matrix(netdata[[noiseparam]]) values <- vector() for(i in 1:dim(noisemat)[1]) { if(netdata[[3]][[i]][1]==1) { values[i] <- 0.5 } else {values[i] <- 0} } step <- function(x) { out <- vector() for(i in 1:dim(noisemat)[1]) { out[i] <- x[i]+(sum(noisemat[-i,i] * x[-i])*0.01) } return(out) } values2 <- step(values) values3 <- step(values2) values4 <- step(values3) values5 <- step(values4) values6 <- step(values5) values7 <- step(values6) values8 <- step(values7) values9 <- step(values8) values10 <- step(values9) res1 <- data.frame(values, values2, values3, values4, values5, values6, values7, values8, values9, values10) outnodes <- colSums(res1[unlist(netdata[[3]]) != 1,]) return(outnodes) } ## Generate networks: undir_simulations <- mainSimU(500) dir_simulations <- mainSimD(500) ## Noiseless nobridges_undir_noiseless <- lapply(undir_simulations, FUN = tensteps_bridge, noiseparam=1) nostrength_undir_noiseless <- lapply(undir_simulations, FUN = tensteps_str, noiseparam=1) nobet_undir_noiseless <- lapply(undir_simulations, FUN = tensteps_bet, noiseparam=1) free_undir_noiseless <- lapply(undir_simulations, FUN = tensteps_free, noiseparam=1) nobridges_dir_noiseless <- lapply(dir_simulations, FUN = tensteps_bridge, noiseparam=1) nostrength_dir_noiseless <- lapply(dir_simulations, FUN = tensteps_str, noiseparam=1) nobet_dir_noiseless <- lapply(dir_simulations, FUN = tensteps_bet, noiseparam=1) free_dir_noiseless <- lapply(dir_simulations, FUN = tensteps_free, noiseparam=1) ## Find averages stepmeans <- function(x) { m <- do.call(rbind, x) out <- vector() for(i in 1:dim(m)[2]) { out[i] <- mean(m[,i]) } return(out) } nobridges_undir_mean <- stepmeans(nobridges_undir) nostrength_undir_mean <- stepmeans(nostrength_undir) nobet_undir_mean <- stepmeans(nobet_undir) free_undir_mean <- stepmeans(free_undir) nobridges_dir_mean <- stepmeans(nobridges_dir) nostrength_dir_mean <- stepmeans(nostrength_dir) nobet_dir_mean <- stepmeans(nobet_dir) free_dir_mean <- stepmeans(free_dir) nobridges_undir_noiseless_mean <- stepmeans(nobridges_undir_noiseless) nostrength_undir_noiseless_mean <- stepmeans(nostrength_undir_noiseless) nobet_undir_noiseless_mean <- stepmeans(nobet_undir_noiseless) free_undir_noiseless_mean <- stepmeans(free_undir_noiseless) nobridges_dir_noiseless_mean <- stepmeans(nobridges_dir_noiseless) nostrength_dir_noiseless_mean <- stepmeans(nostrength_dir_noiseless) nobet_dir_noiseless_mean <- stepmeans(nobet_dir_noiseless) free_dir_noiseless_mean <- stepmeans(free_dir_noiseless) ## Table table4 <- rbind(nobridges_undir_noiseless_mean, nostrength_undir_noiseless_mean, nobet_undir_noiseless_mean, nobridges_dir_noiseless_mean, nostrength_dir_noiseless_mean, nobet_dir_noiseless_mean, nobridges_undir_mean, nostrength_undir_mean, nobet_undir_mean, nobridges_dir_mean, nostrength_dir_mean, nobet_undir_mean) dirs <- c(rep("Undirected", 3), rep("Directed", 3), rep("Undirected", 3), rep("Directed", 3)) noise <- c(rep("Noiseless", 6), rep("Noisey", 6)) table5 <- cbind(noise, dirs, table4) table5 <- table5[,c(1:2, 5, 7, 12)] ## Plots x11() plot(1:10, free_undir_mean, ylim=c(0,1.5)) lines(free_undir_mean, col="green") lines(nobridges_undir_mean, col="purple") lines(nostrength_undir_mean, col="red") lines(nobet_undir_mean, col="blue") x11() plot(0, type="n", xlim=c(0, 10), ylim=c(0,2.5)) lines(free_dir_mean, col="green") lines(nobridges_dir_mean, col="purple") lines(nostrength_dir_mean, col="red") lines(nobet_dir_mean, col="blue") x11() plot(0, type="n", xlim=c(0, 10), ylim=c(0,0.1)) lines(nobridges_undir_noiseless_mean, col="green") lines(nostrength_undir_noiseless_mean, col="purple") lines(nobet_undir_noiseless_mean, col="red") x11() plot(0, type="n", xlim=c(0, 10), ylim=c(0,0.1)) lines(nobridges_dir_noiseless_mean, lty="solid") lines(nostrength_dir_noiseless_mean, lty="dashed") lines(nobet_dir_noiseless_mean, lty="dotted") x11() plot(0, type="n", xlim=c(1, 10), ylim=c(0,0.09), xlab="Iterations", ylab="Total contagion") lines(nobridges_undir_noiseless_mean, lty="solid") lines(nostrength_undir_noiseless_mean, lty="dashed") lines(nobet_undir_noiseless_mean, lty="dotted") legend(x=1, y=0.09, legend=c("Strength", "Betweenness", "Bridge strength"), lty=c("dashed", "dotted", "solid")) ################################################################################################### # Study 3: empirical analyses ################################################################################################### ## Note: we cannot provide full code, as not all of the original data is ours to share. ## Therefore, we provide the code for the general process we applied require(networktools) require(igraph) require(qgraph) ## Generic function for producing plot and bridges ## Index is 1 strength, 2 between, 3 close, 4 EI, 5 EI2 full_bridge_workup <- function(data, index, communities, file){ w1 <- bridge(data, communities=communities) pdf(paste("BridgeCentrality_", file), width=18, height=12) plot(bridge(data, communities=communities), order="value") dev.off() bridge_w1 <- names(w1[[index]][w1[[index]]>quantile(w1[[index]], probs=0.80, na.rm=TRUE)]) bridge_num_w1 <- which(names(w1[[index]]) %in% bridge_w1) new_communities <- vector() for(i in 1:length(w1[[index]])) { if(i %in% bridge_num_w1) { new_communities[i] <- "Bridge" } else {new_communities[i] <- communities[i]} } new_groups <- list() for(i in 1:length(unique(new_communities))){ new_groups[[i]] <- which(new_communities %in% unique(new_communities)[i]) } names(new_groups) <- unique(new_communities) pdf(paste("BridgePlot_", file), width=24, height=24) qgraph(data, layout="spring", groups = new_groups, colors = c("lightsalmon", "lightblue", "indianred1", "lightgreen", 1:30)) dev.off() return(bridge_w1) } bridges <- full_bridge_workup(YOURNETWORKHERE) bridges ## ---- Robustness of networks ---- require(qgraph) ## Again we cannot provide full code, as not all of the original data is ours to share. ## Networks were simulated based on data from 18 networks in the literature. ## Therefore, we provide the code for the general process we applied ## Combine into lists adjs <- list(borsboom_1, borsboom_2, RobinaughWave1, RobinaughWave2, RobinaughWave3, RuzzanoClinicianPCR, RuzzanoSelfReport, RuzzanoPCalg, boschloo_network1, beard_network, boschloo_network2, adult_OCD_net1, adult_OCD_net3, adolescent_OCD_net1, adolescent_OCD_net2, adolescent_OCD_net3, levinson_adj1, levinson_adj2, levinson_adj3) comms <- list(comm_b_u, comm_b_u, Rob_comm, Rob_comm, Rob_comm_Wave3, Ruzzano_comm, Ruzzano_comm, Ruzzano_comm, comm_bn1, beard_comm, comm_bn2, OCD_comm, OCD_comm, OCD_comm, OCD_comm, OCD_comm, levinson_comm1, levinson_comm2, levinson_comm3) type <- list("Association", "PCalg", "Association", "Association", "Association", "Clinician", "Association", "PCalg", "eLASSO", "GLASSO", "eLASSO", "GLASSO","DAG", "Association", "GLASSO", "DAG", "GLASSO", "GLASSO", "GLASSO") adjs <- adjs[type %in% c("Association", "GLASSO", "eLASSO")] comms <- comms[type %in% c("Association", "GLASSO", "eLASSO")] type2 <- type[type %in% c("Association", "GLASSO", "eLASSO")] # Convert all to partial correlations diag(adjs[[5]]) <- 1 adjs_p <- adjs for(i in which(type2=="Association")){ adjs_p[[i]] <- cor2pcor(adjs[[i]]) } ## Use network structures to generate data of varying forms # n = 50, 100, 300, 500, 1000, 10000 # Partial correlations or Graphical LASSO require(bootnet) require(corpcor) require(networktools) generator <- ggmGenerator() ## this will generate data from a pcor matrix bridge_recovery_pcor <- function(pcornet, comms, n, method=c("GLASSO", "pcor")){ # Generate data data <- try(generator(n, pcornet)) # Re-estimate networks if(method[1]=="GLASSO"){ rec_net <- estimateNetwork(data, default="EBICglasso")$graph }else{ rec_net <- cor2pcor(cor(data)) } # Calculate "true" bridge centrality bridge_true <- bridge(pcornet, communities=comms) # Calculate recovered bridge centrality bridge_rec <- bridge(rec_net, communities=comms) # Correlate corvec <- vector() for(b in 1:5){ corvec[b] <- try(cor(bridge_true[[b]], bridge_rec[[b]])) } return(corvec) } bridge_recovery_pcor(adjs_p[[1]], comms=comms[[1]], n=50, method="GLASSO") ## Select appropriate psychometric networks adjs_p <- adjs_p[c(1:5,7,9:14)] comms <- comms[c(1:5,7,9:14)] ns <- rep(c(50, 100, 300, 500, 1000, 10000), 1200) networks <- rep(rep(1:12, 6)[order(rep(1:12, 6))], 100) it <- rep(1:100, 72) it <- it[order(it)] df2 <- data.frame(ns, networks, it, NA, NA, NA, NA, NA) colnames(df2) <- c("n", "network", "iteration", "strength", "betweenness", "closeness", "EI1", "EI2") set.seed(101) for(i in 1:dim(df2)[1]){ out <- try(bridge_recovery_pcor(adjs_p[[df2[i,"network"]]], comms= comms[[df2[i,"network"]]], n=df2[i,"n"], method="pcor")) if(class(out)!="try-error"){ df2[i,4:8] <- out } print(i) } df_nona <- na.omit(df2) boxplot_func <- function(type, color, title){ boxplot(df2[df2$n == 50, type], df2[df2$n == 100, type], df2[df2$n == 300, type], df2[df2$n == 500, type], df2[df2$n == 1000, type], df2[df2$n == 10000, type], names=c(50,100,300,500,1000,10000), col=color, main=title, outline=F, xlab="n", ylab="Correlation with True Network", ylim=c(-0.7,1)) } op <- par(mfrow=c(1,4)) boxplot_func("strength", "lightblue", "Bridge Strength") boxplot_func("betweenness", "indianred", "Bridge Betweenness") boxplot_func("closeness", "lightgreen", "Bridge Closeness") boxplot_func("EI1", "lightgrey", "Bridge Expected Influence") par(op)