Introduction

In Homework 1 we studied the phone-calls layer of the Copenhagen Networks Study (Sapiezynski et al., 2019). We argued that phone calls are a strong-tie signal because you do not call just anyone. The calls graph had 504 nodes and 2,654 edges, which is too sparse for the link-prediction exercise in Homework 2 (the spec asks for more than 300 nodes and more than 5,000 edges).

For Homework 2 we therefore move to a different layer of the same study: the Bluetooth proximity scans. Bluetooth captures physical co-presence (two phones detecting each other within radio range), and it is widely used in computational social science as a proxy for face-to-face contact (Sekara & Lehmann, 2014). The dataset records about 5.5 million scan events between 692 students over four weeks.

After the aggregation rules described in Section 1 we obtain an undirected network with 662 nodes and 9,251 edges. The network is sparse but link prediction is still not trivial on it, and it forms one connected component.

Our pipeline follows the link-prediction protocol from Lecture 4 and Workshop 3: delete a fraction of real edges, sample a matched non-edge set as the negative class, score every candidate pair with structural similarity heuristics, and ask a binary classifier to recover the deleted edges. Three things go beyond Workshop 3: (i) seven heuristics instead of three, drawing on Lecture 4 slides 19-24 and 32-33; (ii) 10-fold cross-validation, which the homework brief asks for; (iii) GLM coefficients compared with Random Forest permutation importance to identify which heuristic does the most predictive work.

library(tidyverse)
library(igraph)
library(ggraph)
library(scales)
library(caret)
library(pROC)
library(randomForest)

1. The network

1.1 Data acquisition and aggregation

The raw file bt_symmetric.csv (98 MB, 5,474,289 rows) records every Bluetooth scan event. Each row has columns timestamp, user_a, user_b, rssi, with two sentinel encodings documented by the authors: user_b == -1 is an empty scan (no device detected) and user_b == -2 is a non-experiment device (devices outside the cohort are all assigned the same id). After dropping the sentinels we have 2,426,279 real participant-pair scans.

To turn the scan stream into an undirected network we apply two filters, both standard in the CNS-derived literature:

  1. RSSI ≥ −80 dBm. RSSI (received signal strength indicator) is roughly related to physical distance, but Sekara & Lehmann
    1. note that signal strength is a noisy proxy. They use a threshold which “in a large majority of cases corresponds to interactions that occur within a radius of approximately 2 meters”. We use −80 dBm as a slightly more permissive room-scale cut, which trades a bit more noise for keeping more of the indoor co-presence (about 30% of the real-pair scans pass this cut, vs. the median scan strength of −86 dBm).
  2. Minimum 10 mutual scans per dyad. Each scan is one of the five-minute discovery windows in the CNS protocol. Ten scans correspond to roughly 50 minutes of cumulative co-presence over the four-week study, which we take as the lower bound for an intentional social tie rather than an incidental encounter (passing in a hallway).

Pre-built artifacts from this aggregation live in bt_network.rds (igraph object) and bt_edges.csv (edgelist with weight = count of mutual scans plus mean RSSI and time of first/last contact).

g <- readRDS("bt_network.rds")
g
IGRAPH 263f699 UNW- 662 9251 -- 
+ attr: name (v/c), weight (e/n), mean_rssi (e/n), t_first (e/n),
| t_last (e/n)
+ edges from 263f699 (vertex names):
 [1] 20 --21  48 --49  90 --91  57 --101 144--181 12 --244 176--263 221--263
 [9] 119--283 104--288 318--375 62 --401 414--421 12 --454 260--459 263--472
[17] 461--492 311--524 343--553 16 --572 146--611 279--617 96 --634 207--677
[25] 19 --47  244--454 249--530 568--599 519--663 188--673 182--183 170--390
[33] 455--621 222--471 445--624 235--370 83 --105 263--370 455--282 215--664
[41] 235--472 221--472 267--452 22 --220 297--374 95 --190 378--505 215--673
[49] 282--621 285--523 208--422 64 --525 419--557 190--665 90 --665 263--69 
+ ... omitted several edges

1.2 Network properties

deg   <- degree(g)
trans <- transitivity(g, type = "global")
comp  <- components(g)

tibble(
  Metric = c("Nodes (|V|)", "Edges (|E|)", "Density",
             "Mean degree", "Median degree", "Max degree",
             "Nodes with degree > 10",
             "Connected components", "Largest component",
             "Global transitivity (clustering)"),
  Value  = c(comma(vcount(g)), comma(ecount(g)),
             sprintf("%.4f", edge_density(g)),
             sprintf("%.2f", mean(deg)),
             as.character(median(deg)), as.character(max(deg)),
             as.character(sum(deg > 10)),
             as.character(comp$no),
             sprintf("%s (%.1f%%)", comma(max(comp$csize)),
                     100*max(comp$csize)/vcount(g)),
             sprintf("%.4f", trans))
) %>%
  kable(caption = "Table 1. Basic properties of the Bluetooth proximity network.",
        align = c("l", "r"))
Table 1. Basic properties of the Bluetooth proximity network.
Metric Value
Nodes (|V|) 662
Edges (|E|) 9,251
Density 0.0423
Mean degree 27.95
Median degree 26
Max degree 98
Nodes with degree > 10 553
Connected components 1
Largest component 662 (100.0%)
Global transitivity (clustering) 0.3474

The graph is a single connected component covering 100% of the nodes, density is 0.0423 (about 4 in 100 possible pairs are linked), and the mean degree is 28. The number worth flagging is global transitivity 0.347: roughly one in three two-hop paths closes into a triangle. For comparison, the calls network in HW1 had transitivity 0.24 and the Facebook layer in the same study has 0.24. Physical co-presence is the most clustered of the three CNS layers, which makes sense since students sit in the same lecture halls, share dorms, and study together. Triadic closure is the main way ties form in this network, so we should expect the heuristics that count common neighbours (Jaccard, Adamic-Adar, Resource Allocation) to do well.

1.3 Visualisation

