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

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

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

2 Software requirements

The following packages are used throughout the tutorial.

install.packages(c("MRStdLCRT", "dplyr", "tidyr", "ggplot2", "knitr"))

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)

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

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

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

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

5.1 Load and prepare the data

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

#> Number of rows: 10409

cat("Number of clusters:", n_distinct(sw_c$cluster), "\n")

#> Number of clusters: 30

cat("Number of periods:", n_distinct(sw_c$period), "\n")

#> Number of periods: 6

The first rows show the expected long-format structure.

head(sw_c) %>%
  kable()

cluster	period	id	x1	x2	y
1	1	1	1	0.4236987	0.9202646
1	1	2	0	0.2536911	2.3420711
1	1	3	0	0.4517473	-0.3612097
1	1	4	0	0.4314234	-1.0275083
1	1	5	1	0.1000356	1.2406739
1	1	6	0	0.1546604	0.4158517

5.2 Verify the stepped-wedge treatment structure

The first diagnostic confirms that the treatment indicator is constant within each cluster-period cell.

check_cluster_period_treatment(sw_c)

#> PASS: treatment assignment is constant within every cluster-period cell.

The next table counts treated and control clusters by period. Periods with both treated and control clusters are informative for marginal treatment-effect aggregation.

period_treatment_table(sw_c) %>%
  kable()

period	control_clusters	treated_clusters	total_clusters	informative
1	30	0	30	FALSE
2	24	6	30	TRUE
3	18	12	30	TRUE
4	12	18	30	TRUE
5	6	24	30	TRUE
6	0	30	30	FALSE

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.

cluster_period_size_summary(sw_c) %>%
  kable(digits = 2)

cluster_periods	minimum	first_quartile	median	mean	third_quartile	maximum
180	20	39.75	55	57.83	77.25	100

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

summarize_outcome_by_treatment(sw_c, outcome = "y") %>%
  kable()

Treatment	N	Mean	SD	Median	Q1	Q3
Control	4935	2.175	2.463	1.553	0.306	3.778
Intervention	5474	12.490	6.851	12.898	7.275	17.725

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

formula_sw_c_gee <- y ~ trt + factor(period) + x1 + x2
formula_sw_c_gee

#> y ~ trt + factor(period) + x1 + x2

5.5 Fit the model-robust standardized estimator

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

5.6 Review the treatment-effect estimates

summary(
  fit_sw_c_gee,
  level       = 0.95,
  show_counts = TRUE,
  ics         = "global",
  digits      = 4
)

#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: gee   Family: gaussian   Scale: RD
#> GEE corstr: exchangeable
#> Clusters: 30   Periods kept: 4 (of 6 total)
#> Kept periods: 2, 3, 4, 5
#> CIs: 95.0% (t, df = 29)
#> Stepped-wedge edge pattern: YES (first period 1 all control; last period 6 all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> --------------------------------------------------
#> # A tibble: 6 × 4
#>   period  minZ  maxZ is_mixed
#>    <int> <int> <int> <lgl>   
#> 1      1     0     0 FALSE   
#> 2      2     0     1 TRUE    
#> 3      3     0     1 TRUE    
#> 4      4     0     1 TRUE    
#> 5      5     0     1 TRUE    
#> 6      6     1     1 FALSE   
#> 
#> Counts used in aggregation (kept periods only)
#> --------------------------------------------------
#> # A tibble: 4 × 4
#>   period n_obs n_clusters n_cp_cells
#>    <int> <int>      <int>      <int>
#> 1      2  1767         30         30
#> 2      3  1555         30         30
#> 3      4  1808         30         30
#> 4      5  1826         30         30
#> 
#> Per-cluster counts (kept periods only)
#> --------------------------------------------------
#> # A tibble: 30 × 3
#>    cluster n_obs n_periods
#>      <int> <int>     <int>
#>  1       1   164         4
#>  2       2   265         4
#>  3       3   120         4
#>  4       4   277         4
#>  5       5   225         4
#>  6       6   223         4
#>  7       7   319         4
#>  8       8   290         4
#>  9       9   250         4
#> 10      10   285         4
#> # ℹ 20 more rows
#> 
#> ICS hypothesis test (linear-contrast F test)
#> --------------------------------------------------
#> Contrasts (rows of C):
#>      h-iATE h-cATE v-iATE v-cATE
#> [1,]      1     -1      0      0
#> [2,]      0      0      1     -1
#> [3,]      1      0     -1      0
#> 
#>   method_type df1 df2       F p_value note
#> 1  unadjusted   3  29 13.1406       0     
#> 2    adjusted   3  29 13.4264       0     
#> 
#> h-iATE
#>            Estimate     SE    LCL     UCL
#> unadjusted   8.2021 1.2483 5.6489 10.7552
#> adjusted     8.2697 1.2488 5.7157 10.8238
#> 
#> h-cATE
#>            Estimate     SE    LCL     UCL
#> unadjusted   7.3619 1.3456 4.6099 10.1138
#> adjusted     7.4190 1.3321 4.6946 10.1434
#> 
#> v-iATE
#>            Estimate     SE    LCL     UCL
#> unadjusted   8.1319 1.2381 5.5996 10.6642
#> adjusted     8.2006 1.2391 5.6662 10.7349
#> 
#> v-cATE
#>            Estimate     SE    LCL    UCL
#> unadjusted   5.4994 1.1269 3.1945 7.8042
#> adjusted     5.5845 1.1083 3.3176 7.8513

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

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

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

