1 Overview

MRStdLCRT implements model‑robust standardization for longitudinal cluster randomized trials (CRTs), including 2×2 crossover designs. The package computes:

  • Estimands supported by the package
    • h‑iATE: hospital‑weighted individual average treatment effect (horizontal aggregation across periods)
    • h‑cATE: hospital‑weighted cluster average treatment effect (horizontal)
    • v‑iATE and v‑cATE: vertical estimands (not emphasized here, but supported by the API)
  • Estimators
    • Unadjusted (purely design/weight-based)
    • Adjusted (model‑robust standardization using a working outcome model: lmer, glmer, or gee)
  • Uncertainty
    • Delete‑1 cluster jackknife standard errors (SEs) for each estimand and their difference (i‑ATE − c‑ATE)

Focus of this tutorial: we keep your preferred style — four fits for each outcome (continuous & binary) — and emphasize h‑iATE and h‑cATE with SW = FALSE (i.e., using all periods). The package itself is designed for longitudinal CRT analysis and provides all four estimands.

1.1 What does “model‑robust standardization” mean?

We fit a working model for the outcome (GLMM via lmer/glmer or GEE via gee). We then standardize predicted potential outcomes within CPs and aggregate with design‑driven weights (i‑ATE vs c‑ATE). Because the standardization step is design‑based, the effect estimands remain well‑defined even if the working model is misspecified.

2 Installation

The package can be downloaded from https://github.com/fancy575/MRStd_XO, and then install from your local source and load:

install.packages("~/path/to/MRStdLCRT_0.1.0.tar.gz", repos = NULL, type = "source")
suppressPackageStartupMessages({
  library(MRStdLCRT)

})

3 Example data

Two simulated datasets are included:

  • dat_c: continuous outcome (y_cont)
  • dat_b: binary outcome (y_bin)

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

data("dat_c", package = "MRStdLCRT")
data("dat_b", package = "MRStdLCRT")

3.1 Peek at the data

datatable(
  dat_c,
  rownames = FALSE,
  options = list(pageLength = 10, scrollX = TRUE),
  caption = htmltools::tags$caption("dat_c: continuous outcome")
)
datatable(
  dat_b,
  rownames = FALSE,
  options = list(pageLength = 10, scrollX = TRUE),
  caption = htmltools::tags$caption("dat_b: binary outcome")
)

3.2 Cluster‑period sizes (Nij)

cp_sizes_c <- count(dat_c, h, p, name = "Nij")
cp_sizes_b <- count(dat_b, h, p, name = "Nij")

datatable(cp_sizes_c, rownames = FALSE, options = list(dom = 'tp', pageLength = 7),
          caption = "dat_c: cluster‑period sizes (Nij)")
datatable(cp_sizes_b, rownames = FALSE, options = list(dom = 'tp', pageLength = 7),
          caption = "dat_b: cluster‑period sizes (Nij)")

3.3 Data overview

A quick sense of marginal distributions by treatment

3.4 Continuous outcome (dat_c)

3.5 Binary outcome (dat_b)

4 API at a glance

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

Notes

  • For GEE, omit random effects from the formula (no (1|...)). The package internally strips them if present, but it’s clearer to leave them out in your code.
  • We always show h‑iATE and h‑cATE. The package can also compute v‑iATE and v‑cATE; set estimand = NULL to print all four.

5 Continuous outcome — four models

We keep the same covariates and change only the working model.

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: (1|h) + (1|h:p)

