Tractable Bayesian Density Regression via Logit Stick-Breaking Priors

TL;DR

This paper introduces a computationally efficient Bayesian density regression method using logit stick-breaking priors, enabling scalable inference with various algorithms and practical application in toxicology.

Contribution

It proposes a new predictor-dependent mixture model with a simple logistic representation and Pólya-gamma augmentation, improving computational efficiency over existing methods.

Findings

Algorithms successfully applied to a toxicology dataset.

Scalable inference methods demonstrated with Gibbs sampling, EM, and variational Bayes.

Enhanced practical implementation of Bayesian density regression.

Abstract

There is a growing interest in learning how the distribution of a response variable changes with a set of predictors. Bayesian nonparametric dependent mixture models provide a flexible approach to address this goal. However, several formulations require computationally demanding algorithms for posterior inference. Motivated by this issue, we study a class of predictor-dependent infinite mixture models, which relies on a simple representation of the stick-breaking prior via sequential logistic regressions. This formulation maintains the same desirable properties of popular predictor-dependent stick-breaking priors, and leverages a recent P\'olya-gamma data augmentation to facilitate the implementation of several computational methods for posterior inference. These routines include Markov chain Monte Carlo via Gibbs sampling, expectation-maximization algorithms, and mean-field variational Bayes for scalable inference, thereby stimulating a wider implementation of Bayesian density regression by practitioners. The algorithms associated with these methods are presented in detail and tested in a toxicology study.

Equations

$$f(y)=\int K(y;\mathbf{\theta })p(\text{mbox}d\mathbf{\theta })=\sum _{h=1}^{\infty }{\pi }_{h}K(y;{\mathbf{\theta }}_h),$$

$$f_{\mathbf{x}}(y)=\int K_{\mathbf{x}}(y;\mathbf{\theta })p_{\mathbf{x}}(\text{mbox}d\mathbf{\theta })=\sum _{h=1}^{\infty }{\pi }_h(\mathbf{x})K_{\mathbf{x}}(y;{\mathbf{\theta }}_h),$$

$$(y_i\mid G_i=h,{\mathbf{x}}_i)\sim K_{{\mathbf{x}}_i}(y_i;{\mathbf{\theta }}_h),\text{mbox}with\text{mbox}pr(G_i=h\mid {\mathbf{x}}_i)={\pi }_h({\mathbf{x}}_i)={\nu }_h({\mathbf{x}}_i)\prod _{l=1}^{h-1}\{1-{\nu }_l({\mathbf{x}}_i)\},$$

$${\nu }_h({\mathbf{x}}_i)=\frac{{\pi }_h({\mathbf{x}}_i)}{1-{\sum }_{l=1}^{h-1}{\pi }_l({\mathbf{x}}_i)}=\frac{\text{mbox}pr(G_i=h\mid {\mathbf{x}}_i)}{\text{mbox}pr(G_i>h-1\mid {\mathbf{x}}_i)},h\in \mathbb{N},$$

$${\zeta }_{ih}=z_{ih}\prod _{l=1}^{h-1}(1-z_{il}),h\in \mathbb{N},$$

$${\pi }_h({\mathbf{x}}_i)=\frac{\exp \{{\eta }_h({\mathbf{x}}_i)\}}{1+\exp \{{\eta }_h({\mathbf{x}}_i)\}}\prod _{l=1}^{h-1}\left[\frac{1}{1+\exp \{{\eta }_l({\mathbf{x}}_i)\}}\right],h\in \mathbb{N},$$

$${\eta }_h({\mathbf{x}}_i)=\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\alpha }}_h,\text{mbox}with{\mathbf{\alpha }}_h\sim \text{mbox}{N}_R({\mathbf{\mu }}_{\alpha },{\mathbf{\Sigma }}_{\alpha }),h\in \mathbb{N}.$$

$$f_{\mathbf{x}}(y)=\int \sqrt{\tau }\varphi [\sqrt{\tau }\left\{y-\mathbf{\lambda }(\mathbf{x}{)}^{⊺}\mathbf{\beta }\right\}]p_{\mathbf{x}}(\text{mbox}d\mathbf{\beta },\text{mbox}d\tau )=\sum _{h=1}^{\infty }{\pi }_h(\mathbf{x})\sqrt{{\tau }_h}\varphi [\sqrt{{\tau }_h}\left\{y-\mathbf{\lambda }(\mathbf{x}{)}^{⊺}{\mathbf{\beta }}_h\right\}],$$

$$f_{\mathbf{X}}^{(H)}(\mathbf{y})=\text{mbox}E\left(\prod _{i=1}^n\sum _{h=1}^H{\pi }_h({\mathbf{x}}_i)\sqrt{{\tau }_h}\varphi [\sqrt{{\tau }_h}\left\{y_i-\mathbf{\lambda }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\beta }}_h\right\}]\right),$$

$$f(\mathbf{\alpha })f(\mathbf{\beta })f(\mathbf{\tau })\prod _{i=1}^n\left[\prod _{h=1}^H(\sqrt{{\tau }_h}\varphi [\sqrt{{\tau }_h}\{y_i-\mathbf{\lambda }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\beta }}_h\}]{)}^{\text{mathbbm}1(G_i=h)}\prod _{h=1}^{H-1}{\nu }_h({\mathbf{x}}_i{)}^{\text{mathbbm}1(G_i=h)}\{1-{\nu }_h({\mathbf{x}}_i){\}}^{\text{mathbbm}1(G_i>h)}\right],$$

$$f_{{\mathbf{x}}_i}(z_{ih})=\frac{0.5\exp \{(z_{ih}-0.5)\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\alpha }}_h\}}{\text{mbox}cosh\{0.5\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\alpha }}_h\}},f_{{\mathbf{x}}_i}({\omega }_{ih})=\frac{\exp [-0.5\{\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\alpha }}_h{\}}^2{\omega }_{ih}]f({\omega }_{ih})}{[\text{mbox}cosh\{0.5\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\alpha }}_h\}{]}^{-1}},$$

$$\text{mbox}pr(G_i=h\mid -)=\frac{\left[{\nu }_h({\mathbf{x}}_i){\prod }_{l=1}^{h-1}\{1-{\nu }_l({\mathbf{x}}_i)\}\right]\sqrt{{\tau }_h}\varphi [\sqrt{{\tau }_h}\{y_i-\mathbf{\lambda }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\beta }}_h\}]}{{\sum }_{q=1}^H\left[{\nu }_q({\mathbf{x}}_i){\prod }_{l=1}^{q-1}\{1-{\nu }_l({\mathbf{x}}_i)\}\right]\sqrt{{\tau }_q}\varphi [\sqrt{{\tau }_q}\{y_i-\mathbf{\lambda }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\beta }}_q\}]},$$

$${\hat{\zeta }}_{ih}=\frac{\left[{\nu }_h^{(t)}({\mathbf{x}}_i){\prod }_{l=1}^{h-1}\{1-{\nu }_l^{(t)}({\mathbf{x}}_i)\}\right]\sqrt{{\tau }_h^{(t)}}\varphi [\sqrt{{\tau }_h^{(t)}}\{y_i-\mathbf{\lambda }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\beta }}_h^{(t)}\}]}{{\sum }_{q=1}^H\left[{\nu }_q^{(t)}({\mathbf{x}}_i){\prod }_{l=1}^{q-1}\{1-{\nu }_l^{(t)}({\mathbf{x}}_i)\}\right]\sqrt{{\tau }_q^{(t)}}\varphi [\sqrt{{\tau }_q^{(t)}}\{y_i-\mathbf{\lambda }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\beta }}_q^{(t)}\}]},$$

