# Improved inference for vaccine-induced immune responses via shape-
constrained methods

Nilanjana Laha  Zoe Moodie  Ying Huang  Alex Luedtke

###### Abstract

We study the performance of shape-constrained methods for evaluating immune
response profiles from early-phase vaccine trials. The motivating problem for
this work involves quantifying and comparing the IgG binding immune responses
to the first and second variable loops (V1V2 region) arising in HVTN 097 and
HVTN 100 HIV vaccine trials. We consider unimodal and log-concave shape-
constrained methods to compare the immune profiles of the two vaccines, which
is reasonable because the data support that the underlying densities of the
immune responses could have these shapes. To this end, we develop novel shape-
constrained tests of stochastic dominance and shape-constrained plug-in
estimators of the squared Hellinger distance between two densities. Our
techniques are either tuning parameter free, or rely on only one tuning
parameter, but their performance is either better (the tests of stochastic
dominance) or comparable with the nonparametric methods (the estimators of the
squared Hellinger distance). The minimal dependence on tuning parameters is
especially desirable in clinical contexts where analyses must be prespecified
and reproducible. Our methods are supported by theoretical results and
simulation studies.

62G07, 62G10, 62G05, 62P10,

vaccine trial,

shape-constraint,

test of stochastic dominance,

Hellinger distance,

logconcave,

unimodal,

###### keywords:

[class=MSC]

###### keywords:

newfloatplacementnewfloatnamenewfloatfileextnewfloatwithin

,
and

##  1 Introduction

To date, the RV144 trial conducted in Thailand is the only vaccine efficacy
trial to show a signal of efficacy (31%) against HIV infection (Rerks-Ngarm et
al., 2009). RV144 inspired a phase 1b trial, named HIV Vaccine Trials Network
(HVTN) 097, which evaluated the safety and immunogenicity of the same regimen
in a South African population (Gray et al., 2019). The predominant subtype of
HIV in South Africa is Clade C. Therefore, in an effort to increase the
potential for high efficacy against clade C infections, scientists modified
the HVTN 097 regimen to include HIV strains matched to the South African clade
C infections (Bekker et al., 2018). A phase 1/2 trial named HVTN 100 assessed
the safety and immunogenicity of the modified regimen in South Africa.
However, the HVTN 702 phase 2B/3 trial of the HVTN 100 regimen in South Africa
met its non-efficacy criteria in February 2020 at a planned interim analysis.

In light of the above, the comparison between the immune response profiles of
HVTN 097 and HVTN 100 trials becomes important because the latter can shed
some light on why the HVTN 100 regimen lacks efficacy. To fix ideas, here we
focus on one class of immune responses, namely the binding of IgG antibodies
to the first and second variable loops (V1V2 region) of the HIV envelope. This
immune response is of particular interest because the RV144 trial revealed an
inverse association between HIV infection and this immune response among
vaccinees (Haynes et al., 2012). Using data from HVTN 097 and HVTN 100, we
focus on answering the following three questions:

  * Q1.

How can we estimate the densities of the aggregated IgG binding immune
responses (to HIV-1 envelope proteins)?

  * Q2.

Is there any ordering between the distributions of the immune responses from
the two trials?

  * Q3.

How can we measure the discrepancy between the densities of the immune
responses from the two trials?

Although Q1 is not directly related with comparison of the two vaccine trials,
answering Q1 is important because the resulting density estimators can help in
designing subsequent vaccines. To answer Q2, we resort to testing for
stochastic dominance, which addresses the ordering of the underlying
distribution functions. To answer Q3, we rely on the squared Hellinger
distance as a measure of discrepancy between two densities, where, for
densities $f$ and $g$, the Hellinger distance ${\mathcal{D}}(f,g)$ is defined
by

${\mathcal{D}}(f,g)=\sqrt{\frac{1}{2}\int_{\mathbb{R}}\bigg{(}\sqrt{f(x)}-\sqrt{g(x)}\bigg{)}^{2}dx}.$

The significance of Q3 may not be immediately obvious. However, the squared
Hellinger distance between the immune responses of HVTN 097 and HVTN 100
trials can serve as a benchmark when new pairs of vaccines are compared in
future vaccine trials. The reason behind choosing the squared Hellinger
distance as the measure of discrepancy in particular is discussed in Section
5.

There are numerous nonparametric methods that can be implemented to carry out
the aforementioned steps. However, traditional nonparametric methods do not
exploit any information on the shape of the underlying densities. The
uniformity of the trial population and exploratory analyses bear evidence that
the underlying densities can be unimodal. The data are also consistent with
the possibility that the densities are log-concave. The latter is an important
subclass of unimodal densities, often advocated for use in modelling because
it contains most of the well-known subexponential unimodal densities and
allows powerful density estimation tools (Walther et al., 2009). Of late,
shape-constrained density estimation has gained much attention. Among other
reasons, the reduced burden of external tuning parameters (Samworth and Sen,
2018; Johnson et al., 2018) makes shape-constrained density estimation an
attractive alternative to traditional nonparametric approaches like kernel-
based or basis expansion type methods, which are known to be sensitive to the
choice of the tuning parameters (cf. pp. 327, Efromovich, 2008; Laha, 2021).
Also, shape-constrained density estimation methods require weaker smoothness
assumptions for asymptotic consistency than do the nonparametric methods.
Shape-constrained techniques are widely used in economics and operations
research (Johnson et al., 2018), and have seen application in other domains
such as circuit design (Hannah and Dunson, 2012). For a detailed account on
the recent development of the shape constraint literature, we refer the survey
articles Samworth and Sen (2018) and Samworth (2018).

Although leveraging shape information can potentially increase efficiency,
there is little to no literature on the application of shape-constrained tools
in vaccine trials. To answer our motivating questions, therefore, we develop
new methods using shape-constrained tools. The application of the new methods
is not limited to vaccine trials.  For example, our tests can be applied to
other areas of medical research where tests of stochastic dominance are
relevant. See Leshno and Levy (2004) for an in-depth discussion of potential
uses of tests of stochastic dominance in medical research. In particular, our
methods are applicable to the mortality data considered in Leshno and Levy
(2004) to infer whether a surgery increases the mortality of patients with
abdominal aortic aneurysma. See also Stinnett and Mullahy (1998) and DeFauw
(2011) for the use of stochastic dominance in cost effective analysis of
healthcare. Outside medical research, our tests have substantial applicability
in finance, economics, social welfare, and operations research, where
stochastic dominance is a popular tool to compare portfolios, income, utility,
poverty, opportunity etc.; see Levy (1992), Sriboonchita et al. (2009), Le
Breton (1991), among others, for a detailed account. On the other hand,
Hellinger distance has also seen successful application as a measure of
discrepancy in various disciplines ranging from machine learning (Cieslak and
Chawla, 2009; González-Castro et al., 2013, 2010) to ecology (Rao, 1995) to
fraud detection (Yamanishi et al., 2004). However, for income data, log-
concavity based methods should be used with caution since distributions with
log-concave density are always sub-exponential (Cule and Samworth, 2010),
whereas income data can have heavier tails. Diagnostic tools such as those in
Asmussen and Lehtomaa (2017) can be used to enquire if the data has heavier
tail, e.g. regular varying tails, in these cases.

###  1.1 Organization of article and main contributions

Although the methods developed in this paper are general, central to our
application lies the HVTN data, which we describe in Section 2. Below we
briefly discuss our methods and main contributions. We write $f_{097}$ and
$f_{100}$ for the densities of the immune responses in the HVTN 097 and HVTN
100 trials respectively, and $F_{097}$ and $F_{100}$ for the corresponding
distribution functions.

### Estimating $f_{097}$ and $f_{100}$:

In vaccine trials, traditionally a kernel density estimator (KDE) is used for
the purpose of density estimation (cf. Miladinovic et al., 2014). However,
using cross-validation, in Section 3, we show that the log-concave maximum
likelihood estimator (MLE) based estimators of Dümbgen and Rufibach (2009) and
Chen and Samworth (2013) minimize the estimated mean integrated squared error
(MISE) among a class of shape-constrained density estimators and KDEs.

### Shape-constrained tests of stochastic dominance:

The claim of a stochastic ordering between two samples is made stronger when
it is backed by a test of stochastic dominance. We say a distribution function
$F$ stochastically dominates another distribution function $G$ in first order
($F\succeq G$) if $F(x)\leq G(x)$ for all $x\in\mathbb{R}$. If $X$ and $Y$ are
two random variables with distribution functions $F$ and $G$, respectively,
then in this case, we say $X$ stochastically dominates $Y$, and write
$X\succeq Y$. The dominance is regarded as “strict” $(F\succ G$ or $X\succ Y$)
if, in addition, there exists $x\in\mathbb{R}$ such that $F(x)<G(x)$. If $F$
does not strictly stochastically dominate $G$, then this event is defined as
the non-dominance ($F\nsucceq G$) of $F$ over $G$ (Whang, 2019, p. 25).

The rejection of the null of non-dominance against the alternative of
stochastic dominance makes the strongest case for ranking one distribution
over the other (Davidson and Duclos, 2013; Álvarez-Esteban et al., 2016;
Ledwina and Wyłupek, 2012). However, the resulting test suffers from lack of
power because the overlap of distribution functions at the tails make
unrestricted stochastic dominance almost impossible to establish via
hypothesis testing (cf. Davidson and Duclos, 2013; Whang, 2019). There are
many ways to deal with this difficulty. For example, Ledwina and Wyłupek
(2012) exploits the fact that the dominance between two distribution functions
employs a dominance between their Fourier coefficients under a carefully
chosen basis. However, their method rejects the null of non-dominance in some
scenarios when the distribution functions cross each other. Another remedy
discussed in the literature (Kaur et al., 1994; Davidson and Duclos, 2013;
Álvarez-Esteban et al., 2016) involves excluding the tail region because it
does not contain enough reliable information for the problem at hand. This
type of tests focus only on a compact set $D$ inside the interior of the
combined support $\\{0<F+G<2\\}$. In Section 4, we follow the latter strategy
because it allows for many ways of incorporating shape-constraints with
desirable asymptotic properties.

To the best of our knowledge, we are the first to introduce the use of shape
constraints in the context of testing the null of non-dominance against
stochastic dominance. Moreover, one of our nonparametric test statistics in
Section 4.1 has not, to our knowledge, previously been studied in the context
of testing the null of non-dominance. In Section 4.2, we show that the shape-
constrained and nonparametric versions of our tests control the asymptotic
type I error at any null configuration under reasonable conditions. We also
show that our tests are asymptotically unbiased, and consistent against all
alternatives lying in the interior of the class of alternative distributions.
In Section 4.3, we empirically show that the shape-constrained tests have
better power than their nonparametric counterparts although they have the same
asymptotic critical values. Section 4.4 analyses the application of these
tests to our data. The proofs are deferred to Appendix A.

### Shape-constrained plug-in estimators of the squared Hellinger distance:

In Section 5, we construct plug-in estimators of the squared Hellinger
distance. It is well known that, unless bias-corrected, plug-in estimators
based on the KDE generally have a first-order bias (cf. Section 2 of Robins et
al., 2009). In contrast, for some smooth functionals, shape constrained MLE
based plug-in estimators do not require further bias correction (cf.
Jankowski, 2010; Mukherjee et al., 2019). The results in Lopuhaä and Musta
(2019) indicate that the squared Hellinger distance is an example of such a
smooth functional. When the underlying density is unimodal, we show that our
unimodal density based plug-in estimator enjoys the same asymptotic guarantees
as that of the bias corrected KDE based estimators (Kandasamy et al., 2015)
under some regularity conditions. The simulation studies in Section 5.3
suggest that similar results hold for our log-concave MLE based plug in
estimators as well. In fact, our smooth log-concave MLE based estimator shows
stable performance across all settings, where even the bias-corrected version
of the KDE based estimator struggles in some cases.

In the process, we develop theoretical tools for analyzing the asymptotic
behavior of plug-in estimators based on the unimodal density estimator of
Birgé, which may of independent interest. We defer the latter analysis to
Appendix B.  The methods developed in this paper are implementable using the R
package SDNNtests (Laha and Luedtke, 2020), which is available on GitHub.

###  1.2 Notations and terminologies

Before proceeding further, we introduce some notation that will be used
throughout this paper. We consider two independent samples
$X_{1},\ldots,X_{m}$ and $Y_{1},\ldots,Y_{n}$ drawn from distributions with
densities $f$ and $g$. We denote the corresponding distribution functions by
$F$ and $G$, respectively. The respective empirical distribution functions
will be denoted by $\mathbb{F}_{m}$ and $\mathbb{G}_{n}$. The pooled sample
$(X_{1},\ldots,X_{m},Y_{1},\ldots,Y_{n})$ has size $N=m+n$. We denote the
corresponding empirical distribution function by $\mathbb{H}_{N}$.

For $k\geq 1$, we let $||\cdot||_{k}$ denote the usual $L_{k}$ norm, i.e.
$||\mu||_{k}=\nolinebreak(\int_{\mathbb{R}}|\mu(x)|^{k}dx)^{1/k}$, where $\mu$
is a function supported on the real line. Also, we denote
$||\mu||_{\infty}=\sup_{x\in\mathbb{R}}\mu(x)$. For a density $f$, denote by
$\text{supp}(f)$ the set $\\{x\ :\ f(x)>0\\}$. For a concave function
$f:\mathbb{R}\mapsto\mathbb{R}$, the domain $\text{dom}(f)$ will be defined as
in (Rockafellar, 1970, p. 40), that is, $\text{dom}(f)=\\{x\in\mathbb{R}\ :\
f(x)>-\infty\\}$. For a sequence of measures $\\{P_{n}\\}_{n\geq 1}$, we say
$P_{n}$ converges weakly to $P$, and write $P_{n}\to_{d}P$, if
$\lim\limits_{n\to\infty}\int\mu dP_{n}=\int\mu dP$ holds for any bounded
continuous function $\mu:\mathbb{R}\mapsto\mathbb{R}$. For any two sets
$A,B\subset\mathbb{R}$, we denote by $\text{dist}(A,B)$ the quantity
$\min\\{|x-y|\ :\ x\in A,y\in B\\}$. We let $\text{int}(A)$ denote the
interior of the set $A$.

##  2 Background: HVTN 097 and HVTN 100

This section presents an exploratory analysis of the dataset. For both trials,
we consider the magnitude of IgG binding to the V1V2 region of seven clade C
glycoprotein $70$ antigens. The immune responses were measured by an HIV-1
binding antibody multiplex assay (BAMA). Following Haynes et al. (2012), Gray
et al. (2019), and Bekker et al. (2018), we use the log-transformed net median
fluorescence intensity (MFI) as the measure of immune response for statistical
analysis.

In HVTN 097 and HVTN 100, four injections of HIV vaccines were given at months
0, 1, 3, and 6. In this study, we only consider the responses measured two
weeks after the month six vaccination, which is considered to be the peak
immune response time point.  We include only those vaccinees in this study who
(a) completed the first four scheduled vaccinations and provided samples at
two weeks after the month six vaccination (known as vaccinated per-protocol
participants), and (b) developed a positive immune response for at least one
of the seven clade C V1V2 antigens. There are $68$ and $180$ vaccinees in the
HVTN 097 and HVTN 100 trial, respectively, who satisfy the above criteria. We
base our analysis on the aggregated response, averaged over the seven clade C
antigens mentioned above, and refer only to the latter when we say “immune
response”. We let $F_{097}$ and $F_{100}$ denote the distribution functions
corresponding to the aggregated response from the two trials.

i ECDF

ii Histogram

iii Boxplot

iv KDE

Figure 1:  Plots (i), (ii), (iii), (iv) display the empirical distribution
functions (ECDF), the histogram, the boxplot, and the KDEs of the average IgG
binding responses corresponding to the HVTN 097 and HVTN 100 regimens.

Figure 1 illustrates the empirical CDFs, the histogram, the boxplot, and the
KDES of the immune responses from the two trials. Figure 1i illustrates that
$F_{100}$ is always greater than $F_{097}$ except at the tails, hinting at the
stochastic dominance of the HVTN 097 immune response over the HVTN 100 immune
response. The histogram in Figure 1ii, the boxplot in Figure 1iii, and the
plot of the KDEs in Figure 1iv also suggest that the HVTN 097 trial induces
higher immune response. Bekker et al. (2018) also indicated that the magnitude
of positive responses in the HVTN 100 trial is lower than that of the RV144
trial, which, on the other hand, is reported to be slightly lower than the
responses in the HVTN 097 trial (Gray et al., 2019). Therefore, it makes sense
to posit the null of non-dominance of $F_{097}$ over $F_{100}$ against the
alternative that $F_{097}$ stochastically dominates $F_{100}$. Note also that
the comparisons between the two sets of immune responses is reasonable because
the trials were conducted on similar populations and share approximately the
same support (cf. Figure 1). See also Table 6 in Appendix D for a comparison
between the two trials.

##  3 Density estimation

This section compares different estimators of $f_{097}$ and $f_{100}$. Since
our study includes some unimodal and log-concave density estimators, we begin
by presenting some observations in support of the shape-restriction
assumptions.
As mentioned previously, the unimodality assumption is not unreasonable owing
to the homogeneity of the trial populations. The histogram in Figure 1ii and
the KDEs displayed in Figure 1iv both support this claim. Although unimodality
is a naturally occurring shape constraint, the class of all unimodal densities
is too large to admit an MLE (Birgé, 1997). The class of log-concave densities
is a subclass of the class of unimodal densities, which is small enough to
admit an MLE (Pal et al., 2007), but contains most of the commonly used
subexponetial unimodal densities (Walther et al., 2009). The log-concave MLE
can be computed efficiently using the R package logcondens. Furthermore, there
is a smoothed version of the log-concave MLE, which is also free of tuning
parameters (Chen and Samworth, 2013), and thus can potentially replace the
smoothed unimodal density estimators, which generally depend crucially on
external tuning parameters. In view of the above, many researchers, e.g.
Walther (2002) and Walther et al. (2009), advocate opting for log-concavity
shape constraint in situations where unimodality may seem plausible.

