# Conformal causal inference for cluster randomized trials: model-robust
inference without asymptotic approximations

Bingkai Wang1 Fan Li2,3 and Mengxin Yu4
1Department of Biostatistics  University of Michigan  Ann Arbor  MI  USA
2Department of Biostatistics  Yale School of Public Health  New Haven  CT  USA
3Center for Methods in Implementation and Prevention Science  Yale School of
Public Health
New Haven  CT  USA
4Department of Statistics and Data Science  University of Pennsylvania
Philadelphia  PA  USA

###### Abstract

In the analysis of cluster randomized trials, two typical features are that
individuals within a cluster are correlated and that the total number of
clusters can sometimes be limited. While model-robust treatment effect
estimators have been recently developed, their asymptotic theory requires the
number of clusters to approach infinity, and one often has to empirically
assess the applicability of those methods in finite samples. To address this
challenge, we propose a conformal causal inference framework that achieves the
target coverage probability of treatment effects in finite samples without the
need for asymptotic approximations. Meanwhile, we prove that this framework is
compatible with arbitrary working models, including machine learning
algorithms leveraging baseline covariates, possesses robustness against
arbitrary misspecification of working models, and accommodates a variety of
within-cluster correlations. Under this framework, we offer efficient
algorithms to make inferences on treatment effects at both the cluster and
individual levels, applicable to user-specified covariate subgroups and two
types of test data. Finally, we demonstrate our methods via simulations and a
real data application based on a cluster randomized trial for treating chronic
pain.

Keywords: Cluster randomized experiments; conformal prediction; machine
learning; individual-level treatment effect; cluster-level treatment effect;
finite-sample coverage.

##  1 Introduction

Cluster randomized trials (CRTs, Murray et al.,, 1998) represent a type of
randomized trials where the treatment is assigned at the cluster level. The
clusters can be schools, hospitals, or villages, and the associated treatment
can be new healthcare practices or policy innovations, among others. In recent
years, CRTs have become increasingly popular in pragmatic clinical trials to
evaluate the treatment effect in routine practice conditions. Compared to
trials with individual randomization, CRTs offer enhanced cost-effectiveness,
effectively mitigate the risk of treatment contamination, and facilitate the
evaluation of interventions that are naturally implemented at the cluster
level.

In analyzing CRTs, there has been an increasing interest in pursuing model-
robust inference in lieu of conventional model-based inference, aiming to
provide stronger statistical guarantees for estimating targeted causal
estimands (Su and Ding,, 2021; Balzer et al.,, 2023; Wang et al.,, 2023).
Briefly, _model robustness_ means that the validity of inference does not rely
on the correct specification of the selected working model. For example, when
a linear mixed model is used to construct a confidence interval for the
average treatment effect, we aim to achieve the target coverage probability
even when the working model is arbitrarily misspecified (Wang et al.,, 2021).
Model robustness is a desirable property since the true and potentially
complex data-generating distribution is unknown, and thus, any working model
can be misspecified. For example, the outcome may be a nonlinear function of
covariates, follow a non-normal distribution, and have a complex covariate-
dependent correlation structure. Because these complexities cannot be ruled
out in practice, it is almost impossible to pre-specify a working model that
perfectly captures the true data-generating process, thereby rendering model-
robust inference attractive for improving the rigor and precision in analyzing
CRTs.

For model-robust inference in CRTs, model-assisted regression and permutation
inference represent two typical approaches, but they both have important
limitations. Model-assisted regression requires the asymptotic theory assuming
the number of clusters (i.e., the number of independent units, due to unknown
intra-cluster correlation) approaches infinity (Balzer et al.,, 2016; Su and
Ding,, 2021; Wang et al.,, 2021; Balzer et al.,, 2023; Wang et al.,, 2023).
Nonetheless, its practical usefulness for CRTs may be constrained by the
limited number of clusters typically available in such studies. For an example
context, Ivers et al., (2011) reviewed 300 CRTs published between 2000 and
2008 across disciplines, and found that at least 50% randomized no more than
21 clusters, which can be too small for the asymptotic results to apply.
Alternatively, the permutation test can achieve model-robust finite-sample
hypothesis testing for the sharp null (Small et al.,, 2008; Zhang et al.,,
2012; Ding and Keele,, 2018). However, the sharp null may be overly
restrictive—that is, even when a permutation test rejects the sharp null, the
average treatment effect may still be zero due to treatment effect
heterogeneity. Besides, the computation for permutation-based confidence
intervals is often non-trivial (Rabideau and Wang,, 2021). To this end,
improvement for model-robust inference in CRTs that can achieve finite-sample
validity with minimal assumptions and efficiently leverage covariate
information is much needed.

In this article, we integrate conformal prediction, as conceptualized by Vovk
et al., (2005), to perform model-robust causal inference for CRTs. In the
prediction context, conformal prediction aims to construct a covariate-based
interval (using an arbitrary working model) such that the outcome resides in
this interval with a desired probability, e.g., 95%. Remarkably, this coverage
probability is model-robust even in finite samples, provided that the data
distribution is exchangeable among individuals. Given that its model
robustness is independent of any asymptotic regimes, conformal prediction
effectively solves the analytical challenges with a limited number of
independent units, making it a promising paradigm for analyzing CRTs. Although
conformal prediction has been connected with causal inference for non-
clustered data (Lei and Candès,, 2021; Yang et al.,, 2022; Qiu et al.,, 2022;
Yin et al.,, 2022; Alaa et al.,, 2023; Jin et al.,, 2023), conformal causal
inference for multilevel data has not been previously investigated and poses
new challenges. First, CRT data may not be exchangeable at the individual
level because the intra-cluster correlation can be much larger than the inter-
cluster correlation, causing non-exchangeability for individuals from
different clusters. As a result, the existing theory of conformal causal
inference cannot be directly applied without essential modifications. Second,
the multilevel data structure of CRTs implies more than one treatment effect
at different levels, whose magnitude and interpretation can differ (Kahan et
al.,, 2023), and conformal causal inference requires adaptation to target each
treatment effect. Third, while existing results for conformal causal inference
focus on the marginal interpretation over the entire study population,
inferring treatment effects among certain subpopulations can lead to more
informative and targeted interpretations, representing a valuable extension
for inference in CRTs.

The principal contribution of this article is to propose a conformal causal
inference framework tailored for CRTs that address the aforementioned
challenges. Specifically, given any working model, we provide algorithms to
construct prediction intervals for both the cluster-level and individual-level
treatment effect and prove their finite-sample model robustness under minimal
assumptions. The cluster-level treatment effect assesses the impact of
treatment on a new cluster as a whole, while the individual-level treatment
effect focuses on the impact experienced by any given individual within a new
cluster: both treatment effects can be meaningful for decision-making
depending on the scientific questions. Furthermore, we extend our methods to
perform subgroup analysis, that is, achieving the desired coverage probability
conditioning on any set with a positive measure in the covariate data space.
Throughout, the proposed method allows for the integration of machine learning
algorithms as a base learner for the outcome given baseline covariates, in
order to narrow the prediction intervals for the target treatment effects
without compromising validity in finite samples.

##  2 Notation, assumptions, and targets of inference

###  2.1 Potential outcomes framework

We consider a CRT with $m$ clusters, and each cluster $i,\ i=1,\dots,m$,
includes $N_{i}$ individuals. At the cluster level, we define $A_{i}$ as a
binary treatment indicator ($A_{i}=1$ for treatment and $A_{i}=0$ for control)
and $R_{i}$ as a vector of cluster-level covariates, such as the geographical
location of the cluster. For each individual $j=1,\dots,N_{i}$ in cluster $i$,
we denote $Y_{ij}$ as their observed outcome, $X_{ij}$ as a vector of
individual-level baseline covariates. The observed data for each cluster are
$O_{i}=\\{(Y_{ij},A_{i},X_{ij},R_{i}):j=1,\dots,N_{i}\\}$. To proceed, we
adopt the potential outcomes framework and define $Y_{ij}(a)$ as the potential
outcome for individual $j$ of cluster $i$ had cluster $i$ been assigned
$A_{i}=a$ for $a=0,1$. We assume consistency such that
$Y_{ij}=A_{i}Y_{ij}(1)+(1-A_{i})Y_{ij}(0)$. Denoting the complete data vector
for each cluster as
$W_{i}=\\{(Y_{ij}(1),Y_{ij}(0),X_{ij},R_{i}):j=1,\dots,N_{i}\\}$, we make the
following structural assumptions on $(W_{1},\dots,W_{m})$ and the assignment
vector $(A_{1},\dots,A_{m})$.

