---
title: "Introduction to MRStdCRT"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to MRStdCRT Package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
editor_options:
  markdown:
    wrap: 72
---


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

## 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:

-   The test statistics for testing non-informative cluster sizes,
-   The number of clusters,
-   The cluster sizes, 
-   The outcome regression model used.

User can call the `MRStdCRT_fit` function as follows:

```{r function-signature, eval=FALSE}
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:

-   `formula`: The outcome regression formula, with a single untransformed
    response column and an intercept. Supported terms are untransformed
    covariate main effects and treatment-by-covariate interactions, as
    described below.
-   `data`: The dataset being analyzed.
-   `cluster`: The column name of the cluster identifier. Cluster IDs must
    not be missing.
-   `trt`: The column name of the treatment variable. Values must be
    numeric or logical 0/1, without missing values, and constant within
    each cluster.
-   `trtprob`: A numeric row-level or cluster-level vector of probabilities
    of assignment to treatment, $P(A_i=1)$. Values must be finite, strictly
    between 0 and 1, and constant within each cluster. If `NULL`, a common
    probability is estimated as the proportion
    of clusters assigned to treatment, giving each cluster one vote. This
    probability is re-estimated in each jackknife sample.
    Use known probabilities from the randomization design when available.
-   `method`: The method used for outcome regression model fitting, i.e.
    "GLM", "LMM", "GEE", "GLMM".
-   `family`: The distributional family for the outcome model (default
    is `gaussian(link="identity")`).
-   `corstr`: The working correlation structure used in GEE.
-   `scale`: The scale of estimand including risk difference (RD), risk
    ratio (RR), and odds ratio (OR).
-   `alpha`: The significance level (default is 0.05).

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:

```{r transformed-covariate-formula, eval=FALSE}
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.

```{r treatment-probabilities}
library(MRStdCRT)
library(dplyr)
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.

```{r estimated-treatment-probability}
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.

```{r hypothetical-block-probabilities}
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.

```{r fit-ppact}
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)
## One may also extract specific statistical feature, such as table for point and interval estimates, and standard errors
example$estimate
```