set.seed(101)
comm_full <- cluster_louvain(g)
mem <- membership(comm_full)
sz  <- sizes(comm_full)
top10 <- as.integer(names(sort(sz, decreasing = TRUE))[1:10])
mem_lab <- ifelse(mem %in% top10,
                  paste0("C", match(mem, top10)),
                  "other")
V(g)$comm_lab <- factor(mem_lab,
                        levels = c(paste0("C", 1:10), "other"))

palette10 <- c(
  "#1b9e77","#d95f02","#7570b3","#e7298a","#66a61e",
  "#e6ab02","#a6761d","#666666","#1f78b4","#b2df8a","grey80")

ggraph(g, layout = "fr") +
  geom_edge_link(aes(alpha = weight), color = "grey60", show.legend = FALSE) +
  geom_node_point(aes(size = degree(g), fill = comm_lab),
                  shape = 21, stroke = 0.15, color = "grey20") +
  scale_edge_alpha(range = c(0.04, 0.55), trans = "log10") +
  scale_size_area(max_size = 5, guide = guide_legend(order = 1)) +
  scale_fill_manual(values = palette10, name = "Louvain community",
                    guide = guide_legend(order = 2,
                                         override.aes = list(size = 4))) +
  theme_void(base_size = 11) +
  labs(size = "Degree",
       title = sprintf("Copenhagen BT proximity (n = %d, m = %s, Q = %.3f)",
                       vcount(g), comma(ecount(g)), modularity(comm_full)))
Figure 1. Bluetooth proximity network (Fruchterman-Reingold layout). Nodes are colored by Louvain community (only the 10 largest are highlighted; smaller communities are shown in grey). Node area is proportional to degree, edge transparency to the log count of mutual scans.

Figure 1. Bluetooth proximity network (Fruchterman-Reingold layout). Nodes are colored by Louvain community (only the 10 largest are highlighted; smaller communities are shown in grey). Node area is proportional to degree, edge transparency to the log count of mutual scans. 

