# The role of randomization inference in unraveling individual treatment
effects in clinical trials: Application to HIV vaccine trials

Zhe Chen  Department of Statistics, University of Illinois, Urbana-Champaign
Xinran Li  Department of Statistics, University of Chicago  Bo Zhang Assistant
Professor of Biostatistics, Vaccine and Infectious Disease Division, Fred
Hutchinson Cancer Center. Email<EMAIL_ADDRESS>Vaccine and Infectious
Disease Division, Fred Hutchinson Cancer Center

Abstract: Randomization inference is a powerful tool in early phase vaccine
trials to estimate the causal effect of a regimen against a placebo or another
regimen. Traditionally, randomization-based inference often focuses on testing
either Fisher’s sharp null hypothesis of no treatment effect for any unit or
Neyman’s weak null hypothesis of no sample average treatment effect. Many
recent efforts have explored conducting exact randomization-based inference
for other summaries of the treatment effect profile, for instance, quantiles
of the treatment effect distribution function. In this article, we
systematically review methods that conduct exact, randomization-based
inference for quantiles of individual treatment effects (ITEs) and extend some
results by incorporating auxiliary information often available in a vaccine
trial. These methods are suitable for four scenarios: (i) a randomized
controlled trial (RCT) where the potential outcomes under one regimen are
constant; (ii) an RCT with no restriction on any potential outcomes; (iii) an
RCT with some user-specified bounds on potential outcomes; and (iv) a matched
study comparing two non-randomized, possibly confounded treatment arms. We
then conduct two extensive simulation studies, one comparing the performance
of each method in many practical clinical settings and the other evaluating
the usefulness of the methods in ranking and advancing experimental therapies.
We apply these methods to an early-phase clinical trail, HIV Vaccine Trials
Network Study 086 (HVTN 086), to showcase the usefulness of the methods.

Keywords: Clinical trials, Individual treatment effect, Effect quantile,
Randomization inference, Vaccine

##  1 Introduction

###  1.1 Early phase vaccine trials; vaccine-induced immune responses;
heterogeneity

One central goal in vaccine development and a primary objective in many early
phase clinical trials of experimental medicine is to evaluate the vaccine-
induced immunogenicity. Vaccine-induced immune responses are often
heterogeneous among study participants. To illustrate, Figure 1 exhibits the
observed serum IgG binding antibody multiplex assay (BAMA) responses to two
antigens, Con 6 gp120/B and gp41, among study participants in a phase 1,
multi-arm, placebo-controlled clinical trial conducted via the HIV Vaccine
Trials Network (HVTN). Ditse et al., (2020); Huang et al., (2022) The study
(HVTN 086/SAAVI 103) enrolled a total of $184$ participants into $4$ study
arms; within each study arm, participants were randomized to receive a
candidate vaccine regimen or placebo. It is transparent from Figure 1 that
vaccine recipients’ binding antibody responses ranged from “potent and mostly
homogeneous” (responses to gp41 among recipients of regimens 1 and 4) to “weak
and heterogeneous” (responses to Con 6 gp120/B among recipients of regimen 2).
This within-participant heterogeneity in vaccine-induced immune responses has
been well-documented in many vaccines, including those for Covid-19,
influenza, dengue, and hepatitis B, Huang et al., (2022) and could at least
partially explain the lack of efficacy in phase 3, HIV-1 vaccine trials.

Figure 1: Observed serum IgG BAMA responses to Con6 gp120/B (left panel) and
gp41 (right panel) among study participants in the HIV Vaccine Trials Network
(HVTN) study 086. A small perturbation was added to each observation to aid
data visualization.

Researchers routinely characterize a vaccine regimen’s induced immune
responses (e.g., the BAMA responses in Figure 1) by estimating and reporting
its sample average treatment effect (SATE) against the placebo, which could be
unbiasedly estimated in a randomized clinical trial, and assess and rank
multiple regimens by comparing their estimated SATEs. When developing a
challenging vaccine product like an HIV-1 vaccine, researchers have long
realized that the response rates among study participants are often highly
variable, and a significant proportion of participants could exhibit immune
responses below the assay limit of detection (LOD) or lower limit of
quantification (LLOQ). Hence, summarizing and comparing immune response
profiles based on the mean difference alone could mask significant and perhaps
meaningful heterogeneities. The current standard practice is to complement the
estimated SATE and its 95% confidence interval by further reporting (i)
descriptive statistics and boxplots summarizing the spread of immune responses
elicited by each regimen, (ii) the percentage of positive or high responders
within each group, and (iii) the mean response among the subset of positive or
high responders. Unfortunately, unlike the SATE which is a well-defined,
albeit less than comprehensive, causal estimand, neither the descriptive
statistics nor the mean difference in responses among positive or high
responders constitutes a formal confidence statement about the “treatment
effect” of a regimen versus placebo or another regimen.

###  1.2 Science table; estimands of interest

What is a well-defined causal estimand that captures the treatment effect
heterogeneities in early phase clinical trials? It is helpful to examine what
Rubin Rubin, (2005) refers to as a “science table.” Table 1 summarizes the
potential outcomes and unit-level treatment effects of $N_{1}+N_{0}$ study
participants in a clinical trial, where $N_{1}$ study participants are
randomized to regimen $1$ and the other $N_{0}$ participants regimen $0$. In
the science table 1, $Y_{i}(1)$ and $Y_{i}(0)$ denote study participant $i$’s
potential immune responses of interest under two regimens, and only one of
these two responses is observed depending on the actual regimen assigned to
the participant. The contrast in the two potential outcomes,
$\tau_{i}=Y_{i}(1)-Y_{i}(0)$, denotes the unit-level treatment effect. Rubin,
(2005) We will refer to $\tau_{i}$ as an _individual treatment effect_ (ITE)
in this article,Lipkovich et al., (2023) and the collection of ITEs,
$\mathcal{T}=\\{\tau_{i},~{}i=1,\dots,N_{1}+N_{0}\\}$, is the causal quantity
of _ultimate interest_ , in the sense that any summary treatment effect can be
derived from $\mathcal{T}$ (e.g., the sample average treatment effect
$\overline{\tau}$). Let $\tau_{(i)}$ denote the $i$th largest treatment
effect. The immunogenicity profile of a vaccine regimen (against placebo or a
competing regimen), as revealed in an early phase clinical trial, is
completely characterized by
$\mathcal{T}=\\{\tau_{(i)},~{}i=1,\dots,N_{1}+N_{0}\\}$. Unfortunately,
elements in $\mathcal{T}$ are almost never completely observed, so statistical
inference is needed.

Table 1: Science table of $N_{0}+N_{1}$ study participants. A total of $N_{0}$
are randomized to regimen 0 and the other $N_{1}$ are randomized to regimen
$1$. Each participant is associated with two potential immune responses
$Y_{i}(0)$ and $Y_{i}(1)$ corresponding to regimen $0$ and $1$, though only
the one corresponding to the actual regimen assignment is observed (boldface).
Each participant $i$ is associated with an individual treatment effect
$\tau_{i}$.

Participants |  Regimen Assignment | Potential Outcomes |  Individual Treatment Effects
---|---|---|---
$Y(1)$ | $Y(0)$
$1$ | Regimen 1 | $\mathbf{Y_{1}(1)}$ | $Y_{1}(0)$ | $\tau_{1}=Y_{1}(1)-Y_{1}(0)$
$2$ | Regimen 1 | $\mathbf{Y_{2}(1)}$ | $Y_{2}(0)$ | $\tau_{2}=Y_{2}(1)-Y_{2}(0)$
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
$N_{1}$ | Regimen 1 | $\mathbf{Y_{N_{1}}(1)}$ | $Y_{N_{1}}(0)$ | $\tau_{N_{1}}=Y_{N_{1}}(1)-Y_{N_{1}}(0)$
$N_{1}+1$ | Regimen 0 | $Y_{N_{1}+1}(1)$ | $\mathbf{Y_{N_{1}+1}(0)}$ | $\tau_{N_{1}+1}=Y_{N_{1}+1}(1)-Y_{N_{1}+1}(0)$
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
$N_{1}+N_{0}$ | Regimen 0 | $Y_{N_{1}+N_{0}}(1)$ | $\mathbf{Y_{N_{1}+N_{0}}(0)}$ | $\tau_{N_{1}+N_{0}}=Y_{N_{1}+N_{0}}(1)-Y_{N_{1}+N_{0}}(0)$

###  1.3 Outline of the article

We have four goals in this article. First, we provide a brief overview of
existing methodologies that examine weak null hypotheses regarding the
distribution function of ITEs. Second, we extend some results by incorporating
useful auxiliary information into existing methods and derive both individual
and simultaneous inference. Third, we compare different methods and evaluate
their performance in (i) characterizing ITE profiles and (ii) ranking regimens
under a wide range of practical data-generating processes mimicking early
phase clinical trials. Finally, we present a case study of immunogenicity data
derived from an early phase HIV vaccine trial and explore how different
methods can facilitate optimizing decision-making and improving the evaluation
of vaccine regimens.

##  2 Framework, notation, and different null hypotheses

