---
title: "Published sample-size examples with the gsDesignNB AI skill"
output: rmarkdown::html_vignette
bibliography: gsDesignNB.bib
vignette: >
  %\VignetteIndexEntry{Published sample-size examples with the gsDesignNB AI skill}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r, message=FALSE, warning=FALSE}
library(gsDesignNB)
```

## Purpose

This vignette demonstrates how the `SKILL.md` file in the **gsDesignNB**
repository helps an AI coding assistant turn natural-language trial
descriptions into transparent `sample_size_nbinom()` calls. The examples are
not meant to replace the original protocols or statistical review; they show
how the skill keeps the translation from text to code disciplined. The skill
teaches the AI to:

1. **Extract** stated parameters from the text.
2. **Prompt** for essential information that is missing.
3. **Apply sensible defaults** for non-essential parameters.
4. **Handle ambiguities** (dropout meaning, $k$ parameterization).
5. **Compute and report** the sample size with full transparency.

## What is a skill file?

A *skill file* (`SKILL.md`) is a structured document that teaches an AI
assistant domain-specific knowledge — in this case, how to translate a
clinical trial description into a `sample_size_nbinom()` call. Unlike a
general prompt, a skill file encodes parameter extraction rules,
validation checks, unit-conversion logic, and domain heuristics (e.g.,
"COPD dispersion is typically 0.5–1.0") so that the AI can work
reliably without requiring the user to know the package API.

The main **gsDesignNB** skill covers sample size calculation, group sequential
design, blinded SSR, simulation, and interim data cuts. This vignette focuses
on the sample size workflow.

## Using the skill with an AI assistant

In the source repository, the main skill lives at
`.agents/skills/gsdesignnb/SKILL.md`. A companion simulation entrypoint lives at
`.agents/skills/gsdesignnb-simulation-entrypoint/SKILL.md`.

How you provide the skill depends on the assistant. Workspace-aware coding
assistants may discover the file automatically when the repository is open. For
assistants that accept attached files, project knowledge, or system
instructions, attach or paste `SKILL.md`, then provide the trial description in
natural language.

The skill is model-agnostic — it uses structured Markdown with
extraction tables, decision rules, and worked examples that any
instruction-following LLM can interpret.

## Overview of examples

The examples span three disease areas, four published trial protocols,
and two methodology papers — deliberately chosen to exercise different
challenges the skill must handle:

| Example | Disease | Design | Key skill challenge |
|:--------|:--------|:-------|:--------------------|
| FLAME [@wedzicha2016flame] | COPD | Non-inferiority | Dropout interpretation (inflation vs exposure integral); inferring follow-up from "annual rate" |
| IMPACT [@lipson2018impact] | COPD | 3-arm superiority | Unequal allocation; deriving control rate from stated treatment rate + % reduction; 28-day event gap |
| OPERA I/II [@hauser2017ocrelizumab] | MS | Superiority (paired studies) | Protocol sized by t-test, not NB; dispersion $k$ must be assumed |
| Friede & Schmidli [@friede2010blinded] | COPD | Fixed superiority + NI | Pooled-rate parameterization ($\lambda$ overall, not per-arm); $\phi = k$ notation |
| Mutze et al. [@mutze2019group] | HF + MS | Group sequential | Inverse dispersion ($\kappa = 1/k$); variable follow-up from staggered entry; monthly vs annual units |
| PARADIGM-HF [@mcmurray2014paradigm] | HF | Event-driven GS | Primary endpoint is not NB; recurrent-event parameters must be sourced externally |

The first example (FLAME) walks through the full skill workflow in
detail. Subsequent examples focus on the extraction and computation
steps, highlighting where each trial poses a new challenge.

## FLAME trial

We replicate the sample size from the FLAME COPD exacerbation protocol
[@wedzicha2016flame].

### The user's input

In a real session, the user might paste a protocol paragraph or summarize the
design as follows:

> Size a non-inferiority COPD exacerbation trial with equal annual rates of
> 3.0, non-inferiority margin 1.15, one-sided alpha 0.025, greater than 95%
> power, dispersion k = 0.625, and 30% dropout or major protocol violations.
> The published target was 3332 randomized patients.

The rest of this vignette walks through what the skill instructs the AI
to do with this input.

### Skill Step 1: Extract what is stated

Following the skill's extraction table, the AI parses the text and
classifies each parameter:
 
| Parameter | Found? | Extracted value | Notes |
|:----------|:------:|:----------------|:------|
| Control event rate | Yes | 3.0 events/year | "the underlying annual rate is 3.0" |
| Treatment effect | Yes | RR = 1.0 | "no difference between treatment groups" |
| Dispersion | Yes | $k = 0.625$ | "index parameter k = 0.625" |
| Power | Yes | > 95% | "greater than 95% power" |
| Alpha | Yes | 0.025, one-sided | "at the level of 0.025 (one-sided)" |
| Null hypothesis | Yes | Non-inferiority, $\text{rr}_0 = 1.15$ | "rule out a 1.15-fold increase" |
| Follow-up | **Inferred** | 1 year | Implied by "annual rate" but not explicitly stated |
| Dropout | Yes | 30% | "30% dropped out or major protocol violators" |
| Event gap | **Not stated** | — | |
| Enrollment | **Not stated** | — | |
| Trial duration | **Not stated** | — | |

### Skill Step 2: Prompt for essentials

Everything essential is present or inferable. In a real session, the AI
would confirm the follow-up duration:

> *"The protocol states an annual rate but does not explicitly give the
> follow-up duration. I'll assume 1 year of maximum follow-up. Is that
> correct?"*

If the dispersion value seemed unusual for the disease area, the skill
would also prompt: *"Is k = 0.625 the overdispersion parameter
($\text{Var} = \mu + k\mu^2$) or the Gamma shape parameter ($1/k$)?"*
Here 0.625 is consistent with published COPD dispersion values, so no
flag is raised.

### Skill Step 3: Apply defaults for non-essentials

The skill instructs the AI to default these with a note:

- **Event gap:** "No event gap specified; assuming `event_gap = 0`."
- **Accrual rate and trial duration:** "Accrual pattern does not affect
  the fixed-design N when all patients have the same max follow-up.
  Placeholder values used."

### Skill Step 4: Handle dropout correctly

The protocol says "30% dropped out or major protocol violators." The
skill recognises this as a combined dropout-plus-violator figure that
reduces the evaluable sample. Following the skill's guidance, the AI
uses the inflation-factor approach: compute N assuming full follow-up,
then divide by 0.70.

### Skill Step 5: Confirm time units

All inputs are in **years**, matching the protocol's "annual rate."

### Skill Step 6: Compute the sample size

With all parameters extracted and defaults applied, the AI translates
directly to R code:

```{r}
# All rates and durations in years
lambda1         <- 3.0                  # control exacerbation rate (events/year)
lambda2         <- 3.0                  # experimental rate (no difference under H1)
dispersion      <- 0.625                # NB dispersion k
rr0             <- 1.15                 # non-inferiority margin
event_gap       <- 0                    # no gap specified in protocol
max_followup    <- 1                    # 1-year follow-up
dropout_retained <- 0.70                # 30% dropout/violator -> inflate by 1/0.70
alpha           <- 0.025
power           <- 0.95

