Package {MRStdCRT}

Model-robust Standardization Estimators for the Cluster Randomized Trials

Description

This function performs cluster randomized trials (CRT) analysis using model-robust standardization estimators to estimate the cluster-average and individual-average treatment effect. It handles different outcome mean models (GLM, LMM, GEE, GLMM) and supports continuous, binary, and count outcomes with options for different correlation structures and scales (risk difference, risk ratio and odds ratio).

Usage

MRStdCRT_fit(
  formula,
  data,
  cluster,
  trt,
  trtprob = rep(0.5, nrow(data)),
  method,
  family = gaussian(link = "identity"),
  corstr,
  scale,
  alpha = 0.05,
  model_weights = c("none", "inverse_cluster_size", "estimand_specific"),
  inference_scale = c("original", "log")
)

Arguments

Details

With model_weights = "inverse_cluster_size", one weighted outcome model supplies predictions for both estimands. With "estimand_specific", two models of the specified method and family are fitted to the same data and formula: the inverse-cluster-size weighted model supplies cATE, and the unweighted model supplies iATE. Cluster sizes are computed from observed membership. Outcome values and trtprob are unchanged, and cluster-size terms are not added to the formula.

Weighted GEE uses geepack::geeglm() with individual fitting weights 1/N_i. Gaussian identity, binomial logit, and Poisson log are supported with independence or exchangeable working correlation.

Weighted LMM and GLMM use WeMix::mix() with one cluster random intercept, conditional individual weights of 1, and cluster weights of 1/N_i. The criterion weights each integrated cluster log-likelihood contribution by 1/N_i; it does not weight individual conditional likelihood contributions before integration. LMM uses a Gaussian identity working model, including when numeric outcomes are binary or counts. GLMM supports Gaussian identity, binomial logit, and Poisson log. The Gaussian GLMM and LMM use the same weighted linear-model backend. Additional WeMix centering is disabled. Full-data and jackknife fits use max_iteration = 100 and nQuad = 13. Gaussian fits use the WeMix linear-model interface (family = NULL) and maximum likelihood. The default unweighted LMM retains its original REML procedure. Nonlinear WeMix models with a random-intercept variance near zero may fail numerically. Unverifiable fits stop; the random intercept is not removed, and no positive variance floor is imposed.

Predictions exclude fitted random effects. Gaussian predictions use the fixed-effects linear predictor. Logit GLMM predictions retain the approximation \operatorname{expit}\{X\widehat\beta / \sqrt{1 + 3\widehat\tau^2/\pi^2}\}, and Poisson GLMM predictions use \exp\{X\widehat\beta + \widehat\tau^2/2\}, where \widehat\tau^2 is the fitted random-intercept variance. GEE uses response-scale marginal predictions.

Both MRS estimands retain the residual correction. The cluster-average estimator gives equal weight to clusters; the individual-average estimator weights clusters by their sizes. The "none" and "inverse_cluster_size" options use one fitted model for both estimands. Every delete-one-cluster sample refits the selected outcome model or models. Jackknife variance is centered at the mean of the deletion estimates, with m-1 degrees of freedom. By default, inference is on the original effect scale. The ICS test uses paired cATE-minus-iATE deletion estimates and always remains on the original effect scale.

With model_weights = "estimand_specific", both models are refitted after deleting the same cluster. The paired deletion contrast is D_{-i}=\widehat\theta_{C,-i}^{w}-\widehat\theta_{I,-i}^{u}, where w and u denote weighted and unweighted fitting. Its jackknife variance is (m-1)/m\sum_i(D_{-i}-\overline D_{(-)})^2, retaining the covariance between the estimates. The procedure requires two model fits per full or deleted sample. Each component retains its existing fitting conventions, including unweighted LMM REML and weighted LMM maximum likelihood. The unweighted component also retains the existing nonparametric fallback and warnings when its outcome-model fit fails. Invalid selected contrasts stop the analysis with the component and full-data or deletion context. No deletion is discarded. This option does not guarantee greater efficiency than common fitting weights.

