---
title: "brood(): Vaccine Coverage Data Structures for All Study Designs"
subtitle: "Pre-aggregated and cohort models including time-series for interrupted time series analysis"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{brood(): Vaccine Coverage Data Structures}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  eval     = TRUE
)
```

## Overview

```{r load}
library(mudnester)
```

`brood()` constructs a `brood_df` object for vaccine coverage analysis. Unlike
`roost()` which counts events over time, vaccine coverage requires an explicitly
documented eligible-population denominator and a declared assessment window.
`brood()` makes both explicit in every object it produces.

The output `brood_df` is designed for `bowerbird::brood_plot()`.

### Why "brood"?

> *A brood is the full clutch under a parent bird's care -- every egg counted,
> every hatchling tracked, none overlooked. The gap between the clutch and the
> hatchlings is not an absence; it is information.*

### Two population models

| Model | Use case | Input |
|---|---|---|
| `"pre_aggregated"` | Already-computed counts per stratum | Summary data frame |
| `"cohort"` | Record-level, one row per person | Wide or long dose data |

The `"cohort"` model covers single time-point snapshots, birth cohort designs
with person-time, and time-series sweeps for interrupted time series analysis.

---

## Model 1: Pre-aggregated denominator

Supply a data frame with `n_vaccinated` and `n_eligible` already computed.
No per-person dose logic is performed.

```{r pre-aggregated}
cov_static <- brood(
  data.frame(
    stratum      = c("0-17", "18-49", "50-64", "65+"),
    n_vaccinated = c(480L,   3550L,   2370L,   1760L),
    n_eligible   = c(1000L,  5000L,   3000L,   2000L)
  ),
  denominator_notes = "ABS ERP 2024, Sunshine Coast LGA."
)
print(cov_static)
```

### Multi-vaccine product comparison

```{r multi-vaccine}
cov_multi <- brood(
  data.frame(
    stratum      = rep(c("18-49", "50-64", "65+"), each = 2),
    vaccine_type = rep(c("Moderna XBB.1.5", "Pfizer XBB.1.5"), 3),
    n_vaccinated = c(1800L, 1750L, 1100L, 1270L, 800L, 960L),
    n_eligible   = c(5000L, 5000L, 3000L, 3000L, 2000L, 2000L)
  ),
  vaccine_type_col  = "vaccine_type",
  denominator_notes = "AIR enrolled persons, SC HHS 2024."
)
print(cov_multi)
```

---

## Model 2: Cohort -- single time point

One row per person in `data`. Use `entry_date_col` and one of the five
assessment windows. Inputs come from `starling::murmuration()` in wide format.

```{r cohort-data}
set.seed(99)
n <- 120L
dialysis_cohort <- data.frame(
  patient_id          = paste0("D", seq_len(n)),
  dialysis_start_date = seq(as.Date("2022-01-01"), by = "week", length.out = n),
  dialysis_end_date   = seq(as.Date("2022-01-01"), by = "week", length.out = n) +
    sample(90:900, n, replace = TRUE),
  age_group           = sample(c("45-64","65-74","75+"), n, replace = TRUE,
                               prob = c(0.25, 0.45, 0.30)),
  vax_date_1          = as.Date(ifelse(
    sample(c(TRUE, FALSE), n, replace = TRUE, prob = c(0.65, 0.35)),
    as.character(seq(as.Date("2022-01-01"), by = "week", length.out = n) +
                   sample(-400:400, n, replace = TRUE)),
    NA_character_)),
  vax_type_1          = sample(c("Influenza","COVID-19",NA), n, replace = TRUE,
                               prob = c(0.5, 0.3, 0.2)),
  stringsAsFactors = FALSE
)
```

### Window 2a: Vaccination status at exit

```{r cohort-at-exit}
cov_exit <- brood(
  dialysis_cohort,
  data_format      = "wide",
  population_model = "cohort",
  vax_date_cols    = "vax_date_1",
  entry_date_col   = "dialysis_start_date",
  exit_date_col    = "dialysis_end_date",
  window           = "at_exit",
  stratum_col      = "age_group",
  denominator_notes = "Dialysis cohort SC HHS. At-exit = vaccinated by dialysis end date."
)
print(cov_exit)
```

### Window 2b: Baseline (pre-entry) vaccination status

```{r cohort-pre-entry}
cov_pre <- brood(
  dialysis_cohort,
  data_format       = "wide",
  population_model  = "cohort",
  vax_date_cols     = "vax_date_1",
  vax_type_cols     = "vax_type_1",
  vax_target        = "Influenza",
  entry_date_col    = "dialysis_start_date",
  window            = "pre_entry",
  stratum_col       = "age_group",
  denominator_notes = "Dialysis cohort. Pre-entry = vaccinated before dialysis start."
)
print(cov_pre)
```

### Window 2c: Vaccination within 1 year of cohort entry

```{r cohort-post-entry}
cov_post <- brood(
  dialysis_cohort,
  data_format       = "wide",
  population_model  = "cohort",
  vax_date_cols     = "vax_date_1",
  vax_type_cols     = "vax_type_1",
  vax_target        = "Influenza",
  entry_date_col    = "dialysis_start_date",
  window            = "post_entry_days",
  window_days       = 365L,
  stratum_col       = "age_group",
  denominator_notes = "Dialysis cohort. Post-entry window = 365 days from dialysis start."
)
print(cov_post)
```

### Window 2d: Post-intervention vaccination

```{r cohort-post-intervention}
cov_intervention <- brood(
  dialysis_cohort,
  data_format       = "wide",
  population_model  = "cohort",
  vax_date_cols     = "vax_date_1",
  vax_type_cols     = "vax_type_1",
  vax_target        = "Influenza",
  entry_date_col    = "dialysis_start_date",
  exit_date_col     = "dialysis_end_date",
  window            = "post_intervention",
  intervention_date = as.Date("2023-06-01"),
  stratum_col       = "age_group",
  denominator_notes = paste0(
    "Dialysis cohort. Post-intervention = vaccinated on or after 2023-06-01."
  )
)
print(cov_intervention)
```

---

## Model 3: Cohort -- birth cohort with person-time

Birth cohort is a special case of the cohort model. Set
`entry_date_col = "dob"` and supply `eligibility_days` to define the age-based
eligibility window. Person-time (days at risk from birth until vaccinated or
window expired) is computed automatically.

```{r birth-cohort-data}
set.seed(42)
n <- 200L
birth_records <- data.frame(
  baby_id    = paste0("B", seq_len(n)),
  dob        = seq(as.Date("2024-01-01"), by = "day", length.out = n),
  gestation  = sample(c("Term","Preterm"), n, replace = TRUE, prob = c(0.85, 0.15)),
  vax_date_1 = as.Date(ifelse(
    sample(c(TRUE, FALSE), n, replace = TRUE, prob = c(0.55, 0.45)),
    as.character(seq(as.Date("2024-01-01"), by = "day", length.out = n) +
                   sample(30:200, n, replace = TRUE)),
    NA_character_)),
  vax_type_1 = "nirsevimab",
  stringsAsFactors = FALSE
)
```

```{r birth-cohort-brood}
cov_birth <- brood(
  birth_records,
  data_format       = "wide",
  population_model  = "cohort",
  vax_date_cols     = "vax_date_1",
  vax_type_cols     = "vax_type_1",
  vax_target        = "nirsevimab",
  entry_date_col    = "dob",       # DOB is the cohort entry date
  eligibility_days  = 180L,        # eligible from birth until 6 months
  censor_date       = as.Date("2024-12-31"),
  stratum_col       = "gestation",
  denominator_notes = paste0(
    "SCPHU birth cohort 2024-01-01 to 2024-12-31. ",
    "Eligible = live births in SC LGA. ",
    "Eligibility window = 180 days from birth. ",
    "Administrative censor date = 2024-12-31."
  )
)
print(cov_birth)
```

The `person_time` column shows days at risk. `coverage` = proportion vaccinated
among those who had at least one eligible day.

---

## Model 4: Cohort -- time series (interrupted time series analysis)

This is the backbone of `brood()`. Set `time_series = TRUE` and supply
`ts_start`, `ts_end`, and `ts_by` to sweep `window = "current_at_date"` across
monthly (or other) intervals automatically. Returns one row per stratum per
time point -- the right shape for `bowerbird::brood_plot()`.

When `intervention_date` is supplied, an `intervention_period` column is
automatically added, labelling each row as `"Pre-intervention"` or
`"Post-intervention"`.

```{r time-series-data}
set.seed(7)
n <- 120L
jak_cohort <- data.frame(
  patient_id          = paste0("P", seq_len(n)),
  entry_date_nominal  = as.Date("2000-01-01"),  # placeholder entry date
  age_group           = sample(c("18-49","50-64","65+"), n, replace = TRUE),
  vax_date_1          = as.Date(ifelse(
    sample(c(TRUE, FALSE), n, replace = TRUE, prob = c(0.45, 0.55)),
    as.character(as.Date("2024-01-01") + sample(0:730, n, replace = TRUE)),
    NA_character_)),
  vax_type_1          = "Shingrix",
  vax_date_2          = as.Date(ifelse(
    sample(c(TRUE, FALSE), n, replace = TRUE, prob = c(0.25, 0.75)),
    as.character(as.Date("2024-06-01") + sample(0:365, n, replace = TRUE)),
    NA_character_)),
  vax_type_2          = "Shingrix",
  stringsAsFactors = FALSE
)
```

```{r time-series-brood}
cov_ts <- brood(
  jak_cohort,
  data_format       = "wide",
  population_model  = "cohort",
  vax_date_cols     = c("vax_date_1", "vax_date_2"),
  vax_type_cols     = c("vax_type_1", "vax_type_2"),
  vax_target        = "Shingrix",
  validity_days     = Inf,
  entry_date_col    = "entry_date_nominal",
  time_series       = TRUE,
  ts_start          = as.Date("2024-12-01"),
  ts_end            = as.Date("2026-06-01"),
  ts_by             = "month",
  intervention_date = as.Date("2025-12-01"),
  denominator_notes = "JAK-inhibitor cohort. Shingrix only, no expiry."
)
print(cov_ts)
```

The result is `r nrow(cov_ts)` rows: one per month (
`r attr(cov_ts, "n_time_points")` time points) across
`r attr(cov_ts, "n_strata")` stratum ("Overall"). With `stratum_col` supplied,
each stratum gets its own row at every time point.

### Stratified time series

```{r time-series-stratified}
cov_ts_strat <- brood(
  jak_cohort,
  data_format       = "wide",
  population_model  = "cohort",
  vax_date_cols     = c("vax_date_1", "vax_date_2"),
  vax_type_cols     = c("vax_type_1", "vax_type_2"),
  vax_target        = "Shingrix",
  validity_days     = Inf,
  entry_date_col    = "entry_date_nominal",
  stratum_col       = "age_group",
  time_series       = TRUE,
  ts_start          = as.Date("2024-12-01"),
  ts_end            = as.Date("2026-06-01"),
  ts_by             = "month",
  intervention_date = as.Date("2025-12-01"),
  denominator_notes = "JAK-inhibitor cohort by age group. Shingrix only."
)
cat("Rows:", nrow(cov_ts_strat),
    "= ", attr(cov_ts_strat, "n_time_points"), "months x",
    attr(cov_ts_strat, "n_strata"), "age groups\n")
