MRStdLCRT implements model‑robust standardization for longitudinal cluster randomized trials (CRTs), including 2×2 crossover designs. The package computes:
lmer
, glmer
, or
gee
)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.
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.
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)
})
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")
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")
)
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)")
A quick sense of marginal distributions by treatment
dat_c
)dat_b
)mrstdlcrt_fit(data, formula, cluster_id, period, trt, method, family, corstr = "independence")
summary(object, level = 0.95, estimand = c("h-iATE","h-cATE"))
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.
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
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"))
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
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
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
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
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
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"))
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
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
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
Each summary()
prints two blocks (because we asked for
estimand = c("h-iATE","h-cATE")
):
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.
(1|...)
. The tutorial examples do this explicitly.summary(fit, level = 0.90)
or 0.99
.summary(fit)
(no estimand=...
) to print h‑iATE,
h‑cATE, v‑iATE,
v‑cATE.plot(fit, estimand = c("h-iATE","h-cATE"))
draws a simple
comparison for the two estimands.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