Model-Assisted Analysis of Covariance for Stepped Wedge Cluster Randomized Experiments

A tutorial using the swcrtcov source code

1 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:

ls()  # functions made available by swcrtcov.R

## [1] "swcrt_ate"            "swcrt_fit"            "swcrt_plugin_cov"     "swcrt_run_simulation"
## [5] "swcrt_simulate"       "swcrt_treat"          "swcrt_true_ate"       "swcrt_vcovCR0"

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.

2 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.

3 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().

4 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.

dat <- swcrt_simulate(I = 30, J = 7, rho0 = 0.1, seed = 2024)
dim(dat)

## [1] 14041    21

head(dat[, c("cluster", "period", "schedule", "treatment",
             "x1", "x2", "y_obs", "w1", "w2", "w3")], 5)

cluster	period	schedule	treatment	x1	x2	y_obs	w1	w2	w3
1	2	7	0	-0.5666667	0.1448924	-2.3869946	1	0.0005777	0.0217391
1	2	7	0	-0.5666667	-0.4338905	0.1295666	1	0.0005777	0.0217391
1	2	7	0	-0.5666667	-0.6267966	1.5701729	1	0.0005777	0.0217391
1	2	7	0	-0.5666667	-0.1671781	-0.8758256	1	0.0005777	0.0217391
1	2	7	0	-0.5666667	-1.3046724	6.0011129	1	0.0005777	0.0217391

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:

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).")

Treatment status by cluster and rollout period (1 = treated).

cluster	schedule	period_2	period_3	period_4	period_5	period_6
5	2	1	1	1	1	1
9	2	1	1	1	1	1
10	2	1	1	1	1	1
19	2	1	1	1	1	1
23	2	1	1	1	1	1
14	3	0	1	1	1	1
15	3	0	1	1	1	1
20	3	0	1	1	1	1
22	3	0	1	1	1	1
24	3	0	1	1	1	1
6	4	0	0	1	1	1
12	4	0	0	1	1	1

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.

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)

##    ind period   cell 
##  0.684  0.588  0.492

5 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.

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.")

Point estimates and standard errors for one simulated trial.

Estimand	Estimator	Estimate	SE (DB)	SE (CRSE)
ind	unadjusted	0.751	0.167	0.150
ind	ANCOVA1	0.650	0.142	0.126
ind	ANCOVA2	0.640	0.137	0.121
ind	ANCOVA3	0.670	0.127	0.115
ind	ANCOVA4	0.676	0.099	0.088
period	unadjusted	0.633	0.165	0.148
period	ANCOVA1	0.540	0.149	0.134
period	ANCOVA2	0.529	0.141	0.127
period	ANCOVA3	0.553	0.131	0.120
period	ANCOVA4	0.557	0.104	0.091
cell	unadjusted	0.505	0.164	0.146
cell	ANCOVA1	0.431	0.146	0.132
cell	ANCOVA2	0.420	0.144	0.132
cell	ANCOVA3	0.446	0.128	0.116
cell	ANCOVA4	0.475	0.108	0.100

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:

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)

## lower   est upper 
## 0.421 0.670 0.919

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.

swcrt_ate(dat, weights = c(ind = "w1"),
          methods = c("unadjusted", "ANCOVA1", "ANCOVA3"),
          covariates = c("x1", "x2", "x3", "x4", "x5"))

estimand	method	estimate	se_db	se_crse
ind	unadjusted	0.7511653	0.1673459	0.1496787
ind	ANCOVA1	0.6450304	0.1391762	0.1238342
ind	ANCOVA3	0.6338824	0.1221130	0.1116755

6 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.

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).

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."))
}

6.1 Individual-average treatment effect

pretty_table(sim_all, "ind", "individual-average treatment effect")

Results for the individual-average treatment effect (\(\tau^{ind}\)). 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.

