Introduction to MRStdCRT

Introduction

The MRStdCRT package provides tools for computing the model-robust standardization estimator with jackknife variance estimator for the cluster average treatment effect(c-ATE) and individual average treatment effect(i-ATE) in clustered randomized trials (CRTs).

Model robust standardization

In CRT, assume the cluster size \(N_{i}\) is the natural cluster panel size. The total sample size of the study is \(N=\sum_{i=1}^m N_{i}\). Let \(A_i\in\{0,1\}\) be the randomized cluster-level treatment indicator, with \(A_i=1\) indicating the assignment to the treatment condition and \(A_i=0\) to usual care. The potential outcomes framework and define \(\{Y_{ij}(1),Y_{ij}(0)\}\) as a pair of potential outcomes for each individual \(j \in \{1,\dots, N_i\}\) under the treatment and usual care conditions, respectively. Denote \(\boldsymbol{X}_i = {\boldsymbol{X}_{i1},\dots,\boldsymbol{X}_{iN_i}}^\top\) as the collection of baseline covariates across all individual, and \(\boldsymbol{H}_i\) as the collection of cluster-level covariates. Writing \(f(a,b)\) as a pre-specified contrast function, a general class of weighted average treatment effect in CRTs is defined as \[\Delta_{\omega}=f(\mu_\omega(1),\mu_\omega(0)),\] where the weighted average potential outcome under treatment condition \(A_i=a\) is \[\begin{align*} \mu_\omega(a)=\frac{E\left((\omega_i/N_i)\sum_{j=1}^{N_i}Y_{ij}(a)\right)}{E(\omega_i)}. \end{align*}\]

In this formulation, \(\omega_i\) is a pre-specified cluster-specific weight determining the contribution of each cluster to the target estimand, and can be at most a function of the cluster size \(N_i\), or additional cluster-level covariates \(\boldsymbol{H}_i\). In CRTs, two typical estimands of interest arise from different specifications of \(\omega_i\). First, setting \(\omega_i=1\) gives each cluster equal weight and leads to the cluster-average treatment effect, \(\Delta_C=f(\mu_C(1),\mu_C(0))\), with \[\begin{align*} \mu_C(a)=E\left(\frac{\sum_{j=1}^{N_i}Y_{ij}(a)}{N_i}\right). \end{align*}\] Second, setting \(\omega_i=N_i\) gives equal weight to each individual in the study regardless of their cluster membership and leads to the individual-average treatment effect, \(\Delta_I=f(\mu_I(1),\mu_I(0))\), with \[\begin{align*} \mu_I(a)=\frac{E\left(\sum_{j=1}^{N_i}Y_{ij}(a)\right)}{E(N_i)}. \end{align*}\] Under the general setup, the average potential outcomes can be estimated by \[\begin{align*} \widehat{\mu}_\omega(a)=\sum_{i=1}^m \frac{\omega_i}{\omega_{+}}\left\{\underbrace{\widehat{E}(\overline{Y}_{i}|A_i=a,\boldsymbol{X}_i,\boldsymbol{H}_i,N_i)}_{\text{regression prediction}}+\underbrace{\frac{I(A_i=a)\left(\overline{Y}_i-\widehat{E}(\overline{Y}_{i}|A_i=a,\boldsymbol{X}_i,\boldsymbol{H}_i,N_i)\right)}{\pi(\boldsymbol{X}_i,\boldsymbol{H}_i,N_i)^a\left(1-\pi(\boldsymbol{X}_i,\boldsymbol{H}_i,N_i)\right)^{1-a}}}_{\text{weighted cluster-level residual}}\right\}, \end{align*}\] where \(\omega_{+}=\sum_{i=1}^m \omega_i\) is the sum of weights across clusters, \(\pi(\boldsymbol{X}_i,\boldsymbol{H}_i,N_i)=P(A_i=1|\boldsymbol{X}_i,\boldsymbol{H}_i,N_i)\) is the conditional probability that each cluster is assigned to treatment given baseline information, and \(\widehat{E}(\overline{Y}_{i}|A_i=a,\boldsymbol{X}_i,\boldsymbol{H}_i,N_i)\) is the conditional mean of the cluster average outcome, \(\overline{Y_i}=N_i^{-1}\sum_{j=1}^{N_i}Y_{ij}\), given baseline covariates and cluster size, which could be estimated via any sensible outcome regression model.

Data Structure and Description

In the context of cluster-randomized trials (CRT), we observe the following data vector for each subject \(j\) in cluster \(i\): \(\{Y_{ij}, A_{i}, \boldsymbol{X}_{ij}, \boldsymbol{H}_i, N_i\}\), where:

Different working outcome mean models are proposed based on either individual-level observations or cluster-level summarizes (means), which can be further used to construct the model-robust standardization estimator for either the weighted cluster-averaged treatment effect (c-ATE) or weighted individual-averaged treatment effect (i-ATE).