######  Assumption 1 (_Random sampling of clusters_).

Each $W_{i}$ is an independent, identically distributed sample from an unknown
distribution $\mathcal{P}^{W}$ on
$W=\\{(Y_{\bullet,j}(1),Y_{\bullet,j}(0),X_{\bullet,j},R):j=1,\dots,N\\}$,
where the subscript ‘${\bullet,j}$’ represents the $j$-th individual in $W$.

######  Assumption 2 (_Cluster randomization_).

Each $A_{i}$ is an independent, identically distributed sample from a
Bernoulli distribution $\mathcal{P}^{A}$ with $P(A=1)=\pi$ for a known
constant $\pi\in(0,1)$. Furthermore, $\mathcal{P}^{A}$ is independent of
$\mathcal{P}^{W}$.

######  Assumption 3 (_Within-cluster exchangeability_).

For any permutation map $\sigma$ on the index set $\\{1,\dots,N\\}$, the
vector $\\{W_{\bullet,1},\dots,W_{\bullet,N}\\}$ has the same distribution as
$\\{W_{\bullet,\sigma(1)},\dots,W_{\bullet,\sigma(N)}\\}$ conditioning on $N$,
where $W_{\bullet,j}=(Y_{\bullet,j}(1),Y_{\bullet,j}(0),X_{\bullet,j},R)$.

Assumption 1 states that each observed cluster is an independent draw from the
same distribution and is standard for drawing inference under a sampling-based
framework. To accommodate the difference in the dimension of $W_{i}$ across
clusters, we can decompose
$\mathcal{P}^{W}=\mathcal{P}^{W|N}\times\mathcal{P}^{N}$ such that $N_{i}$ is
first sampled from $\mathcal{P}^{N}$ and $W_{i}$ is then generated following
$\mathcal{P}^{W|N}$. Assumption 2 often holds by design. Assumption 3 states
that the joint distribution of individuals in a cluster is invariant to
permutations of the index of individuals. Importantly, this assumption still
allows heterogeneous outcome correlation after conditioning on covariates,
which means the conditional correlation of $Y_{\bullet,j}$ and
$Y_{\bullet,j^{\prime}}$ can arbitrarily depend on $X_{\bullet,j}$ and
$X_{\bullet,j^{\prime}}$ through a shared covariance function across
$(j,j^{\prime})$ pairs; see Crespi et al., (2009) for an example data-
generating process in the context of second-order generalized estimating
equations. As a result, Assumption 3 encompasses a wide variety of within-
cluster correlation structures and is not to be confused with the restrictive
random-intercept model assumption (which implies an overly simplistic constant
intracluster correlation and is a special case of Assumption 3). Of note,
Assumption 3 is weaker than those in Balzer et al., (2016) and Benitez et al.,
(2023), which further assumes conditional independence among individual-level
observations within each cluster for asymptotic inference.

###  2.2 Conventional target of inference goals in CRTs

Conventionally, statistical inference for CRTs focuses on the average
treatment effect. For example, when the interest lies in the cluster-average
treatment effect $\Delta_{C}=E[\overline{Y}(1)-\overline{Y}(0)]$ with
$\overline{Y}(a)=\frac{1}{N}\sum_{j=1}^{N}Y_{\bullet,j}(a)$, one important
inferential target is to construct a confidence interval
$\widehat{\textup{CI}}_{m}$ satisfying
$P\left(\Delta_{C}\in\widehat{\textup{CI}}_{m}\right)=1-\alpha$ with $\alpha$
being a pre-specified type I error rate. Often in practice,
$\widehat{\textup{CI}}_{m}$ leverages baseline covariates to sharpen the
interval and reduce the type II error rate. Despite its simple form, this
objective may be difficult without parametric modeling assumptions, such as
linear mixed models or Bayesian methods. However, parametric assumptions
implicitly put strong constraints on the space of possible data-generating
distributions and, therefore, can be violated under non-linearity or non-
normality, among others. To overcome this challenge, Su and Ding, (2021) and
Wang et al., (2023) developed asymptotic theory that guarantees a weaker
result,
$\lim_{m\rightarrow\infty}P\left(\Delta_{C}\in\widehat{\textup{CI}}_{m}\right)=1-\alpha$,
which states that the target coverage probability of a model-assisted
confidence interval is achievable as the number of clusters approaches
infinity, even under arbitrary working model misspecification. Nonetheless,
this modified inferential target may not be fully satisfactory since the
number of clusters $m$ is often limited in CRTs.

###  2.3 Conformal causal inference in CRTs

For finite-sample inference that is robust to model misspecification, we
propose an alternative inference target defined as follows. Given any
$\alpha\in(0,1)$, we aim to construct a conformal interval
$\widehat{\textup{C}}_{C}(\overline{B})\subset\mathbb{R}$ for the cluster-
level treatment effect such that

$P\bigg{\\{}\overline{Y}(1)-\overline{Y}(0)\in\widehat{\textup{C}}_{C}(\overline{B})\bigg{\\}}\geq
1-\alpha,$ (1)

where $\overline{B}=(\overline{X},R,N)$ is the cluster-level covariate
information with $\overline{X}=\frac{1}{N}\sum_{j=1}^{N_{i}}X_{\bullet,j}$
being the cluster-aggregate covariates. Equation (1) has the following
interpretation: for any cluster sampled from $\mathcal{P}^{W}$, its cluster-
level treatment effect, $\overline{Y}(1)-\overline{Y}(0)$, is guaranteed to be
contained in $\widehat{\textup{C}}_{C}(\overline{B})$ with probability at
least $1-\alpha$, without requiring any asymptotic regimes. Therefore,
$\widehat{\textup{C}}_{C}(\overline{B})$ naturally quantifies the cluster-
level treatment effect by characterizing the most likely range of
$\overline{Y}(1)-\overline{Y}(0)$. Distinct from the confidence interval
$\widehat{\textup{CI}}_{m}$, $\widehat{\textup{C}}_{C}(\overline{B})$ is
covariate-based and targets the random variable
$\overline{Y}(1)-\overline{Y}(0)$. Then,
$\widehat{\textup{C}}_{C}(\overline{B})$ essentially provides the
$(\alpha/2,1-\alpha/2)$-quantile for the distribution of
$\overline{Y}(1)-\overline{Y}(0)$. Of note, while the definition in (1) is
directly motivated by the cluster-average estimand, the same framework can be
directly applied to other cluster summaries such as cluster totals (Su and
Ding,, 2021).

Our target of inference also naturally extends to the construction of a
conformal interval $\widehat{\textup{C}}_{I}(B)$ for the individual-level
treatment effect such that

$P\bigg{\\{}Y(1)-Y(0)\in\widehat{\textup{C}}_{I}(B)\bigg{\\}}\geq 1-\alpha.$
(2)

Here, $B=(X,R,N)$ is the individual-level covariate data, and
$\\{Y(1),Y(0),X\\}$ can represent any individual
$\\{Y_{\bullet,j}(1),Y_{\bullet,j}(0),X_{\bullet,j}\\}$ since they are
identically distributed given $N$ by Assumption 3. The interpretation of
Equation (2) is similar to Equation (1): for any cluster sampled from
$\mathcal{P}^{W}$, the individual-level treatment effect $Y(1)-Y(0)$ is
contained in $\widehat{\textup{C}}_{I}(B)$ with probability at least
$1-\alpha$. Similarly, $\widehat{\textup{C}}_{I}(B)$ characterizes the
$(\alpha/2,1-\alpha/2)$-quantile of $Y(1)-Y(0)$ and describes a personalized
range. In contrast to the cluster-level treatment effect, which only pertains
to independent cluster-level data, accomplishing the objectives in Equation
(2) is more challenging due to the unknown intracluster correlations between
observations inherent in CRTs.

Beyond (1) and (2), we further study conformal causal inference on pre-defined
covariate subgroups. Specifically, let $\Omega_{C}$ and $\Omega_{I}$ be
arbitrary sets with positive measures in the support of $\overline{B}$ and
$B$, respectively. For example, if $X_{1}$ represents age, we can define
$\Omega_{C}=\\{\overline{B}:\overline{X}_{1}\geq 70\\}$ and
$\Omega_{I}=\\{B:X_{1}\geq 70\\}$ to focus on the treatment effects on a
specific age group. Given $\Omega_{C}$ and $\Omega_{I}$, we aim to construct
conformal intervals $\widetilde{\textup{C}}_{C}(\overline{B})$ and
$\widetilde{\textup{C}}_{I}(B)$ such that