With inference_scale = "log", let z_{-i}=\log(\widehat\theta_{-i}), where \widehat\theta_{-i} is the RR or OR after deleting cluster i and refitting the selected model. The log-scale standard error is \{(m-1)/m\sum_i(z_{-i}-\overline z_{(-)})^2\}^{1/2}. Confidence limits are \exp\{\log(\widehat\theta)\mathbin{\pm} t_{m-1,1-\alpha/2}\widehat{\mathrm{SE}}_{\log}\}. The two-sided p-value tests a ratio of 1 using \log(\widehat\theta)/\widehat{\mathrm{SE}}_{\log} and a t_{m-1} distribution. This directly jackknifes the log ratios; it does not divide the original-scale standard error by the estimate.

Log inference requires valid arm-specific MRS means for each reported estimand in the full data and every deletion. For "estimand_specific", these are the weighted cATE means and unweighted iATE means. OR requires means strictly between 0 and 1. RR requires finite, positive means; for binary outcomes (all full-data values are 0/1), they must also be at most 1. This outcome classification is retained in all deletion samples. Ratios and log-scale jackknife standard errors must be finite and positive, and alpha must be strictly between 0 and 1. Invalid estimates and confidence limits outside the finite, positive numerical range stop the analysis. Values are not truncated and deletions are not omitted. Positive confidence limits do not guarantee improved small-sample coverage.

Weighted fitting requires complete, finite outcome and covariate data. Weighted binomial outcomes must be 0/1, and Poisson outcomes must be nonnegative integers. Failed weighted fits or invalid predictions stop the analysis and identify the full-data fit or the deleted cluster. Weighted fitting is an additional working-model option; improved efficiency or equality with g-computation is not guaranteed. See the package vignette for examples and a specific GEE comparison.

Value

A list with the following components:

• estimate: A table of point estimates, original-scale standard errors, confidence limits, and p-values. With log inference, only confidence limits and p-values change; the standard-error column remains on the original effect scale.

• m: Number of clusters.

• N: A vector of observed cluster sizes.

• family: The family used for the model.

• model: The method used for the outcome mean model.

• model_weights: The selected outcome-model fitting option.

• inference_scale: The selected inference option.

• log_se: A named numeric vector of log-scale jackknife standard errors for cATE and iATE when log inference is selected; otherwise NULL.

• ics: A named numeric vector containing statistic, the original-scale paired ICS t statistic, and p.value, its two-sided p-value.

• alpha: The specified type-I error rate.

• scale: The requested effect measure.

References

Li F, Tong J, Fang X, Cheng C, Kahan BC, Wang B (2025). Model-Robust Standardization in Cluster-Randomized Trials. Statistics in Medicine, 44, e70270. doi:10.1002/sim.70270.

Lee KM, Harhay MO, Li F (2026). Saturation in G: simple & robust causal inference in cluster randomized trials with informative cluster sizes. arXiv preprint. doi:10.48550/arXiv.2608.28943.

Bailey P, Webb B, Kelley C, Nguyen T, Huo H (2023). WeMix: Weighted Mixed-Effects Models Using Multilevel Pseudo Maximum Likelihood Estimation. R package version 4.0.3. doi:10.32614/CRAN.package.WeMix.

Examples

utils::data("ppact", package = "MRStdCRT")

fit <- MRStdCRT_fit(
  formula = PEGS ~ AGE + FEMALE + comorbid + Dep_OR_Anx + pain_count + PEGS_bl +
    BL_benzo_flag + BL_avg_daily + satisfied_primary + n,
  data     = ppact,
  cluster  = "CLUST",
  trt      = "INTERVENTION",
  trtprob  = NULL,
  method   = "GEE",
  corstr   = "independence",
  scale    = "RR"
)
summary(fit)

weighted_fit <- MRStdCRT_fit(
  PEGS ~ AGE + PEGS_bl, data = ppact,
  cluster = "CLUST", trt = "INTERVENTION",
  method = "GEE", family = gaussian(), corstr = "independence",
  scale = "RD", model_weights = "inverse_cluster_size"
)
weighted_fit$estimate

estimand_specific_fit <- MRStdCRT_fit(
  PEGS ~ AGE + PEGS_bl, data = ppact,
  cluster = "CLUST", trt = "INTERVENTION",
  method = "GEE", family = gaussian(), corstr = "independence",
  scale = "RD", model_weights = "estimand_specific"
)
estimand_specific_fit$estimate
estimand_specific_fit$ics