Although it is difficult to provide visual evidence in favor of the assumption
of log-concavity, the KDE plot in Figure 1iv does not indicate a departure
from log-concavity either. Using the test of log-concavity in Chen and
Samworth, we test the null of log-concavity against the alternative of
violation of log-concavity for $f_{097}$ and $f_{100}$. The corresponding
p-values for $f_{097}$ and $f_{100}$ are $0.4890$ and $0.4631$, respectively,
which implies that our data does not have enough evidence for rejecting the
null of log-concavity.

Figure 2: Birgé’s estimator (left), log-concave MLE (middle), and the smoothed
log-concave MLE (right) based on the immune responses in the HVTN 097 and the
HVTN 100 trial.

We briefly describe below the density estimators that we consider in this
section.
Grenander-type unimodal density estimator (Birgé, 1997): Although the class of
all unimodal densities does not permit an MLE (Birgé, 1997), when the mode of
a unimodal density $f$ is known, the MLE $\widehat{f}^{0}_{m}$ exists (Rao,
1969), and it is a piecewise constant unimodal density with mode at the true
mode. However, the MLE $\widehat{f}^{0}_{m}$ is generally not useful due to
the lack of knowledge on the location of the true mode. The estimator
$\widehat{f}_{m}$ presented in Birgé (1997) is a piecewise constant estimator
of $f$, constructed in such a way so that the Kolmogorov-Smirnov distance
between the corresponding distribution functions can be made arbitrarily
small, in particular, smaller than a pre-fixed number $\eta>0$. Although this
$\eta>0$ is an external parameter, unlike the kernel bandwidth, a smaller
$\eta$ always leads to a more accurate estimation (Birgé, 1997), and hence it
does not require actual tuning. Indeed, if we choose the parameter
$\eta=o(m^{-1})$, then our Lemma A.5 in Appendix A ensures that the total
variation distance between $\widehat{f}_{m}$ and $\widehat{f}^{0}_{m}$ is
$o(m^{-1/2})$ with probability one, where $\widehat{f}^{0}_{m}$ was defined to
be the MLE of $f$ had the true mode been known. Therefore, we choose $\eta$ to
be the inverse of the combined sample size of the two trials.

Smooth unimodal estimator: There is a substantial body of literature on smooth
unimodal density estimators. See for instance, Eggermont and LaRiccia (2000),
Mammen et al. (2001), Hall and Huang (2002), Wolters (2012), Meyer (2012),
Turnbull and Ghosh (2014), and Wolters and Braun (2018), among others. The
smooth unimodal density estimators generally depend crucially on external
tuning parameters, as mentioned previously. Since we opt for shape constraints
mainly to avoid tuning parameters and we already have at our disposal the
tuning free smooth log-concave MLE estimator, the smooth unimodal density
estimators are not particularly attractive to us. Still, since the current
section solely focuses on density estimation, we include some smooth unimodal
estimators for comparison here, namely the estimators of Turnbull and Ghosh
(2014), Hall and Huang (2002), Wolters (2012), and Wolters and Braun (2018).

Turnbull and Ghosh (2014) approximates the unknown unimodal density using
Bernstein polynomial. We use the condition number approach of Turnbull and
Ghosh (2014) to select tuning parameters related to the degree of the
polynomial.  The estimators in Hall and Huang (2002), Wolters (2012), and
Wolters and Braun (2018) are kernel based. To compute them, we use the R
package scdensity with the default choice of bandwidth.

Log-concave estimators: Our first log-concave estimator is the MLE among the
class of all log-concave densities. The MLE is continuous but non-smooth
(Dümbgen and Rufibach, 2009). The second estimator is a smoothed version of
the MLE (Dümbgen and Rufibach, 2009; Chen and Samworth, 2013). The smoothing
parameter for the latter is data dependent and has a closed form formula, and
hence it does not require external tuning. Figure 2 displays these two density
estimators. For more on the properties of the log-concave density estimators,
see, e.g., Balabdaoui et al. (2009), Dümbgen et al. (2011), Cule and Samworth
(2010), and Doss and Wellner (2016).
Kernel density estimators: We consider kernel density estimators (KDE) with
Gaussian kernel. The optimal bandwidth was chosen either by the univariate
plug-in selector of Wand and Jones (1994), or the univariate least square
cross-validation (LSCV) selector of Bowman (1984) and Rudemo (1982) (see
Figure 1iv).

We prefer the estimator $\widehat{f}_{m}$ of a density $f$ with the smallest
mean integrated squared error (MISE), which is given by
$E\int_{\mathbb{R}}(\widehat{f}_{m}(x)-f(x))^{2}dx$. Noting minimizing the
MISE with respect to $\widehat{f}_{n}$ is equivalent to minimizing

$\text{MISE-
err}=E\int_{\mathbb{R}}\widehat{f}_{m}^{2}(x)dx-2E\int_{\mathbb{R}}\widehat{f}_{m}(x)f(x)dx,$
(3.1)

we estimate the latter quantity using a ten folds cross-validation. We also
estimate the negative log-likelihood

$-l_{m}=-m^{-1}\sum_{i=1}^{n}\log\widehat{f}_{m}(X_{i}).$ (3.2)

Table 1: Table of the estimated risks for different density estimators of the aggregated immune response in trials HVTN 097 and HVTN 100 | HVTN 097 |  | HVTN 100 |
---|---|---|---|---
Estimators | MISE-err | $-l_{m}$ | MISE-err | $-l_{m}$
Unimodal (Birgé’s estimator) | -0.172 | 1.662 | -0.126 | 2.071
Unimodal (Bernstein) | -0.191 | 1.797 | -0.130 | 2.133
Unimodal (Wolters, 2012) | -0.190 | 1.878 | -0.102 | 2.549
Unimodal (Wolters and Braun, 2018) | -0.190 | 1.795 | -0.124 | 2.300
Unimodal (Hall and Huang, 2002) | -0.190 | 1.965 | -0.124 | 2.288
Smooth log-concave MLE | -0.196 | 1.758 | -0.129 | 2.127
Log-concave MLE | -0.193 | 1.758 | -0.130 | 2.127
KDE (plug-in bandwidth selector) | -0.189 | 1.750 | -0.128 | 2.140
KDE (LSCV bandwidth selector) | -0.189 | 1.777 | -0.128 | 2.140

Based on these risks (see Table 1), our recommended estimators are the log-
concave estimators which exhibit the lowest risk in an overall sense. Table 1
indicates that Birgé’s estimator excels in minimizing the $-l_{m}$ risk but it
has higher estimated MISE when compared to the other estimators. This can be
attributed to its spikes at the mode (see Figure 2), which contributes large
positive terms to $l_{n}$ and MISE. Grenander type unimodal estimators are
known to exhibit such “spike-problem” at the mode (Walther et al., 2009),
which is caused by the inconsistency of the density estimator at the mode (for
more details, see Woodroofe and Sun, 1993; Balabdaoui et al., 2009).

##  4 Test of stochastic dominance

To provide an answer to Q2, we construct tests for the null of non-dominance
against that of stochastic dominance using the log-concave MLEs and the
unimodal estimator of Birgé. We compare the resulting shape-constrained tests
with their nonparametric counterparts.

Our shape-restricted methods rely on estimating the densities $f$ and $g$. We
denote the corresponding unimodal estimators of Birgé by $\widehat{f}_{m}$ and
$\widehat{g}_{n}$, respectively. The construction of Birgé’s estimators
requires a tuning parameter $\eta$, which we set to be $N^{-1}$ where $N=m+n$.
We let $\tilde{f}_{m}$ and $\tilde{g}_{n}$ denote the log-concave MLEs of $f$
and $g$ (Dümbgen and Rufibach, 2009), and write $\tilde{f}^{sm}_{m}$ and
$\tilde{g}^{sm}_{n}$ for their respective smoothed versions (Chen and
Samworth, 2013). The corresponding distribution functions will be denoted by
$\widehat{F}_{m}$, $\widehat{G}_{n}$, $\tilde{F}_{m}$, $\tilde{G}_{n}$,
$\tilde{F}^{sm}_{m}$, and $\tilde{G}^{sm}_{n}$, respectively.

As $m$ and $n$ approach $\infty$, we assume that $m/N\to\lambda\in(0,1)$.

Letting $H=\lambda F+(1-\lambda)G$, for $p\in(0,1/2)$, we also define the sets

$D_{p}(F,G):=D_{p}=[H^{-1}(p),H^{-1}(1-p)]\quad\text{ and }\quad
D_{p,m,n}:=[{\mathcal{H}}^{-1}(p),{\mathcal{H}}^{-1}(1-p)].$ (4.1)

###  4.1 Construction of the tests

Suppose $D\subset\text{supp}(f)\cup\text{supp}(g)$ is compact. Following Kaur
et al. (1994) and Davidson and Duclos, we formulate the hypotheses as follows:

${\mathcal{H}}_{0}:F(z)\geq G(z)\text{ for some }z\in
D\quad\text{vs.}\quad{\mathcal{H}}_{1}:F(z)<G(z)\text{ for all }z\in D.$ (4.2)

The null configuration $\mathcal{H}_{0}$ occurs if $F=G$, or $G$
stochastically dominates $F$, or if $F$ and $G$ touch or cross each other on
$D$. Thus our formulation is unable to reject the null when $F$ and $G$ touch
at a point in $D$ even if $F\succ G$. This limitation seems to be unavoidable
because such a configuration $(F,G)$ lies on the common boundary shared by
$\\{(F,G):F\succeq G\\}$ and $\\{(F,G):F\npreceq G\\}$, and hence can not be
discriminated from the null without sacrificing control over the size of test.
Notably, our exploratory analysis (see Fig 1i) suggests that it is unlikely
that $F_{097}$ and $F_{100}$ fall in this category. See Figure 3 for examples
of different scenarios associated with our hypotheses.

Regarding the choice of $D$, we need to ensure that $D$ is inside the combined
support of $f$ and $g$ because otherwise, ${\inf_{z\in D}[G(z)-F(z)]}$ will
always be $0$. The set $D_{p}$ defined in (4.1) satisfies this criterion. In
practice, we replace this unknown $D_{p}$ by $D_{p,m,n}$ defined in (4.1),
which always utilizes $100(1-2p)\%$ of the combined data. Naturally, if $p$ is
too small, rejection of the null will be difficult, where a large $p$ will
exclude a large portion of the data, which might be unnecessary. If only some
particular interval of the data is of practical interest (e.g. some particular
range of immune responses or biomarkers), we suggest setting $D$ to be the
smallest superset of that interval. In the absence of such prior knowledge, we
suggest choosing the largest $p$ so that $D_{p,m,n}$ excludes the tail region
where empirical distribution functions overlap. We will return to this issue
later in Section 4.4, with a demonstration on our motivating dataset.

Now we are in a position to introduce our test statistics.

Minimum t-statistic: This statistic was first introduced by Kaur et al. (1994)
in context of second order stochastic dominance, and then extended to the
first order by Davidson and Duclos. For distribution functions $F_{1}$ and
$F_{2}$, this statistic is given by

$T_{m,n}^{\text{min}}(F_{1},F_{2})=\inf_{x\in
D_{p,m,n}}\dfrac{\Big{(}F_{2}(x)-F_{1}(x)\Big{)}}{\sqrt{\dfrac{F_{1}(x)\Big{(}1-F_{1}(x)\Big{)}}{m}+\dfrac{F_{2}(x)\Big{(}1-F_{2}(x)\Big{)}}{n}}}.$
(4.3)

Our tests reject the ${\mathcal{H}}_{0}$ for large values of
$T_{m,n}^{\text{min}}(\widehat{F}_{m},\widehat{G}_{n})$,
$T_{m,n}^{\text{min}}(\tilde{F}_{m},\tilde{G}_{n})$, and
$T_{m,n}^{\text{min}}(\mathbb{F}_{m},\mathbb{G}_{n})$.  The last test-
statistic, which is nonparametric, equals the minimum t-statistic of the Kaur
et al. in context of first order stochastic dominance.

Two sample empirical process (TSEP) type test statistic:Our second test
rejects the ${\mathcal{H}}_{0}$ for large values of
$T_{m,n}^{\text{tsep}}(\widehat{F}_{m},\widehat{G}_{n})$,
$T_{m,n}^{\text{tsep}}(\tilde{F}_{m},\tilde{G}_{n})$, or
$T_{m,n}^{\text{tsep}}(\mathbb{F}_{m},\mathbb{G}_{n})$, where for distribution
functions $F_{1}$ and $F_{2}$, $T_{m,n}^{\text{tsep}}$ is defined by

$T_{m,n}^{\text{tsep}}(F_{1},F_{2})=\sqrt{\dfrac{mn}{N}}\inf_{z\in[p,1-p]}\dfrac{F_{2}({\mathcal{H}}^{-1}(z))-F_{1}({\mathcal{H}}^{-1}(z))}{\sqrt{z(1-z)}}.$
(4.4)

Ledwina and Wyłupek (2013) uses a test statistic similar to
$T_{m,n}^{\text{tsep}}(\mathbb{F}_{m},\mathbb{G}_{n})$ (the second test
statistic in Section 2.2 of Ledwina and Wyłupek, 2013) for testing the null of
stochastic dominance against non-dominance. The pivotal distribution of their
test statistic is completely different from ours because they based their
critical values on the configuration $F=G$, which is very different from what
we will consider. We are not aware of any existing test which uses
$T_{m,n}^{\text{tsep}}(\mathbb{F}_{m},\mathbb{G}_{n})$ for testing non-
dominance against stochastic dominance.

Wilcoxon rank sum (WRS) type test statistic: Wilcoxon rank sum (WRS) test is
widely used for comparing two vaccines (cf. Miladinovic et al., 2014) although
the WRS test is actually designed for testing location shift. It is a popular
choice for testing the null $F=G$ against the alternative $F=G(\cdot-\delta)$
for $\delta>0$ (Lee and Wolfe, 1976). The WRS test is the most powerful
nonparametric test for testing the following hypotheses (cf. Example 25.46 of
Van der Vaart, 1998):

${\mathcal{H}}_{0}^{a}:\int_{\mathbb{R}}G(z)dF(z)\geq 1/2\quad
vs\quad{\mathcal{H}}_{1}^{a}:\int_{\mathbb{R}}G(z)dF(z)<1/2.$ (4.5)

Although the WRS test is not designed to test the null of non-dominance, we
include this test to demonstrate its failure to control the type I error at
some null configurations.

The one-sided WRS test rejects ${\mathcal{H}}_{0}$ for large values of
$T_{m,n}^{\text{wrs}}(\mathbb{F}_{m},\mathbb{G}_{n})$, where, for distribution
functions $F_{1}$ and $F_{2}$,

$T_{m,n}^{\text{wrs}}(F_{1},F_{2})=\sqrt{\dfrac{12mn}{N+1}}\bigg{(}\int_{\mathbb{R}}F_{2}(x)dF_{1}(x)-1/2\bigg{)}.$
(4.6)

$T_{m,n}^{\text{wrs}}(\mathbb{F}_{m},\mathbb{G}_{n})$ is the Mann-Whitney form
of the two-sample WRS statistic. The corresponding shape-constraint versions
are given by $T_{m,n}^{\text{wrs}}(\widehat{F}_{m},\widehat{G}_{n})$ and
$T_{m,n}^{\text{wrs}}(\tilde{F}_{m},\tilde{G}_{n})$.

We excluded tests based on the smoothed log-concave MLE because rigorous
asymptotic analysis of the corresponding tests is out of the scope of the
present paper. However, our empirical study in Section 4.3 includes minimum
t-test and TSEP test based on the smoothed log-concave MLE. Our simulations
indicate that the asymptotic critical values of the tests based on
$T_{m,n}^{\text{min}}(\tilde{F}_{m},\tilde{G}_{n})$ and
$T_{m,n}^{\text{tsep}}(\tilde{F}_{m},\tilde{G}_{n})$ are valid for the
corresponding smoothed log-concave tests. Our simulations also indicate that
the finite sample performance of the tests based on the log-concave MLE and
the smoothed log-concave MLE are quite similar. We leave the rigorous analysis
of the tests based on the smoothed log-concave MLE for future study.

######  Remark 1.

Since the nonparametric tests use the empirical distribution function, shape
constrained methods do not gain any advantage in terms of tuning parameters.
Also, we will see that the nonparametric and shape constrained tests are
asymptotically equivalent. For moderate sized samples, however, our
simulations in Section 4.3 show that the shape-constrained tests exhibit
better performance.

######  Remark 2.

Davidson and Duclos proposed an empirical likelihood ratio approach to test
${\mathcal{H}}_{0}$ vs ${\mathcal{H}}_{1}$. We did not appeal to this
nonparametric approach in this paper because this approach does not extend
easily to shape-constrained scenarios. Regardless, we point out that Davidson
and Duclos showed that their empirical likelihood ratio test is asymptotically
equivalent to the minimum t-test.

In the sequel, we may use the terms “log-concave” or “unimodal” to refer to
the test statistics based on the log-concave or unimodal density estimators.
For example, we may refer to
$T_{m,n}^{\text{min}}(\widehat{F}_{m},\widehat{G}_{n})$ and
$T_{m,n}^{\text{min}}(\tilde{F}_{m},\tilde{G}_{n})$ as the unimodal minimum
t-statistic and the log-concave minimum t-statistic, respectively. Also,
unless otherwise specified, the terms “null” and “alternative” will refer to
the ${\mathcal{H}}_{0}$ and ${\mathcal{H}}_{1}$ defined in (4.2),
respectively.

Figure 3:  This figure displays plots of two distribution functions $F$ (red)
and $G$ (blue). The range of $x$ in these plots correspond to
$D_{p}:=D_{p}(F,G)$. (i) $F=G$ on $D_{p}$. This is a null configuration. (ii)
$F$ and $G$ cross each other on $D_{p}$. This is also a null configuration.
(iii) $F$ strictly stochastically dominates $G$. In fact $G(x)>F(x)$ for all
$x\in D_{p}$. This is an alternative configuration. (iv) $F$ and $G$ touch
each other at the endpoint of $D_{p}$. This is a null configuration.

###  4.2 Asymptotic distribution

