library(igraph)
library(ggplot2)
library(network)
library(intergraph)
library(RColorBrewer)
library(dplyr)
library(visNetwork)
library(htmlwidgets)
library(greed)nodes <- read.csv("Facebook_att.csv", header = TRUE)
links <- read.csv("Facebook_edge.csv", header = TRUE)nodes[1:5,]## friend_count group mutual_friend_count relationship_status sex
## 1 73 Family 37 Married female
## 2 53 GraduateSchool 3 Married male
## 3 363 GraduateSchool 3 Married female
## 4 64 GraduateSchool 4 female
## 5 19 Family 10 female
## vertex.names NodeID
## 1 SE 1
## 2 PH 2
## 3 PD 3
## 4 LN 4
## 5 MB 5
# links?graph_from_data_frame
?asNetwork
facebook <- graph_from_data_frame(d = links, vertices = nodes$NodeID, directed = FALSE)
facebook_net <- asNetwork(facebook)
facebook_net## Network attributes:
## vertices = 93
## directed = FALSE
## hyper = FALSE
## loops = FALSE
## multiple = FALSE
## bipartite = FALSE
## total edges= 322
## missing edges= 0
## non-missing edges= 322
##
## Vertex attribute names:
## vertex.names
##
## No edge attributes
facebook_net %v% "vertex.names" <- nodes$vertex.names
facebook_net %v% "sex" <- nodes$sex
facebook_net %v% "friend_count" <- nodes$friend_count
facebook_net %v% "group" <- nodes$group
facebook_net %v% "mutual_friend_count" <- nodes$mutual_friend_count
facebook_net %v% "relationship_status" <- nodes$relationship_status#summary(facebook_net)plot(facebook_net, vertex.cex = 1.2, main = "Basic plot of Douglas’s Facebook friends")
Question 1.1: Check out the Summary and Plot, how many friends do Douglas have on Facebook? Is this a directed or undirected graph and why? What is the meaning of the link between nodes in the plot? [Question 1, 2 points].
Douglas has 93 friends.
This is a undirected graph because the edges do not have a direction.
It basically shows your friends, who are also friends between them
length(nodes$friend_count)## [1] 93
Question 1.2: Compare the degree, closeness and betweenness of Vertex 1 to the values of other nodes in the network. How will you evaluate the role of Vertex 1 in this network? [Question 2, 3 points].
The nod has 32 edges/friends.
The closeness centrality is low, meaning that that nod/person is not the closest one to all other persons, it costs more edges. Some other nod in the network has a mean shorter distance to other nodes in the network.
Has low betweenness centrality, meaning does not often form a bridge between nods/persons (does not frequently lie on the shortest path between nods in the network).
vertex_ID = 1
df_metrics = data.frame(
degree = degree(facebook, v = vertex_ID, mode = "total"),
closeness = closeness(facebook, v = vertex_ID, normalized = TRUE),
betweenness = betweenness(facebook, v = vertex_ID, directed = FALSE, normalized = TRUE)
)## Warning in closeness(facebook, v = vertex_ID, normalized = TRUE): At
## centrality.c:2874 :closeness centrality is not well-defined for disconnected
## graphs
df_metrics## degree closeness betweenness
## 1 32 0.02120793 0.1411131
Question 1.3: In the lecture, we introduced a few measures of centrality: degree, betweenness, eigenvector. Try to find top 5 nodes according to a) degree, b) betweenness, c) closeness and d) eigenvector. And develop the scatter plot between different metrics, you can refer to the code below. Discuss with your group and describe to your teacher, 1) how well does the top 5 nodes by different metrics overlap with each other? 2) why we need more than one metric to define centrality? [Question 3, 5 points]
library(scales)
vertex_ID = 1
df_metrics = data.frame(
nodeID = nodes$NodeID,
degree = degree(facebook, mode = "all"),
closeness = closeness(facebook, normalized = TRUE),
betweenness = betweenness(facebook, directed = FALSE, normalized = TRUE),
eig = eigen_centrality(facebook)$vector,
lcl = transitivity(facebook, type = "local")
)## Warning in closeness(facebook, normalized = TRUE): At
## centrality.c:2874 :closeness centrality is not well-defined for disconnected
## graphs
top_5 = df_metrics[order(-df_metrics$degree, -df_metrics$closeness, -df_metrics$betweenness, -df_metrics$eig),][1:5,]
coloring = ifelse(df_metrics$nodeID %in% top_5$nodeID, 'red', 'blue')
par(mfrow = c(2,2))
plot(df_metrics$degree, df_metrics$closeness, main = "Degree versus closeness", xlab = "Degree", ylab = "Closeness", col=alpha(coloring, 0.4))
plot(df_metrics$degree, df_metrics$betweenness, main = "Degree versus betweenness", xlab = "Degree", ylab = "Betweenness", col=alpha(coloring, 0.4))
plot(df_metrics$degree, df_metrics$eig, main = "Degree versus eigenvector", xlab = "Degree", ylab = "Eigenvector", col=alpha(coloring, 0.4))
plot(df_metrics$degree, df_metrics$lcl, main = "Degree versus local clustering", xlab = "Degree", ylab = "Local clustering", col=alpha(coloring, 0.4))
top_5## nodeID degree closeness betweenness eig lcl
## 1 1 32 0.02120793 0.141113068 1.0000000 0.3024194
## 12 12 21 0.02110576 0.013670660 0.9012656 0.5809524
## 23 23 19 0.02110576 0.015852965 0.8305933 0.6140351
## 31 31 19 0.02108641 0.002138961 0.9016029 0.7192982
## 26 26 18 0.02108158 0.001726770 0.8678475 0.7385621
- In degree nods 23 and 31 overlap.
- In closeness nods 23 and 31 overlap.
- In betweenness Node 1 has the highest and there is no overlap.
- In eigenvector nods don’t overlap.
Because Centrality describes who is a central actor and it might be handy to analyse in which context a actor has a central behavior.
Question 1.4: Discuss within your group on how you understand each of these measures. And describe to your teacher, 1) why diameter should be larger than 1, and other metrics such as edge density and transitivity are smaller than 1? 2) is this a tightly knitted network? [Question 4, 3 point].
deg_dist <- degree_distribution(facebook)
barplot(deg_dist)
tibble(
edge_density = edge_density(facebook),
diameter = diameter(facebook, directed = FALSE),
transitivity = transitivity(facebook, type = "global"),
centr_degree = centr_degree(facebook)$centralization
)## # A tibble: 1 × 4
## edge_density diameter transitivity centr_degree
## <dbl> <dbl> <dbl> <dbl>
## 1 0.0753 4 0.665 0.273
Because its a count of edges between the shortest path, and the number of edges is always N-1. So having diameter of 1 would lead to 1-1=0. Meaning no connection.
In edge density, knowing that degree lies between 0 and 1 which can be seen as a probability that a pair of nodes, picked uniformly at random from the whole network, is connected by an edge. Furthermore, there is low cohesiveness, explaining the quality of cooperation in the network.
In transitivity, 1 would imply that all nodes in a component are connected to all others which is not the case with 0.66.
In central degree, shows the centrality score based on node-level centrality. In this case the score is low meaning in the network there is less node-level centrality.
There is some degree of being tightly knit network due to components having closely linked nods.
Question 1.5: Reflect on what we have discussed about the Facebook network on Slide 42 of the lecture. Do you think this small network of Douglas resonates some general patterns of the entire Facebook network in terms of the components size and number? How can you explain such an observation [Question 5, 3 point].
barplot(components(facebook)$csize, names.arg=c(1:components(facebook)$no))
nodes['compID'] <- components(facebook)$membershipfactor(nodes[nodes$compID == 1, ]$group)['levels']## [1] <NA>
## Levels: BookClub Family GraduateSchool Music Work
nodes[nodes$compID == 7, ]$group## [1] "HighSchool" "HighSchool" "HighSchool" "HighSchool" "Music"
Question 1.6: From the plot you produced, do you find instances of intermingling of Douglas friends (i.e., belong to different groups but end up in the same component)? Do you find any isolated groups here? What can you conclude about the mixing of Douglas’ Facebook friends? [Question 6, 3 point].
group <- as.factor(get.vertex.attribute(facebook_net, attrname = 'group'))
pal <- brewer.pal(nlevels(group), "Set1")
plot(facebook_net, vertex.col = pal[group], vertex.cex = 1.5, main = "Plot of Facebook Data colored by friend type")
legend(x = "bottomleft", legend = levels(group), col = pal, pch = 19, pt.cex = 1.2, bty = "n", title = "Friend type")
Yes, in component 1 there are groups from family, book club, graduate school and music (don’t know maybe play in a band). Work is not isolated, because it is also occurs in component 1. The group spiel (brown) is isolated, because it exists in one component. The group college (blue) is isolated in a sense, it is isolated in different components but doesn’t overlap.
Yes, brown and blue are isolated.
About mixing we can conclude that he has a separated work life and personal life.
Question 1.7: Using the above code as a reference, check out the attributes of other factors (e.g., sex, relationship_status) in terms of people in the same components. Note that you can specify the ‘attrname’ parameter within the function ‘get.vertex.attribute’. Discuss with your group and describe to your teacher whether or not these factors are the keys in determining the formation of components [Question 7, 3 point].
group <- as.factor(get.vertex.attribute(facebook_net, attrname = 'sex'))
pal <- brewer.pal(nlevels(group), "Set1")## Warning in brewer.pal(nlevels(group), "Set1"): minimal value for n is 3, returning requested palette with 3 different levels
plot(facebook_net, vertex.col = pal[group], vertex.cex = 1.5, main = "Plot of Facebook Data colored by sex type")
legend(x = "bottomleft", legend = levels(group), col = pal, pch = 19, pt.cex = 1.2, bty = "n", title = "Sex type")
So looking at sex - it is not a key factor, it has not an effect on to determine the components. We have two levels, were the majority of the components have both sexes.
group <- as.factor(get.vertex.attribute(facebook_net, attrname = 'relationship_status'))
pal <- brewer.pal(nlevels(group), "Set1")
plot(facebook_net, vertex.col = pal[group], vertex.cex = 1.5, main = "Plot of Facebook Data colored by relationship_status type")
legend(x = "bottomleft", legend = levels(group), col = pal, pch = 19, pt.cex = 1.2, bty = "n", title = "relationship_status type")
The same idea as above applies as well for “relationship status”.
Question 1.8: From this analysis, what do you observe from the density values? Are they similar across different groups? What is the minimum and maximum value you observed here and how do you explain that? [Question 8, 4 point].
groups <- as.factor(get.vertex.attribute(facebook_net, attrname = 'group'))
for (x in levels(groups)){
y <- get.inducedSubgraph(facebook_net, which(facebook_net %v% "group" == x))
print(paste0("Density for ", x, " friends is ", edge_density(asIgraph(y))))
}## [1] "Density for BookClub friends is 1"
## [1] "Density for College friends is 0.19047619047619"
## [1] "Density for Family friends is 0.624505928853755"
## [1] "Density for GraduateSchool friends is 0.266666666666667"
## [1] "Density for HighSchool friends is 0.3"
## [1] "Density for Music friends is 0.316666666666667"
## [1] "Density for Spiel friends is 1"
## [1] "Density for Work friends is 0.3"
There are various densities across groups. Some groups are similar such as Work and Music and HighSchool, Bookclub and Spiel. The maximum is 1 and minimum 0.19.
Knowingly that density lies between the range 0 and 1, the sub-networks varies in everyone knows each other and only few know each other. For example, density of BookClub and Spiel is 1 because everyone knows each other whereas example Work and Music friends is 0.3 because only few of the friends know each other.
Question 1.9: To answer this question, search and discuss within your group on the theory of ‘community detection’. Describe to your teacher what community detection is, and why it is useful to understand complex networks [Question 9, 8 point].
Community detection is a analysis technique, sometimes called clustering. The goal of community detection is to find the natural divisions of a network into groups of nodes such that there are many edges within groups (tightly connected) and few edges between them (loosely connected). Or in other words, it helps us to reveal the hidden relations among the nodes in the network. Furthermore, community detection methods can be categorized into two types; Agglomerative Methods and Divisive Methods. In Agglomerative methods, edges are added one by one to a graph which only contains nodes. Edges are added from the stronger edge to the weaker edge. Whereas, Divisive methods follow the opposite of agglomerative methods. In there, edges are removed one by one from a complete graph. There exist several community detection algorithms.
Well, a complex network e.g. with a million or more nodes can rarely be visualized in its entirety. Such networks are simply too big to be represented usefully on the screen. If the nodes of the network divide naturally into groups, however, then we can make a simpler but still useful picture by representing each group as a single node and the connections between groups as edges. This simplified representation allows us to see the large-scale structure of the network without getting bogged down in the details of the individual nodes. In other words, it gives the possibility to analyse big networks from a bird-eye view by dividing it into groups naturally, to understand the connections within the groups and between the groups. For instance, applying to social media algorithms to discover people with common interests and keep them tightly connected.
Question 10: Compare the plots that you generate from the different algorithms; do you find them similar and Why? Do you think this observation (similar or not similar by using different community detection algorithms) holds in more complex network? [Question 10, 5 point].
cw <- cluster_walktrap(facebook)
plot(cw, facebook, vertex.label = V(facebook)$group,
main = "Walktrap")
ceb <- cluster_edge_betweenness(facebook)
plot(ceb, facebook, vertex.label = V(facebook)$group,
main = "Edge Betweenness")
cfg <- cluster_fast_greedy(facebook)
plot(cfg, facebook, vertex.label = V(facebook)$group,
main = "Fast Greedy")
clp <- cluster_label_prop(facebook)
plot(clp, facebook, vertex.label = V(facebook)$group,
main = "Label Prop")
cle <- cluster_leading_eigen(facebook)
plot(cle, facebook, vertex.label = V(facebook)$group,
main = "Leading Eigen")
- Label Prop and Edge Betweenness have both the same amount of communities.
- Leading Eigen and Edge Betweenness have both 4 communities with connections in between.
- Leading Eigen, Edge Betweenness, Walktrap, Fast Greedy and Label Prop have nods overlapping multi communities
- Leading Eigen and Walktrap have similar set of nodes in several communities
- Fast Greedy and Label Prop have similar set of nodes in several communities
Do you think this observation (similar or not similar by using different community detection algorithms) holds in more complex network?
No, the number of communities and overlap of nodes can differ when the current network complexity changes. It might be that new communities occur or nodes with duplicate communities characteristics leading to overlap.
Question 11: Discuss with your group, then describe to your teacher, under which circumstances, we might need to work on a synthetic network. [Question 11, 3 point].
- Mimic biological computation to obtain knowledge about how neurons process information
- To investigate communities on certain variables
- To investigate the growth of cities depending on certain factors
- To investigate the behavior of a malicious software within a network
Making simplified network and shifting the parameters and creating complex network by using simpler rules.
##Erdos-Renyi Model
Question 12: Plot your network (with n=100, and p=1/100) and compare with those with your group members. Are they identical? Explain why they are/aren’t [Question 12, 2 point].
ER <- sample_gnp(100, 1/100)
plot(ER, vertex.label= NA, edge.arrow.size=0.02,vertex.size = 0.5, xlab = "ER Random Network: G(N,p) model")
No, because the models are Random Graphs which means on each execution the position of the edges changes.
Question 13: Develop three networks with the same number of vertices (n), but different probability (p); Name them as ER1, ER2, and ER3. Develop the plots of ER1, ER2 and ER3, describe how these three graphs look differently as p increase and explain why. [Question 13, 3 point].
ER1 = sample_gnp(100, 1/100)
ER2 = sample_gnp(100, 1/50)
ER3 = sample_gnp(100, 1/25)
par(mfrow = c(2,2))
plot(ER1, vertex.label= NA, edge.arrow.size=0.02,vertex.size = 0.5, xlab = "ER1 Random Network: G(100, 1/100) model")
plot(ER2, vertex.label= NA, edge.arrow.size=0.02,vertex.size = 0.5, xlab = "ER2 Random Network: G(100, 1/50) model")
plot(ER3, vertex.label= NA, edge.arrow.size=0.02,vertex.size = 0.5, xlab = "ER2 Random Network: G(100, 1/25) model")
The higher the probability the more edges in the network. Knowing that p represents the fraction of the edges that are present from the possible N(N-1)/2, so the higher the p the higher the fraction of the edges that are present in the network.
Question 14: If p<1, for n great enough, what happens to the clustering coefficient of an ER random graph and why? (You can use the ‘transitivity’ function to test your guess). Discuss with your group and describe the answer to your teacher. [Question 14, 2 point].
proba_seq = c(1/100, 1/50, 1/25, 1)
N = 100
par(mfrow = c(2,2))
for(p in proba_seq){
ER4 = sample_gnp(N, p)
lcl = transitivity(ER4, type = "localaverage")
plot(ER4, vertex.label= NA, edge.arrow.size=0.02,vertex.size = 0.5, xlab = paste0("ER: G(", N, ", ", p, ") , mean lcl: ", round(lcl, 3)))
}
n_nodes = c(5, 10, 50, 100, 200, 300)
p_values = c(1/5, 1/10, 1/100)
mu_arr = p_values
avg_path_arr = p_values
for(idx in 1:length(p_values)){
# n = n_nodes[idx]
p = p_values[idx]
ER = sample_gnp(100, p)
mu_arr[idx] = transitivity(ER, type = "global")
avg_path_arr[idx] = mean_distance(ER)
}
par(mfrow = c(1,2))
plot(p_values, mu_arr, xlab='n nodes', ylab='mean local clustering coefficient', main='C over n nodes', type="b")
plot(p_values, avg_path_arr, xlab='n nodes', ylab='mean avg path length', main='L over n nodes', type="b")
Because, p represents the fraction of the edges, so the higher the p the more connections will be in the network which automatically leads to higher clustering coefficient. The expected value of a number of neighbors of a node i is equal to p*(N-1), meaning the proportion of the neighbors which are connected to a node is just p.
##Watts and Strotgatz Model (Small-world Random Graph Model)
Question 15: Check the clustering coefficient and average path length of the Regular, SW1, SW2 and SW3. Describe the trend of clustering coefficient and average path length as p increase. Does any of these graphs show the desirable attributes that you are looking for a small world network? [Question 15, 3 point].
Regular <- watts.strogatz.game(dim=1,size=300,nei=6, p=0)
plot(Regular, layout=layout.circle, vertex.label=NA, vertex.size=5, main= "Network with zero rewiring probability ")
SW1 <- watts.strogatz.game(dim=1,size=300,nei=6, p=0.001)
plot(SW1, layout=layout.circle, vertex.label=NA, vertex.size=5, main="Network with 0.001 rewiring probability ")
SW2 <- watts.strogatz.game(dim=1,size=300,nei=6, p=0.01)
plot(SW2, layout=layout.circle, vertex.label=NA, vertex.size=5, main="Network with 0.01 rewiring probability ")
SW3 <- watts.strogatz.game(dim=1,size=300,nei=6, p=0.1)
plot(SW3, layout=layout.circle, vertex.label=NA, vertex.size=5, main= "Network with 0.1 rewiring probability ")
par(mfrow = c(1,2))
wsgame = list(Regular, SW1, SW2, SW3)
lcl = sapply(wsgame, function(x){
transitivity(x, type = "global")
})
avgpath = sapply(wsgame, mean_distance)
p_values = c(0, 0.001, 0.01, 0.1)
plot(y=lcl, x=p_values, main='global coefficient', type="b")
plot(y=avgpath, x=p_values, main='average path length', type="b")
As p increases lcl and average path length decreases - lcl linearly and average path length exponentially.
Does any of these graphs show the desirable attributes that you are looking for a small world network?
Yes, knowingly that small-world has high clustering and small average shortest path length (low diameter) over p, the average shortest path decreases while clustering remains the same for a long time (but it decreases as well). So, clustering is hard to reduce but not the average shortest path length.
Question 16: For the same setting, i.e., size =300, nei=6, what are the range of p that you will suggest to build a small-world network and why? One solution you can consider is to refer to the Figure 2 in Watts and Strogatz (1998), which explains how small-world networks can be created from randomly rewired graphs (similar figure also appeared in Lecture 5).
Reproduce Figure 2 in the current context (i.e., size=300, nei=6) and find the right range of p in order to build a small-world network. Discuss with your group and show the answer to your teacher. [Question 16, 8 point].
proba_arr = seq(0.0001, 1, 0.01)
lcl_arr = proba_arr
avg_path_arr = proba_arr
Regular = watts.strogatz.game(dim=1,size=300,nei=6, p=0.0001)
r_lcl = transitivity(Regular, type="global")
r_avg_path = mean_distance(Regular)
for(idx in 1:length(proba_arr)){
proba = proba_arr[idx]
wsgame <- watts.strogatz.game(dim=1,size=300,nei=6, p=proba)
lcl_arr[idx] = transitivity(wsgame, type="global") / r_lcl
avg_path_arr[idx] = mean_distance(wsgame) / r_avg_path
}plot(proba_arr, lcl_arr, type='b', col='blue', ylab='relative change', xlab='probability')
lines(proba_arr, avg_path_arr, type='b', col='red')
legend("topright", cex=1, lty=1:2, legend=c("Global clustering","avg path length"), col=c("blue","red"))
x_range = c(0.001, 0.20)
plot(proba_arr, lcl_arr, type='b', col='blue', ylab="relative change", xlab='probability', xlim=x_range, main="C & L")
lines(proba_arr, avg_path_arr, type='b', col='red', xlim=x_range, main="L")
legend("topright", cex=1, lty=1:2, legend=c("C(p) / C(0)", "L(p) / L(0)"), col=c("blue","red"))
We assume that the right range of p values is from 0.02 until 0.05 because after 0.1 the avg path length remains constant. We don’t want clustering to decrease to much but the avg path length should be low. And within the chosen range clustering is still high and avg path length remaining at the lowest and for that reason we assume this range of p values is the right one.
Question 17: As a follow-up question of Q16, do you find the p values of very large or relatively small? Do you need to rewire significant amount of connections to make a network small-world-like? [Question 17, 2 point].
The p values are large which leads to difficulty understanding the right range of p values. Fortunately, the function xlim gives the possibility to zoom in an area of p values to closer evaluate.
According to the result no significant amount of connections required to make a network small-world like. In our range for instance we can use 5% of rewiring to result in a small-world like network.
In the paper of Watts and Strogatz (1998), they pointed out that the value of p has two important implications:
-“The idealized construction above reveals the key role of short cuts. It suggests that the small-world phenomenon might be common in sparse networks with many vertices, as even a tiny fraction of short cuts would suffice.”
-“Thus, infectious diseases are predicted to spread much more easily and quickly in a small world; the alarming and less obvious point is how few short cuts are needed to make the world small.”
Use your own words to explain these two implications to your teachers (you are suggested to read the paper for a better understanding). For the second implication, try to explain it in the context of COVID. [Question 18, 4 point].
-
It implies that by introducing shortcuts, referring the few long-range edges, it makes the hops less to connect to certain nods that would otherwise take more hops. This results in immediate drop in average path length L(P).
-
It implies that due to the low average path length L(P), referring to the short distance to neighbors, the spread of disease is faster due to the short distance needed to infect someone.
##Barabasi-Albert Model
Question 19: What does the power in the above function mean? How can it govern the structure of the network (e.g., the formulation of hubs)? (Hint: Change the value of power from 0.05, 0.5, 1, 1.5; See how the plot evolves; if you still fail to see the difference, visualize the vertex size according to the edge number, you can consider the code below.) [Question 19, 3 point].
?barabasi.game
par(mfrow = c(1,2))
power_values = c(0.05, 0.5, 1, 1.5)
for(power in power_values){
g0 <- barabasi.game(100, power = power, m = NULL, out.dist = NULL, out.seq = NULL, out.pref = FALSE, zero.appeal = 1, directed = FALSE,algorithm ="psumtree", start.graph = NULL)
title = paste0("SF network model, power=", power)
plot(g0, vertex.label= NA, edge.arrow.size=0.02,vertex.size =5, main = title)
}
par(mfrow = c(1,2))
power_values = c(0.05, 0.5, 1, 1.5)
for(power in power_values){
g0 <- barabasi.game(100, power = power, m = NULL, out.dist = NULL, out.seq = NULL, out.pref = FALSE, zero.appeal = 1, directed = FALSE,algorithm ="psumtree", start.graph = NULL)
title = paste0("SF Network model, power=", power)
g1Net<-asNetwork(g0)
VS = 3+ 0.5*degree(g0)
plot(g0, vertex.label= NA, edge.arrow.size=0.02,vertex.size =VS, main = title)
}
P[i] ~ k[i]^alpha + a Power: The power (alpha) of the preferential attachment, the default is one, ie. linear preferential attachment.
It influences the probability that an old nod will be chosen to get a edge created by the new nod. Basically, the probability that a edge of a new nod will connect to a old nod.
Gives the ability to decide how strong the old nods (with higher degree) are able to grab edges added to the network. So, the higher the power the stronger they get to grab the edges created by the new nods which leads to more hubs.
Question 20: Discuss with your groups, if you are maintaining a network with a power of 0.5 and 1.5, respectively, what will be your plans to build up resilience for 1) random attack, and 2) targeted attack? (the meanings of ‘random attack’ and ‘targeted attack’ are the same as what is mentioned in Lecture 6, scale-free network) [Question 20, 5 point].
For the network with both power, we would apply a algorithm which is capable of recovering the network using the alive nodes by connecting the links among the nodes in a self-organizing manner of communication.
Question 21: Download the csv data of these network from the BB, import the data to R and build the network. Check out their clustering coefficients, path lengths and degree distribution. Do you find these real networks show some attributes of the synthetic architectures we studied above (e.g., random ER graph, small-world, and scale-free network)? Show your teachers some numbers, plots and how you interpret the results. [Question 21, 8 point].
bnodes <- read.csv("Ex2_brightkite network.csv", header = TRUE)
gnodes <- read.csv("Ex2_network of General Relativity.csv", header = TRUE)# brightkite
uniqueID = seq(min(bnodes$FromNodeID), max(bnodes$ToNodeID), 1)
brightkite <- graph_from_data_frame(d = bnodes, vertices = uniqueID, directed = FALSE)
# general relative
uniqueID = seq(min(gnodes$FromNodeID), max(gnodes$ToNodeID), 1)
general_relative <- graph_from_data_frame(d = gnodes, vertices = uniqueID, directed = FALSE)Check out their clustering coefficients, path lengths and degree distribution.
# brightkite
diameters = c(16, 17)
n_nodes = c(58228, 26184)
n_name = c('brightkite', 'general relative')
networks = list(brightkite, general_relative)
C_list = sapply(networks, function(x){
transitivity(x, type = "global")
})
Degree_list = sapply(networks, function(x){
mean(degree_distribution(x))
})
tibble(
Name = n_name,
Nodes = n_nodes,
C = C_list,
Diameter = diameters,
Degree_D = Degree_list
)## # A tibble: 2 × 5
## Name Nodes C Diameter Degree_D
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 brightkite 58228 0.111 16 0.000441
## 2 general relative 26184 0.630 17 0.00613
Simulation to compare with real-networks!
nodes = c(1000, 1500, 5000, 8000)
for(n in nodes){
ER4 = sample_gnp(n, 1/25)
C = transitivity(ER4, type = "global")
Degree_D = mean(degree_distribution(ER4))
Diameter = diameter(ER4)
title = paste0("ER: nodes: ", n, " Cran: ", round(C, 4), " Diameterran: ", round(Diameter, 4), " Dran: ", round(Degree_D, 6))
print(title)
}## [1] "ER: nodes: 1000 Cran: 0.0392 Diameterran: 3 Dran: 0.016667"
## [1] "ER: nodes: 1500 Cran: 0.0399 Diameterran: 3 Dran: 0.011765"
## [1] "ER: nodes: 5000 Cran: 0.04 Diameterran: 3 Dran: 0.004098"
## [1] "ER: nodes: 8000 Cran: 0.04 Diameterran: 3 Dran: 0.002551"
# Small World
nodes = c(1000, 1500, 5000, 8000)
for(n in nodes){
SF = watts.strogatz.game(dim=1, size=n, nei=6, p=0.005)
C = transitivity(SF, type = "global")
Degree_D = mean(degree_distribution(SF))
Diameter = diameter(SF)
title = paste0("Small-World: nodes: ", n, " Cran: ", round(C, 4), " Diameterran: ", round(Diameter, 4), " Dran: ", round(Degree_D, 6))
print(title)
}## [1] "Small-World: nodes: 1000 Cran: 0.6604 Diameterran: 16 Dran: 0.066667"
## [1] "Small-World: nodes: 1500 Cran: 0.6614 Diameterran: 19 Dran: 0.0625"
## [1] "Small-World: nodes: 5000 Cran: 0.6607 Diameterran: 23 Dran: 0.066667"
## [1] "Small-World: nodes: 8000 Cran: 0.6593 Diameterran: 21 Dran: 0.066667"
# Scale Free
nodes = c(1000, 1500, 5000, 8000)
for(n in nodes){
BA <- barabasi.game(n, power = 1.4, m = NULL, out.dist = NULL, out.seq = NULL, out.pref = FALSE, zero.appeal = 1, directed = FALSE,algorithm ="psumtree", start.graph = NULL)
C = transitivity(BA, type = "global")
Degree_D = mean(degree_distribution(BA))
Diameter = diameter(BA)
title = paste0("Scale Free: nodes: ", n, " C: ", round(C, 4), " Diameterran: ", round(Diameter, 4), " Dran: ", round(Degree_D, 6))
print(title)
}## [1] "Scale Free: nodes: 1000 C: 0 Diameterran: 12 Dran: 0.003344"
## [1] "Scale Free: nodes: 1500 C: 0 Diameterran: 16 Dran: 0.001093"
## [1] "Scale Free: nodes: 5000 C: 0 Diameterran: 24 Dran: 0.000881"
## [1] "Scale Free: nodes: 8000 C: 0 Diameterran: 20 Dran: 0.000185"
The network General Relative is more close to Small-World network due to Cluster Coefficient and Diameter. And Brightkite network is according to the simulation results close to Scale Free network due to Diameter and Degree Distribution.
class2022 = readRDS("classnetwork_2022.rds")
class2022_matrix <- as_edgelist(class2022, names = TRUE)
edges <- data.frame(from = class2022_matrix[,1], to = class2022_matrix[,2])
classnetwork <- graph_from_data_frame(d = edges, vertices = 1:140, directed = FALSE)
nodes <- tibble(id = 1:140, C=transitivity(classnetwork, 'local'), color=ifelse(C > 0.6, "red", "lightblue"))
# saveWidget(visNetwork(nodes, edges, main = "Red nods have C > 0.6") %>%
# visIgraphLayout(layout = "layout_nicely", smooth = TRUE) %>%
# visNodes(shape="circle") %>%
# visOptions(highlightNearest = list(enabled = T, hover = T), nodesIdSelection = T), file = "class2022.html")
visNetwork(nodes, edges, main = "Red nods have C > 0.6") %>%
visIgraphLayout(layout = "layout_nicely", smooth = TRUE) %>%
visNodes(shape="circle")Question 22: Open the html file to view the network of this class. Do you find a lot of strong ties in this network? Highlight three of them, and show to your teachers. Then check the clustering coefficient of this network, and does it match your observation (i.e., this network has a few/ many strong ties)? [Question 22, 5 points].
classnetwork <- graph_from_data_frame(d = edges, vertices = nodes, directed = FALSE)
tibble(id = 1:140, C=transitivity(classnetwork, 'local'), D=degree(classnetwork)) %>%
arrange(-C, -D) %>% head(10)## # A tibble: 10 × 3
## id C D
## <int> <dbl> <dbl>
## 1 88 1 5
## 2 47 1 4
## 3 35 1 3
## 4 43 1 3
## 5 78 1 3
## 6 52 1 2
## 7 63 1 2
## 8 109 0.905 7
## 9 102 0.857 7
## 10 11 0.833 4
Do you find a lot of strong ties in this network? Highlight three of them, and show to your teachers.
Yes, there are certain amount of strong ties in the network, such as nodes: 88, 109 and 102, see html file or table above. Furthermore, this network has few strong ties when we look at C only above 0.7.
classnetwork <- graph_from_data_frame(d = edges, vertices = nodes, directed = FALSE)
c_table = tibble(id = 1:140, C=transitivity(classnetwork, 'local'))
hist(c_table$C, main='Distribution of C values', xlab = 'C value')
In this plot we can see that the majority is C < 0.5
Question 23: Find out the weak ties in this network. Highlight three of them, and show to your teachers.
[Question 23, 2 points].
nodes <- tibble(id = 1:140, C=transitivity(classnetwork, 'local'), color=ifelse(C < 0.3, "red", "lightblue"))
visNetwork(nodes, edges) %>%
visIgraphLayout(layout = "layout_nicely", smooth = TRUE) %>%
visNodes(shape="circle")tibble(id = 1:140, C=transitivity(classnetwork, 'local'), D=degree(classnetwork)) %>%
arrange(C, D) %>% head()## # A tibble: 6 × 3
## id C D
## <int> <dbl> <dbl>
## 1 72 0 2
## 2 60 0 4
## 3 105 0.143 7
## 4 130 0.143 7
## 5 96 0.167 4
## 6 38 0.179 8
Question 24: Following the above question, can you find out the nodes with many weak ties? One possible way is to find the nodes with high betweenness, using the below code. Find out 3 nodes with highest betweenness. Check their positions in the network again to see if they have many weak ties as predicted. Show your teacher the IDs of these 3 nodes, and explain why betweenness can be used as a proxy to find out nodes with many weak ties? [Question 24, 5 points].
nodes <- tibble(id = 1:140,
B=betweenness(class2022, v=V(class2022), directed = FALSE, weights =NULL),
C=transitivity(classnetwork, 'local'),
color=ifelse(B > 500, "red", "lightblue"))
nodes$color[nodes$C > 0.6] = "orange"
nodes['label'] = ''
nodes$label[nodes$color == 'red'] = 'B'
nodes$label[nodes$color == 'orange'] = 'C'
visNetwork(nodes, edges) %>%
visIgraphLayout(layout = "layout_nicely", smooth = TRUE) %>%
visNodes(shape="circle", value=nodes)nodes %>% arrange(-B) %>% head()## # A tibble: 6 × 5
## id B C color label
## <int> <dbl> <dbl> <chr> <chr>
## 1 20 550. 0.275 red "B"
## 2 71 546. 0.233 red "B"
## 3 80 541. 0.233 red "B"
## 4 106 494. 0.181 lightblue ""
## 5 2 483. 0.209 lightblue ""
## 6 8 419. 0.25 lightblue ""
hist(nodes$B, main='Distribution of Betweeness values', xlab = 'betweenness value')
Knowingly that high betweenness is about the occurrence of a node between others (has a lot of influence, e.g. control over the flow) and weak ties is about e.g. a friend which does not belong to close friends nor family but still able to pass on information then betweenness can be useful for finding those nods which are middle nods controlling the flow, these middle nods can e.g. pass on information to other group of nods.
Question 25: Next, we will test “the strength of weak ties in simple contagion”. First, build a simple independent cascade (IC) model with the following characteristics:
After you build the model, apply it to the network of this class, and develop a plot with x-axis as Day, y-axis as the Cumulative Infected Percentage by Day. (You can use the code below to build the IC model or develop your own codes.) [Question 25, 10 points]
stopifnot(require(data.table))## Loading required package: data.table
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
stopifnot(require(Matrix))## Loading required package: Matrix
#search the neighbours of contangious node, and their neighbours have a chance of Pprob being infected
calculate_value2 <- function(node, each_neighbors, Pprob){
return(each_neighbors[[node]][ which(runif(length(each_neighbors[[node]]), 0, 1) <= Pprob)])
}
### function to model simple contagion, return can be modified if you want to see other outputs#
IC <- function(i, seeds, nNode, each_neighbors, Pprob){
node_seed=seeds[i]
node_value <- rep.int(0, nNode) ### node_value
infection <- rep.int(0, nNode) ### infection status 0 uninfected 1 infected
new_infected <- list() # to store the ID of nodes being activated each time step
day <- 1
day_infected <- 1 ### infected num daily
new_infected[[day]] <- node_seed
node_value[as.numeric(node_seed)] <- 1 ### assign value for seed nodes
infection[as.numeric(node_seed)] <- 1
day <- 2
old_infected <- which(infection == 1)
ContagiousID <- which(infection == 1) #nodes that are contagious
activeID <- unlist(lapply(ContagiousID,calculate_value2,each_neighbors,Pprob))
newactiveID<- setdiff(activeID, old_infected)
if (length(activeID)!=0){
new_infected[[day]] <- newactiveID
infection[as.numeric(newactiveID)] <- 1
}else {
new_infected[[day]]<-integer(0)
}
day_infected[day] <- sum(infection == 1)
day<-3
old_infected <- which(infection == 1)
ContagiousID <- which(infection == 1) #nodes that are contagious
activeID <- unlist(lapply(ContagiousID,calculate_value2,each_neighbors,Pprob))
newactiveID <- setdiff(activeID, old_infected)
if (length(activeID)!=0){
new_infected[[day]]<-newactiveID
infection[as.numeric(newactiveID)] <- 1
}else {
new_infected[[day]]<-integer(0)
}
day_infected[day] <- sum(infection == 1)
for (day in c(4:30)){
old_infected<-which(infection == 1)
Within3days <- list(new_infected[[day-3]],new_infected[[day-2]],new_infected[[day-1]])
ContagiousID <- unlist(Within3days, use.names = FALSE)
activeID <- unlist(lapply(ContagiousID, calculate_value2, each_neighbors, Pprob))
newactiveID <- setdiff(activeID, old_infected)
if (length(activeID)!=0){
new_infected[[day]] <- newactiveID
infection[as.numeric(newactiveID)] <- 1
}else {
new_infected[[day]] <- integer(0)
}
day_infected[day] <- sum(infection == 1)
day=day+1
}
return(day_infected)
}
##### function of IC, to be used multiply according to different seeds (i.e., who in the network get the virus first)
IC2 <- function(i, seeds, nNode, each_neighbors, Pprob) {
day_infected_matrix <- do.call('rbind', lapply(1:i, IC, seeds, nNode, each_neighbors, Pprob))
return(day_infected_matrix)
}# For the original network :
# prepare the seed matrix:
#introduce a patient to an otherwise healthy population, this person should be randomly selected to mimic the contagion in real-world
#since we have only 140 nodes, we can select them one by one; therefore, we will have 140 different staring points to kick-off the congation process
seeds<-c(1:140)
nNode=140 #size of the original network
Pprob=0.1 #the probability for an infected person to pass down the virus to his/her neighbours
adj_matrix <- igraph::as_adjacency_matrix(class2022, type = 'both')
each_neighbors <- which(adj_matrix > 0, arr.ind = TRUE)
each_neighbors <- split(each_neighbors[, 2], each_neighbors[, 1]) #get the neigbhour list of each node
result1<- replicate(100, IC2(140, seeds, nNode, each_neighbors, Pprob), simplify=FALSE) #run at least 100 times since each IC model run has its own randomness (see the "calculate_value2 function")
result1 <- do.call(rbind, result1)
infections1 = colMeans(result1)/(140) # percentage of infected person each day
###
all_days = 1:length(infections1)
CIP_result_table = tibble(day=all_days, infection1=infections1)CIP_result_table %>%
ggplot(aes(x=day, y=infection1, color='fullNods')) +
geom_line() +
geom_point() +
labs(title = "Distribution of CIP", x= "Day", y='Cumulative Infected Percentage')
Please see the answer of Question 26, answers of questions 25 and 26 are together in there.
Question 26: Revised the original network of “class2022” by deleting the 3 vertices with highest betweenness (edges directly connecting to these 3 vertices are deleted). Apply the same IC model to this new network, how does the daily cumulative infected percentage change? What if you delete the 3 vertices with lowest betweenness? Show the results to your teacher by plots, and explain the results by “the strength of weak ties”. [Question 26, 15 points]
simulate_spread <- function(to_remove, class2022){
class2022_2 = delete_vertices(class2022, to_remove)
nNode=length(V(class2022_2)) #size of the original network
seeds<-c(1:nNode)
Pprob=0.1 #the probability for an infected person to pass down the virus to his/her neighbours
adj_matrix <- igraph::as_adjacency_matrix(class2022_2, type = 'both')
each_neighbors <- which(adj_matrix > 0, arr.ind = TRUE)
each_neighbors <- split(each_neighbors[, 2], each_neighbors[, 1]) #get the neigbhour list of each node
result1 <- replicate(100, IC2(nNode, seeds, nNode, each_neighbors, Pprob), simplify=FALSE) #run at least 100 times since each IC model run has its own randomness (see the "calculate_value2 function")
result1 <- do.call(rbind, result1)
infections = colMeans(result1)/(nNode) # percentage of infected person each day
infections[is.na(infections)] = 0
return(infections)
}class2022 = readRDS("classnetwork_2022.rds")
# remove nodes
nodes <- tibble(id = 1:140, B=betweenness(class2022, v=V(class2022), directed = FALSE, weights =NULL))
to_remove = nodes %>% arrange(-B) %>% slice_head(n=3) %>% select(id)
infections2 = simulate_spread(to_remove$id, class2022)CIP_result_table['infection2'] = infections2to_remove = nodes %>% arrange(B) %>% slice_head(n=3) %>% select(id)
infections3 = simulate_spread(to_remove$id, class2022)CIP_result_table['infection3'] = infections3colors <- c("fullNods" = "red", "removedHighest" = "green", "removedLowest" = "blue")
CIP_result_table %>%
ggplot(aes(x=day)) +
geom_line(aes(y=infection1, color='fullNods')) +
geom_point(aes(y=infection1, color='fullNods')) +
geom_line(aes(y=infection2, color='removedHighest')) +
geom_point(aes(y=infection2, color='removedHighest')) +
geom_line(aes(y=infection3, color='removedLowest')) +
geom_point(aes(y=infection3, color='removedLowest')) +
labs(title = "Distribution of CIP", x= "Day", y='Cumulative Infected Percentage') +
scale_color_manual(values=colors)
From the fullNods (RED) perspective (Question 25) we see an increase of CIP over days meaning a increase of strength of weak tie spreading the disease faster.
After removing the highest 3 betweenness (GREEN) the strength of weak tie decreases, meaning it has lower diffusion of contagion compared to fullNods (RED). Removing the lowest 3 betweenness slightly increases the strength of weak tie (BLUE) leading to slightly higher CIP.
Question 27: Inspired by your observations from Question 26, can you relate it to some of the COVID measures implemented in the real world? [Question 27, 3 points]
So we had severer measures for weak ties during corona infection. For example - we could have a limited number of visitors during the restrictions period of Covid-19 (strong ties e.g. family, friends), but we had to keep one and half meter distance in public and wear a mask when interacting with weak ties. Moreover, one had to work from home.
Question 28: Recall the fourth and fifth questions in the survey as shown below. They are designed to collect the parameters that we will need to build 1) an independent cascade (IC) model and 2) a threshold model. Can you see which question is for the IC model and which one is for the threshold model? [Question 28, 2 point].
Q 4 is related to the IC model and Q 5 is related to the threshold model. Q 5 is asking about how many people you need to be convinced of adopting a novel behavior, and it is obviously related to threshold model. Q 4 is about trying something new and want to share it. It is related to IC model.
Question 29: We asked your response for three types of behaviors (i.e., Share/watch Youtube video, Share/try a vegetarian recipe, Share/read a paper related to the lecture). The distribution of your answers are shown in the above figure. Among these three types of behaviors, which one is the least contagious? Which one is the most contagious? And why? [Question 29, 3 point].
Least contagious
Vegetarian recipe, just because on average one will probably share it with fewer people, hence the answer 1 dominates in that response (very unlikely to share)
Most contagious
A paper related to the lecture. Firstly, just by looking at the responses - one is very likely to share that (high percentage of fours and fives in that response). Secondly, the likelihood of contagion probably increases because of the context where the survey was conducted. Since, it is in a university, just makes sense that sharing a paper is more relevant than sharing something random with classmates.
Question 30: Based on the model that you’ve built, if you can seed only one person, who will you choose? Explain to your teacher 1) how do you find out this node, 2) the ID of this node 3) its attributes (e.g, degree, betweenness, probability to share), and compared such attributes to other nodes in the network to explain why it is a good candidate to start the diffusion [Question 30, 10 point].
stopifnot(require(data.table))
stopifnot(require(Matrix))
#So added infection to the calculate_value function
calculate_value4 <- function(nodeId, infection, each_neighbors, Pprob){
# get neighbors of the nod
k = each_neighbors[[nodeId]]
# check who watched and not watched
already_infected = which(infection[as.numeric(k)] == 1)
already_not_infected = which(infection[as.numeric(k)] == 0)
# if NONE has already watched then apply to all K otherwise only apply to the not watched neighbors
ran_length = ifelse(length(already_infected) == 0, length(k), length(already_not_infected))
ran_nums = runif(ran_length, 0, 1)
infected = list()
not_infected = list()
# if NONE has already watched then apply to all K
if(length(already_infected) == 0){
# get all watched according to condition
infected = k[which(ran_nums <= Pprob[nodeId])]
# get all none watched according to condition
not_infected = k[which(ran_nums > Pprob[nodeId])]
}else{
# only apply to the not watched neighbors
# get all watched according to condition then concatenate with the already watched ones
infected = unlist(list(already_infected, already_not_infected[which(ran_nums <= Pprob[nodeId])]))
# get all none watched according to condition for later analysis
not_infected = already_not_infected[which(ran_nums > Pprob[nodeId])]
}
# not_infected contains nods with status 0
# if there are more then 0 then proceed the process of checking their neighbors if they are activated
if(length(not_infected) > 0){
# loop over each nod with status 0
for(targetId in not_infected){
# get the neighbors
target_k = each_neighbors[[targetId]]
# get count of the neighbors who watched
k_states = sum(infection[as.numeric(target_k)] == 1)
# when more than 1 neighbor has seen the video
if(k_states > 1){
# generate sequential increasing probabilities (depending on watched neighbors) examples: 0.1, 0.2, 0.3 when
# 3 neighbors has seen the video
probas = seq(0.1, (length(k_states)/10), 0.1)
# generate for each a random number
random_number = runif(length(k_states), 0, 1)
# validate per neighbor is the nod with state 0 should be activated or not
for(idx in 1:length(probas)){
ran_num = random_number[idx]
prob = probas[idx]
# when condition is met then save the nods ID as activated
if(ran_num < prob){
infected = c(infected, nodeId)
}
}
}
}
}
# return all activated nod IDs
return(infected)
}
### FUNCTION TO MODEL SIMPLE CONTAGION, return can be modified if you want to see other outputs
IC_watched <- function(seeds, nNode, each_neighbors, Pprob, max_days){
node_seed = seeds
node_value <- rep.int(0, nNode) ### node_value
infection <- rep.int(0, nNode) ### infection status 0 uninfected 1 infected
new_infected <- list() # to store the ID of nodes being activated each time step
day <- 1
day_infected <- 1 ### infected num daily
new_infected[[day]] <- node_seed
node_value[as.numeric(node_seed)] <- 1 ### assign value for seed nodes
infection[as.numeric(node_seed)] <- 1
day_infected[day] <- sum(infection == 1)
for (day in c(2:max_days)){
old_infected <- which(infection == 1)
prev_day_nodes <- unlist(new_infected[[day-1]], use.names = FALSE)
activeID = list()
for(nodeId in prev_day_nodes){
tmp = calculate_value4(nodeId, infection, each_neighbors, Pprob)
activeID = c(activeID, tmp)
}
activeID = unlist(activeID)
newactiveID <- setdiff(activeID, old_infected)
if (length(activeID) != 0){
new_infected[[day]] <- activeID
infection[as.numeric(newactiveID)] <- 1
}else {
new_infected[[day]] <- list()
}
day_infected[day] <- sum(infection == 1, na.rm = T)
}
return(list(day_infected, infection))
}
##### FUNCTION OF IC, to be used multiply according to different seeds (i.e., who in the network get the virus first)
IC3 <- function(seeds, nNode, each_neighbors, Pprob, max_days) {
day_infected_matrix = IC_watched(seeds, nNode, each_neighbors, Pprob, max_days)
return(day_infected_matrix)
}#Load the data of who have watched the YouTube video with probabilities
class2022 <- readRDS("classnetwork_2022.rds")
class2022_attributes <- readRDS("class2022_attributes.rds")
class2022_matrix <- as_edgelist(class2022, names = TRUE)
edges <- data.frame(from = class2022_matrix[,1], to = class2022_matrix[,2])# Seed nod with the highest betweenness
seed_nod = tibble(id = 1:140, B=betweenness(class2022), D=degree(class2022), C=transitivity(class2022, 'local'), TH=class2022_attributes$pYoutube) %>% arrange(-B)
seed_nod## # A tibble: 140 × 5
## id B D C TH
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 20 550. 18 0.275 0.225
## 2 71 546. 16 0.233 0.225
## 3 80 541. 16 0.233 0.225
## 4 106 494. 15 0.181 0.325
## 5 2 483. 14 0.209 0.225
## 6 8 419. 8 0.25 0.325
## 7 104 398. 10 0.267 0.325
## 8 66 382. 10 0.289 0.325
## 9 122 363. 14 0.352 0.425
## 10 65 344. 13 0.231 0.425
## # … with 130 more rows
Q: If you can seed only one person to spread the recipe, who will you choose?
A: We decided to seed the nod with the highest betweenness. Obviously, there are other approaches, but we chose betweenness.
Q: How do you find out this node?
A: See above
Q: The ID of the node?
A: See above, ID is 20.
Q: Its attributes (e.g, degree, betweenness, probability to share)
A: See above (degree is 18, betweenness is 549.7434 and the probability to share is 0.225).
Even though the probability is low, but having high betweeness is more important in regards to the question at hand (you don’t have time and energy), you do not want to infect everyone, but want to infect “everyone” as fast as you can. So betweenness is a more important metric in our opinion and is the basis of why we chose it.
nNode = length(class2022_attributes$ID) #size of the original network
Pprob = class2022_attributes$pYoutube #the probability for an infected person to pass down the virus to his/her neighbors
adj_matrix <- igraph::as_adjacency_matrix(class2022, type = 'both')
each_neighbors <- which(adj_matrix > 0, arr.ind = TRUE)
each_neighbors <- split(each_neighbors[, 2], each_neighbors[, 1]) #get the neighbor list of each
seeds = c(seed_nod$id) %>% head(1)
computed <- IC3(seeds, nNode, each_neighbors, Pprob, max_days = 8)result1 = computed[1]
infected_status = computed[[2]]
result1 <- do.call(rbind, result1)
infections1 = colMeans(result1, na.rm = T) / 140 # percentage of infected person each day
all_days = 1:length(infections1)
CIP_result_youtube_table = tibble(day=all_days, watched=infections1)CIP_result_youtube_table %>%
ggplot(aes(x=day, y=watched, color='fullNods')) +
geom_line() +
geom_point() +
labs(title = "Distribution of CIP, YouTube", x= "Day", y='Cumulative Infected Percentage')
Moreover, here in the graph one can see, that on the basis on our assumption the video will be sufficiently spread in five days. Also, here below is a visualization of how many nods will be infected. Considering the time span one could argue that is pretty aggressive “infection” spread.
Visualizing the nods who watched the YouTube video
nodes <- tibble(
id = 1:length(infected_status),
watched=infected_status,
color=ifelse(watched == 1, "red", "lightblue"))
visNetwork(nodes, edges) %>%
visIgraphLayout(layout = "layout_nicely", smooth = TRUE) %>%
visNodes(shape="circle")If you can seed only one person to spread the recipe, who will you choose? Explain to your teacher 1) how do you find out this node, 2) the ID of this node 3) its attributes (e.g, degree, betweenness, threshold), and compared such attributes to other nodes in the network to explain why it is a good candidate to start the diffusion [Question 31, 10 point].
stopifnot(require(data.table))
stopifnot(require(Matrix))
#Modified the functions
calculate_value5 <- function(nodeId, infection, each_neighbors, Pprob){
if(infection[as.numeric(nodeId)] == 1){
return(nodeId)
}else{
k = each_neighbors[[nodeId]]
infected_neighbours = which(infection[as.numeric(k)] == 1)
proba_neighbours = Pprob[as.numeric(infected_neighbours)]
total_convincing = sum(proba_neighbours, na.rm = T)
infected = 0
if(total_convincing > Pprob[nodeId]){
infected = nodeId
}
return(infected)
}
}
###Threshold model for the vegetarian recipe
TH_vega <- function(seeds, nNode, each_neighbors, Pprob, max_days){
node_seed = seeds
infection <- rep.int(0, nNode) ### infection status 0 uninfected 1 infected
new_infected <- list() # to store the ID of nodes being activated each time step
day <- 1
day_infected <- 1 ### infected num daily
new_infected[[day]] <- node_seed
infection[as.numeric(node_seed)] <- 1
day_infected[day] <- sum(infection == 1)
for (day in c(2:max_days)){
old_infected <- which(infection == 1)
activeID = list()
for(nodeId in 1:nNode){
tmp = calculate_value5(nodeId, infection, each_neighbors, Pprob)
if(tmp != 0){
activeID = c(activeID, tmp)
}
}
newactiveID <- setdiff(activeID, old_infected)
if (length(activeID) != 0){
new_infected[[day]] <- activeID
infection[as.numeric(newactiveID)] <- 1
}else {
new_infected[[day]] <- activeID
}
day_infected[day] <- sum(infection == 1, na.rm = T)
}
return(list(day_infected, infection))
}
##### Threshold function, to be used multiply according to different seeds (i.e., who in the network get the virus first)
TH4 <- function(seeds, nNode, each_neighbors, Pprob, max_days) {
day_infected_matrix = TH_vega(seeds, nNode, each_neighbors, Pprob, max_days)
return(day_infected_matrix)
}#Again load data
class2022 <- readRDS("classnetwork_2022.rds")
class2022_attributes <- readRDS("class2022_attributes.rds")
class2022_matrix <- as_edgelist(class2022, names = TRUE)
edges <- data.frame(from = class2022_matrix[,1], to = class2022_matrix[,2])# Seed nod 109
seed_nod = tibble(id = 1:140, B=betweenness(class2022), D=degree(class2022), C=transitivity(class2022, 'local'), Recipy_TH=class2022_attributes$tRecipe) %>% arrange(-C, -D)
seed_nod## # A tibble: 140 × 5
## id B D C Recipy_TH
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 88 0 5 1 0.92
## 2 47 0 4 1 0.67
## 3 35 0 3 1 0.0125
## 4 43 0 3 1 0.67
## 5 78 0 3 1 0.92
## 6 52 0 2 1 0.0125
## 7 63 0 2 1 0.33
## 8 109 1.61 7 0.905 0.0125
## 9 102 4.72 7 0.857 0.0125
## 10 11 2.52 4 0.833 0.0925
## # … with 130 more rows
If you can seed only one person to spread the recipe, who will you choose? Explain to your teacher 1) how do you find out this node, 2) the ID of this node 3) its attributes (e.g, degree, betweenness, threshold), and compared such attributes to other nodes in the network to explain why it is a good candidate to start the diffusion.
Compared to the IC model beforehand, this model is adjusted to capture complex contagion.
- See above
- The ID of the node is 109
- It’s betweenness is 1.6, degree is 7 and threshold is 0.125
It’s a good candidate since it has an initial low threshold and high degree and clustering coefficient. Again, considering the question at hand, here is an example. So let’s call nod 109 Karen. Okay, Karen has an initial low threshold but a high degree and clustering coefficient, meaning, that she has a lot of connections. Even though she might not convince someone to try the vegetarian recipe on her own, many triads could and will form and eventually after a couple of time steps the total threshold will be so high, that it could convince an individual with a really high threshold.
Please see the graph below.
nNode = length(class2022_attributes$ID) #size of the original network
Pprob = class2022_attributes$tRecipe #the probability for an infected person to pass down the virus to his/her neighbors
adj_matrix <- igraph::as_adjacency_matrix(class2022, type = 'both')
each_neighbors <- which(adj_matrix > 0, arr.ind = TRUE)
each_neighbors <- split(each_neighbors[, 2], each_neighbors[, 1]) #get the neighbor list of each
seeds = c(seed_nod[8,]$id)
computed <- TH4(seeds, nNode, each_neighbors, Pprob, max_days = 12)result1 = computed[1]
infected_status = computed[[2]]
result1 <- do.call(rbind, result1)
infections1 = colMeans(result1, na.rm = T) / 140 # percentage of infected person each day
all_days = 1:length(infections1)
CIP_result_youtube_table = tibble(day=all_days, watched=infections1)CIP_result_youtube_table %>%
ggplot(aes(x=day, y=watched)) +
geom_line() +
geom_point() +
labs(title = "Distribution of CIP, VEGA", x= "Day", y='Cumulative Infected Percentage')
The adopters of the recipe follow an empirical distribution of a
successful diffusion study. It follows a stretched out S shape which
corresponds to: 1) Up until approximately day 2.5 there are there are
the “innovators” (initial diffusion) 2) Followed by “opinion leaders”
(early adapters) till about day 5 (θ ≤ 0.16), and the diffusion starts
to take off 3) From day 5 to day 7 we have the early majority (0.16 ≤
θ ≤ 0.5) 4) From approximately day 7 to day 9 we have the late
majority (0.5 ≤ θ ≤ 0.84) 5) From approximately day 9 to day 12 we
have the “laggards”.
At day 12 everyone should have heard about the vegeterian recipe.
nodes <- tibble(
id = 1:length(infected_status),
watched=infected_status,
color=ifelse(watched == 1, "red", "lightblue"))
visNetwork(nodes, edges) %>%
visIgraphLayout(layout = "layout_nicely", smooth = TRUE) %>%
visNodes(shape="circle")One can see that the network is fully infected after 12 days.
Question 32
We didn’t need to remove the 140 nodes one by one, hence our function.
If you want to find out 5 people to seed, can you apply such explicit search again?
Yes, again based on the degree and clustering coefficient from the tibble.
Question 33
Now you have a little more time or budget to seed more than 1 persons in this class. To promote the Youtube video, you can seed 3 people. To promote the vegetarian recipe, you can seed 5 people. Use degree heuristics, betweenness heuristics and greedy algorithm to find out 1) 3 seeds to spread Youtube video and 2) 5 seeds to spread the vegetarian recipe. Show your teachers 1) the ID of these seeds, 2) Among degree heuristics, betweenness heuristics and greedy algorithm, which method provides the best result? 3) Do you find some common properties of the seeds provided by different methods? [Question 33, 20 point].
lala = tibble(
ID=1:140,
B=betweenness(class2022)
)
comm = cluster_walktrap(class2022)
lolz = function(x, comm){
for(clusterIdx in 1:length(comm)){
if(x %in% unlist(comm[clusterIdx])){
return(clusterIdx)
}
}
return(NA)
}
lala = lala %>% mutate(commID = sapply(FUN = lolz, ID, comm), D=degree(class2022))communities_grouped = lala %>% arrange(-commID, -D) %>% group_by(commID) %>% summarise(maxD = max(D), B=max(B)) %>% arrange(B, -maxD)
communities_grouped## # A tibble: 11 × 3
## commID maxD B
## <int> <dbl> <dbl>
## 1 9 11 218.
## 2 1 9 240.
## 3 6 14 277.
## 4 8 9 330.
## 5 3 11 334.
## 6 7 12 382.
## 7 10 8 419.
## 8 2 15 494.
## 9 5 16 541.
## 10 4 16 546.
## 11 11 18 550.
justin = function(bieber_ID, each_neighbors, lala2, communities_grouped){
k = unlist(each_neighbors[bieber_ID])
highest_d_communities = 1:length(k)
for(idx in 1:length(k)){
k_id = k[idx]
k_comm = lala2[lala2$ID == k_id, ]$commID
highest_d_communities[idx] = communities_grouped[communities_grouped$commID == k_comm, ]$maxD
}
return(sum(highest_d_communities))
}
adj_matrix <- igraph::as_adjacency_matrix(class2022, type = 'both')
each_neighbors <- which(adj_matrix > 0, arr.ind = TRUE)
each_neighbors <- split(each_neighbors[, 2], each_neighbors[, 1])
lala2 = lala %>% arrange(-B)
candidates = lala2 %>% mutate(total_comm_D=sapply(FUN = justin, ID, each_neighbors, lala2, communities_grouped))candidates %>% arrange(-total_comm_D, -D)## # A tibble: 140 × 5
## ID B commID D total_comm_D
## <int> <dbl> <int> <dbl> <dbl>
## 1 20 550. 11 18 254
## 2 80 541. 5 16 243
## 3 71 546. 4 16 228
## 4 106 494. 2 15 219
## 5 122 363. 5 14 208
## 6 2 483. 4 14 205
## 7 65 344. 4 13 203
## 8 100 335. 4 13 203
## 9 50 279. 11 12 197
## 10 138 266. 6 14 196
## # … with 130 more rows
youtube_seeds = candidates %>% arrange(-total_comm_D) %>% head(3)
vega_seeds = candidates %>% arrange(-total_comm_D) %>% head(5)vega_model <- TH4(vega_seeds$ID, nNode, each_neighbors, Pprob, max_days = 12)result1 = vega_model[1]
infected_status = vega_model[[2]]
result1 <- do.call(rbind, result1)
infections1 = colMeans(result1, na.rm = T) / 140 # percentage of infected person each day
all_days = 1:length(infections1)
tibble(day=all_days, watched=infections1) %>%
ggplot(aes(x=day, y=watched)) +
geom_line() +
geom_point() +
labs(title = "Distribution of CIP, VEGA", x= "Day", y='Cumulative Infected Percentage')
youtube_model <- IC3(youtube_seeds$ID, nNode, each_neighbors, Pprob, max_days = 6)result1 = youtube_model[1]
infected_status = youtube_model[[2]]
result1 <- do.call(rbind, result1)
infections1 = colMeans(result1, na.rm = T) / 140 # percentage of infected person each day
all_days = 1:length(infections1)
tibble(day=all_days, watched=infections1) %>%
ggplot(aes(x=day, y=watched)) +
geom_line() +
geom_point() +
labs(title = "Distribution of CIP, YouTube", x= "Day", y='Cumulative Infected Percentage')
ABOVE we combine Degree and Betweenness:
See the plots and the tibble above for the results for choosing the set of seeds (The Seed IDs) for YouTube and Vegetarian recipe. We have chosen to combine the heuristics Degree and Betweenness to find a set of seeds who have high influence for diffusion. The idea is to sum the degrees of all neighbouring communities of node i to give it a score as the strength of influence. Further, the scores of each node is arranged from high to low to pick the best 3 or 5.
class2022_attributes <- readRDS("class2022_attributes.rds")
youtube_seeds_D = candidates %>% arrange(-D) %>% mutate(Pprob=class2022_attributes$pYoutube) %>% head(3)
youtube_seeds_B = candidates %>% arrange(-B) %>% mutate(Pprob=class2022_attributes$pYoutube) %>% head(3)
vega_seeds_D = candidates %>% arrange(-D) %>% mutate(Pprob=class2022_attributes$tRecipe) %>% head(5)
vega_seeds_B = candidates %>% arrange(-B) %>% mutate(Pprob=class2022_attributes$tRecipe) %>% head(5)Here we compare Degree with Betweenness
compute_stuff = function(FUN, seeds, nNode, each_neighbors, Pprob, max_days, attribute, type_context){
model <- FUN(seeds, nNode, each_neighbors, Pprob, max_days = max_days)
result1 = model[1]
infected_status = model[[2]]
result1 <- do.call(rbind, result1)
infections1 = colMeans(result1, na.rm = T) / 140
all_days = 1:length(infections1)
return(tibble(day=all_days, infected=infections1, attribute=attribute, type_context=type_context))
}
all_combi_seeds = list(youtube_seeds_D, vega_seeds_D, youtube_seeds_B, vega_seeds_B)
all_combi_proba = list(class2022_attributes$pYoutube, class2022_attributes$tRecipe, class2022_attributes$pYoutube, class2022_attributes$tRecipe)
all_combi_attr = list('D', 'D', 'B', 'B')
all_combi_type = list('Y', 'V', 'Y', 'V')
all_func = list(IC3, TH4)
result_tibble = list(tibble(day=numeric(), infected=numeric(), attribute=character(), type_context=character()),
tibble(day=numeric(), infected=numeric(), attribute=character(), type_context=character()))
days = c(6, 12)
for(idx in 1:length(all_func)){
fun = all_func[[idx]]
fun_tibble = result_tibble[[idx]]
total_day = days[[idx]]
for(idy in 1:length(all_combi_seeds)){
combi = all_combi_seeds[[idy]]
attribute = all_combi_attr[[idy]]
type_context = all_combi_type[[idy]]
proba = all_combi_proba[[idy]]
result = compute_stuff(FUN=fun, combi$ID, nNode, each_neighbors, proba, max_days = total_day, attribute=attribute, type_context=type_context)
result_tibble[[idx]] = bind_rows(result_tibble[[idx]], result)
}
}result_tibble[[1]] %>%
ggplot(aes(x=day, y=infected)) +
geom_line() +
labs(title = "Distribution of CIP, IC Model", x= "Day", y='Cumulative Infected Percentage') +
facet_wrap(type_context ~ attribute)
Using the IC model we can see from the result that for Vegatarian (V) and YouTube (T) the attributes Betweenness (B) and Degree (D) provide not similar result, they differ in the amount of CIP over time steps.
result_tibble[[2]] %>%
ggplot(aes(x=day, y=infected)) +
geom_line() +
labs(title = "Distribution of CIP, Threshold Model", x= "Day", y='Cumulative Infected Percentage') +
facet_wrap(type_context ~ attribute)
- Among degree heuristics, betweenness heuristics and greedy algorithm, which method provides the best result?
Using the Threshold model we can see from the result that for Vegatarian (V) and YouTube (Y) the attributes betweenness and degree have similar result. Whereas using the IC model we can see from the result that for Vegatarian (V) and YouTube (T) the attributes Betweenness (B) and Degree (D) provide not similar result, they differ in the amount of CIP over time steps.
- Do you find some common properties of the seeds provided by different methods?
Yes there are common properties specially looking at the result of Threshold model where for YouTube and Vega the Betweenness and Degree attributes are similar.
# fun_tibble = tibble(day=numeric(), infected=numeric(), attribute=character(), type_context=character(), reference=numeric())
# total_day = 12
#
# num_seeds = 1:140
#
# vega_probas = class2022_attributes$tRecipe
#
# for(idy in 1:140){
#
# curr_seeds = num_seeds[1:(length(num_seeds)+1)-idy]
# curr_seeds = sample(curr_seeds)
#
# combi = curr_seeds
#
# attribute = 'not important'
# type_context = 'Vega'
# result = compute_stuff(FUN=TH4, combi$ID, nNode, each_neighbors, vega_probas, max_days = total_day, attribute=attribute, type_context=type_context)
#
# sorty = result %>% arrange(-infected) %>% head(1)
# sorty = sorty %>% mutate((reference=length(num_seeds)+1)-idy)
#
# fun_tibble = bind_rows(fun_tibble, sorty)
# }peep = function(x, t){
input = (x - t) / 1
to_return = 1 / (1+exp(-input))
return(to_return)
}
result = sapply(1:14, peep, 7)
plot(result, type='b', ylab='Probability', xlab = 'total convincing neigh')