print(head(cov_ts_strat[, c("stratum","reference_date","intervention_period",
                              "n_vaccinated","n_eligible","coverage")], 9))
```

---

## Passing to bowerbird::brood_plot()

All models return `brood_df` objects accepted by `bowerbird::brood_plot()`.
The time-series output auto-detects as a line chart:

```{r bowerbird, eval = FALSE}
# Pre-aggregated: bar chart with target line
bowerbird::brood_plot(cov_static, target_line = 0.80,
                       title = "COVID-19 booster coverage by age group")

# Time series: line chart with vertical intervention line (auto-detected)
bowerbird::brood_plot(cov_ts,
                       title = "Monthly Shingrix coverage, JAK-inhibitor cohort")

# Stratified time series: one line per stratum
bowerbird::brood_plot(cov_ts_strat,
                       title = "Shingrix coverage by age group")
```

---

## The denominator_notes requirement

Every `brood()` call should include `denominator_notes` documenting:

- Who counts as eligible (age, disease, program criteria)
- The source of the denominator (ABS ERP, ESKD registry, AIR cohort)
- The observation window (dates, censor date)
- What intervention or event defines the window anchor

`brood()` stores these notes as a metadata attribute. `print.brood_df()`
displays them and `bowerbird::brood_plot()` can use them as a plot caption.

---

## Session info

```{r session}
sessionInfo()
```