# Arbitrary accrual (does not affect N when all patients get the same max FU)
accrual_rate    <- 100                  # subjects/year (placeholder)
accrual_dur     <- 1                    # 1-year enrollment (placeholder)
trial_duration  <- 2                    # years
```

#### Skill-compatible replication

Because the FLAME wording combines dropout and major protocol violations, the
skill treats the 30% as an inflation factor rather than as pure time-to-dropout.
First compute the evaluable sample size under full follow-up, then divide by
0.70.

```{r}
wald_evaluable <- sample_size_nbinom(
  lambda1          = lambda1,
  lambda2          = lambda2,
  dispersion       = dispersion,
  power            = power,
  alpha            = alpha,
  rr0              = rr0,
  accrual_rate     = accrual_rate,
  accrual_duration = accrual_dur,
  trial_duration   = trial_duration,
  dropout_rate     = 0,
  max_followup     = max_followup,
  event_gap        = event_gap,
  test_type        = "wald"
)

n_inflated <- ceiling(ceiling(wald_evaluable$n_total) / dropout_retained)
```

If the 30% represented true dropout only, an alternative is to model it through
the exposure integral by converting cumulative dropout to an exponential hazard.
That gives partial credit for follow-up before dropout, so it is smaller than a
sample-size inflation.

```{r}
wald_exposure <- sample_size_nbinom(
  lambda1          = lambda1,
  lambda2          = lambda2,
  dispersion       = dispersion,
  power            = power,
  alpha            = alpha,
  rr0              = rr0,
  accrual_rate     = accrual_rate,
  accrual_duration = accrual_dur,
  trial_duration   = trial_duration,
  dropout_rate     = -log(dropout_retained) / max_followup,
  max_followup     = max_followup,
  event_gap        = event_gap,
  test_type        = "wald"
)

data.frame(
  Method           = c("Full follow-up + inflate / 0.70",
                        "Dropout via exposure integral",
                        "Protocol target"),
  N                = c(n_inflated, ceiling(wald_exposure$n_total), 3332)
) |> knitr::kable(caption = "Sample size comparison: two dropout approaches vs protocol.")
```

The protocol's 3332 falls between the two approaches and is closer to the
inflation-factor calculation. That is expected because the protocol quantity
combines true dropout with major protocol violations. For a design where dropout
is purely loss of follow-up, modeling dropout through the exposure integral is
usually more informative; for a replication of this protocol statement, the
inflation calculation is the cleaner match.

The object summary below is the evaluable Wald calculation before the 30%
randomized-sample inflation is applied.

```{r}
summary(wald_evaluable)
```

### Extensions beyond the core skill

The skill's extraction workflow ends with the computed sample size above.
The remaining sections show optional extensions that a user or AI might
explore in a follow-up conversation.

#### Score sizing (comparison)

The skill recommends comparing Wald and score sizing. For
non-inferiority where the null and alternative rates differ, the
score formula's null-variance calibration can produce a meaningfully
different sample size.

```{r}
score_evaluable <- sample_size_nbinom(
  lambda1          = lambda1,
  lambda2          = lambda2,
  dispersion       = dispersion,
  power            = power,
  alpha            = alpha,
  rr0              = rr0,
  accrual_rate     = accrual_rate,
  accrual_duration = accrual_dur,
  trial_duration   = trial_duration,
  dropout_rate     = 0,
  max_followup     = max_followup,
  event_gap        = event_gap,
  test_type        = "score"
)