$\displaystyle
P\bigg{\\{}\overline{Y}(1)-\overline{Y}(0)\in\widetilde{\textup{C}}_{C}(\overline{B})\bigg{|}\overline{B}\in\Omega_{C}\bigg{\\}}$
$\displaystyle\geq 1-\alpha,$ (3) $\displaystyle
P\bigg{\\{}Y(1)-Y(0)\in\widetilde{\textup{C}}_{I}(B)\bigg{|}B\in\Omega_{I}\bigg{\\}}$
$\displaystyle\geq 1-\alpha.$ (4)

Equations (3) and (4) characterize “local coverage” that is distinct from
“marginal coverage” characterized by Equations (1) and (2). With local
coverage, the resulting conformal intervals provide more targeted information
for specific sub-populations of interest (rather than the entire study
population) and carry a finer resolution (in the same spirit of subgroup
analyses). The marginal coverage definition is then viewed as a special case
of local coverage by setting $\Omega_{C}$ and $\Omega_{I}$ to be the entire
covariate support. Notably, local coverage differs from conditional coverage,
where the probability measures are further conditioned on $\overline{B}$ or
$B$. Compared to local coverage, conditional coverage is a much stronger
objective, whereas it has been proved impossible for continuous covariates
(Lei and Wasserman,, 2014).

Before moving to the technical development, we introduce some additional
notation. Let $\overline{Y}=A\overline{Y}(1)+(1-A)\overline{Y}(0)$, $I\\{G\\}$
denote the indicator function of event $G$, $||\cdot||_{2}$ is the $\ell_{2}$
norm, $\delta_{s}$ as a point mass at $s$ in the Euclidean space, and
$\delta_{+\infty}$ is a point mass at infinity.

##  3 Conformal causal inference for cluster-level treatment effects

###  3.1 Inference for an observed cluster

For inferring cluster-level treatment effects, we aim to construct a conformal
interval $\widetilde{\textup{C}}_{C}$ given a test point
$\overline{B}_{\textup{test}}$, i.e., the cluster-aggregate covariates of a
new cluster sampled from the target population. We first consider a basic
scenario, where we observe the complete information
$\overline{O}_{\textup{test}}=(\overline{Y}_{\textup{test}},A_{\textup{test}},\overline{B}_{\textup{test}})$
for the test point. In other words, this test point corresponds to an
“observed cluster”, i.e., its current intervention $A_{\textup{test}}$ (likely
not randomized, often $A_{\textup{test}}=0$ and not included in the CRT) and
current average outcome $\overline{Y}(A_{\textup{test}})$ are also recorded.
Such additional data beyond $\overline{B}_{\textup{test}}$ will simplify the
inference procedure since one of the two potential outcomes
$\\{\overline{Y}_{\textup{test}}(1),\overline{Y}_{\textup{test}}(0)\\}$ is
directly observed, and the conformal causal inference only requires
constructing the conformal interval for the unobserved average potential
outcome.

Under this basic scenario, Algorithm 1 outlines the steps to compute the
conformal interval $\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})$.
In the input phase, the prediction model $f_{a},a\in\\{0,1\\}$ can be an
arbitrary map from covariates to the outcome, e.g., a linear model, or random
forest (Breiman,, 2001). Different choices of $f_{a}$ will not impact the
coverage validity, but can result in conformal intervals with different
lengths and thus affect precision. In addition, we need to specify the
covariate subgroup of interest, e.g., $\Omega_{C}$ as the entire covariate
support or
$\Omega_{C}=\\{\overline{B}_{i}:\overline{X}_{i1}>70,R_{i1}=\textup{``Rural''}\\}$
as a sub-population, and a level $\alpha$, e.g., $\alpha=0.1$. While a more
restrictive $\Omega_{C}$ can target a smaller sub-population, the resulting
conformal interval may be wider due to fewer data points; we further discuss
this trade-off in Section 7.

Algorithm 1 Computing the conformal interval
$\widetilde{\textup{C}}_{C}(\overline{O})$ for cluster-level treatment
effects.

Input: Cluster-level data
$\\{(\overline{Y}_{i},A_{i},\overline{B}_{i}):i=1,\dots,m\\}$, a test point
$\overline{O}_{\textup{test}}=(\overline{Y}_{\textup{test}},A_{\textup{test}},\overline{B}_{\textup{test}})$,
an arbitrary prediction model $f_{a}(\overline{B})$ for $\overline{Y}(a)$,
$a\in\\{0,1\\},$ a covariate subgroup of interest $\Omega_{C}$, and level
$\alpha$.

Step 1 (Constructing the conformal interval
$\widetilde{\textup{C}}_{C,a}(\overline{B})$ for $\overline{Y}(a)$.)

For $a=0,1$,

1\. Randomly split the arm-$a$ covariate subgroup data
$\\{(\overline{Y}_{i},\overline{B}_{i}):i=1,\dots,m,A_{i}=a,\overline{B}_{i}\in\Omega_{C}\\}$
into a training fold $\mathcal{O}_{tr}(a)$ and a calibration fold
$\mathcal{O}_{ca}(a)$ with index set $\mathcal{I}_{ca}(a)$.

2\. Train the prediction model $f_{a}(\overline{B})$ using the training fold
$\mathcal{O}_{tr}(a)$, and obtain the estimated model
$\widehat{f}_{a}(\overline{B})$.

3\. For each $i\in\mathcal{I}_{ca}(a)$, compute the non-conformity score
$s(\overline{B}_{i},\overline{Y}_{i})=|\overline{Y}_{i}-\widehat{f}_{a}(\overline{B}_{i})|$.

4\. Compute the $1-\alpha$ quantile $\widehat{q}_{1-\alpha}(a)$ of the
distribution

$\widehat{F}=\frac{1}{|\mathcal{I}_{ca}(a)|+1}\left\\{\sum_{i\in\mathcal{I}_{ca}(a)}\delta_{s(\overline{B}_{i},\overline{Y}_{i})}+\delta_{+\infty}\right\\}.$

5\. Obtain
$\widetilde{\textup{C}}_{C,a}(\overline{B})=\\{y\in\mathbb{R}:|y-\widehat{f}_{a}(\overline{B})|\leq\widehat{q}_{1-\alpha}(a)\\}$.

Step 2 (Constructing the conformal interval
$\widetilde{\textup{C}}_{C}(\overline{O})$ for
$\overline{Y}(1)-\overline{Y}(0)$.)

If $A_{\textup{test}}=1,$ then set
$\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})=\overline{Y}_{\textup{test}}-\widetilde{\textup{C}}_{C,0}(\overline{B}_{\textup{test}})$;

if $A_{\textup{test}}=0,$ then set
$\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})=\widetilde{\textup{C}}_{C,1}(\overline{B}_{\textup{test}})-\overline{Y}_{\textup{test}}$.

Output: $\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})$.

Given the input, we take two steps to construct the conformal interval. In the
first step, we apply the split conformal prediction (Papadopoulos et al.,,
2002) to construct a conformal interval for one potential outcome
$\overline{Y}(a)$, based on arm-$a$ data in the covariate subgroup
characterized by $\Omega_{C}$. In split conformal prediction, the data are
randomly partitioned into two halves, one used to train the prediction model
$f_{a},$ $a\in\\{0,1\\},$ (Steps 1.1-1.2) and the other used to construct the
conformal interval (Steps 1.3-1.5). Compared to alternative methods such as
full conformal prediction (Lei et al.,, 2018) and cross-validated conformal
prediction (Barber et al.,, 2021), split conformal prediction uses random data
partition to substantially reduce computational complexity. It therefore
embraces practicality, despite its less efficient use of the data (Tibshirani
et al.,, 2019). enGiven the model fit $\widehat{f}_{a}$, Step 1.3 computes the
non-conformity score $s(b,y)$, which is the absolute value of prediction error
on calibration data. Intuitively, a large non-conformity score indicates an
abnormal data point (i.e., it does not “conform”) with respect to
$\widehat{f}_{a}$. In Step 1.5, we construct
$\widetilde{\textup{C}}_{C,a}(\overline{B})$ by excluding $y$ with a large
non-conformity score $s(\overline{B},y)$. The threshold,
$\widehat{q}_{1-\alpha}$, as the exclusion criterion is defined as the
$(1-\alpha)$-quantile of the empirical distribution on the non-conformity
score (plus a point mass at infinity). The underlying idea of this threshold
comes from a neat observation: for exchangeable variables
$(V_{1},\dots,V_{n},V_{n+1})$, $V_{n+1}$ is smaller than the
$(1-\alpha)$-quantile of their empirical distribution with probability at
least $1-\alpha$. Following this logic, the constructed
$\widetilde{\textup{C}}_{C,a}(\overline{B})$ has probability no less than
$1-\alpha$ to contain $\overline{Y}(a)$. With
$\widetilde{\textup{C}}_{C,a}(\overline{B})$ obtained from Step 1, we proceed
to Step 2, where
$\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})=(-1)^{A_{\textup{test}}+1}\\{\overline{Y}_{\textup{test}}-\widetilde{\textup{C}}_{C,1-A_{\textup{test}}}(\overline{B}_{\textup{test}})\\}$
is the output. In other words, the final conformal interval is a contrast
between the observed potential outcome and the interval constructed for the
unobserved potential outcome.

