Case study: an Lp(a) testing stepped-wedge trial

This vignette works through the study that the package originally grew out of: a stepped-wedge trial of an intervention intended to raise lipoprotein(a) (Lp(a)) testing rates among physicians.

It is worth being explicit about what is application-specific here and what is not. The package knows nothing about Lp(a), physicians, or specialties. It knows about clusters, sequences, periods, and aggregated binary outcomes. In this study the clusters happen to be physicians and the sequences happen to be specialties – but a different trial could map schools onto clusters and districts onto sequences and use exactly the same functions.

library(stepwedgepower)

The applied vocabulary

Study concept	Package concept
Physician	cluster
Specialty (Cardiology, …)	sequence
Study step	period
Patients seen by a physician in a step	cluster-period sample size
Patients with an Lp(a) test ordered	events

Baseline data

The package ships a small synthetic physician-level dataset. Real analyses would substitute their own extract with the same shape.

physicians <- read_example_physician_data()
str(physicians, max.level = 1)
#> 'data.frame':    16 obs. of  7 variables:
#>  $ prov_id      : int  1 2 3 4 5 6 7 8 9 10 ...
#>  $ PROV_NAME    : chr  "A. Cardio" "B. Cardio" "C. Cardio" "D. Cardio" ...
#>  $ specialty    : chr  "CARDIOLOGY" "CARDIOLOGY" "CARDIOLOGY" "CARDIOLOGY" ...
#>  $ n_total_pat  : int  950 880 730 1100 650 740 520 900 410 350 ...
#>  $ n_lpa_pat    : int  72 61 41 79 16 19 12 24 15 11 ...
#>  $ n_ldl_pat    : int  260 245 210 310 290 300 250 330 155 130 ...
#>  $ n_ldl_lpa_pat: int  30 24 18 33 10 11 8 13 9 7 ...

prepare_cluster_data() is the generic filter: keep the groups of interest and trim clusters that are implausibly small or large. prepare_physician_data() is a thin application wrapper around it with the study’s defaults.

analysis_set <- prepare_physician_data(
  physicians,
  min_patients = 100,
  max_patients = 10000
)

summarize_by_specialty(analysis_set, vars = "n_total_pat")
#>      variable         specialty n min    q1 median  mean    q3  max
#> 1 n_total_pat        CARDIOLOGY 4 730 842.5    915 915.0 987.5 1100
#> 2 n_total_pat   FAMILY MEDICINE 4 520 617.5    695 702.5 780.0  900
#> 3 n_total_pat INTERNAL MEDICINE 4 350 395.0    455 470.0 530.0  620
#> 4 n_total_pat         NEUROLOGY 4 380 410.0    465 480.0 535.0  610

Baseline testing rates by specialty

Baseline rates differ substantially across specialties, which matters because the power of the trial depends on where each sequence starts. fit_specialty_rate_model() wraps the generic fit_grouped_rate_model().

rate_model <- fit_specialty_rate_model(
  analysis_set,
  successes = "n_lpa_pat",
  trials = "n_total_pat",
  specialty_var = "specialty",
  provider_var = "prov_id",
  random_intercept = TRUE
)
#> boundary (singular) fit: see help('isSingular')

baseline_rates <- estimate_specialty_rates(rate_model, specialty_var = "specialty")
baseline_rates[, c("specialty", "probability")]
#>           specialty probability
#> 1        CARDIOLOGY  0.06912568
#> 2   FAMILY MEDICINE  0.02526690
#> 3 INTERNAL MEDICINE  0.03776596
#> 4         NEUROLOGY  0.02083333

The physician-level random-effect variance from this model is what we will carry into the power calculation as the cluster random effect. Expressed as an ICC:

tau2 <- as.numeric(lme4::VarCorr(rate_model)$prov_id[1, 1])
cluster_sd_to_icc(sqrt(tau2))
#> [1] 0

Powering the trial

Now the design. Four specialties cross over one at a time after a single baseline period, giving five periods.

design <- sw_design(
  clusters_per_sequence = c(30, 30, 30, 30),
  crossover_period = c(2, 3, 4, 5),
  n_periods = 5,
  sequence_names = c("Cardiology", "Internal Medicine",
                     "Family Medicine", "Neurology")
)
design
#> <sw_design>
#>   4 sequences, 5 periods, 120 clusters total
#>   clusters per sequence: 30, 30, 30, 30 
#>   treatment schedule (rows = sequences, cols = periods):
#>                   Period 1 Period 2 Period 3 Period 4 Period 5
#> Cardiology               0        1        1        1        1
#> Internal Medicine        0        0        1        1        1
#> Family Medicine          0        0        0        1        1
#> Neurology                0        0        0        0        1

The assumptions carry the specialty-specific baseline rates estimated above, the target effect, and the ICC.

assumptions <- sw_assumptions(
  baseline_prob = c(0.06, 0.04, 0.03, 0.02),
  treatment_or = 1.5,
  icc = 0.05,
  n_per_cluster_period = 20
)
assumptions
#> <sw_assumptions>
#>   outcome: binary 
#>   treatment: OR = 1.500 (log-odds 0.405)
#>   cluster_sd = 0.416 (icc = 0.050)
#>   baseline (prob scale): 0.060, 0.040, 0.030, 0.020

power <- power_swcrt(design, assumptions, nsim = 1000, seed = 1, n_cores = 4)
summary(power)

Sensitivity

A single power number is rarely what a protocol needs. The effect size and the ICC are both assumptions, so vary them.

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)
  ),
  nsim = 500, n_cores = 4, seed = 1
)
summary(grid)
plot(grid)

A note on the analysis model

fit_stepwedge_model() does not adjust for sequence by default. In a stepped-wedge trial each cluster sits in exactly one sequence for the whole study, so sequence is a cluster-level attribute that the cluster random intercept already absorbs. Adding factor(sequence) as a fixed effect is close to redundant, and in this design it measurably lowers power and raises the rate of singular fits. The baseline rates still differ by specialty – that information enters through baseline_prob, where it belongs, rather than through an extra fixed effect in the analysis.