---
title: "Designing a stepped-wedge trial"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Designing a stepped-wedge trial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(stepwedgepower)
```

A stepped-wedge trial assigns clusters to *sequences*. Every sequence starts
under control; sequences then cross over to the intervention at different
periods, and once a sequence has crossed it stays crossed. By the end of the
trial every cluster has received the intervention.

The package encodes exactly that and nothing more. A cluster can be a hospital
ward, a clinic, a school, a general practice, a village, or an individual
clinician; a sequence can be a region, a wave, a specialty, or an arbitrary
randomisation group. The software does not care.

Design and assumptions are kept apart, because they answer different questions:

* [`sw_design()`] -- **who** crosses over, and **when**.
* [`sw_assumptions()`] -- **what** the outcome is expected to do.

## The design

The most common design gives one baseline period and then crosses one sequence
per period.

```{r}
design <- sw_design(
  clusters_per_sequence = c(8, 8, 8, 8),
  crossover_period = c(2, 3, 4, 5),
  n_periods = 5
)
design
```

Designs are not restricted to that shape. Sequences may hold unequal numbers of
clusters, several sequences may cross at once, and there may be more than one
baseline period. Anything expressible as a 0/1 schedule is expressible here:

```{r}
schedule <- matrix(
  c(0, 0, 1, 1, 1,
    0, 0, 0, 1, 1,
    0, 0, 0, 0, 1),
  nrow = 3, byrow = TRUE
)

custom <- sw_design(
  clusters_per_sequence = c(8, 10, 12),  # unequal
  treatment = schedule                   # two baseline periods
)
custom
```

## The assumptions

```{r}
assumptions <- sw_assumptions(
  baseline_prob = 0.05,     # control-arm risk
  treatment_or = 1.6,       # target effect
  icc = 0.05,               # cluster correlation
  period_effects = c(0, 0.05, 0.10, 0.15, 0.20),  # secular trend, log-odds
  n_per_cluster_period = 20
)
assumptions
```

Three points worth drawing out.

**Specify the ICC, not the random-effect SD.** Most people can defend an ICC and
almost nobody can defend a random-intercept standard deviation. Supply `icc` and
the conversion is handled; supply `cluster_sd` if you prefer. Supplying both is
an error unless they agree.

```{r}
icc_to_cluster_sd(0.05)
cluster_sd_to_icc(0.4161)
```

**Model the secular trend explicitly.** Outcomes often drift over time for
reasons unrelated to the intervention. In a stepped-wedge design time and
treatment are partially confounded by construction, so a trend you do not
simulate is a trend you have not accounted for. `period_effects` puts it in the
generating model, where the analysis has to contend with it.

**Cluster sizes need not be equal.** Supply a single number, a
sequence-by-period matrix, or a function:

```{r}
variable_size <- sw_assumptions(
  baseline_prob = 0.05, treatment_or = 1.6, icc = 0.05,
  n_per_cluster_period = function(n) rpois(n, lambda = 20) + 1
)
trial <- simulate_swcrt(design, variable_size, seed = 1)
range(trial$n)
```

## Simulating and analysing one trial

```{r}
trial <- simulate_swcrt(design, assumptions, seed = 123)
head(trial)
```

The columns are deliberately neutral: `cluster_id`, `sequence`, `period`,
`intervention`, `n`, `events`. Any aggregated cluster-period dataset in this
shape can be fed to the analysis function, whether it came from
`simulate_swcrt()` or from a real trial.

```{r}
fit <- analyze_swcrt(trial)
fit$estimate    # log-odds
fit$p_value
```

The default analysis is the standard Hussey-Hughes model: treatment, a
categorical period effect, and a cluster random intercept. You can rename any
column, or hand over a formula outright:

```{r, eval = FALSE}
fit_stepwedge_model(
  trial,
  formula = cbind(events, n - events) ~ intervention + factor(period) +
    (1 | cluster_id)
)
```

## Power

```{r, eval = FALSE}
power <- power_swcrt(design, assumptions, nsim = 1000, seed = 1, n_cores = 4)
summary(power)
plot(power)
```

The returned object reports more than a single number, and the extras are not
decoration. The Monte Carlo interval tells you how much of the estimate is
simulation noise. The failure and singular-fit rates tell you whether the
analysis model is actually estimable at this design size -- a "power" computed
from replicates that mostly failed to converge is not a power. Bias and coverage
tell you whether the estimator is behaving.

The interval for power is exact (Clopper-Pearson) rather than Wald. A Wald
interval collapses to zero width when the estimate hits 0 or 1, which is exactly
the high-power regime a power analysis is trying to land in.

## Sensitivity

Power depends on assumptions you cannot verify in advance, so vary them.

```{r, eval = FALSE}
grid <- power_grid(
  design, assumptions,
  vary = list(
    treatment_or = c(1.3, 1.5, 1.75, 2.0),
    icc = c(0.02, 0.05, 0.10),
    clusters_per_sequence = list(c(6, 6, 6, 6), c(8, 8, 8, 8), c(10, 10, 10, 10))
  ),
  nsim = 500, n_cores = 4, seed = 1
)
summary(grid)
plot(grid)
```

Simulations are reproducible under a given `seed` regardless of `n_cores`: each
replicate draws from its own random-number stream, so adding cores makes the run
faster without moving the answer.

## Checking the engine

`validate_engine()` re-runs the power calculation with an independently written
reference implementation that shares no code with the main engine, and reports
whether the two agree within Monte Carlo error.

```{r, eval = FALSE}
validate_engine(design, assumptions, nsim = 500, seed = 1)
```

`compare_with_swdpwr()` will also line the simulation up against a closed-form
calculation, if you have that package installed. Treat the result as information
rather than as a target: closed-form stepped-wedge power is asymptotic in the
number of clusters, and with the small cluster counts these trials usually run
it can report power well above what a simulated trial achieves. Where the two
disagree, the simulation is the one that reflects the model you will actually
fit.
