1 Overview

MRStdLCRT implements model-robust standardization for longitudinal cluster randomized trials (LCRTs), with emphasis on crossover cluster-randomized trials (CRTs) in this tutorial. The package computes four estimands using either unadjusted or augmented (model-robust standardized) estimators.

1.1 Estimands

h-iATE: horizontal individual average treatment effect
h-cATE: horizontal cluster average treatment effect
v-iATE: vertical individual average treatment effect
v-cATE: vertical cluster average treatment effect

In crossover CRTs, we will emphasize h-iATE and h-cATE only, but the API returns all four.

1.2 Estimators and variance

Unadjusted: non-parametric estimators
Adjusted: augmented estimators based on model-robust standardization with covariate adjustment. In this pacakge, the working model can be:
- Linear mixed effects model (lmer) for continuous outcomes
- Generalized linear mixed effects model (glmer) for binary outcomes
- Generalized estimating equations (gee) for continuous or binary outcomes
Variance: delete-1 cluster jackknife SEs and covariance

1.3 Type of the design

This pacakge can identify the design type based on the treatment assignment pattern. In this tutorial, we focus on crossover CRTs (each cluster receives both treatment and control at least once). The package also supports stepped-wedge CRTs and parallel CRTs (no crossover).

The working model is based on all observations in the dataset, but the aggregation of cluster-period means differs by different designs.

A period is kept if it has a mixture of treatment assignments across clusters (both 0 and 1 occur among cluster-period assignments in that period).
Periods with all clusters treated or all clusters control are excluded from aggregation, but the working model is still fit on the full data.

2 Installation

Install the development version directly from GitHub:

# install.packages("remotes")
remotes::install_github("fancy575/MRStdLCRT")

Load the package:

suppressPackageStartupMessages({
  library(MRStdLCRT)
  library(dplyr)
  library(ggplot2)
})

3 Example data with crossover CRT

Two example datasets are included for this tutorial:

xo_c: continuous outcome dataset (outcome column: y_cont)
xo_b: binary outcome dataset (outcome column: y_bin)

Both contain:

cluster ID: h
period: p
treatment indicator: trt (0/1)
covariates: x_c01, x_c02, x_b01, x_cat1_2, x_cat1_3

data("xo_c", package = "MRStdLCRT")
data("xo_b", package = "MRStdLCRT")

3.1 Design overview:

The heatmap below shows the treatment assignment for each cluster-period cell.
Cells in blue indicate treatment (trt=1); cells in gray indicate control (trt=0).

3.2 Outcome overview

3.2.1 Continuous outcome (`xo_c`)

3.2.2 Binary outcome (`xo_b`)

4 API at a glance

Fit
mrstdlcrt_fit(data, formula, cluster_id, period, trt, method, family, corstr, scale)
Summarize
summary(object, level = 0.95, estimand = c("h-iATE","h-cATE"))
Plot
plot(object, estimand = c("h-iATE","h-cATE"))

4.1 Notes on formulas

For GEE, it is recommended to omit random effects terms like (1|h); the package strips them internally if present.
The package supports complex fixed effects, including interactions such as factor(p) * trt, etc.

5 Continuous outcome (crossover) — four working models

We keep the same covariates and vary only the working model specification.

fml_c_lmer_full <- y_cont ~ 0 + factor(p) + trt +
  x_c01 + x_c02 + x_b01 + x_cat1_2 + x_cat1_3 +
  (1 | h) + (1 | h:p)

fml_c_lmer_h <- y_cont ~ 0 + factor(p) + trt +
  x_c01 + x_c02 + x_b01 + x_cat1_2 + x_cat1_3 +
  (1 | h)

fml_c_gee <- y_cont ~ 0 + factor(p) + trt +
  x_c01 + x_c02 + x_b01 + x_cat1_2 + x_cat1_3

5.1 lmer: random intercepts (1|h) + (1|h:p)

fit_c_lmer_full <- mrstdlcrt_fit(
  data       = xo_c,
  formula    = fml_c_lmer_full,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "lmer",
  family     = "gaussian"
)

