Mixture modeling on related samples by $\psi$-stick breaking and kernel perturbation

There has been great interest recently in applying nonparametric kernel mixtures in a hierarchical manner to model multiple related data samples jointly. In such settings several data features are commonly present: (i) the related samples often share some, if not all, of the mixture components but with differing weights, (ii) only some, not all, of the mixture components vary across the samples, and (iii) often the shared mixture components across samples are not aligned perfectly in terms of their location and spread, but rather display small misalignments either due to systematic cross-sample difference or more often due to uncontrolled, extraneous causes. Properly incorporating these features in mixture modeling will enhance the efficiency of inference, whereas ignoring them not only reduces efficiency but can jeopardize the validity of the inference due to issues such as confounding. We introduce two techniques for incorporating these features in modeling related data samples using kernel mixtures. The first technique, called $\psi$-stick breaking, is a joint generative process for the mixing weights through the breaking of both a stick shared by all the samples for the components that do not vary in size across samples and an idiosyncratic stick for each sample for those components that do vary in size. The second technique is to imbue random perturbation into the kernels, thereby accounting for cross-sample misalignment. These techniques can be used either separately or together in both parametric and nonparametric kernel mixtures. We derive efficient Bayesian inference recipes based on MCMC sampling for models featuring these techniques, and illustrate their work through both simulated data and a real flow cytometry data set in prediction/estimation, cross-sample calibration, and testing multi-sample differences.


Introduction
Kernel mixtures are a powerful tool for modeling a variety of data sets, especially in the presence of a natural clustering structure (Escobar and West, 1995;MacEachern and Müller, 1998). A good portion of the rapidly expanding literature on Bayesian nonparametrics is aimed at building effective mixture models. A recent focus of the literature is on how to jointly model in a hierarchical manner data samples that are similar or otherwise related, the main objective being effective borrowing of strength across samples, thereby substantially enhancing inference on the underlying data generative mechanisms as well as prediction. This is particularly important for complex data sets, for which each individual sample may only contain very limited information regarding the underlying probability distribution. Among many notable efforts in this direction, Lopes et al. (2003) proposed a hierarchical model for multiple finite mixtures. Müller et al. (2004) proposed a nonparametric extension of Lopes et al. (2003)'s model by replacing finite mixtures with Dirichlet process (DP) mixtures. In a different vein, Cron et al. (2013) proposed to use the hierarchical DP, or HDP, (Teh et al., 2006) as the mixing distribution to characterize variation across multiple mixture distributions. Rodríguez et al. (2008) proposed the nested DP (NDP) mixture, which is an infinite mixture of DP mixtures that induces an additional level of clustering among multiple mixture distributions themselves (to be distinguished from the clustering within each mixture distribution).
While applicable to a variety of mixture modeling contexts, our work is motivated during our attempt to apply existing hierarchical mixture models to the analysis of data collected from flow cytometry experiments. Flow cytometry is a laser-based technology that measures biomarkers on a large number of cells, so each cell is an observation from a distribution in R p , where p is the number of biomarkers measured. The cell population typically comes from a blood sample in immunological studies, and it consists of cells of various subtypes-e.g., T cells, B cells, etc.-with each subtype forming a "cluster" in the sample space. Because each cell subtype has a specific function in the immune system, inference on the abundance of the various subtypes across blood samples of a patient under different stimulating conditions, for instance, is of interest. Mixture models are natural tools for characterizing such data as the data is indeed a mixture of various cell types (Chan et al., 2008), and because a typical flow cytometry study will involve multiple samples collected under different conditions, the need for jointly modeling to achieve effective borrowing of strength also naturally arises (Cron et al., 2013).
During the analysis of flow cytometry experiments using mixtures, we encountered a number of important challenges that we believe are present in numerous (if not most of) other applications involving mixture modeling of related samples (not only with locationscale kernels but beyond). Below we summarize the three main data features/challenges that motivate the current work: I. Samples often share clusters but with differing weights. Related samples tend to share some (even most) of their clusters, and these common clusters vary across related samples in their weights. In flow cytometry, for instance, data samples often share a vast majority of the cell subtypes, and the most common type of variation across samples is the differences in the relative sizes of the subtypes.
II. Only some, not all, clusters vary. Often, only a fraction, not all, of the clusters vary across samples. In flow cytometry, not all cell subtypes are affected by the experimental conditions of interest. Very often only one or two cell types are affected and thus vary across the samples while the rest do not.
III. Misalignment across samples in shared clusters. Even the same cluster shared among samples is often not perfectly aligned across samples, either due to actual systematic difference across the samples, or very often due to the presence of extraneous, uncontrolled additional sources of variation, i.e., some "random" effect. This is easily seen in mixtures of location-scale families, where the location and spread of some shared clusters differ to various extent across samples. Such misalignment is ubiquitous in flow cytometry data, with numerous potential causes. For example even tiny differences in the chemical concentrations applied in the experimental protocol across experiments can cause noticeable "perturbations" in the cell subtypes.
As far as we know, none of the existing hierarchical approaches satisfactorily address all of these issues in a single coherent framework. Table 1 provides a summary of these data features and the extent to which some of the state-of-the-art methods (along with the method we propose herein) address each of them.

