Bayesian Distance Weighted Discrimination

Treating $\mathbf{y}$ as a discrete random vector, it follows directly from (4) that its entries $y_{i}$ are conditionally independent and that

$$P(y_{i}\mid{\bm{\beta}},\beta_{0},{\bf x}_{i})\propto P(y_{i})e^{-V\left(y_{n}(\beta_{0}+{\bf x}_{i}^{T}{\bm{\beta}})\right)}.$$

The conditional distribution for $y_{i}\in\{-1,1\}$ thus depends on its unconditional class probability $P(y_{i}=1)=1-P(y_{i}=-1)$, and its corresponding score $u_{i}:=\beta_{0}+{\bf x}_{i}^{T}{\bm{\beta}}$:

$\displaystyle P(y_{i}=1\mid u_{i})=\frac{P(y_{i}=1)e^{-V\left(u_{i}\right)}}{P(y_{i}=1)e^{-V\left(u_{i}\right)}+P(y_{i}={-1})e^{-V\left(-u_{i}\right)}}.$

(5)

Choice of $P(y_{i}=1)$ does not influence posterior inference for $\{\beta_{0},{\bm{\beta}}\}$ in (4), but it is nevertheless important for the calibration of the predictive model and is discussed further in Section 6.1. Given $P(y_{i}=1)$, (5) defines a link function from $u_{i}\in\mathbb{R}$ to a probability in $[0,1]$; this function is smooth and comparable to a probit or logit link. Figure 1 plots this function under equal class proportions, $P(y_{i}=1)=1/2$.

Estimating $\beta_{0}$ and ${\bm{\beta}}$ via maximizing the likelihood under model (5) would perform poorly in the HDLSS context, analogous to the poor performance of unconstrained probit or logistic regression models. However, the posterior distribution (4) also implies a “prior” distribution for $\{\beta_{0},{\bm{\beta}}\}$ that is not conditional on $\mathbf{y}$:

	$\displaystyle p({\bm{\beta}},\beta_{0}\mid\mathbf{X},\mathbf{y},\lambda)$	$\displaystyle\propto P(\mathbf{y}\mid{\bm{\beta}},\beta_{0},\mathbf{X})\,p({\bm{\beta}},\beta_{0}\mid\mathbf{X},\lambda)$
	$\displaystyle\rightarrow p({\bm{\beta}},\beta_{0}\mid\mathbf{X},\lambda)$	$\displaystyle\propto A\cdot B$		(6)

where

$\displaystyle A=\prod_{i=1}^{n}\left(P(y_{i}=1)e^{-V\left(\beta_{0}+{\bf x}_{n}^{T}{\bm{\beta}}\right)}+P(y_{i}=-1)e^{-V\left(-\beta_{0}-{\bf x}_{n}^{T}{\bm{\beta}}\right)}\right)$

(7)

and

$\displaystyle B=e^{-(\lambda N/2)||{\bm{\beta}}||^{2}_{2}}.$

(8)

This prior facilitates shrinkage in two important ways that improve performance in HDLSS settings. Term $B$ (8) gives the kernel of a multivariate normal distribution for ${\bm{\beta}}$, in which the $\beta_{i}$’s are independent with mean $0$ and variance $1/(\lambda N)$. This derives from the $l_{2}$ penalty in (2). There is a vast literature on the connection between $l_{2}$-regularized regression (e.g., ridge regression) and using Gaussian priors on the coefficients in a Bayesian framework, as mentioned in Section 1. This shrinks the inferred coefficients toward $\mathbf{0}$, which improves performance in the HDLSS setting or in the presence of multicollinearity. Term $A$ (7) favors ${\bm{\beta}}$ that correspond to directions in $\mathbf{X}$ with high variability or bimodality. To illustrate, Figure 2 plots a single term of the product $A$ as a function of the DWD score $u_{i}=\beta_{0}+{\bf x}_{i}^{T}{\bm{\beta}}$ when $P(y_{i}=1)=1/2$,

$\displaystyle f(u_{i})=e^{-V(u_{i})}+e^{-V(-u_{i})}.$

(9)

This term gives higher prior density to directions in $\mathbf{X}$ for which the corresponding scores $\{u_{i}\}_{i=1}^{N}$ are not concentrated near $0$, and are therefore more useful for discriminating a negative class from a positive class.