score_inflated <- ceiling(ceiling(score_evaluable$n_total) / dropout_retained)

summary(score_evaluable)
```

#### Compare the two sizing rules

```{r}
comparison <- data.frame(
  Sizing        = c("Wald", "Score", "Protocol target"),
  evaluable_N   = c(ceiling(wald_evaluable$n_total), ceiling(score_evaluable$n_total), NA),
  randomized_N  = c(n_inflated, score_inflated, 3332),
  events        = round(c(wald_evaluable$total_events, score_evaluable$total_events, NA), 1)
)

knitr::kable(comparison,
  caption = "Wald vs score sizing with randomized-sample inflation."
)
```

For this non-inferiority design with equal planned rates ($\lambda_1 =
\lambda_2 = 3.0$), score sizing is slightly smaller than Wald sizing. That is
not a general rule; it depends on the null margin, rates, dispersion, and the
variance reference used by the formula. The skill asks the analyst to compare
both rules rather than assume one is always more conservative.

#### What if the protocol required an event gap?

FLAME does not specify a minimum gap between exacerbations, but many COPD
protocols require 28 symptom-free days to count a "new" exacerbation. Here
we show how the Jensen correction would change the sample size if a 28-day
gap were imposed. (Since rates are in years, we convert: 28/365.25 years.)

```{r}
wald_gap_evaluable <- sample_size_nbinom(
  lambda1          = lambda1,
  lambda2          = lambda2,
  dispersion       = dispersion,
  power            = power,
  alpha            = alpha,
  rr0              = rr0,
  accrual_rate     = accrual_rate,
  accrual_duration = accrual_dur,
  trial_duration   = trial_duration,
  dropout_rate     = 0,
  max_followup     = max_followup,
  event_gap        = 28 / 365.25,       # 28 days in years
  test_type        = "wald"
)

n_gap_inflated <- ceiling(ceiling(wald_gap_evaluable$n_total) / dropout_retained)

data.frame(
  Formula   = c("No gap (FLAME protocol)", "28-day gap (Jensen-corrected)"),
  N         = c(n_inflated, n_gap_inflated),
  delta_n   = c(0, n_gap_inflated - n_inflated)
) |> knitr::kable(
  caption = "Impact of a hypothetical 28-day event gap on the FLAME sample size."
)
```

#### Verify by simulation (fixed design)

A simulation can check whether the randomized sample size delivers the target
power under the alternative ($\lambda_1 = \lambda_2 = 3.0$, so
$\text{rr} = 1$) when the analysis tests non-inferiority against
$\text{rr}_0 = 1.15$. The code below is intentionally not run during vignette
building; use thousands of replicates for a real operating-characteristic
check.

Since `mutze_test()` tests $H_0\!: \log(\text{rr}) = 0$, we shift the
z-statistic to test the NI null:

$$z_{\text{NI}} = \frac{\log(\text{rr}_0) - \hat{\beta}}{\widehat{\text{SE}}}$$

and reject when $z_{\text{NI}} \geq z_{1-\alpha}$.

```{r, eval=FALSE}
n_sims <- 5000
n_per_arm <- ceiling(n_inflated / 2)

enroll_rate <- data.frame(rate = n_per_arm * 2, duration = accrual_dur)
fail_rate <- data.frame(
  treatment  = c("Control", "Experimental"),
  rate       = c(lambda1, lambda2),
  dispersion = c(dispersion, dispersion)
)
dropout_df <- data.frame(
  treatment = c("Control", "Experimental"),
  rate      = c(-log(dropout_retained), -log(dropout_retained)),
  duration  = c(100, 100)
)

set.seed(20260506)
rejections <- 0L
for (i in seq_len(n_sims)) {
  dat <- nb_sim(
    enroll_rate  = enroll_rate,
    fail_rate    = fail_rate,
    dropout_rate = dropout_df,
    max_followup = max_followup
  )
  dat <- cut_data_by_date(dat, trial_duration)
  tst <- mutze_test(dat, test_type = "wald")
  # NI z-statistic: shift from H0: beta=0 to H0: rr >= rr0
  z_ni <- (log(rr0) - tst$estimate) / tst$se
  if (z_ni >= qnorm(1 - alpha)) rejections <- rejections + 1L
}

empirical_power <- rejections / n_sims

data.frame(
  Metric = c("Target power", "Empirical power", "Replicates"),
  Value  = c(sprintf("%.1f%%", power * 100),
             sprintf("%.1f%%", empirical_power * 100),
             n_sims)
) |> knitr::kable(caption = "Fixed-design simulation verification.")
```

At 5,000 replicates, the Monte Carlo standard error near 95% power is about
0.3 percentage points.

## IMPACT trial

The IMPACT trial [@lipson2018impact] compared triple therapy
(fluticasone furoate/umeclidinium/vilanterol, FF/UMEC/VI) to two dual
therapies (FF/VI and UMEC/VI) in COPD patients with a history of
exacerbations. The protocol (Section 8.2.1) gives exact sample size
assumptions.

### The user's input

The user might summarize the protocol assumptions this way:

> Size the two pairwise IMPACT contrasts from a 2:2:1 trial: 4000 on triple
> therapy, 4000 on FF/VI, and 2000 on UMEC/VI. Use 90% power, two-sided 1%
> alpha, a 0.80/year triple-therapy event rate, 15% reduction versus UMEC/VI,
> 12% reduction versus FF/VI, dispersion k = 0.75, 52 weeks of follow-up, and a
> 28-day event gap.

### Extraction and computation

The skill extracts the FF/UMEC/VI vs UMEC/VI comparison:

- Control rate (UMEC/VI): $0.80 / (1 - 0.15) \approx 0.941$/year
- Treatment rate (FF/UMEC/VI): 0.80/year → RR = 0.85
- Dispersion: $k = 0.75$
- Power: 90%, Alpha: 0.01 two-sided (→ 0.005 one-sided)
- Allocation: 2:1 (FF/UMEC/VI : UMEC/VI)
- Follow-up: 52 weeks ≈ 1 year
- Event gap: 28 days = 28/365.25 years

```{r}
# IMPACT: FF/UMEC/VI vs UMEC/VI
lambda_triple <- 0.80                     # FF/UMEC/VI rate (events/year)
lambda_umec   <- lambda_triple / (1 - 0.15)  # UMEC/VI rate (~0.941)