Shared clusters
Only a subset Misalignment with varying weights of clusters differ in kernels Lopes et al. (2003); Müller et al. (2004) Not allowed Allowed Not allowed Teh et al. (2006); Cron et al. (2013) Allowed Not allowed Not allowed Rodríguez et al. (2008) Not allowed Not allowed Not allowed This work Allowed Allowed Allowed Table 1: Comparison of hierarchical mixture models in terms of how they cope with the three common data features/challenges in modeling multiple related data samples.
Specifically, the existing approaches exploit some aspects of these features but do not fully take them into account. By introducing a cluster-specific hierarchical relationship among the samples, Lopes et al. (2003) and Müller et al. (2004) allow some clusters to be shared among the samples. However, their models require that the kernel parameters and the mixture weight for each cluster be either both shared across samples or both different, without the option to decouple these two different types of variations. In particular, no clusters are allowed to have only one type of variation-e.g., mixing weights-under these models. In the context of flow cytometry, for instance, this would mean that cell subtypes cannot change just in abundance across the samples but not in their location and spread, clearly an unrealistic assumption. On the other hand, by using the hierarchical DP (Teh et al., 2006) as the mixing distribution, Cron et al. (2013) does allow variations to exist in weights alone, but enforces the constraint that all clusters must all vary across samples, excluding the common situation in applications such as flow cytometry that only some clusters (e.g., subtypes) vary while others remain unchanged across conditions. Finally, under the nested DP mixture (Rodríguez et al., 2008) In Section 3 we compare our method to current methods through simulation studies that cover prediction/estimation, cross-sample calibration, and testing multi-sample differences, and finally use it to analyze two flow cytometry data sets.

Background: Dirichlet process mixtures and stick breaking
While our techniques can be embedded into mixture models with various weight generating mechanisms and kernel families, we shall introduce and illustrate them in the context of DP mixtures of Gaussians, which is the most widely adopted nonparametric mixture model.
Suppose n observations y = (y 1 , y 2 , . . . , y n ) are from a mixture model: where f denotes the probability density function of F , g(·|λ) is a kernel distribution parametrized by λ, π k the associated (mixture) weight, and K the countable (possibly infinite) index set of the mixture components (or clusters). Location-scale families are commonly adopted as the kernel distribution, in which case λ k specifies the location and spread of the kth cluster. By definition the weights satisfy π k ≥ 0 and k π k = 1. An alternative and computationally attractive formulation utilizes a latent cluster membership label Z i ∈ K for each observation, such that y i | Z i = k ∼ g(·|λ k ) and Pr(Z i = k) = π k for i = 1, 2, . . . , N and k ∈ K.
Bayesian inference under mixture models can proceed after specifying prior distributions on the weights and the kernel parameters {(π k , λ k ) : k ∈ K} (Marin et al., 2005). A flexible and convenient choice on the prior for the mixing weights is a generative procedure called the stick breaking process (SBP) (Sethuraman, 1994;Ishwaran and James, 2001).
The general scheme of SBP starts with the drawing of a sequence of independent random variables v 1 , v 2 , . . . supported on (0, 1). Then the weight for the kth cluster is given as A popular two-parameter specification is the Poisson-Dirichlet process (Kingman, 1975;Pitman and Yor, 1997), corresponding to v i ∼ Beta(1 − γ, α + γ) for some parameters α and γ. In particular, when γ = 0, this boils down to the weight generative mechanism from a Dirichlet process (Ferguson, 1973;Sethuraman, 1994), which we shall refer to as the SBP(α) process.
By adopting the SBP(α) prior on the weights, along with a prior H on the kernel parameters, we obtain a Dirichlet process mixture (DPM) model: The most commonly adopted kernel distributions are location-scale families such as the (multivariate) Gaussian family, i.e., g(·|λ k ) = N (·|µ k , Σ k ). In this case, H is often chosen to be the corresponding conjugate prior such as a normal-inverse-Wishart (NIW) prior on (µ k , Σ k ).