###  2.1 Randomization inference; Fisher’s sharp hypothesis; Neyman’s weak
null hypothesis

We consider a two-arm randomized experiment on $N$ study participants under
Neyman-Rubin’s potential outcomes framework. Neyman, (1923); Rubin, (1974) Let
$Y_{i}(1)$ and $Y_{i}(0)$ denote potential outcomes of unit $i$ under two
competing regimens and $\tau_{i}=Y_{i}(1)-Y_{i}(0)$ unit $i$’s individual
treatment effect comparing Regimen 1 to Regimen 0. We collect the set of $N$
potential outcomes under Regimen 1 in
$\mathbf{Y}(1)=(Y_{1}(1),Y_{i}(1),\dots,Y_{N}(1))^{\intercal}$, those under
Regimen 0 in $\mathbf{Y}(0)=(Y_{1}(0),Y_{i}(0),\dots,Y_{N}(0))^{\intercal}$,
and the set of all ITEs in
$\boldsymbol{\tau}=(\tau_{1},\tau_{2},\dots,\tau_{N})^{\intercal}$. Let
$Z_{i}$ denote the treatment assignment to study participant $i$ such that
$Z_{i}=1$ if participant $i$ is assigned Regimen 1 and $0$ otherwise. The set
of treatment assignments is collected in
$\mathbf{Z}=(Z_{1},Z_{2},\dots,Z_{N})^{\intercal}$. For each study participant
$i$, the observed outcome $Y_{i}$ satisfies $Y_{i}=Z_{i}\cdot
Y_{i}(1)+(1-Z_{i})\cdot Y_{i}(0)$. We collect $N$ observed outcomes in the
vector $\mathbf{Y}=(Y_{1},Y_{2},\dots,Y_{N})^{\intercal}$. In randomization
inference, potential outcomes of study participants are viewed as fixed
quantities and researchers rely _solely_ on the treatment assignment mechanism
to draw valid causal conclusions. In other words, randomization forms what Sir
Ronald Fisher referred to as the “reasoned basis” for causal inference.
Fisher, (1935) In the current paper, we will focus primarily on the completely
randomized experiment (CRE), where $N_{1}$ study participants are randomly
assigned to Regimen 1 and the other $N_{0}=N-N_{1}$ to Regimen 0. In his
seminal work _Design of Experiments_ , FisherFisher, (1935) advocated testing
the following _sharp_ null hypothesis:

$H_{0,\text{sharp}}\mathrel{\mathop{\ordinarycolon}}\mathbf{Y}(1)=\mathbf{Y}(0),$
(1)

which states that the treatment has no effect whatsoever on any study
participant, or equivalently $\boldsymbol{\tau}=0$. Under this null
hypothesis, the missing potential outcome of each study participant can then
be imputed because $Y_{i}(1)=Y_{i}(0)=Y_{i}$ for $i=1,\dots,N$, and any test
statistic has its null distribution completely specified, leading to an
_exact_ $p$-value upon comparing the observed value of the chosen test
statistic to this null distribution. The Fisher randomization test (FRT)
scheme works for a generic sharp null hypothesis of the following form:

$H_{\boldsymbol{\delta}}\mathrel{\mathop{\ordinarycolon}}\boldsymbol{\tau}=\boldsymbol{\delta},$
(2)

where
$\boldsymbol{\delta}=\left(\delta_{1},\delta_{2},\ldots,\delta_{N}\right)^{\intercal}\in\mathbb{R}^{N}$
is a prespecified vector of constants. Let $t(\boldsymbol{z},\boldsymbol{y})$
denote a generic function of the assignment vector
$\boldsymbol{z}\in\\{0,1\\}^{N}$ and outcome vector
$\boldsymbol{y}\in\mathbb{R}^{N}$. Following Rosenbaum,Rosenbaum, (2002) we
consider test statistics of the form
$t\left(\boldsymbol{Z},\boldsymbol{Y}_{\boldsymbol{Z},\boldsymbol{\delta}}(0)\right)$,
where $\boldsymbol{Z}$ is the vector of observed treatment assignment and
$\boldsymbol{Y}_{\boldsymbol{Z},\boldsymbol{\delta}}(0)$ is the vector of
potential outcomes under Regimen 0. Given the observed data $\boldsymbol{Y}$
and under $H_{\boldsymbol{\delta}}$, we have
$\boldsymbol{Y}_{\boldsymbol{Z},\boldsymbol{\delta}}(0)=\boldsymbol{Y}-\boldsymbol{Z}\circ\boldsymbol{\delta}$,
where $\circ$ represents element-wise multiplication. In a CRE design, the
null distribution of the test statistic has the following tail probability:

$\displaystyle G_{\boldsymbol{Z},\boldsymbol{\delta}}(c)$
$\displaystyle\equiv\Pr\\{t(\boldsymbol{Z},\boldsymbol{Y}_{\boldsymbol{Z},\boldsymbol{\delta}}(0))\geq
c\\}=\binom{N}{N_{1}}^{-1}\cdot\sum_{\boldsymbol{z}\in\\{0,1\\}^{N}\mathrel{\mathop{\ordinarycolon}}\sum_{i=1}^{N}z_{i}=N_{1}}\mathbbm{1}\left\\{t(\boldsymbol{z},\boldsymbol{Y}_{\boldsymbol{Z},\boldsymbol{\delta}}(0))\geq
c\right\\},$ (3)

and the corresponding randomization-based exact $p$-value is obtained by
evaluating the tail probability in (3) at the observed value of the test
statistic:

$\displaystyle p_{\boldsymbol{Z},\boldsymbol{\delta}}\equiv
G_{\boldsymbol{Z},\boldsymbol{\delta}}\left\\{t\left(\boldsymbol{Z},\boldsymbol{Y}_{\boldsymbol{Z},\boldsymbol{\delta}}(0)\right)\right\\}.$
(4)

The sharp null $H_{0,\text{sharp}}$ specifies the entire ITE vector
$\boldsymbol{\tau}$ and often equates it to a constant (e.g.,
$\boldsymbol{\tau}=c\boldsymbol{1}$). In some circumstances, treatment effect
heterogeneity is likely to exist although its details (e.g., how the treatment
effect varies across subgroups) may be unknown. In these cases, the sharp null
hypothesis of no effect whatsoever, i.e.,
$H_{0,\text{sharp}}\mathrel{\mathop{\ordinarycolon}}\mathbf{Y}(1)=\mathbf{Y}(0)$,
is still relevant and a meaningful first step of data analysis; Imbens and
Rubin, (2015)however, the null hypothesis of a constant additive effect may
not be of primary interest.

Testing a weak null hypothesis emerges as one alternative to testing a sharp
null hypothesis. A weak null hypothesis only hypothesizes one aspect or
property of the collection of ITEs (e.g., the average treatment effect)
without putting any restriction on each component ITE. Much of the literature
on testing a weak null hypothesis, dating back to Neyman and
Iwaszkiewicz,Neyman and Iwaszkiewicz, (1935) has focused almost exclusively on
testing the null hypothesis of no average treatment effect or its variation.
Unlike the FRT that is _exact_ , tests for sample average treatment effects in
general require large-sample approximation, although some of them also enjoy
finite-sample validity for a certain sharp null hypothesis. For instance, Ding
and Dasgupta Ding and Dasgupta, (2018) showed that, under a CRE design, the
FRT for the sharp null hypothesis of no effect using the studentized
$t$-statistic is also asymptotically valid for the weak null hypothesis of no
average treatment effect; see also Cohen and Fogarty Cohen and Fogarty, (2022)
and Wu and Ding Wu and Ding, (2021) for extensions to more general test
statistics and experimental designs.

As discussed in Section 1.1, the SATE is one important, albeit not
comprehensive, assessment of a regimen’s treatment effect profile. Moreover,
when the outcomes have heavy tails, the SATE could be sensitive to the
outliers and the finite population asymptotic approximation tends to work
poorly in these cases.Caughey et al., (2021) These aspects are particularly
relevant in many pre-clinical studies (e.g., nonhuman primates studies) and
early phase clinical trials, where the sample size could be as small as $10$
to $20$ per arm and treatment effect heterogeneity is often expected; see,
e.g., HVTN 086 in Figure 1.

###  2.2 Beyond the SATE: quantiles and proportions