Theorem 1 proves that
$\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})$ achieves finite-
sample coverage guarantee for the cluster-level treatment effect, without
requiring Assumption 3 on within-cluster exchangeability.

######  Theorem 1.

Under Assumptions 1-2 and assuming that
$\overline{O}_{\textup{test}}=(\overline{Y}_{\textup{test}},A_{\textup{test}},\overline{B}_{\textup{test}})$
is an independent sample from the distribution $\mathcal{P}^{\overline{O}}$
induced by $\mathcal{P}^{W}\times\widetilde{\mathcal{P}}^{A|W}$ and
$\overline{Y}_{\textup{test}}=A_{\textup{test}}\overline{Y}_{\textup{test}}(1)+(1-A_{\textup{test}})\overline{Y}_{\textup{test}}(0)$,
where $\widetilde{\mathcal{P}}^{A|W}$ is an arbitrary unknown distribution for
$A_{\textup{test}}$. Then, the conformal interval
$\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})$ output by Algorithm
1 satisfies

$P\bigg{\\{}\overline{Y}_{\textup{test}}(1)-\overline{Y}_{\textup{test}}(0)\in\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})\bigg{|}\overline{B}_{\textup{test}}\in\Omega_{C}\bigg{\\}}\geq
1-2\alpha$ (5)

for any set $\Omega_{C}$ in the support of $\overline{B}_{\textup{test}}$ with
a positive measure.

Theorem 1 is a direct generalization of the conformal prediction theory to
conformal causal inference on covariate subgroups. Since we observe
$\overline{Y}_{\textup{test}}$ and $A_{\textup{test}}$, the valid coverage for
the treatment effect is straightforward once we achieve valid coverage for the
potential outcome. However, since we need to account for arbitrary
distribution shift on $A_{\textup{test}}$ from $\mathcal{P}^{A}$ to
$\widetilde{\mathcal{P}}^{A|W}$, the non-coverage rate $\alpha$ is doubled in
Equation (5) by applying the union bound. If $\widetilde{P}^{A|W}$ is
independent of $W$ (e.g., $P(A_{\textup{test}}=1)=0$ when studying a new
intervention or $A_{\textup{test}}$ is randomized), then the coverage
probability achieves $1-\alpha$ in Equation (5). For local coverage, we
consider a simple yet effective approach that only clusters with
$\overline{B}_{i}\in\Omega_{C}$ are included. Technically, it changes the
covariate distribution but not the outcome distribution given covariates,
thereby not changing the coverage result. Recently, Hore and Barber, (2023)
proposed an approach to construct the conformal interval that is independent
of $\Omega_{C}$, whereas this uniformity property can lead to poor coverage
(e.g., zero) in many covariate subgroups and, more importantly, discrete
covariates cannot be directly handled by their localization kernels.
Furthermore, it is worth mentioning two important special cases of Theorem 1.
When $P(\overline{B}\in\Omega_{C})=1$, Equation (5) implies marginal coverage
in Equation (1). Setting $N_{i}=1$, Theorem 1 also applies to conformal causal
inference for individually randomized trials.

In practice, how to choose $\alpha$ and the size of the training fold depends
on the number of clusters $m$. To obtain a non-trivial conformal interval such
that $\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})\neq\mathbb{R}$,
the definition of $\widehat{q}_{1-\alpha}(a)$ requires that
$\alpha\geq|\mathcal{I}_{ca}(a)|^{-1}$. For example, if we set $\alpha=0.1$,
then the calibration fold for each arm should contain at least 10 clusters,
and the rest $m-20$ clusters can be used as the training fold. For stability
of model fitting, we recommend using at least 20 clusters for training, while
fewer clusters are also acceptable if parsimonious models such as linear
regression are used. When the number of clusters is small, e.g., $m=20$,
Algorithm 1 based on split conformal prediction cannot construct non-trivial
conformal intervals with $90\%$ coverage probability. However, this goal is
achievable with full conformal prediction or Jackknife+, but demands a
substantially higher computation cost (Barber et al.,, 2021).

###  3.2 Inference based on cluster-level covariates

When the target of inference concerns a new cluster that only has covariate
information $\overline{B}_{\textup{test}}$ (e.g., a cluster having not taken
either treatment or control studied in the current CRT), a direct approach
based on Algorithm 1 is to combine $\widetilde{\textup{C}}_{C,1}$ and
$\widetilde{\textup{C}}_{C,0}$, for which Corollary 1 characterizes its local
coverage property.

######  Corollary 1.

Under Assumptions 1-2 and assuming that
$(\overline{Y}_{\textup{test}}(1),\overline{Y}_{\textup{test}}(0),\overline{B}_{\textup{test}})$
is an independent sample from the distribution $\mathcal{P}^{W}$. Then,
$\widetilde{\textup{C}}_{C,1}(\overline{B}_{\textup{test}})$ and
$\widetilde{\textup{C}}_{C,0}(\overline{B}_{\textup{test}})$ output by
Algorithm 1 satisfy

$P\bigg{\\{}\overline{Y}_{\textup{test}}(1)-\overline{Y}_{\textup{test}}(0)\in\widetilde{\textup{C}}_{C,1}(\overline{B}_{\textup{test}})-\widetilde{\textup{C}}_{C,0}(\overline{B}_{\textup{test}})\bigg{|}\overline{B}_{\textup{test}}\in\Omega_{C}\bigg{\\}}\geq
1-\alpha$ (6)

for any set $\Omega_{C}$ in the support of $\overline{B}_{\textup{test}}$ with
a positive measure.

Compared to Equation (5), Equation (6) provides a stronger coverage
probability $1-\alpha$, but the length of the conformal interval can be
approximately doubled since we no longer observe
$\overline{Y}_{\textup{test}}$ and $A_{\textup{test}}$. As a result, the
conformal interval is based on $\overline{B}_{\textup{test}}$ tends to be less
informative than $\widetilde{\textup{C}}_{C}(\overline{O}_{\textup{test}})$, a
natural result of less observed information.

Beyond this direct approach, Lei and Candès, (2021) provided a more flexible
nested approach to construct
$\widetilde{\textup{C}}_{C}(\overline{B}_{\textup{test}})$. This method first
computes $\overline{C}_{i}$, the $(1-\alpha)$-conformal interval for
$\overline{Y}(1)-\overline{Y}(0)$, using $\widetilde{\textup{C}}_{C,a}$ from
Algorithm 1, and then run split conformal prediction again on
$(\overline{C}_{i},\overline{B}_{i})$ with level $\gamma$ and prediction
models $(m^{L},m^{R})$. Since we run split conformal prediction twice, the
resulting conformal interval
$\widetilde{\textup{C}}_{C}(\overline{B}_{\textup{test}})$ has non-coverage
probability up to $\alpha+\gamma$. For completeness, we provide the detailed
algorithm and its theoretical guarantee in the Supporting Information.
Finally, we remark that the nested approach becomes similar to the direct
approach with $\gamma$ set equal to $\alpha$. In practice, both the direct and
nested approaches tend to be conservative since
$\\{\overline{Y}(1),\overline{Y}(0)\\}$ are not observed, but the nested
approach can alleviate conservativeness by specifying a larger $\gamma$ and
trading the theoretical guarantee for better empirical performance.

##  4 Conformal causal inference for individual-level treatment effects