impact_wald <- sample_size_nbinom(
  lambda1          = lambda_umec,           # control = UMEC/VI
  lambda2          = lambda_triple,         # experimental = FF/UMEC/VI
  dispersion       = 0.75,
  power            = 0.90,
  alpha            = 0.005,                 # two-sided 1% → one-sided 0.5%
  sided            = 1,
  ratio            = 2,                     # 2:1 (experimental:control)
  accrual_rate     = 100,
  accrual_duration = 1,
  trial_duration   = 2,
  max_followup     = 1,                     # 52-week follow-up
  event_gap        = 28 / 365.25,           # 28-day gap
  test_type        = "wald"
)

data.frame(
  Metric = c("N control (UMEC/VI)", "N experimental (FF/UMEC/VI)",
             "N total (this comparison)", "Protocol target (comparison)"),
  Value  = c(ceiling(impact_wald$n1), ceiling(impact_wald$n2),
             ceiling(impact_wald$n_total), "2,000 + 4,000 = 6,000")
) |> knitr::kable(caption = "IMPACT: FF/UMEC/VI vs UMEC/VI.")
```

The protocol's 10,000-patient total also powers the FF/UMEC/VI vs FF/VI
comparison with a smaller 12% effect:

```{r}
# IMPACT: FF/UMEC/VI vs FF/VI
lambda_ffvi <- lambda_triple / (1 - 0.12)  # FF/VI rate (~0.909)

impact_ffvi <- sample_size_nbinom(
  lambda1          = lambda_ffvi,
  lambda2          = lambda_triple,
  dispersion       = 0.75,
  power            = 0.90,
  alpha            = 0.005,
  sided            = 1,
  ratio            = 1,                     # 1:1 (4,000 : 4,000)
  accrual_rate     = 100,
  accrual_duration = 1,
  trial_duration   = 2,
  max_followup     = 1,
  event_gap        = 28 / 365.25,
  test_type        = "wald"
)

data.frame(
  Metric = c("N per arm (FF/UMEC/VI vs FF/VI)", "N total", "Protocol target"),
  Value  = c(ceiling(impact_ffvi$n1), ceiling(impact_ffvi$n_total), "4,000 + 4,000 = 8,000")
) |> knitr::kable(caption = "IMPACT: FF/UMEC/VI vs FF/VI (12% reduction).")
```

### Sensitivity to dispersion

The protocol states $k = 0.75$ but acknowledges uncertainty. How does
sample size vary with $k$?

```{r, message=FALSE, warning=FALSE, fig.width=6, fig.height=4}
library(ggplot2)

k_vals <- seq(0.3, 1.5, by = 0.05)
n_by_k <- vapply(k_vals, function(k) {
  res <- sample_size_nbinom(
    lambda1 = lambda_umec, lambda2 = lambda_triple, dispersion = k,
    power = 0.90, alpha = 0.005, sided = 1, ratio = 2,
    accrual_rate = 100, accrual_duration = 1, trial_duration = 2,
    max_followup = 1, event_gap = 28 / 365.25, test_type = "wald"
  )
  ceiling(res$n_total)
}, numeric(1))

sens_df <- data.frame(k = k_vals, N = n_by_k)

p <- ggplot(sens_df, aes(x = k, y = N)) +
  geom_line(linewidth = 0.8, colour = "#2166AC") +
  geom_point(
    data = sens_df[sens_df$k == 0.75, ],
    colour = "#B2182B", size = 3
  ) +
  annotate("text", x = 0.75, y = sens_df$N[sens_df$k == 0.75],
           label = "Protocol assumption", vjust = -1.2,
           colour = "#B2182B", size = 3.5) +
  scale_y_continuous(labels = scales::comma) +
  labs(
    x = "Dispersion parameter k",
    y = "Total sample size N",
    title = "IMPACT: sample size sensitivity to dispersion"
  ) +
  theme_minimal(base_size = 12)

p
```

The sample size increases substantially with $k$ — motivating blinded
SSR (see the `ssr-example` vignette).

### Impact of the 28-day event gap

The Jensen correction reduces the effective rate (more events are
"absorbed" by the gap), which inflates the sample size:

```{r}
impact_nogap <- sample_size_nbinom(
  lambda1 = lambda_umec, lambda2 = lambda_triple, dispersion = 0.75,
  power = 0.90, alpha = 0.005, sided = 1, ratio = 2,
  accrual_rate = 100, accrual_duration = 1, trial_duration = 2,
  max_followup = 1, event_gap = 0, test_type = "wald"
)