An exclusive focus on the SATE could arrive at unsatisfactory conclusions in
face of large treatment effect heterogeneity. For example, a treatment that is
harmful to a majority of individuals could still have a positive treatment
effect on average, due to possibly some large ITEs on a small proportion of
the cohort; see Figure S1 in Supplemental Material B for an illustration. In
the following, we will study finite-sample valid, randomization-based tests
for inferring the distribution of ITEs. Let $F(c)\equiv
N^{-1}\sum_{i=1}^{N}\mathbbm{1}(\tau_{i}\leq c)$, for $c\in\mathbb{R}$, denote
the distribution function of ITEs and
$F^{-1}(\beta)\equiv\inf\\{c\mathrel{\mathop{\ordinarycolon}}F(c)\geq\beta\\}$,
for $\beta\in(0,1)$, denote the corresponding quantile function. We focus on
the $N$ experimental units at hand so that the quantile function can take at
most $N$ values, i.e., the sorted ITEs
$\tau_{(1)}\leq\tau_{(2)}\leq\ldots\leq\tau_{(N)}$. More precisely, we have
$F^{-1}(\beta)=\tau_{(k)}$, with $k=\lceil N\beta\rceil$ denoting the ceiling
of $N\beta$, for $\beta\in(0,1]$. In addition, the distribution function can
be equivalently written as $F(c)=1-N(c)/N$, for $c\in\mathbb{R}$, where
$N(c)\equiv\sum_{i=1}^{N}\mathbbm{1}(\tau_{i}>c)$ denotes the number of units
with ITEs exceeding a threshold $c$. We consider the following null hypothesis
for any $0\leq k\leq N$ and $c\in\mathbb{R}$:

$\displaystyle H_{k,c}\mathrel{\mathop{\ordinarycolon}}\tau_{(k)}\leq c\
\Longleftrightarrow\ N(c)\leq N-k.$ (5)

For descriptive simplicity, we define $\tau_{(0)}=-\infty$. Here we focus on
one-sided testing with alternatives favoring large treatment effects, based on
which we may then obtain lower confidence bounds for effect quantiles
$F^{-1}(\beta)$’s (or equivalently $\tau_{(k)}$’s) and proportions of units
whose ITEs exceeding any thresholds $1-F(c)$’s (or equivalently $N(c)$’s). To
obtain one-sided tests with alternatives favoring small treatment effects, one
may use the same procedure but flip observed outcomes’ signs or switch the
labels for treatment and control. Two-sided tests can be constructed by
combining two one-sided tests using, say, Bonferroni correction.

##  3 Placebo-controlled trials with highly specific endpoints

###  3.1 Inference for a single effect quantile or proportion

In a placebo-controlled vaccine trial, participants’ potential outcomes under
the placebo are often known _a priori_ , and these controls are nonetheless
included in the study primarily for blinding purposes (e.g., preventing
treatment arm information from being revealed to lab technicians). For
instance, in many early phase HIV vaccine trials, the endpoints of interest
are vaccine-induced, antigen-specific immune responses, and healthy and naïve
study participants receiving the placebo are not expected to exhibit any of
these highly specific immune responses to the antigens. Hence, we often have
auxiliary information $Y_{i}(0)=\textsf{LOD}$ for all
$i=1,\cdots,N_{1}+N_{0}$, where LOD represents an assay-specific limit of
detection. In this stylistic case, entries left unknown in Table 1 are
potential outcomes under treatment for $N_{0}$ placebo recipients. Without
loss of generality, we assume $\textsf{LOD}=0$ so $Y_{i}(0)=0$ for all $i$.
The observed data inform us of ITEs of $N_{1}$ treated units, based on which
we then infer the entire distribution of ITEs.

Consider the null hypothesis $H_{k,c}$ in (5) for any $0\leq k\leq N$ and
$c\in\mathbb{R}$. Under the CRE design, treated units are a simple random
sample of size $N_{1}$ from a total of $N$ experimental units. Therefore,
$n(c)\equiv\sum_{i=1}^{N}Z_{i}\mathbbm{1}(Y_{i}>c)=\sum_{i=1}^{N}Z_{i}\mathbbm{1}(\tau_{i}>c)$
follows a Hypergeometric distribution with parameters $(N,N(c),N_{1})$. The
Hypergeometric distribution with parameters $(N,n,N_{1})$ becomes
stochastically larger as $n$ increases111Specifically, if $X$ follows a
Hypergeometric distribution with parameters $(N,n,N_{1})$, $Y$ follows a
Hypergeometric distribution with parameters $(N,n^{\prime},N_{1})$ such that
$n^{\prime}\geq n$, then $\mathbb{P}\left(X\leq
t\right)\geq\mathbb{P}\left(Y\leq t\right)$ for every $t\in\mathbb{R}$..
Together, these facts imply a finite-sample valid $p$-value for testing
$H_{k,c}$, as summarized in the following proposition.

######  Proposition 1.

For any $1\leq k\leq N$ and $c\in\mathbb{R}$, $p_{k,c}^{\textup{H}}\equiv
G_{\textup{H}}(n(c);N,N-k,N_{1})$ is a valid $p$-value for testing the null
hypothesis $H_{k,c}$ in (5), where $G_{\textup{H}}(x;N,n,N_{1})\equiv\Pr(X\geq
x)$ denotes the tail probability of a Hypergeometric random variable $X$ with
parameters $(N,n,N_{1})$. Specifically, under $H_{k,c}$,
$\Pr(p_{k,c}^{\textup{H}}\leq\alpha)\leq\alpha$ for any $\alpha\in(0,1)$.

###### Proof.

All proofs in the article can be found in Supplemental Material A. ∎

From Proposition 1, we can then conduct Lehmann-style test inversion to
construct confidence intervals for $\tau_{(k)}$ and $n(c)$, for any $1\leq
k\leq N$ and $c\in\mathbb{R}$. Moreover, due to the monotonicity of the
$p$-value $p_{k,c}^{\textup{H}}$ in $k$ and $c$, the resulting confidence sets
are intervals and have simpler forms that facilitate their computation.
Proposition 2 summarizes these results. For descriptive convenience, we let
$y_{(1)}\leq y_{(2)}\leq\ldots\leq y_{(N_{1})}$ denote the sorted observed
outcomes for treated units and further define $y_{(0)}=-\infty$.

######  Proposition 2.

  * (a)

$p_{k,c}^{\textup{H}}$ is monotone increasing in $c$ and decreasing in $k$.

  * (b)

For any $1\leq k\leq N$ and $\alpha\in(0,1)$,
$\\{c\mathrel{\mathop{\ordinarycolon}}p_{k,c}^{\textup{H}}>\alpha,c\in\mathbb{R}\\}$
is a $1-\alpha$ confidence interval for $\tau_{(k)}$, and it can be
equivalently written as $[y_{(k(\alpha))},\infty)$ with $k(\alpha)\equiv
N_{1}-Q_{\textup{H}}(1-\alpha;N,N-k,N_{1})$, where
$Q_{\textup{H}}(\theta;N,n,N_{1})$ denotes the $\theta$-th quantile function
of a Hypergeometric random variable $X$ with parameters $(N,n,N_{1})$.

  * (c)

For any $c\in\mathbb{R}$ and $\alpha\in(0,1)$,
$\\{N-k\mathrel{\mathop{\ordinarycolon}}p_{k,c}^{\textup{H}}>\alpha,0\leq
k\leq N\\}$ is a $1-\alpha$ confidence set for $N(c)$, and it can be
equivalently written as $\\{n_{c,\alpha},n_{c,\alpha}+1,\ldots,N\\}$ with
$n_{c,\alpha}=N-\max\\{k\mathrel{\mathop{\ordinarycolon}}G_{\textup{H}}(n(c);N,N-k,N_{1})>\alpha,0\leq
k\leq N\\}$.

######  Remark 1.

Confidence interval results for $\tau_{(k)}$ can also be derived by applying
results in Sedransk and Meyer, Sedransk and Meyer, (1978) who studied
inference for quantiles of a finite population based on simple random
sampling. In Supplemental Material A.3, we prove the equivalence between the
confidence interval results in Proposition 2 and those derived by applying
results in Sedransk and Meyer. Sedransk and Meyer, (1978) In Supplemental
Material A.3, we also prove that the confidence interval results for $N(c)$ in
Proposition 2 are equivalent to results in Wang Wang, (2015) who studied
inference for Hypergeometric distribution parameters. Hence, by the optimality
results in Wang, Wang, (2015) confidence intervals for $N(c)$ in Proposition
2(c) are optimal one-sided confidence intervals.

###  3.2 Simultaneous confidence region for the quantile function and
distribution function

Evaluating the ITE profile of a treatment regimen goes beyond constructing
confidence intervals for each individual ITE. Proposition 3 derives the
simultaneous confidence region for multiple ITE quantiles, which provides a
valid summary of treatment effect profile and is of primary interest in
practice.

######  Proposition 3.

Consider a CRE design and suppose that we are interested in multiple quantiles
of ITEs: $\tau_{(k_{1})},\tau_{(k_{2})},\ldots,\tau_{(k_{J})}$, where $1\leq
k_{1}\leq\ldots\leq k_{J}\leq N$. Under a CRE design, the simultaneous
coverage probability for quantiles of ITEs is

$\displaystyle\Pr\left(\bigcap_{j=1}^{J}\\{\tau_{(k_{j})}\geq
y_{(k_{j}(\alpha))}\\}\right)$ $\displaystyle\geq
1-\Pr\left(\bigcup_{j=1}^{J}\left\\{\sum_{i=1}^{N}Z_{i}\mathbbm{1}(i>k_{j})>Q_{\textup{H}}(1-\alpha;N,N-k_{j},N_{1})\right\\}\right),$
(6)

