A Nearest-Neighbor Based Nonparametric Test for Viral Remodeling in Heterogeneous Single-Cell Proteomic Data

An important problem in contemporary immunology studies based on single-cell protein expression data is to determine whether cellular expressions are remodeled post infection by a pathogen. One natural approach for detecting such changes is to use non-parametric two-sample statistical tests. However, in single-cell studies, direct application of these tests is often inadequate because single-cell level expression data from uninfected population often contains attributes of several latent sub-populations with highly heterogeneous characteristics. As a result, viruses often infect these different sub-populations at different rates in which case the traditional nonparametric two-sample tests for checking similarity in distributions are no longer conservative. We propose a new nonparametric method for Testing Remodeling Under Heterogeneity (TRUH) that can accurately detect changes in the infected samples compared to possibly heterogeneous uninfected samples. Our testing framework is based on composite nulls and is designed to allow the null model to encompass the possibility that the infected samples, though unaltered by the virus, might be dominantly arising from under-represented sub-populations in the baseline data. The TRUH statistic, which uses nearest neighbor projections of the infected samples into the baseline uninfected population, is calibrated using a novel bootstrap algorithm. We demonstrate the non-asymptotic performance of the test via simulation experiments and derive the large sample limit of the test statistic, which provides theoretical support towards consistent asymptotic calibration of the test. We use the TRUH statistic for studying remodeling in tonsillar T cells under different types of HIV infection and find that unlike traditional tests, TRUH based statistical inference conforms to the biologically validated immunological theories on HIV infection.


Introduction
In many contemporary scientific methodologies, it is extremely difficult, even in well-regulated laboratory experiments, to simultaneously control the multitude of factors that give rise to heterogeneity in the population (Chapter 3 of Holmes and Huber (2018)). Nevertheless, these experiments are very powerful, and are often our only recourse to study several interesting biological phenomena. For example, in singlecell proteomic and genomic studies (Jia et al. (2017), Jiang et al. (2018), Shi and Huang (2017), Wang et al. (2018)), it is now well understood that there is high heterogeneity in cellular responses from controlled cell population. Statistical tests are often used on these datasets to determine differences between the case and control samples. The presence of heterogeneity greatly complicates statistical inference, and direct application of existing two-sample testing methods, without modulating for the latent heterogeneity in the samples, may lead to erroneous statistical decisions and scientific consequences. The problem of testing similarity in the distributions of two samples under heterogeneity arises in a host of modern immunology research set-ups where heterogeneous protein expression datasets collected at single-cell resolution are analyzed to detect viral perturbation. To provide a rigorous statistical hypothesis testing framework for these immunology studies, we develop a new non-parametric testing procedure based on nearest-neighbor distances, that can accurately detect if there are differences between the case and control samples in the presence of unknown heterogeneity in the data-generation process. We next provide the background of the problem through an immunology study on human immunodeficiency virus (HIV) infection in tonsillar cells.

