---
title: "Social bubbles on a small-world network"
author:
  - George Vega Yon
date: "`r Sys.Date()`"
vignette: >
  %\VignetteIndexEntry{Social bubbles on a small-world network}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

# Introduction

Social-bubble policies preserve contacts within small, fixed groups while
cutting down transmission between groups. How much of the outside contact is
cut is controlled by `transmission_factor`: `0` is a perfectly observed bubble
(no contact outside it at all), `1` leaves the population untouched, and values
in between represent partial compliance. This vignette compares an unmitigated
outbreak with three household-based policies created by `bubbles()`:

1. one household per bubble,
2. up to three households per bubble, and
3. up to three households per bubble, whose members all meet while the policy
   is in force (`ties = "complete"`).

All three policies operate from day 20 through day 59. In `bubbles()`, `end_day` is
exclusive, so `start_day = 20` and `end_day = 60` give that interval. The following diagram illustrates three potential bubbles:

![Three different bubbles that can be accomplished using the `bubbles()` function](img/bubble-small-diagram.svg){ width="80%"}

::: {.callout-important}
## Why this example uses SEIR

`ModelSIR()` has no exposed state and therefore cannot represent the requested
three-day incubation period. We use `ModelSEIR()`, the direct extension of the
SIR model, so that all the requested disease timings are represented. For a
strict SIR example, replace `ModelSEIR()` with `ModelSIR()` and omit
`incubation_days`.
:::

# Baseline outbreak

We simulate 10,000 agents for 100 days, beginning with 10 cases. The
small-world network has mean degree 40 and rewires 20% of its edges. With a
seven-day mean infectious period, the approximation

$$
R_0 \approx k \times \beta \times D
$$

gives a per-contact transmission probability of
$\beta = 2 / (40 \times 7)$. This is only an approximation on a clustered
network, which is why we describe the target $R_0$ as "close to 2."

```{r}
#| label: baseline-setup
library(epiworldR)

n_agents          <- 10000L
initial_cases      <- 10L
mean_degree        <- 40L
rewire_probability <- 0.20
incubation_days    <- 3
infectious_days    <- 7
target_r0          <- 2
ndays              <- 100L
nsims              <- 200L
simulation_seed    <- 2026L

transmission_rate <- target_r0 / (mean_degree * infectious_days)

baseline_model <- ModelSEIR(
  name              = "Example infection",
  prevalence        = initial_cases / n_agents,
  transmission_rate = transmission_rate,
  incubation_days   = incubation_days,
  recovery_rate     = 1 / infectious_days
)

agents_smallworld(
  baseline_model,
  n = n_agents,
  k = mean_degree,
  d = FALSE,
  p = rewire_probability
)
```

The network and disease parameters will be held fixed across policies. We also
use the same `simulation_seed` in each call to `run_multiple()`. The outbreak
size saver records the cumulative number infected, including the 10 initial
cases.

```{r}
#| label: helpers
run_scenario <- function(model, policy) {
  saver <- make_saver("outbreak_size")

  run_multiple(
    model,
    ndays    = ndays,
    nsims    = nsims,
    seed     = simulation_seed,
    saver    = saver,
    verbose  = FALSE,
    nthreads = 2L
  )

  outbreak <- run_multiple_get_results(
    model,
    nthreads = 2L
  )$outbreak_size

  final <- outbreak[outbreak$date == ndays, ]
  data.frame(
    policy        = policy,
    replicate     = final$sim_num,
    outbreak_size = final$outbreak_size
  )
}

summarize_outbreaks <- function(x) {
  pieces <- split(x$outbreak_size, x$policy)

  rows <- lapply(names(pieces), function(policy) {
    z <- pieces[[policy]]
    data.frame(
      policy        = policy,
      mean          = round(mean(z), 1),
      median        = unname(median(z)),
      lower_95_pct  = unname(quantile(z, 0.025)),
      upper_95_pct  = unname(quantile(z, 0.975)),
      mean_percent  = round(100 * mean(z) / n_agents, 1),
      row.names     = NULL
    )
  })

  do.call(rbind, rows)
}
```

# One household per bubble

First we generate household sizes from 2 through 6 with equal probability.
The only exception is at the end of the vector, where a draw that would leave
one unassigned agent is rejected. Household IDs are then assigned sequentially
in agent order.