Two techniques for hierarchically modeling related samples
Now assume J samples of observations y j = (y 1,j , . . . , y n j ,j ) for j = 1, . . . , J have been collected, and the observations in each sample are modeled by a mixture: where f j is the probability density function of F j , and λ j,k represent the kernel parameter for the kth cluster in the jth sample. To characterize potential relationship across the samples, let us assume that the kth component under each sample represent the same cluster (e.g., cell subtype). Note that this does not exclude the possibility of having novel clusters that appear in only one or some of the samples, in which case the weights π j,k = 0 if cluster k is absent in the jth sample. Again we let K be the collection of all cluster indices over all the samples. Let Z i,j be a latent variable indicating that the data point y i,j belongs to the kth cluster with k ∈ K. Then the model can be equivalently written as We next introduce techniques for prior choices on the weights and on the kernel parameters by extending the stick breaking prior and the kernel respectively, which will address the three data features and challenges described in the Introduction.
ψ-stick breaking for weights We consider a generative stick breaking procedure called "ψ-stick breaking" (for reasons to be explained below), which breaks J sticks of unit length-one for each sample-in a dependent manner to generate the mixing weights {π j,k : k = 1, 2, . . .} for j = 1, 2, . . . , J. We start by observing that each cluster falls into one of two categories K 0 and K 1 , that is K = K 0 ∪ K 1 with K 0 ∩ K 1 = ∅: those in K 0 have weights that do not vary across the J samples (e.g., cell types whose abundance is constant across experimental conditions), i.e., π j,k = π j ,k for j, j = 1, . . . , J for k ∈ K 0 , whereas those in K 1 have varying weights across samples.
The generative process proceeds in two steps and is illustrated in Figure 1. In the first step, we break the J sticks at exactly the same spot into two pieces of length ρ and 1 − ρ  respectively, where ρ ∈ (0, 1) is drawn as a Beta random variable. Then in the second step, we use the J pieces of length ρ to generate the weights for the components in K 0 , and the J pieces of length 1 − ρ for the subtypes in K 1 . Hence the parameter ρ is interpreted as the overall proportion of the clusters with constant weights across samples.
Specifically, one can imagine that we tie the J sticks of length ρ together and break them using a single SBP as if they were a single stick-always at the same locations. For this reason, we shall refer to the common stick formed by tying the J sticks of length ρ as the "shared" stick, or the s-stick. Let {w 0,k : k ∈ K 0 } with k∈K 0 w 0,k = 1 be the randomly generated relative sizes of the components in K 0 in terms of the proportions of the s-stick. So the absolute size of each cluster that does not change across samples is given by π j,k = ρw 0,k for all j = 1, 2, . . . , J and k ∈ K 0 .
On the other hand, we break the J sticks of length 1 − ρ independently using separate independent SBPs, each generating the weights for one of the J samples, corresponding to the sizes of clusters that vary across samples. For this reason, we shall refer to the J sticks of length 1 − ρ as the "idiosyncratic" sticks, or the i-sticks. We let {w j,k : k ∈ K 1 } for j = 1, 2, . . . , J with k∈K 1 w j,k = 1 be the randomly generated lengths of the components as proportions of the corresponding i-stick. So for the kth cluster, its weight in the jth sample is given by π j,k = (1 − ρ)w j,k .
Using SBP(α) processes for breaking each of the s-and i-sticks, we arrive at a joint gen-erative model for the weights in all of the J samples, which we call "shared/idiosyncratic" (si or ψ) stick breaking. Specifically, with a Beta prior on the length of the shared stick, we arrive at the following hierarchical model for weights See Figure 1 for a visualization of the hierarchical prior on the mixture weights.
The hyperparameter α specifies the size of the clusters as well as the number of clusters (in K 0 and K 1 respectively), with a smaller α corresponding to a small number of large clusters and a larger α corresponding to a large number of small clusters. We infer on α in a hierarchical Bayesian paradigm by placing Gamma hyperprior on it: α ∼ Gamma(τ α,1 , τ α,2 ).