where the equality holds when all ITEs $\tau_{i}$’s are distinct. Recall that
$y_{(k_{j}(\alpha))}$ is the observed outcome of rank $k_{j}(\alpha)\equiv
N_{1}-Q_{\textup{H}}(1-\alpha;N,N-k_{j},N_{1})$ in the treatment group.

The lower bound of the simultaneous coverage probability in (6) is sharp and
is obtained when all ITEs are distinct across all $N$ units. More importantly,
the lower bound does not depend on any unknown quantities and can be
efficiently approximated using the Monte Carlo method.

Technically, the assay limit of detection specifies the lowest level of immune
response to be detected and merely places an upper bound on $Y_{i}(0)$ rather
than specifies its precise value. Proposition 4 extends previous results and
suggests that $p$-values and confidence statements derived from Propositions
1-3 remain valid when $Y_{i}(0)\leq 0$ rather than $Y_{i}(0)=0$.

######  Proposition 4.

Propositions 1-3 hold when $Y_{i}(0)\leq 0$.

##  4 Treatment effect against a second regimen

In a head-to-head comparison of two active regimens, we do not have the _a
priori_ knowledge that the unit-level potential outcomes under a second
regimen are always equal to a constant. In a recent article, Caughey et al.
Caughey et al., (2021) extended the FRT to test the null hypothesis $H_{k,c}$
in (5) without imposing any restriction on potential outcomes. Because
$H_{k,c}$ is a composite null hypothesis and permits infinitely many
imputation schemes, FRT is not directly applicable. Nevertheless, a valid
$p$-value for testing $H_{k,c}$ can be obtained by maximizing the
randomization $p$-value $p_{\boldsymbol{Z},\boldsymbol{\delta}}$ in (4) over
$\boldsymbol{\delta}\in\mathcal{H}_{k,c}$, where $\mathcal{H}_{k,c}$ denotes
the set of vectors whose elements of rank $k$ are bounded by $c$, i.e.,
$\mathcal{H}_{k,c}\equiv\left\\{\boldsymbol{\delta}\in\mathbb{R}^{N}\mathrel{\mathop{\ordinarycolon}}\delta_{(k)}\leq
c\right\\}\subset\mathbb{R}^{N}$. However, optimizing
$\sup_{\boldsymbol{\delta}\in\mathcal{H}_{k,c}}p_{\boldsymbol{Z},\boldsymbol{\delta}}$
is computationally challenging and can even be NP hard. To address this
challenge, Caughey et al. Caughey et al., (2021) proposed to utilize the class
of rank score statistics of the following form:

$t(\boldsymbol{z},\boldsymbol{y})=\sum_{i=1}^{N}z_{i}\phi\left(\mathrm{r}_{i}(\boldsymbol{y})\right),$
(7)

where $\phi(\cdot)$ is a monotone increasing function and
$\mathrm{r}_{i}(\boldsymbol{y})$ denotes the rank of the $i$th coordinate of
$\boldsymbol{y}$ using index ordering to break ties, assuming that the
ordering of units has been randomly permuted before analysis.

Under a CRE design, the rank score statistic $t(\cdot,\cdot)$ defined in (7)
is distribution free, in the sense that for any
$\boldsymbol{y},\boldsymbol{y^{\prime}}\in\mathbb{R}^{N}$,
$t(\boldsymbol{Z},\boldsymbol{y})$ and
$t(\boldsymbol{Z},\boldsymbol{y^{\prime}})$ follow the same distribution,
where $\boldsymbol{Z}$ follows the treatment assignment mechanism. Because of
the distribution free property of rank score statistics, the imputed
randomization distribution in (3) reduces to a distribution that does not
depend on the observed treatment assignment $\boldsymbol{Z}$ or the
hypothesized treatment effect $\boldsymbol{\delta}$, i.e.,

$G_{\boldsymbol{Z},\boldsymbol{\delta}}(c)\equiv\operatorname{Pr}\left\\{t\left(\boldsymbol{Z},\boldsymbol{Y}_{\boldsymbol{Z},\boldsymbol{\delta}}(0)\right)\geq
c\right\\}=\operatorname{Pr}\\{t(\boldsymbol{Z},\boldsymbol{y})\geq c\\}\equiv
G(c),$ (8)

where $\boldsymbol{y}$ can be any constant vector in $\mathbb{R}^{N}$.
Consequently, the valid $p$-value
$\sup_{\boldsymbol{\delta}\in\mathcal{H}_{k,c}}p_{\boldsymbol{Z},\boldsymbol{\delta}}$
for testing $H_{k,c}$ in (5) simplifies to

$p_{k,c}^{\textup{R}}\equiv\sup_{\boldsymbol{\delta}\in\mathcal{H}_{k,c}}p_{\boldsymbol{Z},\boldsymbol{\delta}}=\sup_{\boldsymbol{\delta}\in\mathcal{H}_{k,c}}G\\{t(\boldsymbol{Z},\boldsymbol{Y}-\boldsymbol{Z}\circ\boldsymbol{\delta})\\}=G\left\\{\inf_{\boldsymbol{\delta}\in\mathcal{H}_{k,c}}t(\boldsymbol{Z},\boldsymbol{Y}-\boldsymbol{Z}\circ\boldsymbol{\delta})\right\\},$
(9)

where the last equality holds because $G$ is monotone decreasing and
$t(\boldsymbol{Z},\boldsymbol{Y}-\boldsymbol{Z}\circ\boldsymbol{\delta})$
achieves its infimum over $\boldsymbol{\delta}\in\mathcal{H}_{k,c}$. That is,
to obtain a valid $p$-value for testing $H_{k,c}$ based on a distribution free
test statistic, it suffices to minimize the value of the test statistic
$t\left(\boldsymbol{Z},\boldsymbol{Y}_{\boldsymbol{Z},\boldsymbol{\delta}}(0)\right)$
over $\boldsymbol{\delta}\in\mathcal{H}_{k,c}$. It can be shown that the
worst-case $p$-value has a closed-form solution when using the rank score
statistics in (7);Caughey et al., (2021) however, the solution places
$Y_{i}(0)$’s of a subset of treated units at $-\infty$ and this can be unduly
conservative in many practical scenarios. For instance, in most vaccine
studies, we could know _a priori_ that $Y_{i}(0)$’s are at least greater than
some pre-specified level (e.g., in a study of non-naïve population) and this
auxiliary information $Y_{i}(0)\geq C$ for some pre-specified constant $C$,
can be leveraged to improve the power of the method in Caughey et al.Caughey
et al., (2021) Without loss of generality, we let $C=0$. Proposition 5 extends
the results in Caughey et al.Caughey et al., (2021) by incorporating this
auxiliary information.

######  Proposition 5.

Define
$\mathcal{H}_{k,c}^{0}\equiv\left\\{\boldsymbol{\delta}\in\mathbb{R}^{N}\mathrel{\mathop{\ordinarycolon}}\delta_{(k)}\leq
c,\delta_{i}\leq Y_{i}\text{ for
}i\in\mathcal{T}\right\\}\subset\mathbb{R}^{N}$,
$\mathcal{T}_{c}=\left\\{i\mathrel{\mathop{\ordinarycolon}}Z_{i}=1,Y_{i}-c>0,1\leq
i\leq\right.$ $N\\}$ and $l=\min(N-k,\left|\mathcal{T}_{c}\right|)$. Let
$J_{k}$ be an index set such that $\mathcal{J}_{k}\subset\mathcal{T}_{c}\text{
and }|\mathcal{J}_{k}|=l$. Under a CRE, for any rank score statistic
$t(\cdot,\cdot)$, any $1\leq k\leq N$ and any constant $c\in\mathbb{R}$,

$p_{k,c}^{{\textup{R},0}}\equiv\sup_{\boldsymbol{\delta}\in\mathcal{H}_{k,c}^{0}}p_{\boldsymbol{Z},\boldsymbol{\delta}}=G\left\\{\inf_{\boldsymbol{\delta}\in\mathcal{H}_{k,c}^{0}}t\left(\boldsymbol{Z},\boldsymbol{Y}_{\boldsymbol{Z},\boldsymbol{\delta}}(0)\right)\right\\}=G\left\\{\inf_{\mathcal{J}_{k}\subset\mathcal{T}_{c}\text{
and
}|\mathcal{J}_{k}|=l}t(\boldsymbol{Z},\boldsymbol{Y}-\boldsymbol{Z}\circ\boldsymbol{\eta})\right\\}$
(10)

is a valid $p$-value for testing the null hypothesis $H_{k,c}$ in (5) with the
auxiliary information $Y_{i}(0)\geq 0$ for all $1\leq i\leq N$, where
$\boldsymbol{\eta}=\left(\eta_{1},\eta_{2},\ldots,\eta_{N}\right)^{\intercal}\in\mathbb{R}^{N}$
satisfies $\eta_{i}=c$ if $i\in\mathcal{T}_{c}\backslash\mathcal{J}_{k}$ and
$Y_{i}$ otherwise. When $k=N$, $\mathcal{J}_{N}$ is an empty set.
Specifically, under $H_{k,c}$,
$\operatorname{Pr}\left(p_{k,c}^{{\textup{R},0}}\leq\alpha\right)\leq\alpha$
for any $\alpha\in(0,1)$.