6.1 Specify the mixed-effects working model

A random intercept for cluster is included to represent between-cluster outcome heterogeneity.

formula_sw_c_lmer <- y ~ trt + factor(period) + x1 + x2 + (1 | cluster)
formula_sw_c_lmer

#> y ~ trt + factor(period) + x1 + x2 + (1 | cluster)

6.2 Fit the estimator

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

6.3 Review the treatment-effect estimates

summary(
  fit_sw_c_lmer,
  level       = 0.95,
  show_counts = FALSE,
  ics         = "global",
  digits      = 4
)

#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: lmer   Family: gaussian   Scale: RD
#> Clusters: 30   Periods kept: 4 (of 6 total)
#> Kept periods: 2, 3, 4, 5
#> CIs: 95.0% (t, df = 29)
#> Stepped-wedge edge pattern: YES (first period 1 all control; last period 6 all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> --------------------------------------------------
#> # A tibble: 6 × 4
#>   period  minZ  maxZ is_mixed
#>    <int> <int> <int> <lgl>   
#> 1      1     0     0 FALSE   
#> 2      2     0     1 TRUE    
#> 3      3     0     1 TRUE    
#> 4      4     0     1 TRUE    
#> 5      5     0     1 TRUE    
#> 6      6     1     1 FALSE   
#> 
#> ICS hypothesis test (linear-contrast F test)
#> --------------------------------------------------
#> Contrasts (rows of C):
#>      h-iATE h-cATE v-iATE v-cATE
#> [1,]      1     -1      0      0
#> [2,]      0      0      1     -1
#> [3,]      1      0     -1      0
#> 
#>   method_type df1 df2       F p_value note
#> 1  unadjusted   3  29 13.1406       0     
#> 2    adjusted   3  29 13.4264       0     
#> 
#> h-iATE
#>            Estimate     SE    LCL     UCL
#> unadjusted   8.2021 1.2483 5.6489 10.7552
#> adjusted     8.2697 1.2488 5.7157 10.8238
#> 
#> h-cATE
#>            Estimate     SE    LCL     UCL
#> unadjusted   7.3619 1.3456 4.6099 10.1138
#> adjusted     7.4190 1.3321 4.6946 10.1434
#> 
#> v-iATE
#>            Estimate     SE    LCL     UCL
#> unadjusted   8.1319 1.2381 5.5996 10.6642
#> adjusted     8.2006 1.2391 5.6662 10.7349
#> 
#> v-cATE
#>            Estimate     SE    LCL    UCL
#> unadjusted   5.4994 1.1269 3.1945 7.8042
#> adjusted     5.5845 1.1083 3.3176 7.8513

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

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