data.frame(
  Design      = c("No event gap", "28-day gap (Jensen-corrected)"),
  N           = c(ceiling(impact_nogap$n_total), ceiling(impact_wald$n_total)),
  Delta       = c(0, ceiling(impact_wald$n_total) - ceiling(impact_nogap$n_total))
) |> knitr::kable(caption = "Impact of the 28-day event gap on IMPACT sample size.")
```


## OPERA I/II trials

The OPERA I and II trials [@hauser2017ocrelizumab] evaluated ocrelizumab
vs interferon beta-1a (Rebif) in relapsing multiple sclerosis. The primary
endpoint was the annualized relapse rate (ARR) at 96 weeks, analyzed with a
negative binomial model.

### The user's input

The user might summarize the published design as:

> Each OPERA study randomized about 400 patients per arm. Use a Rebif ARR of
> 0.33/year, a 50% relative reduction for ocrelizumab, two-sided alpha 0.05,
> 84% power, 96 weeks of follow-up, and 20% dropout. The protocol does not give
> an NB dispersion value for sample-size replication.

### Extraction and computation

The protocol used a t-test for sizing but the analysis used NB regression.
We can check what `sample_size_nbinom()` gives with plausible dispersion
values for MS relapses ($k \in [0.5, 1.5]$):

```{r}
# OPERA: ocrelizumab vs Rebif
opera_k_vals <- c(0.5, 0.75, 1.0, 1.5)
opera_results <- lapply(opera_k_vals, function(k) {
  res <- sample_size_nbinom(
    lambda1          = 0.33,                 # Rebif ARR (events/year)
    lambda2          = 0.165,                # ocrelizumab ARR
    dispersion       = k,
    power            = 0.84,
    alpha            = 0.025,                # two-sided 5% → one-sided 2.5%
    sided            = 1,
    ratio            = 1,
    accrual_rate     = 100,
    accrual_duration = 1,
    trial_duration   = 3,
    max_followup     = 96 / 52,              # 96 weeks ≈ 1.846 years
    dropout_rate     = -log(0.80) / (96/52), # 20% dropout over 96 weeks
    test_type        = "wald"
  )
  data.frame(k = k, N_per_arm = ceiling(res$n1), N_total = ceiling(res$n_total))
})

opera_df <- do.call(rbind, opera_results)
opera_df$Protocol <- 400

knitr::kable(opera_df,
  caption = "OPERA sample size for varying dispersion k (84% power, RR = 0.50)."
)
```

With a 50% rate reduction, this NB calculation gives fewer than 400 subjects per
arm across the assumed dispersion range. That does not mean the protocol was
over-sized: the published calculation used a different sizing model, and this
NB replication requires an assumed $k$ that was not stated in the excerpt. The
skill should therefore present this as a sensitivity analysis, not an exact
replication.


## Friede & Schmidli - formula validation

Friede & Schmidli [-@friede2010blinded] provide sample sizes for balanced COPD
superiority trials using the same Wald formula implemented in
`sample_size_nbinom()`. Their Table 1 gives fixed-design $n_0$ per group
for specific combinations of rate ratio $\theta$, overall event rate
$\lambda$, shape parameter $\phi$, and target power.

### Published $n_0$ values

For their base scenario ($\lambda^* = 1.5$/year, $\phi^* = 0.5$, equal
allocation, one-sided $\alpha = 0.025$, follow-up $t_i = 1$ year):

| $\theta$ | Power | Published $n_0$ |
|-----------|-------|-----------------|
| 0.7       | 80%   | 147             |
| 0.7       | 90%   | 196             |
| 0.8       | 80%   | 370             |
| 0.8       | 90%   | 496             |

Note that $\lambda$ in Friede & Schmidli is the overall (pooled) rate,
$\lambda = (\lambda_0 + \lambda_1) / 2$, and $\phi$ is our $k$ (dispersion
such that $\text{Var} = \mu + \phi\mu^2$). Given $\lambda = 1.5$ and
$\theta = \lambda_1/\lambda_0$:

$$\lambda_0 = \frac{2\lambda}{1 + \theta}, \quad \lambda_1 = \theta \cdot \lambda_0$$

### Replication

```{r}
friede_scenarios <- expand.grid(
  theta = c(0.7, 0.8),
  power = c(0.80, 0.90)
)
lambda_overall <- 1.5
phi <- 0.5

friede_results <- mapply(function(theta, pwr) {
  lam0 <- 2 * lambda_overall / (1 + theta)
  lam1 <- theta * lam0
  res <- sample_size_nbinom(
    lambda1          = lam0,
    lambda2          = lam1,
    dispersion       = phi,
    power            = pwr,
    alpha            = 0.025,
    sided            = 1,
    ratio            = 1,
    accrual_rate     = 100,
    accrual_duration = 1,
    trial_duration   = 2,
    max_followup     = 1,               # ti = 1 year (equal FU for all)
    dropout_rate     = 0,
    test_type        = "wald"
  )
  ceiling(res$n1)
}, friede_scenarios$theta, friede_scenarios$power)

published_n0 <- c(147, 370, 196, 496)

