--- title: "Model-Assisted Analysis of Covariance for Stepped Wedge Cluster Randomized Experiments" subtitle: "A tutorial using the `swcrtcov` source code" author: "Companion tutorial to Chen and Li (2025, *Scandinavian Journal of Statistics*)" date: "`r format(Sys.Date(), '%B %d, %Y')`" output: html_document: toc: true toc_float: true toc_depth: 3 number_sections: true theme: flatly highlight: tango df_print: kable --- ```{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, fig.align = "center", cache = TRUE ) options(width = 100) ## Load the source code that implements every calculation in this tutorial. source("../R/swcrtcov.R") set.seed(2024) ``` # Overview Stepped wedge cluster randomized experiments (SW-CREs) are a class of unidirectional crossover designs in which all clusters begin under the usual-care condition and the intervention is rolled out across the follow-up periods in a staggered fashion until every cluster is exposed. Although these designs have become popular in public health, medicine, and implementation science, clearly defined estimands and robust, estimand-aligned methods of analysis remain comparatively underdeveloped. This tutorial illustrates the model-assisted analysis of covariance (ANCOVA) framework of Chen and Li (2025) using the accompanying source code `swcrtcov.R`. The framework (i) defines a class of weighted average treatment effect (WATE) estimands that acknowledge the multilevel data structure of a stepped wedge design, (ii) introduces four ANCOVA working models that achieve covariate adjustment while remaining consistent even under model misspecification, and (iii) provides two complementary variance estimators---a design-based (DB) plug-in estimator motivated by the staggered rollout randomization scheme, and the cluster-robust standard error (CRSE). Throughout, we generate a single stepped wedge data set, target the three interpretable estimands, compare the five estimators, and then study their finite-sample operating characteristics through a small simulation. The source code depends only on base R: ```{r show-source} ls() # functions made available by swcrtcov.R ``` The user-facing functions cover the entire workflow---data generation, estimation, variance, and Monte Carlo evaluation---as summarized below. | Function | Purpose | Key arguments | Returns | |:--|:--|:--|:--| | `swcrt_treat()` | Monotone 0/1 treatment indicator across periods for a cluster that first adopts at a given time. | `j` (adoption time), `J` | Length-`J` integer vector. | | `swcrt_simulate()` | Simulate one individual-level balanced stepped wedge data set with informative cluster-period sizes, two structured covariates plus auxiliary noise covariates, known potential outcomes, and the three estimand weight columns. | `I`, `J`, `rho0` (ICC), `seed` | Analysis-period data frame (periods $2,\ldots,J-1$). | | `swcrt_true_ate()` | Compute the true weighted ATE from the known potential outcomes (for simulation benchmarking). | `data`, `weight`, `y1`, `y0` | Numeric scalar. | | `swcrt_plugin_cov()` | Design-based (DB) plug-in covariance of the period-specific treatment effects under the staggered rollout scheme. | `data`, `residual`, `weights`, `I`, `J` | $(J-2)\times(J-2)$ covariance matrix. | | `swcrt_vcovCR0()` | Cluster-robust (CR0) sandwich covariance for a weighted `lm`; reproduces `clubSandwich::vcovCR(., "CR0")`. | `fit` (an `lm`), `cluster` | $p\times p$ covariance matrix. | | `swcrt_fit()` | Fit **one** estimator for **one** estimand by weighted least squares and return the point estimate with both standard errors. | `data`, `outcome`, `weight`, `method`, `covariates` | List: `estimate`, `se_db`, `se_crse`, plus the fitted model. | | `swcrt_ate()` | Convenience wrapper that runs **all** estimators across **all** estimands. | `data`, `outcome`, `weights`, `methods`, `covariates` | Tidy data frame, one row per (estimand, estimator). | | `swcrt_run_simulation()` | Monte Carlo driver: repeatedly generates data, records the truth, and estimates every combination, summarizing bias, RMSE, average SEs, and CI coverage. | `n_rep`, `gen_fun`, `weights`, `methods`, `covariates`, `seed` | List: `estimates` (raw draws) and `summary` (metrics). | The `method` argument of `swcrt_fit()` and `swcrt_ate()` accepts `"unadjusted"`, `"ANCOVA1"`, `"ANCOVA2"`, `"ANCOVA3"`, and `"ANCOVA4"`, corresponding to the five working models. Helper functions prefixed with a dot (e.g. `.swcrt_formula()`, `.swcrt_design()`) are internal and are not intended to be called directly. # Notation and estimands We adopt a finite population framework in which all potential outcomes are fixed quantities and the only source of variability is the randomization distribution. Consider a stepped wedge design with $J$ periods indexed $j = 1,\ldots,J$. The first period (all clusters under control) and the last period (all clusters under treatment) carry no treatment contrast and are excluded from the analysis owing to the absence of positivity; the analysis therefore proceeds over the $J-2$ rollout periods $2,\ldots,J-1$. For cluster $i$, period $j$, and individual $k$, let $Y_{ijk}(z)$ denote the potential outcome under treatment status $z \in \{0,1\}$, with cluster-period size $N_{ij}$. The general class of weighted average treatment effect estimands during the rollout is $$ \tau \;=\; \frac{\sum_j \sum_i \sum_k w_{ijk}\{Y_{ijk}(1) - Y_{ijk}(0)\}} {\sum_j \sum_i \sum_k w_{ijk}}, $$ where $w_{ijk} \ge 0$ is an individual-specific weight. Different choices of weight yield different interpretable members of this family. We focus on the three estimands summarized below, which weigh the population of all individuals, all rollout periods, or all cluster-period cells, respectively. | Estimand | Weight $w_{ijk}$ | Weight column | Interpretation | |:--|:--:|:--:|:--| | $\tau^{\text{ind}}$ (individual-average) | $1$ | `w1` | Average of individual-level counterfactual contrasts over all individuals during the rollout; equal weight to each individual. | | $\tau^{\text{period}}$ (period-average) | $N_j^{-1}$ | `w2` | Average of individual-level contrasts per period, further averaged over rollout periods; equal weight to each period. | | $\tau^{\text{cell}}$ (cell-average) | $N_{ij}^{-1}$ | `w3` | Average of cluster-period cell-level mean contrasts over all cells during the rollout; equal weight to each cluster-period cell. | When cluster-period sizes are heterogeneous, these estimands need not coincide, especially under informative cluster size---that is, when the cluster-period size $N_{ij}$ is associated with the within-cluster counterfactual contrast. They coincide only when the sizes are homogeneous or the treatment effect is constant across cells. The choice among them should be guided by the study context and objectives. # The ANCOVA estimators To target the estimands above, we fit a working ANCOVA model by weighted least squares (WLS), using the weights $w_{ijk}$ that correspond to the estimand of interest. Writing $Z_{ij}$ for the treatment indicator and $\widetilde{\mathbf{X}}_{ijk}$ for the period-mean-centered covariate vector, the four model formulations differ only in whether period-specific covariate main effects and treatment-by-covariate interactions are included: | Estimator | Period-specific covariate effects | Treatment-by-covariate interactions | |:--|:--:|:--:| | **Unadjusted** | --- | --- | | **ANCOVA I** | $\times$ | $\times$ | | **ANCOVA II** | $\checkmark$ | $\times$ | | **ANCOVA III** | $\times$ | $\checkmark$ | | **ANCOVA IV** | $\checkmark$ | $\checkmark$ | Each estimator is *model-assisted*: the period-specific effect $\widehat{\tau}_j$, and hence the overall estimator $\widehat{\tau} = \sum_j w_j \widehat{\Delta}_j / \sum_j w_j$, remains consistent for the target estimand even when the working model is incorrectly specified. The unadjusted estimator is the special case of ANCOVA I with no covariate adjustment. The variance of $\widehat{\tau}$ is obtained from the covariance matrix of the period-specific effects, sandwiched by the (known) normalized weight vector $\boldsymbol{\varpi}$. We consider two estimators of this covariance: * the **design-based (DB)** plug-in estimator, derived from the randomization-based variance expression under the staggered rollout scheme and known to be conservative because a non-negative between-sequence term is generally inestimable; and * the **cluster-robust standard error (CRSE)** of Liang and Zeger (1986), which assumes working independence between clusters but permits arbitrary within-cluster correlation. In `swcrtcov.R`, the DB covariance is computed by `swcrt_plugin_cov()` and the CRSE (type CR0) by `swcrt_vcovCR0()`; both are assembled, together with the point estimate, inside `swcrt_fit()`. # Generating a stepped wedge data set The function `swcrt_simulate()` generates one individual-level data set from a balanced stepped wedge design. We set the total number of periods to $J = 7$, which corresponds to five rollout periods together with one pre-rollout and one post-rollout period, and use $I = 30$ clusters so that $I/(J-1) = 5$ previously untreated clusters adopt the intervention at each rollout time. The cluster-period sizes increase with the period and are deliberately informative; the potential outcomes are generated to be nonlinear in the covariates so as to probe the model-robustness of the estimators. ```{r generate} dat <- swcrt_simulate(I = 30, J = 7, rho0 = 0.1, seed = 2024) dim(dat) head(dat[, c("cluster", "period", "schedule", "treatment", "x1", "x2", "y_obs", "w1", "w2", "w3")], 5) ``` The returned data frame is already restricted to the rollout periods $2,\ldots,6$ and contains the period indicator columns `period_2`,...,`period_6`, the adoption time `schedule`, the treatment indicator `treatment`, the centered covariates, the observed outcome `y_obs`, the (unobservable) potential outcomes `y0` and `y1`, and the three weight columns. The staggered rollout is most transparent in the cluster-by-period treatment matrix, in which clusters are ordered by their adoption time: ```{r design-matrix} cp <- unique(dat[, c("cluster", "schedule", "period", "treatment")]) wide <- reshape(cp[, c("cluster", "schedule", "period", "treatment")], idvar = c("cluster", "schedule"), timevar = "period", direction = "wide") wide <- wide[order(wide$schedule, wide$cluster), ] colnames(wide) <- sub("treatment.", "period_", colnames(wide), fixed = TRUE) knitr::kable(head(wide, 12), row.names = FALSE, caption = "Treatment status by cluster and rollout period (1 = treated).") ``` Because the potential outcomes are known in this simulated trial, we can compute the true value of each estimand directly from the individual-level counterfactual contrasts using `swcrt_true_ate()`. The three estimands differ in value, a direct consequence of the informative cluster-period size built into the data-generating process. ```{r true-values} true_vals <- c( ind = swcrt_true_ate(dat, "w1"), period = swcrt_true_ate(dat, "w2"), cell = swcrt_true_ate(dat, "w3") ) round(true_vals, 3) ``` # Estimating the three estimands on a single trial A single call to `swcrt_ate()` fits all five estimators for each of the three estimands and returns the point estimate together with the design-based and cluster-robust standard errors. Following the primary analysis in Chen and Li (2025), we adjust for the two baseline covariates $\mathbf{X}_{ijk} = (X_{ij1}, X_{ijk2})$, supplied here as `x1` and `x2`. ```{r single-fit} fit_tab <- swcrt_ate(dat, outcome = "y_obs", weights = c(ind = "w1", period = "w2", cell = "w3"), covariates = c("x1", "x2")) fit_tab[, c("estimate", "se_db", "se_crse")] <- round(fit_tab[, c("estimate", "se_db", "se_crse")], 3) knitr::kable(fit_tab, col.names = c("Estimand", "Estimator", "Estimate", "SE (DB)", "SE (CRSE)"), caption = "Point estimates and standard errors for one simulated trial.") ``` Two patterns, both anticipated by the theory, are already visible. First, within each estimand the covariate-adjusted estimators return standard errors that are smaller than those of the unadjusted estimator, illustrating the efficiency gain from covariate adjustment. Second, the design-based standard errors are systematically at least as large as the cluster-robust ones, consistent with the conservativeness of the DB estimator. A 95% Wald confidence interval for any row is formed in the usual way. For the individual-average treatment effect estimated by ANCOVA III, using the design-based standard error: ```{r ci-example} row <- subset(swcrt_ate(dat, weights = c(ind = "w1"), methods = "ANCOVA3", covariates = c("x1", "x2"))) ci <- c(lower = row$estimate - 1.96 * row$se_db, est = row$estimate, upper = row$estimate + 1.96 * row$se_db) round(ci, 3) ``` The framework places no restriction on the covariate set. As discussed in Chen and Li (2025), the cluster-period size $N_{ij}$ may itself be adjusted for as an additional covariate; here we illustrate by appending the three auxiliary covariates available in the generated data. ```{r more-covariates} swcrt_ate(dat, weights = c(ind = "w1"), methods = c("unadjusted", "ANCOVA1", "ANCOVA3"), covariates = c("x1", "x2", "x3", "x4", "x5")) ``` # Finite-sample operating characteristics A single trial cannot speak to bias or coverage. We therefore evaluate the estimators over repeated stepped wedge experiments using `swcrt_run_simulation()`, which for each replication generates a new data set, records the true estimand from the potential outcomes, and estimates every method-by-estimand combination. We summarize performance through the relative bias (BIAS), the root mean square error (RMSE), the empirical standard error (ESE), the average standard errors (ASE) from the DB and CRSE approaches, and the empirical coverage of 95% confidence intervals. To keep the runtime modest we use 100 replications at two cluster sizes, $I \in \{30, 60\}$; the original study used 1{,}000 replications and $I \in \{18, 30, 60, 120\}$. The qualitative conclusions are unchanged. ```{r run-sim} run_one <- function(I, n_rep = 100) { sim <- swcrt_run_simulation( n_rep = n_rep, gen_fun = function() swcrt_simulate(I, 7, rho0 = 0.1), weights = c(ind = "w1", period = "w2", cell = "w3"), covariates = c("x1", "x2"), seed = 2024, progress = FALSE) cbind(I = I, sim$summary) } sim30 <- run_one(30) sim60 <- run_one(60) sim_all <- rbind(sim30, sim60) ``` The helper below renders the results in the layout of the simulation tables in Chen and Li (2025). ```{r format-sim} pretty_table <- function(df, estimand_label, estimand_name) { d <- df[df$estimand == estimand_label, ] meth_map <- c(unadjusted = "UN", ANCOVA1 = "AN I", ANCOVA2 = "AN II", ANCOVA3 = "AN III", ANCOVA4 = "AN IV") out <- data.frame( I = d$I, Estimator = meth_map[d$method], BIAS = round(d$rel_bias, 3), RMSE = round(d$rmse, 3), ESE = round(d$emp_sd, 3), `ASE.DB` = round(d$mean_se_db, 3), `ASE.CRSE`= round(d$mean_se_crse, 3), `Cov.DB` = round(d$cover_db, 3), `Cov.CRSE`= round(d$cover_crse, 3), check.names = FALSE, row.names = NULL) out$I[duplicated(out$I)] <- "" knitr::kable(out, row.names = FALSE, caption = paste0( "Results for the ", estimand_name, " ($\\", "tau^{", estimand_label, "}$). BIAS = relative bias; RMSE = root mean squared error; ", "ESE = empirical standard error; ASE = average standard error ", "(DB = design-based plug-in, CRSE = cluster-robust); ", "Coverage = empirical coverage of 95% CIs over 200 replications.")) } ``` ## Individual-average treatment effect ```{r tab-ind} pretty_table(sim_all, "ind", "individual-average treatment effect") ``` ## Period-average treatment effect ```{r tab-period} pretty_table(sim_all, "period", "period-average treatment effect") ``` ## Cell-average treatment effect ```{r tab-cell} pretty_table(sim_all, "cell", "cell-average treatment effect") ``` # Interpretation and practical recommendations The simulation reproduces the central messages of Chen and Li (2025). For every estimand, the relative bias and RMSE decrease as the number of clusters $I$ grows, corroborating the consistency of all five estimators. The covariate-adjusted estimators attain smaller RMSE than the unadjusted estimator, confirming that adjusting for baseline covariates improves efficiency, although the most efficient working model depends on the data-generating process. The two variance estimators exhibit complementary finite-sample behavior. For the unadjusted, ANCOVA I, and ANCOVA II estimators, the design-based approach tends to be conservative when the number of clusters is small, whereas the cluster-robust approach yields average standard errors closer to the empirical standard error. For ANCOVA III, which includes treatment-by-covariate interactions but no period-specific covariate main effects, the design-based approach yields average standard errors close to the empirical standard error with near-nominal coverage across cluster sizes. For ANCOVA IV, which additionally includes period-specific interactions and therefore estimates the largest number of parameters, the design-based variance can be slightly underestimated when $I$ is small but improves as $I$ increases. These observations lead to a consistent practical recommendation. When a sufficient number of clusters are available, ANCOVA estimators with treatment-by-covariate interactions (ANCOVA III or IV) are preferred for their efficiency. When only a limited number of clusters are present, the ANCOVA III estimator---with treatment-by-covariate interactions but without period-specific covariate main effects---offers the best balance of efficiency and finite-sample stability, and the design-based standard error is the more reliable choice for ANCOVA III and IV, whereas the cluster-robust standard error is preferable for the unadjusted, ANCOVA I, and ANCOVA II estimators. # Reproducibility ```{r session} sessionInfo() ``` # Reference Chen, X. and Li, F. (2025). Model-assisted analysis of covariance estimators for stepped wedge cluster randomized experiments. *Scandinavian Journal of Statistics*.