bind_rows(
  extract_estimates(fit_sw_c_gee,  "Gaussian GEE"),
  extract_estimates(fit_sw_c_lmer, "Linear mixed-effects model")
) %>%
  kable(digits = 4)

working_model	method_type	h-iATE	h-cATE	v-iATE	v-cATE
Gaussian GEE	unadjusted	8.2021	7.3619	8.1319	5.4994
Gaussian GEE	adjusted	8.2697	7.4190	8.2006	5.5845
Linear mixed-effects model	unadjusted	8.2021	7.3619	8.1319	5.4994
Linear mixed-effects model	adjusted	8.2697	7.4190	8.2006	5.5845

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.

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

7.1 Load and prepare the data

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

#> Number of rows: 3386

cat("Number of clusters:", n_distinct(sw_b$cluster), "\n")

#> Number of clusters: 30

cat("Number of periods:", n_distinct(sw_b$period), "\n")

#> Number of periods: 4

head(sw_b) %>%
  kable()

cluster	period	id	x2	y
1	1	1	-0.0825638	1
1	1	2	-0.0001809	1
1	1	3	-0.1694425	0
1	1	4	-0.0588241	0
1	1	5	0.0374889	1
1	1	6	0.0682931	1

7.2 Verify the stepped-wedge treatment structure

check_cluster_period_treatment(sw_b)

#> PASS: treatment assignment is constant within every cluster-period cell.

period_treatment_table(sw_b) %>%
  kable()

period	control_clusters	treated_clusters	total_clusters	informative
1	30	0	30	FALSE
2	20	10	30	TRUE
3	10	20	30	TRUE
4	0	30	30	FALSE

7.3 Inspect cluster-period sizes

cluster_period_size_summary(sw_b) %>%
  kable(digits = 2)

cluster_periods	minimum	first_quartile	median	mean	third_quartile	maximum
120	5	16	28.5	28.22	41	50

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

summarize_outcome_by_treatment(sw_b, outcome = "y") %>%
  kable()

Treatment	N	Mean	SD	Median	Q1	Q3
Control	1671	0.543	0.498	1	0	1
Intervention	1715	0.898	0.303	1	1	1

7.5 Specify the working model

formula_sw_b_gee <- y ~ trt + factor(period) + x1 + x2
formula_sw_b_gee

#> y ~ trt + factor(period) + x1 + x2

7.6 Fit the estimator on the risk-difference scale

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

7.7 Review the treatment-effect estimates

summary(
  fit_sw_b_gee_rd,
  level       = 0.95,
  show_counts = TRUE,
  ics         = "global",
  digits      = 4
)

#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: gee   Family: binomial   Scale: RD
#> GEE corstr: exchangeable
#> Clusters: 30   Periods kept: 2 (of 4 total)
#> Kept periods: 2, 3
#> CIs: 95.0% (t, df = 29)
#> Stepped-wedge edge pattern: YES (first period 1 all control; last period 4 all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> --------------------------------------------------
#> # A tibble: 4 × 4
#>   period  minZ  maxZ is_mixed
#>    <int> <int> <int> <lgl>   
#> 1      1     0     0 FALSE   
#> 2      2     0     1 TRUE    
#> 3      3     0     1 TRUE    
#> 4      4     1     1 FALSE   
#> 
#> Counts used in aggregation (kept periods only)
#> --------------------------------------------------
#> # A tibble: 2 × 4
#>   period n_obs n_clusters n_cp_cells
#>    <int> <int>      <int>      <int>
#> 1      2   894         30         30
#> 2      3   823         30         30
#> 
#> Per-cluster counts (kept periods only)
#> --------------------------------------------------
#> # A tibble: 30 × 3
#>    cluster n_obs n_periods
#>      <int> <int>     <int>
#>  1       1    35         2
#>  2       2    44         2
#>  3       3    85         2
#>  4       4    28         2
#>  5       5    86         2
#>  6       6    68         2
#>  7       7    57         2
#>  8       8    94         2
#>  9       9    62         2
#> 10      10    33         2
#> # ℹ 20 more rows
#> 
#> ICS hypothesis test (linear-contrast F test)
#> --------------------------------------------------
#> Contrasts (rows of C):
#>      h-iATE h-cATE v-iATE v-cATE
#> [1,]      1     -1      0      0
#> [2,]      0      0      1     -1
#> [3,]      1      0     -1      0
#> 
#>   method_type df1 df2      F p_value note
#> 1  unadjusted   3  29 2.6359  0.0685     
#> 2    adjusted   3  29 2.5067  0.0786     
#> 
#> h-iATE
#>            Estimate     SE    LCL    UCL
#> unadjusted   0.3773 0.0378 0.3000 0.4547
#> adjusted     0.3754 0.0361 0.3015 0.4494
#> 
#> h-cATE
#>            Estimate     SE    LCL    UCL
#> unadjusted   0.3496 0.0459 0.2557 0.4436
#> adjusted     0.3489 0.0442 0.2584 0.4393
#> 
#> v-iATE
#>            Estimate     SE    LCL    UCL
#> unadjusted   0.3758 0.0377 0.2987 0.4529
#> adjusted     0.3739 0.0359 0.3004 0.4474
#> 
#> v-cATE
#>            Estimate     SE    LCL    UCL
#> unadjusted   0.3206 0.0468 0.2248 0.4163
#> adjusted     0.3206 0.0451 0.2284 0.4127

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

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