Syntax

The primary data fitting function is MRStdCRT_fit, which generates a summary for the target estimands, including both c-ATE and i-ATE. In particular, the output provides:

User can call the MRStdCRT_fit function as follows:

MRStdCRT_fit(
  formula, data, cluster, trt, trtprob = rep(0.5, nrow(data)),
  method, family = gaussian(link = "identity"), corstr, scale,
  alpha = 0.05
)

with the following arguments:

The treatment main effect is added automatically. Each interaction must involve treatment and a single covariate, and that covariate must also appear as a main effect. For example, with trt = "A", both y ~ x + A:x and y ~ A * x are supported. Transformations inside the formula, interactions between covariates, higher-order interactions, dot notation (.), no-intercept formulas, and term subtraction are not supported. To include a transformed covariate, first create a separate column and then refer to that column in the formula:

data$x_sq <- data$x^2
formula <- y ~ x + x_sq

summary() displays Estimate, SE, CI, and p-value for both estimands. The displayed confidence level is \(100(1-\alpha)\%\), using the alpha specified in MRStdCRT_fit().

Illustrative example with PPACT data sets

Illustrative Example

In this example, we demonstrate how to use the MRStdCRT_fit function to estimate treatment effects in a CRT using the ppact dataset. The goal is to estimate the cluster-averaged treatment effect (c-ATE) and the individual-averaged treatment effect (i-ATE) using the marginal model fitted by generalized estimating equation (GEE).

Step 1: Prepare the Treatment Assignment Probabilities

Before fitting the model, specify each cluster’s probability of assignment to treatment, \(P(A_i=1)\), regardless of its observed treatment assignment. Use known probabilities from the randomization design when available. For illustration, the following analysis uses a known probability of 0.5 for every cluster. See Section 3.2 in the main manuscript for details on randomization probabilities under other designs.

library(MRStdCRT)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
data(ppact)

prob <- rep(0.5, nrow(ppact))

When estimating a common treatment probability is appropriate, trtprob = NULL calculates the proportion of treated clusters. The following code shows the same calculation explicitly; individuals in larger clusters do not receive extra weight in this calculation.

cluster_assignments <- ppact %>%
  distinct(CLUST, INTERVENTION)

estimated_prob <- rep(
  mean(cluster_assignments$INTERVENTION == 1),
  nrow(ppact)
)

As another example, the following hypothetical treatment probability vector represents a blocked CRT with three blocks, assigned probabilities of 0.3, 0.5, and 0.6, respectively.

clusters <- sort(unique(ppact$CLUST))
n_clusters <- length(clusters) 
block_info <- data.frame(
  CLUST = clusters
) %>%
  mutate(
    # Create a 'block' identifier for each cluster
    block = case_when(
      row_number() <= 25 ~ 1,
      row_number() <= 66 ~ 2,
      TRUE ~ 3
    )
  ) %>%
  mutate(
    # Assign the desired hypothetical treatment probability to each block
    hypothetical_prob = case_when(
      block == 1 ~ 0.3,
      block == 2 ~ 0.5,
      block == 3 ~ 0.6
    )
  )
prob_lookup <- setNames(block_info$hypothetical_prob, block_info$CLUST)
probs_vector <- prob_lookup[as.character(ppact$CLUST)]

Step 2: Fit the MRStdCRT_fit Model

The argument trtprob accepts a vector with one entry per data row or one entry per cluster. For a common known treatment probability of 0.5, use rep(0.5, nrow(data)); a scalar value is not supported. Row-level vectors must follow the order of the input data, with the same probability for every individual in a cluster. An unnamed cluster-level vector follows the order of unique(data[[cluster]]); a named vector uses cluster IDs as its names, with unique names matching all cluster IDs exactly. The example below uses the known probability vector prob defined above.

example <- 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 = prob,
  method = "GEE",
  corstr = "independence",
  scale = "RR"
)
## To view the summary, use the following command
summary(example)
## 
## Model-robust Standardization
## =========================================
##   Method   : GEE
##   Family   : gaussian (link = identity)
##   Clusters : 106
##   Scale    : Risk ratio
## 
## Estimates:
##       Estimate Std. Error         95% CI p-value
## c-ATE    0.907      0.028 (0.852, 0.962) 0.001**
## i-ATE    0.926      0.024 (0.879, 0.973) 0.002**
## 
## Test for no informative cluster size:
##   Statistic: -1.7204
##   p-value  : 0.0883
## One may also extract specific statistical feature, such as table for point and interval estimates, and standard errors
example$estimate
##       Estimate Std. Error  CI lower  CI upper     p-value
## cATE 0.9070587 0.02760612 0.8523209 0.9617965 0.001063919
## iATE 0.9262683 0.02369230 0.8792909 0.9732458 0.002393551