The function $A\cdot B$ in (6) defines a proper density for ${\bm{\beta}}$ for any fixed intercept $\beta_{0}$, as stated in Theorem 2.

The prior (6) does not define a proper density for $\beta_{0}$. Improper priors are commonly used for a variety of Bayesian applications, and still provide a valid framework for inference if the posterior is proper (Taraldsen and Lindqvist, 2010; Sun et al., 2001). Nevertheless, it would be straightforward to modify the model to allow for shrinkage of $\beta_{0}$, if one desired it to have a proper probability density without conditioning on $\mathbf{y}$.

No explicit probabilistic assumptions on $\mathbf{X}$ are needed for the validity of the posterior and prior distributions in the previous two sections. However, it is informative to consider the broader implications of the model with respect to the distribution for $\mathbf{X}$. Without loss of generality we set $\beta_{0}=0$ and assume $P(y_{i}=1)=0.5$. Note that

$$p({\bm{\beta}}\mid{\bf x}_{i},\lambda)=\left(e^{-V\left({\bf x}_{i}^{T}{\bm{\beta}}\right)}+e^{-V\left(-{\bf x}_{i}^{T}{\bm{\beta}}\right)}\right)e^{-(\lambda/2)||{\bm{\beta}}||^{2}_{2}}/\phi({\bf x}_{i},\lambda),$$

where

$$\phi({\bf x}_{i},\lambda)=\int_{{\bm{\beta}}}\left(e^{-V\left({\bf x}_{i}^{T}{\bm{\beta}}\right)}+e^{-V\left(-{\bf x}_{i}^{T}{\bm{\beta}}\right)}\right)e^{-(\lambda/2)||{\bm{\beta}}||^{2}_{2}}d{\bm{\beta}},$$

Given any marginal density for ${\bf x}$ (which may depend on $\lambda$), ${\bf x}_{i}\sim p({\bf x}\mid\lambda)$,

	$\displaystyle p({\bm{\beta}},{\bf x}_{i}\mid\lambda)=p({\bm{\beta}}\mid{\bf x}_{i},\lambda)p({\bf x}_{i}),$
	$\displaystyle p({\bm{\beta}}\mid\lambda)=\int p({\bm{\beta}},{\bf x}_{i}\mid\lambda)d{\bf x}_{i},$
	$\displaystyle p({\bf x}_{i}\mid{\bm{\beta}},\lambda)=p({\bm{\beta}},{\bf x}_{i}\mid\lambda)/p({\bm{\beta}}\mid\lambda).$

Thus, a marginal prior for ${\bm{\beta}}$ exists but it depends on the distribution ${\bf x}_{i}$. Moreover,

$\displaystyle p({\bf x}_{i}\mid{\bm{\beta}},\lambda)=p({\bf x}_{i}\mid{\bm{\beta}},y_{i}=1,\lambda)+p({\bf x}_{i}\mid{\bm{\beta}},y_{i}=-1,\lambda)$

implies that

$$p({\bf x}_{i}\mid{\bm{\beta}},y_{i})\propto e^{-V(y_{i}{\bf x}_{i}^{T}{\bm{\beta}})}\phi({\bf x}_{i},\lambda)p({\bf x}_{i}\mid\lambda)/\phi({\bf x}_{i},\lambda).$$

Thus, given any marginal distribution for the data ${\bf x}_{i}$, it is possible to obtain the conditional distribution of ${\bf x}_{i}$ given its class membership $\mathbf{y}_{i}$ and the DWD vector ${\bm{\beta}}$.

The prior class probabilities $P(y_{i}=1)$ and $P(y_{i}=-1)$ in (5) do not affect the posterior for $\{\beta_{0},{\bm{\beta}}\}$ in (4). Thus, they need not be specified for tasks that do not require a predictive model for $y_{i}$, such as exploratory visualization or batch adjustment (Huang et al., 2012). However, their specification is important to obtain appropriately calibrated probabilities for classification. In particular, they must be considered carefully when the proportions in the training data differ from that of the target population. The ability to adjust the marginal class probabilities for out-of-sample prediction is useful, as the decision boundary for DWD ($\beta_{0}$) is sensitive to the class proportions used for training (Qiao et al., 2010; Qiao and Zhang, 2015). We let each observation have the same class membership probability a priori:

$$P(y_{i}=1)=P_{1}\;\text{ for }i=1,\ldots,n,$$

and by default we set $P_{1}=0.5$. In other situations $P_{1}$ may be known and fixed at another value, or estimated as the sample proportion from class $1$. Alternatively, one can give $P_{1}$ its own prior (e.g., $P_{1}\sim\text{Uniform}(0,1)$) and infer it with the rest of the Bayesian model; such an approach may be useful in (e.g.,) situation for which the $y_{i}$’s are only partially observed and the population class proportions are of interest.

Rather than fixing the penalty parameter $\lambda$ a priori, we recommend inferring its value within the Bayesian framework. For a given prior density $\lambda\sim p(\lambda)$, the conditional density for $\lambda$ is

$\displaystyle\begin{split}p(\lambda\mid\mathbf{X},\mathbf{y},{\bm{\beta}},\beta_{0})&\propto p(\lambda)P(\mathbf{y}\mid\mathbf{X},{\bm{\beta}},\beta_{0})p({\bm{\beta}},\mathbf{X}\mid\beta_{0},\lambda)\\ &\propto\prod_{i=1}^{n}p(\lambda)\left(P_{1}e^{-V\left(\beta_{0}+{\bf x}_{i}^{T}{\bm{\beta}}\right)}+(1-P_{1})e^{-V\left(-\beta_{0}-{\bf x}_{i}^{T}{\bm{\beta}}\right)}\right)/\phi(\mathbf{X},\beta_{0},\lambda)\end{split}$

(10)

where

$\displaystyle\phi(\mathbf{X},\lambda,\beta_{0})=\int_{{\bm{\beta}}}\prod_{i=1}^{n}\left(P_{1}e^{-V\left(\beta_{0}+{\bf x}_{i}^{T}{\bm{\beta}}\right)}+(1-P_{1})e^{-V\left(-\beta_{0}-{\bf x}_{i}^{T}{\bm{\beta}}\right)}\right)e^{-(\lambda/2)||{\bm{\beta}}||^{2}_{2}}d{\bm{\beta}}.$

(11)

In practice we use a prior for $\lambda$ that is uniform over a large range, by default, $\lambda\sim\text{Uniform}(1/128,128)$. The impact of $\lambda$ depends on the scale of $\mathbf{X}$ and the number of observations $n$, so the prior may need to be suitably modified for other scenarios.

The conditional distribution of ${\bm{\beta}}$ on $\mathbf{X}$ in (6) is not flat, and therefore observations for which $y_{i}$ are not observed can still inform ${\bm{\beta}}$. Consider the context in which $n_{o}$ observations are fully observed $\{({\bf x}_{i},\mathbf{y}_{i})\}_{i=1}^{n_{o}}$, $n_{u}$ have missing outcome $\{({\bf x}_{i}^{u},y_{i})\}_{i=1}^{n_{u}}$, and $n=n_{o}+n_{u}$. The full posterior distribution for $\{\beta_{0},{\bm{\beta}}\}$ is

	$\displaystyle p({\bm{\beta}},\beta_{0}\mid\mathbf{X},\mathbf{y},\lambda)$	$\displaystyle\propto\left[\prod_{i=1}^{n_{u}}\left(P_{1}e^{-V\left(\beta_{0}+{\bf x}_{n}^{T}{\bm{\beta}}\right)}+(1-P_{1})e^{-V\left(-\beta_{0}-{\bf x}_{n}^{T}{\bm{\beta}}\right)}\right)\right]$
		$\displaystyle\;\;\;\;\;\;\left[\prod_{i=1}^{n_{o}}e^{-V\left(y_{i}(\beta_{0}+{\bf x}_{i}^{T}{\bm{\beta}})\right)}\right]e^{-(\lambda n/2)||{\bm{\beta}}||^{2}_{2}}.$

This full posterior can be used for semi-supervised learning.

To simulate from the posterior distribution for $\{\beta,{\bm{\beta}}_{0}\}$ (4), we implement a “Metropolis-in-Gibbs” algorithm (Carlin and Louis, 2008). That is, we iteratively draw each $\beta_{i}$ from its conditional distribution given the current values of the rest of the parameters ${\bm{\beta}}[-i]$, $p(\beta_{i}\mid{\bm{\beta}}[-i],\mathbf{X},\mathbf{y},\lambda)$; each conditional draw is accomplished via a Metropolis sub-step. For the Metropolis step, we generate proposals from a normal distribution with mean given by the previous value for $\beta_{i}$ and a variance that is adaptively scaled over the course of the algorithm. We initialize the Markov chain at the posterior mode, which can be identified using existing software for DWD; we use the sdwd package for this initialization.