$${\mathbf{\alpha }}_h^{(t+1)}=\{\mathbf{\Psi }(\mathbf{x}{)}^{⊺}\text{mbox}diag({\hat{\bar{\omega }}}_{1h},\dots ,{\hat{\bar{\omega }}}_{nh})\mathbf{\Psi }(\mathbf{x})+{\mathbf{\Sigma }}_{\mathbf{\alpha }}^{-1}{\}}^{-1}\{\mathbf{\Psi }(\mathbf{x}{)}^{⊺}({\hat{\bar{\kappa }}}_{1h},\dots ,{\hat{\bar{\kappa }}}_{nh}{)}^{⊺}+{\mathbf{\Sigma }}_{\mathbf{\alpha }}^{-1}{\mathbf{\mu }}_{\mathbf{\alpha }}\},$$

$${\mathbf{\beta }}_h^{(t+1)}$$

$${\tau }_h^{(t+1)}$$

$$\sum _{i=1}^n{\ell }_{{\mathbf{x}}_i}(\mathbf{\alpha },\mathbf{\beta },\mathbf{\tau };y_i,{\mathbf{\zeta }}_i,{\bar{\mathbf{\omega }}}_i)+\sum _{h=1}^{H-1}\log f({\mathbf{\alpha }}_h)+\sum _{h=1}^H\log f({\mathbf{\beta }}_h)+\sum _{h=1}^H\log f({\tau }_h)+\text{mbox}const,$$

$$\sum _{h=1}^H{\zeta }_{ih}\left[-\frac{{\tau }_h\{y_i-\mathbf{\lambda }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\beta }}_h{\}}^2}{2}+\frac{1}{2}\log ({\tau }_h)\right]+\sum _{h=1}^{H-1}\left[{\bar{\kappa }}_{ih}\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\alpha }}_h-{\bar{\omega }}_{ih}\frac{\{\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\alpha }}_h{\}}^2}{2}\right]+\text{mbox}const,$$

$$\text{mbox}logit({\rho }_{ih})=\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}\text{mbox}E({\mathbf{\alpha }}_h)+\sum _{l=h}^H{\zeta }_{il}^{(h)}\left[0.5\cdot \text{mbox}E(\log {\tau }_l)-0.5\cdot \text{mbox}E({\tau }_l)\text{mbox}E\{(y_i-\mathbf{\lambda }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\beta }}_l{)}^2\}\right],$$

$$f(\mathbf{\alpha })f(\mathbf{\beta })f(\mathbf{\tau })\prod _{i=1}^n\left[\prod _{h=1}^H(\sqrt{{\tau }_h}\varphi [\sqrt{{\tau }_h}\{y_i-\mathbf{\lambda }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\beta }}_h\}]{)}^{z_{ih}{\prod }_{l=1}^{h-1}(1-z_{il})}\prod _{h=1}^{H-1}\frac{f({\omega }_{ih})}{2}\frac{\exp \{(z_{ih}-0.5)\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\alpha }}_h\}}{\exp \{0.5{\omega }_{ih}(\mathbf{\psi }({\mathbf{x}}_i{)}^{⊺}{\mathbf{\alpha }}_h{)}^2\}}\right],$$

$$\sum _{\mathbf{z}}\int q_{\mathbf{x}}(\mathbf{\alpha },\mathbf{\beta },\mathbf{\tau },\mathbf{z},\mathbf{\omega })\left[\log \left\{\frac{f_{\mathbf{x}}(\mathbf{y}\mid \mathbf{z},\mathbf{\beta },\mathbf{\tau })f_{\mathbf{x}}(\mathbf{z}\mid \mathbf{\alpha })f_{\mathbf{x}}(\mathbf{\omega }\mid \mathbf{\alpha })f(\mathbf{\alpha })f(\mathbf{\beta })f(\mathbf{\tau })}{q_{\mathbf{x}}(\mathbf{\alpha },\mathbf{\beta },\mathbf{\tau },\mathbf{z},\mathbf{\omega })}\right\}\right]\text{mbox}d(\mathbf{\alpha },\mathbf{\beta },\mathbf{\tau },\mathbf{\omega }).$$

$$q_{\mathbf{x}}(\mathbf{\alpha },\mathbf{\beta },\mathbf{\tau },\mathbf{z},\mathbf{\omega })=\prod _{h=1}^{H-1}{q}_{\mathbf{x}}({\mathbf{\alpha }}_h)\prod _{h=1}^H{q}_{\mathbf{x}}({\mathbf{\beta }}_h)\prod _{h=1}^H{q}_{\mathbf{x}}({\tau }_h)\prod _{h=1}^{H-1}\prod _{i=1}^n{q}_{{\mathbf{x}}_i}(z_{ih})\prod _{h=1}^{H-1}\prod _{i=1}^n{q}_{{\mathbf{x}}_i}({\omega }_{ih}).$$

$$\begin{array}{rl}& \log q_{\mathbf{x}}^{\ast }({\mathbf{\beta }}_h)=\text{mbox}{E}_{\mathbf{\tau },\mathbf{z}}[\log \{f_{\mathbf{x}}(\mathbf{y}\mid \mathbf{z},\mathbf{\beta },\mathbf{\tau })f({\mathbf{\beta }}_h)\}]+c_{{\mathbf{\beta }}_h},h=1,\dots ,H,\\ & \log q_{\mathbf{x}}^{\ast }({\tau }_h)=\text{mbox}{E}_{\mathbf{\beta },\mathbf{z}}[\log \{f_{\mathbf{x}}(\mathbf{y}\mid \mathbf{z},\mathbf{\beta },\mathbf{\tau })f({\tau }_h)\}]+c_{{\tau }_h},h=1,\dots ,H,\\ & \log q_{\mathbf{x}}^{\ast }({\mathbf{\alpha }}_h)=\text{mbox}{E}_{\mathbf{z},\mathbf{\omega }}[\log \{f_{\mathbf{x}}(\mathbf{z},\mathbf{\omega }\mid \mathbf{\alpha })f({\mathbf{\alpha }}_h)\}]+c_{{\mathbf{\alpha }}_h},h=1,\dots ,H-1,\\ & \log q_{{\mathbf{x}}_i}^{\ast }(z_{ih})=\text{mbox}{E}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\tau },{\mathbf{z}}_{i,-h}}[\log f_{\mathbf{x}}(\mathbf{y},\mathbf{z}\mid \mathbf{\beta },\mathbf{\tau },\mathbf{\alpha })]+c_{z_{ih}},h=1,\dots ,H-1,i=1,\dots ,n,\\ & \log q_{{\mathbf{x}}_i}^{\ast }({\omega }_{ih})=\text{mbox}{E}_{\mathbf{\alpha }}[\log f_{\mathbf{x}}({\omega }_{ih}\mid \mathbf{\alpha })]+c_{{\omega }_{ih}},h=1,\dots ,H-1,i=1,\dots ,n,\end{array}$$

$$||f_{\mathbf{X}}^{(H)}(\mathbf{y})-f_{\mathbf{X}}^{(\infty )}(\mathbf{y})|{|}_1\le 4\left[1-\text{mbox}E\left\{\prod _{i=1}^n\sum _{h=1}^{H-1}{\pi }_h({\mathbf{x}}_i)\right\}\right]=4\cdot \text{mbox}E\left[1-\prod _{i=1}^n\sum _{h=1}^{H-1}{\pi }_h({\mathbf{x}}_i)\right].$$

Full text

1 Introduction