Example Dataset: Simulated CRT (binary outcome)

Description

A simulated dataset for demonstrating MRStdCRT with a binary outcome. Treatment is assigned at the cluster level and is constant within cluster.

Usage

data(data_sim_binary)

Format

A data frame with the following variables (10 columns):

A

Cluster-level treatment assignment (0/1), constant within cluster.

H1

Cluster-level covariate 1.

H2

Cluster-level covariate 2.

N

Cluster size recorded on each row (repeats within cluster).

X1

Individual-level covariate 1 (numeric).

X2

Individual-level covariate 2 (numeric or binary coded 0/1).

Y

Observed binary outcome (0/1).

Y0

Potential outcome under control (0/1).

Y1

Potential outcome under treatment (0/1).

cluster_id

Cluster identifier (integer or factor), constant within cluster.

Source

Simulated data included with the package for examples.

Examples

data(data_sim_binary)
head(data_sim_binary)
with(data_sim_binary, table(A, Y))

Example Dataset: Simulated CRT (continuous outcome)

Description

A simulated dataset for demonstrating MRStdCRT with a continuous outcome. Treatment is assigned at the cluster level and is constant within cluster.

Usage

data(data_sim_continuous)

Format

A data frame with the following variables (10 columns):

A

Cluster-level treatment assignment (0/1), constant within cluster.

H1

Cluster-level covariate 1.

H2

Cluster-level covariate 2.

N

Cluster size recorded on each row (repeats within cluster).

X1

Individual-level covariate 1 (numeric).

X2

Individual-level covariate 2 (numeric or binary coded 0/1).

Y

Observed continuous outcome.

Y0

Potential outcome under control (continuous).

Y1

Potential outcome under treatment (continuous).

cluster_id

Cluster identifier (integer or factor), constant within cluster.

Source

Simulated data included with the package for examples.

Examples

data(data_sim_continuous)
head(data_sim_continuous)
table(data_sim_continuous$cluster_id)

Example Dataset: PPACT

Description

The Pain Program of Active Coping and Training(PPACT) is a large-scale, mixed methods, cluster-randomized trial (CRT) to compare the effectiveness of an integrated, interdisciplinary program versus usual care in treating patients with chronic pain on long-term opioid treatment (CP-LOT). The primary outcome is the impact of pain (assessed using the PEGS)

Usage

ppact

Format

A data frame with primary outcome, cluster-level, individual level covariates:

SID

Study ID

CLUST

Cluster

INTERVENTION

Study arm

AGE

Patient age at randomization

FEMALE

Participant gender

comorbid

Diagnosis of 2 or more of the chronic medical conditions in 6 month prior to randomization

Dep_OR_Anx

Anxiety and/or depression diagnosis in 6 months prior to randomization

pain_count

Number of different pain types from which participants have diagnoses in 12 months prior to randomization

BL_benzo_flag

Benzodiazepine dispensed in 6 months prior to randomization

BL_avg_daily

Average morphine miligram equivalents dose per day in 6 month prior to randomization

PEGS_bl

PEGS score at baseline

satisfied_primary

Satisfaction with primary care services in prior 3 months

PEGS

PEGS score

n

cluster size

Source

ClinicalTrials.gov: NCT02113592, The manuscript of the study's main outcomes is published in the Annals of Internal Medicine (https://doi.org/10.7326/M21-1436).

Summarize a MRS_obj Fit

Description

Print a concise summary of a model-robust standardization CRT fit, including the c-ATE and i-ATE estimates with SEs and CIs.

Usage

## S3 method for class 'MRS_obj'
summary(object, ...)

Arguments

Value

Invisibly returns the original MRS_obj object, after printing:

• Fitting method and family, with the weighting option and the estimand-specific assignments when selected,

• Number of clusters; per-cluster sizes are available in object$N,

• A four-column table (Estimate, SE, CI, p-value) with rownames c-ATE and i-ATE. The confidence level is 100(1-\alpha)\%, using the alpha stored in the fit. SE is always on the original effect scale. With log inference, CI and p-value use the log scale, with confidence limits transformed back to the ratio scale; the log-scale SE is available separately in object$log_se,

• The NICS test statistic and p-value.