7.9 Optional alternative binary scales

The package also allows alternative binary-outcome scales.

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
)

#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: gee   Family: binomial   Scale: OR
#> GEE corstr: exchangeable
#> Clusters: 30   Periods kept: 2 (of 4 total)
#> Kept periods: 2, 3
#> CIs: 95.0% (t, df = 29)
#> Stepped-wedge edge pattern: YES (first period 1 all control; last period 4 all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> --------------------------------------------------
#> # A tibble: 4 × 4
#>   period  minZ  maxZ is_mixed
#>    <int> <int> <int> <lgl>   
#> 1      1     0     0 FALSE   
#> 2      2     0     1 TRUE    
#> 3      3     0     1 TRUE    
#> 4      4     1     1 FALSE   
#> 
#> Counts used in aggregation (kept periods only)
#> --------------------------------------------------
#> # A tibble: 2 × 4
#>   period n_obs n_clusters n_cp_cells
#>    <int> <int>      <int>      <int>
#> 1      2   894         30         30
#> 2      3   823         30         30
#> 
#> Per-cluster counts (kept periods only)
#> --------------------------------------------------
#> # A tibble: 30 × 3
#>    cluster n_obs n_periods
#>      <int> <int>     <int>
#>  1       1    35         2
#>  2       2    44         2
#>  3       3    85         2
#>  4       4    28         2
#>  5       5    86         2
#>  6       6    68         2
#>  7       7    57         2
#>  8       8    94         2
#>  9       9    62         2
#> 10      10    33         2
#> # ℹ 20 more rows
#> 
#> ICS hypothesis test (linear-contrast F test)
#> --------------------------------------------------
#> Contrasts (rows of C):
#>      h-iATE h-cATE v-iATE v-cATE
#> [1,]      1     -1      0      0
#> [2,]      0      0      1     -1
#> [3,]      1      0     -1      0
#> 
#>   method_type df1 df2      F p_value note
#> 1  unadjusted   3  29 7.4666   8e-04     
#> 2    adjusted   3  29 7.1961   9e-04     
#> 
#> h-iATE
#>            Estimate     SE    LCL    UCL
#> unadjusted    2.296 0.2660 1.7519 2.8401
#> adjusted      2.283 0.2598 1.7517 2.8144
#> 
#> Exponentiated (OR)
#>             Ratio Ratio_LCL Ratio_UCL
#> unadjusted 9.9348    5.7658   17.1180
#> adjusted   9.8064    5.7644   16.6826
#> 
#> h-cATE
#>            Estimate     SE    LCL    UCL
#> unadjusted   2.0032 0.3273 1.3337 2.6727
#> adjusted     1.9970 0.3228 1.3368 2.6572
#> 
#> Exponentiated (OR)
#>             Ratio Ratio_LCL Ratio_UCL
#> unadjusted 7.4127    3.7952   14.4784
#> adjusted   7.3669    3.8067   14.2569
#> 
#> v-iATE
#>            Estimate     SE    LCL    UCL
#> unadjusted   2.2861 0.2648 1.7445 2.8277
#> adjusted     2.2734 0.2585 1.7448 2.8020
#> 
#> Exponentiated (OR)
#>             Ratio Ratio_LCL Ratio_UCL
#> unadjusted 9.8365    5.7232   16.9060
#> adjusted   9.7127    5.7249   16.4782
#> 
#> v-cATE
#>            Estimate     SE    LCL    UCL
#> unadjusted   1.7628 0.2873 1.1751 2.3504
#> adjusted     1.7608 0.2819 1.1843 2.3373
#> 
#> Exponentiated (OR)
#>             Ratio Ratio_LCL Ratio_UCL
#> unadjusted 5.8285    3.2385   10.4897
#> adjusted   5.8172    3.2684   10.3537