So far, we have extended conformal causal inference for the cluster-level
treatment effects by treating cluster-level data as independent units.
However, for individual-level treatment effects, the existing theory for
conformal inference cannot be directly applied since the individual-level data
are clustered. This feature presents apparent challenges in obtaining
conformal intervals with either marginal or local coverage guarantees. Lee et
al., (2023) developed recent results for conformal inference with hierarchical
data, focusing on the setting where individuals in the same cluster are
conditionally independent. We consider this setting overly restrictive because
it rules out any unmeasured source of within-cluster correlations and is not
even compatible with a random-intercept linear mixed model. Below, we provide
the algorithm and theory for conformal causal inference on the individual-
level treatment effect estimand under a more general exchangeability setup
characterized by Assumption 3.

###  4.1 Inference for an observed individual

We first consider a basic setting, where the test point has the complete
information
$O_{\textup{test}}=(Y_{\textup{test}},A_{\textup{test}},B_{\textup{test}})$.
This setting corresponds to an observed individual in an observed cluster of
interest, whose current treatment $A_{\textup{test}}$ may not be randomized.
Similar to our development in Section 3.1, we only need to construct the
conformal interval for the unobserved potential outcome, leading to higher
chances of narrower and more informative conformal intervals.

In Algorithm 2, we show how to construct the conformal interval
$\widetilde{\textup{C}}_{I}(O)$ for $Y(1)-Y(0)$. Compared to Algorithm 1, the
major difference is in Step 1.4: the new empirical distribution function
$\widehat{F}$ only involves individuals meeting the criteria
$B_{ij}\in\Omega_{I}$, and each individual is further weighted by
$\left(\sum_{j=1}^{N_{i}}I\\{B_{ij}\in\Omega_{I}\\}\right)^{-1}$; all other
steps remain similar but now apply to the individual data (rather than cluster
aggregate). Since we target inference on specific covariate subgroups
characterized by $\Omega_{I}$ (which again can be the entire support), we
prove a Lemma in the Supporting Information to show that Assumptions 1-3 holds
when restricting to any subgroup $\Omega_{I}$ with a positive measure. This is
a critical step to preserve the exchangeable structure within clusters and
develop our theoretical results. Furthermore, the clustered data structure is
handled by inverse weighting: with these weights, point masses within a
cluster are symmetric, and the cluster-average point masses are also symmetric
across clusters. Although the technical details are more complex and explained
in the Supporting Information, intuitively, the multilevel symmetry of
$\widehat{F}$ links the test data and calibration data as in the setting of
exchangeable data, based on which we can derive finite-sample robustness. The
theoretical guarantee of the resulting conformal interval is formally stated
in Theorem 2.

Algorithm 2 Computing the conformal interval $\widetilde{\textup{C}}_{I}({O})$
for individual-level treatment effects.

Input: Individual-level data
$\\{(Y_{ij},A_{i},B_{ij}):i=1,\dots,m;j=1,\dots,N_{i}\\}$, a test point
$O_{\textup{test}}=(Y_{\textup{test}},A_{\textup{test}},B_{\textup{test}})$, a
prediction model $f_{a}(B)$ for $Y(a)$, $a\in\\{0,1\\},$ a covariate subgroup
of interest $\Omega_{I}$, and level $\alpha$.

Step 1 (Constructing the conformal interval $\widetilde{\textup{C}}_{I,a}(B)$
for $Y(a)$.)

For $a=0,1$,

1\. Randomly split the arm-$a$ covariate subgroup data
$\\{(Y_{ij},A_{i},B_{ij}):i=1,\dots,m;j=1,\dots,N_{i};A_{i}=a;B_{ij}\in\Omega_{I}\\}$
into a training fold $\mathcal{O}_{tr}(a)$ and a calibration fold
$\mathcal{O}_{ca}(a)$ with index set $\mathcal{I}_{ca}(a)$.

2\. Train the prediction model $f_{a}(B)$ using the training fold
$\mathcal{O}_{tr}(a)$, and obtain the estimated model $\widehat{f}_{a}(B)$.

3\. For each $(i,j)\in\mathcal{I}_{ca}(a)$, compute the non-conformity score
$s(B_{ij},Y_{ij})=|Y_{ij}-\widehat{f}_{a}(B_{ij})|$.

4\. Compute the $1-\alpha$ quantile $\widehat{q}_{1-\alpha}(a)$ of the
distribution

$\widehat{F}=\frac{1}{|\mathcal{I}_{ca}(a)|+1}\left\\{\sum_{i\in\mathcal{I}_{ca}(a)}\frac{1}{\sum_{j=1}^{N_{i}}I\\{B_{ij}\in\Omega_{I}\\}}\sum_{j=1}^{N_{i}}I\\{B_{ij}\in\Omega_{I}\\}\delta_{s(B_{ij},Y_{ij})}+\delta_{+\infty}\right\\}.$

5\. Obtain
$\widetilde{\textup{C}}_{I,a}(B)=\\{y\in\mathbb{R}:|y-\widehat{f}_{a}(B)|\leq\widehat{q}_{1-\alpha}(a)\\}$.

Step 2 (Constructing the conformal interval $\widetilde{\textup{C}}_{I}(O)$
for $Y(1)-Y(0)$.)

If $A_{\textup{test}}=1,$ then set
$\widetilde{\textup{C}}_{I}(O_{\textup{test}})=Y_{\textup{test}}-\widetilde{\textup{C}}_{I,0}(B_{\textup{test}})$;

if $A_{\textup{test}}=0,$ then set
$\widetilde{\textup{C}}_{I}(O_{\textup{test}})=\widetilde{\textup{C}}_{I,1}(B_{\textup{test}})-Y_{\textup{test}}$.

Output: $\widetilde{\textup{C}}_{I}(O_{\textup{test}})$.

######  Theorem 2.

Under Assumptions 1-3 and assuming that
$O_{\textup{test}}=(Y_{\textup{test}},A_{\textup{test}},B_{\textup{test}})$ is
an arbitrary individual in a new cluster independently sampled from
$\mathcal{P}^{W}\times\widetilde{\mathcal{P}}^{A|W}$ with
$Y_{\textup{test}}=A_{\textup{test}}Y_{\textup{test}}(1)+(1-A_{\textup{test}})Y_{\textup{test}}(0)$,
where $\widetilde{\mathcal{P}}^{A|W}$ is an arbitrary unknown distribution for
$A_{\textup{test}}$. Then, the $\widetilde{\textup{C}}_{I}(O_{\textup{test}})$
output by Algorithm 1 satisfies

$P\bigg{\\{}Y_{\textup{test}}(1)-Y_{\textup{test}}(0)\in\widetilde{\textup{C}}_{I}(O_{\textup{test}})\bigg{|}B_{\textup{test}}\in\Omega_{I}\bigg{\\}}\geq
1-2\alpha$ (7)

for any set $\Omega_{I}$ in the support of $B_{\textup{test}}$ with a positive
measure.

Theorem 2 is the counterpart of Theorem 1 characterizing the local or marginal
coverage property for individual-level treatment effects, depending on the
choice of $\Omega_{I}$. Due to the distribution shift on $A_{\textup{test}}$,
the coverage probability is $1-2\alpha$ instead of $1-\alpha$, while the
latter can be achieved if $\widetilde{\mathcal{P}}^{A|W}$ is independent of
$W$ (e.g., $A_{\textup{test}}\equiv 0$). When each cluster only has one
individual, Algorithm 1 and Algorithm 2 coincide, and their resulting coverage
guarantees also become identical.

From the perspective of implementation, conformal causal inference for
individual-level treatment effects is advantageous in model training over that
for cluster-level treatment effects. As a concrete context, given $10$
clusters, learning the cluster-level model for
$\overline{Y}(1)-\overline{Y}(0)$ includes $10$ data points as its effective
sample size, which may potentially lead to overfitting with a handful of
covariates. In contrast, the individual-level model for $Y(1)-Y(0)$ often has
a larger effective sample size for training, e.g., 300 data points given 30
individuals per cluster. Since the within-cluster correlation is usually small
(Campbell et al.,, 2005), the increased effective sample size in the training
data can substantially improve the numerical stability. As a result, conformal
causal inference for individual-level treatment effects could be more
informative in CRTs with a small $m$ compared to that for cluster-level
treatment effects. As we demonstrate in the ensuing simulations, with $m=30$
clusters, we are able to perform conformal subgroup analysis for the
individual-level treatment effect but not the cluster-level treatment effect.