In Proposition 5, the worse-case $p$-value is obtained when the largest $N-k$
ITEs are equal to the observed outcomes (i.e., $\tau_{i}$ is maximized subject
to $Y_{i}(0)\geq 0$), and these $N-k$ maximal ITEs are assigned to the treated
units with $Y_{i}>c$ in such a way that the test statistic is minimized. To
derive such a worst-vase $p$-value, one needs to enumerate
$\left|\mathcal{T}_{c}\right|\choose N-k$ possibilities and calculate the
corresponding test statistics. This is in general a mild task for a small-size
clinical trial where randomization inference is most appropriate. The solution
is different from the original solution in Caughey et al. Caughey et al.,
(2021) in that, in the absence of the auxiliary information, they let ITEs of
the largest $N-k$ equal infinity and simply assign these infinity values to
the $\min(N_{1},N-k)$ treated units with the largest observed outcomes. In
other words, the worst-case solution in Caughey et al. Caughey et al., (2021)
admits a more explicit, closed-form solution. Nevertheless, incorporating the
auxiliary information, Proposition 5 derives a less conservative $p$-value and
ultimately delivers valid but more powerful statistical inference.

Proposition 6 builds upon Proposition 5 and derives simultaneous inference for
multiple ITE quantiles with the auxiliary information $Y_{i}(0)\geq 0$.

######  Proposition 6.

Suppose $Y_{i}(0)\geq 0$ for all $1\leq i\leq N$. Under a CRE and using the
$p$-value $p_{k,c}^{{\textup{R},0}}$ in (10) with any rank score statistic
$t(\cdot,\cdot)$, for any $\alpha\in(0,1)$ and any $1\leq k\leq N$, a
$1-\alpha$ confidence set for $\tau_{(k)}$ is
$\mathcal{S}_{k,c}\equiv\left\\{c\mathrel{\mathop{\ordinarycolon}}p_{k,c}^{{\textup{R},0}}>\alpha,c\in\mathbb{R}\right\\}$.
Moreover, these confidence sets cover the corresponding individual treatment
effects simultaneously with at least $1-\alpha$ probability, i.e.,
$\operatorname{Pr}\left(\tau\in\bigcap_{k=1}^{n}\mathcal{S}_{k,c}\right)\geq
1-\alpha.$

More recently, Chen and LiChen and Li, (2023) enhanced the approach in Caughey
et al. Caughey et al., (2021) and obtained improved statistical power. Instead
of presuming that all the largest individual treatment effects are assigned to
the treated units with the largest outcomes, they consider using the $p$-value
when the number of treated units with effects greater than $c$ is bounded by
$l$. This approach better captures the underlying assignment mechanism of a
randomized experiment, in the sense that it allows units with the upper $k$th
quantile effects to be randomly assigned to each regimen. As shown by Chen and
Li,Chen and Li, (2023) this approach delivers more powerful statistical
inference compared to the method in Caughey et al.,Caughey et al., (2021)
while it can still preserve the testing level by adding a correction term
which bounds, from the above, the probability that the number of treated units
with effects greater than $c$ exceeds $l$. We summarize their results in
Supplemental Material C.1 for completeness.

##  5 Simulation study

###  5.1 Goal and structure; methods to be compared; measurement of success

One primary objective in a typical experimental vaccine trial is to
_characterize_ the immunogenicity profile of each candidate regimen. The
primary goal of the simulation studies is to assess and compare the power of
several competing methodologies in characterizing one regimen. Our simulation
studies can be compactly summarized by a $3\times 2\times 2\times 2\times 3$
factorial design. Factors 1 and 2 specify the clinical trial design and are
meant to mimic the size of a typical early phase or experimental medicine
trial:

Factor 1:

sample size, $N$: $40$, $100$ and $200$;

Factor 2:

proportion assigned to treatment, $p\equiv N_{1}/N$: $0.5$ (balanced) and
$0.8$ (unbalanced).

Factor 3 specifies potential outcomes under control:

Factor 3:

potential outcome under control, $Y_{i}(0)$: $Y_{i}(0)=2$ (corresponding to
LOD in our case study) and $Y_{i}(0)\sim\text{Normal}(0,1)$.

We will refer to the $Y_{i}(0)=2$ case as Scenario I and the
$Y_{i}(0)\sim\text{Normal}(0,1)$ case as Scenario II.

We consider a simple pattern for the ITE profile: $\alpha\times N$ units have
a null ITE ($\tau_{i}=0$) and the other $(1-\alpha)\times N$ positive
responders have a positive ITE ($\tau_{i}>0$). The distribution of $\tau_{i}$
among $(1-\alpha)\times N$ positive responders follows a truncated normal
distribution with parameters $\mu$, $\sigma=1.5$, and truncation at $0.1$ and
$4$. This simple ITE pattern is meant to mimic a practical setting in HIV
vaccine development where the experimental vaccine only induces an antigen-
specific immune response for a subset of study participants. Factors 4 and 5
specify $\alpha$ and $\mu$.

Factor 4:

proportion of units exhibiting a null ITE, $\alpha$: $0.2$ and $0.5$;

Factor 5:

mean parameter of truncated normal among positive responders, $\mu$: $1.5$,
$2.0$ and $2.5$.

Each unit has its potential outcome $Y_{i}(1)=Y_{i}(0)+\tau_{i}$, and the
observed outcome $Y_{i}$ satisfies $Y_{i}=Z_{i}\cdot Y_{i}(1)+(1-Z_{i})\cdot
Y_{i}(0)$.

For each of the $72$ data generating processes considered in Section 5.1, we
perform randomization inference for the distribution function of ITEs (and
equivalently the proportion of units with ITEs exceeding different
thresholds). We construct (i) individual confidence intervals for selected
quantiles; and (ii) simultaneous confidence region for selected quantiles
using the following methods: M1: The method outlined in Section 3; M2: The
method in Caughey et al. Caughey et al., (2021) as outlined in Section 4; in
particular, we considered the Wilcoxon signed rank test (M2-Wilcox) and
Stephenson rank sum test with $s=4$ (M2-S4), $10$ (M2-S10) and $20$ (M2-S20);
M3: The method in Chen and Li Chen and Li, (2023) as reviewed in Supplemental
Material C.1; analogous to the method M2, we considered the Wilcoxon signed
rank test (M3-Wilcox) and three types of Stephenson rank sum test (M3-S4,
M3-S10 and M3-S20). Method M1 is valid only in the constant $Y(0)$ case
(Scenario I) while methods M2 and M3 yield valid inference irrespective of the
nature of $Y(0)$. We repeat simulation 1000 times.

We measure the success of each method in recovering the ITE profile according
to the following criteria. First, for a selected quantile $k$ and the target
estimand $\tau_{(k)}$, we report the proportion of times out of MC simulations
that the one-sided confidence limit derived by each method excludes $c$. This
proportion is denoted by $\textsf{PP}_{k,c}$ and we report $\textsf{PP}_{k,c}$
over a range of $c$ for $k=\lceil 0.95\times N\rceil$, $\lceil 0.9\times
N\rceil$ and $\lceil 0.75\times N\rceil$. For a fixed $k$ and $c$, a larger
$\textsf{PP}_{k,c}$ corresponds to a method being more powerful.

Second, we compare the simultaneous confidence region to the true ITE profile
$\mathcal{T}=\\{\tau_{(i)},~{}i=1,\dots,N\\}$ by calculating
$\textsf{SS}=|\mathcal{K}|^{-1}\sum_{k\in\mathcal{K}}(\hat{L}_{(k)}-\tau_{(k)})^{2}$,
where $\mathcal{K}$ is a collection of quantiles of interest, and
$\hat{L}_{(k)}$ denotes the lower confidence limit delivered by the method. We
consider $\mathcal{K}=\\{\lceil 0.5\times N\rceil,\lceil 0.75\times
N\rceil,\lceil 0.8\times N\rceil,\lceil 0.85\times N\rceil,\lceil 0.9\times
N\rceil,\lceil 0.95\times N\rceil\\}$. We will report SS averaged over MC
simulations. A smaller SS value corresponds to a method being more powerful.

###  5.2 Simulation results

Figure 2 summarizes the performance of each method in Scenario I when
$(n,p,\alpha,\mu)=(100,0.5,0.2,2)$. Results for other data-generating
processes can be found in Figures S2 and S3 in the Supplemental Material D.
Each panel of Figure 2 plots the cumulative distribution functions of the
lower confidence limits of $95\%$ confidence intervals for selected quantiles
of ITEs. The Wilcoxon rank sum statistic and the Stephenson rank sum statistic
with $s=4$ have little power detecting ITEs and are omitted from plots. Recall
that the average SS summarizes a method’s power to identify the true ITE
profile. Table 2 reports the average SS of each method under different data-
generating processes in Scenario I. We identified several consistent trends in
Scenario I. First, each of the three methods had improved power, reflected by
larger $\textsf{PP}_{k,c}$ and smaller average SS, as the sample size $N$
increases. Second, M1 provided the most informative confidence intervals for
effect quantiles with the largest $\textsf{PP}_{k,c}$ and smallest average SS
in all simulation settings. This is as expected because, unlike M2 or M3, M1
incorporates the auxiliary information that $Y_{i}(0)$ is constant. Third, for
M2 and M3, $\textsf{PP}_{k,c}$ increases with $s$ and the average SS decreases
with $s$, that is, larger $s$ tends to yield more informative lower confidence
limits for effect quantiles.