Phenotypic Profiling of T Cells Under HIV Infection
In single-cell immunology, phenotypic profiling of immune cells under the influence of a target virus, such as the HIV (Cavrois et al., 2017), the varicella zoster virus (VZV) (Sen et al., 2014), or the rotavirus (RV) (Sen et al., 2012), is a critical research endeavor. It enhances understanding of which subsets of cells are most or least susceptible to infection, leading to new insights regarding the magnitude of viral persistence, which is crucial in the development of life saving drugs (Sen et al., 2015). Mass cytometry based techniques (Bendall et al., 2011, Giesen et al., 2014 are popularly used for generating proteomic datasets for such phenotypic analysis. These techniques can simultaneously measure around fifty protein expressions on individual cells. In this paper we provide a rigorous statistical analysis for testing if there are any HIV induced changes in the proteomic expressions of tonsillar T cells, which are a type of lymphocyte that plays a central role in the immune response, based on the dataset generated in (Cavrois et al., 2017). technique. Tonsillar T cells from 4 healthy donors were infected with two variants of a HIV viral strain: Nef rich HIV and Nef deficient HIV. Nef (Negative Regulatory Factor) is a protein encoded by HIV which enhances virus replication in the host cell by protecting infected cells from immune surveillance. We study the differential impact of these two variants on the immune cells. The healthy cells were cultured, processed and batched into three identically distributed populations for each donor. For each patient, one among the three batches were randomly selected and phenotyped to generate the expression data of the uninfected population, while the other two batches were contaminated with the Nef rich

Viral Remodeling
If the virus changes the expression of any of the surface markers, which are proteins attached to the cell membrane, then the cell is said to have undergone viral remodeling of its phenotypic characteristics (Sen et al., 2014). A virally remodelled cell will have aberrant inter-cellular activities, therefore, detecting the presence of remodeling is a fundamental step towards understanding the mechanism of pathogenesis and disease progression. Detecting remodeling translates to testing if there is enough evidence in the data to reject the null hypothesis that the joint distribution of all the surface proteins is same between the uninfected and virus infected sample. A natural approach for this problem is to invoke non-parametric two-sample testing methods to see if there is enough evidence to support the alternative hypothesis that the virus has changed the distribution of least one of the sub-populations. However, for single-cell level expression data, the hypothesis test described above is particularly difficult because of the following two reasons: (a) the presence of heterogeneity in the uninfected population, and (b) due to the phenomenon of preferential infection. Single-cell resolution expression data from processed uninfected population often contains attributes from several latent sub-populations with highly heterogeneous characteristics. This sub-population level heterogeneity in the uninfected (also referred to as the control or baseline) samples can arise from varied attributes that cannot be controlled in experiments, such as differences in the cell effector functions, trafficking and longevity (Cavrois et al., 2017). Viruses often infect these different sub-populations at different rates. If a virus infects different sub-population at different rates, but does not alter the marker expressions for any of the sub-populations, still the distribution of the overall viral sample will be different from the uninfected samples. In these situations, the difference in distribution between the infected and the uninfected samples is not due to viral remodeling but due to preferential infection (for a detailed biological explanation see Figures 2A and 2B of (Cavrois et al., 2017)) of the uninfected sub-populations by the virus. Figure 2 presents two scenarios that may arise when the cloud of infected and uninfected cells are analyzed with respect to a single marker A. In this toy example, Panel 1 in Figure 2 shows that the uninfected T cells arise from three sub-populations with varying expression levels for marker A which may reflect their inherent heterogeneity with respect to cell longevity. The scenario of preferential infection is depicted in Panel 2 where the HIV preferentially infects the T cell subpopulation that has a lower expression level for marker A amongst the uninfected cells. Moreover, the virus does not alter the expression levels of these infected cells when compared to Panel 1. In Panel 3, which represents HIV remodeling, the virus targets those uninfected cells that have low to medium expression for marker A amongst the uninfected cells and alters their original expression levels upon infection. The distinct pink and yellowish shade of the infected cells in panel 3 depicts their phenotypic change associated with infection. Here, we have described the phenomenon of viral remodeling only for the HIV. However, remodeling analysis is widely conducted across virology for understanding mechanism of other pathogens also. For correct scientific understanding of the viral mechanism, it is extremely important to accurately detect the instances of viral remodeling from mere preferential infection. However, popular single-cell based segmentation and classification algorithms (Amir et al., 2013, Bruggner et al., 2014, Linderman et al., 2012, Qiu, 2012) lack a rigorous statistical hypothesis testing framework for conducting two-sample inference and can greatly suffer in testing problems, particularly if there is high imbalance in the sizes of the uninfected (control) Figure 2: Schematic representation of HIV remodeling of T cells with respect to a single marker A. Panel 1 shows that the uninfected T cells arise from three sub-populations with varying expression levels for marker A. Panel 2 depicts preferential infection where the HIV preferentially infects the T cell subpopulation that has a lower expression level for marker A amongst the uninfected cells and the infection does not alter the expression levels of the T cells when compared to Panel 1. Panel 3 represents HIV remodeling where the HIV targets those uninfected cells that have low to medium expression for marker A amongst the uninfected cells and alters their original expression levels upon infection. and infected (case) samples, which is often the situation in virology.

Testing Procedures in Existing Literature and Statistical Challenges
The statistical framework for testing remodeling falls under the realm of non-parametric two-sample testing. For univariate data, non-parametric two-sample tests like the Kolmogorov-Smirnov test, the Wilcoxon rank-sum test, and the Wald-Wolfowitz runs test are extremely popular and find a place in every practitioner's toolkit. Multidimensional versions of these widely used tests date back to the randomization tests of Chung and Fraser (1958) and to the generalized Kolmogorov-Smirnov test of Bickel (1969). Friedman and Rafsky (Friedman and Rafsky, 1979) proposed the first computationally efficient non-parametric 2-sample test, which applies to high-dimensional data. The Friedman-Rafsky edgecount test, which can be viewed as a generalization of the univariate runs test, computes the Euclidean minimal spanning tree (MST) 1 of the pooled sample, and rejects the null if the number of edges with endpoints in different samples is small. Many variants of the edgecount test, based on nearest-neighbor distances and geometric graphs have been proposed over the years Hall and Tajvidi (2002), Henze (1984), Rosenbaum (2005), Schilling (1986), Weiss (1960). Recently, Chen and Friedman (2017) suggested novel modifications of the edge-count test for high-dimensional and object data, and Chen et al. (2018) proposed new and powerful tests to deal with the issue of sample size imbalance. Asymptotic properties of two-sample tests based on geometric graphs can be studied in the general framework described in Bhattacharya (2019). Other popular two-sample tests include the test of Baringhaus and Franz (2004), the energy distance test of Aslan and Zech (2005), and the kernel based test using maximum mean discrepancy of Gretton et al. (2007). More recently, Chen et al. (2013) address the problem of sample size imbalances in the two-sample problem by constructing an ensemble subsampling scheme for the nearestneighbor tests Henze (1984), Schilling (1986). Very recently, Deb and Sen (2019) and Ghosal and Sen (2019) proposed distribution-free two-sample tests based on the concept of multivariate ranks, defined using optimal transport.
One of the main challenges for devising a statistically correct test to detect viral remodeling from preferential infection is that the virus may infect different sub-populations at different rates. In Section 2, we show that even in very large sample sizes direct application of existing non-parametric two-sample tests can lead to erroneous inference. We expound this phenomenon by exhibiting explicit scenarios of preferential infection and remodeling where traditional tests fail in a simple setting of d = 2 markers.
In Figure 3, the green triangles correspond to a sample of uninfected (UI) cells that arise from three different subpopulations while the red dots reflect the infected (VI) cells. The leftmost panel presents a setting where the virus has infected all the three cellular subpopulations and the overlap of the UI and VI cells indicate no remodeling. The middle panel presents a scenario where the cells have undergone remodeling under the influence of the virus as is evident through a shift in the location of the VI cells.
The rightmost panel reflect no remodeling but preferential infection. The different g-tests , Chen and Friedman, 2017, Friedman and Rafsky, 1979, the cross-match test (Rosenbaum, 2005), and the energy test (Aslan and Zech, 2005) reject the null hypothesis of no remodeling in all the three cases, in each of the 100 simulation replications (see Table 3 in Section 3.1). This is not surprising, because these tests are designed to test the simple null hypothesis of equality of the two distributions.
Due to the presence of sub-population level heterogeneity the problem of testing for remodeling warrants testing a composite null hypothesis. To this end, note that under preferential infection, the two-samples arise from the mixture distribution with identical component distributions but with different mixing weights. This is the case for the right most sub-plot in Figure 3. In this paper, we formulate the problem of testing for preferential infection versus remodeling as a composite two-sample hypothe- sis with mixture distributions, and develop a new nearest-neighbor based test that can consistently and efficiently detect the differences between the two-samples.

The TRUH Testing Framework: Novel attributes and Our Contributions
In this article, we propose a novel procedure for Testing Remodeling Under Heterogeneity (TRUH), that effectively incorporates the underlying heterogeneity and imbalance in the samples, and provides a conservative test for the composite null hypothesis that the two-samples arise from the same mixture distribution but may differ with respect to the mixing weights. We summarize its key attributes below: • The TRUH statistic is based on a nearest neighbor approach (Cover andHart, 1967, Devroye et al., 2013) that relies on first identifying for every infected cell a predictive precursor cell which is the phenotypically closest cell in the uninfected population. It then measures the relative dissimilarities between the infected cells and their predictive precursors and the predictive precursors to their most phenotypically similar uninfected cells. A large relative dissimilarity between the infected cells and their predictive precursors indicates surface protein regulation or remodeling by the virus, while a small relative dissimilarity provides evidence for preferential infection or no remodeling.
• We describe an efficient bootstrap based approach for calibrating the TRUH test statistic, and evaluate its performance in finite-sample simulations. We then use this method to test for viral remodeling in tonsillar T cells under different types of HIV infection, corroborating the efficacy of our proposed procedure.
• We provide an extensive theoretical understanding of the large sample characteristics of our proposed test statistics. We establish the L 2 -limit of our proposed statistic using asymptotic properties of functionals of random geometric graphs Penrose and Yukich (2003). The limit can be expressed in terms of the densities of the uninfected and infected populations and dimension dependent constants obtained from nearest-neighbor distances defined on a homogeneous Poisson process. Using these properties, we can select a cut-off for the TRUH statistic that is asymptotically consistent against biologically-relevant location alternatives. Traditional non-parametric tests enjoy these consistency properties in homogeneous populations but not under heterogeneity. We show that using a nearest neighbor based approach this inefficiency of existing non-parametric tests in heterogeneous data can be mitigated.
The rest of the paper is organized as follows: In Section 2, we formulate the problem of testing for remodeling in single-cell virology as a heterogeneous two-sample problem, describe the TRUH framework, and show how it can be calibrated using the bootstrap. Numerical experiments demonstrating the non-asymptotic performance of our testing procedure are given in Section 3. In Section 4 we use TRUH for studying remodeling in tonsillar T cells under different types of HIV infection. The asymptotic properties of the test statistic are discussed in Section 5. We conclude the paper in Section 6 with a discussion.
The technical details and proofs of the theoretical results are given in the supplementary materials.

Statistical Framework and the Proposed TRUH Statistic
In this section we formulate the problem of testing for remodeling in single-cell virology as a heterogeneous two-sample problem (Section 2.1), introduce the TRUH statistic (Section 2.2), and discuss how to calibrate it using the bootstrap (Section 2.3).

The Heterogeneous Two-Sample Problem
In our virology example, the baseline constitutes the m uninfected cells. For each cell, i ∈ {1, . . . m}, we denote by U i a d-dimensional vector of cellular characteristics typically measuring expressions corresponding to different genes or proteins. Denote the uninfected/baseline population by U m = {U 1 , . . . , U m }.
Let F 0 be the cumulative distribution function (cdf) of the baseline population with the heterogeneity in the population being reflected by K different sub-groups each having unimodal distributions with distinct modes and cdfs F 1 , . . . , F K , and mixing proportions w 1 , . . . , w K , such that where w a ∈ (0, 1) and K a=1 w a = 1.
Note that, the number of components K, the mixing distributions F 1 , . . . , F K , and the mixing weights w 1 , . . . , w K are fixed (non-random) attributes, which are unknown. Also, as F 1 , . . . , F K are cdfs from unimodal distributions with distinct modes, F 0 is well-defined with a unique specification. In addition to the uninfected population, we observe n i.i.d. infected observations V n = {V 1 , . . . , V n } from a distribution function G in R d . Note that, the infected and uninfected samples U m and V n are collected from separate experiments and are independent of each other.
Simple versus Composite Null In single-cell virology when an uninfected population is exposed to a pathogen, the virus may infect the different sub-populations at different rates. Therefore, even if the virus does not cause any change in the cellular characteristics, the virus infected sample might have different representations of the uninfected sub-populations than the uninfected mixing proportions {w 1 , . . . , w K }.
As such, it is quite possible that a few of uninfected sub-populations are completely absent in the viral population, which biologically implies that the virus preferentially targets few cellular sub-populations.
Thus, if the virus does not induce any change in the cellular characteristics, then the distribution of the infected population G lies in a class of distributions F(F 0 ) that contains any convex combination of {F 1 , . . . , F K } including the boundaries, that is, Note that the uninfected cdf F 0 is a particular member of the class F(F 0 ). If the virus induces changes in the cellular characteristics, then the viral population distribution would contain at least one non-trivial sub-population with distribution substantially different from {F 1 , F 2 , . . . , F K } or their linear combinations. Thus, the test for viral remodeling tantamounts to testing the following composite null hypothesis: If the null hypothesis is accepted, we say the virus exhibits preferential infection, otherwise we say the virus exhibits remodeling (see Figure 6 below), and the hypothesis testing problem (3) will be referred to as the problem of testing remodeling under heterogeneity (TRUH). Later on, to facilitate proofs of the theoretical properties of our proposed method, we will assume that the baseline cdfs F 1 , . . . , F K have unimodal densities f 1 , . . . , f K (with respect to Lebesgue measure). In this case, the baseline uninfected population will have density f 0 = K a=1 w a f a , and the set of distributions in (2) can be represented in terms of the densities f 1 , . . . , f K , and will be denoted by F(f 0 ).

Inefficiency of Existing Tests
Traditional non-parametric graph-based two-sample tests, such as the edgecount (EC) test of Friedman and Rafsky (1979) or the crossmatch (CM) test of Rosenbaum (2005), are tailored for the null hypothesis H 0 : F 0 = G, that is, testing whether the distributions of the uninfected samples U m and the infected samples V n are the same. However, not surprisingly, direct application of these tests to the composite hypothesis testing problem described in (3) above gives nonconservative procedures. To see this, consider the EC test. Recall that the EC test is based on the statistic R(U m , V n ) which counts the number of edges in the minimal spanning tree (MST) of the pooled sample {U 1 , . . . , U m , V 1 , . . . , V n } that connect points from different samples. Then, the EC test rejects the null Henze and Penrose (1999) in the usual limiting regime where m, n → ∞ and n/m → ρ ∈ (0, ∞). In particular, it follows from Theorem 1 of Henze and Penrose (1999) that , and δ d is a constant depending only on dimension d. More precisely, δ d is the variance of the degree of the origin 0 ∈ R d in the minimal spanning tree built on a homogeneous Poisson process of rate 1 in R d with the origin added to it. Note that (4) shows that the test with rejection region {R m,n (U m , V n ) < C m,n (α)} is asymptotically level α for the null hypothesis The following proposition shows that direct application of the EC test as described above, will not be conservative for testing the hypothesis (3) of viral remodeling. In fact, for cases of preferential infection but no remodeling the EC test will produce undesired false discoveries.
Proposition 1. Fix α ∈ (0, 1/2). Then for F 0 as in (1) and for any G ∈ F(F 0 ) \ {F 0 } in the usual limiting regime, from g, where f 0 and g are the densities (with respect to the Lebesgue measure) of F 0 and G, respectively.
The proof of the above result is given in the supplementary materials (Section A). This shows that for any level α, the EC test will be asymptotically inconsistent as it would reject with certainty all cases of preferential infection but no remodeling. This phenomenon is demonstrated in Figure 4 through a simple univariate simulation experiment. Here, we consider m = 1000, n = 50, and d = 1. The true distribution of the uninfected and infected sub-populations are gaussian mixtures. We consider two cases: In the top row, we describe the density of the true uninfected F0 (in continuous line) and the density of the infected G (in dotted line) for the two cases. In both cases, F0 and G are mixtures of normal distributions. In the first case, all the three equiprobable sub-populations in F0 have undergone a discernible location change in G. In case B, F0 again has three equiprobable sub-groups while G has two of those three sub-groups. Thus, while case A signifies viral remodeling, there is no remodeling but only preferential infection in Case B. In the bottom row, we have the histogram of the values of the TRUH statistic in the left (defined below in (7)) and the edgecount statistic in the right, respectively, under the two cases.
• Case A: Here, F 0 and G are mixtures of 3 equi-likely Gaussians, with each sub-population in . In this case, there is preferential infection, but no remodeling, that is, G ∈ F(F 0 ) with the middle population in F 0 being resistant to viral infection.
Any test for the hypothesis (3) should ideally reject Case A and fail to reject Case B. However, Figure 4 shows that the histogram of EC test statistic values across 500 replications under cases A and B have a significant overlap. Table 1 shows the rejection rate (proportion of false discoveries) in Case B and power  our proposed test statistic (TRUH), described in the following section, entertains possibilities where both the rejection rate and the power attain the desired limit.

Proposed Test Statistic: TRUH
In this section we describe a nearest-neighbor based statistic for testing the hypothesis of remodeling.
To this end, recall that the Euclidean distance of V i to its nearest point in the uninfected sample U m . The point in U m which attains this minimum will be denoted by N (V i , U m ) 3 and constitutes a key quantity in measuring the relative phenotypic difference between the infected cells and their closest uninfected counterparts. In which is the distance of N (V i , U m ) to its nearest neighbor in U m . Our proposed test statistic for testing (3), hereafter referred to as the TRUH statistic, is whereD m,n = 1 n n i=1 D i andC m,n = 1 n n i=1 C i . Note that the TRUH statistic above is an aggregated measure of how far apart each viral cell is from the uninfected sample compared to the local distance between uninfected sample points in its vicinity. Consider, for example, panel A in Figure 6 that represents a schematic for remodeling, while panel B depicts preferential infection. Here, the three infected cells (in red) in Panel A are phenotypically different than their uninfected counterparts and thus the average gap |D m,n −C m,n | in Panel A, averaged over the three infected cells, is relatively larger than what is observed under preferential infection in Panel B. Therefore, we develop a test to reject the null hypothesis of no remodeling for large values of T m,n . The cut-off for T m,n can be chosen based on a bootstrap calibration (Section 2.3) or using the asymptotic limit of T m,n (Section 5). Note that, since the nearest neighbor of a point, in a cloud of n random points in R d , typical lies within a ball of radius n − 1 d centered at that point, the TRUH statistic is scaled by n 1− 1 d , which makes T m,n bounded in probability.
One of the interesting properties of the quantity T m,n is that it only involves enumeration of distance based features for the viral sample, unlike classical graph-based two-sample tests Rafsky, 1979, Rosenbaum, 2005) which are built using the inter-point distances of the pooled sample. As a consequence, the TRUH test statistic is not symmetric in its usage of the uninfected and infected samples, even when the sample sizes are equal and the two-samples were actually generated from the same population distribution. This asymmetric sample usage of TRUH helps in tackling possibly different heterogeneity levels in the two-samples. Finally, note that even though the quantities D i and C i are defined above using the Euclidean distance, they can be easily generalized to any arbitrary distance function, and the statistic T m,n can potentially be used in non-Euclidean data spaces, such as graph data or functional data, as well.