In terms of the length of conformal intervals,
$\widetilde{\textup{C}}_{C}(\overline{B}_{\textup{test}})$ tends to be more
informative than $\widetilde{\textup{C}}_{I}(B_{\textup{test}})$ given a
sufficient number of clusters. This is because
$\overline{Y}(1)-\overline{Y}(0)$ often has smaller variance than $Y(1)-Y(0)$,
especially if $N_{i}$ is large. As a result,
$\widetilde{\textup{C}}_{C}(\overline{B}_{\textup{test}})$ is more likely to
exclude zero than $\widetilde{\textup{C}}_{I}(B_{\textup{test}})$ given the
same treatment effect size. In practice, choosing between the two types of
inferential targets requires a case-by-case evaluation. Although the
scientific question should drive the target of inference, from a statistical
perspective, conformal causal inference for the cluster-level treatment
effects with $\alpha=0.1$ is typically more informative and efficient when $m$
is large, e.g., $m\geq 80$. Given a small to moderate number of clusters, the
individual-level treatment effects are preferable, and $\alpha$ can be
adjusted accordingly to meet the practical needs.

###  4.2 Inference based on individual-level covariates

When the test point only contains individual-level covariates
$B_{\textup{test}}$, we follow the same strategy as in Section 3.2 to
construct conformal intervals. Corollary 2 characterizes the local coverage
property of the direct approach as an application of Theorem 2.

######  Corollary 2.

Under Assumptions 1-3 and assuming that
$(Y_{\textup{test}}(1),Y_{\textup{test}}(0),B_{\textup{test}})$ is an
arbitrary individual from a new cluster independently sampled from
$\mathcal{P}^{W}$. Then $\widetilde{\textup{C}}_{I,1}(B_{\textup{test}})$ and
$\widetilde{\textup{C}}_{I,0}(B_{\textup{test}})$ output by Algorithm 2
satisfy

$P\bigg{\\{}Y_{\textup{test}}(1)-Y_{\textup{test}}(0)\in\widetilde{\textup{C}}_{I,1}(B_{\textup{test}})-\widetilde{\textup{C}}_{I,0}(B_{\textup{test}})\bigg{|}B_{\textup{test}}\in\Omega_{I}\bigg{\\}}\geq
1-\alpha$ (8)

for any set $\Omega_{I}$ in the support of $B_{\textup{test}}$ with a positive
measure.

In parallel to Corollary 1, Corollary 2 establishes the coverage guarantee on
conformal intervals for the individual-level treatment effect based on
covariates. In addition, the resulting conformal interval enjoys the same
benefit of stability and small-sample compatibility as discussed in Section
4.1. In addition to the direct approach, we provided the nested approach with
$1-\alpha-\gamma$ coverage guarantee in the Supporting Information.

##  5 Simulations

Through a simulation study, we aim to demonstrate our algorithms and verify
our finite-sample theoretical results for both cluster-level and individual-
level treatment effects. We consider the combination of the following
settings: CRTs with a large ($m=100$) versus small ($m=30$) number of
clusters, and full-data analysis versus covariate subgroup analysis.

###  5.1 Simulation design

For $i=1,\dots,m$, we independently generate the cluster size
$N_{i}\sim\mathcal{U}([10,50])$ and two cluster-level covariates
$R_{i1}|N_{i}\sim\mathcal{N}(N_{i}/10,1)$,
$R_{i2}|(N_{i},R_{i1})\sim\mathcal{B}\\{(1+e^{-R_{i1}/2})^{-1}\\}$, where
$\mathcal{U},\mathcal{N},\mathcal{B}$ represent the uniform, normal, and
Bernoulli distribution, respectively. For each individual $j=1,\dots,N_{i}$,
we independently generate covariates
$X_{ij1}|(N_{i},R_{i1},R_{i2})\sim\mathcal{B}(0.3+0.4R_{i2})$ and
$X_{ij2}=(2I\\{R_{i1}>0\\}-1)\overline{X}_{i1}+\mathcal{N}(0,1)$, and
potential outcomes

$\displaystyle Y_{ij}(0)$
$\displaystyle=\sin(R_{i1})(2R_{i2}-1)+|X_{ij1}X_{ij2}|+\gamma_{i}+\varepsilon_{ij}$
$\displaystyle Y_{ij}(1)$
$\displaystyle=\sin(R_{i1})(2R_{i2}-1)+|X_{ij1}X_{ij2}|+N_{i}/50+\varepsilon_{ij},$

where $\gamma_{i}\sim\mathcal{N}(0,0.5^{2})$ is the random intercept and
$\varepsilon_{ij}\sim N(0,1)$ is the random noise. Then we independently
generate the treatment indicator $A_{i}\sim\mathcal{B}(0.5)$, and obtain
$Y_{ij}=A_{i}Y_{ij}(1)+(1-A_{i})Y_{ij}(0)$. The simulated observed data are
$\\{(Y_{ij},A_{i},X_{ij1},X_{ij2},R_{i1},R_{i2}):i=1,\dots,N_{i}\\}_{i=1}^{m}$.
In addition, we generate 1,000 new clusters as the test data, following the
same data-generating distribution. We repeat the above procedure to generate
1,000 data replicates for evaluation of the conformal interval estimators.

For each simulated data set, we first construct the conformal interval for the
cluster-average treatment effect. Given complete test data (with treatment and
outcomes), we run Algorithm 1, and refer to this approach as “O”. Given
cluster-level covariates only, we refer to the direct approach as “B-direct”
and the nested approach as “B-nested”. For covariate subgroup analysis, we
consider $\Omega_{C}=\\{\overline{B}_{i}:R_{i1}\geq 2,R_{i2}=1\\}$, which
contains 60% of all clusters. For inferring the individual-average treatment
effect, we adopt the same names “O”, ”B-direct”, “B-nested” to refer to the
output of Algorithm 2, the direct approach and the nested approach, and the
covariate subgroup is defined as $\Omega_{I}=\\{B_{ij}:|X_{ij2}|<0.5\\}$,
which includes $30\%$ of all individuals. For each method, the training model
uses an ensemble learner of linear regression and random forest through the
SuperLearner package (van der Laan et al.,, 2007) in the R software. We
consider two metrics of performance based on the test data: the probability
that the conformal interval contains the true treatment effect, and the
average length of the conformal interval. When $m=100$, we set $\alpha=0.1$,
while we consider both $\alpha=0.1$ and $\alpha=0.2$ given $m=30$. For the
nested approach, we set $\gamma=0.5$ to improve the informativeness of the
resulting conformal interval.

###  5.2 Simulation results

Figure 1 summarizes the simulation results for the marginal and local cluster-
level treatment effects given $m=100$. In Figure 1, the upper panels show that
all three methods achieve the target 90% coverage probability, thereby
confirming our theoretical results. Comparing the three methods, the “O”
method has a coverage probability close to 0.9, whereas the “B-direct” method
is the most conservative one (coverage probability near 1). This difference
can be explained by the lower panels, where the “B-direct” method yields wider
conformal intervals, whose length nearly doubles the oracle length. In
contrast, due to leveraging the complete information in the test data, the “O”
method achieves the near-optimal length of the conformal interval. The
performance of the “B-nested” method lies between the other two since we set a
loose parameter $\gamma=0.5$. If we use $\gamma=0.1$ instead, it will perform
similarly to the “B-direct” method as demonstrated in Lei and Candès, (2021)
for non-clustered data. Comparing marginal versus local treatment effects, the
coverage probability and length of conformal intervals are both similar, while
the results for the local treatment effect are slightly less stable with fewer
data points. In the Supporting Information, we reproduce Figure 1 with
$m=500$, where the span of the boxplot substantially decreases and the “O”
method is nearly optimal under both metrics, confirming the large-sample
stability.

Figure 1: Simulation results (boxplot) for the marginal (left column) and
local (right column, conditioning on $\\{R_{i1}\geq 2,R_{i2}=1\\}$) cluster-
level treatment effects with $m=100$. In the upper panels, the dashed line is
the target 90% coverage probability. In the lower panels, the solid line is
the oracle length of conformal intervals, computed as the average length
between the $(\alpha/2,1-\alpha/2)$-quantiles of
$\overline{Y}(1)-\overline{Y}(0)$ among test data.