In Scenario II, $Y_{i}(0)$ is not constant so M1 is not applicable and we
focus on investigating the power of M2 and M3. Figure S4 in Supplemental
Material D shows the cumulative distribution function of the $95\%$ lower
confidence limits for selected ITE quantiles when
$(n,p,\alpha,\mu)=(100,0.5,0.2,2.5)$. Table S1 in Supplemental Material D
further summarizes the average SS. Interestingly, when $Y(0)$ is not constant,
a large value of $s$ may deteriorate power and increase the average SS. Based
on these simulation results, we recommend a moderately large $s$, say $s=10$.
With a fixed $s$, M3 appeared to be superior to M2 in most simulation
settings. In particular, our simulation studies suggest that M3 outperforms M2
for both pointwise and simultaneous inference in the following settings: (i)
Inferring lower quantiles. As shown in Figure 2, M2 has negligible power when
making inference for the median individual treatment effect. On the other
hand, M3 (with $s=10$) detects a positive median ITE 59% times, i.e.,
$\text{PP}_{0.5,0}=0.59$. For quantiles larger than the median, both M2 and M3
yield informative confidence intervals, and their performance becomes more
similar when inferring larger quantiles. For instance, when $k=\lceil
0.95\times N\rceil$, i.e., the $95\%$ quantile, M2 and M3 provide exactly the
same distribution of the lower confidence limits: the CDFs, green for M2 and
blue for M3, completely overlap in the $95\%$ quantile panel in Figure 2. (ii)
With a larger sample size $N$. As shown in Figure S2 in Supplemental Material
D, with a small sample size $N=40$, the lower confidence limits for the median
ITE of treatment effects are not informative (almost all zeros) for both M3
and M2; on the other hand, with larger sample size $N=100$, M3 is notably
better compared to M2, as shown in Figure 2. Moreover, upon examining Table 2
and Table S1, in a balanced design with a sample size of $N=100$, M3 appears
to be more favorable than M2 with a smaller average SS. (iii) In a balanced
design ($p=0.5$). In a balanced design, M3 outperforms M2 in both the
individual and simultaneous inference for quantiles of treatment effects; see
Figures S2 and S3. Moreover, the overall average SS of M3 is also smaller than
that of M2 with $p=0.5$; see Table 2 and Table S1. On the other hand, in an
unbalanced design with $p=0.8$, M3 is slightly inferior to M2 in making
simultaneous inference: M3 has a slightly higher average SS. This implies
that, in an unbalanced design, the improvement of M3 over M2 in making
inference for individual ITEs may fail to compensate the cost of multiple
testing corrections.

###  5.3 Additional simulation studies

A second objective in vaccine studies is to rank multiple regimens and
recommend the most promising ones to enter the next phase of clinical trials.
In an additional simulation study, we compared ranking two regimens. Regimen I
induced a small and homogeneous treatment effect on a majority of units and a
large and heterogeneous treatment effect on a small fraction of units. Regimen
II induced a large and heterogeneous treatment effect on a larger fraction of
units. We found that, by making inference for the proportion of units with
ITEs exceeding a large threshold, the method M3 (chosen because of its
superior performance in Section 5.2) would recommend Regimen II over Regimen I
in most scenarios; in the remaining scenarios, the method has no power
differentiating between the two regimens. See Supplemental Material E for
details.

(a) Inference for a single quantile of treatment effects

(b) Simultaneous inference for multiple quantiles of treatment effects

Figure 2: Cumulative distribution functions of the (a) individual and (b) simultaneous lower confidence limits of $\tau_{(k)}$’s with $Y_{i}(0)=2$, $N=100$, $p=0.5$, $\alpha=0.2$, $\mu=2.5$ for $k=\lceil 0.5\times N\rceil,\lceil 0.75\times N\rceil,\lceil 0.8\times N\rceil,\lceil 0.85\times N\rceil,\lceil 0.9\times N\rceil,\text{and }\lceil 0.95\times N\rceil$. Table 2:  SS averaged over 1000 simulations with $Y_{i}(0)=2$ for various sample size $n$, treatment proportion $p$, proportion of units with a null treatment effect $\alpha$, mean of truncated normal among positive responders $\mu$, and inference methods for ITE distributions. $\alpha$ | $\mu$ | M1 | M2 | M3
---|---|---|---|---
M2-Wilcox | M2-S4 | M2-S10 | M2-S20 | M3-Wilcox | M3-S4 | M3-S10 | M3-S20
|  |  | n = 40, p = 0.5
0.2 | 1.50 | 0.37 | 6.38 | 4.36 | 2.77 | 2.23 | 6.55 | 4.62 | 2.56 | 2.09
2.00 | 0.37 | 7.58 | 5.01 | 3.05 | 2.44 | 7.80 | 5.26 | 2.76 | 2.22
2.50 | 0.36 | 8.63 | 5.50 | 3.24 | 2.55 | 8.93 | 5.76 | 2.82 | 2.26
0.5 | 1.50 | 0.59 | 4.75 | 4.54 | 3.58 | 3.07 | 4.75 | 4.69 | 3.76 | 3.25
2.00 | 0.61 | 5.85 | 5.58 | 4.34 | 3.68 | 5.85 | 5.77 | 4.58 | 3.90
2.50 | 0.65 | 6.75 | 6.43 | 4.96 | 4.23 | 6.75 | 6.65 | 5.23 | 4.45
|  |  | n = 40, p = 0.8
0.2 | 1.50 | 0.08 | 6.48 | 4.53 | 2.71 | 2.39 | 6.52 | 4.75 | 2.78 | 2.46
2.00 | 0.08 | 7.53 | 5.05 | 2.86 | 2.50 | 7.58 | 5.30 | 2.93 | 2.57
2.50 | 0.07 | 8.50 | 5.41 | 2.87 | 2.49 | 8.59 | 5.73 | 2.94 | 2.56
0.5 | 1.50 | 0.12 | 4.77 | 4.77 | 4.17 | 3.84 | 4.77 | 4.77 | 4.25 | 3.97
2.00 | 0.13 | 5.82 | 5.82 | 5.07 | 4.65 | 5.82 | 5.82 | 5.17 | 4.83
2.50 | 0.12 | 6.88 | 6.87 | 5.86 | 5.32 | 6.88 | 6.88 | 6.01 | 5.57
|  |  | n = 100, p = 0.5
0.2 | 1.50 | 0.18 | 6.21 | 4.24 | 2.22 | 1.53 | 6.10 | 3.52 | 1.66 | 1.01
2.00 | 0.17 | 7.20 | 4.78 | 2.38 | 1.64 | 7.06 | 3.82 | 1.74 | 1.03
2.50 | 0.14 | 8.07 | 5.13 | 2.42 | 1.67 | 7.86 | 3.90 | 1.72 | 0.99
0.5 | 1.50 | 0.28 | 4.97 | 4.70 | 3.16 | 2.30 | 4.97 | 4.69 | 2.80 | 1.70
2.00 | 0.28 | 6.05 | 5.70 | 3.75 | 2.70 | 6.05 | 5.68 | 3.27 | 1.92
2.50 | 0.27 | 7.08 | 6.63 | 4.26 | 3.05 | 7.08 | 6.61 | 3.62 | 2.04
|  |  | n = 100, p = 0.8
0.2 | 1.50 | 0.04 | 5.91 | 3.66 | 1.75 | 1.16 | 6.03 | 3.69 | 1.81 | 1.13
2.00 | 0.04 | 6.83 | 3.99 | 1.83 | 1.23 | 6.98 | 4.00 | 1.89 | 1.14
2.50 | 0.03 | 7.53 | 4.08 | 1.83 | 1.27 | 7.71 | 4.08 | 1.88 | 1.09
0.5 | 1.50 | 0.06 | 4.98 | 4.91 | 3.01 | 1.83 | 4.98 | 4.95 | 3.17 | 1.99
2.00 | 0.05 | 6.01 | 5.93 | 3.46 | 2.03 | 6.01 | 5.98 | 3.66 | 2.20
2.50 | 0.05 | 7.08 | 6.96 | 3.90 | 2.21 | 7.08 | 7.03 | 4.13 | 2.36
|  |  | n = 200, p = 0.5
0.2 | 1.50 | 0.09 | 6.02 | 3.97 | 1.98 | 1.34 | 5.66 | 2.92 | 1.25 | 0.64
2.00 | 0.09 | 7.05 | 4.51 | 2.14 | 1.46 | 6.59 | 3.14 | 1.30 | 0.64
2.50 | 0.08 | 7.83 | 4.81 | 2.17 | 1.52 | 7.23 | 3.13 | 1.28 | 0.58
0.5 | 1.50 | 0.14 | 5.04 | 4.66 | 2.94 | 2.04 | 5.04 | 4.53 | 2.17 | 1.07
2.00 | 0.13 | 6.09 | 5.58 | 3.45 | 2.40 | 6.09 | 5.43 | 2.45 | 1.13
2.50 | 0.12 | 7.17 | 6.49 | 3.92 | 2.71 | 7.17 | 6.29 | 2.63 | 1.11
|  |  | n = 200, p = 0.8
0.2 | 1.50 | 0.02 | 5.52 | 3.19 | 1.50 | 0.83 | 5.53 | 3.14 | 1.38 | 0.69
2.00 | 0.02 | 6.31 | 3.40 | 1.55 | 0.86 | 6.31 | 3.33 | 1.40 | 0.65
2.50 | 0.02 | 6.89 | 3.42 | 1.57 | 0.86 | 6.88 | 3.33 | 1.37 | 0.58
0.5 | 1.50 | 0.03 | 5.05 | 4.77 | 2.54 | 1.26 | 5.05 | 4.86 | 2.48 | 1.16
2.00 | 0.03 | 6.12 | 5.75 | 2.90 | 1.34 | 6.12 | 5.86 | 2.80 | 1.20
2.50 | 0.03 | 7.16 | 6.66 | 3.16 | 1.34 | 7.16 | 6.81 | 3.02 | 1.16