knitr::kable(
  data.frame(
    theta     = friede_scenarios$theta,
    power     = friede_scenarios$power,
    published = published_n0,
    gsDesignNB = friede_results,
    difference = friede_results - published_n0,
    within_one = ifelse(abs(friede_results - published_n0) <= 1, "Yes", "No")
  ),
  caption = "Friede & Schmidli (2010) Table 1: fixed-design validation."
)
```

The replicated values match the published table exactly or within one subject
per arm, which is the expected resolution for small differences in rounding
conventions.

### Non-inferiority extension

Friede & Schmidli also consider a non-inferiority design using the
COPD combination therapy (Calverley et al.) as active control, with a
15% non-inferiority margin:

- Control (combination therapy): $\lambda_0 = 1.16$/year
- Shape: $\phi = 0.46$
- Rate ratio: $\theta = 1.0$ (no difference under $H_1$)
- Margin: $\text{rr}_0 = 1.15$

```{r}
friede_ni <- sample_size_nbinom(
  lambda1          = 1.16,                # active control rate
  lambda2          = 1.16,                # experimental (no difference)
  dispersion       = 0.46,
  power            = 0.80,
  alpha            = 0.025,
  sided            = 1,
  ratio            = 1,
  rr0              = 1.15,               # NI margin
  accrual_rate     = 100,
  accrual_duration = 1,
  trial_duration   = 2,
  max_followup     = 1,
  dropout_rate     = 0,
  test_type        = "wald"
)

summary(friede_ni)
```

This parallels the FLAME NI design above, but with a lower event rate
and different dispersion — both are COPD trials with identical 15%
non-inferiority margins.


## Mutze et al. - group sequential parameterization examples

Mutze et al. [-@mutze2019group] extended the fixed-design NB sample size formula
to group sequential designs and provide Tables 4 and 5 with information
and sample size calculations under O'Brien-Fleming and Pocock spending
functions.

### Fixed-design information

The tables report the fixed-design statistical information $I_{\text{fix}}$
and corresponding $n_1$ per arm. Their shape parameter $\kappa$ is
**the inverse** of the `gsDesignNB` dispersion $k$: $\kappa = 1/k$.

**Heart failure scenarios** (Table 4): enrollment 15 months, study
48 months, variable follow-up (33–48 months), annualized rates.

```{r}
# Mutze Table 4: HF scenarios, κ=2 (k=0.5), θ=0.70
# Their λ₁=λ₂ under H₀ gives the "overall" rate; under H₁ λ₂=θ·λ₁
# We verify a subset of fixed-design n₁ values
mutze_hf <- expand.grid(
  lambda_annual = c(0.08, 0.10, 0.12, 0.14),
  theta         = c(0.70, 0.80),
  kappa         = 2,
  power         = c(0.80, 0.90)
)

mutze_hf_results <- mapply(function(lam, theta, kappa, pwr) {
  k <- 1 / kappa
  # Rates are annualized; follow-up is variable (15-mo enrollment, 48-mo study)
  # Convert to monthly for enrollment/trial_duration consistency
  lam_mo <- lam / 12
  res <- sample_size_nbinom(
    lambda1          = lam_mo,              # control rate (monthly)
    lambda2          = theta * lam_mo,      # experimental rate
    dispersion       = k,
    power            = pwr,
    alpha            = 0.025,
    sided            = 1,
    ratio            = 1,
    accrual_rate     = 100,
    accrual_duration = 15,                  # 15 months enrollment
    trial_duration   = 48,                  # 48 months study
    dropout_rate     = 0,
    test_type        = "wald"
  )
  ceiling(res$n1)
}, mutze_hf$lambda_annual, mutze_hf$theta, mutze_hf$kappa, mutze_hf$power)

mutze_hf$k <- 1 / mutze_hf$kappa
mutze_hf$n1_gsDesignNB <- mutze_hf_results

knitr::kable(
  mutze_hf[, c("lambda_annual", "theta", "k", "power", "n1_gsDesignNB")],
  caption = "Mutze et al. (2019): HF fixed-design n₁ (κ=2, k=0.5)."
)
```

**MS scenarios** (Table 5): fixed 6-month follow-up, monthly CUAL rates.

```{r}
# Mutze Table 5: MS scenarios, κ=2 (k=0.5), fixed 6-month FU
mutze_ms <- expand.grid(
  lambda_annual = c(6, 8, 10),   # annualized monthly CUAL rates
  theta         = c(0.50, 0.70),
  kappa         = 2,
  power         = c(0.80, 0.90)
)

mutze_ms_results <- mapply(function(lam, theta, kappa, pwr) {
  k <- 1 / kappa
  lam_mo <- lam / 12          # monthly rate
  res <- sample_size_nbinom(
    lambda1          = lam_mo,
    lambda2          = theta * lam_mo,
    dispersion       = k,
    power            = pwr,
    alpha            = 0.025,
    sided            = 1,
    ratio            = 1,
    accrual_rate     = 100,
    accrual_duration = 1,
    trial_duration   = 7,       # 6-month FU + enrollment
    max_followup     = 6,       # fixed 6 months
    dropout_rate     = 0,
    test_type        = "wald"
  )
  ceiling(res$n1)
}, mutze_ms$lambda_annual, mutze_ms$theta, mutze_ms$kappa, mutze_ms$power)

mutze_ms$k <- 1 / mutze_ms$kappa
mutze_ms$n1_gsDesignNB <- mutze_ms_results