In this section, we explore the asymptotic distribution of our test
statistics. We show that the minimum t-test and the two sample empirical
process (TSEP) test asymptotically control type I error for each null
configuration and they are asymptotically consistent against each
$(F,G)\in{\mathcal{H}}_{1}$. We also show that, with the exception of the test
based on log-concave MLE, the WRS type tests control the type I error for
distributions in ${\mathcal{H}}_{0}^{a}$ and are consistent against
$(F,G)\in{\mathcal{H}}_{1}^{a}$. We first prove the asymptotic results on the
nonparametric test statistics. Then, we show that the shape-constrained test
statistics are equivalent to their nonparametric counterparts up to a
$o_{p}(N^{-1/2})$ term, which implies that the same critical values can be
used for them.

We exclude the log-concave MLE based WRS statistic from our discussion because
we are unable to infer on its asymptotic limit. The difficulty arises due to
our inability to track the asymptotic behaviour of
$\sqrt{m}||\tilde{F}_{m}-\mathbb{F}_{m}||_{\infty}$. In the remainder of this
section, by “shape-constrained test statistics”, we will therefore refer to
$T_{m,n}^{\text{min}}(\widehat{F}_{m},\widehat{G}_{n})$,
$T_{m,n}^{\text{tsep}}(\widehat{F}_{m},\widehat{G}_{n})$,
$T_{m,n}^{\text{wrs}}(\widehat{F}_{m},\widehat{G}_{n})$,
$T_{m,n}^{\text{min}}(\tilde{F}_{m},\tilde{G}_{n})$, and
$T_{m,n}^{\text{tsep}}(\tilde{F}_{m},\tilde{G}_{n})$ only. For the sake of
clarity, in Table 2, we summarize the current state of results on the
different tests discussed in this paper.

Test | Asymptotic results | Empirical results
---|---|---
Nonparametric
Minimum t-test |  previously known (Kaur et al., 1994) | included in this paper
TSEP | we derived | $---^{\prime\prime}---$
WRS |  previously known (Dwass, 1956) | $---^{\prime\prime}---$
Unimodal
Minimum t-test | we derived | included in this paper
TSEP | $---^{\prime\prime}---$ | $---^{\prime\prime}---$
WRS | $---^{\prime\prime}---$ | $---^{\prime\prime}---$
Log-concave
Minimum t-test | we derived | included in this paper
TSEP | $---^{\prime\prime}---$ | $---^{\prime\prime}---$
WRS | unknown | $---^{\prime\prime}---$
Smoothed log-concave
Minimum t-test | unknown | included in this paper
TSEP | $---^{\prime\prime}---$ | $---^{\prime\prime}---$
WRS | $---^{\prime\prime}---$ | not included
Table 2: Table summarizing different tests used in this paper and current
results on them.

Before going into further details, we state a technical condition that will be
required by all of our Theorems.

######  Condition N.

$F$ and $G$ are continuous. Also, $F$ and $G$ have densities $f$ and $g$,
respectively, such that $\text{supp}(f)\cup\text{supp}(g)$ contains an open
neighborhood of $D_{p}$.

The first requirement of Condition Condition N, i.e., the continuity of $F$
and $G$, is necessary for the weak convergence of the empirical processes to
Brownian bridges. The second requirement ensures that
$\text{dist}(D_{p,m,n},D_{p})$ approaches $0$ as $m,n\to\infty$ with
probability one. These assumptions are likely to be satisfied by our immune
response data provided $p$ is not too small (see Figure 1i and Figure 1ii).
For the rest of the paper, we restrict our attention to $F$ and $G$ that
satisfy Condition Condition N.

####  4.2.1 Asymptotic critical values of the nonparametric tests

We begin our discussion with the minimum t-statistic and the TSEP statistic.
Our first objective is to identify the null configurations that lead to the
highest asymptotic type I error. Here we remind the readers that $(F,G)$ is a
null configuration if there exists $x\in D_{p}$ so that $G(x)\leq F(x)$. One
may guess that the interesting cases appear on the boundary of
${\mathcal{H}}_{0}$. However, to formally discuss the boundary of
${\mathcal{H}}_{0}$, we need to equip it with a suitable topology. To
formalize our discussion, we consider the space $\mathcal{F}$ of all
continuous distribution functions on $\mathbb{R}$, and equip it with the
uniform metric ${d(F,F^{\prime})=\sup_{x\in\mathbb{R}}|F(x)-F^{\prime}(x)|.}$
Consider the product space $\mathcal{F}\times\mathcal{F}$ with the metric

$d_{2}\Big{(}(F,G),(F^{\prime},G^{\prime})\Big{)}=\max\bigg{\\{}d(F,F^{\prime}),d(G,G^{\prime})\bigg{\\}}.$

By an abuse of notation, we denote by ${\mathcal{H}}_{0}$ and
${\mathcal{H}}_{1}$ the set of all combinations
$(F,G)\in\mathcal{F}\times\mathcal{F}$ that satisfy the hypotheses
${\mathcal{H}}_{0}$ and ${\mathcal{H}}_{1}$, respectively. For $i=0,1$, we
denote the closure of ${\mathcal{H}}_{i}$ in $\mathcal{F}\times\mathcal{F}$ by
$\text{cl}({\mathcal{H}}_{i})$. Then the boundary of ${\mathcal{H}}_{i}$ is
given by $\text{cl}({\mathcal{H}}_{i})\setminus\text{int}({\mathcal{H}}_{i})$.
The following lemma characterizes $\text{bd}({\mathcal{H}}_{0})$ and
$\text{int}({\mathcal{H}}_{0})$.

######  Lemma 1.

${\mathcal{H}}_{0}$ is a closed subset of $\mathcal{F}\times\mathcal{F}$ with
boundary

$\text{bd}({\mathcal{H}}_{0})=\bigg{\\{}(F,G)\in\mathcal{F}\times\mathcal{F}\
:\ \sup_{x\in D_{p}}\Big{(}F(x)-G(x)\Big{)}=0\bigg{\\}}.$

Moreover, $\text{bd}({\mathcal{H}}_{0})=\text{bd}({\mathcal{H}}_{1})$. Also,
the interior of ${\mathcal{H}}_{0}$ is given by

$\text{int}({\mathcal{H}}_{0})=\bigg{\\{}(F,G)\in\mathcal{F}\times\mathcal{F}\
:\ F(z)>G(z)\text{ for some }z\in D_{p}\bigg{\\}}.$

Figure 3(ii) gives an example of an $(F,G)$ pair in the interior of
${\mathcal{H}}_{0}$. The following lemma entails that the minimum t-statistic
and the TSEP statistic are asymptotically degenerate on
$\text{int}({\mathcal{H}}_{0})$.

######  Lemma 2.

Suppose $(F,G)\in\text{int}({\mathcal{H}}_{0})$ satisfies Condition Condition
N, $m/N\to\lambda$, and $\text{supp}(f)\cup\text{supp}(g)$ contains an open
neighborhood of $D_{p}$. Then,