Local kernel perturbation
We utilize a hierarchical setup to incorporate local perturbation in the kernel parameters, thereby adjusting for the misalignment and allowing more effective borrowing of information across the samples on each cluster. Specifically, we model the kernel parameters {λ j,k } as follows where λ 0,k represent the cross-sample "centroid" kernel parameters for the kth cluster, with a hyperprior H 0 specified by hyperparameter φ 0 . Given λ 0,k , the sample-specific kernel parameters for the kth cluster λ j,k is drawn from H with additional hyperparameter , which specifies the dispersion of cluster k among the samples around the "centroid".
The above specification enforces that each cluster k will have misalignment. More generally, in some problems misalignment may exist in only a subset of the clusters. To allow for such cases, again appeal to a "spike-and-slab" setup by introducing an additional Bernoulli latent indicator S k for each cluster, such that S k = 1 if there is misalignment in cluster k whereas S k = 0 if otherwise. That is, where δ · represents a point mass.
Putting the pieces together in the context of Gaussian kernels, we arrive at the following spike-and-slab version of the locally perturbed kernel model: This model is illustrated in Figure 2. The hyperparameter specifies the total amount of local variation between the means of each group µ j,k and the grand mean µ 0,k , and ϕ specifies the proportion of clusters that have misalignment. The hyperparameters m 1 , Ψ 1 , k 0 , , and ϕ are all characterizing "global" features of the data that pertain to all of the clusters and samples. We can reliably infer them by pooling information through hierarchical Bayes. In particular, in our numerical examples we adopt the following hy-  (Ishwaran and James, 2001), which is more efficient in terms of mixing and computational speed, and hence is what we recommend. Figure 2: A locally perturbed Gaussian kernel with a spike-and-slab setup. When S k = 0, all kernels for the kth cluster are identical across samples. When S k = 1, the kernel is centered around a common mean but are not identical.

Posterior inference based on MCMC sampling
To this end, two different finite approximation strategies are commonly adopted for DPMs and other stick breaking mixtures: (i) truncating the stick breaking at some maximum number of components and (ii) using finite-dimensional symmetric Dirichlet distribution. These two approximations might look very different at first, but the main difference between the two is in the induced stochastic ordering of the weights, which is irrelevant in mixture models. In fact, as Kurihara et al. (2007) points out, one can apply a size-biased permutation to the order of the weights of a finite symmetric Dirichlet distribution and obtain a distribution which is practically identical to the truncated SBP. However, the two strategies are not computationally equivalent for mixture models. The weights under the symmetric finite-Dirichlet approximation are exchangeable, which results in substantially improved mixing over truncating the SBP. Therefore we opt for the symmetric finite Dirichlet approximation in our implementation. This approximation has been studied and used by many authors in a variety of contexts. See Neal (2000), Green and Richardson (2001) and Ishwaran and Zarepour (2002), among others. Specifically, under this approximation, the infinite sequences of mixture weights in Eq. (1) are replaced by: where K 0 and K 1 represent the numbers of mixture components that are shared and differential across the groups, respectively. In the nonparametric case, both K 0 and K 1 are infinite, while in the finite approximation we need to choose K 0 and K 1 . A simple choice is to set K 0 = K 1 = K for some large K which represents an upperbound to the a priori expected number of mixture components.
3. Latent perturbation state variables for k ∈ K:

7.
A Metropolis step to explore different modes of the posterior distribution by swapping an index from K 0 with an index from K 1 . The proposal distribution is defined as follows. An initial index k is drawn proportionally to √ n j,k for k ∈ K, where n j,k = |(i, j) : Z i,j = k|, and a second index k is drawn uniformly from K 0 if k ∈ K 1 and uniformly from K 1 if k ∈ K 0 . Since the proposal is symmetric, the swap is accepted with probability: where Z and Z new represent the vectors of the latent assignments before and after the swap. Since the mixture components are exchangeable within K 0 and K 1 , the acceptance probability depends only on the swapped indices. Similar strategies to improve the exploration of the sample space have been proposed by Porteous et al. (2012) and Papaspiliopoulos and Roberts (2008).
8. The Dirichlet pseudo-count parameter α is updated using a Metropolis-Hastings step with the following proposal: where is a is a tuning parameter calibrated in the burn-in. 9. Mean shrinkage parameter 12. The perturbation parameter is updated using a Metropolis step with the following proposal: 14. The "length" of the shared stick [ρ| . . .] ∼ Beta(a ρ + n 0 , b ρ + j n j ), where n j = k n j,k .