Bootstrap Based Calibration for TRUH
In this section, we present a bootstrap based procedure to determine the cut-off t m,n,α for a level α test using T m,n . To this end, recall that F(F 0 ) contains any convex combination of the baseline distribution functions {F 1 , . . . , F K }. Therefore, the proposed bootstrap procedure relies on random sampling of the mixing proportions a large number of times to construct surrogate samples from F(F 0 ), and then computing the TRUH statistic in each case to generate a pseudo null distribution which is finally used to estimate the level α cut-offs.
Our algorithm leverages the fact that in our virology example the number m of uninfected samples is much larger than the size n of the infected samples. Therefore, we can use the prediction strength approach of Tibshirani and Walther (2005) on the uninfected samples to obtain an estimateK of the unknown number of heterogeneous subgroups K. We then use this value ofK to estimate the class memberships of the baseline samples U m using aK-means algorithm. For 1 ≤ a ≤K, denote bŷ J a ⊆ {1, 2, . . . , m} the subset of indices which belong to class a in the output of theK-means algorithm.
Let UĴ a = {U i : i ∈Ĵ a } be the subset of the baseline samples estimated to be in the a-th class by thê K-means algorithm. Note that U m = {UĴ a : a = 1, 2, . . . ,K} and K }, we construct a surrogate infected sample from F(F 0 ) as follows: for 1 ≤ a ≤K, randomly sample nλ (b) a elements without replacement from UĴ a . Denote the chosen elements by and set the remaining m a − nλ Note that under the null hypothesis of no remodeling (G ∈ F(F 0 )), the bootstrapped samples in the b-th (which are surrogates for U m and V n , respectively), can be used to compute the This quantity is the surrogate of the TRUH statistic in the b-th bootstrap round. Observe that compared to (7), we have introduced a tuning parameter τ f c in (8) above. We define it as the fold change (fc) hyper-parameter and will consider values of τ f c ≥ 1. Biologically relevant remodeling corresponds to significant fold change increase or decrease in the magnitude of cellular expressions between the infected and the uninfected cells. As we test the global null hypothesis of no change in any of the concerned genes, alternative hypothesis of remodeling with meager fold changes, if accepted, will only lead to biologically uninteresting discoveries. For discovering virologically interesting alternatives, it is natural to set τ f c slightly larger than 1. (Note that τ f c = 1 corresponds to the bootstrapped version of the TRUH statistic in (7).) In the simulation experiments presented later in Section 3 we set τ f c = 1 whereas in Section 4 τ f c is fixed at 1.1 as we study a real-world virology dataset.
≤ m a then STEP 2: Draw a simple random sample V the baseline residual sample in class a; else Stop: Go to STEP 1; : a = 1, . . . ,K}; the maximum number of iterations over which the K-means algorithm runs before stopping are both fixed and thus has O(md) complexity. Therefore the overall computational complexity of Algorithm 1 is O(nmd).