There is a growing interest in density regression models which allow the entire distribution of a univariate response variable $y\in\mathcal{Y}$ to be unknown and changing with a vector of predictors $\bf{x}\in\mathcal{X}$. Indeed, the increased flexibility provided by these procedures allows improvements in inference and prediction compared to classical regression frameworks, as seen in applications [e.g. Dunson and Park, 2008, Griffin and Steel, 2011, Wade et al., 2014].

Within the Bayesian nonparametric framework, there is a wide variety of alternative methods to provide flexible inference for conditional distributions. Most of these strategies represent generalizations of the marginal density estimation problem for $f(y)$, which is commonly addressed via Bayesian nonparametric mixture models of the form $f(y)=\int K(y;\bm{\theta})p(\mbox{d}\bm{\theta})$, where $K(y;\bm{\theta})$ is a known parametric kernel indexed by $\bm{\theta}\in\bm{\Theta}$, and $p$ is a random probability measure which is assigned a flexible prior $\Pi$. Popular choices for $\Pi$ are the Dirichlet process [Ferguson, 1973], the two-parameter Poisson-Dirichlet process [Pitman and Yor, 1997], and other random measures having a stick-breaking representation [Ishwaran and James, 2001]. This choice leads to the infinite mixture model

[TABLE]

where $\pi_{h}=\nu_{h}\prod_{l=1}^{h-1}(1-\nu_{l})$ for every $h\geq 1$, with $\pi_{1}=\nu_{1}$. In equation (1), the kernel parameters $(\bm{\theta}_{h})_{h\geq 1}$ are distributed according to a diffuse base measure $P_{0}$, whereas the stick-breaking weights $(\nu_{h})_{h\geq 1}$ have independent $\mbox{Beta}(a_{h},b_{h})$ priors, so that $\sum_{h=1}^{\infty}\pi_{h}=1$ almost surely.

Model (1) has key computational benefits in allowing the implementation of simple Markov chain Monte Carlo methods for posterior inference [e.g. Escobar and West, 1995, Neal, 2000], and provides a consistent strategy for density estimation [e.g. Ghosal et al., 1999, Tokdar, 2006, Ghosal and Van Der Vaart, 2007]. This has motivated different generalizations of (1) to incorporate the conditional density inference problem for $f(y\mid{\bf{x}})=f_{{\bf x}}(y)$, by allowing the random mixing measure $p_{\bf{x}}$ to change with $\bf{x}\in\mathcal{X}$, under a dependent stick-breaking characterization [MacEachern, 1999, 2000]. Popular representations consider predictor-independent mixing weights $\pi_{h}$, and incorporate changes with $\bf{x}\in\mathcal{X}$ in the atoms $\bm{\theta}_{h}(\bf{x})$; see for instance De Iorio et al. [2004], Gelfand et al. [2005], De la Cruz-Mesía et al. [2007]. As noted in MacEachern [2000] and Griffin and Steel [2006], the predictor-independent assumption for the mixing weights might have limited flexibility in practice. This has motivated more general formulations allowing also $\pi_{h}(\bf{x})$ to change with the predictors. Relevant examples include the order-based dependent Dirichlet process [Griffin and Steel, 2006], the kernel stick-breaking process [Dunson and Park, 2008], the infinite mixture model with predictor-dependent weights [Antoniano-Villalobos et al., 2014], and more recent representations for Bayesian dynamic inference [Gutiérrez et al., 2016]. These formulations provide a broader class of priors for Bayesian density regression, but their flexibility comes at a computational cost. In particular, the availability of simple algorithms for tractable posterior inference is limited by the specific construction of these representations.

The above issues motivate alternative formulations which preserve theoretical properties, but facilitate tractable posterior computation under a broader variety of algorithms. We aim to address this goal via a logit stick-breaking prior (lsbp), which relates each stick-breaking weight $\nu_{h}({\bf{x}})\in(0,1)$ to a function $\eta_{h}(\bf{x})\in\Re$ of the covariates, using the logit link. The proposed formulation is closely related to the probit stick-breaking prior (psbp) of Rodriguez and Dunson [2011]. Indeed, as we will discuss in Section 2, both lsbp and psbp are characterized by a continuation-ratio representation [Tutz, 1991], which allows to express the underlying clustering assignment in terms of independent and sequential binary regressions. This representation has key computational benefits and has been exploited by Rodriguez and Dunson [2011] to derive a Markov chain Monte Carlo (mcmc) algorithm for posterior inference. However, while the mcmc for psbp relies on the truncated Gaussian data augmentation for probit regression [Albert and Chib, 1993], the one for lsbp exploits the recent Pólya-gamma data augmentation for logistic regression [Polson et al., 2013], which might improve mixing compared to the psbp, especially in imbalanced situations [Johndrow et al., 2019]. As we will clarify in Section 2, these imbalanced settings can also occur in our case, since the binary regressions are associated to latent clustering allocations. We illustrate the mcmc algorithm for the lsbp in Section 3.

Besides developing tractable Gibbs sampling methods, we further derive alternative computational routines which address the scalability and mixing issues of mcmc in high-dimensional studies. Specifically, in Section 3 we illustrate a tractable expectation-maximization (em) routine for point estimation, and a simple variational Bayes (vb) algorithm for scalable inference. Both strategies leverage again the sequential representation of the lsbp and the associated Pólya-gamma data augmentation. Note that a vb routine for lsbp is also presented in Ren et al. [2011], but it is based on the bound of Jaakkola and Jordan [2000]. As a consequence of the recent theoretical findings in Durante and Rigon [2019], it can be shown that our approach is intimately related to the one of Ren et al. [2011], although being developed by means of seemingly unrelated strategies. Finally, while tractable algorithms such as em or vb could be possibly obtained also for psbp, we are not aware of any actual discussion or implementation. Indeed, the analytical derivations might be slightly more complex in the psbp case compared to the lsbp, as discussed in Section 3.

We shall emphasize that the overarching focus of our contribution is not on developing a novel methodological framework for Bayesian density regression, but on deriving a broad set of routine-use computational strategies under a suitable and tractable representation. To our knowledge this goal remains partially unaddressed, but represents a fundamental condition to facilitate routine implementation of Bayesian density regression by practitioners. The three proposed algorithms are empirically compared in Section 4 using a real data toxicology study, previously considered in Dunson and Park [2008]. Section 5 provides concluding remarks.

2 Logit stick-breaking prior

This section presents a formal construction of the lsbp via continuation-ratio logistic regressions. As a natural extension of model (1), we consider the general class of predictor-dependent infinite mixture models

[TABLE]

where $\pi_{h}({\bf x})=\nu_{h}({\bf x})\prod_{l=1}^{h-1}\{1-\nu_{l}({\bf x})\}$ are predictor-dependent mixing probabilities having a stick-breaking representation, whereas $K_{{\bf x}}(y;\bm{\theta})$ denotes a predictor-dependent kernel, indexed by parameters $\bm{\theta}$ and covariates ${\bf x}$.

To highlight the continuation-ratio representation of the lsbp, let us first consider an equivalent formulation of the predictor-dependent mixture model in (2). In particular, following standard hierarchical representations of mixture models, independent samples $y_{1},\ldots,y_{n}$ of the variable with density function displayed in (2), can be obtained from

[TABLE]

for each unit $i=1,\ldots,n$, where $\bm{\theta}_{h}\sim P_{0}$ independently for $h\in\mathbb{N}$, whereas $G_{i}\in\mathbb{N}$ is the categorical variable denoting the mixture component associated with the $i$th unit. According to (3), every $G_{i}$ has probability mass function $f(G_{i}\mid{\bf x}_{i})=\prod_{h=1}^{\infty}\pi_{h}({\bf x}_{i})^{\mathbbm{1}(G_{i}=h)}$, where $\mathbbm{1}(\cdot)$ denotes the indicator function. Hence, re-writing $\{\nu_{h}({\bf x}_{i})\}_{h\in\mathbb{N}}$ as a function of the mixing probabilities $\{\pi_{h}({\bf x}_{i})\}_{h\in\mathbb{N}}$ via

