--- title: "A Tutorial for Model-Robust Standardization in Stepped-Wedge Cluster-Randomized Trials Using `MRStdLCRT`" date: "`r format(Sys.Date(), '%B %d, %Y')`" output: html_document: toc: true toc_float: true toc_depth: 3 number_sections: true theme: readable highlight: tango pdf_document: toc: true toc_depth: 3 number_sections: true keep_tex: false latex_engine: xelatex geometry: margin=1in fontsize: 11pt urlcolor: blue linkcolor: black --- ```{r setup, include = FALSE} knitr::opts_chunk$set( echo = TRUE, warning = FALSE, message = FALSE, comment = "#>", tidy = FALSE, fig.align = "center", fig.width = 10, fig.height = 6, out.width = "100%" ) options(width = 78) ``` # Introduction This tutorial presents a workflow for applying model-robust standardization to stepped-wedge cluster-randomized trials using the `MRStdLCRT` R package. The tutorial introduces the required data structure, the available estimands, the model-fitting syntax, and the interpretation of output returned by `mrstdlcrt_fit()`. The tutorial uses the example datasets distributed with the `MRStdLCRT` package: - `sw_c`, an illustrative stepped-wedge dataset with a continuous outcome; - `sw_b`, an illustrative stepped-wedge dataset with a binary outcome. The numerical output produced by this tutorial will be interpreted as instructional output from example data. A stepped-wedge cluster-randomized trial is a longitudinal cluster-randomized design in which clusters transition from control to intervention according to a pre-specified rollout schedule. Outcomes are observed over multiple periods, and the timing of intervention adoption varies across clusters. A conventional regression coefficient for treatment may not directly correspond to a marginal causal estimand of scientific interest. Model-robust standardization addresses this issue by defining explicit marginal estimands and estimating them through a combination of nonparametric standardization and working-model-based augmentation. The `MRStdLCRT` package returns four treatment-effect estimands. These estimands differ according to two design choices: - whether aggregation is conducted horizontally or vertically over the stepped-wedge rollout; - whether the target population weights individuals or clusters. ## Estimands reported by `MRStdLCRT` | Estimand | Description | |----------|-------------| | `h-iATE` | Horizontal individual-average treatment effect. The estimand targets the average treatment effect for individuals across the rollout periods. | | `h-cATE` | Horizontal cluster-average treatment effect. The estimand gives equal weight to clusters before averaging across the rollout periods. | | `v-iATE` | Vertical individual-average treatment effect. The estimand first forms period-specific individual-average contrasts and then averages across periods. | | `v-cATE` | Vertical cluster-period-average treatment effect. The estimand gives equal weight to eligible cluster-period cells in the period-level aggregation. | Individual-weighted and cluster-weighted estimands may differ when cluster-period sizes vary. Horizontal and vertical estimands may differ when the treatment effect varies over time or when the distribution of cluster-period sizes changes over the rollout. ## Estimators and inference For each estimand, `MRStdLCRT` reports two estimators: - an unadjusted estimator based on observed cluster-period outcome summaries; - an adjusted model-robust standardized estimator that uses a working outcome model for augmentation and covariate adjustment. The working model is not itself the final estimand. Instead, it is used to improve precision and to adjust for pre-specified covariates. The resulting treatment-effect estimates remain marginal standardized estimands. Variance is evaluated using delete-one-cluster jackknife inference. Confidence intervals are reported using a reference distribution with degrees of freedom determined by the number of clusters. --- # Software requirements The following packages are used throughout the tutorial. ```{r installation, eval = FALSE} install.packages(c("MRStdLCRT", "dplyr", "tidyr", "ggplot2", "knitr")) ``` ```{r packages} required_packages <- c("MRStdLCRT", "dplyr", "tidyr", "ggplot2", "knitr") missing_packages <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)] if (length(missing_packages) > 0) { stop( "Install the following package(s) before knitting this tutorial: ", paste(missing_packages, collapse = ", ") ) } library(MRStdLCRT) library(dplyr) library(tidyr) library(ggplot2) library(knitr) ``` --- # Data structure required for analysis `MRStdLCRT` expects a long-format analysis dataset. In the example datasets used here, each row represents one individual observed within a specific cluster and period. | Variable role | Example variable | Description | |---------------|------------------|-------------| | Cluster identifier | `cluster` | Randomization unit, such as a site, clinic, school, or hospital. | | Period identifier | `period` | Ordered time period in the stepped-wedge rollout. | | Individual identifier | `id` | Identifier for each individual record. | | Treatment assignment | `trt` | Binary cluster-period treatment indicator, coded 0 for control and 1 for intervention. | | Outcome | `y` | Continuous in `sw_c`; binary in `sw_b`. | | Covariates | `x1`, `x2` | Auxiliary covariates used for illustration. | The treatment variable must be constant within each cluster-period cell. This requirement follows from the cluster-period assignment structure in a stepped-wedge design. > **Analytic principle.** The treatment variable supplied to `mrstdlcrt_fit()` should represent cluster-period treatment assignment. If a dataset records individual-level intervention receipt, that variable should not be used as the cluster-period treatment indicator unless it is constant within every cluster-period cell. --- # Helper functions The following helper functions can be used to inspect the design, summarize the outcome, and create compact reporting tables. They are not required by the `MRStdLCRT` package, but can provide useful checks before fitting a model. ```{r helper-functions} check_cluster_period_treatment <- function(dat) { violations <- dat %>% group_by(cluster, period) %>% summarise(n_treatment_values = n_distinct(trt), .groups = "drop") %>% filter(n_treatment_values > 1) if (nrow(violations) == 0) { cat("PASS: treatment assignment is constant within every cluster-period cell.\n") } else { cat("FAIL: at least one cluster-period cell contains more than one treatment value.\n") print(violations) } invisible(violations) } period_treatment_table <- function(dat) { dat %>% distinct(cluster, period, trt) %>% group_by(period) %>% summarise( control_clusters = sum(trt == 0), treated_clusters = sum(trt == 1), total_clusters = n(), informative = control_clusters > 0 & treated_clusters > 0, .groups = "drop" ) } summarize_outcome_by_treatment <- function(dat, outcome = "y") { dat %>% group_by(trt) %>% summarise( N = n(), Mean = round(mean(.data[[outcome]], na.rm = TRUE), 3), SD = round(sd(.data[[outcome]], na.rm = TRUE), 3), Median = round(median(.data[[outcome]], na.rm = TRUE), 3), Q1 = round(quantile(.data[[outcome]], 0.25, na.rm = TRUE), 3), Q3 = round(quantile(.data[[outcome]], 0.75, na.rm = TRUE), 3), .groups = "drop" ) %>% mutate(Treatment = ifelse(trt == 1, "Intervention", "Control")) %>% select(Treatment, N, Mean, SD, Median, Q1, Q3) } cluster_period_size_summary <- function(dat) { dat %>% count(cluster, period, trt, name = "n_individuals") %>% summarise( cluster_periods = n(), minimum = min(n_individuals), first_quartile = quantile(n_individuals, 0.25), median = median(n_individuals), mean = mean(n_individuals), third_quartile = quantile(n_individuals, 0.75), maximum = max(n_individuals) ) } extract_estimates <- function(fit, model_label) { if (is.null(fit$estimates)) { stop("The fitted object does not contain an `estimates` element.") } as_tibble(fit$estimates) %>% mutate(working_model = model_label, .before = 1) } ``` --- # Example 1: Continuous outcome with a GEE working model This section analyzes the continuous-outcome example dataset `sw_c`. The working model is a Gaussian generalized estimating equation with period fixed effects and an exchangeable working correlation structure. ## Load and prepare the data ```{r swc-load} data("sw_c", package = "MRStdLCRT") sw_c <- as_tibble(sw_c) %>% mutate( cluster = as.integer(cluster), period = as.integer(period), id = as.integer(id), trt = as.integer(trt), x1 = as.integer(x1), x2 = as.numeric(x2), y = as.numeric(y) ) %>% drop_na(cluster, period, id, trt, x1, x2, y) required_variables <- c("cluster", "period", "id", "trt", "x1", "x2", "y") stopifnot(all(required_variables %in% names(sw_c))) cat("Number of rows:", nrow(sw_c), "\n") cat("Number of clusters:", n_distinct(sw_c$cluster), "\n") cat("Number of periods:", n_distinct(sw_c$period), "\n") ``` The first rows show the expected long-format structure. ```{r swc-preview} head(sw_c) %>% kable() ``` ## Verify the stepped-wedge treatment structure The first diagnostic confirms that the treatment indicator is constant within each cluster-period cell. ```{r swc-treatment-check} check_cluster_period_treatment(sw_c) ``` The next table counts treated and control clusters by period. Periods with both treated and control clusters are informative for marginal treatment-effect aggregation. ```{r swc-period-treatment-table} period_treatment_table(sw_c) %>% kable() ``` In a stepped-wedge design, some edge periods may contain only control clusters or only intervention clusters. Such periods can still contribute information to the working outcome model, but they do not provide direct treated-versus-control contrasts for aggregation. ```{r swc-cluster-period-size} cluster_period_size_summary(sw_c) %>% kable(digits = 2) ``` ## Summarize the outcome descriptively The following table summarizes the continuous outcome by observed treatment assignment. This table is descriptive and does not adjust for period, cluster, or covariates. ```{r swc-descriptive-summary} summarize_outcome_by_treatment(sw_c, outcome = "y") %>% kable() ``` ## Specify the working model The working model includes treatment assignment, period indicators, and the two illustrative covariates. The treatment coefficient in this working model should not be interpreted as the final marginal estimand. ```{r swc-gee-formula} formula_sw_c_gee <- y ~ trt + factor(period) + x1 + x2 formula_sw_c_gee ``` ## Fit the model-robust standardized estimator ```{r swc-fit-gee} fit_sw_c_gee <- mrstdlcrt_fit( data = sw_c, formula = formula_sw_c_gee, cluster_id = "cluster", period = "period", trt = "trt", method = "gee", family = "gaussian", corstr = "exchangeable" ) ``` ## Review the treatment-effect estimates ```{r swc-summary-gee} summary( fit_sw_c_gee, level = 0.95, show_counts = TRUE, ics = "global", digits = 4 ) ``` ```{r swc-plot-gee} plot( fit_sw_c_gee, estimand = c("h-iATE", "h-cATE", "v-iATE", "v-cATE") ) + labs( title = "Continuous outcome example", subtitle = "Gaussian GEE working model with exchangeable working correlation", caption = "h = horizontal; v = vertical; i = individual-weighted; c = cluster-weighted" ) + theme_bw() ``` ## Interpretation The summary table reports the four estimands under both unadjusted and adjusted estimation. The adjusted estimates are standardized marginal treatment-effect estimates. They use the GEE working model for covariate adjustment and precision improvement, while retaining the marginal estimand definitions specified by the package. Differences between `h-iATE` and `h-cATE` illustrate the consequence of weighting individuals rather than clusters. Differences between `v-iATE` and `v-cATE` similarly illustrate the effect of the aggregation scheme over eligible cluster-period cells. --- # Example 2: Continuous outcome with a linear mixed-effects working model This section analyzes the same continuous-outcome dataset using a linear mixed-effects working model. The purpose is to demonstrate how the working model can be changed while preserving the same four marginal estimands. ## Specify the mixed-effects working model A random intercept for cluster is included to represent between-cluster outcome heterogeneity. ```{r swc-lmer-formula} formula_sw_c_lmer <- y ~ trt + factor(period) + x1 + x2 + (1 | cluster) formula_sw_c_lmer ``` ## Fit the estimator ```{r swc-fit-lmer} fit_sw_c_lmer <- mrstdlcrt_fit( data = sw_c, formula = formula_sw_c_lmer, cluster_id = "cluster", period = "period", trt = "trt", method = "lmer", family = "gaussian" ) ``` ## Review the treatment-effect estimates ```{r swc-summary-lmer} summary( fit_sw_c_lmer, level = 0.95, show_counts = FALSE, ics = "global", digits = 4 ) ``` ```{r swc-plot-lmer} plot( fit_sw_c_lmer, estimand = c("h-iATE", "h-cATE", "v-iATE", "v-cATE") ) + labs( title = "Continuous outcome example", subtitle = "Linear mixed-effects working model with random cluster intercept", caption = "The estimands are unchanged; the working model used for adjustment is changed." ) + theme_bw() ``` ## Compare the two continuous-outcome working models The table below places the GEE and mixed-effects results in a common display. This comparison is a working-model sensitivity assessment, not a comparison of different estimands. ```{r swc-working-model-comparison} bind_rows( extract_estimates(fit_sw_c_gee, "Gaussian GEE"), extract_estimates(fit_sw_c_lmer, "Linear mixed-effects model") ) %>% kable(digits = 4) ``` When the adjusted estimates are stable across working models, the display provides reassurance that the standardized estimates are not strongly dependent on the particular working-model specification. When the estimates differ meaningfully, the analyst should inspect covariate coding, time trends, cluster heterogeneity, convergence behavior, and the plausibility of the model structure. --- # Example 3: Binary outcome with a logistic GEE working model This section analyzes the binary-outcome example dataset `sw_b`. The working model is a logistic generalized estimating equation, and the marginal treatment-effect estimates are reported on the risk-difference scale. ## Load and prepare the data ```{r swb-load} data("sw_b", package = "MRStdLCRT") sw_b <- as_tibble(sw_b) %>% mutate( cluster = as.integer(cluster), period = as.integer(period), id = as.integer(id), trt = as.integer(trt), x1 = as.integer(x1), x2 = as.numeric(x2), y = as.integer(y) ) %>% drop_na(cluster, period, id, trt, x1, x2, y) stopifnot(all(required_variables %in% names(sw_b))) stopifnot(all(sw_b$y %in% c(0, 1))) cat("Number of rows:", nrow(sw_b), "\n") cat("Number of clusters:", n_distinct(sw_b$cluster), "\n") cat("Number of periods:", n_distinct(sw_b$period), "\n") ``` ```{r swb-preview} head(sw_b) %>% kable() ``` ## Verify the stepped-wedge treatment structure ```{r swb-treatment-check} check_cluster_period_treatment(sw_b) ``` ```{r swb-period-treatment-table} period_treatment_table(sw_b) %>% kable() ``` ## Inspect cluster-period sizes ```{r swb-cluster-period-size} cluster_period_size_summary(sw_b) %>% kable(digits = 2) ``` ## Summarize the outcome descriptively For a binary outcome, the mean of `y` is the observed event proportion. The following table is descriptive and does not adjust for the stepped-wedge design. ```{r swb-descriptive-summary} summarize_outcome_by_treatment(sw_b, outcome = "y") %>% kable() ``` ## Specify the working model ```{r swb-gee-formula} formula_sw_b_gee <- y ~ trt + factor(period) + x1 + x2 formula_sw_b_gee ``` ## Fit the estimator on the risk-difference scale ```{r swb-fit-gee-rd} fit_sw_b_gee_rd <- mrstdlcrt_fit( data = sw_b, formula = formula_sw_b_gee, cluster_id = "cluster", period = "period", trt = "trt", method = "gee", family = "binomial", corstr = "exchangeable", scale = "RD" ) ``` ## Review the treatment-effect estimates ```{r swb-summary-gee-rd} summary( fit_sw_b_gee_rd, level = 0.95, show_counts = TRUE, ics = "global", digits = 4 ) ``` ```{r swb-plot-gee-rd} plot( fit_sw_b_gee_rd, estimand = c("h-iATE", "h-cATE", "v-iATE", "v-cATE") ) + labs( title = "Binary outcome example", subtitle = "Logistic GEE working model with estimates reported as risk differences", caption = "For binary outcomes, scale = 'RD' requests marginal risk-difference estimates." ) + theme_bw() ``` ## Interpretation The risk-difference scale expresses the marginal treatment effect as an absolute difference in event probability. For example, an estimate of 0.05 would correspond to a five-percentage-point higher event probability under intervention than under control, while an estimate of -0.05 would correspond to a five-percentage-point lower event probability under intervention. Because the example is based on illustrative data, the numerical value of the estimate is used only to demonstrate interpretation and reporting. ## Optional alternative binary scales The package also allows alternative binary-outcome scales. ```{r swb-alternative-scales, eval = T} fit_sw_b_gee_rr <- mrstdlcrt_fit( data = sw_b, formula = formula_sw_b_gee, cluster_id = "cluster", period = "period", trt = "trt", method = "gee", family = "binomial", corstr = "exchangeable", scale = "RR" ) fit_sw_b_gee_or <- mrstdlcrt_fit( data = sw_b, formula = formula_sw_b_gee, cluster_id = "cluster", period = "period", trt = "trt", method = "gee", family = "binomial", corstr = "exchangeable", scale = "OR" ) summary( fit_sw_b_gee_or, level = 0.95, show_counts = TRUE, ics = "global", digits = 4 ) ``` The selected scale should correspond to the scientific question, the outcome frequency, and the reporting conventions of the relevant field. ## Working-model selection The working model should be selected before examining treatment-effect results. A typical formula includes treatment assignment, period indicators, and pre-specified covariates: ```{r working-model-template, eval = FALSE} y ~ trt + factor(period) + baseline_covariate_1 + baseline_covariate_2 ``` Possible working-model choices include the following. | Outcome type | Method | Working-model description | |--------------|--------|---------------------------| | Continuous | `gee` | Marginal Gaussian working model with a selected working correlation structure. | | Continuous | `lmer` | Linear mixed-effects working model with random effects. | | Binary | `gee` | Marginal logistic working model with risk-difference, risk-ratio, or odds-ratio scale selected through `scale`. | | Binary | `glmer` | Logistic mixed-effects working model with random effects. | The working model should be viewed as an adjustment model. The treatment-effect estimates reported by `MRStdLCRT` are standardized marginal estimands, not unstandardized regression coefficients. --- # References - Fang X, Li F. *MRStdLCRT: Model-Robust Standardization for Longitudinal Cluster-Randomized Trials.* R package. - Fang X, Wang X, Heagerty PJ, Wang B, Li F. *Model-robust standardization in stepped wedge cluster randomized trials.* - Hussey MA, Hughes JP. Design and analysis of stepped wedge cluster randomized trials. *Contemporary Clinical Trials.* 2007;28(2):182--191.