##  6 Application to early phase HIV-1 vaccine trials

###  6.1 Characterizing immunogenicity profiles

Heterogeneity in vaccine-induced immune responses has been widely observed
among vaccinees. For instance, as shown in Figure 1, vaccine-induced immune
responses are heterogeneous among participants receiving 4 vaccine regimens in
the HVTN 086 study, spanning from no response (“nonresponders”) to individuals
with exceptional responses (“high-responders”). Recently, Huang et al. Huang
et al., (2022) conducted a comprehensive meta-analysis exploring the
variations in immune responses induced by $26$ vaccine regimens. Here we focus
on the immunogenicity data from HVTN 086, a phase-1 clinical trial studying 4
vaccine regimens: SAAVI MVA-C priming with sequential or concurrent Novartis
subtype C gp140/MF59 vaccine boost; SAAVI DNA-C2 priming with SAAVI MVA-C
boosting, with or without Novartis subtype C gp140/MF59 vaccine. The HVTN 086
study has 4 study arms. Each arm enrolled $46$ HIV-negative, healthy, vaccine-
naïve adult participants between 18 and 45 in the Republic of South Africa,
among whom 38 were randomly assigned to one candidate vaccine regimen and the
other 8 to the placebo. Our analysis followed Huang et al. Huang et al.,
(2022) and focused on data from study participants who successfully completed
all study visits and received all scheduled vaccination.

Our first goal is to characterize each vaccine regimen’s immunogenicity
profile of antigen-specific immune responses. The primary outcome of interest
is the serum IgG response to the antigen Con6 gp120/B measured by a validated
binding antibody multiplex assay (BAMA) $2$ weeks post the last vaccination.
Because study participants are all naïve to the antigen, we assume that their
potential outcomes under placebo are always less than or equal to $100$, the
limit of detection of the BAMA assay, and use method M1 to conduct both
individual and simultaneous inference for all quantiles of ITEs.

Figure 3 plots the $95\%$ lower confidence limits for selected quantiles of
ITEs. The vertical dashed red lines in Figure 3 show the $95\%$ lower
confidence limits of the mean difference between each regimen and placebo
derived from the standard t-test. Inference for the effect quantiles largely
enriches the data summary based on the t-test; for instance, in addition to
stating that the lower confidence limit of the mean treatment effect of the
vaccine regimen T1 (comparing to the placebo) is $2.15$, researchers could
further state with confidence that the largest individual treatment effect is
at least $2.49$, the top $25\%$ ITE is at least $2.40$, the median ITE is at
least $2.30$, etc. Note that the t-test requires additional distributional
assumptions on the potential outcomes and large sample approximation, which
may be violated due to heavy-tailed outcome distribution and insufficient
sample size. Figure 3 also suggests that the simultaneous confidence intervals
are very similar to the point-wise confidence intervals, which agrees well
with our simulation studies and confirms again that conducting simultaneous
inference according to Proposition 3 does not sacrifice much power.

ITE profiles also facilitate comparing and ranking 4 candidate regimens. From
Figure 3, the lower confidence limits of the median ITE for all 4 regimens,
except T3, exceed 0. If the median ITE is deemed as an important summary of
immunogenicity, then T3 may be the least favorable. Compared to those of T2
and T3, the lower confidence limits of effect quantiles for T1 and T4 are
significantly larger, which suggests that T1 and T4 elicit strong immune
response among more participants and are more favorable compared to T2 or T3.
We may rank 4 regimens via making inference for $N(c)$, the number of
participants with individual treatment effects exceeding $c$. Table S6 in
Supplemental Material F summarizes the $95\%$ lower confidence limits of
$N(c)$ for each regimen and selected values of $c$. According to Table S6, a
$95\%$ confidence interval for $N(2)$ of T1 is $[31,41]$, suggesting that we
are $95\%$ confident that at least $31/41=75.6\%$ participants had a treatment
effect at least as large as $2$ (in the $\log_{10}$-scale) when receiving T1.
Similarly, we are $95\%$ confident that at least $17/31=54.8\%$ participants
had a treatment effect as large as $2$ (in the $\log_{10}$-scale) when
receiving T4. On the other hand, the lower confidence limits of $N(2)$ equal
$0$ for T2 and T3. These results suggest that if researchers are interested in
advancing a vaccine regimen that can elicit high immune responses among most
participants, then T1 is the most promising based on this early phase clinical
trial.

(a) Inference for a single quantile of treatment effects

(b) Simultaneous inference for multiple quantiles of treatment effects

Figure 3: Vaccine regimens versus placebo: (a) individual and (b) simultaneous
lower confidence limits of $\tau_{(k)}$’s for $k=\lceil 0.5\times
N\rceil,\lceil 0.55\times N\rceil,\lceil 0.6\times N\rceil,\lceil 0.65\times
N\rceil,\lceil 0.7\times N\rceil,\lceil 0.75\times N\rceil,\lceil 0.8\times
N\rceil,\lceil 0.85\times N\rceil,\lceil 0.9\times N\rceil,\text{and }\lceil
0.95\times N\rceil$.

###  6.2 Comparison of two non-randomized treatment regimens

Next, we consider a head-to-head comparison of T1 and T3. Comparing the
treatment effect of T1 versus T3 is more subtle because we no longer have
randomization by design. In fact, a naïve comparison of T1 versus T3 could
suffer from a version of the sample selection bias, Stuart et al., (2011);
Dahabreh et al., (2019) which we refer to as “study arm selection bias,” if
participants were self-selected into the study arms T1/C1 or T3/C3. To
alleviate this selection bias, we used statistical matching to control for
observed baseline characteristics. Specifically, we used an optimal tripartite
matching methodZhang et al., (2023) and constructed 24 matched pairs, each
consisting of one study participant receiving vaccine regimen T1 and the other
receiving T3. The matching algorithm matched exactly on sex, closely matched
on age and body mass index (BMI), and minimized the earth mover’s distance
between propensity score distributions of two groups. Table S7 in Supplemental
Material F summarizes the covariate balance; after matching, two groups are
well balanced in age, sex, and BMI, with the standardized mean differences all
below 0.1, or one-tenth of one pooled standard deviation.

(a) Randomization inference with the stratified Wilcoxon rank sum statistic
(b) Sensitivity analysis

Figure 4: $95\%$ lower confidence limits for quantiles of ITEs under
randomization inference and sensitivity analysis.

We first assume that, after matching, two study participants in each pair have
the same probability of receiving T1 or T3; in other words, the matched
observational study reduces to a finely stratified randomized
experiment.Rosenbaum, (2002); Chen et al., (2023) We conducted randomization
inference for quantiles of ITEs using the stratified Wilcoxon rank sum
statistic;Su and Li, (2023) see Supplemental Material C.2 for details. Figure
4(a) shows the simultaneous lower confidence limits for selected quantiles as
follows: for each point in Figure 4(a), its x-axis value denotes the $95\%$
lower confidence limit for the ITE at its y-axis rank. According to Figure
4(a), the $95\%$ lower confidence limit for ITE at rank $40$ barely exceeds
$0$. This implies that, at $95\%$ confidence level, at least $18\%$ (rank $40$
to $48$) of all study participants have a positive effect, or equivalently T1
would induce higher MFI level than T3 for at least $18\%$ of study
participants in the study. The $95\%$ individual lower confidence limits of
$N(c)$ for T1 (n = $24$) versus T3 (n = $24$), with $c=0,0.5,1$ and $1.5$, are
$9,7,6$ and $4$, respectively.