summary(fit_c_lmer_full, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: lmer   Family: gaussian   Scale: RD
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159120 0.132915 -1.484352 -0.833889
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179611 0.130738 -1.499516 -0.859707
plot(fit_c_lmer_full, estimand = c("h-iATE","h-cATE"))

5.2 lmer: random intercept (1|h)

fit_c_lmer_h <- mrstdlcrt_fit(
  data       = xo_c,
  formula    = fml_c_lmer_h,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "lmer",
  family     = "gaussian"
)

summary(fit_c_lmer_h, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: lmer   Family: gaussian   Scale: RD
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159147 0.132820 -1.484145 -0.834150
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179642 0.130622 -1.499264 -0.860021

5.3 GEE: independence

fit_c_gee_ind <- mrstdlcrt_fit(
  data       = xo_c,
  formula    = fml_c_gee,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "gaussian",
  corstr     = "independence"
)

summary(fit_c_gee_ind, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: gee   Family: gaussian   Scale: RD
#> GEE corstr: independence
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159131 0.132818 -1.484124 -0.834138
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179568 0.130649 -1.499254 -0.859883

5.4 GEE: exchangeable

fit_c_gee_exch <- mrstdlcrt_fit(
  data       = xo_c,
  formula    = fml_c_gee,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "gaussian",
  corstr     = "exchangeable"
)

summary(fit_c_gee_exch, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: gee   Family: gaussian   Scale: RD
#> GEE corstr: exchangeable
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159147 0.132820 -1.484145 -0.834150
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179642 0.130623 -1.499264 -0.860021

5.5 Interaction example (continuous outcome)

This example illustrates fixed-effect interactions that frequently arise in crossover analyses.

fml_c_lmer_int <- y_cont ~ 0 + factor(p) * trt +
  factor(p) * x_c01 + x_c02 + x_b01 + x_cat1_2 + x_cat1_3 +
  (1 | h) + (1 | h:p)

fit_c_lmer_int <- mrstdlcrt_fit(
  data       = xo_c,
  formula    = fml_c_lmer_int,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "lmer",
  family     = "gaussian"
)

summary(fit_c_lmer_int, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: lmer   Family: gaussian   Scale: RD
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159155 0.132941 -1.484450 -0.833860
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179632 0.130749 -1.499563 -0.859702

6 Binary outcome (crossover) — four working models and reporting scale

For binary outcomes, you can specify the reporting scale via scale:

scale = "RD": risk difference
scale = "RR": log risk ratio (summary also prints exponentiated RR)
scale = "OR": log odds ratio (summary also prints exponentiated OR)

fml_b_glmer_full <- y_bin ~ 0 + factor(p) + trt +
  x_c01 + x_c02 + x_b01 + x_cat1_2 + x_cat1_3 +
  (1 | h) + (1 | h:p)

fml_b_glmer_h <- y_bin ~ 0 + factor(p) + trt +
  x_c01 + x_c02 + x_b01 + x_cat1_2 + x_cat1_3 +
  (1 | h)

fml_b_gee <- y_bin ~ 0 + factor(p) + trt +
  x_c01 + x_c02 + x_b01 + x_cat1_2 + x_cat1_3

6.1 glmer: (1|h) + (1|h:p), scale = OR (log OR)

fit_b_glmer_full_or <- mrstdlcrt_fit(
  data       = xo_b,
  formula    = fml_b_glmer_full,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "glmer",
  family     = "binomial",
  scale      = "OR"
)

summary(fit_b_glmer_full_or, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: glmer   Family: binomial   Scale: OR
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.769585 0.132807 -1.094553 -0.444618
#> adjusted   -0.818665 0.099803 -1.062873 -0.574457
#> 
#> Exponentiated (OR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.463205  0.334689  0.641069
#> adjusted   0.441020  0.345462  0.563011
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.780413 0.131000 -1.100960 -0.459867
#> adjusted   -0.840393 0.099997 -1.085076 -0.595709
#> 
#> Exponentiated (OR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.458217  0.332552  0.631367
#> adjusted   0.431541  0.337876  0.551171
plot(fit_b_glmer_full_or, estimand = c("h-iATE","h-cATE"))

6.2 glmer: (1|h), scale = OR (log OR)

fit_b_glmer_h_or <- mrstdlcrt_fit(
  data       = xo_b,
  formula    = fml_b_glmer_h,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "glmer",
  family     = "binomial",
  scale      = "OR"
)

summary(fit_b_glmer_h_or, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: glmer   Family: binomial   Scale: OR
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.769585 0.132807 -1.094553 -0.444618
#> adjusted   -0.818440 0.099916 -1.062927 -0.573954
#> 
#> Exponentiated (OR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.463205  0.334689  0.641069
#> adjusted   0.441119  0.345443  0.563294
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.780413 0.131000 -1.100960 -0.459867
#> adjusted   -0.840183 0.100066 -1.085036 -0.595330
#> 
#> Exponentiated (OR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.458217  0.332552  0.631367
#> adjusted   0.431631  0.337890  0.551380

6.3 GEE: independence, scale = OR (log OR)

fit_b_gee_ind_or <- mrstdlcrt_fit(
  data       = xo_b,
  formula    = fml_b_gee,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "binomial",
  corstr     = "independence",
  scale      = "OR"
)

summary(fit_b_gee_ind_or, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: gee   Family: binomial   Scale: OR
#> GEE corstr: independence
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.769585 0.132807 -1.094553 -0.444618
#> adjusted   -0.817738 0.100249 -1.063038 -0.572438
#> 
#> Exponentiated (OR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.463205  0.334689  0.641069
#> adjusted   0.441429  0.345405  0.564148
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.780413 0.131000 -1.100960 -0.459867
#> adjusted   -0.839421 0.100221 -1.084653 -0.594189
#> 
#> Exponentiated (OR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.458217  0.332552  0.631367
#> adjusted   0.431961  0.338019  0.552010

6.4 GEE: exchangeable, scale = OR (log OR)

fit_b_gee_exch_or <- mrstdlcrt_fit(
  data       = xo_b,
  formula    = fml_b_gee,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "binomial",
  corstr     = "exchangeable",
  scale      = "OR"
)

summary(fit_b_gee_exch_or, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: gee   Family: binomial   Scale: OR
#> GEE corstr: exchangeable
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.769585 0.132807 -1.094553 -0.444618
#> adjusted   -0.818096 0.100348 -1.063639 -0.572552
#> 
#> Exponentiated (OR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.463205  0.334689  0.641069
#> adjusted   0.441271  0.345197  0.564084
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.780413 0.131000 -1.100960 -0.459867
#> adjusted   -0.839700 0.100487 -1.085583 -0.593816
#> 
#> Exponentiated (OR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.458217  0.332552  0.631367
#> adjusted   0.431840  0.337705  0.552216

6.5 Alternative scales (RR and RD)

fit_b_gee_ind_rr <- mrstdlcrt_fit(
  data       = xo_b,
  formula    = fml_b_gee,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "binomial",
  corstr     = "independence",
  scale      = "RR"
)

fit_b_gee_ind_rd <- mrstdlcrt_fit(
  data       = xo_b,
  formula    = fml_b_gee,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "binomial",
  corstr     = "independence",
  scale      = "RD"
)

summary(fit_b_gee_ind_rr, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: gee   Family: binomial   Scale: RR
#> GEE corstr: independence
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.370624 0.051181 -0.495859 -0.245390
#> adjusted   -0.394445 0.037920 -0.487231 -0.301659
#> 
#> Exponentiated (RR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.690303  0.609047   0.78240
#> adjusted   0.674054  0.614325   0.73959
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.379936 0.053852 -0.511707 -0.248164
#> adjusted   -0.408698 0.042047 -0.511582 -0.305813
#> 
#> Exponentiated (RR)
#>               Ratio Ratio_LCL Ratio_UCL
#> unadjusted 0.683906  0.599471  0.780232
#> adjusted   0.664515  0.599546  0.736524
summary(fit_b_gee_ind_rd, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-LCRT Summary
#> ========================================================================
#> Method: gee   Family: binomial   Scale: RD
#> GEE corstr: independence
#> Clusters: 7   Periods kept: 2 (of 2 total)
#> Kept periods: 1, 2
#> CIs: 95.0% (t, df = 6)
#> Stepped-wedge edge pattern: NO (first period not all control and/or last period not all treatment)
#> 
#> Period mixture table (min/max of Z_ij by period)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p  minZ  maxZ is_mixed
#>   <int> <int> <int> <lgl>   
#> 1     1     0     1 TRUE    
#> 2     2     0     1 TRUE    
#> 
#> Counts used in aggregation (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 2 × 4
#>       p n_obs n_clusters n_cp_cells
#>   <int> <int>      <int>      <int>
#> 1     1  4455          7          7
#> 2     2  4535          7          7
#> 
#> Per-cluster counts (kept periods only)
#> ------------------------------------------------------------------------
#> # A tibble: 7 × 3
#>       h n_obs n_periods
#>   <int> <int>     <int>
#> 1     1   901         2
#> 2     2  1329         2
#> 3     3  1028         2
#> 4     4  1632         2
#> 5     5  1782         2
#> 6     6  1690         2
#> 7     7   628         2
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.189803 0.031290 -0.266367 -0.113238
#> adjusted   -0.201386 0.023287 -0.258368 -0.144404
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.192533 0.031128 -0.268701 -0.116365
#> adjusted   -0.206688 0.023582 -0.264391 -0.148985

7 Interpreting the output

Each summary() call prints:

Method / family / scale
Which periods were kept by the mixture rule
A period mixture table (min/max treatment by period)
Whether the dataset matches the stepped-wedge CRT pattern

For each requested estimand (e.g., h-iATE, h-cATE), it reports:

Unadjusted and Adjusted estimates
Jackknife SE
t-based CI with df = clusters − 1
If scale is "RR" or "OR", an exponentiated ratio table is also printed

8 Session info

sessionInfo()
#> R version 4.5.2 (2025-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.1   dplyr_1.1.4     MRStdLCRT_0.1.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] Matrix_1.7-4       gtable_0.3.6       jsonlite_2.0.0     compiler_4.5.2    
#>  [5] Rcpp_1.1.0         tidyselect_1.2.1   gee_4.13-29        tidyr_1.3.2       
#>  [9] jquerylib_0.1.4    splines_4.5.2      scales_1.4.0       boot_1.3-32       
#> [13] yaml_2.3.12        fastmap_1.2.0      lattice_0.22-7     R6_2.6.1          
#> [17] labeling_0.4.3     generics_0.1.4     knitr_1.50         rbibutils_2.4     
#> [21] MASS_7.3-65        tibble_3.3.0       nloptr_2.2.1       minqa_1.2.8       
#> [25] bslib_0.9.0        pillar_1.11.1      RColorBrewer_1.1-3 rlang_1.1.6       
#> [29] cachem_1.1.0       xfun_0.55          sass_0.4.10        S7_0.2.1          
#> [33] cli_3.6.5          withr_3.0.2        magrittr_2.0.4     Rdpack_2.6.4      
#> [37] digest_0.6.39      grid_4.5.2         rstudioapi_0.17.1  lme4_1.1-38       
#> [41] nlme_3.1-168       lifecycle_1.0.4    reformulas_0.4.3   vctrs_0.6.5       
#> [45] evaluate_1.0.5     glue_1.8.0         farver_2.1.2       purrr_1.2.0       
#> [49] rmarkdown_2.30     tools_4.5.2        pkgconfig_2.0.3    htmltools_0.5.9

Model-Robust Standardization for Crossover Trials (MRStdLCRT) — Tutorial

Xi Fang

1 Overview

1.1 Estimands

1.2 Estimators and variance

1.3 Type of the design

2 Installation

3 Example data with crossover CRT

3.1 Design overview:

3.2 Outcome overview

3.2.1 Continuous outcome (`xo_c`)

3.2.2 Binary outcome (`xo_b`)

4 API at a glance

4.1 Notes on formulas

5 Continuous outcome (crossover) — four working models

5.1 lmer: random intercepts (1|h) + (1|h:p)

5.2 lmer: random intercept (1|h)

5.3 GEE: independence

5.4 GEE: exchangeable

5.5 Interaction example (continuous outcome)

6 Binary outcome (crossover) — four working models and reporting scale

6.1 glmer: (1|h) + (1|h:p), scale = OR (log OR)

6.2 glmer: (1|h), scale = OR (log OR)

6.3 GEE: independence, scale = OR (log OR)

6.4 GEE: exchangeable, scale = OR (log OR)

6.5 Alternative scales (RR and RD)

7 Interpreting the output

8 Session info

Model-Robust Standardization for Crossover Trials (MRStdLCRT) — Tutorial

Xi Fang

1 Overview

1.1 Estimands

1.2 Estimators and variance

1.3 Type of the design

2 Installation

3 Example data with crossover CRT

3.1 Design overview:

3.2 Outcome overview

3.2.1 Continuous outcome (xo_c)

3.2.2 Binary outcome (xo_b)

4 API at a glance

4.1 Notes on formulas

5 Continuous outcome (crossover) — four working models

5.1 lmer: random intercepts (1|h) + (1|h:p)

5.2 lmer: random intercept (1|h)

5.3 GEE: independence

5.4 GEE: exchangeable

5.5 Interaction example (continuous outcome)

6 Binary outcome (crossover) — four working models and reporting scale

6.1 glmer: (1|h) + (1|h:p), scale = OR (log OR)

6.2 glmer: (1|h), scale = OR (log OR)

6.3 GEE: independence, scale = OR (log OR)

6.4 GEE: exchangeable, scale = OR (log OR)

6.5 Alternative scales (RR and RD)

7 Interpreting the output

8 Session info

3.2.1 Continuous outcome (`xo_c`)

3.2.2 Binary outcome (`xo_b`)