[TABLE]

allows to interpret each $\nu_{h}({\bf x}_{i})$ as the probability of being allocated to component $h$, conditionally on the event of surviving to the previous $1,\ldots,h-1$ components, namely $\nu_{h}({\bf x}_{i})=\mbox{pr}(G_{i}=h\mid G_{i}>h-1,{\bf x}_{i})$. This result provides a formal characterization of the stick-breaking construction in (3) as the continuation-ratio parameterization [Tutz, 1991] of the probability mass function for each component membership variable $G_{i}$. This connection with the literature on sequential inference for categorical data is common to all the stick-breaking priors—as mentioned also by Rodriguez and Dunson [2011] in the probit case.

As we will describe in Section 3, the above result facilitates the implementation of different routine-use algorithms in Bayesian inference, and provides a simple generative process for each $G_{i}$. In particular, as illustrated in Figure 1, in the first step of this continuation-ratio generative mechanism, unit $i$ is either assigned to the first component with probability $\nu_{1}({\bf{x}}_{i})$ or to one of the others with complement probability. If $G_{i}=1$ the process stops, otherwise it continues considering the reduced set $\{h:h>1\}$. A generic step $h$ is reached if $i$ has not been assigned to $1,\ldots,h-1$, and the decision at this step will be to either allocate $i$ to component $h$ with probability $\nu_{h}({\bf{x}}_{i})$, or to one of the subsequent components with probability $1-\nu_{h}({\bf{x}}_{i})$, conditioned on $G_{i}>h-1$. Based on this representation, the assignment indicator $\zeta_{ih}=\mathbbm{1}(G_{i}=h)$ can be expressed, for every unit $i=1,\ldots,n$, as

[TABLE]

where the generic $z_{ih}$, $h\in\mathbb{N}$, is a Bernoulli variable $(z_{ih}\mid{\bf{x}}_{i})\sim\mbox{Bern}\{\nu_{h}({\bf{x}}_{i})\}$ denoting the decision at the $h$th step to either allocate $i$ to component $h$ or to one of the subsequents. Hence, according to (5), the sampling of each $G_{i}$, under the predictor-dependent stick-breaking representation for each $\pi_{h}({\bf{x}}_{i})$ in (3), can be reformulated as a set of sequential Bernoulli choices with natural parameters $\eta_{h}({\bf x}_{i})=\mbox{logit}\{\nu_{h}({\bf{x}}_{i})\}=\log[\nu_{h}({\bf x}_{i})/\{1-\nu_{h}({\bf x}_{i})\}]$ under an exponential family representation. Hence, we can write

[TABLE]

allowing each $\eta_{h}({\bf x}_{i})$ to be explicitly interpreted as the log-odds of the probability of being allocated to component $h$, conditionally on the event of surviving to the first $1,\ldots,h-1$ components. This result might be helpful in driving prior specification for the stick-breaking weights, while allowing recent computational advances in Bayesian logistic regression [Polson et al., 2013] to be inherited in our density regression problem.

To conclude our Bayesian representation, we require priors for the log-odds $\eta_{h}({\bf{x}}_{i})$, $h\in\mathbb{N}$ in the continuation-ratio logistic regressions. A natural choice, which is consistent with classical generalized linear models [e.g. Nelder and Wedderburn, 1972], is to define $\eta_{h}({\bf{x}}_{i})$ as a linear combination of selected functions of the covariates $\bm{\psi}({\bf{x}}_{i})=\{\psi_{1}({\bf{x}}_{i}),\ldots,\psi_{R}({\bf{x}}_{i})\}^{\intercal}$ and consider Gaussian priors for the coefficients, thus obtaining

[TABLE]

Although the linearity assumption in (7) may seem restrictive, note that flexible formulations for $\eta_{h}({\bf{x}}_{i})$, including regression via splines and Gaussian processes, induce linear relations in the coefficients. Moreover, as we will outline in Section 3, the linearity assumption simplifies computations, while inducing a logistic-normal prior for each $\nu_{h}({\bf{x}}_{i})$. Although such a prior can closely approximate Dirichlet distributions [Aitchison and Shen, 1980], the logit stick-breaking does not induce beta distributed stick-breaking weights, and therefore it cannot be included in the class discussed by Ishwaran and James [2001]. However, one can easily adapt the theoretical results in Rodriguez and Dunson [2011] to our logit link. For example, the infinite summation of the mixing weights is such that $\sum_{h=1}^{\infty}\pi_{h}({\bf x}_{i})=1$ almost surely for any ${\bf x}\in\mathcal{X}$; see the Appendix for details. Moreover, the lsbp is highly similar in its probabilistic nature and properties to other popular predictor-dependent stick-breaking constructions. In particular, psbp can be approximated by lsbp, and viceversa, up to a simple transformation of the prior for each $\bm{\alpha}_{h}$. This is a natural consequence of the well known relationship between the probit and the logit function [Amemiya, 1981], since the mapping $\{1+\exp(-\bm{\psi}({\bf{x}})^{\intercal}{\bm{\alpha}}_{h})\}^{-1}$ can be roughly approximated by $\Phi\{\bm{\psi}({\bf{x}})^{\intercal}{\bm{\alpha}}_{h}\sqrt{\pi/8}\}$. This is summarized in Remark 1.

Remark 1**.**

The logit stick-breaking prior in (6)–(7), can be approximated by a probit stick-breaking process $\nu_{h}({\bf{x}})\approx\Phi\{\bm{\psi}({\bf{x}})^{\intercal}\bar{\bm{\alpha}}_{h}\}$, with $\bar{\bm{\alpha}}_{h}={\bm{\alpha}}_{h}\sqrt{\pi/8}\sim\mbox{N}_{R}\{\sqrt{\pi/8}\bm{\mu}_{\alpha},(\pi/8)\bm{\Sigma}_{\alpha}\}$, for every ${\bf{x}}\in\mathcal{X}$ and $h\in\mathbb{N}$.

Hence, a researcher considering a psbp could perform approximate inference leveraging our algorithms, after a suitable rescaling of the prior for each $\bm{\alpha}_{h}$. Moreover, this link suggests that the $O(\log n)$ growth of the number of clusters found in empirical studies on the psbp, should hold also for lsbp.

3 Bayesian computational methods

Although the lsbp and the associated computational procedures apply to a wider set of dependent mixture models and kernels, we focus, for the sake of clarity, on the general class of predictor-dependent infinite mixtures of Gaussians

[TABLE]

where $\tau_{h}=\sigma^{-2}_{h}$ is the precision parameter, whereas $\bm{\beta}_{h}=(\beta_{1h},\ldots,\beta_{Mh})^{\intercal}$ denotes a vector of coefficients linearly related to selected functions of the observed predictors $\bm{\lambda}({\bf x})=\{\lambda_{1}({\bf x}),\ldots,\lambda_{M}({\bf x})\}^{\intercal}$. Formulation (8) provides a flexible construction [Barrientos et al., 2012, Pati et al., 2013], and is arguably the most widely used in Bayesian density regression. As mentioned in Section 1, we provide here a detailed derivation of three computational methods for Bayesian density regression under model (8), with logit stick-breaking prior (6)–(7) for the mixing weights. In particular, we consider a Gibbs sampler converging to the exact posterior, an expectation-maximization (em) algorithm for point estimation, and a mean-field variational Bayes (vb) approximation for scalable posterior inference. The algorithms associated with these methods are available at https://github.com/tommasorigon/LSBP, along with the code to reproduce the application in Section 4.