Inference for $\lambda$ can be accomplished by iteratively simulating from its full conditional distribution (10) to update its value within the Gibbs sampling algorithm in Section 7.1. However, this is not straightforward, because the integral $\phi(\mathbf{X},\beta_{0},\lambda)$ (11) does not have a closed form. In our implementation, we first estimate $\phi(\mathbf{X},\lambda,\beta_{0})$ over a grid of values $\{\lambda_{j}:j=1,\ldots,J\}$ that span the range of the prior for $\lambda$. For this, we utilize the fact that $e^{-(\lambda n/2)||{\bm{\beta}}||^{2}_{2}}$ resembles a Gaussian kernel to approximate the integral via Monte Carlo with simulated random normal vectors. Given ${\bm{\beta}}^{(t)}\sim\mbox{MVN}(\mathbf{0},1/(\lambda n)\mathbf{I})$ for $t=1,\ldots,T$,

$$\phi(\mathbf{X},\lambda,\beta_{0})\approx\left(\frac{2\pi}{n\lambda}\right)^{p/2}\frac{1}{T}\sum_{t=1}^{T}\prod_{i=1}^{n}\left(P_{1}e^{-V\left(\beta_{0}+{\bf x}_{i}^{T}{\bm{\beta}}^{(t)}\right)}+(1-P_{1})e^{-V\left(-\beta_{0}-{\bf x}_{i}^{T}{\bm{\beta}}^{(t)}\right)}\right).$$

Using this approximation, we estimate $\hat{\phi}(\mathbf{X},\lambda_{j},\beta_{0})$ for each $j$, and then linearly interpolate on the log-scale to approximate the full function $\hat{\phi}(\mathbf{X},\lambda,\beta_{0})$ over the range of the prior for $\lambda$. Posterior estimation then proceeds as in Section 7.1, with an additional Metropolis sub-step to update the value of $\lambda$. The full algorithm for posterior sampling is given in Appendix B.

In practice, we find that $1000$ MCMC cycles yield appropriate coverage of the posterior when $\lambda$ is fixed (see Section 8.1) and $10000$ are sufficient when $\lambda$ is inferred. On a single CPU, $10000$ MCMC iterations takes approximately $5$ minutes for a dataset of size $\mathbf{X}:500\times 500$ (comparable to the data considered in Section 9) and computing time scales linearly with the dimension $d$. Thus, full posterior inference via MCMC is feasible for most applications (with patience), but computing time can be a limitation.

Here we fix $\beta_{0}=0$ and simulate data from the assumed model as follows:

1. 1.

   Generate $\mathbf{X}:d\times n$ according to a specified probability distribution,

2. 2.

   Generate ${\bm{\beta}}_{\text{true}}$ from its conditional prior $\eqref{bprior}$ (via a Metropolis-Hastings algorithm), given $\lambda$

3. 3.

   For $i=1,\ldots,n$, compute the score $u_{i}={\bf x}_{i}^{T}{\bm{\beta}}_{\text{true}}$ and generate $y_{i}$ according to its class probabilities defined by (5).

As manipulated conditions, we consider three values for $\lambda$ ($\lambda=0.1,1,$ or $10$), two configurations of $d$ and $n$ ($\mathbf{X}:20\times 100$ or $\mathbf{X}:100\times 20$), and three different probability distributions for $\mathbf{X}$: (1) entries are simulated independently from a Uniform$(-1,1)$ distribution, (2) entries are simulated independently from an Exponential$(1)$ distribution centered to have mean $0$, and (3) each column ${\bf x}_{i}$ is simulated from the following mixture of multivariate normal distributions:

$\displaystyle x_{i}\sim\begin{cases}\mbox{MVN}(\bm{\mu}_{0},\mathbf{I})\text{ with probability }1/2\\ \mbox{MVN}(\bm{\mu}_{1},\mathbf{I})\text{ with probability }1/2,\end{cases}$

