Study the operating characteristics of stepped wedge cluster-randomised trials using the {simswjm} R package

simswjm-example.R

# Note:
# This code is released under the MIT license.
# Copyright (c) 2025, Alessandro Gasparini
# More details on the license here: https://choosealicense.com/licenses/mit

# Load packages
library(tidyverse)
library(simswjm) # from: https://github.com/RedDoorAnalytics/simswjm
library(lme4)
library(broom)
library(broom.mixed)
library(ragg)
library(ggtext)
library(mgcv)

# Set random seed, for reproducibility
set.seed(1138746)

# Function to fit a mixed effects model to either complete or observed data
fitmod <- function(data, model) {
  model <- match.arg(arg = model, choices = c("observed", "target"))
  if (model == "observed") {
    fit <- tryCatch(
      {
        lmer(yobs ~ 0 + factor(j) + factor(x) + (1 | id) + (1 | i), data = data, REML = TRUE, control = lmerControl(restart_edge = TRUE))
      },
      error = function(e) NULL
    )
  } else {
    fit <- tryCatch(
      {
        lmer(y ~ 0 + factor(j) + factor(x) + (1 | id) + (1 | i), data = data, REML = TRUE, control = lmerControl(restart_edge = TRUE))
      },
      error = function(e) NULL
    )
  }
  if (is.null(fit)) {
    data.frame(model = model, estimate = NA, std.error = NA)
  } else {
    tidy(fit) |>
      filter(term == "factor(x)1") |>
      mutate(model = model) |>
      select(model, estimate, std.error)
  }
}

# Function to simulate a trial and analyse the data under both scenarios
simfun <- function(b, i, n_min = 15, n_max = 150) {
  n_per_cluster <- runif(n = 1, min = n_min, max = n_max)
  n_per_cluster <- round(n_per_cluster)
  simdt <- swtrial_inf(
    repn = 1,
    k = n_per_cluster,
    i = i,
    j = 5,
    intervention_seq = 4,
    deltas = rep(5, 4),
    betas = rep(30, 5),
    family = "gaussian",
    nu = 0.0,
    lambda = 0.1,
    p = 1.0,
    sigma_epsilon = sqrt(40),
    sigma_alpha = sqrt(2),
    sigma_gamma = 0,
    sigma_phi = sqrt(55),
    omega1 = 0,
    omega3 = 0
  )
  res_obs <- fitmod(data = simdt, model = "obs")
  res_target <- fitmod(data = simdt, model = "target")
  bind_rows(res_obs, res_target) |>
    mutate(b = b, n = n_per_cluster, i = i)
}

# Simulating N = 10,000 trials
res <- map(.x = seq(10000), .f = function(b) {
  simfun(b = b, i = 8) # 8 clusters per trial, i.e. 2 per intervention sequence
}, .progress = TRUE)

# Combine the results and calculate the SE inflation in the presence of dropout
# The SE inflation is the ratio of the SE with dropout over the SE with complete data
resb <- list_rbind(res) |>
  select(-estimate) |>
  pivot_wider(names_from = "model", values_from = "std.error") |>
  mutate(ratio = observed / target)

# Plot it
pl <- ggplot(resb, aes(x = n, y = ratio)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_point(size = 1) +
  geom_smooth(method = "gam", formula = y ~ s(x), color = "#6667AB") +
  scale_x_continuous(breaks = seq(15, 150, by = 15)) +
  labs(
    x = "Participants per Cluster",
    y = "SE Inflation",
    caption = "Each dot represents a simulated trial. The <span style = 'color:#6667AB;'>purple line</span> is a smoother using thin plate regression splines."
  ) +
  theme_bw(base_size = 12, base_family = "Atkinson Hyperlegible") +
  theme(
    plot.caption = element_textbox_simple(
      lineheight = 1,
      padding = margin(5, 0, 0, 0)
    )
  )
pl
ggsave(plot = pl, filename = "blog/2025/02/23/simswjm.png", device = agg_png, dpi = 600, width = 5, height = 3)