Numerical Experiments
In this section we evaluate the numerical performance of the TRUH procedure across a wide range of simulation experiments. We consider the following six competing testing procedures that use different methodologies to conduct a non-parametric two-sample test: (i) Energy test (Energy) of Aslan and Zech (2005)  To assess the performance of the competing testing procedures, we simulate U m and V n from F 0 and G, the cdf of the baseline population and the case population respectively, and for each testing procedure, we measure the proportion of rejections across 100 repetitions of the composite null hypothesis test described in (3) at 5% level of significance. For TRUH, we fix the fold change constant τ f c = 1 and take B = 200 bootstrap repetitions. The R code that reproduces our simulation results can be downloaded from the following link: https://github.com/trambakbanerjee/TRUH paper and the associated R package is available at https://github.com/trambakbanerjee/TRUH.

Experiment 1
In the setup of experiment 1, we consider testing H 0 : G ∈ F(F 0 ) versus H A : G / ∈ F(F 0 ), when F 0 is the cdf of a d dimensional Gaussian mixture distribution with three components: where µ 1 = 0 d , µ 2 = −31 d , µ 3 = −µ 2 , and Σ K , for K = 1, 2, 3, are d dimensional positive definite matrices with eigenvalues randomly generated from the interval [1, 10]. To simulate V n from G, we consider two scenarios as follows:  (3), direct application of traditional two-sample tests such as those considered here, is no longer conservative while TRUH is both conservative and powerful against departures from H 0 : G ∈ F(F 0 ). In Table 3 we present the results of the simulation exercise that correspond to the three scenarios described in Figure 3. The two dimensional uninfected marker expressions (X 1 , X 2 ) are randomly sampled For the panel on the left of Figure 3, infected marker expressions arise from F 0 but with sample size n = 500, while for the center panel the infected marker expressions represent a random sample of size n from G = 0.5N 2 (µ 4 , I 2 ) + 0.5 N 2 (µ 5 , I 2 ), where µ 4 = 0.25µ 2 + 0.5µ 3 and µ 5 = (3/4)µ 2 + (9/8)µ 3 .
Clearly in this case G / ∈ F(F 0 ). For the right most panel, infected marker expressions are again a random sample of size n from F(F 0 ) with mixing weights given by the vector (w 1 , w 2 , w 3 ) = (0.8, 0.1, 0.1).
Under this setting, the three g-tests , Chen and Friedman, 2017, Friedman and Rafsky, 1979, the cross-match test of Rosenbaum (2005), and the energy test of Aslan and Zech (2005) infer G / ∈ F(F 0 ) in all of the 100 repetitions of the experiment thus suggesting their inability to tackle sub-population level heterogeneity.