I	Estimator	BIAS	RMSE	ESE	ASE.DB	ASE.CRSE	Cov.DB	Cov.CRSE
30	UN	0.010	0.142	0.152	0.163	0.146	0.97	0.94
	AN I	0.009	0.083	0.102	0.130	0.117	0.99	0.99
	AN II	0.011	0.096	0.110	0.130	0.117	0.98	0.97
	AN III	0.007	0.081	0.099	0.109	0.098	0.98	0.98
	AN IV	0.002	0.086	0.103	0.097	0.085	0.97	0.94
60	UN	0.009	0.091	0.113	0.112	0.106	0.98	0.97
	AN I	0.013	0.060	0.091	0.088	0.084	1.00	0.99
	AN II	0.017	0.064	0.095	0.090	0.085	1.00	1.00
	AN III	0.008	0.059	0.092	0.074	0.071	0.98	0.98
	AN IV	0.006	0.061	0.095	0.071	0.067	0.98	0.97

6.2 Period-average treatment effect

pretty_table(sim_all, "period", "period-average treatment effect")

Results for the period-average treatment effect (\(\tau^{period}\)). 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.

I	Estimator	BIAS	RMSE	ESE	ASE.DB	ASE.CRSE	Cov.DB	Cov.CRSE
30	UN	0.005	0.132	0.145	0.159	0.142	0.97	0.94
	AN I	0.002	0.082	0.104	0.128	0.115	0.99	0.99
	AN II	0.005	0.094	0.112	0.127	0.115	0.98	0.98
	AN III	0.000	0.078	0.100	0.106	0.095	1.00	0.97
	AN IV	-0.015	0.084	0.103	0.095	0.083	0.96	0.93
60	UN	-0.001	0.088	0.109	0.108	0.103	0.99	0.98
	AN I	0.006	0.055	0.086	0.087	0.082	1.00	1.00
	AN II	0.011	0.059	0.089	0.089	0.084	1.00	1.00
	AN III	-0.001	0.056	0.087	0.072	0.069	0.99	0.99
	AN IV	-0.009	0.059	0.091	0.069	0.065	0.99	0.97

6.3 Cell-average treatment effect

pretty_table(sim_all, "cell", "cell-average treatment effect")

Results for the cell-average treatment effect (\(\tau^{cell}\)). 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.

I	Estimator	BIAS	RMSE	ESE	ASE.DB	ASE.CRSE	Cov.DB	Cov.CRSE
30	UN	0.014	0.132	0.149	0.160	0.143	0.96	0.95
	AN I	0.015	0.086	0.110	0.129	0.116	1.00	0.99
	AN II	0.022	0.097	0.118	0.129	0.117	0.98	0.97
	AN III	0.016	0.080	0.105	0.108	0.097	0.99	0.98
	AN IV	0.026	0.086	0.107	0.097	0.085	0.96	0.92
60	UN	0.008	0.085	0.107	0.108	0.102	1.00	0.97
	AN I	0.013	0.057	0.088	0.087	0.083	0.99	0.99
	AN II	0.022	0.060	0.091	0.089	0.084	0.99	0.99
	AN III	0.008	0.056	0.088	0.073	0.069	0.99	0.99
	AN IV	0.013	0.058	0.091	0.070	0.066	0.99	0.97

7 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.

8 Reproducibility

sessionInfo()

## R version 4.6.0 (2026-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8        LC_COLLATE=C.UTF-8    
##  [5] LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8    LC_PAPER=C.UTF-8       LC_NAME=C             
##  [9] LC_ADDRESS=C           LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39    R6_2.6.1         codetools_0.2-20 fastmap_1.2.0    xfun_0.59       
##  [6] cachem_1.1.0     knitr_1.51       htmltools_0.5.9  rmarkdown_2.31   lifecycle_1.0.5 
## [11] cli_3.6.6        sass_0.4.10      jquerylib_0.1.4  compiler_4.6.0   tools_4.6.0     
## [16] evaluate_1.0.5   bslib_0.11.0     yaml_2.3.12      rlang_1.2.0      jsonlite_2.0.0

9 Reference

Chen, X. and Li, F. (2025). Model-assisted analysis of covariance estimators for stepped wedge cluster randomized experiments. Scandinavian Journal of Statistics.