(12)

where the entries of $\bm{\mu}_{0}$ and $\bm{\mu}_{1}$ are generated independently from a N$(0,0.5)$ distribution.

For each set of manipulated conditions, we generate $100$ datasets $\{\mathbf{X},{\bm{\beta}}_{\text{true}},\mathbf{y}\}$ under the assumed model and approximate the posterior distribution for ${\bm{\beta}}$ using the MCMC algorithm in Section 7.1. We construct 95% Bayesian credible intervals for ${\bm{\beta}}_{\text{true}}$ and the scores $\mathbf{u}=\mathbf{X}^{T}{\bm{\beta}}_{\text{true}}$ based on the posterior samples. We also consider credible intervals constructed using the posterior mode and the asymptotic approximation in Theorem 3 (i.e., the CLT approximation), and confidence intervals based on the percentiles of $200$ non-parametric bootstrap resamples as in Lyu et al. (2017).

Table 1 gives the average coverage rates for the true values under each set of conditions and the three different approaches to inference. The credible intervals using MCMC all have nominal coverage rates, which validates the algorithm. The CLT approximation gives appropriate coverage in several instances, even when $p>n$. However, it performs less well when $\lambda$ is smaller and has relatively lower coverage rates for the scores $\mathbf{u}=\mathbf{X}^{T}{\bm{\beta}}_{\text{true}}$. The bootstrap percentile approach yields more narrow intervals and tends to undercover in all cases. As a representative example, Figure 3 plots the mode and 95% intervals for the coefficients for an instance under the Uniform scenario with $\lambda=0.1$, $n=100$ and $p=20$; this shows close agreement between the MCMC and CLT intervals, and the bootstrap intervals are universally more narrow.

To assess the calibration of estimated class probabilities, we generate a new set of $1000$ observations under the true model as a “test” set for each instance of the simulation. We compute the class probabilities for each test case, given the posterior fit to the $n=20$ or $n=100$ training cases. We group the estimated probabilities into narrow bins of length $0.02$ ($[0,0.02),[0.02,0.04)$, etc.) and consider the proportion of cases with $y=1$ in each bin. Figure 4 plots these proportions, aggregated across all conditions, vs. the midpoint of each bin. Using the posterior mean, the empirical proportions match the estimated probabilities perfectly, confirming that the model is well-calibrated and again validating the MCMC procedure. Using the posterior mode, the relationship deviates slightly but still tends to provide a reasonable fit in most cases.

Here we generate data according to a two-class multivariate normal distribution, in which each class has a different mean vector, as follows:

1. 1.

   Generate $d$-dimensional mean vectors $\bm{\mu}_{0}$ and $\bm{\mu}_{1}$, with independent entries from a N$(0,\tau^{2})$ distribution.

2. 2.

   Set $y_{i}=-1$ for $i=1,\ldots,n/2$ and $y_{i}=1$ for $i=n/2+1,\ldots,n$.

3. 3.

   For $i=1,\ldots,n$, generate ${\bf x}_{i}$ via

   $$x_{i}\sim\begin{cases}\mbox{MVN}(\bm{\mu}_{0},\mathbf{I})\text{ if }y_{0}=-1\\ \mbox{MVN}(\bm{\mu}_{1},\mathbf{I})\text{ if }y_{i}=1.\end{cases}$$

Note that this scenario differs from the ‘Bimodal’ case in Section 8.1 (12); for that scenario, membership in a given mixture component (defined by $\bm{\mu}_{0}$ and $\bm{\mu}_{1}$) does not necessarily correspond to the supervised class labels $y_{i}$. The present scenario does not precisely match our assumed model. The “true” (i.e., oracle) probabilistic link between $y_{i}$ and ${\bf x}_{i}$ in this scenario corresponds to a probit link with coefficients that are proportional to the mean difference vector $\bm{\mu}_{0}-\bm{\mu}_{1}$; as shown in Figure 1 the link for Bayesian DWD is not a probit, but it is close.