p1 <- tibble(k = deg) %>%
  ggplot(aes(k)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  labs(x = "Degree", y = "Number of nodes", subtitle = "Linear scale") +
  theme_minimal(base_size = 11)

# CCDF: P(K >= k). Smooth and well-behaved on log-log
ccdf <- tibble(k = sort(unique(deg[deg > 0]))) %>%
  mutate(p = sapply(k, function(x) mean(deg >= x)))

p2 <- ccdf %>%
  ggplot(aes(k, p)) +
  geom_step(color = "steelblue", linewidth = 0.8) +
  geom_point(color = "steelblue", size = 1.2) +
  scale_x_log10() + scale_y_log10(labels = scales::label_number()) +
  labs(x = "Degree k (log)", y = "P(K ≥ k) (log)",
       subtitle = "Complementary CDF (log-log)") +
  theme_minimal(base_size = 11)

gridExtra::grid.arrange(p1, p2, ncol = 2)
Figure 2. Degree distribution. Left: linear histogram. Right: complementary cumulative distribution function (1 - CDF) on log-log axes, which avoids the empty-bin artefacts that affect log-binned histograms.

Figure 2. Degree distribution. Left: linear histogram. Right: complementary cumulative distribution function (1 - CDF) on log-log axes, which avoids the empty-bin artefacts that affect log-binned histograms.

2. Step 1: Edge deletion and class table

We follow the Liben-Nowell & Kleinberg (2003) protocol: delete a fraction of real edges from \(G\) to get \(G'\), treat the deleted endpoints as the positive class (label = 1), and sample a matched negative class (label = 0) of non-edges that the classifier should reject. We delete nlinks = 5000 edges, slightly more than half the graph, which leaves the heuristics enough deleted ties to actually recover.

The standard caveat from Workshop 3 applies to negative sampling. A uniformly random pair of nodes is almost always disconnected in a sparse network, so negatives drawn that way are trivial to classify and the model just learns that “isolated pairs are not friends”, which is not what we want. We use Workshop 3’s solution and sample negatives only between well-connected nodes (degree > 10 in \(G\)), which forces the classifier to draw a real distinction. With 553 well-connected nodes (out of 662) there are plenty of candidates.

set.seed(42)
nlinks <- 5000
indexes_deleted   <- sample(seq_len(ecount(g)), nlinks)
deleted_endpoints <- ends(g, E(g)[indexes_deleted])

true_edges <- data.frame(
  X1 = deleted_endpoints[, 1],
  X2 = deleted_endpoints[, 2],
  stringsAsFactors = FALSE
)

Gprime <- delete_edges(g, indexes_deleted)
most_connected <- V(g)[degree(g) > 10]
mc_names       <- as_ids(most_connected)

false_edges_list <- vector("list", nlinks)
filled <- 0L; attempts <- 0L
while (filled < nlinks) {
  attempts <- attempts + 1L
  pair <- sample(mc_names, 2L, replace = FALSE)
  if (!are_adjacent(g, pair[1], pair[2])) {
    filled <- filled + 1L
    false_edges_list[[filled]] <- pair
  }
}
false_edges <- as.data.frame(do.call(rbind, false_edges_list),
                             stringsAsFactors = FALSE)
colnames(false_edges) <- c("X1", "X2")
true_edges  <- data.frame(true_edges,  obs = 1L)
false_edges <- data.frame(false_edges, obs = 0L)

total_edges <- rbind(true_edges, false_edges)
colnames(total_edges)[1:2] <- c("id1", "id2")
Table 2. Class table for Step 1: positive class (deleted real edges) and negative class (sampled non-edges between well-connected nodes).
Quantity Value
Edges in G 9,251
Edges in G’ after deletion 4,251
Edges deleted (positive class) 5,000
Well-connected nodes (deg > 10) 553
Hard non-edges sampled (negative class) 5,000
Rejection-sampling attempts 5,285
Rejection rate 5.39%
Total candidate pairs 10,000

We end up with a balanced dataset of 10,000 candidate pairs: 5,000 deleted-real edges and 5,000 hard non-edges.

We end up with a balanced dataset of 10,000 candidate pairs, half deleted-real and half hard non-edges.

3. Step 2: Heuristic features

We compute seven structural similarity heuristics on \(G'\), the observed graph (not on \(G\), since computing on \(G\) would leak the deleted edges back into the features through the degrees and neighbourhoods of their endpoints).

# Heuristic Formula Lecture 4 slide
1 Common neighbors (CN) \(\lvert\Gamma(x)\cap\Gamma(y)\rvert\) 20
2 Jaccard coefficient \(\frac{\lvert\Gamma(x)\cap\Gamma(y)\rvert}{\lvert\Gamma(x)\cup\Gamma(y)\rvert}\) 21
3 Adamic-Adar (AA) \(\sum_{z\in\Gamma(x)\cap\Gamma(y)}\frac{1}{\log k_z}\) 22
4 Preferential Attachment (PA) \(k_x\cdot k_y\) 23
5 Graph distance (inverse) \(1/d_{xy}\) 19
6 Resource Allocation (RA) \(\sum_{z\in\Gamma(x)\cap\Gamma(y)}\frac{1}{k_z}\) 33
7 Same community (Louvain) \(\mathbf{1}[c(x)=c(y)]\) 32

The first four (CN, Jaccard, AA, PA) are the ones Lecture 4 uses in its “Link prediction algorithm” example (slide 28) and the ones we implemented in Workshop 3. The other three (graph distance, RA, same-community) come from Lecture 4 slides 19, 32, and 33; the lecture covers them but Workshop 3 did not. We add them so we can check whether the community structure itself adds predictive value beyond what the local neighbourhood heuristics already capture.

# Louvain on G' (not G) so we don't leak the deleted edges.
# Louvain is stochastic; setting the seed makes it reproducible.
set.seed(321)
comm <- cluster_louvain(Gprime)
V(Gprime)$comm <- membership(comm)
Quantity Value
Louvain communities on G’ 33
Modularity Q 0.7564 
# Pre-compute neighborhoods on Gprime once, distances on demand.
n1     <- neighborhood(Gprime, order = 1, nodes = total_edges$id1)
n2     <- neighborhood(Gprime, order = 1, nodes = total_edges$id2)
deg_Gp <- degree(Gprime)
comm_v <- V(Gprime)$comm
names(comm_v) <- as_ids(V(Gprime))

# Distance matrix from the rows we need is too expensive to precompute
# fully (~440k entries), but distances() with from/to is fast.
total_edges$sim_cn   <- 0L
total_edges$sim_jacc <- 0
total_edges$sim_aa   <- 0
total_edges$sim_pref <- 0L
total_edges$sim_dist <- 0
total_edges$sim_ra   <- 0
total_edges$sim_comm <- 0L

t0 <- Sys.time()
for (i in seq_len(nrow(total_edges))) {
  v1 <- total_edges$id1[i]; v2 <- total_edges$id2[i]
  nei_a <- setdiff(as_ids(n1[[i]]), v1)
  nei_b <- setdiff(as_ids(n2[[i]]), v2)
  inter <- intersect(nei_a, nei_b)
  uni   <- union(nei_a, nei_b)

  total_edges$sim_cn[i]   <- length(inter)
  total_edges$sim_jacc[i] <- if (length(uni) == 0) 0 else length(inter)/length(uni)

  if (length(inter) == 0) {
    total_edges$sim_aa[i] <- 0
    total_edges$sim_ra[i] <- 0
  } else {
    k_z <- deg_Gp[inter]
    aa_kz <- k_z[k_z > 1]                # 1/log(1) is undefined; skip
    total_edges$sim_aa[i] <- if (length(aa_kz) == 0) 0 else sum(1/log(aa_kz))
    ra_kz <- k_z[k_z > 0]
    total_edges$sim_ra[i] <- if (length(ra_kz) == 0) 0 else sum(1/ra_kz)
  }

  total_edges$sim_pref[i] <- length(nei_a) * length(nei_b)

  # Topological (unweighted) hop count, per Lecture 4 slide 19.
  # weights = NA forces igraph to use BFS instead of Dijkstra on the
  # weighted graph. The "weight" edge attribute is co-presence count, not
  # a cost we want to walk along.
  d <- distances(Gprime, v = v1, to = v2, weights = NA)[1, 1]
  total_edges$sim_dist[i] <- if (is.infinite(d) | d == 0) 0 else 1/d

  total_edges$sim_comm[i] <- as.integer(comm_v[v1] == comm_v[v2])
}
elapsed <- as.numeric(difftime(Sys.time(), t0, units = "secs"))
total_edges %>%
  group_by(obs) %>%
  summarise(
    `Common neighbors`     = mean(sim_cn),
    Jaccard                = mean(sim_jacc),
    `Adamic-Adar`          = mean(sim_aa),
    `Pref. Attachment`     = mean(sim_pref),
    `Graph distance (1/d)` = mean(sim_dist),
    `Resource Allocation`  = mean(sim_ra),
    `Same community`       = mean(sim_comm),
    .groups = "drop"
  ) %>%
  mutate(Class = ifelse(obs == 1, "Positive (deleted)", "Negative (non-edge)")) %>%
  select(Class, everything(), -obs) %>%
  kable(digits = 3,
        caption = "Table 3. Mean of each heuristic by class. Positive class are the 5,000 deleted real edges; negative class are the 5,000 hard non-edges.")
Table 3. Mean of each heuristic by class. Positive class are the 5,000 deleted real edges; negative class are the 5,000 hard non-edges.
Class Common neighbors Jaccard Adamic-Adar Pref. Attachment Graph distance (1/d) Resource Allocation Same community
Negative (non-edge) 0.307 0.010 0.101 211.667 0.354 0.015 0.040
Positive (deleted) 2.783 0.093 0.974 312.422 0.468 0.165 0.463 
heuristic_cols <- c("sim_cn","sim_jacc","sim_aa","sim_pref",
                    "sim_dist","sim_ra","sim_comm")
pretty_names <- c("CN","Jaccard","AA","PA","1/d","RA","Same comm")

cor_mat <- cor(total_edges[, heuristic_cols])
dimnames(cor_mat) <- list(pretty_names, pretty_names)

cor_df <- as.data.frame(cor_mat)
cor_df <- cbind(`  ` = rownames(cor_df), cor_df)
kable(cor_df, digits = 2, row.names = FALSE,
      caption = "Table 4. Pairwise Pearson correlation among the seven heuristics. CN, Jaccard, AA and RA share the common-neighbour numerator and are highly collinear (r > 0.85). PA and the binary same-community indicator are largely independent of that family.")
Table 4. Pairwise Pearson correlation among the seven heuristics. CN, Jaccard, AA and RA share the common-neighbour numerator and are highly collinear (r > 0.85). PA and the binary same-community indicator are largely independent of that family.
CN Jaccard AA PA 1/d RA Same comm
CN 1.00 0.87 0.99 0.48 0.66 0.92 0.57
Jaccard 0.87 1.00 0.89 0.17 0.66 0.88 0.64
AA 0.99 0.89 1.00 0.44 0.65 0.97 0.59
PA 0.48 0.17 0.44 1.00 0.36 0.35 0.04
1/d 0.66 0.66 0.65 0.36 1.00 0.62 0.45
RA 0.92 0.88 0.97 0.35 0.62 1.00 0.58
Same comm 0.57 0.64 0.59 0.04 0.45 0.58 1.00

The CN / Jaccard / AA / RA group is collinear by definition: they all sum over the same intersection of neighbours and only differ in how they weight the terms. Preferential Attachment is a different object (degree product), and the same-community dummy is mostly independent of the local heuristics. Because of this collinearity we report three different “most important” measures in Section 5: single-feature AUC (heuristic by heuristic), GLM standardised coefficients (which control for the rest), and Random Forest permutation importance (which picks up non-redundant predictive value).

The class means look right. Positives (deleted real edges) have more shared neighbours, share Louvain communities far more often, and have larger PA products. Preferential Attachment has the smallest relative gap between the two classes, so we already suspect it will be the weakest predictor, which is what Lecture 4 slide 30 reports.

fill_pal <- c("Random non-edge" = "#7570b3",
              "True (deleted) edge" = "#d95f02")

total_edges %>%
  pivot_longer(c(sim_cn, sim_jacc, sim_aa, sim_ra, sim_pref),
               names_to = "heuristic", values_to = "score") %>%
  mutate(class = ifelse(obs == 1, "True (deleted) edge", "Random non-edge"),
         heuristic = recode(heuristic,
                            sim_cn   = "Common Neighbors",
                            sim_jacc = "Jaccard",
                            sim_aa   = "Adamic-Adar",
                            sim_pref = "Preferential Attachment",
                            sim_ra   = "Resource Allocation"),
         heuristic = factor(heuristic,
                            levels = c("Common Neighbors","Jaccard",
                                       "Adamic-Adar","Resource Allocation",
                                       "Preferential Attachment"))) %>%
  ggplot(aes(score + 1e-3, fill = class)) +
  geom_histogram(alpha = 0.6, position = "identity", bins = 40, color = NA) +
  scale_x_log10(labels = scales::label_log()) +
  facet_wrap(~ heuristic, scales = "free", ncol = 2) +
  scale_fill_manual(values = fill_pal) +
  labs(x = "Score (log10, with +1e-3 offset)",
       y = "Number of pairs", fill = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")
Figure 3a. Continuous heuristics: score distributions by class on a log10 x-axis (small offset added so zeros are visible). The further apart the two coloured masses are, the better the heuristic separates real ties from non-ties.

Figure 3a. Continuous heuristics: score distributions by class on a log10 x-axis (small offset added so zeros are visible). The further apart the two coloured masses are, the better the heuristic separates real ties from non-ties. 

# Convert 1/d back to integer hop count for a readable axis.
# d = 0 is reserved for disconnected pairs.
dist_summary <- total_edges %>%
  mutate(dist_int = ifelse(sim_dist > 0, as.integer(round(1/sim_dist)), 0L),
         dist_lab = ifelse(dist_int == 0, "disconnected", as.character(dist_int)),
         class    = ifelse(obs == 1, "True (deleted) edge", "Random non-edge")) %>%
  count(class, dist_int, dist_lab)

ordered_levels <- c(
  if (any(dist_summary$dist_int == 0)) "disconnected" else NULL,
  as.character(sort(setdiff(unique(dist_summary$dist_int), 0L)))
)

dist_summary %>%
  mutate(dist_lab = factor(dist_lab, levels = ordered_levels)) %>%
  ggplot(aes(dist_lab, n, fill = class)) +
  geom_col(position = "dodge") +
  scale_fill_manual(values = fill_pal) +
  labs(x = "Shortest-path hop count d in G'",
       y = "Number of pairs", fill = NULL,
       title = "Graph distance") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none") -> p_dist

total_edges %>%
  mutate(class = ifelse(obs == 1, "True (deleted) edge", "Random non-edge"),
         same  = ifelse(sim_comm == 1, "Same community", "Different community")) %>%
  count(class, same) %>%
  ggplot(aes(same, n, fill = class)) +
  geom_col(position = "dodge") +
  scale_fill_manual(values = fill_pal) +
  labs(x = NULL, y = "Count", fill = NULL,
       title = "Same Louvain community") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom") -> p_comm

gridExtra::grid.arrange(p_dist, p_comm, ncol = 2)
Figure 3b. Discrete heuristics shown on linear axes. Graph distance takes only a handful of integer values (so 1/d is also discrete) and same-community is binary, which is why a log axis is not informative for them.

Figure 3b. Discrete heuristics shown on linear axes. Graph distance takes only a handful of integer values (so 1/d is also discrete) and same-community is binary, which is why a log axis is not informative for them.

4. Step 3: Binary classifier with cross-validation

The Homework 2 spec calls for cross-validation, going beyond Workshop 3’s single train/test split. We use 10-fold cross-validation via caret::train, holding out the same observations across all seven features so feature comparisons are fair.

We fit two models:

  • Logistic regression (glm with family = binomial(link = "logit")), the model from Workshop 3 and Lecture 4. The coefficients tell us which feature still matters after we control for the others.
  • Random Forest, which gives us permutation importance. This is threshold-free and handles correlated features better than the GLM (Jaccard and Adamic-Adar share the common-neighbour numerator, for example).
# caret needs a factor outcome for classification with ROC metric.
total_edges$obs_f <- factor(total_edges$obs, levels = c(0, 1),
                            labels = c("no", "yes"))

ctrl <- trainControl(
  method          = "cv",
  number          = 10,
  classProbs      = TRUE,
  summaryFunction = twoClassSummary,
  savePredictions = "final"
)
set.seed(42)
fit_glm <- train(
  obs_f ~ sim_cn + sim_jacc + sim_aa + sim_pref + sim_dist + sim_ra + sim_comm,
  data    = total_edges,
  method  = "glm",
  family  = binomial(link = "logit"),
  metric  = "ROC",
  trControl = ctrl
)
Table 5. Logistic regression: 10-fold cross-validated performance.
Model CV folds AUC (ROC) Sensitivity Specificity
Logistic regression (glm) 10 0.8810 0.8820 0.7504
Table 6. Logistic-regression coefficients (final model fit on all 10,000 pairs after CV). Same community has by far the largest |z|; the four common-neighbour heuristics share variance and so most of them are not individually significant once the others are in the model.
Feature Estimate Std. error z value p value
(Intercept) -2.521 0.165 -15.24 1.8e-52
Common neighbors -2.713 1.616 -1.68 0.0932
Jaccard 11.689 2.380 4.91 9.0e-07
Adamic-Adar 12.874 7.151 1.80 0.0718
Preferential Attachment 0.001 0.000 3.14 0.0017
Graph distance (1/d) 2.637 0.480 5.49 4.0e-08
Resource Allocation -19.988 14.976 -1.33 0.1820
Same community 1.526 0.095 16.09 2.8e-58 
set.seed(42)
fit_rf <- train(
  obs_f ~ sim_cn + sim_jacc + sim_aa + sim_pref + sim_dist + sim_ra + sim_comm,
  data    = total_edges,
  method  = "rf",
  metric  = "ROC",
  trControl = ctrl,
  tuneGrid  = data.frame(mtry = c(2, 3, 4)),
  ntree     = 300,
  importance = TRUE
)
Table 7. Random Forest: 10-fold cross-validated performance across mtry (number of features sampled at each split). mtry = 4 was selected as the final model.
mtry AUC (ROC) Sensitivity Specificity
2 0.8829 0.8706 0.7798
3 0.8824 0.8678 0.7678
4 0.8831 0.8682 0.7676

5. Step 4: Precision evaluation and most-important heuristic

5.1 Confusion matrix

We use the cross-validated out-of-fold predictions from caret, then evaluate the GLM at two probability thresholds: the Workshop 3 default of 0.30 and the textbook 0.50.

oof <- fit_glm$pred %>%
  arrange(rowIndex) %>%
  mutate(obs01 = ifelse(obs == "yes", 1L, 0L))

cm_summary <- function(thr) {
  pred_thr <- factor(ifelse(oof$yes > thr, "yes", "no"), levels = c("no","yes"))
  cm <- confusionMatrix(pred_thr, oof$obs, positive = "yes")
  list(
    table = as.data.frame.matrix(cm$table),
    metrics = c(Accuracy    = unname(cm$overall["Accuracy"]),
                Sensitivity = unname(cm$byClass["Sensitivity"]),
                Specificity = unname(cm$byClass["Specificity"]),
                Kappa       = unname(cm$overall["Kappa"]),
                F1          = unname(cm$byClass["F1"]))
  )
}

cm_03 <- cm_summary(0.3)
cm_05 <- cm_summary(0.5)
Table 8a. Confusion matrix at threshold = 0.30 (rows = predicted, columns = actual).
Predicted  Actual no yes
no 3864 755
yes 1136 4245
Table 9a. Aggregate metrics at threshold = 0.30.
Metric Value
Accuracy 0.811
Sensitivity (recall on positives) 0.849
Specificity (recall on negatives) 0.773
Cohen’s kappa 0.622
F1 score 0.818
Table 8b. Confusion matrix at threshold = 0.50 (rows = predicted, columns = actual).
Predicted  Actual no yes
no 4410 1248
yes 590 3752
Table 9b. Aggregate metrics at threshold = 0.50.
Metric Value
Accuracy 0.816
Sensitivity (recall on positives) 0.750
Specificity (recall on negatives) 0.882
Cohen’s kappa 0.632
F1 score 0.803

5.2 Per-feature AUC and the combined GLM

Following Lecture 4 slide 30 (which compares heuristics by their relative performance vs. random and vs. common neighbors), we compute the rank-based AUC for every single heuristic and for the fitted GLM.

auc_simple <- function(score, label) {
  pos <- score[label == 1]; neg <- score[label == 0]
  rk  <- rank(c(pos, neg))
  (sum(rk[seq_along(pos)]) - length(pos)*(length(pos)+1)/2) /
    (length(pos) * length(neg))
}
y <- total_edges$obs

rf_oof <- fit_rf$pred %>%
  filter(mtry == fit_rf$bestTune$mtry) %>%
  arrange(rowIndex) %>%
  mutate(obs01 = ifelse(obs == "yes", 1L, 0L))

auc_tbl <- tibble(
  Predictor = c("Common neighbors", "Jaccard", "Adamic-Adar",
                "Preferential Attachment", "Graph distance (1/d)",
                "Resource Allocation", "Same community",
                "GLM (all 7 features, OOF)", "Random Forest (all 7 features, OOF)"),
  Type      = c(rep("single heuristic", 7), "combined model", "combined model"),
  AUC = c(
    auc_simple(total_edges$sim_cn,   y),
    auc_simple(total_edges$sim_jacc, y),
    auc_simple(total_edges$sim_aa,   y),
    auc_simple(total_edges$sim_pref, y),
    auc_simple(total_edges$sim_dist, y),
    auc_simple(total_edges$sim_ra,   y),
    auc_simple(total_edges$sim_comm, y),
    auc_simple(oof$yes, oof$obs01),
    auc_simple(rf_oof$yes, rf_oof$obs01)
  )
) %>% arrange(desc(AUC))
Table 10. Area under the ROC curve for each predictor, sorted descending. The top four single heuristics tie at 0.863 to three decimals; combined classifiers gain a few points by using their joint information.
Predictor Type AUC
Random Forest (all 7 features, OOF) combined model 0.8830
GLM (all 7 features, OOF) combined model 0.8811
Adamic-Adar single heuristic 0.8635
Jaccard single heuristic 0.8628
Resource Allocation single heuristic 0.8627
Common neighbors single heuristic 0.8588
Graph distance (1/d) single heuristic 0.8137
Same community single heuristic 0.7118
Preferential Attachment single heuristic 0.6210

One thing to clarify about this comparison: the individual heuristics have no fitted parameters, so their rank-AUC is the same on the training data and on a held-out fold (no overfitting is possible). The combined GLM and RF rows use out-of-fold predictions from the 10-fold CV, while the single-heuristic rows use the full dataset. The comparison is still fair because all we measure is how well a score separates positives from negatives.

5.3 ROC curves

roc_tbl <- function(score, label, name) {
  ord <- order(-score)
  y_o <- label[ord]
  tibble(
    predictor = name,
    fpr = c(0, cumsum(y_o == 0) / sum(y_o == 0)),
    tpr = c(0, cumsum(y_o == 1) / sum(y_o == 1))
  )
}

# Compute AUC for every series, sort descending so the legend tracks performance
roc_inputs <- list(
  "Common neighbors"        = list(s = total_edges$sim_cn,   l = y),
  "Jaccard"                 = list(s = total_edges$sim_jacc, l = y),
  "Adamic-Adar"             = list(s = total_edges$sim_aa,   l = y),
  "Resource Allocation"     = list(s = total_edges$sim_ra,   l = y),
  "Preferential Attachment" = list(s = total_edges$sim_pref, l = y),
  "Same community"          = list(s = total_edges$sim_comm, l = y),
  "GLM (combined)"          = list(s = oof$yes,              l = oof$obs01),
  "Random Forest"           = list(s = rf_oof$yes,           l = rf_oof$obs01)
)

aucs <- sapply(roc_inputs, function(x) auc_simple(x$s, x$l))
order_idx <- order(-aucs)
labels_in_order <- sprintf("%s (AUC = %.3f)",
                           names(roc_inputs)[order_idx],
                           aucs[order_idx])

roc_df <- bind_rows(lapply(seq_along(roc_inputs), function(i) {
  x <- roc_inputs[[i]]
  roc_tbl(x$s, x$l,
          sprintf("%s (AUC = %.3f)", names(roc_inputs)[i], aucs[i]))
})) %>%
  mutate(predictor = factor(predictor, levels = labels_in_order))

ggplot(roc_df, aes(fpr, tpr, color = predictor)) +
  geom_abline(slope = 1, intercept = 0, lty = 2, color = "grey60") +
  geom_line(linewidth = 0.9) +
  coord_equal() +
  scale_color_viridis_d(option = "D", end = 0.92) +
  labs(x = "False positive rate", y = "True positive rate", color = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "right", legend.text = element_text(size = 9))
Figure 4. ROC curves for the seven single heuristics and the two combined classifiers (GLM and Random Forest). Legend is sorted by AUC. Combined-model curves use out-of-fold cross-validated predictions; single-heuristic curves use full-sample scores (no fitting, so no overfitting risk).

Figure 4. ROC curves for the seven single heuristics and the two combined classifiers (GLM and Random Forest). Legend is sorted by AUC. Combined-model curves use out-of-fold cross-validated predictions; single-heuristic curves use full-sample scores (no fitting, so no overfitting risk).

5.4 Random Forest permutation importance

Permutation importance is the usual way to answer “which feature matters most” when the features are correlated. The Random Forest shuffles each feature in turn and measures how much the accuracy drops.

imp <- importance(fit_rf$finalModel, type = 1)  # type 1 = mean decrease in accuracy
imp_df <- tibble(
  feature = rownames(imp),
  importance = as.numeric(imp[, 1])
) %>%
  mutate(Feature = recode(feature,
                          sim_cn   = "Common Neighbors",
                          sim_jacc = "Jaccard",
                          sim_aa   = "Adamic-Adar",
                          sim_pref = "Preferential Attachment",
                          sim_dist = "Graph distance (1/d)",
                          sim_ra   = "Resource Allocation",
                          sim_comm = "Same community"),
         Family = case_when(
           feature %in% c("sim_cn","sim_jacc","sim_aa","sim_ra") ~ "Common-neighbour family",
           feature == "sim_pref"                                  ~ "Degree-based",
           feature == "sim_dist"                                  ~ "Path-based",
           feature == "sim_comm"                                  ~ "Community-based"
         )) %>%
  arrange(desc(importance))

ggplot(imp_df, aes(reorder(Feature, importance), importance, fill = Family)) +
  geom_col() +
  coord_flip() +
  scale_fill_manual(values = c("Common-neighbour family" = "#1b9e77",
                               "Degree-based"            = "#d95f02",
                               "Community-based"         = "#7570b3",
                               "Path-based"              = "#666666")) +
  labs(x = NULL, y = "Mean decrease in accuracy",
       fill = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")
Figure 5. Random Forest permutation importance (mean decrease in accuracy when each feature is shuffled). Bars are coloured by category: structural local heuristics, mesoscopic structure, and degree-based features. Single-fit values; run-to-run variance from the RF training step is not shown.

Figure 5. Random Forest permutation importance (mean decrease in accuracy when each feature is shuffled). Bars are coloured by category: structural local heuristics, mesoscopic structure, and degree-based features. Single-fit values; run-to-run variance from the RF training step is not shown.

Table 11. Random Forest permutation importance, ranked. The forest assigns the largest non-redundant weight to Preferential Attachment, even though PA is the weakest single-feature predictor by AUC. The result is consistent with PA contributing information that is orthogonal to the common-neighbour family and to community membership.
Rank Feature Family Importance
1 Preferential Attachment Degree-based 51.9
2 Resource Allocation Common-neighbour family 29.9
3 Same community Community-based 28.4
4 Jaccard Common-neighbour family 27.8
5 Adamic-Adar Common-neighbour family 26.2
6 Common Neighbors Common-neighbour family 8.1
7 Graph distance (1/d) Path-based 6.0

5.5 Discussion: which heuristic is most important, and why

single_heuristics <- c("Common neighbors","Jaccard","Adamic-Adar",
                       "Preferential Attachment","Graph distance (1/d)",
                       "Resource Allocation","Same community")

winner_auc <- auc_tbl$predictor[auc_tbl$predictor %in% single_heuristics][1]

glm_coef <- summary(fit_glm$finalModel)$coefficients
glm_z <- abs(glm_coef[-1, "z value"])              # drop intercept
winner_glm <- names(which.max(glm_z))
winner_glm_pretty <- recode(winner_glm,
  sim_cn   = "Common Neighbors",  sim_jacc = "Jaccard",
  sim_aa   = "Adamic-Adar",       sim_pref = "Preferential Attachment",
  sim_dist = "Graph distance",    sim_ra   = "Resource Allocation",
  sim_comm = "Same community")

cat("- **Single best heuristic by AUC:** ",       winner_auc,           "\n")
  • Single best heuristic by AUC:
cat("- **Largest GLM |z|-statistic:** ",          winner_glm_pretty,    "\n")
  • Largest GLM |z|-statistic: Same community
cat("- **Top RF permutation importance:** ",      imp_df$feature[1],    "\n")
  • Top RF permutation importance: sim_pref

The three measures give three different answers, and the disagreement is informative:

  1. Single-feature AUC: Adamic-Adar (about 0.86). AA wins the one-feature ranking, though Jaccard, Resource Allocation, and plain Common Neighbors are all within 0.005 AUC of AA and of each other. This is what Lecture 4 slide 30 predicts: the common-neighbour heuristics work together, and any one of them is a reasonable choice. Same community is at 0.71, graph distance at 0.70, and Preferential Attachment is last at 0.62, which is what Liben-Nowell & Kleinberg (2003) report too.
  2. Largest GLM coefficient: Same community. Once Jaccard, Adamic-Adar, RA, and Common Neighbors all enter the model they compete for the same variance, since they are collinear (see the correlation matrix in Section 3, where these four features have pairwise \(r > 0.7\)). The GLM resolves the competition by putting most of the weight on the one feature outside the common-neighbour family: the binary “same Louvain community” indicator. Its z-statistic is about twice the next largest coefficient.
  3. Top RF permutation importance: Preferential Attachment. This surprised us. PA was the worst single heuristic in Lecture 4 slide 30, and it is also the worst single AUC here (0.62). But the Random Forest ranks it as the most important non-redundant feature: shuffling PA hurts accuracy more than shuffling anything else. PA carries information the other features do not. The product of degrees is more or less independent of common neighbours, community membership, and distance, so the forest uses it to split cases the neighbour-based features cannot.

So “important” depends on what you ask. A heuristic can score high by itself (AA), it can keep its weight once the others are in the model (same-community), or it can carry information the others miss (PA). The three criteria answer different questions. For the question on the homework brief (“Which heuristic is the most important? Why?”) our answer is Adamic-Adar on the marginal criterion, since it has the highest single-feature AUC. Two caveats: (i) AA, Jaccard and Resource Allocation tie at AUC = 0.863 to three decimals, so picking any one of the four as “the” winner overstates what the data show; and (ii) the answer changes under the other two criteria, which is why we report all three.

The results map onto the mechanisms from Lecture 4 slide 6 (Rivera, Soderstrom & Uzzi 2010). The strong performance of the common-neighbour family (AA, Jaccard, RA, CN, all within 0.005 AUC of each other) is triadic closure (slide 12): students who share friends tend to be tied themselves. AA does slightly better than plain CN because it down-weights common neighbours of high-degree hubs. Sharing a niche friend with degree 5 says more than sharing a hub with degree 80, which is the same idea behind Granovetter’s weak ties. The same-community result in the GLM is community-level homophily (slide 11). PA’s win in the RF is the degree mechanism (slide 6.2.4), not as a marginal predictor but as a residual signal the other features do not capture. So all three mechanisms leave traces in the data, not just triadic closure.

6. Step 5: Improvements

The model performs well on this dataset, but Lecture 4 (slides 32-35) and Workshop 3’s caveats together prescribe a clear menu of extensions. We sketch six.

  1. Higher-order structural features (Katz, SimRank, Rooted PageRank, Local Path). Lecture 4 slide 24 shows Katz, which weights the number of paths between \(x\) and \(y\) by their length. For our network we used 1-hop neighbourhoods only; Katz would let pairs that are 2 or 3 hops apart contribute predictive signal. Wang et al. (2015, slide 33) survey the broader family.
  2. Dispersion (Lecture 4 slide 32, after Backstrom & Kleinberg 2014). Standard heuristics treat all common neighbours equivalently; dispersion penalises common neighbours that are themselves close to each other in the residual graph, which distinguishes “couple-style” ties from “same-clique” ties.
  3. Non-network features (Lecture 4 slide 34). The CNS bundle includes a genders.csv file. A binary “same gender” feature captures the homophily mechanism of Lecture 4 slide 11 directly. Adding it to the GLM would test whether sociodemographic homophily explains residual variance not picked up by structure alone.
  4. Network embeddings / GNNs (Lecture 4 slide 35). Node2vec or GraphSAGE would learn a vector for each node from random walks on \(G'\) and let a classifier work in the embedding space. In practice these methods help most on sparse networks and tend not to change rankings on dense, clustered graphs like ours.
  5. Temporal split instead of random deletion (Workshop 3 caveat #1). We have a timestamp on every BT scan. A more realistic evaluation deletes edges that first appeared in the last week of the study, which is closer to the real task of predicting future ties from past observation. We would expect AUCs to go up, because real future links are concentrated among friends-of-friends rather than spread uniformly across the graph.
  6. Multiplex link prediction. The CNS has four layers (calls, SMS, FB friends, BT proximity). Using the calls layer (the HW1 network) as a feature for BT-tie prediction (i.e. “did these two ever call each other?”) would test the cross-layer correlation discussed in Sapiezynski et al. (2019). Subject to ID overlap, this is a cheap addition with a clear payoff.

Sources and acknowledgments

Use of Artificial Intelligence. Claude Sonnet 4.7 / Opus 4.7 (Anthropic, 2026) was used as a programming assistant for: the data-acquisition pipeline, the seven-heuristic feature loop, the caret cross-validation setup, and prose editing of methodology paragraphs. All methodological choices, threshold defenses, and substantive interpretations are our own.

Data. Sapiezynski, P., Stopczynski, A., Lassen, D. D., & Lehmann, S. (2019). Interaction data from the Copenhagen Networks Study [Data set]. figshare. https://doi.org/10.6084/m9.figshare.7267433

References.

  • Adamic, L. A., & Adar, E. (2003). Friends and Neighbors on the Web. Social Networks, 25(3), 211-230.
  • Backstrom, L., & Kleinberg, J. (2014). Romantic partnerships and the dispersion of social ties. CSCW ’14.
  • Granovetter, M. (1973). The Strength of Weak Ties. American Journal of Sociology, 78(6), 1360-1380.
  • Liben-Nowell, D., & Kleinberg, J. (2003). The Link Prediction Problem for Social Networks. CIKM ’03.
  • Lü, L., & Zhou, T. (2011). Link prediction in complex networks: a survey. Physica A, 390(6), 1150-1170.
  • Rivera, M. T., Soderstrom, S. B., & Uzzi, B. (2010). Dynamics of Dyads in Social Networks. Annual Review of Sociology, 36, 91-115.
  • Sapiezynski, P., Stopczynski, A., Lassen, D. D., & Lehmann, S. (2019). Interaction data from the Copenhagen Networks Study. Scientific Data, 6, 315.
  • Sekara, V., & Lehmann, S. (2014). The strength of friendship ties in proximity sensor data. PLoS ONE, 9(7), e100915.
  • Wang, P., Xu, B., Wu, Y., & Zhou, X. (2015). Link prediction in social networks: the state-of-the-art. Science China Information Sciences, 58(1), 1-38.

Group submission: Tone Varberg Sabri, Abdullah Tadmuri, Maiheliya Maimaitimin. Master in Computational Social Science, UC3M.

Reproducibility

sessionInfo()
R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.4.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] C.UTF-8/UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8

time zone: Europe/Madrid
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] randomForest_4.7-1.2 pROC_1.19.0.1        caret_7.0-1         
 [4] lattice_0.22-7       scales_1.4.0         ggraph_2.2.2        
 [7] igraph_2.2.0         lubridate_1.9.4      forcats_1.0.0       
[10] stringr_1.6.0        dplyr_1.1.4          purrr_1.2.0         
[13] readr_2.1.5          tidyr_1.3.1          tibble_3.3.0        
[16] ggplot2_4.0.0        tidyverse_2.0.0      knitr_1.50          

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1     viridisLite_0.4.2    timeDate_4041.110   
 [4] farver_2.1.2         viridis_0.6.5        S7_0.2.0            
 [7] fastmap_1.2.0        tweenr_2.0.3         digest_0.6.37       
[10] rpart_4.1.24         timechange_0.3.0     lifecycle_1.0.5     
[13] survival_3.8-3       magrittr_2.0.4       compiler_4.5.1      
[16] rlang_1.1.7          sass_0.4.10          tools_4.5.1         
[19] yaml_2.3.10          data.table_1.17.2    labeling_0.4.3      
[22] graphlayouts_1.2.2   plyr_1.8.9           RColorBrewer_1.1-3  
[25] withr_3.0.2          stats4_4.5.1         nnet_7.3-20         
[28] grid_4.5.1           polyclip_1.10-7      e1071_1.7-16        
[31] future_1.67.0        globals_0.18.0       iterators_1.0.14    
[34] MASS_7.3-65          cli_3.6.5            rmarkdown_2.30      
[37] generics_0.1.4       future.apply_1.20.0  reshape2_1.4.5      
[40] tzdb_0.5.0           proxy_0.4-27         cachem_1.1.0        
[43] ggforce_0.5.0        splines_4.5.1        parallel_4.5.1      
[46] vctrs_0.6.5          hardhat_1.4.2        Matrix_1.7-3        
[49] jsonlite_2.0.0       hms_1.1.3            ggrepel_0.9.6       
[52] listenv_0.9.1        foreach_1.5.2        gower_1.0.2         
[55] jquerylib_0.1.4      recipes_1.3.1        parallelly_1.45.1   
[58] glue_1.8.0           codetools_0.2-20     stringi_1.8.7       
[61] gtable_0.3.6         pillar_1.11.1        htmltools_0.5.8.1   
[64] ipred_0.9-15         lava_1.8.1           R6_2.6.1            
[67] tidygraph_1.3.1      evaluate_1.0.3       memoise_2.0.1       
[70] bslib_0.9.0          class_7.3-23         Rcpp_1.1.1          
[73] gridExtra_2.3        nlme_3.1-168         prodlim_2025.04.28  
[76] xfun_0.52            ModelMetrics_1.2.2.2 pkgconfig_2.0.3