In the classical predictor-independent mixture of Gaussians framework, these computational methods are closely related, and relevant connections can be drawn also with $k$-means and Bayesian $k$-means algorithms [Bishop, 2006, Kurihara and Welling, 2009]. A summary of these relations is depicted in Figure 1 of Kurihara and Welling [2009]. Broadly speaking, these strategies differ in how they handle unknown parameters and the involved latent quantities, either through maximization or by taking expectations. These connections are paralleled in the lsbp model, although our focus is mainly on Gibbs sampling, em and vb.

Before providing a detailed derivation of these different algorithms, we first study a truncated version of the random probability measure $p_{\bf{x}}$, which will be employed as an approximation of the infinite process. Indeed, although Gibbs samplers for infinite representations are available [Kalli et al., 2011], developing em and vb algorithms is not straightforward. In line with Rodriguez and Dunson [2011] and Ren et al. [2011], we develop detailed routines based on a finite representation. In particular, we model the first $H-1$ weights $\nu_{1}({\bf{x}}),\ldots,\nu_{H-1}({\bf{x}})$ and let $\nu_{H}({\bf{x}})=1$ for any ${\bf{x}}\in\mathcal{X}$, so that $\sum_{h=1}^{H}\pi_{h}({\bf{x}})=1$. Based on Theorem 1 below, this choice provides an accurate approximation of the infinite representation for sufficiently large truncations $H$.

Theorem 1**.**

For a sample ${\bf{y}}=(y_{1},\dots,y_{n})^{\intercal}$ with covariates ${\bf{X}}=({\bf{x}}_{1},\ldots,{\bf{x}}_{n})^{\intercal}$, let

[TABLE]

be the marginal joint density arising from a truncated lsbp prior with $H$ components, and define with $f_{\bf{X}}^{(\infty)}({\bf{y}})$ the same quantity in the infinite case. Moreover, let $\mu_{\nu}({\bf x})=\mbox{E}\{\nu_{h}({\bf x})\}$ be the expected value of a generic stick-breaking weight. Then $||f_{\bf{X}}^{(H)}({\bf{y}})-f_{\bf{X}}^{(\infty)}({\bf{y}})||_{1}\leq 4\sum_{i=1}^{n}\{1-\mu_{\nu}({\bf x}_{i})\}^{H-1},$ where $||\cdot||_{1}$ denotes the $L^{1}$–norm. Note that in the above formula the expectation is taken with respect to the lsbp prior law.

According to Theorem 1, for fixed sample size $n$ and covariates ${\bf{X}}$, the $L^{1}$ distance between $f_{\bf{X}}^{(H)}({\bf{y}})$ and $f_{\bf{X}}^{(\infty)}({\bf{y}})$ vanishes as $H\rightarrow\infty$, implying that the marginal density $f_{\bf{X}}^{(H)}({\bf{y}})$ converges to $f_{\bf{X}}^{(\infty)}({\bf{y}})$. This rate of decay is exponential in $H$, and therefore the number of components does not have to be very large in practice to accurately approximate the infinite representation, thus motivating computational methods based on truncated versions.

3.1 MCMC via Gibbs sampling

In deriving a Gibbs sampler for model (8) we focus on a dependent mixture of Gaussians with fixed $H$, and exploit the hierarchical representation (3) along with the continuation-ratio characterization of the logit stick-breaking prior, given in Section 2. Under these constructions, the joint law for the augmented model (3) and its parameters becomes

[TABLE]

with $\nu_{h}({\bf x}_{i})=\exp\{\psi({\bf x}_{i})^{\intercal}\bm{\alpha}_{h}\}/[1+\exp\{\psi({\bf x}_{i})^{\intercal}\bm{\alpha}_{h}\}]$, whereas $f(\bm{\alpha})f(\bm{\beta})f(\bm{\tau})=\prod_{h=1}^{H-1}f(\bm{\alpha}_{h})\prod_{h=1}^{H}f(\bm{\beta}_{h})f({\tau}_{h})$ denote the prior laws of the parameters comprising $\bm{\alpha}$, $\bm{\beta}$ and $\bm{\tau}$. As is clear from (9), given ${\bf{G}}=(G_{1},\ldots,G_{n})$, sampling of $\bm{\beta}_{h}$ and $\tau_{h}$, for $h=1,\ldots,H$, requires standard methods for Gaussian linear regression within each mixture component, as long as conditionally conjugate priors $\bm{\beta}_{h}\sim\mbox{N}_{M}(\bm{\mu}_{\beta},\bm{\Sigma}_{\beta})$ and $\tau_{h}\sim\mbox{Ga}(a_{\tau},b_{\tau})$, or normal-gammas for the pair $(\bm{\beta}_{h},\tau_{h})$, are employed. Here we focus on the first choice to keep notation more compact.

The updating of the $\bm{\alpha}_{h}$ parameters, for $h=1,\ldots,H-1$, relies instead on a set of separate Bayesian logistic regressions with responses $z_{ih}=1$ when $G_{i}=h$ and $z_{ih}=0$ if $G_{i}>h$, for those units $i$ having $G_{i}>h-1$, thus allowing parallel sampling from the full-conditional of each $\bm{\alpha}_{h}$. Adapting results from the recent Pólya-gamma data augmentation scheme [Polson et al., 2013] to our statistical model, these updatings can be easily accomplished by noticing that $\nu_{h}({\bf x}_{i})^{z_{ih}}\{1-\nu_{h}({\bf x}_{i})\}^{1-z_{ih}}=\int f_{{\bf{x}}_{i}}(z_{ih})f_{{\bf{x}}_{i}}(\omega_{ih})\mbox{d}\omega_{ih},$ with laws $f_{{\bf{x}}_{i}}(z_{ih})$ and $f_{{\bf{x}}_{i}}(\omega_{ih})$ defined as

[TABLE]

for every $i:G_{i}>h-1$ and $h=1,\ldots,H-1$. In (10), $f_{{\bf{x}}_{i}}(\omega_{ih})$ and $f(\omega_{ih})$ are the density functions of the Pólya-gamma random variables $\mbox{PG}\{1,\bm{\psi}({\bf{x}}_{i})^{\intercal}{\bm{\alpha}}_{h}\}$, and $\mbox{PG}(1,0)$, respectively. Hence, based on (10), the contribution to the augmented likelihood for each pair $(z_{ih},\omega_{ih})$ is proportional to a Gaussian kernel for transformed data $(z_{ih}-0.5)/\omega_{ih}$, provided that $f_{{\bf{x}}_{i}}(z_{ih})f_{{\bf{x}}_{i}}(\omega_{ih})$ $\propto\exp[(z_{ih}-0.5)\bm{\psi}({\bf{x}}_{i})^{\intercal}{\bm{\alpha}}_{h}-0.5\{\bm{\psi}({\bf{x}}_{i})^{\intercal}{\bm{\alpha}}_{h}\}^{2}\omega_{ih}]$. This allows conditionally conjugate updating steps for each $\bm{\alpha}_{h}$ under a classical Bayesian linear regression framework. Refer to Choi and Hobert [2013], Wang and Roy [2018a, b] for further theoretical properties of the Pólya-gamma scheme. Finally, note that in (10), the latent indicators $z_{ih}$ and the Pólya-gamma random variables $\omega_{ih}$ are conditionally independent given the coefficients $\bm{\alpha}_{h}$ for $i:G_{i}>h-1$. This is in contrast with the data augmentation underlying the psbp, which would lead to more complex calculations, especially in the em and vb algorithms discussed in Sections 3.2 and 3.3.