knitr::kable(
  mutze_ms[, c("lambda_annual", "theta", "k", "power", "n1_gsDesignNB")],
  caption = "Mutze et al. (2019): MS fixed-design n₁ (κ=2, k=0.5)."
)
```

### Real-data reverse engineering: CHARM-Preserved

Mutze et al. Table 1 reports results from the CHARM-Preserved trial
(candesartan vs placebo in heart failure with preserved ejection fraction):

- Placebo: $n = 1509$, 4374 person-years, 547 HF admissions
- Candesartan: $n = 1514$, 4425 person-years, 392 admissions
- Estimated NB rate ratio: $\hat{\theta} = 0.71$

What sample size would `sample_size_nbinom()` recommend if we were
designing a new trial with these observed parameters?

```{r}
# Estimate rates from CHARM-Preserved
rate_placebo   <- 547 / 4374     # events per person-year
rate_treatment <- 392 / 4425
rr_observed    <- rate_treatment / rate_placebo

charm <- sample_size_nbinom(
  lambda1          = rate_placebo,
  lambda2          = rate_treatment,
  dispersion       = 0.5,           # typical HF dispersion
  power            = 0.90,
  alpha            = 0.025,
  sided            = 1,
  ratio            = 1,
  accrual_rate     = 100,
  accrual_duration = 12,
  trial_duration   = 48,            # ~4-year study
  max_followup     = 36,            # ~3-year median FU
  dropout_rate     = 0,
  test_type        = "wald"
)

data.frame(
  Metric = c("Control rate (events/person-year)", "Rate ratio",
             "Required N (90% power)", "CHARM-Preserved actual N"),
  Value  = c(round(rate_placebo, 3), round(rr_observed, 3),
             ceiling(charm$n_total), 3023)
) |> knitr::kable(caption = "Design from CHARM-Preserved observed data.")
```


### Small-sample regime: BOLD trial

The BOLD study (Mutze et al. Table 2) was a dose-ranging trial of
siponimod in relapsing MS with very small sample sizes — a regime where
the Wald test's Type I error inflation is known to be problematic:

- Placebo ($n = 61$): monthly CUALs = 1.39
- Siponimod 2 mg ($n = 45$): monthly CUALs = 0.42, RR = 0.303

```{r}
bold_wald <- sample_size_nbinom(
  lambda1          = 1.39,                # monthly CUAL rate, placebo
  lambda2          = 0.42,                # siponimod 2mg
  dispersion       = 0.5,                # assumed
  power            = 0.80,
  alpha            = 0.025,
  sided            = 1,
  ratio            = 1,
  accrual_rate     = 100,
  accrual_duration = 1,
  trial_duration   = 4,
  max_followup     = 3,                  # 3 months
  dropout_rate     = 0,
  test_type        = "wald"
)

bold_score <- sample_size_nbinom(
  lambda1          = 1.39,
  lambda2          = 0.42,
  dispersion       = 0.5,
  power            = 0.80,
  alpha            = 0.025,
  sided            = 1,
  ratio            = 1,
  accrual_rate     = 100,
  accrual_duration = 1,
  trial_duration   = 4,
  max_followup     = 3,
  dropout_rate     = 0,
  test_type        = "score"
)

data.frame(
  Test  = c("Wald", "Score"),
  N_per_arm = c(ceiling(bold_wald$n1), ceiling(bold_score$n1)),
  N_total   = c(ceiling(bold_wald$n_total), ceiling(bold_score$n_total))
) |> knitr::kable(
  caption = "BOLD trial sizing: Wald vs Score (RR = 0.30, high rate)."
)
```

At these small sample sizes, the score test is recommended for the
final analysis to maintain proper Type I error control (see Mutze et
al., 2019, Section 5.2).


## PARADIGM-HF: recurrent-event re-analysis

The PARADIGM-HF trial [@mcmurray2014paradigm] was designed with a
time-to-first-event primary endpoint (CV death or HF hospitalization) analyzed
by Cox regression. The protocol specified:

- N = 7,980 (3,990 per arm)
- CV death hazard reduction: 15%, event rate 7%/year in enalapril
- Primary composite rate: 14.5%/year → power ≥ 97%
- O'Brien-Fleming spending, 2 interim analyses at 1/3 and 2/3 information

But recurrent HF hospitalizations were also analyzed. What if the
trial had been sized for recurrent hospitalizations using an NB model?
Using rates from the gsDesignNB paper ($\sim 0.3$ events/month,
$k \approx 0.5$) and a 20% rate reduction:

```{r}
paradigm_nb <- sample_size_nbinom(
  lambda1          = 0.3,                  # HF hospitalization rate, monthly
  lambda2          = 0.3 * 0.80,           # 20% reduction
  dispersion       = 0.5,
  power            = 0.90,
  alpha            = 0.025,
  sided            = 1,
  ratio            = 1,
  accrual_rate     = 362,                  # ~7,980 over 22 months
  accrual_duration = 22,
  trial_duration   = 43,                   # 22 enrollment + 21 minimum FU
  dropout_rate     = 0,
  test_type        = "wald"
)

