library(tidyverse)
library(igraph)
library(ggraph)
library(scales)
library(ggthemes) # theme_hc(): the workshop plotting theme
library(caret) # 10-fold CV, mirrors Homework 2
library(randomForest) # RF, mirrors Homework 2
library(cowplot)
library(gridExtra)
library(data.table) # fast read for the arXiv ca-GrQc edgelist (§8)
theme_set(theme_minimal(base_size = 11))
Abstract
We study SIR spreading on the Copenhagen Bluetooth proximity network
(662 nodes, 9,251 edges) and replicate every analysis on the arXiv
ca-GrQc collaboration network (4,158 nodes, 13,422 edges). The
heterogeneous-mean-field threshold is \(\beta_c=0.00265\) at \(\mu=0.1\). At a data-chosen operating point
\(\beta^*=2\beta_c\), closeness-based
1% seeds infect 74 more nodes than random (\(p<0.001\)) and bring the prevalence peak
forward by 4 steps; an 8-step quarantine of 20% of susceptibles,
triggered at the detected epidemic-doubling time, produces no
significant change at this \(\beta\).
Removing the 5% highest-closeness nodes cuts the outbreak by 107
infections, twice the random-removal effect. A two-part ML model
(logistic AUC 0.97; random-forest \(R^2=0.38\)) chooses blockers that suppress
spread below any single centrality. The same closeness-central set acts
as both the best accelerant (as seeds) and the best brake (as blockers);
this duality replicates on ca-GrQc.
Introduction
The objective of this work is to study the problem of
spreading on a social network with the SIR epidemiological
model, following the framework taught in Lecture 5 (Spreading)
and Lecture 6 (Robustness) and the dynamics workshop.
Network choice. We reuse the network we analysed in
Homework 1 and 2: the Bluetooth proximity layer of the
Copenhagen Networks Study (Sapiezynski et al., 2019). Bluetooth
co-presence is a standard proxy for face-to-face social contact in
computational social science (Sekara & Lehmann, 2014), so it is a
genuine social network. After the Homework 2 aggregation it is a single
connected component with 662 nodes and 9,251 edges.
This satisfies the brief’s \(>5{,}000\) edges. It does
not reach the \(>1{,}000\) nodes guideline, and no
Copenhagen layer can: the study is a fixed cohort of 845
participants, so no union of layers can exceed 845 nodes (we verified
this on the raw files). Rather than abandon the network we have analysed
all year, we keep Copenhagen as the primary network
and, in Section 8, re-run the entire analysis on a second,
larger social network (the arXiv General Relativity
collaboration network, \(>5{,}000\)
nodes). This covers the size guideline and, more usefully, lets us test
whether our spreading conclusions generalise across two
independent social networks.
Use of Artificial Intelligence. AI tools were used
to assist with coding and editing. All methodological choices,
derivations, and interpretations are our own.
The SIR engine
We implement the SIR model exactly as defined in Lecture 5 (slides
40–41) and the dynamics workshop. The model has three states,
Susceptible, Infected,
Recovered, and two parameters: \(\beta\) (infection probability from an
infected to a susceptible neighbour, per time step) and \(\mu\) (recovery probability of an infected
node, per time step). At each time step every infected node infects each
susceptible neighbour with probability \(\beta\) and recovers with probability \(\mu\).
# run_sir(): SIR engine (Workshop "Spreading"; Lecture 5 slide 41).
# Within a step: infected nodes first RECOVER (prob mu), then those still
# infected INFECT susceptible contacts (prob beta). Returns BOTH per-step
# compartment counts (for the prevalence curve I(t), Lecture 5 slide 43)
# and per-node infection times (used in Section 7).
# NOTE: uses integer indices 1:vcount(g); Copenhagen vertex names are
# CNS ids 0..845, so `seeds` must be indices, never V(g)$name values.
run_sir <- function(g, beta, mu, seeds) {
n <- vcount(g)
state <- integer(n) # 0 = S, 1 = I, 2 = R
state[seeds] <- 1L
infections <- data.frame(t = 0L, inf = seeds) # when each node was infected
counts <- data.frame(t = 0L,
S = sum(state == 0L),
I = sum(state == 1L),
R = sum(state == 2L))
t <- 0L
while (any(state == 1L)) {
t <- t + 1L
# I -> R : each infected recovers with probability mu
inf_idx <- which(state == 1L)
state[inf_idx] <- ifelse(runif(length(inf_idx)) < mu, 2L, 1L)
# S -> I : nodes still infected expose their susceptible neighbours
inf_idx <- which(state == 1L)
sus <- which(state == 0L)
contacts <- as.numeric(unlist(adjacent_vertices(g, inf_idx)))
contacts <- contacts[contacts %in% sus]
new_inf <- unique(contacts[runif(length(contacts)) < beta])
if (length(new_inf) > 0) {
state[new_inf] <- 1L
infections <- rbind(infections, data.frame(t = t, inf = new_inf))
}
counts <- rbind(counts, data.frame(t = t,
S = sum(state == 0L),
I = sum(state == 1L),
R = sum(state == 2L)))
}
list(infections = infections, counts = counts)
}
# Convenience extractors -----------------------------------------------------
# Final epidemic size = everyone who was ever infected. At the end all I -> R,
# so this is the last R count (== nrow(infections)).
final_size <- function(sir) tail(sir$counts$R, 1)
# "Peak of infection" = peak PREVALENCE: the time step where the number of
# simultaneously-infected I(t) is maximal (Lecture 5 slide 43).
peak_time <- function(sir) sir$counts$t[which.max(sir$counts$I)]
peak_prev <- function(sir) max(sir$counts$I)
# 1% seed count. The course idiom sample(1:V, V*0.01) truncates 6.62 -> 6;
# we fix ONE constant and reuse it everywhere for consistency.
n_seed <- function(g, frac = 0.01) max(1L, floor(frac * vcount(g)))
Why the engine returns compartment counts. The
original course sim_sir returns only infection
times. The “peak of infection” is defined on the
prevalence curve \(I(t)\) (Lecture 5 slide 43; the “flatten
the curve” framing of slide 6), which requires the recovery process.
Returning counts per step gives \(I(t)\) directly, while
infections still provides the per-node time-to-infection
used in Section 7.
Network-agnostic design
Every analysis below (Sections 1–7) is written as a function
of an igraph object. This is what makes the Section 8
generality claim credible: the identical code runs on both
networks, so any difference in results is a property of the networks,
not of the implementation.
# One place that turns a graph into the metrics the brief and lectures use.
network_summary <- function(g, label) {
deg <- degree(g); comp <- components(g)
tibble(
Network = label,
`Nodes |V|` = comma(vcount(g)),
`Edges |E|` = comma(ecount(g)),
Density = sprintf("%.4f", edge_density(g)),
`Mean degree <k>` = sprintf("%.2f", mean(deg)),
`<k^2>` = sprintf("%.1f", mean(deg^2)),
`Max degree` = max(deg),
Components = comp$no,
`Largest comp.` = sprintf("%s (%.1f%%)", comma(max(comp$csize)),
100 * max(comp$csize) / vcount(g)),
`Transitivity` = sprintf("%.3f", transitivity(g, "global"))
)
}
# Giant component, exactly the workshop's get_gc() (Lecture 6 / Robustness).
get_gc <- function(g) {
cc <- components(g)
induced_subgraph(g, which(cc$membership == which.max(cc$csize)))
}
mu_global <- 0.1 # recovery probability, the course default (Lecture 5)
0. The networks
0.1 Primary network: Copenhagen Bluetooth proximity
g_cph <- readRDS("bt_network.rds")
stopifnot(vcount(g_cph) == 662, ecount(g_cph) == 9251,
components(g_cph)$no == 1) # fail loudly if the data drifts
g_cph
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
Primary network: the Homework 2 Copenhagen Bluetooth proximity
graph.
| Nodes |V| |
662 |
| Edges |E| |
9,251 |
| Density |
0.0423 |
| Mean degree |
27.95 |
| <k^2> |
1083.9 |
| Max degree |
98 |
| Components |
1 |
| Largest comp. |
662 (100.0%) |
| Transitivity |
0.347 |
The network is a single connected component (good: SIR needs a path
for the infection to travel), it is strongly clustered (transitivity
\(\approx 0.35\)), and its degree
distribution is right-skewed (\(\langle
k^2\rangle \gg \langle k\rangle^2\)). That degree heterogeneity
is exactly what drives a small epidemic threshold, which we make precise
in Section 1.
0.2 Network-choice justification
Compliance note. The brief asks us to reuse the
Homework 1/2 network if it is social, and to email the instructor
only if we change it. The Copenhagen Bluetooth network is a
genuine social network (face-to-face proxy; Sekara & Lehmann, 2014),
is the network from our previous homeworks, and clears the \(>5{,}000\)-edge requirement (9,251
edges). It has 662 nodes, below the \(>1{,}000\) guideline; this cannot be
fixed within the Copenhagen study, whose cohort is fixed at 845 people
(verified on the raw files: 706 Bluetooth, 800 Facebook, 536 calls, 787
roster; union \(=845\)). We therefore
keep Copenhagen as the primary network and
repeat the full analysis on a larger independent social network in
Section 8 (arXiv ca-GrQc, \(>5{,}000\) nodes), which both satisfies
the size guideline and tests the generality of our
conclusions. This is treated as an addition, not a network change; the
choice was communicated in advance.
0.3 Secondary network: arXiv ca-GrQc (loaded in Section 8)
The secondary network is the arXiv General Relativity and
Quantum Cosmology collaboration network (Leskovec, Kleinberg
& Faloutsos, 2007; SNAP). Nodes are authors; an edge means two
authors co-wrote a paper, so it is a social (collaboration) network. It
is downloaded, built, reduced to its giant component, and its size
verified in the document (never hard-coded) in Section
8.
# Defined here, invoked in Section 8. Downloads once, caches locally,
# verifies counts at knit time. Undirected, simplified, giant component.
load_cagrqc <- function(
url = "https://snap.stanford.edu/data/ca-GrQc.txt.gz",
local = "ca-GrQc.txt.gz") {
if (!file.exists(local)) {
utils::download.file(url, local, mode = "wb", quiet = TRUE)
}
el <- read.table(gzfile(local), header = FALSE, comment.char = "#",
col.names = c("from", "to")) # robust: skips # header
g <- graph_from_data_frame(el, directed = FALSE)
g <- igraph::simplify(g) # SNAP file lists each undirected edge twice;
# ::simplify is explicit to avoid purrr::simplify
# masking if library order is ever changed.
get_gc(g) # connected component, as the workshop does
}
1. Theoretical epidemic threshold \(\beta_c\)
For the information to reach a significant fraction of the network
the epidemic must be supercritical. The control
parameter is the basic reproduction number \(R_0\): the average number of new infections
caused by one infected node in an otherwise susceptible population. When
\(R_0>1\) the prevalence grows
exponentially; when \(R_0<1\) the
outbreak dies out (Lecture 5, slides 43–44).
On a network the relevant quantity is the degree of a node’s
neighbours, not its own. Following the friendship-paradox
argument in Lecture 5 (slides 61–62): a susceptible reached along an
edge is attached to a node of degree \(k'\) with probability \(k' p_{k'}/\langle k\rangle\), and
one of that node’s links has already been used to infect it, leaving
\(k'-1\) onward links. Linearising
the early-time dynamics (\(s_k\approx
1\)) gives \(\dfrac{d\Theta}{dt}\approx\Big(\beta\dfrac{\langle
k^2\rangle-\langle k\rangle}{\langle k\rangle}-\mu\Big)\Theta\),
so growth is exponential iff
\[R_0=\frac{\beta}{\mu}\,\frac{\langle
k^2\rangle-\langle k\rangle}{\langle k\rangle}>1
\qquad\Longrightarrow\qquad
\beta_c=\mu\,\frac{\langle k\rangle}{\langle k^2\rangle-\langle
k\rangle}.\]
This is the heterogeneous-population threshold of Lecture 5 slide 43
(Pastor-Satorras & Vespignani, 2001). The factor \(\langle k^2\rangle-\langle k\rangle\) (the
“degree variance” annotation on slide 42) is what makes hubs so
important: the more right-skewed the degree distribution, the smaller
\(\beta_c\). The well-mixed expression
\(\beta_c=\mu/\langle k\rangle\) (same
slide) ignores this and is shown only for contrast.
deg <- degree(g_cph)
kavg <- mean(deg)
ksq <- mean(deg^2)
mu <- mu_global # 0.1, the course default
beta_c <- mu * kavg / (ksq - kavg) # heterogeneous (slide 43)
beta_c_wm <- mu / kavg # well-mixed contrast (slide 43)
tibble(
Quantity = c("$\\langle k\\rangle$", "$\\langle k^2\\rangle$",
"$\\langle k\\rangle^2$",
"$\\langle k^2\\rangle-\\langle k\\rangle$", "$\\mu$",
"$\\beta_c$ (heterogeneous)", "$\\beta_c$ (well-mixed)"),
Value = c(sprintf("%.2f", kavg), sprintf("%.1f", ksq),
sprintf("%.1f", kavg^2), sprintf("%.1f", ksq - kavg),
sprintf("%.2f", mu),
sprintf("%.5f", beta_c), sprintf("%.5f", beta_c_wm))
) %>% kable(caption = "Epidemic threshold for the Copenhagen Bluetooth network ($\\mu=0.1$).")
Epidemic threshold for the Copenhagen Bluetooth network (\(\mu=0.1\)).
| \(\langle
k\rangle\) |
27.95 |
| \(\langle
k^2\rangle\) |
1083.9 |
| \(\langle
k\rangle^2\) |
781.1 |
| \(\langle k^2\rangle-\langle
k\rangle\) |
1056.0 |
| \(\mu\) |
0.10 |
| \(\beta_c\)
(heterogeneous) |
0.00265 |
| \(\beta_c\)
(well-mixed) |
0.00358 |
For this network \(\beta_c \approx
0.00265\). It is tiny precisely because the degree distribution
is heterogeneous: \(\langle k^2\rangle =
1084\) is far larger than \(\langle
k\rangle^2 = 781\), so the friendship paradox makes the network
very easy to ignite. The heterogeneous threshold is about \(1.4\times\) smaller than the well-mixed
value \(\beta_c^{\text{wm}}=0.0036\):
ignoring degree variance would badly over-estimate how infectious the
information has to be. Equivalently, any \(\beta\) we pick maps to a reproduction
number \(R_0=(\beta/\mu)\,(\langle
k^2\rangle-\langle k\rangle)/\langle k\rangle\), so \(\beta=m\,\beta_c\) gives exactly \(R_0=m\), a fact we use in Section 3 to
choose a \(\beta\) “well above”
threshold.
Modelling note. The Bluetooth graph is weighted
(number of mutual scans), but the course \(\beta_c\) and sim_sir use the
unweighted degree and a single \(\beta\) on every edge, so we follow that
convention throughout and treat edge-weighted transmission as a
limitation/extension.
2. SIR sweep over \(\beta\) with 1%
random seeds
We now seed 1% of the nodes at random and simulate
the SIR model for a range of \(\beta\)
below and above \(\beta_c\), with \(\mu=0.1\) fixed. SIR is stochastic, so for
each \(\beta\) we run
N_REP independent realisations and report the mean
final epidemic size (the number of nodes ever infected,
i.e. the \(R\) compartment at the end).
Near \(\beta_c\) realisations show
large fluctuations (Lecture 5; the critical-regime
behaviour of slide 43): some die out, some explode, so we also show the
10–90% band.
set.seed(42)
# Grid design: dense between 1x and 2x beta_c (where the transition
# happens), sparse out on the plateau.
m_grid <- c(0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 5, 6, 8, 10)
beta_grid <- m_grid * beta_c
N_REP <- 50
ns <- n_seed(g_cph) # 1% = 6 nodes, integer indices
sweep <- map_dfr(beta_grid, function(b) {
fin <- replicate(N_REP, {
seeds <- sample(vcount(g_cph), ns) # random 1%, integer indices
final_size(run_sir(g_cph, b, mu_global, seeds))
})
tibble(beta = b, mult = b / beta_c,
mean_final = mean(fin),
lo = unname(quantile(fin, 0.10)),
hi = unname(quantile(fin, 0.90)))
})
Mean final epidemic size over 50 realisations per \(\beta\) (1% random seeds, \(\mu=0.1\)).
| 0.0000 |
0.00 |
6 |
0.9 |
| 0.0007 |
0.25 |
7 |
1.1 |
| 0.0013 |
0.50 |
11 |
1.6 |
| 0.0020 |
0.75 |
15 |
2.2 |
| 0.0026 |
1.00 |
25 |
3.8 |
| 0.0033 |
1.25 |
41 |
6.2 |
| 0.0040 |
1.50 |
118 |
17.7 |
| 0.0046 |
1.75 |
207 |
31.3 |
| 0.0053 |
2.00 |
272 |
41.1 |
| 0.0066 |
2.50 |
364 |
55.0 |
| 0.0079 |
3.00 |
437 |
66.1 |
| 0.0106 |
4.00 |
497 |
75.1 |
| 0.0132 |
5.00 |
544 |
82.2 |
| 0.0159 |
6.00 |
566 |
85.5 |
| 0.0212 |
8.00 |
596 |
90.0 |
| 0.0265 |
10.00 |
609 |
92.0 |
ggplot(sweep, aes(beta, mean_final)) +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, fill = "#1f78b4") +
geom_line(color = "#1f78b4", linewidth = 0.7) +
geom_point(color = "#1f78b4", size = 1.7) +
geom_vline(xintercept = beta_c, linetype = 2) +
annotate("text", x = beta_c, y = max(sweep$hi) * 0.97,
label = "beta[c]", parse = TRUE, hjust = -0.25, size = 3.6) +
scale_y_continuous(
sec.axis = sec_axis(~ . / vcount(g_cph) * 100, name = "% of network")) +
labs(x = expression(beta),
y = "Mean final epidemic size (nodes ever infected)") +
theme_minimal(base_size = 11)
Reading the transition. Below \(\beta_c\) the system is
subcritical: outbreaks die out almost immediately, the
final size stays at roughly the 6 seeds. At \(\beta=\beta_c\) the mean outbreak is still
small (25 nodes, 3.8% of the network) but the 10–90% band is at its
widest — the critical regime’s large fluctuations. The
mean outbreak first reaches 10% of the network at \(\beta \approx 0.0040\) (\(1.5\times\beta_c\)), and by \(4\beta_c\) (\(R_0\approx4\)) it is firmly
supercritical, infecting on average 497 nodes (75% of
the network).
The empirical takeoff sits at \(1.5\,\beta_c\), not exactly at \(\beta_c\), and the transition is smooth
rather than a sharp step. Both are expected from Section 1: the
heterogeneous-mean-field \(\beta_c\) is
a large-network approximation, and because this network is heterogeneous
but not heavy-tailed (\(\beta_c\) is
only \(1.4\times\) below the well-mixed
value), the transition is broad. Rather than fix an operating point by
hand, Section 3 sweeps this supercritical range and data-chooses
\(\beta^*\) at the value where
a better seed set helps most, and separately characterises the saturated
end (\(\approx 4\beta_c\)), where seed
choice no longer changes the final size.
3. A better 1% of seeds
We must “choose a \(\beta\) well
above \(\beta_c\)”. Section 2 showed
the final size is steeply \(\beta\)-dependent up to \(\sim 3\beta_c\) and then saturates. Whether
a better seed set helps therefore also depends on \(\beta\): near the transition there is room
to recruit more of the network, while deep in the supercritical regime
the epidemic fills the reachable component regardless of where
it starts, so a targeted seed set can only bring the peak
forward without changing its height.
So instead of guessing one \(\beta\), we sweep the seeding
advantage across the supercritical regime (\(1.5\)–\(4\,\beta_c\), i.e. \(R_0=1.5\)–\(4\), all well above threshold), report (a)
the difference in total infected and (b) the difference in
peak-prevalence time as functions of \(\beta\), and then fix \(\beta^*\) at the operating point where a
better seed set is most effective — a data-chosen \(\beta^*\) for Sections 4–7 rather than a
heuristic one.
Lecture 5 (slides 51–52) makes degree the natural
first choice (“high-degree nodes tend to be good”) while warning it is
redundant (hubs adjacent to hubs) and ignores other centrality notions.
We compare random 1% seeds against degree, closeness, betweenness,
eigenvector, PageRank, \(k\)-core, and
a community-spread rule (top-degree nodes distributed
across Louvain communities, to beat the redundancy the lecture warns
about). The full greedy optimum (Kempe, Kleinberg & Tardos, 2003;
slide 52) is \(O(k\,N\,R)\)
simulations; we let the learned seeding of Section 7
play that “smarter than one centrality” role at feasible cost.
set.seed(42)
ns <- n_seed(g_cph) # 1% = 6 seeds, integer indices
NREP3 <- 50
mults3 <- c(1.5, 2, 2.5, 3, 4) # R0 = 1.5 .. 4, all well above 1
# Centralities, computed once on integer-indexed vertices.
cen <- list(
degree = degree(g_cph),
closeness = closeness(g_cph),
betweenness = betweenness(g_cph),
eigenvector = eigen_centrality(g_cph)$vector,
pagerank = page_rank(g_cph)$vector,
kcore = coreness(g_cph)
)
top_by <- function(v, k) order(v, decreasing = TRUE)[seq_len(k)]
# Community-spread: Louvain, allocate seeds across communities by size,
# take the top-degree node(s) inside each --- spreads the seeds out
# instead of clustering them among mutually-adjacent hubs (slide 52).
set.seed(42)
comm <- cluster_louvain(g_cph); mem <- membership(comm)
seeds_community <- function(k) {
tab <- sort(table(mem), decreasing = TRUE)
alloc <- floor(as.numeric(tab) / sum(tab) * k)
while (sum(alloc) < k) { # leftovers to the big communities
j <- which.max((as.numeric(tab) / sum(tab) * k) - alloc); alloc[j] <- alloc[j] + 1
}
unlist(lapply(seq_along(tab), function(i) {
idx <- which(mem == as.integer(names(tab)[i]))
if (alloc[i] == 0) return(integer(0))
idx[order(cen$degree[idx], decreasing = TRUE)][seq_len(alloc[i])]
}))
}
strategies <- list(
random = NULL, # re-sampled every realisation
degree = top_by(cen$degree, ns),
closeness = top_by(cen$closeness, ns),
betweenness = top_by(cen$betweenness, ns),
eigenvector = top_by(cen$eigenvector, ns),
pagerank = top_by(cen$pagerank, ns),
kcore = top_by(cen$kcore, ns),
community = seeds_community(ns)
)
run_strategy <- function(seed_set, beta) {
fin <- pk <- numeric(NREP3)
for (r in seq_len(NREP3)) {
s <- if (is.null(seed_set)) sample(vcount(g_cph), ns) else seed_set
sir <- run_sir(g_cph, beta, mu_global, s)
fin[r] <- final_size(sir); pk[r] <- peak_time(sir)
}
list(fin = fin, pk = pk)
}
# Full grid: 5 beta x 8 strategies x 50 realisations.
set.seed(42)
raw3 <- list()
for (m in mults3) for (nm in names(strategies)) {
raw3[[paste(m, nm)]] <- c(list(mult = m, strategy = nm),
run_strategy(strategies[[nm]], m * beta_c))
}
ci <- function(x) qt(0.975, length(x) - 1) * sd(x) / sqrt(length(x))
sum3 <- map_dfr(raw3, function(o) tibble(
mult = o$mult, strategy = o$strategy,
fin_m = mean(o$fin), fin_ci = ci(o$fin),
pk_m = mean(o$pk), pk_ci = ci(o$pk)))
# Difference of the BEST targeted strategy vs random, at each beta.
delta_best <- map_dfr(mults3, function(m) {
rr <- raw3[[paste(m, "random")]]$fin
rpk <- raw3[[paste(m, "random")]]$pk
cand <- setdiff(names(strategies), "random")
best <- cand[ which.max(vapply(cand,
function(nm) mean(raw3[[paste(m, nm)]]$fin), numeric(1))) ]
bf <- raw3[[paste(m, best)]]$fin; bpk <- raw3[[paste(m, best)]]$pk
tt <- t.test(bf, rr)
tibble(mult = m, beta = m * beta_c, best = best,
d_fin = mean(bf) - mean(rr),
d_lo = tt$conf.int[1], d_hi = tt$conf.int[2],
d_p = tt$p.value,
d_pk = mean(bpk) - mean(rpk))
})
# Data-chosen operating point: the beta where a better seed set helps
# total infections the most. Reused as beta_star in Sections 4-7.
star <- delta_best %>% slice_max(d_fin, n = 1, with_ties = FALSE)
beta_star <- star$beta
mult_star <- star$mult
primary_strategy <- star$best
seeds_primary <- strategies[[primary_strategy]]
p_size <- ggplot(delta_best, aes(mult, d_fin)) +
geom_hline(yintercept = 0, color = "grey70") +
geom_ribbon(aes(ymin = d_lo, ymax = d_hi), alpha = 0.15, fill = "#1f78b4") +
geom_line(color = "#1f78b4") + geom_point(color = "#1f78b4", size = 2) +
geom_vline(xintercept = mult_star, linetype = 2) +
labs(x = expression(beta/beta[c]),
y = "Extra infected vs random") +
theme_minimal(base_size = 11)
p_time <- ggplot(delta_best, aes(mult, d_pk)) +
geom_hline(yintercept = 0, color = "grey70") +
geom_line(color = "#d95f02") + geom_point(color = "#d95f02", size = 2) +
geom_vline(xintercept = mult_star, linetype = 2) +
labs(x = expression(beta/beta[c]),
y = "Peak-time shift (steps)") +
theme_minimal(base_size = 11)
gridExtra::grid.arrange(p_size, p_time, ncol = 2)
At the data-chosen operating point \(\beta^* = 0.0053\) (\(2.0\,\beta_c\), \(R_0\approx2.0\)) we report the full
ranking; this \(\beta^*\) is reused in
Sections 4–7.
Seeding strategies at the chosen \(\beta^*=2.0\beta_c\) (6 seeds = 1%, 50
realisations). Δ = targeted minus random.
| closeness |
313 |
[302, 324] |
+74 |
44.6 |
-4.3 |
| degree |
310 |
[299, 322] |
+70 |
37.6 |
-11.3 |
| betweenness |
295 |
[280, 309] |
+55 |
37.6 |
-11.3 |
| community |
292 |
[277, 306] |
+52 |
35.1 |
-13.8 |
| kcore |
287 |
[268, 307] |
+47 |
41.4 |
-7.5 |
| pagerank |
274 |
[251, 297] |
+34 |
41.5 |
-7.4 |
| eigenvector |
252 |
[218, 286] |
+12 |
42.6 |
-6.3 |
| random |
240 |
[204, 275] |
— |
48.9 |
— |
set.seed(7)
cur_rnd <- run_sir(g_cph, beta_star, mu_global,
sample(vcount(g_cph), ns))$counts
cur_best <- run_sir(g_cph, beta_star, mu_global, seeds_primary)$counts
bind_rows(
cur_rnd %>% transmute(t, I, which = "Random"),
cur_best %>% transmute(t, I, which = sprintf("Targeted (%s)", primary_strategy))
) %>%
ggplot(aes(t, I, color = which)) +
geom_line(linewidth = 0.8) +
scale_color_manual(values = c("grey55", "#d95f02"), name = NULL) +
labs(x = "Time step", y = "Infected (prevalence)") +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom")
Result. The two effects the brief asks for behave
very differently. (a) Total infected: a better seed set
helps most near the transition and the gain shrinks toward zero as
\(\beta\) grows: at the saturated
end (\(4\beta_c\)) the best strategy
adds only +21 nodes, while at the chosen \(\beta^*=2.0\beta_c\) it adds +74 (best
strategy: closeness). (b) Peak timing:
targeted seeding brings the peak forward by -4.3 steps at \(\beta^*\) and does so consistently
across the whole supercritical range (right panel).
This reconciles the two parts of the brief: extra total
infections are only achievable near the transition; deep in the
supercritical regime the epidemic percolates the giant component
whatever the seeds, and the only remaining leverage is
timing, which our \(\beta\)-resolved sweep recovers as the
saturation limit of a single curve. The primary strategy is therefore
closeness, the size-winner at the data-chosen \(\beta^*\); it is what Sections 5–7
reuse.
4. Quarantine strategy
Using the same \(\beta^*=2.0\beta_c\), we add a temporary
quarantine during the exponential-growth phase. We model quarantine by
temporarily removing susceptible nodes, but
not with delete_vertices: the course engine is
indexed by \(1{:}|V|\), and deleting
vertices renumbers them, which would silently corrupt the
seed/centrality indices. Instead the engine carries a fourth state,
\(Q\) (a state mask): quarantined nodes
cannot be infected and cannot transmit, and at release they return to
S with their identity intact.
# Course SIR + quarantine. state: 0=S 1=I 2=R 3=Q. At t_q a fraction
# `frac` of the CURRENT susceptibles is moved S->Q; `release` steps
# later Q->S (susceptible again). who = random | degree (which 20%).
run_sir_quar <- function(g, beta, mu, seeds, t_q, frac = 0.20,
release = 8L, who = c("random", "degree")) {
who <- match.arg(who)
n <- vcount(g); state <- integer(n); state[seeds] <- 1L
deg_g <- degree(g)
t_rel <- NA_integer_
rec <- function(tt) data.frame(t = tt, S = sum(state == 0L),
I = sum(state == 1L), R = sum(state == 2L), Q = sum(state == 3L))
counts <- rec(0L); t <- 0L
while (any(state == 1L)) {
t <- t + 1L
if (!is.na(t_rel) && t == t_rel) state[state == 3L] <- 0L # Q -> S
if (t == t_q) { # S -> Q
sus <- which(state == 0L); k <- floor(frac * length(sus))
if (k > 0) {
pick <- if (who == "random") sample(sus, k)
else sus[order(deg_g[sus], decreasing = TRUE)][seq_len(k)]
state[pick] <- 3L; t_rel <- t + release
}
}
inf_idx <- which(state == 1L) # I -> R
state[inf_idx] <- ifelse(runif(length(inf_idx)) < mu, 2L, 1L)
inf_idx <- which(state == 1L); sus <- which(state == 0L) # S -> I
contacts <- as.numeric(unlist(adjacent_vertices(g, inf_idx)))
contacts <- contacts[contacts %in% sus]
new_inf <- unique(contacts[runif(length(contacts)) < beta])
if (length(new_inf) > 0) state[new_inf] <- 1L
counts <- rbind(counts, rec(t))
}
counts
}
# Mean prevalence curve I(t) over `R` reps, padded to a common length
# (I = 0 after an epidemic ends), for no-quarantine vs the two policies.
mean_I <- function(run_one, R, Tmax = 160L) {
M <- matrix(0, nrow = R, ncol = Tmax + 1L)
for (r in seq_len(R)) {
cn <- run_one()
ix <- pmin(cn$t, Tmax) + 1L
M[r, ix] <- cn$I
}
colMeans(M)
}
set.seed(42)
NREP4 <- 50
ns <- n_seed(g_cph) # random 1% seeds = uncontrolled epidemic
# Mean uncontrolled I(t) at beta_star, to locate the exponential phase.
baseI <- mean_I(function() run_sir(g_cph, beta_star, mu_global,
sample(vcount(g_cph), ns))$counts, NREP4)
peakI_t <- which.max(baseI) - 1L # time index (t starts at 0)
# Exponential-phase onset = the EPIDEMIC-DOUBLING time: first step at
# which mean prevalence has doubled relative to the seeds. This is the
# standard descriptor of exponential growth, is robust to the seed
# level (a "% of peak" rule fails when seeds already exceed it), and
# scales with the epidemic's own speed.
t_q <- max(2L, which(baseI >= 2 * baseI[1])[1] - 1L)
The brief suggests \(t\approx 3\)
“but check”. At \(\beta^*=2\beta_c\)
the mean uncontrolled peak sits at \(t\approx
49\), so the exponential phase begins much later in absolute time
than \(t=3\). We quarantine at the
detected epidemic-doubling time \(t_q=13\) (the first step at which mean
prevalence has doubled relative to the seeds: the standard onset of
exponential growth), holding for 8 steps. A fixed \(t\approx3\) would intervene before the
epidemic has even entered exponential growth here.
set.seed(42)
# Paired design: identical seed set per replicate across the three
# scenarios, so each difference is within-replicate. This removes the
# enormous between-run variance of this near-critical regime
# (Section 2's "large fluctuations"), which a bare mean would hide.
pair <- map_dfr(seq_len(NREP4), function(r) {
s <- sample(vcount(g_cph), ns)
c0 <- run_sir(g_cph, beta_star, mu_global, s)$counts
cR <- run_sir_quar(g_cph, beta_star, mu_global, s, t_q, who = "random")
cD <- run_sir_quar(g_cph, beta_star, mu_global, s, t_q, who = "degree")
tibble(rep = r,
f0 = tail(c0$R, 1), fR = tail(cR$R, 1), fD = tail(cD$R, 1),
p0 = max(c0$I), pR = max(cR$I), pD = max(cD$I),
k0 = c0$t[which.max(c0$I)], kR = cR$t[which.max(cR$I)],
kD = cD$t[which.max(cD$I)])
})
# Paired difference vs no quarantine: mean, 95% CI, paired-t p-value.
pt <- function(x, y) { tt <- t.test(x - y)
list(m = mean(x - y), lo = tt$conf.int[1], hi = tt$conf.int[2],
p = tt$p.value, sig = tt$p.value < 0.05) }
dF_R <- pt(pair$fR, pair$f0); dF_D <- pt(pair$fD, pair$f0)
dP_R <- pt(pair$pR, pair$p0); dP_D <- pt(pair$pD, pair$p0)
dT_R <- pt(pair$kR, pair$k0); dT_D <- pt(pair$kD, pair$k0)
base_f <- mean(pair$f0); base_p <- mean(pair$p0); base_k <- mean(pair$k0)
Paired effect of an 8-step quarantine of 20% of susceptibles
started at the epidemic-doubling time \(t_q=13\), \(\beta^*=2.0\beta_c\), 50 paired
realisations. Δ = scenario minus no-quarantine on the same seeds; ‘n.s.’
= 95% CI includes 0.
| Final size (nodes ever infected) |
236.6 |
+15.9 [-28.2, +60.1] n.s. |
+20.1 [-24.1, +64.3] n.s. |
| Peak height (max simultaneous I) |
45.2 |
-1.7 [-9.7, +6.3] n.s. |
+0.3 [-8.0, +8.5] n.s. |
| Peak time (step of max I) |
44.7 |
+7.4 [-5.1, +20.0] n.s. |
+8.7 [-2.9, +20.2] n.s. |
curR <- mean_I(function() run_sir_quar(g_cph, beta_star, mu_global,
sample(vcount(g_cph), ns), t_q, who = "random"), NREP4)
curD <- mean_I(function() run_sir_quar(g_cph, beta_star, mu_global,
sample(vcount(g_cph), ns), t_q, who = "degree"), NREP4)
tibble(t = 0:160,
`No quarantine` = baseI, `20% random` = curR, `20% by degree` = curD) %>%
pivot_longer(-t, names_to = "Scenario", values_to = "I") %>%
filter(t <= max(which(baseI > 0.5)) + 10) %>%
ggplot(aes(t, I, color = Scenario)) +
annotate("rect", xmin = t_q, xmax = t_q + 8, ymin = -Inf, ymax = Inf,
alpha = 0.10, fill = "grey40") +
geom_line(linewidth = 0.8) +
scale_color_manual(values = c("grey55", "#1f78b4", "#d95f02")) +
labs(x = "Time step", y = "Infected (prevalence)", color = NULL) +
theme_minimal(base_size = 11) + theme(legend.position = "bottom")
Result. With paired replicates and 95% confidence
intervals the conclusion is clear. At \(\beta^*=2\beta_c\) this network is
near-critical (Section 2’s large fluctuations), and an 8-step
reversible quarantine of 20% of susceptibles produces no
statistically significant change in any outcome: final size
(random +15.9 (95% CI [-28.2, +60.1], n.s.), by-degree +20.1 (95% CI
[-24.1, +64.3], n.s.)), peak height (random -1.7 (95% CI [-9.7, +6.3],
n.s.), by-degree +0.3 (95% CI [-8.0, +8.5], n.s.)), and even peak timing
(random +7.4 steps (95% CI [-5.1, +20.0], n.s.), by-degree +8.7 steps
(95% CI [-2.9, +20.2], n.s.)) all have 95% CIs that include zero. The
point estimates are mechanistically sensible: a mild peak
delay of roughly the 8-step quarantine duration, with
by-degree directionally stronger than random (the hub leverage of
Sections 3 and 6). None of these signals survives the intrinsic
run-to-run variability of a near-critical epidemic.
That null is the answer to “measure the difference with
respect to no quarantine”: a brief, reversible, partial
quarantine is statistically swamped by the fluctuations of a
near-critical outbreak. Truly “flattening the curve” (Lecture 5)
requires a sustained or larger intervention, not a short pulse
lifted while transmission is still ongoing; this is why real
non-pharmaceutical interventions are held in place rather than released
early. At a faster, well-supercritical \(\beta\) the same pulse bites harder: the
\(\beta\)-dependence mirrors the
seeding result of Section 3.
Which 20% to quarantine (random vs by-degree).
Resolved by the data: at \(\beta^*=2\beta_c\) neither random
nor by-degree quarantine produces a statistically significant change in
final size, peak height, or peak time: every 95% CI includes
zero, and the random-vs-by-degree difference is itself within noise. We
report the null directly, with CIs and the mechanistically sensible but
non-significant point estimates, rather than picking a “winner” from
noise.
5. Convincing 5% not to spread
A person convinced “not to spread at all” is, for the dynamics,
exactly an immunised/removed node: Lecture 5 (slides
47–50) gives the equivalence “vaccinating a node so it stops
spreading … is equivalent to removing it”, and the
immunisation-of-complex-networks literature builds on the same identity
(Pastor-Satorras & Vespignani, 2002). The brief itself seeds “1% of
the remaining nodes”, which only makes sense if the 5%
are taken out. So we delete the 5%, then run the SIR at the
same \(\beta^*=2\beta_c\) on the reduced network,
seeding 1% of the survivors at random.
The 5% are chosen two ways: (A) at random, and
(B) by closeness, the same centrality the data
picked as the Section 3 primary, so Section 6 can relate the two roles
of the identical metric. We measure the A-vs-B difference exactly as in
Section 3 (total infected and peak-prevalence time), and add the
Lecture-5 immunisation theory: the random critical fraction \(g_c\approx 1-\frac{\mu}{\beta}\frac{\langle
k\rangle}{\langle k^2\rangle}\) (slide 48) and the reproduction
number \(R_0\) of each reduced
graph.
set.seed(42)
NREP5 <- 50
nblock <- round(0.05 * vcount(g_cph)) # 5% non-spreaders
clo <- closeness(g_cph) # = the Section 3 primary
# Removal RE-INDEXES vertices to 1:|V'|, so seeds are sampled on the
# *subgraph*. We never map back to names: §5 is purely index-based.
make_reduced <- function(rm)
induced_subgraph(g_cph, setdiff(seq_len(vcount(g_cph)), rm))
set.seed(42)
gr_none <- g_cph
gr_rand <- make_reduced(sample(vcount(g_cph), nblock))
gr_clo <- make_reduced(order(clo, decreasing = TRUE)[seq_len(nblock)])
run_block <- function(gr) {
nseed_r <- max(1L, floor(0.01 * vcount(gr))) # 1% of the REMAINING
fin <- pk <- numeric(NREP5)
for (r in seq_len(NREP5)) {
s <- sample(vcount(gr), nseed_r) # random seeds, subgraph idx
sir <- run_sir(gr, beta_star, mu_global, s)
fin[r] <- final_size(sir); pk[r] <- peak_time(sir)
}
list(fin = fin, pk = pk)
}
set.seed(42)
B <- list(`No blockers` = run_block(gr_none),
`5% random` = run_block(gr_rand),
`5% closeness` = run_block(gr_clo))
R0_of <- function(g) { d <- degree(g) # Lecture 5 slide 42
(beta_star / mu_global) * (mean(d^2) - mean(d)) / mean(d) }
gc_random <- 1 - (mu_global / beta_star) *
(mean(degree(g_cph)) / mean(degree(g_cph)^2)) # slide 48
ci_ <- function(x) qt(0.975, length(x) - 1) * sd(x) / sqrt(length(x))
bsum <- tibble(
Scenario = names(B),
fin_m = sapply(B, function(o) mean(o$fin)),
fin_ci = sapply(B, function(o) ci_(o$fin)),
pk_m = sapply(B, function(o) mean(o$pk)),
R0 = c(R0_of(gr_none), R0_of(gr_rand), R0_of(gr_clo)))
base_fin <- bsum$fin_m[bsum$Scenario == "No blockers"]
# A-vs-B difference, exactly as Section 3 (closeness blockers - random)
dAB_fin <- t.test(B$`5% closeness`$fin, B$`5% random`$fin)
dAB_pk <- t.test(B$`5% closeness`$pk, B$`5% random`$pk)
5% non-spreaders removed, then SIR at \(\beta^*=2\beta_c\) with 1% random seeds
among the survivors (50 realisations). \(R_0\) is the reproduction number of the
reduced graph.
| No blockers |
271 |
[240, 303] |
— |
47.5 |
2.00 |
| 5% random |
220 |
[186, 253] |
-52 |
50.0 |
1.93 |
| 5% closeness |
164 |
[139, 190] |
-107 |
56.9 |
1.64 |
set.seed(7)
cv <- function(gr) mean_I(function() {
run_sir(gr, beta_star, mu_global,
sample(vcount(gr), max(1L, floor(0.01 * vcount(gr)))))$counts }, 40L)
tibble(t = 0:160, `No blockers` = cv(gr_none),
`5% random` = cv(gr_rand), `5% closeness` = cv(gr_clo)) %>%
pivot_longer(-t, names_to = "Scenario", values_to = "I") %>%
filter(t <= 90) %>%
ggplot(aes(t, I, color = Scenario)) + geom_line(linewidth = 0.8) +
scale_color_manual(values = c("grey55", "#1f78b4", "#d95f02")) +
labs(x = "Time step", y = "Infected (prevalence)", color = NULL) +
theme_minimal(base_size = 11) + theme(legend.position = "bottom")
Result. The brief’s A-vs-B comparison: choosing the
5% by closeness rather than at random suppresses the
total infected significantly more, -56 (95% CI [-97, -14],
p<0.05) (closeness minus random); the peak-time difference is smaller
and not significant (+7 (95% CI [-9, +22], n.s.)). Both
removals have a real, partial effect: 5% random cuts the outbreak by -52
(reduced \(R_0\) from 2.00 to 1.93),
and 5% by-closeness by roughly twice as much, -107 (\(R_0\) to 1.64). Neither removal eliminates
the epidemic, and that is exactly slide 48: the random critical fraction
here is \(g_c\approx0.51\), so 5%
(\(\ll g_c\)) cannot drive \(R_0<1\) — even the targeted 5% only
lowers it to \(1.64>1\). The slide
49–50 lesson is efficiency, not elimination: per node
removed, taking the highest-closeness 5% is about twice as effective as
random at shrinking the outbreak (and delays its peak), but stopping the
epidemic would still need a far larger fraction. We carry this A-vs-B
contrast into Section 6.
6. Relating Sections 3 and 5
The brief asks us to relate Sections 3 and 5 using the same type
of centrality. We did exactly that by construction: the data picked
closeness as the Section 3 primary, and Section 5
removed the 5% highest-closeness nodes. So the
identical structural quantity plays two opposite roles, and no
new simulation is needed; the two effects are placed side by side
below.
mr <- mean(degree(g_cph)^2) / mean(degree(g_cph)) # Molloy-Reed ratio (s.25)
fc <- 1 - 1 / (mr - 1) # random-failure f_c (s.28)
# Closeness as 1% SEEDS (Section 3): extra infected vs random seeds.
s3 <- res3 %>% filter(strategy %in% c("closeness", "random"))
seed_gain <- s3$fin_m[s3$strategy == "closeness"] -
s3$fin_m[s3$strategy == "random"]
# Closeness as 5% BLOCKERS (Section 5): change vs random blockers.
blk_drop <- bsum$fin_m[bsum$Scenario == "5% closeness"] -
bsum$fin_m[bsum$Scenario == "5% random"]
n_s <- n_seed(g_cph); n_b <- round(0.05 * vcount(g_cph))
dual <- tibble(
`Role of the closeness-central nodes` =
c(sprintf("As 1%% seeds (%d nodes), Section 3", n_s),
sprintf("As 5%% blockers (%d nodes), Section 5", n_b)),
`Effect vs the random choice` =
c(sprintf("%+.0f infected", seed_gain),
sprintf("%+.0f infected", blk_drop)),
`Per node` =
c(sprintf("%+.1f / seed", seed_gain / n_s),
sprintf("%+.1f / blocker", blk_drop / n_b)))
The same closeness centrality as accelerant and as brake (vs
the random choice, at \(\beta^*=2\beta_c\)). Both effects are
statistically significant (Sections 3 and 5).
| As 1% seeds (6 nodes), Section 3 |
+74 infected |
+12.3 / seed |
| As 5% blockers (33 nodes), Section 5 |
-56 infected |
-1.7 / blocker |
tibble(role = c("Closeness seeds (S3)", "Closeness blockers (S5)"),
d = c(seed_gain, blk_drop)) %>%
ggplot(aes(reorder(role, d), d, fill = d > 0)) +
geom_col(width = 0.55) + geom_hline(yintercept = 0, color = "grey40") +
coord_flip() +
scale_fill_manual(values = c("#1f78b4", "#d95f02"), guide = "none") +
labs(x = NULL, y = "Extra infected vs random choice") +
theme_minimal(base_size = 11)
Relationship. The same closeness-central nodes act
as both accelerant (as 1% seeds, they add +74 infected over random seeds
in Section 3) and brake (as 5% blockers, they remove -56 relative to
random blockers in Section 5). The set of nodes that best accelerates a
cascade when used as seeds is also the set that best suppresses it when
used as blockers.
This is one phenomenon, not two, and Lecture 6 names its cause. The
Molloy–Reed ratio here is \(\langle
k^2\rangle/\langle k\rangle \approx 39
\gg 2\), so the giant component is deeply robust to
random loss (random-failure threshold \(f_c \approx 0.97\) — about 97% of nodes
would have to fail at random to fragment it). That is precisely why
Section 5’s random 5% blockers, and the random epidemic
threshold \(g_c\approx0.51\), are so
weak: 5% is far below either bar. But heterogeneous networks are
fragile to targeted attack (slides 29–30; Albert, Jeong &
Barabási, 2000): removing the highest-closeness hubs collapses \(\langle k^2\rangle\), which is what made
Section 5’s closeness blockers roughly twice as effective. The
common root is the same heavy second moment \(\langle k^2\rangle \gg \langle k\rangle^2\)
that gave Section 1 its tiny \(\beta_c\) and the friendship paradox its
bite: one structural fact ties together the epidemic threshold (S1),
seed leverage (S3), targeted-blocker leverage (S5), and the
robust-to-random / fragile-to-targeted asymmetry (Lecture 6).
7. Predicting time-to-infection and seeding from it
We now learn who gets infected early from structure alone,
then use that model to choose seeds and blockers. This is the
social-network sensors idea of Christakis & Fowler (2010),
recast as a predictive model. Note this is not
circular: the predictors are purely structural (known before
any epidemic), the label is the emergent time-to-infection from the
Section 3 random-seed runs, and every evaluation below uses
fresh simulations with the model-chosen sets.
set.seed(42)
R_TRAIN <- 60L
NV <- vcount(g_cph) # NOT `V` (would shadow igraph::V)
ns_ml <- n_seed(g_cph)
n_inf <- n_ns <- sum_t <- numeric(NV)
for (r in seq_len(R_TRAIN)) {
sd_ <- sample(NV, ns_ml)
inf <- run_sir(g_cph, beta_star, mu_global, sd_)$infections
firstt <- tapply(inf$t, inf$inf, min) # earliest t per node
tvec <- rep(NA_real_, NV)
tvec[as.integer(names(firstt))] <- as.numeric(firstt)
is_seed <- logical(NV); is_seed[sd_] <- TRUE
use <- !is_seed # exclude reps where seed
n_ns[use] <- n_ns[use] + 1L
got <- use & !is.na(tvec) & tvec > 0
n_inf[got] <- n_inf[got] + 1L
sum_t[got] <- sum_t[got] + tvec[got]
}
p_inf <- ifelse(n_ns > 0, n_inf / n_ns, NA_real_)
t_inf <- ifelse(n_inf > 0, sum_t / n_inf, NA_real_)
lc <- transitivity(g_cph, type = "local", isolates = "zero"); lc[is.na(lc)] <- 0
cm <- membership(cluster_louvain(g_cph))
feat <- data.frame(
degree = degree(g_cph), strength = strength(g_cph),
closeness = closeness(g_cph), between = betweenness(g_cph),
eigen = eigen_centrality(g_cph)$vector,
pagerank = page_rank(g_cph)$vector, kcore = coreness(g_cph),
clustering = lc,
comm_size = as.integer(table(cm)[as.character(cm)]))
# gender joined by CNS id = V(g)$name (guarded: NOT by 1:NV position)
gd <- read.table("genders.csv", sep = ",", header = FALSE,
comment.char = "#", col.names = c("user", "female"))
gid <- suppressWarnings(as.integer(V(g_cph)$name))
fem <- gd$female[match(gid, gd$user)]
feat$female <- factor(ifelse(is.na(fem), "U", as.character(fem)))
gender_cov <- mean(!is.na(fem))
feat$t_inf <- t_inf
# Censoring: at beta*=2beta_c the epidemic reaches only a minority,
# so many nodes are seldom infected and t_inf is undefined/noisy for
# them. A single regression would have to drop them or pollute the
# target with a sentinel. We use a TWO-PART model:
# (A) logistic for whether a node is reliably reached (p_inf >= 0.5),
# (B) RF regression of t_inf on the reliably-reached "core".
never <- mean(p_inf == 0, na.rm = TRUE)
core <- !is.na(p_inf) & p_inf >= 0.5
feat$core <- factor(ifelse(core, "yes", "no"), levels = c("no", "yes"))
preds <- c("degree","strength","closeness","between","eigen",
"pagerank","kcore","clustering","comm_size","female")
ctrlC <- trainControl(method = "cv", number = 10, classProbs = TRUE,
summaryFunction = twoClassSummary)
ctrlR <- trainControl(method = "cv", number = 10)
set.seed(42)
fitA <- train(reformulate(preds, "core"), data = feat, method = "glm",
family = binomial(), metric = "ROC", trControl = ctrlC)
dcore <- feat[core, c(preds, "t_inf")]
set.seed(42)
fitB <- train(reformulate(preds, "t_inf"), data = dcore, method = "rf",
metric = "RMSE", trControl = ctrlR, ntree = 300,
tuneGrid = data.frame(mtry = c(3, 5, 7)), importance = TRUE)
set.seed(42)
fitL <- train(reformulate(preds, "t_inf"), data = dcore, method = "lm",
metric = "RMSE", trControl = ctrlR)
aucA <- max(fitA$results$ROC)
rfR2 <- max(fitB$results$Rsquared, na.rm = TRUE)
rfRMSE <- min(fitB$results$RMSE)
lmRMSE <- min(fitL$results$RMSE)
vi <- varImp(fitB)$importance
vtop <- rownames(vi)[order(-vi[, 1])][1:3]
predA <- predict(fitA, feat, type = "prob")[, "yes"]
predT <- predict(fitB, feat)
cand <- which(predA >= 0.5)
need <- round(0.05 * NV) # most we will ever draw (blockers)
rankN <- if (length(cand) >= need) cand[order(predT[cand])] else order(predT)
# predicted earliest-infected (= most influential), robust if 'core' is small
Two-part time-to-infection model, 10-fold CV (Homework-2
style). 3% of nodes are never infected across the 60 random-seed runs
and 35% are ‘core’ (p_inf>=0.5); gender coverage 96%. Top predictors:
degree, comm_size, eigen.
| (A) Logistic: node reliably reached? |
AUC = 0.975 |
| (B) RF: time-to-infection | reached |
R2 = 0.38, RMSE = 4.52 |
| (B’) Linear baseline |
RMSE = 4.84 |
set.seed(42)
NREP7 <- 40L
clo_rank <- order(closeness(g_cph), decreasing = TRUE)
ml_seeds <- rankN[seq_len(ns_ml)]
arm <- function(ss) replicate(NREP7, final_size(run_sir(
g_cph, beta_star, mu_global,
if (is.null(ss)) sample(NV, ns_ml) else ss)))
S <- list(Random = arm(NULL),
`Closeness (S3)` = arm(clo_rank[seq_len(ns_ml)]),
`ML model` = arm(ml_seeds))
# Repeat Section 5 with this knowledge: 5% blockers = ML top influentials
nb <- round(0.05 * NV)
blk <- function(rm) { gr <- induced_subgraph(g_cph, setdiff(seq_len(NV), rm))
nse <- max(1L, floor(0.01 * vcount(gr)))
replicate(NREP7, final_size(run_sir(gr, beta_star, mu_global,
sample(vcount(gr), nse)))) }
set.seed(42)
Bk <- list(`Random 5%` = blk(sample(NV, nb)),
`Closeness 5% (S5)`= blk(clo_rank[seq_len(nb)]),
`ML 5%` = blk(rankN[seq_len(nb)]))
ci7 <- function(x) qt(0.975, length(x) - 1) * sd(x) / sqrt(length(x))
tab7 <- function(L) imap_dfr(L, ~ tibble(
Choice = .y, mean = mean(.x), lo = mean(.x) - ci7(.x),
hi = mean(.x) + ci7(.x)))
seedT <- tab7(S); blokT <- tab7(Bk)
Model-chosen vs random vs the Section 3/5 closeness baseline,
SIR at \(\beta^*=2\beta_c\), 40
realisations.
| Seeds (maximise spread) |
Random |
274 |
[241, 308] |
| Seeds (maximise spread) |
Closeness (S3) |
313 |
[300, 325] |
| Seeds (maximise spread) |
ML model |
303 |
[279, 327] |
| 5% blockers (minimise spread) |
Random 5% |
230 |
[194, 266] |
| 5% blockers (minimise spread) |
Closeness 5% (S5) |
145 |
[110, 180] |
| 5% blockers (minimise spread) |
ML 5% |
113 |
[83, 143] |
Result. The structural model is predictive: the
two-part fit gives a logistic AUC of 0.98 for “is this node reliably
reached?” and a random-forest \(R^2=0.38\) for time-to-infection on the
core (linear baseline RMSE 4.84 vs RF 4.52). Its most important
predictors are degree, comm_size, eigen, with
degree on top — directly recovering Lecture 5 slide 51,
“nodes with a higher degree get infected faster”. As
seeds, the model infects 303 on average, well above
random (274) and statistically on par with the single best centrality of
Section 3 (313; CIs overlap). As 5% blockers (repeating
Section 5 with this learned knowledge) it is the strongest of
all: only 113 infected, below closeness (145) and far below random
(230). The learned multi-feature rule edges out the best single
centrality at minimising spread while matching it at maximising
spread: the same seed/blocker duality as Section 6, now obtained from
structure alone.
Handling censoring in the ML target. Resolved by the
data: across the 60 runs only 35% of nodes are “core” (infected in a
majority of runs, hence a reliable time-to-infection); most of
the rest are reached only occasionally, giving high-variance t_inf
estimates (just 3% are never infected at all). Pooling every node into
one regression would fit that noise. The two-part model (logistic for
whether a node is reliably reached, regression of when
on the core) handles the censoring explicitly, with its real core/never
fractions reported rather than hidden by silently dropping nodes.
8. Generality on a larger network (arXiv ca-GrQc)
Copenhagen has 662 nodes; the brief’s guideline is \(>1{,}000\). As explained in Section 0,
no Copenhagen layer can reach that, so we additionally analyse
a second, larger social network — the arXiv General Relativity
collaboration network (Leskovec, Kleinberg & Faloutsos, 2007; SNAP).
This both satisfies the size guideline and, more usefully, tests
whether our conclusions generalise. Because every routine above
is a function of an igraph object, the identical code runs here
— any difference is a property of the network, not the
implementation.
g_ca <- load_cagrqc() # downloads/caches; giant component
stopifnot(vcount(g_ca) > 1000, ecount(g_ca) > 5000,
components(g_ca)$no == 1) # the brief's size rule, met here
d_ca <- degree(g_ca)
bc_ca <- mu_global * mean(d_ca) / (mean(d_ca^2) - mean(d_ca))
The two networks. ca-GrQc satisfies the >1000-node and
>5000-edge guideline; it is sparser and a different domain, so it is
a genuine generality test, not a clone.
| Nodes |V| |
662 |
4,158 |
| Edges |E| |
9,251 |
13,422 |
| Mean degree |
27.95 |
6.46 |
| <k^2> |
1084 |
116 |
| Molloy-Reed <k^2>/ |
38.8 |
18.0 |
| beta_c (mu=0.1) |
0.00265 |
0.00589 |
set.seed(42)
NREP8 <- 30L
bstar8 <- mult_star * bc_ca # SAME R0 multiple as Section 3
nsc <- n_seed(g_ca)
cca <- list(degree = degree(g_ca), closeness = closeness(g_ca),
betweenness = betweenness(g_ca),
eigenvector = eigen_centrality(g_ca)$vector,
pagerank = page_rank(g_ca)$vector, kcore = coreness(g_ca))
tb <- function(v, k) order(v, decreasing = TRUE)[seq_len(k)]
set.seed(42); cm8 <- membership(cluster_louvain(g_ca))
seeds_comm8 <- function(k) {
tt <- sort(table(cm8), decreasing = TRUE)
al <- floor(as.numeric(tt) / sum(tt) * k)
while (sum(al) < k) { j <- which.max(as.numeric(tt)/sum(tt)*k - al)
al[j] <- al[j] + 1 }
unlist(lapply(seq_along(tt), function(i) {
idx <- which(cm8 == as.integer(names(tt)[i]))
if (al[i] == 0) return(integer(0))
idx[order(cca$degree[idx], decreasing = TRUE)][seq_len(al[i])] }))
}
strat8 <- c(list(random = NULL),
lapply(cca, tb, k = nsc),
list(community = seeds_comm8(nsc)))
arm8 <- function(ss) mean(replicate(NREP8, final_size(run_sir(
g_ca, bstar8, mu_global,
if (is.null(ss)) sample(vcount(g_ca), nsc) else ss))))
set.seed(42)
S8 <- vapply(strat8, arm8, numeric(1))
nm8 <- setdiff(names(S8), "random")
best8 <- nm8[which.max(S8[nm8])]
seed_gain8 <- unname(S8[best8] - S8["random"])
set.seed(42)
nb8 <- round(0.05 * vcount(g_ca))
mkr <- function(rm) induced_subgraph(g_ca, setdiff(seq_len(vcount(g_ca)), rm))
blk8 <- function(gr) { nse <- max(1L, floor(0.01 * vcount(gr)))
mean(replicate(NREP8, final_size(run_sir(gr, bstar8, mu_global,
sample(vcount(gr), nse))))) }
set.seed(42)
b_none <- blk8(g_ca)
b_rand <- blk8(mkr(sample(vcount(g_ca), nb8)))
b_clo <- blk8(mkr(order(cca$closeness, decreasing = TRUE)[seq_len(nb8)]))
# Copenhagen reference numbers (from Sections 3 and 5, lazy-loaded)
cph_seed <- res3$fin_m[res3$strategy == "closeness"] -
res3$fin_m[res3$strategy == "random"]
cph_blk <- bsum$fin_m[bsum$Scenario == "5% closeness"] -
bsum$fin_m[bsum$Scenario == "No blockers"]
Replication at the same operating point \(R_0\approx2\) (30 realisations on
ca-GrQc).
| Targeted seeds beat random (best strategy, Δ
infected) |
closeness, +74 |
community, +187 |
| 5% closeness blockers vs no blockers (Δ infected) |
-107 |
-128 |
| Seed/blocker duality (same centrality, opposite
sign) |
yes |
yes |
Result. All three main findings replicate on
ca-GrQc. The heterogeneous-mean-field threshold still applies (\(\beta_c=0.00589\), larger than Copenhagen’s
because the collaboration graph is sparser, with \(\langle k\rangle\approx6.5\) vs 28). At the
same \(R_0\approx2\) operating point, a
centrality-targeted seed set again beats random (community, +187
infected vs random), removing the 5% highest-closeness nodes again
suppresses the epidemic far more than no removal (-128), and 5%
random removal again does much less (-31). The accelerant/brake
duality of Section 6 therefore holds on a second, structurally different
network. The single difference is which centrality leads
(community on ca-GrQc, closeness on Copenhagen), which is itself
expected and motivates the Section 7 learned multi-feature rule as the
more portable tool.
Conclusion
One structural fact organises every result. The second moment \(\langle k^2\rangle \gg \langle k\rangle^2\)
gives the network a tiny epidemic threshold (Section 1); large seed
leverage from a 1% targeted set near \(\beta_c\) (Section 3); and a sharp
asymmetry between random and targeted 5% blockers (Section 5). The same
closeness-central nodes act as both the best accelerant and the best
brake (Section 6’s seed/blocker duality), and a learned multi-feature
model recovers the same nodes from structure alone, edging out any
single centrality at minimising spread (Section 7). At the near-critical
operating point \(\beta^*=2\beta_c\) a
short reversible quarantine is statistically swamped by run-to-run
variability (Section 4); flattening the curve requires a sustained or
larger intervention. All findings replicate on a structurally different
second network (Section 8), confirming they are properties of
heterogeneous social networks rather than artefacts of the Copenhagen
cohort.
Sources and acknowledgments
Data. Sapiezynski, P., Stopczynski, A., Lassen, D.
D., & Lehmann, S. (2019). Interaction data from the Copenhagen
Networks Study. Scientific Data, 6, 315. — Leskovec, J.,
Kleinberg, J., & Faloutsos, C. (2007). Graph evolution:
densification and shrinking diameters. ACM TKDD, 1(1) (arXiv
ca-GrQc; SNAP).
References.
- Pastor-Satorras, R., & Vespignani, A. (2001). Epidemic spreading
in scale-free networks. Physical Review Letters, 86,
3200–3203.
- Pastor-Satorras, R., & Vespignani, A. (2002). Immunization of
complex networks. Physical Review E, 65, 036104.
- Kempe, D., Kleinberg, J., & Tardos, É. (2003). Maximizing the
spread of influence through a social network. KDD ’03.
- Albert, R., Jeong, H., & Barabási, A.-L. (2000). Error and
attack tolerance of complex networks. Nature, 406,
378–382.
- Christakis, N. A., & Fowler, J. H. (2010). Social network
sensors for early detection of contagious outbreaks. PLoS ONE,
5(9), e12948.
- Sekara, V., & Lehmann, S. (2014). The strength of friendship
ties in proximity sensor data. PLoS ONE, 9(7), e100915.
Group submission: Tone Varberg Sabri, Abdullah Tadmuri, Maiheliya
Maimaitimin. Master in Computational Social Science, UC3M.
Reproducibility
R version: R version 4.5.1 (2025-06-13)
Platform : aarch64-apple-darwin20
Key packages:
caret 7.0-1
cowplot 1.1.3
data.table 1.17.2
dplyr 1.1.4
forcats 1.0.0
ggplot2 4.0.0
ggraph 2.2.2
ggthemes 5.1.0
gridExtra 2.3
igraph 2.2.0
knitr 1.50
lattice 0.22-7
lubridate 1.9.4
purrr 1.2.0
randomForest 4.7-1.2
readr 2.1.5
scales 1.4.0
stringr 1.6.0
tibble 3.3.0
tidyr 1.3.1
tidyverse 2.0.0