This ordering is useful for a small-world network: before rewiring, nearby node
IDs are connected on the ring lattice. Sequential household members are
therefore likely to already share edges, even after 20% rewiring. This matters
because, by default, `bubbles()` restricts existing network ties; it does not
create new contacts between household members (we relax this in the last
policy).

```{r}
#| label: households
set.seed(2026)
household_sizes <- sample_household_sizes(n_agents)
household_id <- rep(
  seq_along(household_sizes),
  times = household_sizes
)

stopifnot(
  length(household_id) == n_agents,
  all(household_sizes >= 2L & household_sizes <= 6L)
)

table(household_sizes)
```

For the first policy, each household is its own bubble (`group_size = 1`).
Contacts within a household are untouched -- those are what the policy
preserves. Contacts outside the household are reduced by 80% during the policy
window, meaning that 20% of the original per-contact probability remains
(`transmission_factor = 0.20`), which stands in for imperfect compliance rather
than a hermetic bubble.

```{r}
#| label: household-only-model
household_only_model <- clone_epi(baseline_model)

bubbles(
  household_only_model,
  household_id        = household_id,
  flavor              = "household",
  group_size          = 1L,
  transmission_factor = 0.20,
  start_day           = 20L,
  end_day             = 60L
)
```

`bubbles()` does not keep the factor for itself: it stores it as a model
parameter, which the intervention reads on every exposure. Compliance can
therefore be inspected or changed -- between runs, or from a global event in the
middle of one -- without rebuilding the policy.

```{r}
#| label: transmission-factor-param
get_param(household_only_model, "Bubble transmission factor")
```

# Up to three households per bubble

The second policy changes only `group_size`. `bubbles()` joins connected
households into disjoint groups containing at most three households. The same
80% out-of-bubble reduction and day 20--60 policy window apply.

```{r}
#| label: three-household-model
three_household_model <- clone_epi(baseline_model)

bubbles(
  three_household_model,
  household_id        = household_id,
  flavor              = "household",
  group_size          = 3L,
  transmission_factor = 0.20,
  start_day           = 20L,
  end_day             = 60L
)
```

# Completing the bubbles

The two policies above realize a bubble only as a transmission rule: contacts
inside the bubble are left alone and everything else is damped. That
under-states what a bubble is. If two households were merged because one member
of each happened to be tied, that single contact is the only one preserved
between them -- the rest of the two households never meet, although the policy
says they have bubbled. The same gap exists inside a household, whose members
need not all be tied to each other.

With `ties = "complete"`, `bubbles()` completes every bubble to a clique while
the policy is in force: each pair of agents sharing a bubble is tied, so
households become complete and merged households meet in full. Grouping is
unchanged -- the same seed draws the same bubbles -- and the ties are withdrawn
when the policy lifts (and at the end of each run), so the network is left as it
was found. Contacts outside the bubble are still damped by the same 80%.

```{r}
#| label: complete-ties-model
complete_ties_model <- clone_epi(baseline_model)

bubbles(
  complete_ties_model,
  household_id        = household_id,
  flavor              = "household",
  group_size          = 3L,
  transmission_factor = 0.20,
  start_day           = 20L,
  end_day             = 60L,
  ties                = "complete"
)
```

Every tie is an independent daily exposure, so completing a bubble raises the
force of infection inside it. Here a bubble holds up to three households of two
to six members, so an agent may be tied to as many as 17 bubble-mates, all of
them transmitting at the full, undamped rate while the policy lasts.

# Compare outbreak sizes

Finally, we run 200 replicates for each scenario and report the final outbreak
size. Each replicate draws its own bubbles from its own seed, and the partition
belongs to the model rather than to the intervention, so the replicates can be
spread across threads.

```{r}
#| label: run-scenarios
#| message: false
results <- rbind(
  run_scenario(baseline_model, "No bubbles"),
  run_scenario(household_only_model, "One household"),
  run_scenario(three_household_model, "Up to three households"),
  run_scenario(complete_ties_model, "Up to three households, complete")
)

summarize_outbreaks(results)
```

The table reports the mean, median, central 95% simulation interval, and mean
percentage of the population ever infected. Comparing the first two bubble
policies shows the trade-off directly: allowing households to combine preserves
more contacts, but it generally produces larger outbreaks than keeping each
household in its own bubble. Completing the bubbles adds the contacts that the
policy allows but the network lacked, and the outbreak grows accordingly. The
two settings bracket the policy: `ties = "existing"` assumes bubbled households
only meet through the contacts they already had, while `ties = "complete"`
assumes everyone in a bubble meets everyone else in it every day.