data.frame(
  Design      = c("NB recurrent events (hypothetical)",
                   "Cox time-to-first (protocol)"),
  N           = c(ceiling(paradigm_nb$n_total), 7980),
  Endpoint    = c("Recurrent HF hospitalizations", "CV death or first HF hosp.")
) |> knitr::kable(
  caption = "PARADIGM-HF: NB recurrent-event sizing vs actual protocol."
)
```

This illustrates a key trade-off: time-to-first-event designs require
more patients but capture mortality; recurrent-event NB designs use
information from all events and can be more efficient for rate-based
endpoints.


## Summary

This vignette applied the `SKILL.md` extraction workflow to six
published examples spanning COPD, MS, and heart failure. The table
below summarizes what the skill handled well and where challenges arose.

### Successes

| Skill step | Demonstrated in | What worked |
|:-----------|:----------------|:------------|
| Parameter extraction | FLAME, IMPACT | All stated parameters (rates, $k$, power, $\alpha$, allocation, margin) were correctly parsed from protocol text |
| Sensible defaults | FLAME, OPERA | Follow-up inferred from "annual rate"; accrual pattern recognized as irrelevant for fixed-FU sizing |
| Dispersion check | FLAME, IMPACT | Values checked against disease-area ranges (COPD: $k \in [0.5, 1.0]$); no false alarms |
| Event gap handling | IMPACT | 28-day gap automatically triggered Jensen correction; sensitivity table showed practical impact |
| Non-inferiority | FLAME, Friede NI | NI margin ($\text{rr}_0 = 1.15$) correctly distinguished from the treatment effect |
| Formula validation | Friede & Schmidli | Published $n_0$ values matched exactly or within one subject per arm, confirming the Wald formula implementation up to rounding |
| Score vs Wald comparison | FLAME, BOLD | Both sizing rules computed; difference clearly explained in context |
| Dropout modeling | FLAME | Inflation and exposure-integral approaches separated so protocol replication and prospective dropout modeling are not conflated |

### Challenges

| Challenge | Example | Issue | Resolution |
|:----------|:--------|:------|:-----------|
| Unstated dispersion | OPERA | Protocol sized by t-test; $k$ not given | Skill must ask user or sweep over plausible range ($k \in [0.5, 1.5]$) |
| Inverse parameterization | Mutze et al. | Paper uses $\kappa = 1/k$ | Skill's dispersion-check step must detect and convert — easy to miss |
| Pooled-rate formulation | Friede & Schmidli | Paper gives overall $\lambda$, not per-arm | Need algebra: $\lambda_0 = 2\lambda/(1+\theta)$ — skill should prompt for clarification |
| Unequal allocation | IMPACT | 2:2:1 three-arm design | Skill handles pairwise comparisons but user must specify which contrast |
| Wrong endpoint family | PARADIGM-HF | Primary is Cox, not NB | Skill correctly flags that NB is not the primary model; recurrent-event parameters sourced externally |
| Time-unit ambiguity | Mutze HF vs MS | Monthly CUALs vs annualized rates in same paper | Skill must confirm units before computing; mixing units is a common error source |
| Small-sample regime | BOLD | $n \approx 50$/arm — Wald test inflates Type I error | Skill recommends score test and simulation verification |

### Lessons for skill development

1. **Always confirm dispersion parameterization.** Three conventions
   appear in the literature: $k$ (our convention), $\phi = k$
   (Friede & Schmidli), and $\kappa = 1/k$ (Mutze). The skill should
   ask whenever the source notation is ambiguous.
2. **Sweep when $k$ is unknown.** When a protocol does not state the
   dispersion (OPERA, PARADIGM-HF), a sensitivity table over plausible
   values is more useful than a single-point answer.
3. **Flag non-NB primaries.** When the protocol's primary analysis is
   Cox or another model, the skill should state clearly that the NB
   calculation is exploratory or hypothetical.
4. **Dropout handling depends on wording.** If dropout is a true loss of
   follow-up, model it through the exposure integral with `dropout_rate > 0`.
   If the protocol combines dropout with non-evaluable or protocol-violator
   assumptions, an inflation factor may be the better replication.
5. **Event gap matters more at high rates.** The Jensen correction had
   a visible impact for IMPACT (COPD, rate ≈ 0.9/year, 28-day gap) but
   would be negligible for MS relapses at 0.3/year.

### Practical notes for literature replication

- **Dropout handling.** Modeling dropout through the exposure integral
  (`dropout_rate > 0`) gives partial credit for dropout follow-up. Convert
  cumulative dropout incidence $p$ to an instantaneous rate:
  `dropout_rate = -log(1 - p) / T`. Use a simple $N/(1 - p)$ inflation when the
  protocol describes non-evaluable patients, major protocol violators, or a
  conventionally inflated randomized sample size rather than a pure dropout
  process.
- **Accrual pattern.** For a fixed-design replication where all patients
  get the same maximum follow-up, the accrual shape does not affect N.
- **Event gap.** When a protocol does not mention a minimum inter-event
  gap, set `event_gap = 0`. The Jensen correction activates automatically
  when `event_gap > 0`.
- **Dispersion parameterization.** Friede & Schmidli's $\phi$ equals
  our $k$. Mutze et al.'s $\kappa$ equals $1/k$. Always check
  whether a paper uses $\text{Var} = \mu + k\mu^2$ (our convention)
  or the Gamma shape $1/k$.

For group sequential designs and SSR studies, see the
`group-sequential-simulation`, `ssr-example`, and
`ssr-simulation-study` vignettes.

## References