The detailed steps of the Gibbs sampler for the truncated representation of model (8) are outlined in Algorithm 1. In this routine, $\bm{\Lambda}_{h}({\bf{x}})$ and $\bm{\Psi}_{h}({\bf{x}})$ denote the $n_{h}\times M$ and the $\bar{n}_{h}\times R$ predictor matrices in (8) and (7) having row entries $\bm{\lambda}({\bf x}_{i})^{\intercal}$ and $\bm{\psi}({\bf x}_{i})^{\intercal}$, for only those statistical units $i$ such that $G_{i}=h$ and $G_{i}>h-1$, respectively. We shall also emphasize that step [1] can be run in parallel across units $i=1,\ldots,n$, whereas parallel computing for the different mixture components can be easily implemented in steps [2], [3] and [4].

3.2 EM algorithm

In high-dimensional studies, the Gibbs sampler described in Section 3.1 could face computational bottlenecks. If a point estimate of model (8) is the main quantity of interest, for example for prediction purposes, one possibility is to rely on a more efficient procedure specifically designed for this goal, such as the em [Dempster et al., 1977]. The implementation of a simple em for a finite representation of model (8) under the lsbp prior benefits from the Pólya-gamma data augmentation, which has analytical expectation and allows direct maximization within a Gaussian linear regression framework. Note that, although the em algorithm is commonly implemented for maximum likelihood estimation, it can be easily modified to estimate posterior modes [e.g. Dempster et al., 1977].

The proposed em in Algorithm 2 alternates between a maximization step for the parameters $(\bm{\alpha},\bm{\beta},\bm{\tau})$ and an expectation step for the augmented data $(\bm{\zeta}_{i},\bar{\bm{\omega}}_{i})$, $i=1,\ldots,n$, with $\bm{\zeta}_{i}=\{\zeta_{i1}=\mathbbm{1}(G_{i}=1),\ldots,\zeta_{iH}=\mathbbm{1}(G_{i}=H)\}^{\intercal}$ the vector of binary indicators denoting the membership to a mixture component, and $\bar{\bm{\omega}}_{i}=(\bar{\omega}_{i1},\ldots,\bar{\omega}_{iH-1})^{\intercal}$ the corresponding Pólya-gamma augmented data. Although this data augmentation parallels the one described for the Gibbs sampler, we adopt a slightly different notation for the Pólya-gamma random variables $\bar{\omega}_{ih}$, to emphasize that we are considering $n$ units, and not only those for which the cluster indicators $G_{i}>h-1$. Indeed, in line with the em rationale, we do not condition on the membership indicators and on the Pólya-gamma latent random variables, but we rather take expectations with respect to their conditional distributions. For the same reason, in this case we work directly with the component indicator variables $\bm{\zeta}_{i}$ instead of the binary vectors ${\bf{z}}_{i}=(z_{i1},\ldots,z_{iH-1})^{\intercal}$ in (5).

Based on the data augmentations outlined in (3) and (10), the complete log-posterior $\mbox{log}f_{{\bf{x}}}(\bm{\alpha},\bm{\beta},\bm{\tau}\mid{\bf{y}},\bm{\zeta},\bar{\bm{\omega}})$ underlying the proposed em routine, can be written as

[TABLE]

where $\ell_{{\bf{x}}_{i}}(\bm{\alpha},\bm{\beta},\bm{\tau};{y}_{i},\bm{\zeta}_{i},\bar{\bm{\omega}}_{i})$ is the contribution of unit $i$ to the complete log-likelihood. Working on the complete log-likelihood has relevant benefits. Indeed, exploiting equations (3) and (5), and the results in Polson et al. [2013] summarized in (10), the term $\ell_{{\bf{x}}_{i}}(\bm{\alpha},\bm{\beta},\bm{\tau};{y}_{i},\bm{\zeta}_{i},\bar{\bm{\omega}}_{i})=\ell_{{\bf{x}}_{i}}(\bm{\beta},\bm{\tau};{y}_{i},\bm{\zeta}_{i})+\ell_{{\bf{x}}_{i}}{(\bm{\alpha};\bm{\zeta}_{i},\bar{\bm{\omega}}_{i})}$, can be factorized as

[TABLE]

where $\bar{\kappa}_{ih}=\zeta_{ih}-0.5\sum_{l=h}^{H}\zeta_{il}$. Hence, both terms in equation (12) are linear in the augmented data $(\bm{\zeta}_{i},\bar{\bm{\omega}}_{i})$, and represent the sum of Gaussian kernels. This linearity property simplifies computations in the expectation step for the complete log-posterior in equation (11), whereas the Gaussian structure allows simple maximizations. Since the joint maximization of the expected complete log-posterior with respect to $(\bm{\beta},\bm{\tau})$ is intractable, we rely on a conditional maximization procedure [Meng and Rubin, 1993] in the last step of Algorithm 2, which provides analytical solutions.

3.3 Mean-field variational Bayes

Section 3.2 provides a scalable procedure for estimation of posterior modes in large-scale problems. However, an appealing aspect of the Bayesian approach is in allowing uncertainty quantification via inference on the entire posterior. The Gibbs sampler in Section 3.1 represents an appealing procedure which converges to the exact posterior, but faces computational bottlenecks. This motivates scalable variational methods for approximate Bayesian inference [Bishop, 2006, Blei et al., 2017]. Clearly, these computational gains do not come without some drawbacks. For example, variational approximations typically underestimate posterior variability. This issue might be mitigated via a post-processing operation as in Giordano et al. [2015], at the cost of an additional computational step.

Due to the Pólya-gamma data augmentation, our variational strategy is framed within the well-established exponential family setting, for which there exists a closed-form coordinate ascent variational inference algorithm (cavi). Compared to more accurate black-box variational strategies [e.g. Ranganath et al., 2014], the cavi algorithm is appealing because it requires no tuning. Moreover, recent theoretical properties for this class of computational methods [Blei et al., 2017] are inherited by our variational algorithm. This seems in contrast with the variational strategy discussed by Ren et al. [2011], which considers a local approximation based on the lower bound of Jaakkola and Jordan [2000]. However, the recent contribution of Durante and Rigon [2019] allows to draw a sharp connection between the Pólya-gamma data augmentation and the Jaakkola and Jordan [2000] lower bound. As a consequence, the vb approach we propose relies on the same optimization problem considered by Ren et al. [2011].

Compared to the Gibbs sampler in Section 3.1, here we augment the entire model (8) with respect to the binary vectors ${\bf z}_{i}=(z_{i1},\dots,z_{iH-1})^{\intercal}$, $i=1,\dots,n$ comprising ${\bf z}$, rather than using the membership indicators ${\bf G}$. Hence, the joint law $f_{\bf{x}}({\bf{y}},\bm{\alpha},\bm{\beta},\bm{\tau},{\bf{z}},\bm{\omega})=f_{\bf{x}}({\bf{y}}\mid{\bf{z}},\bm{\beta},\bm{\tau})f_{\bf{x}}({\bf{z}}\mid\bm{\alpha})f_{\bf{x}}(\bm{\omega}\mid\bm{\alpha})f(\bm{\alpha})f(\bm{\beta})f(\bm{\tau})$ is equal to

[TABLE]