fit_c_lmer_full <- mrstdlcrt_fit(
  data       = dat_c,
  formula    = fml_c_lmer_full,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "lmer",
  family     = "gaussian"
)
summary(fit_c_lmer_full, level = 0.95, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-XO Summary (all periods)
#> ========================================================================
#> Method: lmer   (family: gaussian)
#> Clusters: 7   Periods: 2   N: 8990
#> Cluster id: h   Period: p   Trt: trt
#> Formula used for fitting: y_cont ~ 0 + factor(p) + trt + x_c01 + x_c02 + x_b01 + x_cat1_2 + 
#>     x_cat1_3 + (1 | h) + (1 | h:p)
#> CIs: 95.0% (t, df = 6)
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159120 0.132915 -1.484352 -0.833889
#> 
#> Coef (model):
#>        TE       SE     Var
#> 1 -1.1622 0.132777 0.01763
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179611 0.130738 -1.499516 -0.859707
#> 
#> Coef (model):
#>        TE       SE     Var
#> 1 -1.1622 0.132777 0.01763
# Optional visual summary for these two estimands
plot(fit_c_lmer_full, estimand = c("h-iATE","h-cATE"))

5.2 lmer: (1|h)

fit_c_lmer_h <- mrstdlcrt_fit(
  data       = dat_c,
  formula    = fml_c_lmer_h,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "lmer",
  family     = "gaussian"
)
summary(fit_c_lmer_h, level = 0.95, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-XO Summary (all periods)
#> ========================================================================
#> Method: lmer   (family: gaussian)
#> Clusters: 7   Periods: 2   N: 8990
#> Cluster id: h   Period: p   Trt: trt
#> Formula used for fitting: y_cont ~ 0 + factor(p) + trt + x_c01 + x_c02 + x_b01 + x_cat1_2 + 
#>     x_cat1_3 + (1 | h)
#> CIs: 95.0% (t, df = 6)
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159147 0.132820 -1.484145 -0.834150
#> 
#> Coef (model):
#>          TE       SE      Var
#> 1 -1.151398 0.027046 0.000731
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179642 0.130622 -1.499264 -0.860021
#> 
#> Coef (model):
#>          TE       SE      Var
#> 1 -1.151398 0.027046 0.000731

5.3 GEE: independence

fit_c_gee_ind <- mrstdlcrt_fit(
  data       = dat_c,
  formula    = fml_c_gee,      # no (1|...) for GEE
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "gaussian",
  corstr     = "independence"
)
summary(fit_c_gee_ind, level = 0.95, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-XO Summary (all periods)
#> ========================================================================
#> Method: gee   (family: gaussian)
#> GEE corstr: independence
#> Clusters: 7   Periods: 2   N: 8990
#> Cluster id: h   Period: p   Trt: trt
#> Formula used for fitting: y_cont ~ 0 + factor(p) + trt + x_c01 + x_c02 + x_b01 + x_cat1_2 + 
#>     x_cat1_3
#> CIs: 95.0% (t, df = 6)
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159131 0.132818 -1.484124 -0.834138
#> 
#> Coef (model):
#>          TE      SE      Var
#> 1 -1.167619 0.10165 0.010333
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179568 0.130649 -1.499254 -0.859883
#> 
#> Coef (model):
#>          TE      SE      Var
#> 1 -1.167619 0.10165 0.010333

5.4 GEE: exchangeable

fit_c_gee_exch <- mrstdlcrt_fit(
  data       = dat_c,
  formula    = fml_c_gee,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "gaussian",
  corstr     = "exchangeable"
)
summary(fit_c_gee_exch, level = 0.95, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-XO Summary (all periods)
#> ========================================================================
#> Method: gee   (family: gaussian)
#> GEE corstr: exchangeable
#> Clusters: 7   Periods: 2   N: 8990
#> Cluster id: h   Period: p   Trt: trt
#> Formula used for fitting: y_cont ~ 0 + factor(p) + trt + x_c01 + x_c02 + x_b01 + x_cat1_2 + 
#>     x_cat1_3
#> CIs: 95.0% (t, df = 6)
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159147 0.132820 -1.484145 -0.834150
#> 
#> Coef (model):
#>          TE       SE      Var
#> 1 -1.151379 0.076398 0.005837
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179642 0.130623 -1.499264 -0.860021
#> 
#> Coef (model):
#>          TE       SE      Var
#> 1 -1.151379 0.076398 0.005837

5.5 Interaction example for continous outcome

Below we add period×treatment interaction using formula
0 + factor(p) * trt. We use the lmer as example:

fml_c_lmer_full_int <- y_cont ~ 0 + factor(p) * trt +
  x_c01 + x_c02 + x_b01 + x_cat1_2 + x_cat1_3 +
  (1 | h) + (1 | h:p)
fit_c_lmer_full_int <- mrstdlcrt_fit(
  data       = dat_c,
  formula    = fml_c_lmer_full_int,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "lmer",
  family     = "gaussian"
)
summary(fit_c_lmer_full_int, level = 0.95, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-XO Summary (all periods)
#> ========================================================================
#> Method: lmer   (family: gaussian)
#> Clusters: 7   Periods: 2   N: 8990
#> Cluster id: h   Period: p   Trt: trt
#> Formula used for fitting: y_cont ~ 0 + factor(p) * trt + x_c01 + x_c02 + x_b01 + x_cat1_2 + 
#>     x_cat1_3 + (1 | h) + (1 | h:p)
#> CIs: 95.0% (t, df = 6)
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.093397 0.163533 -1.493549 -0.693245
#> adjusted   -1.159120 0.132915 -1.484351 -0.833889
#> 
#> Coef (model):
#>                 TE       SE      Var
#> estimate -1.160227 0.132736 0.017619
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -1.104814 0.157185 -1.489433 -0.720196
#> adjusted   -1.179611 0.130738 -1.499514 -0.859708
#> 
#> Coef (model):
#>                 TE       SE      Var
#> estimate -1.163136 0.132713 0.017613

6 Binary outcome — four models

We adopt the same structure with binary outcomes. Results are on the log‑odds‑ratio scale.

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)

fit_b_glmer_full <- mrstdlcrt_fit(
  data       = dat_b,
  formula    = fml_b_glmer_full,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "glmer",
  family     = "binomial"
)
summary(fit_b_glmer_full, level = 0.95, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-XO Summary (all periods)
#> ========================================================================
#> Method: glmer   (family: binomial)
#> Clusters: 7   Periods: 2   N: 8990
#> Cluster id: h   Period: p   Trt: trt
#> Formula used for fitting: y_bin ~ 0 + factor(p) + trt + x_c01 + x_c02 + x_b01 + x_cat1_2 + 
#>     x_cat1_3 + (1 | h) + (1 | h:p)
#> CIs: 95.0% (t, df = 6)
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.769585 0.132807 -1.094553 -0.444618
#> adjusted   -0.818665 0.099803 -1.062874 -0.574457
#> 
#> Coef (model):
#>          TE       SE      Var
#> 1 -1.182661 0.128071 0.016402
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.780413 0.131000 -1.100960 -0.459867
#> adjusted   -0.840393 0.099997 -1.085077 -0.595710
#> 
#> Coef (model):
#>          TE       SE      Var
#> 1 -1.182661 0.128071 0.016402
plot(fit_b_glmer_full, estimand = c("h-iATE","h-cATE"))

6.2 glmer: (1|h)

fit_b_glmer_h <- mrstdlcrt_fit(
  data       = dat_b,
  formula    = fml_b_glmer_h,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "glmer",
  family     = "binomial"
)
summary(fit_b_glmer_h, level = 0.95, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-XO Summary (all periods)
#> ========================================================================
#> Method: glmer   (family: binomial)
#> Clusters: 7   Periods: 2   N: 8990
#> Cluster id: h   Period: p   Trt: trt
#> Formula used for fitting: y_bin ~ 0 + factor(p) + trt + x_c01 + x_c02 + x_b01 + x_cat1_2 + 
#>     x_cat1_3 + (1 | h)
#> CIs: 95.0% (t, df = 6)
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.769585 0.132807 -1.094553 -0.444618
#> adjusted   -0.818440 0.099916 -1.062927 -0.573953
#> 
#> Coef (model):
#>          TE      SE      Var
#> 1 -1.167529 0.06887 0.004743
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.780413 0.131000 -1.100960 -0.459867
#> adjusted   -0.840183 0.100066 -1.085037 -0.595330
#> 
#> Coef (model):
#>          TE      SE      Var
#> 1 -1.167529 0.06887 0.004743

6.3 GEE: independence

fit_b_gee_ind <- mrstdlcrt_fit(
  data       = dat_b,
  formula    = fml_b_gee,    # no (1|...) for GEE
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "binomial",
  corstr     = "independence"
)
summary(fit_b_gee_ind, level = 0.95, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-XO Summary (all periods)
#> ========================================================================
#> Method: gee   (family: binomial)
#> GEE corstr: independence
#> Clusters: 7   Periods: 2   N: 8990
#> Cluster id: h   Period: p   Trt: trt
#> Formula used for fitting: y_bin ~ 0 + factor(p) + trt + x_c01 + x_c02 + x_b01 + x_cat1_2 + 
#>     x_cat1_3
#> CIs: 95.0% (t, df = 6)
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.769585 0.132807 -1.094553 -0.444618
#> adjusted   -0.817738 0.100249 -1.063038 -0.572438
#> 
#> Coef (model):
#>          TE       SE      Var
#> 1 -1.142398 0.105884 0.011211
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.780413 0.131000 -1.100960 -0.459867
#> adjusted   -0.839421 0.100221 -1.084653 -0.594189
#> 
#> Coef (model):
#>          TE       SE      Var
#> 1 -1.142398 0.105884 0.011211

6.4 GEE: exchangeable

fit_b_gee_exch <- mrstdlcrt_fit(
  data       = dat_b,
  formula    = fml_b_gee,
  cluster_id = "h",
  period     = "p",
  trt        = "trt",
  method     = "gee",
  family     = "binomial",
  corstr     = "exchangeable"
)
summary(fit_b_gee_exch, level = 0.95, estimand = c("h-iATE","h-cATE"))
#> 
#> MRS-XO Summary (all periods)
#> ========================================================================
#> Method: gee   (family: binomial)
#> GEE corstr: exchangeable
#> Clusters: 7   Periods: 2   N: 8990
#> Cluster id: h   Period: p   Trt: trt
#> Formula used for fitting: y_bin ~ 0 + factor(p) + trt + x_c01 + x_c02 + x_b01 + x_cat1_2 + 
#>     x_cat1_3
#> CIs: 95.0% (t, df = 6)
#> 
#> h-iATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.769585 0.132807 -1.094553 -0.444618
#> adjusted   -0.818096 0.100348 -1.063639 -0.572552
#> 
#> Coef (model):
#>         TE       SE      Var
#> 1 -1.12929 0.079166 0.006267
#> 
#> h-cATE
#>             Estimate       SE       LCL       UCL
#> unadjusted -0.780413 0.131000 -1.100960 -0.459867
#> adjusted   -0.839700 0.100487 -1.085583 -0.593816
#> 
#> Coef (model):
#>         TE       SE      Var
#> 1 -1.12929 0.079166 0.006267

7 Interpreting the output

Each summary() prints two blocks (because we asked for estimand = c("h-iATE","h-cATE")):

  • Unadjusted: estimates without covariate adjustment
  • Adjusted: model‑robust standardized estimates

For each block we report Estimate, SE, and a t‑based CI using df = (#clusters − 1) (jackknife convention). The difference (i‑ATE − c‑ATE) is also available via object$estimates$diff and object$jk_se$'SE(diff)' if you print all estimands.

8 Tips, pitfalls, and reproducibility

  • GEE formulas: Do not include (1|...). The tutorial examples do this explicitly.
  • Confidence levels: change via summary(fit, level = 0.90) or 0.99.
  • All four estimands: call summary(fit) (no estimand=...) to print h‑iATE, h‑cATE, v‑iATE, v‑cATE.
  • Plotting: plot(fit, estimand = c("h-iATE","h-cATE")) draws a simple comparison for the two estimands.

9 Session info

sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: AlmaLinux 8.10 (Cerulean Leopard)
#> 
#> Matrix products: default
#> BLAS/LAPACK: FlexiBLAS OPENBLAS;  LAPACK version 3.10.1
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: NA
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] DT_0.33         ggplot2_3.5.1   dplyr_1.1.4     MRStdLCRT_0.1.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.9        utf8_1.2.4        generics_0.1.3    tidyr_1.3.1      
#>  [5] lattice_0.22-6    lme4_1.1-35.5     digest_0.6.36     magrittr_2.0.3   
#>  [9] evaluate_0.24.0   grid_4.4.1        fastmap_1.2.0     jsonlite_1.8.8   
#> [13] Matrix_1.7-0      purrr_1.0.2       fansi_1.0.6       crosstalk_1.2.1  
#> [17] scales_1.3.0      jquerylib_0.1.4   cli_3.6.3         rlang_1.1.4      
#> [21] munsell_0.5.1     splines_4.4.1     withr_3.0.0       cachem_1.1.0     
#> [25] yaml_2.3.8        tools_4.4.1       nloptr_2.1.1      minqa_1.2.8      
#> [29] colorspace_2.1-1  boot_1.3-31       vctrs_0.6.5       R6_2.5.1         
#> [33] lifecycle_1.0.4   htmlwidgets_1.6.4 MASS_7.3-61       pkgconfig_2.0.3  
#> [37] pillar_1.9.0      bslib_0.7.0       gee_4.13-27       gtable_0.3.5     
#> [41] glue_1.7.0        Rcpp_1.0.13-1     highr_0.11        xfun_0.45        
#> [45] tibble_3.2.1      tidyselect_1.2.1  rstudioapi_0.16.0 knitr_1.47       
#> [49] farver_2.1.2      htmltools_0.5.8.1 nlme_3.1-166      labeling_0.4.3   
#> [53] rmarkdown_2.27    compiler_4.4.1