Figure 2 presents the simulation results for the marginal and local
individual-level treatment effects given $m=100$. All three methods reach the
target coverage probability, and their comparisons are similar to Figure 1.
However, the conformal intervals for the individual-level treatment effect are
more stable than the cluster-level treatment effect, as reflected by shorter
spans of the boxplot, but have wider lengths due to the increased variance in
treatment effects. The large sample stability is similarly confirmed by
reproducing Figure 2 with $m=500$ in the Supporting Information.

Figure 2: Simulation results (boxplot) for the marginal (left column) and
local (right column, conditioning on $\\{|X_{ij2}|<0.5\\}$) individual-level
treatment effects with $m=100$. In the upper panels, the dashed line is the
target 90% coverage probability. In the lower panels, the solid line is the
oracle length of conformal intervals, computed as the average length between
the $(\alpha/2,1-\alpha/2)$-quantiles of $Y(1)-Y(0)$ among test data.

When $m=30$, we provide the simulation results in Table 1. Due to limited
clusters, we do not implement the “B-nested” method (which requires two sample
splittings) or perform the covariate-subgroup analysis for the cluster-level
treatment effect (which drops 40% of all clusters). For the remaining
scenarios, the “O” and “B-direct” methods can still achieve the target
coverage probability, confirming their finite-sample validity. However, they
become more conservative compared to $m=100$, as reflected by the increased
length of intervals. This is expected given limited clusters used in model
training and quantile estimation. An ad-hoc solution is to consider an 80%
coverage probability, in which case the resulting conformal intervals have a
reduced length and are likely more informative. Finally, it is worth
highlighting that, even with few clusters, we are able to perform subgroup
analysis for the individual-level treatment effect, and its empirical
performance is similar to the marginal analysis that uses all data. This
observation confirms the statistical benefit of conformal causal inference for
individual-level treatment effects in small CRTs.

Table 1: Summary of simulation results for 80% and 90% conformal intervals
given $m=30$ clusters. For both coverage probability and length of intervals,
we report their averages and standard errors. The local treatment effect is
conditioned on $|X_{ij1}|<0.5$.

Treatment effects | Methods | $\alpha=0.2$ | $\alpha=0.1$
---|---|---|---
Coverage probability |  Length of intervals |  Coverage probability |  Length of intervals
Marginal cluster-level | O | 0.816(0.094) | 3.119(1.465) | 0.912(0.072) | 4.211(2.230)
B-direct | 0.998(0.008) | 5.891(2.833) | 1.000(0.002) | 7.970(4.447)
Marginal individual-level | O | 0.868(0.035) | 4.311(0.489) | 0.977(0.014) | 6.637(0.794)
B-direct | 0.989(0.007) | 8.345(0.768) | 1.000(0.000) | 12.905(1.232)
Local individual-level | O | 0.870(0.048) | 4.337(0.675) | 0.977(0.021) | 6.933(1.260)
B-direct | 0.988(0.009) | 8.398(1.138) | 1.000(0.000) | 13.542(2.114)

In the Supporting Information, we provide additional simulations to evaluate
the benefit of machine learning in fitting the prediction model. Specifically,
we repeat our simulation study for $m=30$ with only linear regression as
$f_{a}$ (in contrast to the ensemble method with linear regression and random
forest). We find that, by including random forest, the coverage probability
has negligible differences, but the length of intervals is reduced by
8.2-46.5% across scenarios. This example further demonstrates that our method
is able to leverage the power of machine learning to improve accuracy.

##  6 Data example with the PPACT cluster randomized trial

The Pain Program for Active Coping and Training study (PPACT) is a CRT
evaluating the effect of a care-based cognitive behavioral therapy
intervention for treating long-term opioid users with chronic pain (DeBar et
al.,, 2022). The study equally randomized 106 primary care providers
(clusters) to receive the intervention or usual care, with 1-10 participants
in each cluster. We focus on the primary outcome, the PEGS (pain intensity and
interference with enjoyment of life, general activity, and sleep) score at 12
months, a continuous measure of pain scale ranging from 1 to 10. For more
accurate conformal causal inference, we adjust for 13 individual-level
baseline variables, including the baseline PEGS score, age, gender,
disability, smoking status, body mass index, alcohol abuse, drug abuse,
comorbidity, depression, number of pain types, average morphine dose, and
heavy opioid usage.

Since the underlying truth of PPACT is unknown, we randomly sample 20 clusters
as the test data and use the remaining 86 clusters to construct the conformal
intervals. To account for the uncertainty introduced by random sampling, we
repeat this procedure 100 times and report the average and standard error for
two performance metrics: length of intervals and fraction of negatives. Here,
the fraction of negatives defines the proportion of conformal intervals that
are subsets of $(-\infty,0)$ among the test data. Since negative values are in
the direction of treatment benefits, this metric reveals how many
clusters/individuals are associated with beneficial treatment effects with
probability $1-\alpha$. Because the test data have the complete information,
we directly run the “O” method, i.e., Algorithm 1 and Algorithm 2, with
$f_{a}$ set as the ensemble method of linear regression and random forest.

Table 2 summarizes the results for the marginal treatment effects setting
$\alpha\in\\{0.1,0.2,0.3,0.4\\}$. As $\alpha$ increases, the length of
intervals decreases, and the fraction of negative becomes larger. Since the
treatment effect is small (relative to the variability of treatment effects,
DeBar et al.,, 2022), only a small to moderate proportion of the population
has negative conformal intervals. In practice, these negative conformal
intervals can be informative for new patients with similar baseline
characteristics.

Table 2: Summary results of data application for marginal treatment effects.
For both the length of intervals and the fraction of negatives, we present the
average and standard error over 100 runs.

Coverage probability | Marginal cluster-level treatment effect | Marginal individual-level treatment effect
---|---|---
Length of intervals | Fraction of negatives | Length of intervals | Fraction of negatives
90% | 4.056(0.557) | 0.089(0.069) | 6.908(0.590) | 0.055(0.026)
80% | 2.874(0.412) | 0.173(0.092) | 4.800(0.345) | 0.132(0.044)
70% | 2.233(0.313) | 0.238(0.104) | 3.811(0.220) | 0.199(0.052)
60% | 1.799(0.255) | 0.304(0.117) | 3.108(0.189) | 0.257(0.055)

In addition to the marginal analysis, we performed subgroup analyses on
individuals with severe baseline pain (baseline PEGS score equal to or larger
than 7, 40% participants) and moderate baseline pain (baseline PEGS score
between 4 and 7, 50% participants). We summarize the results in Table 3, which
shows no difference in the length of intervals or the fraction of negatives
between the two subgroups. Overall, these analyses indicate treatment benefits
for a small subset of the clusters/individuals, which drive the overall
treatment effect signal provided in previous analyses (DeBar et al.,, 2022;
Wang et al.,, 2023).

Table 3: Summary results of data application for individual-level treatment
effects on covariate subgroups. For both the length of intervals and the
fraction of negatives, we present the average and standard error over 100
runs.

Coverage probability | Subgroup: severe baseline pain | Subgroup: moderate baseline pain
---|---|---
Length of intervals | Fraction of negatives | Length of intervals | Fraction of negatives
90% | 6.538(0.571) | 0.072(0.042) | 6.508(0.538) | 0.076(0.042)
80% | 4.982(0.465) | 0.120(0.051) | 4.928(0.417) | 0.132(0.053)
70% | 3.935(0.353) | 0.180(0.064) | 3.960(0.349) | 0.185(0.068)
60% | 3.180(0.262) | 0.243(0.067) | 3.221(0.296) | 0.241(0.079)

##  7 Discussion

In this work, we contributed a novel conformal causal inference framework for
finite-sample model-robust inference in CRTs. Distinct from the classical
framework that focuses on average treatment effects and often relies on the
asymptotic theory, conformal causal inference yields covariate-based (thus
personalized) intervals and is robust to arbitrary working model
misspecification without appealing to asymptotic regimes. These new features
address the limitation of classical methods in analyzing CRTs and provide a
fresh perspective on treatment effects with a new inferential target. More
broadly, we expand the conformal inference toolkit to accommodate clustered
data and their covariate subgroups, addressing model-robust finite-sample
validity under a broader set of data-generating distributions.