We fix $n=100$ and consider different levels of signal ($\tau=0,0.1,0.2,0.5$) and different dimensions ($d=10,100,1000$) as manipulated conditions. For each instance, we approximate the posterior separately with $\lambda$ fixed at each of $15$ possible values, that range from $1/128$ to $128$ and are equally spaced on a log scale. For each case we compute the Kullback-Leibler divergence (KL-divergence) of the estimated class probabilities from the oracle model. For each case we also generate a test set of $5000$ samples and evaluate (1) the misclassification rate on the test samples using a posterior predictive probability of $0.5$ as a threshold, and the mean squared error (MSE) of the predicted class probabilities for class 1 ($p_{i}$) from the true class memberships as

$$\text{MSE}=\frac{1}{n_{\text{test}}}\sum_{i=1}^{n_{\text{test}}}(1-p_{i})^{2}\mathbbm{1}_{\{y_{i}=-1\}}+p_{i}^{2}\mathbbm{1}_{\{y_{i}=1\}}.$$

Figure 5 plots these metrics as a function of $\lambda$ for the different conditions. In general, lower values of $\lambda$ perform better when the distinguishing signal is stronger. This is intuitive, as $\lambda$ can be considered an inverse scale parameter for ${\bm{\beta}}$, and larger magnitudes for ${\bm{\beta}}$ correspond to probabilities closer to $0$ and $1$. The misclassification rate is robust to the value of $\lambda$, but its specification is important to achieve well-calibrated probabilities.

Motivated by the sensitivity of the estimated probabilities to $\lambda$, we also approximate the posterior with a uniform prior on $\lambda$ as in Section 6.2. Table 2 summarizes the posterior mean for $\lambda$ under each set of conditions, and compares the performance of the model with uniform prior on $\lambda$ with that of the best performing fixed value of $\lambda$ for each metric. The model with a uniform prior performs comparatively well across the different scenarios, and the posterior mean for $\lambda$ appropriately ranges from very high values when the signal is nonexistent to low values when the signal is strong. However, for the conditions with $p=1000$ and strong signal the posterior mean is substantially larger than the best performing values of $\lambda$; this is likely because the model fits well and perfectly classifies the data with probabilities close to $0$ or $1$ for a wide range of $\lambda$ values (see Figure 5), so there is not strong evidence to distinguish between different values of $\lambda$ in this range.

In Figure 6 we consider the overall calibration of the fitted probabilities across the different conditions, as in Section 8.1. This shows clearly that fixing $\lambda$ at a large value ($\lambda=128$) can yield overly conservative estimated probabilities and a relatively small value ($\lambda=1/128$) can yield overly confident probabilities. The model in which $\lambda$ is inferred with a uniform prior is mostly well-calibrated but does show some deviations from the estimated probabilities and empirical proportions. These deviations are likely in part due to the fact that the assumed probabilistic link does not perfectly match the true generative model.

Here we generate data in which some of the class labels $y_{i}$ may be unknown, to illustrate a semi-supervised approach. We consider two scenarios. For the first scenario, the observations ${\bf x}_{i}$ are generated exactly as in Section 8.2 with $d=100$ and $\tau=0.3$. For the second scenario, $d=100$ and entries of each ${\bf x}_{i}$ are generated independently from a $N(0,1)$ distribution; then class labels $y_{i}$ are defined deterministically by the sign of ${\bf x}_{i}^{T}{\bm{\beta}}_{\text{sep}}$, where the entries of ${\bm{\beta}}_{\text{sep}}$ are generated from a standard normal distribution. The key difference between these two scenarios, for our purpose, is that the first scenario (bimodal) has an underlying cluster structure that distinguishes the two classes whereas the second (unimodal) does not. So, we expect to improve performance by considering unlabeled data for the bimodal scenario but not the unimodal scenario.

As manipulated conditions, we vary the number of observations with $y_{i}$ observed ($n_{o}=10,20,100,1000)$ and with $y_{i}$ unobserved ($n_{u}=0,10,20,100,1000$) in our training set, for each of the bimodal and unimodal scenarios. For each combination of conditions, we generate $20$ datasets, use MCMC to approximate the posterior, and consider the resulting misclassification rates for a test set. The average test misclassification rates under each set of conditions are shown in Table 3. This shows that for the bimodal scenario a semi-supervised analysis with a large number of unlabeled observations can improve performance dramatically; the misclassification rates for $n_{o}=10$ observed and $n_{u}=1000$ unobserved are comparable to that with $n_{o}=1000$ observed. As expected, including data without observed class labels does not help in the unimodal scenario; however, it does not lead to a substantial decrease in performance.