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:

• \(Y_{ij}\) represents the observed outcome for individual \(j\) in cluster \(i\),
• \(A_{i}\) is the treatment assignment for cluster \(i\),
• \(\boldsymbol{X}_{ij}\) denotes the individual-level covariates,
• \(\boldsymbol{H}_i\) denotes the cluster-level covariates,
• \(N_i\) is the number of individuals in cluster \(i\).

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).

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

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

Illustrative example with PPACT data sets

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))

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

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

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)]

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