where $z_{iH}=1$. Our goal is to find a variational distribution $q_{{\bf{x}}}(\bm{\alpha},\bm{\beta},\bm{\tau},{\bf{z}},\bm{\omega})$ that best approximates the joint posterior $f_{{\bf{x}}}(\bm{\alpha},\bm{\beta},\bm{\tau},{\bf z},\bm{\omega}\mid{\bf y})$, while maintaining simple computations. This can be obtained by minimizing the Kullback-Leibler divergence $\textsc{kl}\{q_{{\bf{x}}}(\bm{\alpha},\bm{\beta},\bm{\tau},{\bf{z}},\bm{\omega})\mid\mid f_{{\bf{x}}}(\bm{\alpha},\bm{\beta},\bm{\tau},{\bf{z}},\bm{\omega}\mid{\bf{y}})\}$ between the variational distribution and the full posterior, or, alternatively, by maximizing the evidence lower bound $\textsc{elbo}\{q_{{\bf{x}}}(\bm{\alpha},\bm{\beta},\bm{\tau},{\bf{z}},\bm{\omega})\}$ of the log-marginal density $\log f_{{\bf{X}}}^{(H)}({\bf{y}})$, since $\log f_{{\bf{X}}}^{(H)}({\bf{y}})$ can be analytically expressed as the sum of the elbo and the positive kl divergence. This evidence lower bound to be maximized can be expressed as

[TABLE]

Without further restrictions, the Kullback-Leibler divergence is minimized when the variational distribution is equal to the true posterior, which is intractable. To address this issue, a common strategy is to assume that the variational distribution belongs to a mean-field family [see e.g. Blei et al., 2017]. This incorporates a posteriori independence among distinct groups of parameters, implying that the variational distribution can be expressed as the product of marginal laws. Specifically, we consider the following factorization for the variational distribution

[TABLE]

Note that we are not making specific assumptions about the functional form of the variational distributions. Combining (13) with (14), we obtain a tractable expression for the elbo, which can be easily maximized as in Bishop [2006, Ch. 10]. In particular, the optimal solutions are provided by the following system of equations

[TABLE]

where ${\bf{z}}_{i,-h}$ denotes the vector of binary indicators ${\bf z}_{i}$ without considering the $h$th one, whereas $c_{\bm{\beta}_{h}}$, $c_{{\tau}_{h}}$, $c_{\bm{\alpha}_{h}}$, $c_{z_{ih}}$ and $c_{\omega_{ih}}$, are additive constants with respect to the argument in the corresponding variational distribution. Each expectation in the above equations is evaluated with respect to the variational distribution of the other parameters, and therefore we need to rely on iterative methods to find the optimal solution. We consider the coordinate ascent variational inference (cavi) iterative procedure—described in Algorithm 3—which maximizes the variational distribution of each parameter based on the current estimate for the remaining ones [e.g. Bishop, 2006, Ch. 10]. This procedure generates a monotone sequence for the $\textsc{elbo}\{q_{{\bf{x}}}(\bm{\alpha},\bm{\beta},\bm{\tau},{\bf{z}},\bm{\omega})\}$, which ensures convergence to a local joint maximum. Refer to Blei et al. [2017] for practical guidelines to address issues of local maxima via multiple runs. Finally, note that as shown in Algorithm 3, the normalizing constants in the above equations have not to be computed numerically, since kernels of well known distributions can be recognized.

4 Epidemiology application

We compare the performance of the three computational methods developed in Section 3, in a toxicology study. Consistent with recent interests in Bayesian density regression [e.g. Dunson and Park, 2008, Hwang and Pennell, 2014, Canale et al., 2018], we focus on a dataset aimed at studying the relationship between the DDE concentration in maternal serum, and the gestational days at delivery [Longnecker et al., 2001].

The DDE is a metabolite of ddt, which is still used against malaria-transmitting mosquitoes in certain developing countries—according to the Malaria Report 2015 from the World Health Organization—thus raising concerns about its adverse effects on premature delivery. Popular studies in reproductive epidemiology address this goal by dichotomizing the gestational age at delivery (GAD) with a clinical threshold, so that births occurred before the $37$th week are considered preterm. Although this approach allows for a simpler modeling strategy, it leads to a clear loss of information. In particular, a greater risk of mortality and morbidity is associated with preterm birth, which increases rapidly as the GAD decreases. This has motivated an increasing interest in modeling how the entire distribution of GAD changes with DDE exposure [e.g. Dunson and Park, 2008, Hwang and Pennell, 2014, Canale et al., 2018].

Data are composed by $n=2312$ measurements $(x_{i},y_{i})$, $i=1,\dots,n$, where $x_{i}$ denotes the DDE concentration, and $y_{i}$ is the gestational age at delivery for woman $i$. Our goal is to reproduce the analyses in Dunson and Park [2008] on this dataset, and compare the inference and computational performance of the mcmc via Gibbs sampling, the em algorithm, and the vb routine proposed in Section 3. Note that, consistent with the main novelty of this contribution, we do not attempt to improve the flexibility and the efficiency of the available statistical models for Bayesian density regression—such as the kernel stick-breaking [Dunson and Park, 2008], and the psbp [Rodriguez and Dunson, 2011]. Indeed, as discussed in Sections 1 and 2, these representations are expected to provide a comparable performance to our lsbp in terms of inference. However, unlike current models for Bayesian density regression, inference under the lsbp is available under a broader variety of simple computational methods, thus facilitating implementation of the same model in a wider range of applications—including large $M$, $R$ and $n$ settings. Due to this, the main focus is on providing an empirical comparison of the algorithms in Section 3, while using results in Dunson and Park [2008] as a benchmark to provide reassurance that inference under the lsbp is comparable to alternative representations.

We apply the predictor-dependent mixture of Gaussians (8) with lsbp (6)–(7), to a normalized version of the DDE and GAD $(\bar{x}_{i},\bar{y}_{i})$, $i=1,\ldots,n$, and then show results for $f_{x}(y)$ on the original scale of the data. Consistent with previous works [Dunson and Park, 2008, Canale et al., 2018], we let $M=2$, with $\lambda_{1}(\bar{x}_{i})=1$ and $\lambda_{2}(\bar{x}_{i})=\bar{x}_{i}$, for every $i=1,\ldots,n$, and rely instead on a flexible representation for $\eta_{h}(\bar{x}_{i})$ to characterize changes in the stick-breaking weights with DDE. In particular, each $\eta_{h}(\bar{x}_{i})$ is defined via a natural cubic spline basis $\bm{\psi}(\bar{x}_{i})=\{1,\psi_{1}(\bar{x}_{i}),\dots,\psi_{5}(\bar{x}_{i})\}^{\intercal}$, for every $h=1,\ldots,H-1$. Bayesian posterior inference—under the three computational methods developed in Section 3—is instead performed with default hyperparameters $\bm{\mu}_{\beta}=(0,0)^{\intercal}$, $\bm{\Sigma}_{\beta}=\mbox{I}_{2\times 2}$, $\bm{\mu}_{\alpha}=(0,\dots,0)^{\intercal}$, $\bm{\Sigma}_{\alpha}=\mbox{I}_{6\times 6}$ and $a_{\tau}=b_{\tau}=1$. For the total number of mixture components we consider $H=20$, and allow the shrinkage induced by the stick-breaking prior to adaptively delete redundant components not required to characterize the data. As shown in Figure 2, these choices allows accurate inference on $f_{x}(y)$.

In providing posterior inference under the Gibbs sampling algorithm described in Section 3.1, we rely on $30{,}000$ iterations, after discarding the first $5{,}000$ as a burn-in, and initialize the routine from random starting values sampled from the prior. Analysis of the traceplots for the quantities discussed in Figures 2 and 3 showed that this choice is sufficient for good convergence. The em algorithm and the vb procedures discussed in Sections 3.2 and 3.3, respectively, are instead run until convergence to a modal solution. Since such modes could be local, we run both algorithms for different initial values, and consider the solutions having the highest log-posterior and elbo, respectively. We also controlled the monotonicity of the sequences for these quantities, in order to further validate the correctness of our derivations. In this study, the em and the vb reach convergence in about $2$ and $6$ seconds, respectively, whereas the Gibbs sampler requires $5$ minutes, using a MacBook Air with a Intel Core i5.