Although our data example in Section 6 includes a large number of clusters,
our proposed methods are not confined to a large number of clusters, given
appropriate adjustments of the target coverage probability. As discussed in
Sections 3, for achieving 90% coverage probability of the conformal interval
for the cluster-level treatment effects, we need at least 10 clusters per arm
for numerically stable calibration, and our simulation shows that the
conformal interval is valid but moderately conservative given $m=30$ and
$\pi=0.5$. Given fewer clusters, e.g., $m=20$, targeting 90% coverage
probability will only result in intervals spanning the entire real line, and
informative intervals are possible with a lower target coverage probability
(say 80%). In practice, we recommend setting $\alpha=0.05$ if $m\geq 80$,
$\alpha=0.1$ if $40\leq m<80$, and $\alpha=0.2$ if $20\leq m<40$ given
$\pi=0.5$. If $10<m<20$, we can still construct the 80% conformal interval for
the individual-level treatment effect, while it may be challenging to make
conformal causal inference given $m\leq 10$. Conformal causal inference with
an extremely small number of clusters is a topic of future research.

## References

  * Alaa et al., (2023) Alaa, A., Ahmad, Z., and van der Laan, M. (2023).  Conformal meta-learners for predictive inference of individual treatment effects.  arXiv preprint arXiv:2308.14895.
  * Balzer et al., (2016) Balzer, L. B., Petersen, M. L., van der Laan, M. J., and Collaboration, S. (2016).  Targeted estimation and inference for the sample average treatment effect in trials with and without pair-matching.  Statistics in Medicine, 35(21):3717–3732.
  * Balzer et al., (2023) Balzer, L. B., van der Laan, M., Ayieko, J., Kamya, M., Chamie, G., Schwab, J., Havlir, D. V., and Petersen, M. L. (2023).  Two-stage TMLE to reduce bias and improve efficiency in cluster randomized trials.  Biostatistics, 24(2):502–517.
  * Barber et al., (2021) Barber, R. F., Candès, E. J., Ramdas, A., and Tibshirani, R. J. (2021).  Predictive inference with the jackknife+.  The Annals of Statistics, 49(1):486 – 507.
  * Benitez et al., (2023) Benitez, A., Petersen, M. L., van der Laan, M. J., Santos, N., Butrick, E., Walker, D., Ghosh, R., Otieno, P., Waiswa, P., and Balzer, L. B. (2023).  Defining and estimating effects in cluster randomized trials: A methods comparison.  Statistics in Medicine, 42(19):3443–3466.
  * Breiman, (2001) Breiman, L. (2001).  Random forests.  Machine learning, 45:5–32.
  * Campbell et al., (2005) Campbell, M. K., Fayers, P. M., and Grimshaw, J. M. (2005).  Determinants of the intracluster correlation coefficient in cluster randomized trials: the case of implementation research.  Clinical Trials, 2(2):99–107.
  * Crespi et al., (2009) Crespi, C. M., Wong, W. K., and Mishra, S. I. (2009).  Using second-order generalized estimating equations to model heterogeneous intraclass correlation in cluster-randomized trials.  Statistics in Medicine, 28(5):814–827.
  * DeBar et al., (2022) DeBar, L., Mayhew, M., Benes, L., Bonifay, A., Deyo, R. A., Elder, C. R., Keefe, F. J., Leo, M. C., McMullen, C., Owen-Smith, A., et al. (2022).  A primary care–based cognitive behavioral therapy intervention for long-term opioid users with chronic pain: a randomized pragmatic trial.  Annals of Internal Medicine, 175(1):46–55.
  * Ding and Keele, (2018) Ding, P. and Keele, L. (2018).  Rank tests in unmatched clustered randomized trials applied to a study of teacher training.  The Annals of Applied Statistics, 12(4):2151–2174.
  * Hore and Barber, (2023) Hore, R. and Barber, R. F. (2023).  Conformal prediction with local weights: randomization enables local guarantees.  arXiv preprint arXiv:2310.07850.
  * Ivers et al., (2011) Ivers, N., Taljaard, M., Dixon, S., Bennett, C., McRae, A., Taleban, J., Skea, Z., Brehaut, J., Boruch, R., Eccles, M., et al. (2011).  Impact of consort extension for cluster randomised trials on quality of reporting and study methodology: review of random sample of 300 trials, 2000-8.  Bmj, 343.
  * Jin et al., (2023) Jin, Y., Ren, Z., and Candès, E. J. (2023).  Sensitivity analysis of individual treatment effects: A robust conformal inference approach.  Proceedings of the National Academy of Sciences, 120(6):e2214889120.
  * Kahan et al., (2023) Kahan, B. C., Li, F., Copas, A. J., and Harhay, M. O. (2023).  Estimands in cluster-randomized trials: choosing analyses that answer the right question.  International Journal of Epidemiology, 52(1):107–118.
  * Lee et al., (2023) Lee, Y., Barber, R. F., and Willett, R. (2023).  Distribution-free inference with hierarchical data.  arXiv preprint arXiv:2306.06342.
  * Lei et al., (2018) Lei, J., G’Sell, M., Rinaldo, A., Tibshirani, R. J., and Wasserman, L. (2018).  Distribution-free predictive inference for regression.  Journal of the American Statistical Association, 113(523):1094–1111.
  * Lei and Wasserman, (2014) Lei, J. and Wasserman, L. (2014).  Distribution-free prediction bands for non-parametric regression.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 76(1):71–96.
  * Lei and Candès, (2021) Lei, L. and Candès, E. J. (2021).  Conformal inference of counterfactuals and individual treatment effects.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 83(5):911–938.
  * Murray et al., (1998) Murray, D. M. et al. (1998).  Design and Analysis of Group-Randomized Trials, volume 29.  Oxford University Press, USA.
  * Papadopoulos et al., (2002) Papadopoulos, H., Proedrou, K., Vovk, V., and Gammerman, A. (2002).  Inductive confidence machines for regression.  In Machine Learning: ECML 2002: 13th European Conference on Machine Learning Helsinki, Finland, August 19–23, 2002 Proceedings 13, pages 345–356. Springer.
  * Qiu et al., (2022) Qiu, H., Dobriban, E., and Tchetgen, E. T. (2022).  Prediction sets adaptive to unknown covariate shift.  arXiv preprint arXiv:2203.06126.
  * Rabideau and Wang, (2021) Rabideau, D. J. and Wang, R. (2021).  Randomization-based confidence intervals for cluster randomized trials.  Biostatistics, 22(4):913–927.
  * Small et al., (2008) Small, D. S., Ten Have, T. R., and Rosenbaum, P. R. (2008).  Randomization inference in a group–randomized trial of treatments for depression: covariate adjustment, noncompliance, and quantile effects.  Journal of the American Statistical Association, 103(481):271–279.
  * Su and Ding, (2021) Su, F. and Ding, P. (2021).  Model-assisted analyses of cluster-randomized experiments.  Journal of the Royal Statistical Society, Series B, 83(5):994–1015.
  * Tibshirani et al., (2019) Tibshirani, R. J., Foygel Barber, R., Candes, E., and Ramdas, A. (2019).  Conformal prediction under covariate shift.  Advances in neural information processing systems, 32:1–11.
  * van der Laan et al., (2007) van der Laan, M. J., Polley, E. C., and Hubbard, A. E. (2007).  Super learner.  Statistical Applications in Genetics and Molecular Biology, 6(1):1–21.
  * Vovk et al., (2005) Vovk, V., Gammerman, A., and Shafer, G. (2005).  Algorithmic learning in a random world, volume 29.  Springer.
  * Wang et al., (2021) Wang, B., Harhay, M. O., Small, D. S., Morris, T. P., and Li, F. (2021).  On the mixed-model analysis of covariance in cluster-randomized trials.  arXiv preprint arXiv:2112.00832.
  * Wang et al., (2023) Wang, B., Park, C., Small, D. S., and Li, F. (2023).  Model-robust and efficient inference for cluster-randomized experiments.  Journal of the American Statistical Association.
  * Yang et al., (2022) Yang, Y., Kuchibhotla, A. K., and Tchetgen, E. T. (2022).  Doubly robust calibration of prediction sets under covariate shift.  arXiv preprint arXiv:2203.01761.
  * Yin et al., (2022) Yin, M., Shi, C., Wang, Y., and Blei, D. M. (2022).  Conformal sensitivity analysis for individual treatment effects.  Journal of the American Statistical Association, pages 1–14.
  * Zhang et al., (2012) Zhang, K., Traskin, M., and Small, D. S. (2012).  A powerful and robust test statistic for randomization inference in group-randomized trials with matched pairs of groups.  Biometrics, 68(1):75–84.