Experiment 2
For experiment 2, we consider a more complex setup wherein F 0 is the cdf of a d dimensional mixture distribution which is not necessarily Gaussian. Here, where Gam d and Exp d are d dimensional Gamma and Exponential distributions. For generating correlated Gamma and Exponential variables, we use the Gaussian copula approach based function from the R-package lcmix (Dvorkin, 2012, Xue-Kun Song, 2000. We consider tapering matrices with positive and negative autocorrelations: (Σ 1 ) ij = 0.7 |i−j| and (Σ 2 ) ij = −0.9 |i−j| for 1 ≤ i, j ≤ d. For simulating V n from G, we consider the following two scenarios: • Scenario I: Here, G = Exp d (rate = 1 d , Σ 2 ). In this case, G arises from only one of the components of F 0 , that is, G ∈ F(F 0 ).
In this setting, G / ∈ F(F 0 ) and the composite null H 0 is not true. When the ratio n/m is small, this scenario presents a difficult setting for detecting departures from H 0 as majority of the case samples from V n will arise from Exp d (rate = 1 d , Σ 2 ) and the tests will rely on only a small fraction of samples from Gam d (shape = 101 d , rate = 0.51 d , Σ 1 ) to reject the null hypothesis. Substantial improvement in the rejection rates is evident when m = 2000. However, for both these cases, m = 500 and m = 2000, the Energy test followed by TRUH exhibit the largest rejection rates. Although Energy test rejects H 0 in almost all of the testing instances in scenario II, its performance in scenario I (Table 4) reveals that it can be severely non-conservative when G ∈ F(F 0 ).

Experiment 3
For experiment 3, we introduce zero inflation in both the baseline and case samples to mimic the scenario that is often encountered in virology studies wherein some of the markers exhibit only a small probability of expressing themselves. We let p = (p 1 , . . . , p d ) denote the d dimensional vector of point masses at 0 across dimensions, and consider In the above representation, p regulates the differential zero inflation across the d dimensions. For the purposes of this experiment, we chose the first 0.8d coordinates of p independently from Unif(0.5, 0.6), and the remaining 0.2d coordinates are set to 0. Thus, the zero inflation is encountered only in the first 0.8d coordinates of F 0 . To simulate the baseline sample U m from F 0 , we use the R-package lcmix with Σ 1 , Σ 2 as described in experiment 2 (Section 3.2). For simulating V n from G, we consider the following two scenarios: Here, G arises from only one of the components of F 0 , that is, G ∈ F(F 0 ).
• Scenario II: Here, we let G = 0.5G 1 + 0.5G 2 , where G 1 qδ {0} + (1 d − q) Gam d (shape = 51 d , rate = 0.51 d , Σ 1 ) and G 2 = qδ {0} + (1 d − q) Exp d (rate = 1 d , Σ 2 ), and we set the first 0.8d coordinates of q to 0.3 and the remaining 0.2d coordinates to 0. Note that in this setting, apart from the difference in the rate parameter of the Gamma distribution, we also have differential zero inflation across G and F 0 , as q = p. Thus, G / ∈ F(F 0 ) and the composite null H 0 is not true.
Moreover, when n is small, this scenario presents a challenging setting for detecting departures from H 0 as the tests will have to rely on both the differences in the rate parameter and differential zero expression between U m , V n to reject the null hypothesis.

Remodeling Analysis of HIV-Infected T Cells
In this section, we analyze the data collected in Cavrois et al. (2017). It contains protein expressions of uninfected and HIV infected CD4 (which is a protein found on the surface of immune cells) positive tonsillar T cells. We show that existing two-sample tests might lead to incorrect conclusions, which can be corrected by using our proposed TRUH hypothesis testing framework. As discussed in section 1.1, the goal in Cavrois et al. (2017) was to conduct a mass cytometric assessment of subsets of CD4+ T cells that support HIV entry and viral infection in humans using two variants of the HIV virus: Nef rich HIV and Nef deficient HIV. It is known in the immunology literature, that Nef-rich cells are more prone to viral remodeling (Basmaciogullari and Pizzato, 2014). The data set we analyze here contains uninfected and infected data from two different sets of experiments. Both the experiments have four replications based on tonsillar T cells from 4 healthy donors. In Experiment I, the infection was done by Nef-rich HIV where as in Experiment II the infection was done by Nef-deficient HIV. We expect remodeling, if any, in the infected cells to be higher in experiment I than in experiment II compared to their respective baseline uninfected populations.
The cells in the data were phenotyped in a 38 parameter CyToF   expression values with cdf F 0 . We assume that the heterogeneity in the uninfected population is reflected by K heterogeneous cellular sub-groups each having unimodal probability distribution functions with cdfs F 1 , F 2 , . . . , F K and mixing proportions w 1 , w 2 , . . . , w K , such that F 0 is of the form represented in equation (1). We observe the virus infected sample V n,r = {V 1,r , . . . , V n,r } consisting of n i.i.d. Among the 35 markers considered here it is known that the expressions of the four markers CD4, CCR5, CD28, and CD62L are mainly changed due to HIV infection (Garcia and Miller, 1991, Matheson et al., 2015, Michel et al., 2005, Ross et al., 1999, Swigut et al., 2001, Vassena et al., 2015. • CASE 3 corresponds to Experiment II A where we test viral remodeling on Nef-deficient Figure 7: This is a t-SNE plot (Maaten and Hinton, 2008) of the data where the d dimensional uninfected and infected cellular expression levels are projected to a two dimensional space for each of the four donors across the four cases.

d-dimensional arcsin transformed observations from G and the goal is to test H
infected cells based on all 35 markers including the four which are known to be remodeled.
• CASE 4 corresponds to Experiment II B where we test viral remodeling on Nef-deficient infected cells based on 31 markers which are known to be mainly invariant under HIV infection.
In all the four cases, we have four replications corresponding to four donors. It has been established through validation experiments in Cavrois et al. (2017) that there is no remodeling but only preferential infection in cases 2 and 4 whereas cases 1 and 3 exhibits remodeling with the intensity of remodeling being much higher in the former than the later. Biologically, it corresponds to the fact that there is Nefindependent remodeling but the intensity of remodeling is higher in presence of Nef. Also, remodeling in cellular expressions is confined to the four markers CD4, CCR5, CD28, and CD62L in the set of markers considered in the study. Figure 7 presents a t-SNE plot (Maaten and Hinton, 2008) of the data where the d dimensional uninfected and infected cellular expression levels are projected to a two dimensional space for each of the four donors across the four cases. While these plots exhibit the underlying heterogeneity in the uninfected sample and the sample size imbalance, instances of remodeling are also visible in cases 1 and 3 wherein a relatively large fraction of the infected cells in red occupy a distinct position in the two dimensional space with no overlap with their uninfected counterparts.
For conducting statistical hypothesis test for the above four cases, along with our proposed TRUH procedure, we also use the five other competing tests statistics described in section 3 which are the Energy test (Aslan and Zech, 2005), CrossMatch (Rosenbaum, 2005), E Count (Friedman and Rafsky, 1979), GE Count (Chen and Friedman, 2017), WE Count  and MTE Count (Zhang and Chen, 2017). Figure 8 presents the values of the TRUH statistic and the 2.5-th, 50-th, 97.5-th percentiles of the associated null distribution. From the plots, it is evident that at 5% level our proposed procedure correctly captures the biological phenomena of remodeling or no remodeling across the four cases. The other five tests fail to correctly detect the phenomena correctly in some of the four cases due to heterogeneity in the data. Next, we describe the results in further details. In Tables 8 and 9 we report the p-values of the six competing tests statistics for testing remodeling under HIV infection in Nef-rich environment. In Table 8 all six tests reject the null hypothesis of no remodeling, thus verifying that CD4+ T cells exhibit remodeling under the influence of Nef rich HIV infection. In Table 9, however, we present the p-values of the six tests when the four cell surface markers, CD4, CCR5, CD28, and CD62L, known to be down regulated by Nef, were removed from our analysis (d = 31). Other than donor 1, TRUH indicates no remodeling in this scenario for the remaining three donors which is expected given the mechanism of remodeling that Nef pursues by down-regulating CD4, CCR5, CD28, and CD62L (Swigut et al., 2001). The absence of these four cell markers from the uninfected and infected samples reduces the phenotypic gap between these samples as measured through their surface markers.
The top row in Figure 8 shows that while the null distribution shifts down from CASE 1 (left plot) to CASE 2 (right plot) across all four donors, the drop in the magnitude of the TRUH statistic is far more significant when the four surface markers are excluded. The remaining five test statistics appear to be insensitive to these subtle changes in the uninfected and infected samples across the two scenarios and, continue to detect remodeling in Case 2 which is actually no remodeling but preferential infection. This demonstrates their inability to handle heterogeneity in the data that TRUH tackles via the composite null testing framework of equations (1)-(3).
In Tables 10 and 11, we present the p-values of the six test statistics for testing the null hypothesis H 0 of no remodeling when the HIV infected sample lacks the critical Nef gene (see Construction and validation of reporter viruses in Supplemental Experimental Procedures of Cavrois et al. (2017) for details around the generation of Nef-deficient HIV infected cells). We see that TRUH rejects the null hypothesis of no remodeling in CASE 3 (Table 10) while fails to do so in CASE 4 (Table 11), thus corroborating the biological phenomena that (a) Nef independent remodeling is prevalent in HIV infected cells and, (b) even in the absence of Nef, the down regulation of the four surface markers by other mechanisms contributes to remodeling. The bottom row in Figure 8 presents the values of the TRUH statistic and the 2.5-th, 50-th, 97.5-th percentiles of the associated null distribution. Similar observations from the top row continue to hold for cases 3 and 4 in the bottom row of Figure 8 wherein the drop in the magnitude of TRUH statistic is far more significant when the four surface markers are excluded.
Moreover, from Figure 8, we see that for every donor the TRUH statistic obeys a rank ordering across the scenarios which is of the form TRUH 1 > TRUH 3 > TRUH 2 > TRUH 4 where TRUH s is the magnitude of the TRUH statistic under cases s = 1, . . . , 4. This is not accidental for the relative strength of remodeling is known to be highest under the influence of Nef-rich HIV infection and more so when Nef downregulates the four cell surface markers, CD4, CCR5, CD28, and CD62L. As was seen in cases 1 and 2, the remaining five tests continue to side in favor of remodeling in both cases 3 and 4, thus reflecting  their relative lack of conservatism in detecting remodeling under our composite null testing framework.  The remodeling analysis of the HIV-infected T Cells using TRUH reveals that our proposed testing procedure conforms to the biologically validated phenomenon of remodeling of human tonsillar T cells under both Nef-rich (case 1) and Nef-deficient (case 3) HIV infection. However unlike traditional tests that continue to infer remodeling in cases 2 and 4, TRUH detects preferential infection and concludes that phenotypic differences between the HIV infected and uninfected T cells are primarily driven by variations in the expression levels of CD4, CCR5, CD28, and CD62L across the uninfected and infected cells. Moreover, through cases 1 and 2, TRUH corroborates the findings in Chaudhuri et al. (2007), Michel et al. (2005), Swigut et al. (2001), Vassena et al. (2015) that HIV remodeling of the T cells is driven by Nef dependent down-regulation of CD4, CCR5, CD28, CD62L while through cases 3 and 4 TRUH reveals Nef independent remodeling of T cells as evidenced in Cavrois et al. (2017).

Optimality Properties of the TRUH Statistic
In this section we derive the L 2 -limit of the proposed test statistic T m,n in the usual limiting regime where the sample sizes m, n → ∞, such that n/m → ρ > 0. This can be used to choose a cut-off and construct a test based on T m,n , and show asymptotic consistency for biologically relevant location alternatives.
Recall, that the uninfected and infected samples are denoted as which are i.i.d. samples from two unknown densities f 0 and g in R d , respectively. To derive the limit of T m,n we need certain integrability/moment assumptions on f 0 and g.
Assumption 1. The densities f 0 and g have a common support S ⊆ R d and satisfy either one of the following two assumptions depending on the dimension: 1. For d ≤ 2, the support S is compact (with a non-empty interior) and f 0 and g are bounded away from zero on S.
To describe the limit of T m,n we need a few definitions: For λ > 0, denote by P λ the homogeneous Poisson process of intensity λ ≥ 0 in R d , and P x λ = P λ ∪ {x}, for x ∈ R d . Now, define the following two quantities: that is, the distance from the origin 0 0 0 in R d to its nearest neighbor in the Poisson process P 1 , and the distance of this point to its neighbor in P 1 , respectively.
Theorem 1. Let T m,n be as in (7). Then, for f 0 and g as in Assumption 1 above, as m, n → ∞ such that n/m → ρ, with ∆ d = (ζ 2 − ζ 1 ), where ζ 1 = Eζ 1 (0 0 0, P 1 ), the expected distance from the origin in 0 ∈ R d to its nearest neighbor in P 1 , and ζ 2 = Eζ 2 (0 0 0, P 1 ), the excepted distance between the nearest neighbor of the origin in P 1 to its nearest neighbor in P 1 .
The above theorem gives the L 2 -limit of the test statistic for general distributions f 0 and g. The proof of the theorem, which is given in the supplementary materials (Section B), uses the machinery of geometric stabilization, introduced by Penrose and Yukich (2003), which obtains the asymptotics of nearest neighbor based functionals in terms of functionals defined on a homogeneous Poisson process.
Before we discuss how the result in Theorem 1 can be used to construct a test based on T m,n for the hypothesis (2), we discuss some properties and the consequences of the limit in (11): • Note that the finiteness of the limit in (11) is ensured by Assumption 1. For d ≥ 3, the moment conditions in Assumption 1 are required to establish the L 2 convergence in (11). This assumption can be relaxed to S |y| r f 0 (y)dy < ∞ and S |y| r g(y)dy < ∞, for some r > d/(d − 1), if we are only interested in L 1 convergence (by combining the proof of Theorem 1 with that of (Penrose and Yukich, 2003, Proposition 3.2)). However, this still does not apply for d = 1, where it is necessary to assume the compactness of the support, in order to ensure that the limit in (11) is finite. This is a well-known constraint which arises in a large family of random geometric graphs, while dealing with the asymptotics of edge-lengths (see, for example, (Penrose and Yukich, 2003, Theorem 1.1) and the references therein). Even though the compactness assumption technically rules out some natural distributions, from a practical standpoint, there is no real concern because one can approximate the univariate density by truncating it to a large interval, on which the above result applies.
• Note that ζ 1 and ζ 2 are both constants, which depend only on the dimension d. In fact, ζ 1 has a closed form expression which can be easily derived. To this end, denote by V d and S d the volume and the surface area of the unit ball in R d , respectively. It is easy to verify that S d = dV d .
Moreover, for r > 0 and x ∈ R d , denote by B(x, r) the ball of radius r centered at x ∈ R d . Then, using the observation that a point b is the nearest neighbor of the origin, if are there no points of the Poisson process P 1 in the ball B(0, ||b||), it follows that which, by the change of variable x = V d t d equals where Γ(·) denotes the Gamma function.
Theorem 1 shows that for K fixed densities f 1 , . . . , f K , and where the last step uses the fact that λ a ∈ [0, 1], for 1 ≤ a ≤ K, and K a=1 λ K = 1. Note that the RHS above is unknown, because the densities f 1 , . . . , f K , the weights w 1 , . . . , w K , as well as the number K of mixture components, are all unknown. However, if we can consistently estimate the RHS of (13), then the test which rejects H 0 in (3) when T m,n is greater than the estimated value of (13), would have zero asymptotic Type I error and would be powerful whenever g has some separation from the set F(f 0 ) (recall definition in (2)).
The approach described above is, in general, infeasible because non-parametric estimation of mixture parameters in multivariate problems, especially when the number K is unknown, can often be difficult.
In the following, we show how in location families, one can obtain a slightly weaker upper bound on ϕ(f 0 , g, ρ), which is free of the unknown parameters, that can be used to construct a valid and powerful test for the remodeling hypothesis (3). To this end, consider {p(y|θ) = p(y − θ) : θ ∈ Θ} a family of densities indexed by the parameter space Θ ⊆ R d , where p : R d → R ≥0 such that R d p(y)dy = 1. Throughout we assume that the densities in the family satisfy Assumption 1. Suppose the baseline samples U 1 , U 2 , . . . , U m are i.i.d. from the density f 0 (·) = K a=1 w a p(·|θ a ), where θ 1 , . . . , θ K ∈ Θ are fixed (but unknown), and there exists a known constant L > 0 such that w a ≥ L, for all 1 ≤ a ≤ K. If the infected samples V 1 , V 2 , . . . , V n are i.i.d. from a density g in R d , then the hypothesis of remodeling (2), in this parametric setting, becomes, where θ = (θ 1 , . . . , θ K ) and F(θ) is defined as follows: is the collection of K-mixtures of p(·|θ 1 ), p(·|θ 2 ), . . . , p(·|θ K ). Note that under the null H 0 , g(·) = K a=1 λ a p(·|θ a ), for some λ 1 , λ 2 , . . . , λ K ∈ [0, 1], such that K a=1 λ a = 1. Then using where the last step follows by the change of variable z = y − θ a . Note that the constant γ depends on L (the lower bound on the mixing weights of the baseline population), the dimension d, and the base function p defining the location family (which is assumed to be known); but not on the unknown means (θ 1 , θ 2 , . . . , θ K ), the unknown weights (w 1 , w 2 , . . . , w K ), or the number of components, and hence can be directly calculated.
The corollary below shows how the bound in (15) can be used to construct a test based on T m,n which is powerful for radially symmetric mixtures, such as Gaussian mixtures and t-mixtures, among others.
The proof of the corollary is given in the supplementary materials (Section C). Note that the condition on g(y) in (17) quantifies a natural notion of separation between g and the set F(θ), by assuming that at least one of the mixture means of g is ε-far (in L 2 -distance) from all the unknown null means of the baseline density. Explicit bounds on the separation ε(γ) can be obtained from the proof of Corollary 1, based on the tail decay of the base density p (details given in supplementary materials, Section C).

Discussion
We propose a novel nearest neighbor based two-sample test for detecting changes between the baseline and the case samples, in the presence of heterogeneity, as is often the case in single-cell virology. Our testing procedure is specially designed for mass cytometry based techniques (Bendall et al., 2011, Giesen et al., 2014 which produces moderate dimensional (d ∼ 50) cellular characteristics. In the future, it will be interesting to extend our methodology for dealing with single-cell RNA-seq based techniques , Hwang et al., 2018, Jaitin et al., 2014, Schiffman et al., 2017, which can produce highly multivariate phenotypes (d ∼ 10 4 ). A possible approach can be using the random projection wrapper on our testing procedure. Also, it will be interesting to develop efficient testing procedures where the underlying population contains heterogeneous sub-populations with highly varying sizes including some very rare sub-populations.

B Proof of Theorem 1
The proof of Theorem 1 is an immediate consequence of the following two lemmas. The first lemma computes the limit of n 1 dD m,n . Lemma 1. Let D 1 , D 2 , . . . , D n be as defined in (5). Then, under Assumption 1, as m, n → ∞, where ζ 1 is as defined in the statement of Theorem 1.
The next lemma computes the limit of n 1 dC m,n , which combined with Lemma 1 completes the proof of Theorem 1.
The proofs of Lemma 1 and Lemma 2 are given below in Section B.2 and Section B.3, respectively.
We begin with some preliminaries about Poisson processes and stabilization of geometric functionals, introduced by Penrose and Yukich (2003), in Section B.1 below.

B.1 Preliminaries
Given z ∈ R d , denote by ϕ(z, Z) a measurable R + valued function defined for all locally finite set Z ⊂ R d and z ∈ Z. If z / ∈ Z, then ϕ(z, Z) := ϕ(z, Z ∪ {z}). The function ϕ is said to be translation invariant if ϕ(y + z, y + Z) = ϕ(z, Z). Penrose and Yukich (2003) defined stabilizing functions as follows: Definition 1. (Penrose and Yukich (2003) We will be interested in functionals that stabilize almost surely on P λ , the homogeneous Poisson process with rate λ in R d . Note that with probability 1, ϕ(P λ , M ) is non-increasing in M and ϕ(P λ , M ) is non-decreasing in M , therefore, they both converge. The definition of stabilization in (21) means they converge to the same limit almost surely. Note that any functional ϕ(z, Z) which depends only on the points of Z within a fixed distance of z is stabilizing on P λ . In our proofs, we will consider the following two functionals: • For y ∈ R d , and Z ⊂ R d finite, define which is the distance from y to its nearest neighbor in Z.
using the inequality ab ≤ a 2 +b 2 2 and the fact m j=1 1{U j = N (V 1 , U m )} = 1. Now, for n large enough, where the functional φ(A, B) := a∈A b∈B ||a − b|| 2 1{b = N (a, B)}, where A, B ⊂ R d are finite and disjoint. Note that for any partition {S 0 , S 1 , . . . , } of R d , that is, the functional φ is sub-additive. (Note that the sum above is, in fact, finite because the sets A and B are finite.) Then by a modification of (Yukich, 2006, Lemma 3.3), one can obtain the growth bound d . Now, choosing S 0 to be the ball of radius 2 centered at the origin, and S K to be the annulus centered at the origin with inner radius 2 K and outer radius 2 K+1 , for K ≥ 1, it follows from (28) that . Now, taking expectations above and the Jensen's inequality gives, for n large enough, both of which are finite by the integrality assumptions on f 0 and g (using arguments in (Yukich, 2006, Page 85)). The result now follows by combining the bound above with (26) and (27).
Then, recalling n/m → ρ gives, which establishes the limit in (19) in expectation.
To complete the proof of the lemma we need to show that the variance of the LHS in (19) goes to zero. To this end, note that where, as before, V is a random variable distributed according to the density g, and P f0(V ) is a Cox process with intensity measure f 0 (V ). This combined with (31) and (30), shows that Var 1 This completes the proof of Lemma 1.
which establishes the limit in (20) in expectation.
Finally, similar to the proof of Lemma 1, it can be shown that the variance of the LHS in (20) goes to zero, completing the proof.
However, by an application of the FKG inequality for Poisson processes (Janson, 1984, Last andPenrose, 2017), it can be shown that ζ 2 ≥ ζ 1 . This is described in the following remark.
Note that the separation ε(γ) depends on η(M, d), the rate of decay of the tail of the base density p.
For instance, when p is the standard multivariate normal distribution N (0, I d ), then it suffices to choose η(M, d) = K(d) √ log M , where K(d) is a constant depending on d.