Similarly to Figure 3 in Dunson and Park [2008], Figure 2 provides posterior inference for the conditional density $f_{x}(y)$ evaluated at the $0.1$, $0.6$, $0.9$, $0.99$ quantiles of DDE, for the three algorithms. Histograms for the GAD, are instead obtained by grouping the response data according to a binning of the DDE with cut-offs at the central values of subsequent quantiles, so that the conditional density can be plotted alongside the corresponding histogram. Results in Figure 2 confirm accurate fit to the data and suggest that the left tail of the GAD distribution—associated with preterm deliveries—increasingly inflates as DDE grows. Moreover, as seen in Figure 2, the three algorithms have similar results, thus providing empirical reassurance for the goodness of the proposed routines. Although the em outputs a posterior mode, such an estimate matches closely the posterior mean of the Gibbs sampler, whereas the vb tends to over-smooth some modes of the conditional distribution. This is likely due to the fact that the vb provides an approximation of the posterior, instead of the exact one. However, unlike for the em, this routine allows uncertainty quantification, and provides a much scalable methodology compared to the Gibbs sampler, thus representing a valid candidate in high-dimensional inference when the focus is on specific functionals of $f_{x}(y)$. Indeed, as shown in Figure 3, when the aim is to exploit $f_{x}(y)$ to infer conditional preterm probabilities $\mbox{pr}(y<y^{*}\mid x)=\int_{-\infty}^{y^{*}}f_{x}(y)\mbox{d}y$ with $y^{*}\in(7\times 33,7\times 35,7\times 37,7\times 40)$ denoting a clinical threshold, the vb provides very similar conclusions.

Prior to conclude our analysis, note that the results in Figures 2 and 3 are similar to those obtained under the kernel stick-breaking prior in Dunson and Park [2008]. This provides empirical guarantee that the flexibility characterizing popular Bayesian nonparametric models for density regression is maintained also under lsbp, which has the additional relevant benefit of facilitating computational implementation of these methodologies. Minor differences are found at extreme DDE exposures, but this is mainly due to the sparsity of the data in this subset of the predictor space.

5 Discussion

The focus of this paper has been on providing novel methods to facilitate the computational implementation of Bayesian nonparametric models for density regression in a broad range of applications. To address this goal, we have proposed an alternative reparameterization of the predictor-dependent stick-breaking weights, which relies on a set of sequential logistic regressions. This representation has relevant connections with continuation-ratio logistic regressions and Pólya-gamma data augmentation, thus allowing simple derivation of several algorithms of routine use in Bayesian inference. The proposed computational methods are empirically evaluated in a toxicology study, obtaining good results and reassurance that the lsbp maintains the same flexibility and efficiency properties characterizing popular Bayesian nonparametric models for density regression.

Although our dependent mixture of Gaussians provides a flexible representation, it is worth considering extensions to other kernels. For example, all our algorithms can be easily adapted to predictor-independent kernels coming from an exponential family, when conjugate priors for their parameters are used. Similar derivations are also possible for predictor-dependent kernels within a generalized linear model representation, provided that conjugate priors for the coefficients can be found [e.g. Chen and Ibrahim, 2003].

Acknowledgments

The authors are grateful to the Associate Editor and the four referees for the insightful comments and suggestions that led to a substantial improvement of the paper. This work was partially supported by grant 1R01ES027498 of the National Institutes of Environmental Health Sciences of the United States National Institutes of Health, and by the MIUR–PRIN 2017 grant 20177BRJXS.

Appendix. Proofs and additional properties of the lsbp

Proposition 1**.**

For any fixed ${\bf{x}}\in\mathcal{X}$, $\sum_{h=1}^{\infty}\pi_{h}({\bf{x}})=1$ almost surely, with $\pi_{h}({\bf{x}})$ factorized as in (6) and $\bm{\alpha}_{h}\sim\mbox{N}_{R}(\bm{\mu}_{\alpha},\bm{\Sigma}_{\alpha})$ independently for every $h\in\mathbb{N}$. Hence, the lsbp provides a well defined predictor-dependent random probability measure $p_{\bf x}$ at every ${\bf x}\in\mathcal{X}$.

Proof of Proposition 1. Recalling results in Ishwaran and James [2001], we have that $\sum_{h=1}^{\infty}\pi_{h}({\bf{x}})=1$ almost surely if and only if the equality $\sum_{h=1}^{\infty}\mbox{E}[\log\{1-\nu_{h}({\bf{x}})\}]=-\infty$ holds. Since $\log\{1-\nu_{h}({\bf{x}})\}$ is concave in $\nu_{h}({\bf{x}})$ for every ${\bf{x}}\in\mathcal{X}$ and $h\in\mathbb{N}$, by the Jensen inequality $\mbox{E}[\log\{1-\nu_{h}({\bf{x}})\}]\leq\log[1-\mbox{E}\{\nu_{h}({\bf{x}})\}]$. Therefore, since $\nu_{h}({\bf{x}})\in(0,1)$, we have that $0<\mbox{E}\{\nu_{h}({\bf{x}})\}=\mu_{\nu}({\bf{x}})<1$, thereby providing $\log\{1-\mu_{\nu}({\bf{x}})\}<0$. Leveraging these results, the proof of Proposition 1 follows after noticing that $\sum_{h=1}^{\infty}\mbox{E}[\log\{1-\nu_{h}({\bf{x}})\}]\leq\sum_{h=1}^{\infty}\log\{1-\mu_{\nu}({\bf{x}})\}=-\infty$.

Proof of Theorem 1. Adapting the proof of Theorem 1 in Ishwaran and James [2002] to our representation we have

[TABLE]

Since $\sum_{h=1}^{H-1}\pi_{h}({\bf x}_{i})\leq 1$, and $1=\prod_{i=1}^{n}1$, we can write $1-\prod_{i=1}^{n}\sum_{h=1}^{H-1}\pi_{h}({\bf x}_{i})=\prod_{i=1}^{n}1-\prod_{i=1}^{n}\sum_{h=1}^{H-1}\pi_{h}({\bf x}_{i})\leq\sum_{i=1}^{n}\{1-\sum_{h=1}^{H-1}\pi_{h}({\bf x}_{i})\}$ [Billingsley, 1995, pp. 358]. Hence $||f_{\bf{X}}^{(H)}({\bf{y}})-f_{\bf{X}}^{(\infty)}({\bf{y}})||_{1}\leq 4[n-\sum_{i=1}^{n}\sum_{h=1}^{H-1}\mbox{E}\{\pi_{h}({\bf x}_{i})\}],$ with $\sum_{h=1}^{H-1}\mbox{E}\{\pi_{h}({\bf x}_{i})\}=\sum_{h=1}^{H-1}\mu_{\nu}({\bf x}_{i})\{1-\mu_{\nu}({\bf x}_{i})\}^{h-1}$ $=1-\{1-\mu_{\nu}({\bf x}_{i})\}^{H-1}$. Substituting this quantity in $4[n-\sum_{i=1}^{n}\sum_{h=1}^{H-1}\mbox{E}\{\pi_{h}({\bf x}_{i})\}]$, we obtain the final bound $4\sum_{i=1}^{n}\{1-\mu_{\nu}({\bf x}_{i})\}^{H-1}$. ∎