As is true for all non-randomized studies, the head-to-head comparison between
T1 and T3 could suffer from unmeasured confounding. For instance, this would
be the case if healthier people tend to participate in the study arm T1 rather
than T3, and any treatment effect comparing T1 to T3 could be attributed to
this hidden bias. Rosenbaum, (2002) Figure 4(b) shows the sensitivity analysis
results and plots the lower confidence limits for effect quantiles under
Rosenbaum’s sensitivity analysis model indexed by $\Gamma$. Rosenbaum, (2002)
Values of $1.17,2.08,3.48$ and $8.87$ are the largest $\Gamma$ values such
that the resulting confidence intervals for the $85\%,90\%,95\%$ and $100\%$
ITE quantiles do not cover $0$. For example, when the bias in the treatment
assignment is as large as $\Gamma=2.08$, the ITE at rank $44$, i.e.,
approximately $90\%$ quantile, is positive at significance level $0.05$, which
implies that T1 would induce a higher MFI response than T3 for at least $10\%$
of study participants in the study even under a bias of magnitude $2.08$.
According to the method described in Rosenbaum and Silber, Rosenbaum and
Silber, (2009) an unobserved covariate associated with a $3.5$-fold increase
in the odds of selecting in the study arm T1/C1 as opposed to T3/C3 _and_ a
$4.4$-fold increase in the odds of a positive matched-pair difference in MFI
values is needed in order to explain away the ITE at rank $44$. Such a strong
confounding factor appears unlikely, particularly after we have controlled for
observed covariates using matching. Therefore, we can conclude, with
reasonable confidence, that T1 induces higher MFI response compared to T3 for
at least $10\%$ of the cohort, and this causal conclusion is robust to the
treatment arm selection bias.

##  7 Discussion

Many recent research has centered around moving beyond Fisher’s sharp null
hypothesis and exploring aspects of the collection of individual treatment
effects other than the sample mean. In this article, we provide a systematic
review of relevant methods and closely examine the usefulness of different
methods under a wide range of scenarios. We found that randomization inference
that tests quantile treatment effects could be a useful complement to the
SATE, especially when the scientific interest lies in uncovering and
quantifying the heterogeneous treatment effects. In addition, these methods
hold promise in helping relevant stakeholders advance an experimental therapy
that has a meaningfully large treatment effect on possibly a fraction of study
participants, as opposed to a competing therapy that has a similar SATE but
nevertheless does not show any treatment effect at a magnitude of practical
relevance. Another interesting and practically relevant finding is that,
perhaps contrary to expectation, when constructing a simultaneous confidence
region jointly for many quantiles, or even every quantile, the cost of
multiple hypothesis testing could be minimal.

Overall, our assessment is that the suite of randomization-based inferential
methods, testing Fisher’s sharp null of no effect whatsoever, Neyman’s weak
null of no SATE, and various quantiles of ITEs discussed in the current
article, should always have a place in empirical researchers’ toolbox when
analyzing data derived from clinical trials because these methods are always
reliable and could potentially yield informative results.

In many senses, quantiles of ITEs are the most fine-grained estimands. This
important line of research could be furthered in at least two directions.
First, choosing an appropriate test statistic is a key component. As
demonstrated in simulation studies, the power of the Stephenson rank sum
statistic depends critically on the choice of $s$ and an optimal choice of the
test statistic could depend largely on the data-generating process. One
possibility is to develop a data-driven, adaptive approach that produces an
optimal or near-optimal choice of the test statistic (or the tuning parameter
in the test statistic) or combines several test statistics. Second, we have
reviewed and considered methods that (i) assume a constant $Y(0)$; (ii) do not
place any restriction on any aspect of the potential outcomes; (iii) place a
uniform upper or lower bound on potential outcomes. There are still plenty of
other possibilities to leverage auxiliary information from historical data and
improve the power of statistical inference.

## Acknowledgements

Research reported in this publication was partially supported by the National
Science Foundation (award DMS-2238128 to Xinran Li).

## Data Availability Statement

The data that support the findings of this study are available from the
public-facing HIV Vaccine Trials Network (HVTN) website
(https://atlas.scharp.org/).

## References

  * Caughey et al., (2021) Caughey, D., Dafoe, A., Li, X., and Miratrix, L. (2021).  Randomization inference beyond the sharp null: Bounded null hypotheses and quantiles of individual treatment effects.  arXiv preprint arXiv:2101.09195.
  * Chen et al., (2023) Chen, K., Heng, S., Long, Q., and Zhang, B. (2023).  Testing biased randomization assumptions and quantifying imperfect matching and residual confounding in matched observational studies.  Journal of Computational and Graphical Statistics, 32(2):528–538.
  * Chen and Li, (2023) Chen, Z. and Li, X. (2023).  Distributions of individual treatment effects in sampling-based randomized experiments.  Technical Report.
  * Cohen and Fogarty, (2022) Cohen, P. L. and Fogarty, C. B. (2022).  Gaussian prepivoting for finite population causal inference.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(2):295–320.
  * Dahabreh et al., (2019) Dahabreh, I. J., Robertson, S. E., Tchetgen, E. J., Stuart, E. A., and Hernán, M. A. (2019).  Generalizing causal inferences from individuals in randomized trials to all trial-eligible individuals.  Biometrics, 75(2):685–694.
  * Ding and Dasgupta, (2018) Ding, P. and Dasgupta, T. (2018).  A randomization-based perspective on analysis of variance: a test statistic robust to treatment effect heterogeneity.  Biometrika, 105(1):45–56.
  * Ditse et al., (2020) Ditse, Z., Mkhize, N. N., Yin, M., Keefer, M., Montefiori, D. C., Tomaras, G. D., Churchyard, G., Mayer, K. H., Karuna, S., Morgan, C., et al. (2020).  Effect of hiv envelope vaccination on the subsequent antibody response to hiv infection.  Msphere, 5(1):e00738–19.
  * Fisher, (1935) Fisher, R. A. (1935).  The Design of Experiments. Oliver and Boyd. London and Edinburgh.
  * Huang et al., (2022) Huang, Y., Zhang, Y., Seaton, K. E., De Rosa, S., Heptinstall, J., Carpp, L. N., Randhawa, A. K., McKinnon, L. R., McLaren, P., Viegas, E., et al. (2022).  Baseline host determinants of robust human hiv-1 vaccine-induced immune responses: A meta-analysis of 26 vaccine regimens.  Ebiomedicine, 84:104271.
  * Imbens and Rubin, (2015) Imbens, G. W. and Rubin, D. B. (2015).  Causal Inference in Statistics, Social, and Biomedical Sciences.  Cambridge University Press.
  * Lipkovich et al., (2023) Lipkovich, I., Svensson, D., Ratitch, B., and Dmitrienko, A. (2023).  Overview of modern approaches for identifying and evaluating heterogeneous treatment effects from clinical data.  Clinical Trials, page 17407745231174544.
  * Neyman and Iwaszkiewicz, (1935) Neyman, J. and Iwaszkiewicz, K. (1935).  Statistical problems in agricultural experimentation.  Supplement to the Journal of the Royal Statistical Society, 2(2):107–180.
  * Neyman, (1923) Neyman, J. S. (1923).  On the application of probability theory to agricultural experiments. Essay on principles. Section 9.  Annals of Agricultural Sciences, 10:1–51.
  * Rosenbaum, (2002) Rosenbaum, P. R. (2002).  Observational Studies.  Springer.
  * Rosenbaum and Silber, (2009) Rosenbaum, P. R. and Silber, J. H. (2009).  Amplification of sensitivity analysis in matched observational studies.  Journal of the American Statistical Association, 104(488):1398–1405.
  * Rubin, (1974) Rubin, D. B. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of Educational Psychology, 66(5):688.
  * Rubin, (2005) Rubin, D. B. (2005).  Causal inference using potential outcomes: Design, modeling, decisions.  Journal of the American Statistical Association, 100(469):322–331.
  * Sedransk and Meyer, (1978) Sedransk, J. and Meyer, J. (1978).  Confidence intervals for the quantiles of a finite population: simple random and stratified simple random sampling.  Journal of the Royal Statistical Society: Series B (Methodological), 40(2):239–252.
  * Stuart et al., (2011) Stuart, E. A., Cole, S. R., Bradshaw, C. P., and Leaf, P. J. (2011).  The use of propensity scores to assess the generalizability of results from randomized trials.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 174(2):369–386.
  * Su and Li, (2023) Su, Y. and Li, X. (2023).  Treatment Effect Quantiles in Stratified Randomized Experiments and Matched Observational Studies.  Biometrika.  asad030.
  * Wang, (2015) Wang, W. (2015).  Exact optimal confidence intervals for hypergeometric parameters.  Journal of the American Statistical Association, in press.
  * Wu and Ding, (2021) Wu, J. and Ding, P. (2021).  Randomization tests for weak null hypotheses in randomized experiments.  Journal of the American Statistical Association, 116(536):1898–1913.
  * Zhang et al., (2023) Zhang, B., Small, D. S., Lasater, K. B., McHugh, M., Silber, J. H., and Rosenbaum, P. R. (2023).  Matching one sample according to two criteria in observational studies.  Journal of the American Statistical Association, 118(542):1140–1151.