The selected scale should correspond to the scientific question, the outcome frequency, and the reporting conventions of the relevant field.

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

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.

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

A Tutorial for Model-Robust Standardization in Stepped-Wedge Cluster-Randomized Trials Using `MRStdLCRT`

June 01, 2026

1 Introduction

1.1 Estimands reported by `MRStdLCRT`

1.2 Estimators and inference

2 Software requirements

3 Data structure required for analysis

4 Helper functions

5 Example 1: Continuous outcome with a GEE working model

5.1 Load and prepare the data

5.2 Verify the stepped-wedge treatment structure

5.3 Summarize the outcome descriptively

5.4 Specify the working model

5.5 Fit the model-robust standardized estimator

5.6 Review the treatment-effect estimates

5.7 Interpretation

6 Example 2: Continuous outcome with a linear mixed-effects working model

6.1 Specify the mixed-effects working model

6.2 Fit the estimator

6.3 Review the treatment-effect estimates

6.4 Compare the two continuous-outcome working models

7 Example 3: Binary outcome with a logistic GEE working model

7.1 Load and prepare the data

7.2 Verify the stepped-wedge treatment structure

7.3 Inspect cluster-period sizes

7.4 Summarize the outcome descriptively

7.5 Specify the working model

7.6 Fit the estimator on the risk-difference scale

7.7 Review the treatment-effect estimates

7.8 Interpretation

7.9 Optional alternative binary scales

7.10 Working-model selection

8 References

A Tutorial for Model-Robust Standardization in Stepped-Wedge Cluster-Randomized Trials Using MRStdLCRT

June 01, 2026

1 Introduction

1.1 Estimands reported by MRStdLCRT

1.2 Estimators and inference

2 Software requirements

3 Data structure required for analysis

4 Helper functions

5 Example 1: Continuous outcome with a GEE working model

5.1 Load and prepare the data

5.2 Verify the stepped-wedge treatment structure

5.3 Summarize the outcome descriptively

5.4 Specify the working model

5.5 Fit the model-robust standardized estimator

5.6 Review the treatment-effect estimates

5.7 Interpretation

6 Example 2: Continuous outcome with a linear mixed-effects working model

6.1 Specify the mixed-effects working model

6.2 Fit the estimator

6.3 Review the treatment-effect estimates

6.4 Compare the two continuous-outcome working models

7 Example 3: Binary outcome with a logistic GEE working model

7.1 Load and prepare the data

7.2 Verify the stepped-wedge treatment structure

7.3 Inspect cluster-period sizes

7.4 Summarize the outcome descriptively

7.5 Specify the working model

7.6 Fit the estimator on the risk-difference scale

7.7 Review the treatment-effect estimates

7.8 Interpretation

7.9 Optional alternative binary scales

7.10 Working-model selection

8 References

A Tutorial for Model-Robust Standardization in Stepped-Wedge Cluster-Randomized Trials Using `MRStdLCRT`

1.1 Estimands reported by `MRStdLCRT`