Numerical examples
In this section we provide three numerical examples. In the first example data are simulated under different mixture distributions, and we compare the goodness-of-fit of our method with respect to competing approaches. In the second example we illustrate through a simulated dataset how our model can be used to remove small distributional shifts across related mixture distributions. In the third example we compare the performance of our model to other competing methods in testing and identifying differences across distributions. In the fourth example we analyze two real flow cytometry datasets. In all of the examples, we shall refer to our Dirichlet process mixtures of Gaussians with ψ-stick breaking and kernel perturbation as CREMID, as it models Closely RElated MIxture Distributions.

Example 1: Estimation and predictive performance
In this first example, we investigate how CREMID helps achieve more effective borrowing of information across samples thereby enhancing predictive performance. To this end, we consider four simulation scenarios, representative of a vast variety of real applications. We use the sum of L 1 distances of the estimated univariate predictive densities from the true densities as measure of goodness of fit. (Note that we used this metric instead of the more natural log predictive score or the L 1 distance between the multivariate predictive density from the true density, because at the time of writing, the available software for the competitor HDPM provides the marginal predictive densities but not the other two metrics.) We consider the following multi-sample scenarios in R 4 . In each scenario, there are three data samples (j = 1, 2, 3) and the sample size for each is 100. Below we outline the four different scenarios. Some of the parameters are omitted here, but provided in the Appendix.
We compare our method to Müller et al. (2004)'s hierarchical Dirichlet process mixture (HDPM) method. We use the R package DPpackage (Jara et al., 2011) for fitting HDPM.
In addition, we also compare these to methods to independent finite mixture of Gaussians for each of the three samples, using Mclust (Fraley and Raftery, 2002), available in the R package mclust. In Figure 3 we show the sum of L 1 distances of the estimated univariate predictive densities from the true densities for the three methods. Our approaches outperform HDPM and mclust in the two shift scenarios. CREMID is the most accurate method in the two location shift scenarios as well as in the local weight change scenario. In the global weight change scenario, both our method and HDPM underperforms Mclust. Because the samples are different in all cluster weights, we pay a price for assuming that some cluster weights are shared.

Example 2: Correcting for cross-sample misalignment
A common problem in studies involving data collected from multiple labs or centers is the misalignment of the same clusters across samples due to external confounders, which is what motivated our hierarchical locally perturbed kernel construction. In flow cytometry, for example, misalignment across cell subpopulations can be substantial. An important preprocessing step is cross-sample calibration-that is, to estimate and correct for the misalignment across samples and thereby produce "standardized" data sets for follow up studies. (This shares the registration problem in functional data analysis.) To this end, we note that for each observation y i,j , if Z i,j = k, that is, the observation belongs to cluster k, then we can compute a corrected value by adjusting for the shift in the cluster center across the samples:ỹ where ∆ j,k = µ j,k − µ 0,k is the displacement of cluster k in sample j relative to the centroid. Because Z i,j is unobserved, we can appeal to Bayesian model averaging (BMA) by computing the posterior mean ofỹ i,j is the bth posterior draw on the displacement ∆ Let us consider a numerical example based on mixture of normals in R 4 to illustrate how one can remove cross-sample misalignment. The data are generated as follows: where i = 1, . . . , 1000, µ j,1 = (1, 10 − j, 1, 9), µ 2 = (8, 8,8,8), µ 3 = (1, 1, 1, 1) and µ j,4 = (6 + j, j, 7, 1). The three plots in the first row of Figure 4 show the data pro- In the second row of Figure 4 the three plots show the calibrated data, i.e., after removing the estimated kernel perturbations. The model is able to correctly remove the local distributional shifts across the samples. show the calibrated data, i.e., after removing the estimated kernel perturbations.