$T_{m,n}^{\text{min}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{p}-\infty\quad\text{and}\quad
T_{m,n}^{\text{tsep}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{p}-\infty.$

The proof of Lemma 2 for the minimum t-statistic can be found in Whang (2019)
(see also Davidson and Duclos, 2013). However, we include it in Appendix 2 for
the sake of completeness. Lemma 2 indicates that non-trivial type I errors can
originate only at the boundary of ${\mathcal{H}}_{0}$, which is a subset of
${\mathcal{H}}_{0}$ because the latter is a closed set (see Lemma 1). Lemma 1
also implies that $\text{bd}({\mathcal{H}}_{0})$ consists of all those $F$ and
$G$ that touch each other on $D_{p}$. To concretize this idea, we define the
contact set $C_{p}$ by

$C_{p}=\\{x\in D_{p}\ :\ F(x)=G(x)\\}.$

Note that if $(F,G)\in\text{bd}({\mathcal{H}}_{0})$, then
$C_{p}\neq\emptyset$, where $F=G$ on $C_{p}$, and $F<G$ on $D_{p}\setminus
C_{p}$. Let us also define

$H(C_{p})=\\{t\in[p,1-p]\ |\ H^{-1}(t)\in C_{p}\\}.$ (4.7)

Theorem 1 shows that the asymptotic distribution of the test statistics on
$\text{bd}({\mathcal{H}}_{0})$ crucially depends on this contact set $C_{p}$
and $H(C_{p})$. The proof of Theorem 1 is given in Appendix A.0.2.

######  Theorem 1.

Suppose $(F,G)\in\text{bd}({\mathcal{H}}_{0})$, $m/N\to\lambda$, and $F$ and
$G$ have continuous densities $f$ and $g$ satisfying

$\inf_{x\in D_{p}}\min\\{f(x),g(x)\\}>0,$

where $D_{p}=[H^{-1}(p),H^{-1}(1-p)]$. Let $\mathbb{U}$ denote a standard
Brownian bridge. Then under the stated conditions, the following assertions
hold:

  * A.

$\displaystyle
T_{m,n}^{\text{min}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{d}\inf_{x\in
C_{p}}\dfrac{\mathbb{U}\circ F(x)}{\sqrt{F(x)(1-F(x))}},$ (4.8)

where $C_{p}=\\{x\in D_{p}\ :\ F(x)=G(x)\\}$.

  * B.

$T_{m,n}^{\text{tsep}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{d}\inf_{t\in
H(C_{p})}\sqrt{\dfrac{\lambda}{1-\lambda}}\dfrac{\mathbb{L}_{0}(t)}{\sqrt{t(1-t)}},$

where $\mathbb{L}_{0}$ is the centred Gaussian process given by (A.17) of
Appendix A, and $H(C_{p})$ is as in (4.7).

  * C.

In particular, if $C_{p}\subset\text{int}(D_{p})$, then

$T_{m,n}^{\text{tsep}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{d}\inf_{t\in
H(C_{p})}\dfrac{\mathbb{U}(t)}{\sqrt{t(1-t)}}.$

The Gaussian process $\mathbb{L}_{0}$, which is defined in (A.17) of Appendix
A, depends on $F$ and $G$. We postpone further discussion on the form of
$\mathbb{L}_{0}$ till Appendix A. Next we discuss the implication of Theorem 1
on the minimum t-test. Then we will discuss the case of the TSEP test.

Asymptotic critical value of minimum t-test: Theorem 1 reveals an interesting
fact: the length of $C_{p}$ imposes a stochastic ordering among the limiting
laws of the minimum t-statistic for the boundary configurations. To elaborate,
let us consider $(F_{1},G_{1})$ and
$(F_{2},G_{2})\in\text{bd}({\mathcal{H}}_{0})$, with respective contact sets
$C_{p}^{1}$ and $C_{p}^{2}$. Then, under the conditions of Theorem 1,

$T_{m,n}^{\text{min}}(\mathbb{F}_{m},\mathbb{G}_{n})\xrightarrow{F_{i},G_{i}}_{d}\inf_{x\in
C^{i}_{p}}\dfrac{\mathbb{U}\circ F(x)}{\sqrt{F(x)(1-F(x))}},\quad\text{ for
}i=1,2.$

If the contact sets satisfy the ordering $C_{p}^{1}\subset C_{p}^{2}$, then

$\inf_{z\in C^{2}_{p}}\dfrac{\mathbb{U}\circ
F(z)}{\sqrt{F(z)(1-F(z))}}\preceq\inf_{z\in C^{1}_{p}}\dfrac{\mathbb{U}\circ
F(z)}{\sqrt{F(z)(1-F(z))}},$

implying that the limiting law of the minimum t-statistic under
$(F_{1},G_{1})$ stochastically dominates that under $(F_{2},G_{2})$. The
extreme cases for $C_{p}^{1}$ are the singleton sets $\\{x\\}$, where $x\in
D_{p}$. In this case, the asymptotic distribution of both test statistics is
standard Gaussian. Therefore, we set the critical value of our tests to be
$z_{\alpha}$, the $(1-\alpha)$-th quantile of the standard Gaussian
distribution.  The class of boundary configurations with a singleton contact
set is referred to as “the least favorable class” (LFC) (Davidson and Duclos,
2013). Figure 4 illustrates the difference between an LFC and an ordinary non-
LFC boundary combination. Theorem 2 of Davidson and Duclos shows that there is
no null configuration under which the law of the minimum t-statistic strictly
stochastically dominates that of the LFC configuration. This result, which
holds for any $m$ and $n$, implies that, among the null configurations, the
LFC configurations lead to the greatest dominance of $F$ over $G$. The above
finding, in conjunction with our Theorem 1, imply that our tests, whose
critical values are based on the LFC class, is likely to have asymptotic size
$\alpha$. This being a stronger assertion than the asymptotic control of type
I error can be an interesting topic for further investigation.

i $(F,G)\in\text{bd}({\mathcal{H}}_{0})$ but not LFC

ii $(F,G)$ is LFC

Figure 4: Illustration of $F$ (red) and $G$ (blue) on
$\text{bd}({\mathcal{H}}_{0})$. (i) The contact set $C_{p}$ is an interval
(ii) The contact set is singleton; so this is an LFC configuration.

######  Remark 3.

The asymptotic behavior of the nonparametric minimum t-statistic has been
previously studied (Davidson and Duclos, 2013; Kaur et al., 1994). However,
previous studies focus only on the asymptotic behaviour of the minimum
t-statistic at the LFC configurations and the interior of ${\mathcal{H}}_{0}$,
whereas our results show that there are other classes of boundary
configurations with non-vanishing type I error. Although existing results are
enough for the purpose of constructing critical values, our new results
provide a more complete understanding of the scenario.

Asymptotic critical value of the TSEP test: The asymptotic distribution of the
TSEP statistic under ${\mathcal{H}}_{0}$ also exhibits a monotonocity property
similar to the minimum t-test. To elaborate, suppose the pairs $(F_{1},G_{1})$
and $(F_{2},G_{2})$ have respective contact sets $C_{p}^{1}$ and $C_{p}^{2}$
satisfying $C_{p}^{1}\subset C_{p}^{2}$. Then

$\inf_{t\in
H(C_{p}^{2})}\frac{\mathbb{L}_{0}(t)}{\sqrt{t(1-t)}}\preceq\inf_{t\in
H(C_{p}^{1})}\frac{\mathbb{L}_{0}(t)}{\sqrt{t(1-t)}}.$

Therefore, similar to the case of the minimum t-test, the configurations with
singleton $H(C_{p})$ constitute the class of LFC configurations for the TSEP
test. Under Condition Condition N, $H$ is continuous and strictly increasing
on $C_{p}$. Therefore, $H(C_{p})$ is singletone if and only if $C_{p}$ is
singletone, i.e., $C_{p}$ is of the form $\\{x\\}$, where $x\in D_{p}$. To
find the critical value, it suffices to study the asymptotics in such LFC
cases. If $x\in\text{int}(D_{p})$, the TSEP test statistic weakly converges to
a standard Gaussian distribution by part C of Theorem 1. Till this point,
there has been no difference between the TSEP and the minimum t-test
statistic. If, however, $x\in\text{bd}(D_{p})$ (see Figure 15iv), the
asymptotic distribution of the TSEP test statistic can be different. To this
end, first we state a lemma, and then using this lemma, we explain the
asymptotics of the TSEP test when $C_{p}=\\{x\\}\subset\\{p,1-p\\}$.

######  Lemma 3.

Suppose $F$ and $G$ are as in Theorem 1 and $H(C_{p})=\\{t\\}$ for some
$t\in[p,1-p]$. Then

$\sqrt{\frac{\lambda}{1-\lambda}}\dfrac{\mathbb{L}_{0}(t)}{\sqrt{t(1-t)}}\sim
N(0,\sigma_{TSEP}^{2}),$

where

$\sigma_{TSEP}^{2}=\frac{\lambda(f\circ H^{-1}(t))^{2}+(1-\lambda)(g\circ
H^{-1}(t))^{2}}{\Big{(}\lambda f\circ H^{-1}(t)+(1-\lambda)g\circ
H^{-1}(t)\Big{)}^{2}}.$ (4.9)

Moreover, the following assertions also hold:

  * A.

$\sigma_{TSEP}^{2}=1$ if and only if $f\circ H^{-1}(t)=g\circ H^{-1}(t)$.
Otherwise,

$1<\sigma_{TSEP}^{2}\leq\max\\{(1-\lambda)^{-1},\lambda^{-1}\\}.$ (4.10)

  * B.

Given any $\epsilon\in(0,2(1-\lambda)/\lambda)$, we can find a constant
$C_{\lambda}>0$, depending only on $\lambda>0$, so that whenever
$g(H^{-1}(t))/f(H^{-1}(t))<C_{\lambda}$, then
$\sigma^{2}_{TSEP}>\lambda^{-1}-\epsilon.$

  * C.

Given any $\epsilon\in(0,2\lambda/(1-\lambda))$, we can find a constant
$C^{\prime}_{\lambda}>0$, depending only on $\lambda>0$, so that whenever
$f(H^{-1}(t))/g(H^{-1}(t))<C^{\prime}_{\lambda}$, then
$\sigma^{2}_{TSEP}>(1-\lambda)^{-1}-\epsilon.$

Lemma 3 has some interesting consequences. First, if $C_{p}=\\{H^{-1}(p)\\}$
or $\\{H^{-1}(1-p)\\}$, then using Theorem 1B and Lemma 3 one can show that
$T_{m,n}^{\text{tsep}}$ converges weakly to a centred Gaussian distribution
with variance $\sigma^{2}_{TSEP}$. Part A of Lemma 3 implies that if $F$ and
$G$ touch at $C_{p}$, i.e. if $f=g$ at the point of contact, then
$\sigma^{2}_{TSEP}$ is still one. Hence, the TSEP test statistic is
asymptotically standard Gaussian for this case. However, if $F$ and $G$ cross
at the point of contact instead of touching, i.e. if $f\neq g$ at the point of
contact, then $\sigma^{2}_{TSEP}>1$. The precise value of $\sigma^{2}_{TSEP}$
is given by (4.9). Moreover, the value of $\sigma^{2}_{TSEP}$ increases as the
value of $f$ and $g$ diverges at the point of contact. On one hand, if $f$ is
much larger than $g$, then part B of Lemma 3 implies $\sigma^{2}_{TSEP}$ is
close to $\lambda^{-1}$. On the other hand, if $g$ is much larger than $f$,
then part C of Lemma 3 indicates that $\sigma^{2}_{TSEP}$ is close to
$(1-\lambda)^{-1}$. These bounds are tight because, by part B of Lemma 3,
$\sigma^{2}_{TSEP}$ can not be larger than
$\max\\{\lambda^{-1},(1-\lambda)^{-1}\\}$ under the set up of Theorem 1.

The above discussion leads to the following conclusion for the TSEP test. If
we want to control the asymptotic type I error of the TSEP test at all null
configurations, then we should use the critical value $C_{m,n}z_{\alpha}$
where $C_{m,n}=\max\\{\sqrt{N/m},\sqrt{N/n}\\}$. There is a caveat, however.
To see this, we begin by noting that $C_{m,n}\geq\sqrt{2}$, with equality
holding only when $m=n$. For example, for our motivating dataset,
$C_{m,n}\approx 1.9$. However, if $m/N\to 0$ or $n/N\to 0$, then
$C_{m,n}\to\infty$. Thus if $m$ and $n$ are not close to being equal,
$C_{m,n}$ can be a large quantity. Therefore, using $C_{m,n}z_{\alpha}$ as
critical value yields a conservative test. Hence, we will call the
corresponding TSEP tests the conservative TSEP tests. Our simulations indicate
that conservative TSEP tests have poorer power compared to the minimum t-test
even when $m$ and $n$ are equal, and their power keeps degrading as $C_{m,n}$
increases.

If we use instead use the critical value $z_{\alpha}$, then we control the
asymptotic type I error at the null configurations with
$C_{p}\subset\text{int}(D_{p})$ or those with $f=g$ at $C_{p}$. This only
excludes the null cases where $F$ and $G$ may cross at $\text{bd}(D_{p})$ (see
Figure 3 iv). These type of configurations can be considered pathological
cases. Moreover, our simulations in Section 4.3 (see case b) show that even
when $F$ and $G$ cross at $\text{bd}(D_{p})$, the TSEP tests with critical
value $z_{\alpha}$ control the type I error. The TSEP tests with critical
value $z_{\alpha}$ have decent power and their overall performance is
comparable with the minimum t-tests. In view of the above, we recommend using
the asymptotic critical value $z_{\alpha}$ when using the TSEP test.

Our final result on the minimum t-test and the TSEP test establishes their
asymptotic consistency.

######  Theorem 2.

Suppose $(F,G)\in{\mathcal{H}}_{1}$ satisfy Condition Condition N. Then if
$m/N\to\lambda$, then

$\lim_{m,n\to\infty}T_{m,n}^{\text{min}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{p}\infty\quad\text{and}\quad
T_{m,n}^{\text{tsep}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{p}\infty.$

The asymptotic distribution of the WRS statistic is well-established in the
literature. Suppose $F$ and $G$ are continuous distribution functions. In that
case, it is well known that, when $F=G$, the WRS statistic is asymptotically
distributed as a standard gaussian random variable, i.e.
$T_{m,n}^{\text{wrs}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{d}N(0,1)$ (Dwass,
1956). If $(F,G)\in H^{a}_{0}$ satisfies $\int_{\mathbb{R}}G(z)dF(z)>1/2$,
however, $T_{m,n}^{\text{wrs}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{p}-\infty$,
whereas for $(F,G)\in{\mathcal{H}}_{1}^{a}$, we have
$T_{m,n}^{\text{wrs}}(\mathbb{F}_{m},\mathbb{G}_{n})\to_{p}\infty$.

####  4.2.2 Asymptotic critical values of the shape-constrained tests

We will show that under some additional conditions, the difference between the
nonparametric and the shape-constrained test statistics is $o_{p}(1)$, which
automatically implies that the shape-constrained tests enjoy the same
asymptotic properties as the nonparametric tests.

For the unimodal case, the additional condition is a curvature condition,
which requires $f$ and $g$ to be nowhere flat within their respective domains.

######  Condition A.

For the density $\mu$, the Lebesgue measure of the set
$\\{\mu^{\prime}=0,\mu>0\\}$ is $0$, where $\mu^{\prime}$ is the derivative of
$\mu$.

For densities satisfying Condition Condition A,
$\sqrt{m}(\widehat{F}_{m}^{0}-F)$ almost surely weakly converges to
$\mathbb{V}\circ F$, where $\mathbb{V}$ is a Brownian bridge. In this case, it
can be shown that (Kiefer and Wolfowitz, 1976)
$\sqrt{m}\|\widehat{F}_{m}^{0}-\mathbb{F}_{m}\|_{\infty}=o_{p}(1)$. Our Lemma
A.5 in Appendix A.1 states that
$\sqrt{m}\|\widehat{F}_{m}^{0}-\widehat{F}_{m}\|_{\infty}=o_{p}(1)$, which
implies $\sqrt{m}\|\widehat{F}_{m}-\mathbb{F}_{m}\|_{\infty}=o_{p}(1)$ in this
case. However, under the violation of Condition Condition A, the process
$\sqrt{m}(\widehat{F}_{m}^{0}-F)$ no longer converges to $\mathbb{V}\circ F$
weakly. The densities that violate Condition Condition A form the boundary of
the class of unimodal densities. For these densities,
$\sqrt{m}\|\widehat{F}_{m}^{0}-\mathbb{F}_{m}\|_{\infty}$, and hence,
$\sqrt{m}\|\widehat{F}_{m}-\mathbb{F}_{m}\|_{\infty}$ is non-negligible. The
limiting process of $\sqrt{m}(\widehat{F}_{m}^{0}-F)$ in this case is slightly
convoluted, and we refer to Beare et al. (2017); Carolan and Dykstra (1999);
Carolan (2002) for more details on the limiting process. Just to give an
example, if $f\sim U[0,1]$, then the limiting process is the least concave
majorant of a Brownian bridge (Carolan and Dykstra, 1999). Condition Condition
A is thus required to ensure that neither $f$ nor $g$ is one of these
problematic boundary densities.

######  Lemma 4.

Suppose that $f$ and $g$ are unimodal densities satisfying Condition Condition
A. Further suppose that $f$ and $g$ are bounded away from $0$ on an open set
containing $D_{p}$, and $m,n$ satisfy $m/N\to\lambda$. Then,

$|T_{i}(\mathbb{F}_{m},\mathbb{G}_{n})-T_{i}(\widehat{F}_{m},\widehat{G}_{n})|\to_{a.s.}0,\quad\text{
for }i=1,2,3.$

In case of the log-concave test statistics, however, we require a smoothness
condition as well as a curvature condition. We will quantify smoothness via a
Hölder condition. For a compact set $K\subset\mathbb{R}$, a function $h$ is
said to be in the Hölder class $\mathcal{H}^{\beta,L}(K)$ with exponent
$\beta\in[1,2]$ and constant $L>0$ if, for all $x,y\in K$, $|h(x)-h(y)|\leq
L|x-y|$ if $\beta=1$ and $|h^{\prime}(x)-h^{\prime}(y)|\leq
L|x-y|^{\beta-1}\quad\text{if }\beta>1.$
We say that a density $\mu$ ($\mu=f$ or $g$) satisfies Condition B1 if the
following holds.

######  Condition B1.

There exists $\beta\in[1,2]$, $L>0$, and a compact $K\subset\mathbb{R}$ such
that the density $\mu$ satisfies $\log\mu\in\mathcal{H}^{\beta,L}(K)$.

In addition to Condition Condition B1, we also require $f$ and $g$ to satisfy
a curvature condition.

######  Condition B2.

$\mu$ is a log-concave density with log-density $\phi=\log\mu$. Suppose
$K\subset\text{dom}(\phi)$ is compact. Then there exists $C>0$ such that all
$x,y\in K$ satisfying $x<y$ obeys

$\displaystyle\begin{aligned} \phi^{\prime}(x)-\phi^{\prime}(y)\geq
C(y-x),\end{aligned}$

where $\phi^{\prime}$ is the left derivative or the right derivative of
$\phi$.

Note that, since $\phi$ is concave, its left and right derivatives always
exist. If $\phi^{\prime}$ is differentiable on $K$, Condition Condition B2
reads as $\phi^{\prime\prime}(x)\leq-C$ for $x\in K$. Conditions of type
Condition B1 and Condition B2 also appear in Dümbgen and Rufibach.

######  Lemma 5.

Suppose that $f$ and $g$ are log-concave densities satisfying Conditions
Condition B1 and Condition B2. Suppose, further, $f$ and $g$ are bounded away
from $0$ on an open set containing $D_{p}$, and $m/N\to\lambda$. Then it
follows that
$|T_{i}(\mathbb{F}_{m},\mathbb{G}_{n})-T_{i}(\tilde{F}_{m},\tilde{G}_{n})|=o_{p}(1)$
for $i=1,2$.

######  Remark 4.

All results of Section 4.2 hold if we replace $D_{p}$ by a compact set $D$ as
long as there is a $p>0$ so that $D\subset D_{p}$ and $D_{p}$ satisfies the
conditions of the theorems stated in Section 4.2.

###  4.3 Simulations

This section compares the performance of the shape-constrained tests designed
in Section 4.1 with their nonparametric counterparts. We let $m=n=100$, which
is reflective of the sample sizes anticipated in many phase 1b or phase 2
vaccine trials — for example, our motivating dataset has $m+n=248$. For all
TSEP tests in this section, we use the critical value $z_{\alpha}$. See
Appendix C for the simulations with TSEP tests that have critical value
$C_{m,n}z_{\alpha}$. For the TSEP and the minimum t-test, we also include the
tests based on the smoothed log-concave MLE. We will refer to the
corresponding test as the smoothed log-concave test. Although we do not have
any theoretical result for this test, we use the asymptotic critical value
$z_{\alpha}$.

By the design of our hypotheses, ${\mathcal{H}}_{0}$ encompasses a broad
number of cases ranging from $G\succ F$, $F=G$ to cases where $F$ and $G$
touch or cross each other. We develop simulation schemes so that we can
explore a wide range of scenarios. Our simulation schemes involves a parameter
$\gamma$ varying over the range $[0,1]$. Here $\gamma$ quantifies the
difference between the data generating distribution functions $F_{\gamma}$ and
$G_{\gamma}$. We evaluate the power $\nu(\gamma)$ at a grid of equally spaced
points in $[0,1]$.

For our simulation study, we consider the following cases:

  * (a)

$F_{\gamma}\sim N(\gamma,1)$, and $G_{\gamma}\sim N(0,1)$.

  * (b)

$F_{\gamma}\sim N(3\gamma,1)$, and $G_{\gamma}\sim N(0.5,2)$.

  * (c)

$F_{\gamma}\sim Gamma(2,0.1+0.4\gamma)$ and $G_{\gamma}\sim Gamma(1,0.5)$,
where $Gamma(a,b)$ is a Gamma random variable with shape parameter $a$ and
scale parameter $b$.

  * (d)

$F_{\gamma}\sim Gamma(2,1)$ and $G_{\gamma}\sim Pareto(0.5+2\gamma,1)$. Here
$Pareto(a,b)$ is the Pareto distribution function with shape parameter $a$ and
scale parameter $b$.

  * (e)

$F_{\gamma}\sim N(0,1)$ and $G_{\gamma}\sim
N(2\gamma+4,1)/2+N(2\gamma-2,1)/2$.

Figure 5: Plots of the densities $f_{\gamma}$ and $g_{\gamma}$ corresponding
to the cases (a), (d), and (e) for $\gamma=1$.

Case (a) corresponds to the traditional setting of null of equity against a
shift alternative. The last four cases cover ${\mathcal{H}}_{0}$ combinations
that include crossing. In those cases, $F_{\gamma}$ and $G_{\gamma}$ cross
each other on $\text{int}(D_{p})$ at $\gamma=0$. As $\gamma$ increases,
however, $F_{\gamma}$ and $G_{\gamma}$ eventually touch (cases c, d, and e),
or cross (case b) each other at $\text{bd}(D_{p})$, generating a LFC
configuration. We denote the corresponding $\gamma$ by $\gamma^{*}$. Finally,
at $\gamma=1$, $F_{\gamma}$ strictly dominates $G_{\gamma}$ in the sense of
${\mathcal{H}}_{1}$. Figures 18i, 17i, 17ii, and 18ii in Appendix D display
the plots of $F_{\gamma}$ and $G_{\gamma}$ for several values of $\gamma$ in
cases (b), (c), (d) and (e), respectively.  Figure 5 illustrates the densities
in cases (a), (d), and (e) for $\gamma=1$. Cases (d) and (e) are chosen to
reflect violations of the shape constraints that we consider. The Pareto
density in (d) violates the log-concavity assumption and the normal mixture in
(e) violates both the unimodality and the log-concavity assumptions. All other
densities satisfy both shape constraints.

We evaluate the properties of the tests under consideration using 10,000 Monte
Carlo replicates. We set the $p$ in $D_{p,m,n}$ to be $0.05$, where
$D_{p,m,n}$ was defined to be the set
$[\mathbb{H}_{N}^{-1}(p),\mathbb{H}_{N}^{-1}(1-p)]$. Also, the level of
significance is $0.05$ for all our tests. The tests based on the empirical
cumulative distribution function will be referred to as NP (nonparametric)
tests. For brevity, we will refer to the nonparametric, unimodal, log-concave
and smoothed log-concave tests by NP, UM, LC, and smoothed LC tests,
respectively.

Figure 7 displays the power curves for the minimum t-test and the TSEP test.
In terms of power, the LC and smootheed LC tests generally outperform the UM
tests, which generally outperform the NP tests. In their simulation study,
Davidson and Duclos (2013) reported the NP minimum t-test to be conservative,
which aligns with our observation. Also, the shape-constrained minimum t-tests
have slightly higher power than the shape-constrained TSEP tests, although the
difference is not always significant. For the NP tests, the rejection rates of
the minimum t-tests and the TSEP tests are almost identical.

Except for in case (d), where $G_{\gamma}$ is Pareto, all tests control the
type I error in all null set-ups, including the LFC configuration where
$\gamma=\gamma^{*}$. Since LFC configurations constitute the boundary of
${\mathcal{H}}_{0}$, this observation indicates that our tests have size
$0.05$ for all cases except case (d). Although Pareto density violates only
the log-concavity assumption, apparently no test has the correct size, albeit
the NP and UM tests performing the best in terms of size. Also, for case (d),
the smoothed LC test has lower type I error than LC test. Surprisingly, in
case (e), where both log-concavity and unimodality are violated, the shape-
constrained tests have overall better performance than the NP tests although
all tests exhibit poor power in this case.

All the WRS type tests, including the LC WRS test (whose asymptotic behavior
is yet unknown), exhibit a much larger size than 0.05 in all cases except the
null of equity type case (a), which is the ideal scenario for WRS type tests.
Figure 8 illustrates the power curve $\nu(\gamma)$ in cases (b) and (c), where
all three WRS tests exhibit very high type I error at several null
configurations. Our findings are consistent with the observation in Ledwina
and Wyłupek (2012) that the WRS test can be misleading for testing the null of
non-dominance.

In summary, under a correctly specified model, the shape-constrained minimum
t-test and TSEP test outperform their nonparametric counterparts, with the
log-concave minimum t-test having the best power. When the shape constraints
are violated, the nonparametric tests are not distinguishably better than the
unimodal tests. On the other hand, The WRS tests do not control type I error
for null configurations with crossing, as expected. Another important takeaway
from this section is that with the critical value $z_{\alpha}$, the smoothed
log-concave tests seem to perform as well as the log-concave tests under log-
concavity.

The current article uses asymptotic critical values for performing the above-
mentioned tests, but bootstrap critical values could also be an option.
Bootstrap for the shape-constrained tests, however, is not straightforward.
Generation of observations from shape-constrained LFC configuration poses some
challenges, which involves solving non-trivial optimization problems. Further
discussion in this direction is out of the scope of the present paper.
Therefore, we leave bootstrap tests for future research.

###  4.4 Application to HVTN 097 and HVTN 100 data

Our first task is to select a $p$ for choosing $D_{p,m,n}$. Figure 1 and some
inspection show that the empirical distribution functions are very close on
the sets $(-\infty,{\mathcal{H}}^{-1}(0.072)]$ and
$[{\mathcal{H}}^{-1}(0.975),\infty)$ in that they either cross or touch each
other on these regions. Therefore, this is the problematic region we wish to
exclude from our $D_{p,m,n}$, because clearly there is not enough evidence of
any dominance in this region. Therefore, we set the $p$ in $D_{p,m,n}$ to be
$0.075$. We remark that it is ideal to choose $p$ in a systematic way without
looking at the data. However, constructing a rigorous procedure for choosing
$p$ is out of the scope of the present paper, and we leave it for future
research. Table 3, which tabulates the p-values of the tests, displays that
all tests reject the null at the level of significance $0.05$. The highest
p-value is observed for the NP TSEP test, which is approximately $0.042$.

Tests | Nonparametric | Unimodal | Log-concave | Smoothed log-concave
---|---|---|---|---
Minimum t-test | 0.015 | 0.007 | 0.007 | 0.001
TSEP | 0.042 | 0.037 | 0.035 | 0.030
Table 3: Table of the p-values of different tests applied on our data. Here
the $p$ in $D_{p,m,n}$ is set to be $0.075$. The critical value was
$z_{\alpha}$ for all the tests.

It is natural to ask if we can estimate the power of our tests at
$(F_{097},F_{100})$. Because $(F_{097},F_{100})$ is unavailable, we analyze
the power in a neighborhood of $(\tilde{F}^{sm}_{097},\tilde{F}^{sm}_{100})$
instead, where $\tilde{F}^{sm}_{097}$ and $\tilde{F}^{sm}_{100}$ correspond to
the distributions of the smoothed log-concave MLE (Chen and Samworth, 2013)
estimators of $f_{100}$ and $f_{097}$, respectively.

Let us denote the smoothed log-concave MLE of the pooled sample by
$\tilde{f}_{m,n}^{0}$. Letting $\tilde{F}_{m,n}^{0}$ denote the corresponding
distribution function, we consider the mixture distributions

$\tilde{F}^{sm}_{097}(\gamma)=(1-\gamma)\tilde{F}_{m,n}^{0}+\gamma\tilde{F}^{sm}_{097},\quad\text{
and
}\quad\tilde{F}^{sm}_{100}(\gamma)=(1-\gamma)\tilde{F}_{m,n}^{0}+\gamma\tilde{F}^{sm}_{100},$
(4.11)

where $\gamma\in[0,1]$. Note that, similar to case (a) in Section 4.3, here
also $\gamma$ quantifies the departure of the configuration
$(\tilde{F}^{sm}_{097}(\gamma),\tilde{F}^{sm}_{100}(\gamma))$ from the null of
equality of distributions. Also, the distance between
$\tilde{F}^{sm}_{097}(\gamma)$ and $\tilde{F}^{sm}_{100}(\gamma)$ increases as
$\gamma$ approaches $1$. We denote the densities of
$\tilde{F}^{sm}_{097}(\gamma)$ and $\tilde{F}^{sm}_{100}(\gamma)$ by
$\tilde{f}^{sm}_{097}(\gamma)$ and $\tilde{f}^{sm}_{100}(\gamma)$,
respectively.

Now observe that when $\gamma\in(0,1)$, the mixture densities may not be log-
concave or even unimodal. Hence, we compute the log-concave projections
(Dümbgen et al., 2011) of $\tilde{f}^{sm}_{097}(\gamma)$ and
$\tilde{f}^{sm}_{100}(\gamma)$, respectively. Log-concave projection of a
density $f$ is the log-concave density closest to $f$ in Kullback-Leibler (KL)
distance. Since the log-concave projection of any arbritrary density is not
directly computable, we adopt a two step approach to approximate the log-
concave projections. In the first step, we simulate $1000$ observations from
each of $\tilde{f}^{sm}_{097}(\gamma)$ and $\tilde{f}^{sm}_{100}(\gamma)$. In
the second step, we calculate the smoothed log-concave MLE density estimators
of Chen and Samworth based on the simulated samples in the last step. The
resulting densities are the approximate smoothed log-concave projections of
$\tilde{f}^{sm}_{097}(\gamma)$ and $\tilde{f}^{sm}_{100}(\gamma)$. Finally, we
generate two samples of size $68$ and $180$ from the projected densities using
the methods in Dümbgen and Rufibach (2010) and the R package logcondens, and
replicate this process 10,000 times.

Figure 9 entails that all the tests exhibit decent power for higher values of
$\gamma$. The LC minimum t-test exhibits the highest power, which is
unsurprising since the underlying data is generated from log-concave
densities. Also, since
$\tilde{F}^{sm}_{097}(\gamma)=\tilde{F}^{sm}_{100}(\gamma)$ at $\gamma=0$, the
power curves resemble that of case (a) in our simulation schemes.

##  5 Measures of discrepancy

This section describes an approach for quantifying the difference between $f$
and $g$ via estimates of the squared Hellinger distance
${\mathcal{D}}^{2}(f,g)$, which provides complementary insights to the tests
of stochastic dominance presented in Section 4. The Hellinger distance is an
example of $f$-divergence. These divergences are widely used to measure the
similarity or dissimilarity between two probability measures. Compared to
other commonly used $f$-divergences such as the KL divergence or the
$\chi^{2}$ divergence (Nielsen and Nock, 2014), the Hellinger distance is
appealing due to its symmetry in its arguments, which is a desirable property
for a measure of discrepancy. Another crucial advantage of the Hellinger
distance is that its value is finite for every pair of densities (Gibbs and
Su, 2002). In contrast, the KL and $\chi^{2}$ divergences can both be infinite
if the densities under consideration do not share the same support. It is not
always reasonable to assume that the underlying densities of responses
collected from different vaccine trials will have the same support. Therefore,
the Hellinger distance appeals to us more than the KL divergence or the
$\chi^{2}$ divergence. The Hellinger distance has also seen successful
application in various disciplines ranging from machine learning (Cieslak and
Chawla, 2009; González-Castro et al., 2013, 2010), to ecology (Rao, 1995), to
fraud detection (Yamanishi et al., 2004). Finally, we choose to work with the
squared Hellinger distance instead of the Hellinger distance because due to
its simpler form, the squared version easily lends itself to efficient
estimation procedures. Regardless, a 95% confidence interval for the Hellinger
distance can always be constructed using that of its squared version.

We estimate the squared Hellinger distance between $f$ and $g$ using the same
density estimators involved in the construction of tests of stochastic
dominance, that is, the log-concave MLE (Dümbgen and Rufibach, 2009) and its
smooth version (Chen and Samworth, 2013), or Birgé’s estimator. Recalling the
definition of the density estimators $\widehat{f}_{m}$, $\widehat{g}_{n}$,
$\tilde{f}_{m}$, $\tilde{g}_{n}$, $\tilde{f}^{sm}_{m}$, and
$\tilde{g}^{sm}_{n}$ from Section 4.1, we propose the plug-in estimators
${\mathcal{D}}^{2}(\widehat{f}_{m},\widehat{g}_{n})$,
${\mathcal{D}}^{2}(\tilde{f}_{m},\tilde{g}_{n})$, and
${\mathcal{D}}^{2}(\tilde{f}^{sm}_{m},\tilde{g}^{sm}_{n})$ for the purpose of
estimating ${\mathcal{D}}^{2}(f,g)$. We refer to the resulting estimators as
the “unimodal”, “log-concave”, and the “smoothed log-concave” estimator,
respectively.

###  5.1 Asymptotic properties of the unimodal estimator:

We will show that, under some regularity conditions,
${\mathcal{D}}^{2}(\widehat{f}_{m},\widehat{g}_{n})$ is a
$\sqrt{N}$-consistent estimator of ${\mathcal{D}}^{2}(f,g)$ with asymptotic
variance

$\sigma_{f,g}^{2}=\frac{2{\mathcal{D}}^{2}(f,g)-{\mathcal{D}}^{4}(f,g)}{4\lambda(1-\lambda)}.$
(5.1)

Letting $b,B>0$, we denote by $\mathcal{P}(b,B)$ the class of densities that
are bounded below and above by $b$ and $B$ on their support, that is,

$\mathcal{P}(b,B)=\bigg{\\{}f\in\mathcal{P}\ :\ b\leq f(x)\leq B,\ \text{ for
}x\in\text{supp}(f)\bigg{\\}}.$ (5.2)

We will assume that $f,g\in\mathcal{P}(b,B)$. Simulations suggest that this
condition may not be necessary. However, this condition is required for
technical reasons in our proof. Such technical condition is quite common in
the literature, and has appeared in the analysis of plug-in estimators
(Kandasamy et al., 2015) and functionals of Grenander estimators (Mukherjee et
al., 2019).

######  Theorem 3.

Suppose $f$ and $g$ are unimodal densities in $\mathcal{P}(b,B)$, where $b$,
$B>0$. Further suppose that $f$ and $g$ satisfy condition Condition A and
$m/N\to\lambda$. Then,

$\sqrt{N}[{\mathcal{D}}^{2}(\widehat{f}_{m},\widehat{g}_{n})-{\mathcal{D}}^{2}(f,g)]\to_{d}N(0,\sigma^{2}_{f,g}),$
(5.3)

where $\sigma^{2}_{f,g}$ is as defined in (5.1).

Note that, since $\sigma_{f,g}^{2}$ can be consistently estimated plugging in
the estimator ${\mathcal{D}}^{2}(\widehat{f}_{m},\widehat{g}_{n})$, a Wald
type confidence interval is readily available for ${\mathcal{D}}^{2}(f,g)$.
Also, the asymptotic variance $\sigma^{2}_{f,g}$ equals the lower bound of the
asymptotic variance on a regular estimator of the squared Hellinger distance
under the nonparametric model (Kandasamy et al., 2015). See Van der Vaart
(1998); Birgé and Massart (1995) for more detail on the lower bound and
related theory.

We now give a high-level explanation of the idea behind Theorem 3. It can be
shown that the squared Hellinger distance is a smooth functional of the
underlying distribution functions in some suitable sense. Specifically, we
will show that it allows a first order Von Mises expansion (Fernholz, 2012),
which has the same essence as the Taylor series expansion (Kandasamy et al.,
2015). On the other hand, for $\widehat{F}_{m}^{0}$, the unimodal MLE of $F$
based on the true mode, we show that $\sqrt{m}(\widehat{F}_{m}^{0}-F)$
converges weakly to a Brownian process almost surely under Condition Condition
A. It can then be shown via a delta-method type argument, applied on the
squared Hellinger distance functional, that (5.3) holds for the unimodal MLEs
of $F$ and $G$. The final step in proving Theorem 3 is showing that the
squared Hellinger distance between Birgé’s estimator and the MLE is small.

###  5.2 Asymptotic properties of the log-concave estimators:

Recall that our log-concave estimators of the squared Hellinger distance are
given by ${\mathcal{D}}^{2}(\tilde{f}_{m},\tilde{g}_{n})$ and
${\mathcal{D}}^{2}(\tilde{f}^{sm}_{m},\tilde{g}^{sm}_{n})$, where we remind
the reader that $\tilde{f}_{m}$, $\tilde{g}_{n}$ are the log-concave MLEs and
$\tilde{f}^{sm}_{m}$, $\tilde{g}^{sm}_{n}$ are the smoothed log-concave MLEs.
Lemma 6 implies that these estimators are strongly consistent for
${\mathcal{D}}^{2}(f,g)$ provided the shape constraint holds.

######  Lemma 6.

Suppose that the densities $f$ and $g$ are log-concave and continuous. Then,
as $m,n\to\infty$,
${\mathcal{D}}^{2}(\tilde{f}_{m},\tilde{g}_{n})\to_{a.s.}{\mathcal{D}}^{2}(f,g)$
and
${\mathcal{D}}^{2}(\tilde{f}^{sm}_{m},\tilde{g}^{sm}_{n})\to_{a.s.}{\mathcal{D}}^{2}(f,g)$.

Our simulations indicate that both the log-concave and the smoothed log-
concave plug-in estimators are $\sqrt{N}$-consistent with asymptotic variance
$\sigma^{2}_{f,g}$ when $f$ and $g$ are continuous log-concave densities.
Therefore, in our empirical study, we include Wald type confidence intervals
based on these log-concave plug-in estimators as well. Our simulations in
Section 5.3 indicate that under the violation of the continuity assumption,
the log-concave plug-in estimator may still remain $\sqrt{N}$-consistent, but
the $\sqrt{N}$-consistency of the the smoothed log-concave plug-in estimator
may fail to hold.

It may be possible to analyze the $\sqrt{N}$-consistency of the log-concave
estimators working along the lines of Kulikov and Lopuhaä (2006); Groeneboom
(1984, 1989), which pertain to the shape restriction of monotonicity. However,
we leave this investigation for future research because a detailed treatment
of the $\sqrt{N}$-consistency of the log-concave estimators is out of scope of
the present paper. We remark in passing that proving the
$\sqrt{N}$-consistency of the log-concave plug-in estimator may be easier for
some special cases, e.g., when the logarithm of $f$ and $g$ are linear or
piece-wise affine, by using the results of Kim et al. (2018). We do not pursue
this direction because it is unlikely that, for our data, $f$ and $g$ belong
to such restricted classes.

######  Remark 5.

The case of model misspecification is of natural interest in the study of
shape-constrained estimators. Suppose $f$ and $g$ are not log-concave, but
they are bounded continuous densities with finite first moments. Then it
follows that there exist unique log-concave densities $f^{*}$ and $g^{*}$,
which are almost sure limits of $\tilde{f}_{m}$ and $\tilde{g}_{n}$ in both
uniform and $L_{1}$ metric (Theorem 4 Cule and Samworth, 2010). Here $f^{*}$
and $g^{*}$ are also the log-concave projections of $f$ and $g$, respectively,
in the sense of Dümbgen et al. (2011). Using the same arguments as in the
proof of Lemma 6, it can be shown that the log-concave estimator of the
squared Hellinger distance converges almost surely to
${\mathcal{D}}^{2}(f^{*},g^{*})$. If $f$ and $g$ additionally have finite
second moments, a similar phenomenon takes place for the smoothed log-concave
MLEs as well. In this case, however, $\tilde{f}^{sm}_{m}$ and
$\tilde{g}^{sm}_{n}$ converge uniformly (and also in $L_{1}$) to different
limits $f^{**}$ and $g^{**}$, which can be interprated as the respective
smoothed versions of the log-concave projections $f^{*}$ and $g^{*}$ (cf.
Theorem 1, Chen and Samworth, 2013). In this case also, the smoothed log-
concave estimator of the squared Hellinger distance converges almost surely to
${\mathcal{D}}^{2}(f^{**},g^{**})$. In summary, if the log-concavity
assumption is violated, the log-concave plug-in estimators converge to a
different limit, whose distance from ${\mathcal{D}}^{2}(f,g)$ depends on the
departure of $f$ and $g$ from log-concavity.

KDE based plug-in estimators: The natural non-parametric comparators of the
shape constrained plug-in estimators are the KDE based plug-in estimators.
Though the latter is simple to implement, it can have a bias of order
$||\widehat{f}_{m,k}-f||_{2}+||\widehat{g}_{n,k}-g||_{2}$ (cf. Section 2 of
Robins et al., 2009), where $\widehat{f}_{m,k}$ and $\widehat{g}_{n,k}$ are
the KDEs of $f$ and $g$, respectively. The above bias decreases to zero at a
rate slower than $N^{-1/2}$ (Stupfler, 2014), thereby leading to a suboptimal
performance. See Section 5 of Kandasamy et al. (2015) for more discussion on
the disadvantages of the KDE based naïve plug-in estimators.

One can improve the naïve plug-in estimator, however, using a one-step Newton-
Raphson procedure (cf. Van der Vaart, 1998; Pfanzagl and Wefelmeyer, 1985),
which leads to a bias corrected plug-in estimator. In context of the Hellinger
distance, the bias-corrected estimator
$({\mathcal{D}}^{2})^{*}(\widehat{f}_{m,k},\widehat{g}_{n,k})$ takes the form
(Kandasamy et al., 2015)

$\displaystyle{\mathcal{D}}^{2}(\widehat{f}_{m,k},\widehat{g}_{n,k})+\int_{\mathbb{R}}\psi_{f}(x;\widehat{f}_{m,k},\widehat{g}_{n,k})d\mathbb{F}_{m}(x)+\int_{\mathbb{R}}\psi_{g}(y;\widehat{f}_{m,k},{\widehat{g}_{n,k}\vphantom{\widehat{f}_{m,k}}})d\mathbb{G}_{n}(y)$
$\displaystyle=1-\dfrac{1}{2}\bigg{(}\int_{\mathbb{R}}\sqrt{\widehat{g}_{n,k}(x)/\widehat{f}_{m,k}(x)}d\mathbb{F}_{m}(x)+\int_{\mathbb{R}}\sqrt{\widehat{f}_{m,k}(y)/\widehat{g}_{n,k}(y)}d\mathbb{G}_{n}(y)\bigg{)},$

where $\psi_{f}$ and $\psi_{g}$ are the influence functions corresponding to
the functional $(f,g)\mapsto{\mathcal{D}}^{2}(f,g)$. We will formally
introduce the influence functions in Appendix B. Note that learning the form
of the bias-corrected estimator thus requires the explicit computation of the
influence functions $\psi_{f}$ and $\psi_{g}$. From a broader prospective,
each time one tries to estimate a functional of the underlying distributions
using a bias-corrected plug-in estimator, they have to carry out some extra
analytical calculations that depend on the functional of interest.

In contrast, Theorem 3 shows that when the shape constraint is satisfied, the
unimodal estimator of the squared Hellinger distance does not require any bias
correction for $\sqrt{N}$-consistency. Our analysis in Appendix B indicates
that this is not an artifact of the Hellinger distance, but rather the result
of Birgé’s estimator’s proximity to the unimodal MLE, i.e. the Grenander
estimator based on the true mode. In fact, Theorem 4 in Appendix B shows that,
if $T$ is a smooth functional of the distribution functions $F$, then
$T(\widehat{f}^{0}_{m})$ is $\sqrt{m}$-consistent for $T(F)$ under mild
conditions. Because the total variation distance between $\widehat{f}_{m}$ and
$\widehat{f}^{0}_{m}$ is $o_{p}(m^{-1/2})$ (see our Lemma A.5 in Appendix A),
it is then natural to expect that $T(\widehat{f}_{m})$ would be
$\sqrt{m}$-consistent if $T$ is sufficiently smooth.

Although we do not have any such theoretical evidence for the log-concave
estimators, our simulations and the simulations of Cule et al. (2010) on plug-
in estimators based on log-concave MLEs (see Figure 18 therein) indicate that
these estimators do not have any $O_{p}(N^{-1/2})$ bias term either. This
allows users of correctly specified shape-constrained estimators to avoid
analytic calculation of the influence functions entirely.

In our upcoming simulation study, we use both
${\mathcal{D}}^{2}(\widehat{f}_{m,k},\widehat{g}_{n,k})$ and
$({\mathcal{D}}^{2})^{*}(\widehat{f}_{m,k},\widehat{g}_{n,k})$ as comparators,
where the KDEs are based on the Gaussian kernel. We refer to these estimators
as the KDE estimator and the bias-corrected KDE estimator, respectively. As in
Section 3, the kernel bandwidth is chosen using the univariate least square
cross-validation (LSCV) selector of Bowman (1984) and Rudemo (1982).
Generally, kernel-based bias corrected estimators satisfy $\sqrt{N}$
consistency results of the type (5.3) (cf. Theorem 6, Kandasamy et al., 2015).
Therefore, we use the bias-corrected KDE based confidence intervals to
benchmark the performance of our shape constrained estimators. However, we do
not report any confidence interval based on the naïve plug-in estimator
because our simulations indicate that usually its coverage is much less than
the nominal level. For the sake of clarity, in Table 4, we summarize the
current state of results on the above-mentioned estimators of the squared
Hellinger distance.

Estimator | Consistency |  $\sqrt{N}$-consistency | Showed in |  Confidence interval∗
---|---|---|---|---
Naïve KDE | yes | no |  Kandasamy et al. (2015) | no
Bias corrected KDE | yes | yes | Kandasamy et al. (2015) | yes
Unimodal | yes | yes | current paper | yes
Log-concave | yes |  unknown, but suggested by simulations | current paper | yes
Smoothed log-concave | yes | $---^{\prime\prime}---$ | current paper | yes
Table 4: Table summarizing the current results on different plug-in estimators
of the squared Hellinger distance. The consistency and $\sqrt{N}$-consistency
of the shape constrained estimators are proved under the correct specification
of the shape constraints. Empirical studies are provided on all the estimators
listed above.
* Indicates whether we included confidence intervals for the corresponding estimator in the empirical study of Section 5.3.

###  5.3 Simulations:

Figure 6: Plots of the densities $f$ and $g$ for the cases (a) (left), (e)
(middle), and (f) (right) in the simulation scheme of Section 5.3 Figure 7:
Plot of estimated power $\nu(\gamma)$ vs $\gamma$: the labels (a)-(e)
correspond to simulation schemes (a)-(e). Here MT and TSEP correspond to the
minimum t-test and the TSEP test, respectively. The standard deviation of the
$\nu(\gamma)$ estimate in each case is less than $0.005$. The black horizontal
line corresponds to the level of the test, $\alpha=0.05$. For cases (b)-(e),
the black vertical line represents the LFC configuration $\gamma^{*}$, taking
value 0.70 (b), 0.55 (c), 0.65 (d), and 0.80 (e).  Figure 8: Estimated power
curves of WRS type tests for case (b) (left) and case (c) (right). Here LC.WRS
and UM.WRS correspond to the tests based on
$T_{m,n}^{\text{wrs}}(\tilde{F}_{m},\tilde{G}_{n})$ and
$T_{m,n}^{\text{wrs}}(\widehat{F}_{m},\widehat{G}_{n})$, respectively. The
standard errors of the estimated powers are less than 0.005. The black
horizontal line corresponds to the level of the test, $\alpha=0.05$. The black
vertical lines represent the LFC configuration $\gamma^{*}$, which are 0.70
(b) and 0.55 (c). Figure 9: Estimated power curve applied on our datasets:
Here MT and TSEP correspond to the minimum t-test and the TSEP test,
respectively. The black horizontal line corresponds to the level of the test,
which is $0.05$. The standard deviation of the estimated power does not exceed
$0.005$ in any case.

To compare the performance of different estimators of the squared Hellinger
distance, we consider the following combinations of $f$ and $g$:

  * (a)

$f\sim Gamma(4,1)$ and $g\sim Gamma(3,1)$.

  * (b)

$f\sim N(1,1)$ and $g\sim N(0,1)$, which corresponds to case (a) in Section
4.3 with $\gamma=1$.

  * (c)

$f=\tilde{f}^{sm}_{097}$ and $g=\tilde{f}^{sm}_{100}$, where
$\tilde{f}^{sm}_{097}$ and $\tilde{f}^{sm}_{100}$ are the smoothed log-concave
MLEs of $f_{097}$ and $f_{100}$, respectively (see Figure 2).

  * (d)

$f\sim\text{Exp}(1)$ and $g\sim\text{Exp}(2)$ where $\text{Exp}(t)$ is the
exponential distribution with rate $t$.

  * (e)

$(f,g)$ corresponds to case (e) in Section 4.3 with $\gamma=1$.

  * (f)

$f\sim N(0,1)$, and $g\sim Gamma(3.61,1.41)$.

The plots of the above schemes can be found in Figure 5 and Figure 6. . In
cases (a), (b), (c), (d), and (f), both shape constraints are satisfied, where
in case (e), both shape constraints are violated. We will refer to case (a),
(b), (c), (d), and (f), therefore, as the “correctly specified” cases and case
(e), as the “misspecified” case. In case (d), the logarithm of the densities
are linear on their respective supports. This is also the only case where the
densities are discontinuous (see Figure 6; discontinuity at zero). In case
(f), $g$ was chosen so as to resemble the density estimate of $f_{100}$ (see
Figure 2 and Figure 6). This is the only correctly specified case where the
densities are known to come from different families, and the densities also
have very different shapes.

We generate two samples of same size from each simulation setting. We vary the
common sample size $n$ from $50$ to $500$ in increments of $50$. We do not
consider larger values of $n$ because in our motivating phase 1b and phase 2
vaccine trial applications, the sample sizes are generally no larger than 500.
We consider 10,000 Monte Carlo replications for each sample size.

Figure 10: Plots of the absolute values of the bias (scaled by $\sqrt{n}$) for
cases (a)–(d). Here “KDE (BC)” and “smoothed LC” stand for the bias corrected
KDE estimator and the smoothed log-concave estimator, respectively. For each
estimator, either the standard error is less than 0.006 or the relative
standard error is less than 2%.  Figure 11: Plots of the bias (scaled by $n$)
for cases (e) and (f). Here “KDE (BC)” and “smoothed LC” stand for the bias
corrected KDE estimator and the smoothed log-concave estimator, respectively.
For each estimator, either the standard error is less than 0.006 or the
relative standard error is less than 2%.

i Plots of the MSE (scaled by $n$) for cases (a)–(c).

ii Plots of the MSE (scaled by $n$) for cases (d)–(f).

Figure 12: Plot of MSE: Here “KDE (BC)” and “smoothed LC” stand for the bias
corrected KDE estimator and the smoothed log-concave estimator, respectively.
For each estimator, either the standard error is less than 0.006 or the
relative standard error is less than 2%.

i Coverage probability of the 95$\%$ confidence intervals for cases (a)–(f).

ii Average length of the 95$\%$ confidence intervals for cases (a)–(f).

Figure 13: Here “KDE (BC)” and “smoothed LC” stand for the bias corrected KDE
estimator and the smoothed log-concave estimator, respectively. In all the
cases, the standard error is less than 0.005, and hence, not plotted.

Figures 10 and 11 plot $\sqrt{n}$ times the absolute bias, and Figure 12 plots
$n$ times the mean squared error (MSE). Figure 13i and Figure 13ii display,
respectively, the coverage and the average length of the confidence intervals,
both of which are estimated using 10,000 Monte Carlo samples. In the
misspecified case (e), the performance of the shape-constrained estimators
deteriorate sharply with $n$, which is unsurprising due to the violation of
the shape constraints in this case. In case (d), where the densities are
exponential, the scaled bias of the smoothed log-concave plug-in estimator
increases with $n$, which implies $\sqrt{N}$-consistency does not hold for
this estimator in this case. Closer inspection reveals that although it is not
$\sqrt{N}$-consistent, the smoothed log-concave plug-in estimator is still
consistent in case (d). We found out that the smoothed log-concave density
estimator fails to approximate the exponential densities near zero, their
point of discontinuity (see Figure 19 in Appendix D). It is worth mentioning
that although the curvature Condition Condition B2 is violated in case (d), it
does not affect the performance of the log-concave plug-in estimator.

For all other cases, the smoothed log-concave estimator exhibits the best
performance among the shape-constrained estimators. In all cases, the unimodal
estimator underperforms. Importantly, Figure 13i indicates that the unimodal
estimator would require sample size larger than 500 for the Wald type
confidence interval to be valid, where for the other confidence intervals,
this sample size is sufficient for the asymptotics to kick in.

The KDE-based plug-in estimator experiences an increase in the bias and the
MSE with $n$, which agrees with our previous discussion on KDE based plug-in
estimators. Its bias corrected version performs comparably to the shape-
constrained estimators in all cases except case (f). In case (f), however, the
bias corrected estimator yields a confidence interval with poor coverage,
which is probably due to the large bias incurred by the original KDE-based
plug-in estimator in this case (see panel (f) of Figure 11). Case (f) is a
case where the underlying densities differ by both shape and scale. We suspect
that the cross-validated bandwidth for the KDE estimator does not work in this
case even with bias correction. Finally, the average length of the confidence
intervals do not vary noticeably across different methods.

In summary, the log-concave plug-in estimator exhibits reliable performance
when the shape restriction holds. The smoothed log-concave plug-in estimator
may also require the underlying densities to be continuous, but otherwise it
performs comparably with the bias-corrected KDE-based plug-in estimator. The
latter performs well in many settings, but it is not always reliable. The
lacking performance of the KDE based methods in some settings is probably due
to the variable nature of the optimal bandwidth under different settings.

###  5.4 Application to HVTN 097 and HVTN 100 data:

Table 5 tabulates the point estimates of the squared Hellinger distance
(between $f_{100}$ and $f_{097}$) and the corresponding 95% intervals. Table 5
also displays the 95% confidence intervals of the Hellinger distance, which
are obtained by taking square root of the upper and lower bounds of the
previous confidence intervals.

Estimator |  Naïve KDE |  Bias -corrected KDE | UM | LC |  Smoothed LC
---|---|---|---|---|---
Point estimate of ${\mathcal{D}}^{2}(f_{097},f_{100})$ | 0.16 | 0.19 | 0.18 | 0.15 | 0.21
95% CI for ${\mathcal{D}}^{2}(f_{097},f_{100})$ |  not available | $(0.110,0.274)$ | $(0.128,0.300)$ | $(0.098,0.256)$ | $(0.079,0.228)$
95% CI for ${\mathcal{D}}(f_{097},f_{100})$ | $--^{\prime\prime}--$ | $(0.332,0.523)$ | $(0.358,0.548)$ | $(0.313,0.506)$ | $(0.281,0.447)$
Table 5: Table of the point estimates of ${\mathcal{D}}^{2}(f_{097},f_{100})$,
and 95% confidence intervals of ${\mathcal{D}}^{2}(f_{097},f_{100})$ and
${\mathcal{D}}(f_{097},f_{100})$.

To give the reader some perspective, if $f_{0}$ is a $N(0,1)$ distribution and
$f_{\mu}$ is a $N(\mu,1)$ distribution, then, when $\mu$ is equal to 1.00,
1.25, 1.50, 1.75, and 2.00, ${\mathcal{D}}(f_{0},f_{\mu})$ is equal to 0.346,
0.424, 0.500, 0.566, and 0.624, respectively. Also, in Figure 14, we display
some density pairs $(f,g)$ satisfying $F\succeq G$ with Hellinger distance in
the range $0.35-0.60$, which is similar to our data.

Figure 14: Plot of different density pairs and their distribution functions.
The Hellinger distance $(\mathcal{D})$ is given on top of each pair. Left:
$f\sim Gamma(25,1)$, $g\sim N(0.5,1)$; middle: $f\sim N(1.40,1),$ $g\sim
N(0,1)$; right: $f\sim Gamma(5,1)$, $g\sim Gamma(2,1)$.

##  6 Discussion

The first contribution of our work is a novel analysis of the data from the
HVTN 097 and HVTN 100 trials. All of our tests reject the null of non-
dominance in favor of the strict stochastic dominance of $F_{097}$ over
$F_{100}$. To provide further insight into the discrepancy between the two IgG
binding response distributions, we estimated the squared Hellinger distance
between the corresponding densities, which turns out to be approximately 0.20
(95% CI 0.10-0.30). We remark that our findings are consistent with those of
Bekker et al. (2018), who found that the average magnitude of IgG binding to
V1V2 antigens observed in the HVTN 100 trial is lower than that in the RV 144
trial. Although the latter used the same regimen as HVTN 097, it was conducted
in a different population (Thailand). Our findings indicate that the
difference in the magnitude of IgG binding response between HVTN 100 and RV
144 regimen may be attributable to the HIV clade difference rather than to the
difference in populations.

The outcome of our tests become meaningful when viewed against the lack of
efficacy observed in the phase 2b/3 HVTN 702 trial, which evaluated the safety
and efficacy of the HVTN 100 regimen in South Africa. A possible hypothesis
for why no efficacy was observed when the HVTN 100 regimen was evaluated in
this trial, whereas efficacy was observed when the HVTN 097 regimen was
evaluated in the RV144 trial, is that the HVTN 100 regimen leads to a lower
magnitude of IgG binding to the V1V2 region. This possibility is supported by
the observation made by Haynes et al. (2012) regarding the negative
correlation between rate of infection and the magnitude of IgG binding to V1V2
region. This hypothesis can be tested when the immune profile of the
participants in HVTN 702 trial becomes available.

Another contribution of our work relates to density estimation in the context
of vaccine trials. Based on a cross-validated analysis of the HVTN 097 and
HVTN 100 data, we believe that the log-concave density estimators of Dümbgen
and Rufibach (2009) and Chen and Samworth (2013) may yield improved density
estimation in vaccine studies. In future work, it would be worth further
validating this claim on other vaccine trial datasets.

We also made several methodological contributions. In Section 4, we introduce
three novel shape-constrained tests. These tests have the desirable asymptotic
properties of nonparametric tests, and simulations illustrate that the shape-
constrained tests have better overall performance than the nonparametric
tests. Moreover, even under the violation of the shape-constraints, their
performance is not much worse than the nonparametric tests. We also introduce
shape-constrained plug-in estimators of the squared Hellinger distance and
provide asymptotic consistency and distributional results. Our simulations
suggest that, when the shape constraint is satisfied, the log-concave plug-in
estimators exhibit overall lower MSE and absolute bias than the KDE based
plug-in estimator. In fact, they perform comparably with the bias-corrected
KDE based estimator. However, unlike the bias-corrected KDE plug-in estimator,
the log-concave plug-in estimators require neither selecting a tuning
parameter nor carrying out the analytic calculations needed to derive the
bias-correction term. Therefore, the log-concave plug-in estimators may be
preferred in settings where this shape constraint is plausible.

##  7 Acknowledgements

This work was supported by the National Institutes of Health (NIH) through
award numbers DP2-LM013340 and 5UM1AI068635-09. The content is solely the
responsibility of the authors and does not necessarily represent the official
views of the NIH. GlaxoSmithKline Biologicals SA was provided the opportunity
to review a preliminary version of this manuscript, but the authors are solely
responsible for final content and interpretation.

## References

  * Álvarez-Esteban et al. (2016) Álvarez-Esteban, P., Del Barrio, E., Cuesta-Albertos, J., and Matrán, C. (2016).  A contamination model for the stochastic order.  Test, 25, 751–774.
  * Ambrosio et al. (2000) Ambrosio, L., Fusco, N., and Pallara, D. (2000).  Functions of bounded variation and free discontinuity problems.  Oxford Science Publications.
  * Asmussen and Lehtomaa (2017) Asmussen, S. and Lehtomaa, J. (2017).  Distinguishing log-concavity from heavy tails.  Risks, 5(1), 10.
  * Balabdaoui et al. (2009) Balabdaoui, F., Rufibach, K., and Wellner, J. A. (2009).  Limit distribution theory for maximum likelihood estimation of a log-concave density.  Ann. Statist., 37, 1299.
  * Beare et al. (2017) Beare, B. K., Fang, Z., et al. (2017).  Weak convergence of the least concave majorant of estimators for a concave distribution function.  Electron. J. Stat., 11, 3841–3870.
  * Bekker et al. (2018) Bekker, L.-G., Moodie, Z., Grunenberg, N., Laher, F., Tomaras, G. D., Cohen, K. W., Allen, M., Malahleha, M., Mngadi, K., Daniels, B., et al. (2018).  Subtype c alvac-hiv and bivalent subtype c gp120/mf59 hiv-1 vaccine in low-risk, hiv-uninfected, south african adults: a phase 1/2 trial.  The Lancet HIV, 5(7), e366–e378.
  * Billingsley (2013) Billingsley, P. (2013).  Convergence of probability measures.  John Wiley & Sons.
  * Birgé (1997) Birgé, L. (1997).  Estimation of unimodal densities without smoothness assumptions.  Ann. Statist., 25, 970–981.
  * Birgé and Massart (1995) Birgé, L. and Massart, P. (1995).  Estimation of integral functionals of a density.  Ann. Statist., 23, 11–29.
  * Bobkov and Ledoux (2016) Bobkov, S. and Ledoux, M. (2016).  One-dimensional empirical measures, order statistics and Kantorovich transport distances.  preprint.
  * Bowman (1984) Bowman, A. W. (1984).  An alternative method of cross-validation for the smoothing of density estimates.  Biometrika, 71, 353–360.
  * Carolan and Dykstra (1999) Carolan, C. and Dykstra, R. (1999).  Asymptotic behavior of the grenander estimator at density flat regions.  Can. J. Stat., 27, 557–566.
  * Carolan (2002) Carolan, C. A. (2002).  The least concave majorant of the empirical distribution function.  Canadian Journal of Statistics, 30(2), 317–328.
  * Chen and Samworth (2013) Chen, Y. and Samworth, R. J. (2013).  Smoothed log-concave maximum likelihood estimation with applications.  Statist. Sinica, 23, 1373–1398.
  * Cieslak and Chawla (2009) Cieslak, D. A. and Chawla, N. V. (2009).  A framework for monitoring classifiers’ performance: when and why failure occurs?  Knowl. Inf. Syst., 18, 83–108.
  * Cule and Samworth (2010) Cule, M. and Samworth, R. (2010).  Theoretical properties of the log-concave maximum likelihood estimator of a multidimensional density.  Electron. J. Statist., 4, 254–270.
  * Cule et al. (2010) Cule, M., Samworth, R., and Stewart, M. (2010).  Maximum likelihood estimation of a multidimensional log-concave density.  J. Royal Stat. Soc: Series B, 72, 545–607.
  * Davidson and Duclos (2013) Davidson, R. and Duclos, J.-Y. (2013).  Testing for restricted stochastic dominance.  Econometric Rev., 32, 84–125.
  * DeFauw (2011) DeFauw, M. C. (2011).  Cost Effectiveness Analysis in Healthcare Decision-Making: Stochastic Modeling and Statistical Inference. Ph.D. thesis.
  * Doss and Wellner (2016) Doss, C. R. and Wellner, J. A. (2016).  Global rates of convergence of the MLEs of log-concave and $s$-concave densities.  Ann. Statist., 44, 954–981.
  * Dümbgen and Rufibach (2009) Dümbgen, L. and Rufibach, K. (2009).  Maximum likelihood estimation of a log-concave density and its distribution function: Basic properties and uniform consistency.  Bernoulli, 15, 40–68.
  * Dümbgen and Rufibach (2010) Dümbgen, L. and Rufibach, K. (2010).  logcondens: Computations related to univariate log-concave density estimation.  J. Stat. Softw., 39, 1–28.
  * Dümbgen et al. (2011) Dümbgen, L., Samworth, R., and Schuhmacher, D. (2011).  Approximation by log-concave distributions, with applications to regression.  Ann. Statist., 39, 702–730.
  * Dwass (1956) Dwass, M. (1956).  The large-sample power of rank order tests in the two-sample problem.  Ann. Math. Statist., pages 352–374.
  * Efromovich (2008) Efromovich, S. (2008).  Nonparametric curve estimation: methods, theory, and applications.  Springer Science & Business Media.
  * Eggermont and LaRiccia (2000) Eggermont, P. P. B. and LaRiccia, V. N. (2000).  Maximum likelihood estimation of smooth monotone and unimodal densities.  Ann. Statist., 28, 922 – 947.
  * Fernholz (2012) Fernholz, L. T. (2012).  Von Mises calculus for statistical functionals, volume 19.  Springer Science & Business Media.
  * Gibbs and Su (2002) Gibbs, A. L. and Su, F. E. (2002).  On choosing and bounding probability metrics.  International statistical review, 70(3), 419–435.
  * González-Castro et al. (2010) González-Castro, V., Alaiz-Rodríguez, R., Fernández-Robles, L., Guzmán-Martínez, R., and Alegre, E. (2010).  Estimating class proportions in boar semen analysis using the Hellinger distance.  IEA/AIE, pages 284–293.
  * González-Castro et al. (2013) González-Castro, V., Alaiz-Rodríguez, R., and Alegre, E. (2013).  Class distribution estimation based on the Hellinger distance.  Inf. Sci., 218, 146–164.
  * Gray et al. (2019) Gray, G. E., Huang, Y., Grunenberg, N., Laher, F., Roux, S., Andersen-Nissen, E., De Rosa, S. C., Flach, B., Randhawa, A. K., Jensen, R., Swann, E. M., Bekker, L.-G., Innes, C., Lazarus, E., Morris, L., Mkhize, N. N., Ferrari, G., Montefiori, D. C., Shen, X., Sawant, S., Yates, N., Hural, J., Isaacs, A., Phogat, S., DiazGranados, C. A., Lee, C., Sinangil, F., Michael, N. L., Robb, M. L., Kublin, J. G., Gilbert, P. B., McElrath, M. J., Tomaras, G. D., and Corey, L. (2019).  Immune correlates of the thai rv144 hiv vaccine regimen in south africa.  Sci. Transl. Med., 11.
  * Groeneboom (1984) Groeneboom, P. (1984).  Estimating a monotone density.  Department of Mathematical Statistics.
  * Groeneboom (1989) Groeneboom, P. (1989).  Brownian-motion with a parabolic drift and airy functions.  Probability theory and related fields, 81(1), 79–109.
  * Hall and Huang (2002) Hall, P. and Huang, L.-S. (2002).  Unimodal density estimation using kernel methods.  Statist. Sinica, pages 965–990.
  * Hannah and Dunson (2012) Hannah, L. and Dunson, D. (2012).  Ensemble methods for convex regression with applications to geometric programming based circuit design.  arXiv:1206.4645.
  * Haynes et al. (2012) Haynes, B. F., Gilbert, P. B., McElrath, M. J., Zolla-Pazner, S., Tomaras, G. D., Alam, S. M., Evans, D. T., Montefiori, D. C., Karnasuta, C., Sutthent, R., Liao, H.-X., DeVico, A. L., Lewis, G. K., Williams, C., Pinter, A., Fong, Y., Janes, H., DeCamp, A., Huang, Y., Rao, M., Billings, E., Karasavvas, N., Robb, M. L., Ngauy, V., de Souza, M. S., Paris, R., Ferrari, G., Bailer, R. T., Soderberg, K. A., Andrews, C., Berman, P. W., Frahm, N., De Rosa, S. C., Alpert, M. D., Yates, N. L., Shen, X., Koup, R. A., Pitisuttithum, P., Kaewkungwal, J., Nitayaphan, S., Rerks-Ngarm, S., Michael, N. L., and Kim, J. H. (2012).  Immune-correlates analysis of an HIV-1 vaccine efficacy trial.  N. Engl. J. Med., 366, 1275–1286.
  * Jankowski (2010) Jankowski, H. K. (2010).  Discussion of “estimation of multi-dimensional log-concave density” by Cule, Samworth, and Stewart.
  * Johnson et al. (2018) Johnson, A. L., Jiang, D. R., et al. (2018).  Shape constraints in economics and operations research.  Statist. Sci., 33, 527–546.
  * Kandasamy et al. (2015) Kandasamy, K., Krishnamurthy, A., Póczos, B., Wasserman, L. A., and Robins, J. M. (2015).  Nonparametric von mises estimators for entropies, divergences and mutual informations.  In Adv. Neural. Inf. Process. Syst., pages 397–405.
  * Kaur et al. (1994) Kaur, A., Prakasa Rao, B., and Singh, H. (1994).  Testing for second-order stochastic dominance of two distributions.  Econometric Theory, 10, 849–866.
  * Kiefer and Wolfowitz (1976) Kiefer, J. and Wolfowitz, J. (1976).  Asymptotically minimax estimation of concave and convex distribution functions.  Z. Wahr. verw. Gabiete, 34, 73–85.
  * Kim et al. (2018) Kim, A. K., Guntuboyina, A., Samworth, R. J., et al. (2018).  Adaptation in log-concave density estimation.  Annals of Statistics, 46(5), 2279–2306.
  * Kulikov and Lopuhaä (2006) Kulikov, V. and Lopuhaä, H. (2006).  The limit process of the difference between the empirical distribution function and its concave majorant.  Statistics & Probability Letters, 76, 1781–1786.
  * Laha (2021) Laha, N. (2021).  Adaptive estimation in symmetric location model under log-concavity constraint.  Electron. J. Stat., 15, 2939–3014.
  * Laha and Luedtke (2020) Laha, N. and Luedtke, A. (2020).  SDNNtests.  https://github.com/nilanjanalaha/SDNNtests/.
  * Le Breton (1991) Le Breton, M. (1991).  Stochastic orders in welfare economics.  Lecture Notes-Monograph Series, pages 190–206.
  * Ledwina and Wyłupek (2012) Ledwina, T. and Wyłupek, G. (2012).  Two-sample test against one-sided alternatives.  Scand. J. Stat., 39, 358–381.
  * Ledwina and Wyłupek (2013) Ledwina, T. and Wyłupek, G. (2013).  Tests for first-order stochastic dominance.  preprint.
  * Lee and Wolfe (1976) Lee, Y. J. and Wolfe, D. A. (1976).  A distribution-free test for stochastic ordering.  J. Amer. Statist. Assoc., 71, 722–727.
  * Lehmann (1975) Lehmann, E. L. (1975).  Nonparametrics: Statistical methods based on ranks.  San Francisco : Holden-Day.
  * Leshno and Levy (2004) Leshno, M. and Levy, H. (2004).  Stochastic dominance and medical decision making.  Health. Care. Manag. Sci., 7, 207–215.
  * Levy (1992) Levy, H. (1992).  Stochastic dominance and expected utility: Survey and analysis.  J. Manag. Sci., 38(4), 555–593.
  * Lopuhaä and Musta (2019) Lopuhaä, H. P. and Musta, E. (2019).  A central limit theorem for the hellinger loss of grenander-type estimators.  Statistica Neerlandica, 73(2), 180–196.
  * Mammen et al. (2001) Mammen, E., Marron, J. S., Turlach, B. A., and Wand, M. P. (2001).  A general projection framework for constrained smoothing.  Stat. Sci., 16(3), 232–248.
  * Meyer (2012) Meyer, M. C. (2012).  Nonparametric estimation of a smooth density with shape restrictions.  Stat. Sin., pages 681–701.
  * Miladinovic et al. (2014) Miladinovic, B., Kumar, A., Mhaskar, R., and Djulbegovic, B. (2014).  Benchmarks for detecting ‘breakthroughs’ in clinical trials: empirical assessment of the probability of large treatment effects using kernel density estimation.  BMJ open, 4, e005249.
  * Mukherjee et al. (2019) Mukherjee, R., Sen, B., et al. (2019).  On efficiency of the plug-in principle for estimating smooth integrated functionals of a nonincreasing density.  Electron. J. Stat., 13, 4416–4448.
  * Nielsen and Nock (2014) Nielsen, F. and Nock, R. (2014).  On the chi square and higher-order chi distances for approximating f-divergences.  IEEE Signal Process. Lett., 21, 10–13.
  * Pal et al. (2007) Pal, J. K., Woodroofe, M., and Meyer, M. (2007).  Estimating a pólya frequency function2.  Lecture Notes-Monograph Series, 54, 239–249.
  * Pfanzagl and Wefelmeyer (1985) Pfanzagl, J. and Wefelmeyer, W. (1985).  Asymptotic expansions for general statistical models.  Lecture notes in statistics. Springer-Verlag.
  * Pyke and Shorack (1968) Pyke, R. and Shorack, G. R. (1968).  Weak convergence of a two-sample empirical process and a new approach to Chernoff-Savage theorems.  Ann. Math. Statist., 39, 755–771.
  * Rao (1969) Rao, B. P. (1969).  Estimation of a unimodal density.  Sankhya Ser. A, 31, 23–36.
  * Rao (1995) Rao, C. R. (1995).  A review of canonical coordinates and an alternative to correspondence analysis using Hellinger distance.  Qüestiió: quaderns d’estadística i investigació operativa, 19.
  * Rerks-Ngarm et al. (2009) Rerks-Ngarm, S., Pitisuttithum, P., Nitayaphan, S., Kaewkungwal, J., Chiu, J., Paris, R., Premsri, N., Namwat, C., de Souza, M., Adams, E., et al. (2009).  Vaccination with alvac and aidsvax to prevent hiv-1 infection in thailand.  N. Engl. J. Med., 361, 2209–2220.
  * Robins et al. (2009) Robins, J., Li, L., Tchetgen, E., and van der Vaart, A. W. (2009).  Quadratic semiparametric von mises calculus.  Metrika, 69, 227–247.
  * Rockafellar (1970) Rockafellar, R. T. (1970).  Convex Analysis.  Princeton University Press.
  * Rudemo (1982) Rudemo, M. (1982).  Empirical choice of histograms and kernel density estimators.  Scand. J. Stat., 9, 65–78.
  * Samworth (2018) Samworth, R. J. (2018).  Recent progress in log-concave density estimation.  Statist. Sci., 33, 493–509.
  * Samworth and Sen (2018) Samworth, R. J. and Sen, B. (2018).  Editorial: Special issue on “nonparametric inference under shape constraints”.  Statist. Sci., 33, 469–472.
  * Shorack and Wellner (2009) Shorack, G. and Wellner, J. (2009).  Empirical Processes with Applications to Statistics.  Society for Industrial and Applied Mathematics.
  * Shorack (1984) Shorack, G. R. (1984).  Empirical and rank processes of observations and residuals.  Canad. J. Statist., 12, 319–332.
  * Shorack (2000) Shorack, G. R. (2000).  Probability for Statisticians.  Springer.
  * Sriboonchita et al. (2009) Sriboonchita, S., Nguyen, H. T., Wong, W.-K., and Dhompongsa, S. (2009).  Stochastic dominance and applications to finance, risk and economics.  Chapman and Hall/CRC.
  * Stinnett and Mullahy (1998) Stinnett, A. A. and Mullahy, J. (1998).  Net health benefits: a new framework for the analysis of uncertainty in cost-effectiveness analysis.  Med. Decis. Making., 18, S68–S80.
  * Stupfler (2014) Stupfler, G. (2014).  On the weak convergence of kernel density estimators in lp spaces.  J. Nonparametr. Stat., 26, 721–735.
  * Turnbull and Ghosh (2014) Turnbull, B. C. and Ghosh, S. K. (2014).  Unimodal density estimation using bernstein polynomials.  Comput. Statist. Data Anal., 72, 13 – 29.
  * van de Geer (2000) van de Geer, S. (2000).  Empirical Processes in M-Estimation.  Cambridge University Press.
  * Van der Vaart (1998) Van der Vaart, A. (1998).  Asymptotic Statistics.  Cambridge University Press.
  * Walther (2002) Walther, G. (2002).  Detecting the presence of mixing with multiscale maximum likelihood.  J. Amer. Statist. Assoc., 97, 508–513.
  * Walther et al. (2009) Walther, G. et al. (2009).  Inference and modeling with log-concave distributions.  Statist. Sci., 24(3), 319–327.
  * Wand and Jones (1994) Wand, M. P. and Jones, M. C. (1994).  Multivariate plug-in bandwidth selection.  Comput. Statist., 9, 97–116.
  * Whang (2019) Whang, Y. J. (2019).  Econometric Analysis of Stochastic Dominance: Concepts, Methods, Tools, and Applications.  Cambridge University Press.
  * Wolters (2012) Wolters, M. (2012).  A greedy algorithm for unimodal kernel density estimation by data sharpening.  J. Stat. Softw., 47, 1–26.
  * Wolters and Braun (2018) Wolters, M. A. and Braun, W. J. (2018).  Enforcing shape constraints on a probability density estimate using an additive adjustment curve.  Comm. Statist. Simulation Comput., 47, 672–691.
  * Woodroofe and Sun (1993) Woodroofe, M. and Sun, J. (1993).  A penalized maximum likelihood estimate of f (0+) when f is non-increasing.  Statistica Sinica, 3.
  * Yamanishi et al. (2004) Yamanishi, K., Takeuchi, J. I., Williams, G., and Milne, P. (2004).  On-line unsupervised outlier detection using finite mixtures with discounting learning algorithms.  Data Min. Knowl. Discov., 8, 275–300.

##  Appendix A Proofs for Section 4

Before proceeding any further, we introduce some new notations. We let
$l^{\infty}$ denote the collection of all bounded functions on $\mathbb{R}$,
which we equip with the uniform metric $\|\cdot\|_{\infty}$. For any function
$\mu:\mathbb{R}\mapsto\mathbb{R}$, we define the norm $\|\cdot\|_{a}^{b}$ by
$\|\mu\|_{a}^{b}=\sup_{x\in[a,b]}|\mu(x)|$. Also, we denote the boundary of a
set $A$ by $\text{bd}(A)$.

Suppose $(\Omega,\text{int}({\mathcal{H}}_{0}),\mathcal{P})$ is the common
probability space corresponding to the $X_{i}$’s and $Y_{j}$’s. For the rest
of this section, we let $\to_{p}$ and $\to_{a.s.}$ correspond to this
probability space. Since $F$ and $G$ are continuous, using the construction of
Section 1.1 of Shorack (1984) (see also p.93 of Shorack and Wellner, 2009), we
can show that there exist two independent Brownian bridges $\mathbb{V}_{1}$
and $\mathbb{V}_{2}$ on $(\Omega,\text{int}({\mathcal{H}}_{0}),\mathcal{P})$
such that

$\norm{\sqrt{m}(\mathbb{F}_{m}-F)-\mathbb{V}_{1}\circ
F}_{\infty}\to_{a.s.}0\quad\text{ as }m\to\infty,$ (A.1)

and

$\norm{\sqrt{n}(\mathbb{G}_{n}-G)-\mathbb{V}_{2}\circ
G}_{\infty}\to_{a.s.}0\quad\text{ as }n\to\infty.$ (A.2)

Let us denote

$\mathbb{U}=\lambda^{1/2}\mathbb{V}_{2}-(1-\lambda)^{1/2}\mathbb{V}_{1},$
(A.3)

where $\lambda$ is so that $m/N\to\lambda$. Note that $\mathbb{U}$ is also
distributed as a Brownian bridge. We will show that the asymptotic
distributions of our test statistics depend on $\mathbb{U}$. Also, for ours
$F$ and $G$,

$\|\mathbb{F}_{m}-F\|_{\infty}\to_{a.s.}0,\quad\|\mathbb{G}_{n}-G\|_{\infty}\to_{a.s.}0,\quad\|\mathbb{H}_{N}-H\|_{\infty}\to_{a.s.}0.$

In the sequel, we will also use the following fact on the convergence of the
quantiles of $\mathbb{H}_{N}$ and $D_{p,m,n}$, often without mentioning. As a
corollary to Fact 3.1A, it follows that
$\text{dist}(D_{p,m,n},D_{p})\to_{a.s.}0$ under Condition Condition N.

######  Fact 3.1.

Suppose Condition Condition N holds and $U$ is the open set such that
$D_{p}\subset U\subset\text{supp}(f)\cup\text{supp}(g)$. Further suppose
$p^{\prime}<p$ such that $D_{p^{\prime}}\subset U$. Then

  * A.

For any $t\in[p^{\prime},1-p^{\prime}]$,
$\mathbb{H}_{N}^{-1}(t)\to_{a.s.}H^{-1}(t)$.

  * B.

$D_{p,m,n}\subset D_{p^{\prime}}$ for all sufficiently large $m$ and $n$ with
probability one.

  * C.

For any $q>p$, $D_{q}\subset D_{p,m,n}$ for all sufficiently large $m$ and $n$
with probability one.

###### Proof of Fact 3.1.

Since $h=\lambda f+(1-\lambda)g>0$ on $U$, it follows that $H^{-1}$ is
continuous on $U$ (cf. Proposition A.7, pp. 98 of Bobkov and Ledoux, 2016).
Also, under Condition Condition N, $\|\mathbb{H}_{N}-H\|_{\infty}\to_{a.s.}0$.
Since $D_{p^{\prime}}\subset U$, $H^{-1}(t)\in U$ for any
$t\in[p^{\prime},1-p^{\prime}]$. Hence,
$\mathbb{H}_{N}^{-1}(t)\to_{a.s.}H^{-1}(t)$ by Lemma A5, pp. 96 of Bobkov and
Ledoux (2016). Hence part (A) of Fact 3.1 is proved.

Part (A) of Fact 3.1 implies $\mathbb{H}_{N}^{-1}(p)\to_{a.s.}H^{-1}(p)$ and
$\mathbb{H}_{N}^{-1}(1-p)\to_{a.s.}H^{-1}(1-p)$. To prove part B and C of Fact
3.1, therefore, it suffices to show that $H^{-1}(p^{\prime})<H^{-1}(p)$,
$H^{-1}(1-p^{\prime})>H^{-1}(1-p)$, $H^{-1}(p)<H^{-1}(q)$, and
$H^{-1}(1-q)<H^{-1}(1-p)$. Since $q<p<p^{\prime}$ and $D_{q}\subset
D_{p}\subset D_{p^{\prime}}\subset U$, it is enough to show that $H^{-1}$ is
continuous and strictly increasing on $U$ . The latter holds if $H$ s
continuous and strictly increasing on $U$. The proof now follows since we
already showed that $H$ has a positive density on $U$. ∎

####  A.0.1 Proof of Lemma 1

To prove Lemma 1, we will use an alternative definition of ${\mathcal{H}}_{0}$
and ${\mathcal{H}}_{1}$. Note that continuous distribution functions
$(F,G)\in{\mathcal{H}}_{0}$ if and only if $\sup_{x\in D_{p}}[F(x)-G(x)]\geq
0$. Because $D_{p}$ is compact, such a pair is in ${\mathcal{H}}_{1}$ if and
only if $\sup_{x\in D_{p}}[F(x)-G(x)]<0$. We will prove the current lemma in
some steps.
Step 1: closure of ${\mathcal{H}}_{0}$ Suppose
$(F_{n},G_{n})\in{\mathcal{H}}_{0}$ converges to
$(F,G)\in\mathcal{F}\times\mathcal{F}$ with respect to the metric $d_{2}$.
Then

$\displaystyle\sup_{x\in D_{p}}(F(x)-G(x))\geq-\sup_{x\in
D_{p}}(F_{n}(x)-F(x))$ $\displaystyle\ +\sup_{x\in
D_{p}}(F_{n}(x)-G_{n}(x))-\sup_{x\in D_{p}}(G(x)-G_{n}(x)),$

which is bounded below by $-\|F_{n}-F\|_{\infty}-\|G_{n}-G\|_{\infty}$. Taking
$n\to\infty$ yields $\sup_{x\in D_{p}}(F(x)-G(x))\geq 0$. Therefore,

$\text{cl}({\mathcal{H}}_{0})\subset\bigg{\\{}(F,G)\in\mathcal{F}\times\mathcal{F}:\sup_{x\in
D_{p}}(F(x)-G(x))\geq 0\bigg{\\}}={\mathcal{H}}_{0}.$

Hence, $\text{cl}({\mathcal{H}}_{0})={\mathcal{H}}_{0}$, i.e.
${\mathcal{H}}_{0}$ is closed in $\mathcal{F}\times\mathcal{F}$.

Step 2: interior and boundary of ${\mathcal{H}}_{0}$
If $(F,G)\in{\mathcal{H}}_{0}$, either $\sup_{x\in D_{p}}(F(x)-G(x))>0$ or
$\sup_{x\in D_{p}}(F(x)-G(x))=0$. Let us denote

$\mathcal{A}=\\{(F,G)\in\mathcal{F}\times\mathcal{F}:\sup_{x\in
D_{p}}(F(x)-G(x))=0\\}.$

By Lemma A.1, $\mathcal{A}\subset\text{bd}({\mathcal{H}}_{0})$ because this
set is not a part of the interior. Hence,

$\text{int}({\mathcal{H}}_{0})\subset\\{(F,G)\in\mathcal{F}\times\mathcal{F}:\sup_{x\in
D_{p}}(F(x)-G(x))>0\\}.$

Therefore, to prove

$\text{int}({\mathcal{H}}_{0})=\\{(F,G)\in\mathcal{F}\times\mathcal{F}:\sup_{x\in
D_{p}}(F(x)-G(x))>0\\},$ (A.4)

it suffices to prove

$\\{(F,G)\in\mathcal{F}\times\mathcal{F}:\sup_{x\in
D_{p}}(F(x)-G(x))>0\\}\subset\text{int}({\mathcal{H}}_{0}).$ (A.5)

Suppose $(F,G)\in{\mathcal{H}}_{0}$ satisfies $\sup_{x\in
D_{p}}(F(x)-G(x))>\epsilon$ for some $\epsilon>0$. Consider any
$(\tilde{F},\tilde{G})$ such that
$d_{2}((\tilde{F},\tilde{G}),(F,G))<\epsilon/3$, which means
$\|\tilde{F}-F\|_{\infty}<\epsilon/3$ and
$\|\tilde{G}-G\|_{\infty}<\epsilon/3$. Since $F-G$ is continuous and $D_{p}$
is compact, the supremum of $F-G$ over $D_{p}$ is attained at some $z\in
D_{p}$. We have

$\displaystyle\tilde{F}(z)-\tilde{G}(z)\geq
F(z)-G(z)-\|\tilde{F}-F\|_{\infty}-\|\tilde{G}-G\|_{\infty},$

which is greater than $\epsilon/3$. Thus $\sup_{x\in
D_{p}}(\tilde{F}(x)-\tilde{G}(x))>0$ which implies
$(F,G)\in\text{int}({\mathcal{H}}_{0})$. Therefore, we have established (A.4).
Because $\text{cl}({\mathcal{H}}_{0})={\mathcal{H}}_{0}$, (A.4) also implies
$\text{bd}({\mathcal{H}}_{0})=\mathcal{A}$.

Step 3: boundary of ${\mathcal{H}}_{1}$
Our first step is to show

$\displaystyle\text{cl}({\mathcal{H}}_{1})\subset\bigg{\\{}(F,G)\in\mathcal{F}\times\mathcal{F}:\sup_{x\in
D_{p}}(F(x)-G(x))\leq 0\bigg{\\}}.$ (A.6)

To this end, consider $(F_{n},G_{n})\in{\mathcal{H}}_{1}$ converging to
$(F,G)\in\mathcal{F}\times\mathcal{F}$ in $d_{2}$. Then

$\displaystyle\sup_{x\in D_{p}}(F(x)-G(x))\leq\sup_{x\in
D_{p}}(F(x)-F_{n}(x))+\sup_{x\in D_{p}}(F_{n}(x)-G_{n}(x))+\sup_{x\in
D_{p}}(G_{n}(x)-G(x)),$

which is bounded above by $\|F_{n}-F\|_{\infty}+\|G_{n}-G\|_{\infty}$. Letting
$n\to\infty$, we obtain $\sup_{x\in D_{p}}(F(x)-G(x))\leq 0$ implying
$(F,G)\in{\mathcal{H}}_{1}$. Hence, (A.6) holds, which implies
$\text{cl}({\mathcal{H}}_{1})\subset\mathcal{A}\cup{\mathcal{H}}_{1}$. Using
Lemma A.1 we obtain

$\text{cl}({\mathcal{H}}_{1})\subset\mathcal{A}\cup{\mathcal{H}}_{1}\subset\text{bd}({\mathcal{H}}_{1})\cup{\mathcal{H}}_{1}=\text{cl}({\mathcal{H}}_{1}).$

Hence, the inclusion in (A.6) is actually an equality and
$\mathcal{A}\cup{\mathcal{H}}_{1}=\text{bd}({\mathcal{H}}_{1})\cup{\mathcal{H}}_{1}$.
The proof will be complete if we can show that
$\text{int}({\mathcal{H}}_{1})={\mathcal{H}}_{1}$, because then
$\text{bd}({\mathcal{H}}_{1})=\mathcal{A}$ follows. To that end, consider
$(F,G)\in{\mathcal{H}}_{1}$. Suppose $\sup_{x\in D_{p}}(F(x)-G(x))<-\epsilon$.
Let $(\tilde{F},\tilde{G})\in\mathcal{F}\times\mathcal{F}$ be such that
$d_{2}((\tilde{F},\tilde{G}),(F,G))<\epsilon/3$. Then

$\displaystyle\sup_{x\in D_{p}}(\tilde{F}(x)-\tilde{G}(x))\leq\sup_{x\in
D_{p}}(F(x)-G(x))+\|\tilde{F}-F\|_{\infty}+\|\tilde{G}-G\|_{\infty},$

which is less than $-\epsilon/3$. Therefore,
$(\tilde{F},\tilde{G})\in{\mathcal{H}}_{1}$, which completes the proof. $\Box$

######  Lemma A.1.

Under the set-up of Lemma 1,

$\mathcal{A}:=\bigg{\\{}(F,G)\in\mathcal{F}\times\mathcal{F}:\sup_{x\in
D_{p}}(F(x)-G(x))=0\bigg{\\}}\subset\text{bd}({\mathcal{H}}_{1}).$

###### Proof.

Consider $(F,G)\in\mathcal{A}$. Because
$\mathcal{A}\cap{\mathcal{H}}_{1}=\emptyset$, it suffices to show that $(F,G)$
is a limit point of ${\mathcal{H}}_{1}$. We will show that given $\delta>0$,
there exists $(\tilde{F},\tilde{G})\in{\mathcal{H}}_{1}$ so that
$d_{2}((F,G),(\tilde{F},\tilde{G}))<\delta$.

Since $F-G$ is continuous, its supremum over $D_{p}$ is attained. Thus, there
exists $C_{p}\subset D_{p}$ so that $F=G$ on $C_{p}$. Because $p\leq\lambda
F+(1-\lambda)G\leq 1-p$ on $D_{p}$, $\inf_{z\in C_{p}}F(z)\geq p$ and
$\sup_{z\in C_{p}}F(z)\leq 1-p$. Suppose $a^{\prime}=\inf C_{p}$ and
$\beta^{\prime}=\sup C_{p}$. Since $F$ is continuous, there exists closed
interval $L_{p}=[a,\beta]\supset C_{p}$ such that $L_{p}\subset D_{p}$ and
$p/2<F<1-p/2$ on $[a,\beta]$. We choose $a$ and $\beta$ so that additionally
the followings hold:

  1. 1.

If $a^{\prime}>H^{-1}(p)$ then $a<a^{\prime}$. Thus $a=a^{\prime}$ only if
$a^{\prime}=H^{-1}(p)$, in which case, the only choice for $a$ is $a^{\prime}$
because $[a^{\prime},\beta^{\prime}]\subset[a,\beta]\subset D_{p}$.

  2. 2.

If $\beta^{\prime}<H^{-1}(1-p)$ then $\beta^{\prime}<\beta$. When
$\beta^{\prime}=H^{-1}(1-p)$, the only choice for $\beta$ is $\beta^{\prime}$
because $[a^{\prime},\beta^{\prime}]\subset[a,\beta]\subset D_{p}$.

Letting $\epsilon^{\prime}=\min(\delta/2,p/4)$,
$t_{1}=F^{-1}(1-\epsilon^{\prime})$ and $t_{2}=t_{1}+1$, we define

$\tilde{F}(x)=\begin{cases}\max(0,F(x)-\epsilon^{\prime})&x<a\\\
F(x)-\epsilon^{\prime}&x\in[a,t_{1}]\\\
(x-t_{1})\frac{F(t_{2})-\tilde{F}(t_{1})}{t_{2}-t_{1}}+\tilde{F}(t_{1})&x\in(t_{1},t_{2}]\\\
F(x)&x>t_{2}\end{cases}$ (A.7)

That $\tilde{F}$ is a distribution function is clear from the definition. We
claim that (a) $\tilde{F}$ is continuous, (b)
$\|\tilde{F}-F\|_{\infty}<\delta$, and finally (c)
$(\tilde{F},G)\in{\mathcal{H}}_{1}$. Taking $\tilde{G}=G$, the proof of the
current lemma follows. Hence, it remains to prove Claim (a), (b), and (c).

First, we prove Claim (a). Because $F$ is continuous,
$\max(0,F(x)-\epsilon^{\prime})$ is continuous in $x$ on $(-\infty,a)$. BSince
$F(a)>p/2>\epsilon^{\prime}$, continuity of $F$ also implies
$\tilde{F}(a-)=F(a)-\epsilon^{\prime}=F(a+)$. Therefore, $\tilde{F}$ is
continuous on $(-\infty,a]$. Therefore (A.7) implies $\tilde{F}$ is continuous
on $(-\infty,a]$. The continuity of $F(x)-\epsilon^{\prime}$ on $(a,t_{1}]$
follows because $F$ is continuous. $\tilde{F}$ is linear on $(t_{1},t_{2}]$,
and equals $F$ on $(t_{2},\infty)$ Moreover, the left and right limits of
$\tilde{F}$ agree at $t_{1}$ and $t_{2}$. Therefore $\tilde{F}$ is continuous.

To prove Claim (b), it suffices to show that
$|\tilde{F}(x)-F(x)|\leq\epsilon^{\prime}$ for any
$x\in(-\infty,a)\cup(t_{1},t_{2}]$. Note that if $x<a$ and
$F(x)\leq\epsilon^{\prime}$, then $\tilde{F}(x)=0$. Thus,
$\tilde{F}(x)-F(x)\leq\epsilon^{\prime}$ in this case. If $x<a$ but
$F(x)>\epsilon^{\prime}$, $\tilde{F}(x)-F(x)=\epsilon^{\prime}$. Therefore,
$|\tilde{F}(x)-F(x)|\leq\epsilon^{\prime}$ for $x\in(-\infty,a)$. On the other
hand, since $t_{1}=F^{-1}(1-\epsilon^{\prime})$, $F(x)$ is greater than
$1-\epsilon^{\prime}$ on $(t_{1},t_{2}]$. Also for $x\in[t_{1},t_{2}]$,

$\displaystyle\tilde{F}(x)\geq\tilde{F}(t_{1})=F(t_{1})-\epsilon^{\prime}=F(F^{-1}(1-\epsilon^{\prime}))-\epsilon^{\prime}=1-2\epsilon^{\prime}.$

Thus on $[t_{1},t_{2}]$, $|\tilde{F}(x)-F(x)|$ is bounded by
$2\epsilon^{\prime}$, which is not greater than $\delta$ by our choice of
$\epsilon^{\prime}$. Hence, we have shown that
$\|\tilde{F}-F\|_{\infty}<\delta$.

To prove Claim (c), first note that $\tilde{F}\leq F$ on $\mathbb{R}$. Let us
partition $D_{p}=[H^{-1}(p),a)\cup[a,\beta]\cup(\beta,H^{-1}(1-p)]$. If
$a=H^{-1}(p)$, then we define $[H^{-1}(p),a)$ to be the empty set. Similarly,
$(\beta,H^{-1}(1-p)]$ is non-empty only if $\beta<H^{-1}(1-p)$. Note that
$\beta\leq t_{1}$, implying $[a,\beta]\subset[a,t_{1}]$. Therefore, by (A.7),

$\sup_{x\in[a,\beta]}(\tilde{F}(x)-G(x))=\sup_{x\in[a,\beta]}(F(x)-G(x))-\epsilon^{\prime}\leq\sup_{x\in
D_{p}}(F(x)-G(x))-\epsilon^{\prime},$

which equals $-\epsilon^{\prime}$. Now let us consider the set
$[H^{-1}(p),a)$. Of course if $a=H^{-1}(p)$, there is nothing to prove. So
suppose $a>H^{-1}(p)$. Because $C_{p}\cap[H^{-1}(p),a)=\emptyset$, there is no
$x\in[H^{-1}(p),a)$ such that $F(x)-G(x)=0$. Therefore, if
$\sup_{x\in[H^{-1}(p),a)}(F(x)-G(x))=0$, the set $[H^{-1}(p),a]$ must contain
a limit point of $C_{p}$. However,
$\text{cl}(C_{p})\subset[a^{\prime},\beta^{\prime}]$ where $a<a^{\prime}$
because $a>H^{-1}(p)$. Therefore, $[H^{-1}(p),a]$ can not contain any limit
point of $C_{p}$ either. Thus we must have
$\sup_{x\in[H^{-1}(p),a)}(F(x)-G(x))<0$. Therefore,
$\sup_{x\in[H^{-1}(p),a)}(\tilde{F}(x)-G(x))<0$ as well because $\tilde{F}\leq
F$. On the other hand, using $H^{-1}(1-p)<t_{1}$, we have

$\sup_{x\in(\beta,H^{-1}(1-p)]}(\tilde{F}(x)-G(x))=\sup_{x\in(\beta,H^{-1}(1-p)]}(F(x)-G(x))-\epsilon\leq\sup_{x\in
D_{p}}(F(x)-G(x))-\epsilon^{\prime},$

which equals $-\epsilon^{\prime}$. Combining the above pieces, we obtain
$\sup_{x\in D_{p}}(\tilde{F}(x)-G(x))<0$, which completes the proof of
$(\tilde{F},G)\in{\mathcal{H}}_{1}$.

∎

#### Proof of Lemma 2

Since $F,G\in\text{int}({\mathcal{H}}_{0})$, there exists $x^{\prime}_{0}\in
D_{p}$ and $\delta>0$ such that
$G(x^{\prime}_{0})-F(x^{\prime}_{0})<-3\delta$. Because $F$ and $G$ are
continuous, we can find $x_{0}\in\text{int}(D_{p})$ such that
$G(x_{0})-F(x_{0})<-2\delta$. Hence, (A) indicates that with probability one,
$\mathbb{G}_{n}(x_{0})-\mathbb{F}_{m}(x_{0})<-\delta$ for all sufficiently
large $m$ and $n$, which leads to

$\dfrac{\mathbb{G}_{n}(x_{0})-\mathbb{F}_{m}(x_{0})}{\sqrt{\dfrac{\mathbb{F}_{m}(x_{0})(1-\mathbb{F}_{m}(x_{0}))}{m}+\dfrac{\mathbb{G}_{n}(x_{0})(1-\mathbb{G}_{n}(x_{0}))}{n}}}\leq\dfrac{-\delta}{\sqrt{1/m+1/n}},$