Example 3: Testing cross-sample differences in cluster weights
We consider the same multi-sample scenarios in R 4 used in Example 1. For each dataset we define a corresponding null data set by permuting the labels of the three samples. In Figure 5 we compare the ROC curves of our method and HDPM for testing the hypothesis that the three distributions are identical. Our method is substantially more powerful than HDPM in all four scenarios.
In these simulations, for our method we use E(ρϕ|y) as the test statistic. This quantity goes to zero when there are differences in the mixture weights or in the mixture kernels across samples, and it goes to one when the distributions are identical across samples. One can adopt different test statistics under our method depending on the inference objective.
For instance, if one is interested in testing just the presence of differences in weights then a suitable test statistic is E(ρ|y).
We compare our method only to HDPM since Mclust does not provide a way to test for differences across samples. In HDPM each F j is defined as a mixture of two components: The distribution H 0 represents the common part, and H j represents the idiosyncratic part. The hyperparameter controlling the "degree of similarity" across the F j 's has a beta hyperprior. We use E( |y) as the test statistic.

Application: flow cytometry
In flow cytometry experiments, biomarkers are measured on a large number of blood cells.
Different cell subtypes, i.e., groups of cells sharing similar biomarker's levels, have distinct functions in human immune system. Identifying variations in the abundance of subtypes across multiple samples is an important immunological question. Additionally, the location of a given subtype across samples can slightly change due to both experimental variability and other uncontrolled "random effects".
We analyze two datasets where each one contains three samples of 5,000 blood cells, and for each cell six biomarkers have been measured.

A control study
The blood from a given patient was split in three samples, and each sample went through a separate experimental procedure to generate the data. Since the three samples are essentially biologically identical, one expects no variations in the abundance of the different subtypes or large location shifts of the cell types. Small perturbations of the cell types are likely due to additional variations in the experimental procedures.
In Figure 6 we plot the posterior distributions of ρ and for this data set under our proposed model. The parameter ρ reflects the total mass assigned to mixture components where the mixture weights are identical across groups. In this dataset a posteriori this parameter concentrates around one, indicating that there is no evidence of a difference in the mixture weights across the three replicates. The parameter controls the expected amount of shift in the location of each kernel across samples. Its posterior does not concentrate around zero, indicating the presence of small misalignment among the replicate samples due to uncontrolled sources of variation. It is the decoupling of these two sources of variations that allows us to correctly infer the absence of variations in the mixture weights across the distributions of the three samples.

Samples under different stimulation conditions
In another data set, three blood samples from an individual underwent different stimulation treatments. One sample was left unstimulated, while the two remaining samples were stimulated with CEF and CMV pp65, respectively. The samples underwent separate experimental procedures in data generation. In Figure 7 we plot the posterior distributions of ρ and . The parameter ρ concentrates around 0.6, indicating that there are differences in some of the mixture weights across the three samples. The parameter concentrates around 0.2, either due to effects of the experiment conditions on the locations of the kernels, which is also a systematic cross-sample difference, or substantial additional variations in the experimental procedures in comparison to the control study.
To judge the goodness-of-fit, we also compare the predictive performance of our model with Mclust, evaluated by the log predictive likelihood of the a "test" sample. We randomly select 1,000 data points from the whole data set as a "test" sample, while using 5,000 observations as the "training sample". We had hoped to compare our method to other methods such as Müller et al. (2004) but at the time of writing, the existing software in R (the HDPMdensity function in DPpackage) crashes for the data sets, most probably due to the large sample sizes, and it does not output predictive scores.

Conclusion
In this work we have introduced two useful techniques in modeling related data sets using mixture models-the shared-idiosyncratic stick breaking and the locally perturbed kernel. When used together, they incorporate three common data features observed in real applications-(i) samples often share the same clusters with different weights; (ii) only some clusters vary across samples; (iii) misalignment in the clusters due to extraneous causes.
We have derived Bayesian inference recipe through MCMC sampling and carried out an extensive numerical studies to illustrate the gain in inferential efficiency in both estimation, prediction, and hypothesis testing.
Finally, we note that while the two techniques are introduced and demonstrated in the context of mixtures of location-scale families, they are generally applicable to modeling related mixtures of other forms of kernels as well, such as mixtures of generalized linear models and mixtures of factor models. The computational details will vary but the general ideas remain the same.