Bayesian Variable Selection for Gaussian copula regression models

TL;DR

This paper introduces a Bayesian variable selection method for Gaussian copula regression models that handle multiple response types and complex dependencies, improving model estimation efficiency and applicability.

Contribution

The paper presents a novel Bayesian approach with a specialized MCMC algorithm for variable selection in multivariate regression with diverse response types.

Findings

Effective in simulated data scenarios

Successfully applied to real datasets with mixed response types

Enhanced exploration of predictor model space

Abstract

We develop a novel Bayesian method to select important predictors in regression models with multiple responses of diverse types. A sparse Gaussian copula regression model is used to account for the multivariate dependencies between any combination of discrete and/or continuous responses and their association with a set of predictors. We utilize the parameter expansion for data augmentation strategy to construct a Markov chain Monte Carlo algorithm for the estimation of the parameters and the latent variables of the model. Based on a centered parametrization of the Gaussian latent variables, we design a fixed-dimensional proposal distribution to update jointly the latent binary vectors of important predictors and the corresponding non-zero regression coefficients. For Gaussian responses and for outcomes that can be modeled as a dependent version of a Gaussian response, this proposal leads to a Metropolis-Hastings step that allows an efficient exploration of the predictors' model space. The proposed strategy is tested on simulated data and applied to real data sets in which the responses consist of low-intensity counts, binary, ordinal and continuous variables.

Equations

$$F_{Y_1,\dots ,Y_m}(y_1,\dots ,y_m)=C\{F_{Y_1}(y_1),\dots ,F_{Y_m}(y_m)\}.$$

$$C(u_1,\dots ,u_m;\mathbf{R})={\Phi }_m\{{\Phi }^{-1}(u_1),\dots ,{\Phi }^{-1}(u_m);\mathbf{R}\},$$

$$c(u_1,\dots ,u_m;\mathbf{R})=|\mathbf{R}{|}^{-\frac{1}{2}}\exp \left\{-\frac{1}{2}({\tilde{\mathbf{z}}}^T{\mathbf{R}}^{-1}\tilde{\mathbf{z}}-{\tilde{\mathbf{z}}}^T\tilde{\mathbf{z}})\right\},$$

$$y_{ik}=h_{ik}^{-1}({\tilde{z}}_{ik}),{\tilde{\mathbf{Z}}}_i\sim {\mathrm{N}}_m(0,\mathbf{R}),$$

$${\beta }_{kj}|{\gamma }_{kj}\stackrel{\text{iid}}{\sim }(1-{\gamma }_{kj}){\delta }_0+{\gamma }_{kj}\mathrm{N}(0,v),j=1,\dots ,p_k,$$

$${\gamma }_{kj}|{\pi }_k\stackrel{\text{iid}}{\sim }\mathrm{Ber}({\pi }_k),j=1,\dots ,p_k$$

$${\pi }_k\sim \mathrm{Beta}(a_k,b_k),$$

$$p({\mathbf{\gamma }}_k)=\frac{\mathrm{B}\left(a_k+|{\mathbf{\gamma }}_k|,b_k+p_k-|{\mathbf{\gamma }}_k|\right)}{\mathrm{B}(a_k,b_k)},$$

$$\mathrm{vec}(\tilde{\mathbf{W}})\sim {\mathrm{N}}_{mn}(0,\mathbf{\Sigma }\otimes {\mathbf{I}}_n),$$

$$g_{\ell }\sim \mathrm{Ber}({\pi }_{\mathbf{G}}),{\pi }_{\mathbf{G}}\sim \mathrm{Unif}(0,1),\ell =1,\dots ,m(m-1)/2,\mathbf{G}\in \mathcal{G}.$$

$$p(\mathbf{D},\mathbf{R},\mathbf{G})=p(\mathbf{D}|\mathbf{R})p(\mathbf{R}|\mathbf{G})p(\mathbf{G}),$$

$${\delta }_k^2|r^{kk}\stackrel{\text{iid}}{\sim }\mathrm{IGam}((m+1)/2,r^{kk}/2),k=1,\dots ,m$$

$$p({\mathbf{\beta }}_k,{\mathbf{\gamma }}_k|{\mathbf{y}}_k,{\mathbf{\theta }}_k,{\tilde{\mathbf{Z}}}_{-k},\mathbf{D},\mathbf{R})\propto p({\mathbf{y}}_k|{\mathbf{\beta }}_k,{\mathbf{\gamma }}_k,{\mathbf{\theta }}_k,{\tilde{\mathbf{Z}}}_{-k},\mathbf{D},\mathbf{R})p({\mathbf{\beta }}_k|{\mathbf{\gamma }}_k)p({\mathbf{\gamma }}_k),$$

$$p({\mathbf{y}}_k|{\mathbf{\beta }}_k,{\mathbf{\gamma }}_k,{\mathbf{\theta }}_k,{\tilde{\mathbf{Z}}}_{-k},\mathbf{D},\mathbf{R})$$

$$p({\mathbf{y}}_k|{\mathbf{\beta }}_k,{\mathbf{\gamma }}_k,{\mathbf{\theta }}_k,{\tilde{\mathbf{Z}}}_{-k},\mathbf{D},\mathbf{R})=\prod _{i=1}^n\left\{\Phi \left(\frac{u_{ik}-{\tilde{\mu }}_{ik|-k}}{{\tilde{\sigma }}_{k|-k}}\right)-\Phi \left(\frac{{\ell }_{ik}-{\tilde{\mu }}_{ik|-k}}{{\tilde{\sigma }}_{k|-k}}\right)\right\},$$

$$\{\begin{array}{ll}{Z}_{ik}={\mathbf{x}}_{i,{\mathbf{\gamma }}_k}^T{\mathbf{\beta }}_{{\mathbf{\gamma }}_k}+{\delta }_k{\tilde{Z}}_{ik}& \text{if}{\mathbf{y}}_k\text{is continuous,}\\ Z_{ik}={\mathbf{x}}_{i,{\mathbf{\gamma }}_k}^T{\mathbf{\beta }}_{{\mathbf{\gamma }}_k}+{\tilde{Z}}_{ik}& \text{if}{\mathbf{y}}_k\text{is discrete},\end{array}$$

$$q({\mathbf{\beta }}_k^{\ast }|{\mathbf{\gamma }}_k^{\ast },{\mathbf{z}}_k,{\tilde{\mathbf{Z}}}_{-k},\mathbf{D},\mathbf{R})$$

$$\propto {\mathcal{N}}_{|{\mathbf{\gamma }}_k^{\ast }|}({\mathbf{\beta }}_k^{\ast };{\mathbf{V}}_{{\mathbf{\gamma }}_k^{\ast }}{\mathbf{X}}_{{\mathbf{\gamma }}_k^{\ast }}^T{\sigma }_{k|-k}^{-2}\{{\mathbf{z}}_k-{\sigma }_k{\tilde{\mathbf{\mu }}}_{k|-k}\},{\mathbf{V}}_{{\mathbf{\gamma }}_k^{\ast }}),$$

$$q({\mathbf{\beta }}_k^{\ast }|{\mathbf{\gamma }}_k^{\ast },{\mathbf{z}}_k,\mathbf{D})\propto {\mathcal{N}}_{|{\mathbf{\gamma }}_k^{\ast }|}({\mathbf{\beta }}_k^{\ast };{\mathbf{V}}_{{\mathbf{\gamma }}_k^{\ast }}{\mathbf{X}}_{{\mathbf{\gamma }}_k^{\ast }}^T{\sigma }_k^{-2}{\mathbf{z}}_k,{\mathbf{V}}_{{\mathbf{\gamma }}_k^{\ast }}),$$

$$\alpha =1\wedge \frac{{\left|{\mathbf{V}}_{{\mathbf{\gamma }}_k^{\ast }}\right|}^{1/2}{v}^{|{\mathbf{\gamma }}_k|/2}\exp \left(\frac{1}{2}{\mathbf{m}}_{{\mathbf{\gamma }}_k^{\ast }}^T{\mathbf{V}}_{{\mathbf{\gamma }}_k^{\ast }}^{-1}{\mathbf{m}}_{{\mathbf{\gamma }}_k^{\ast }}\right)p({\mathbf{\gamma }}_k^{\ast })q({\mathbf{\gamma }}_k|{\mathbf{\gamma }}_k^{\ast })}{{\left|{\mathbf{V}}_{{\mathbf{\gamma }}_k}\right|}^{1/2}{v}^{|{\mathbf{\gamma }}_k^{\ast }|/2}\exp \left(\frac{1}{2}{\mathbf{m}}_{{\mathbf{\gamma }}_k}^T{\mathbf{V}}_{{\mathbf{\gamma }}_k}^{-1}{\mathbf{m}}_{{\mathbf{\gamma }}_k}\right)p({\mathbf{\gamma }}_k)q({\mathbf{\gamma }}_k^{\ast }|{\mathbf{\gamma }}_k)},$$

$$p({\tilde{z}}_{ik}|{\mathbf{y}}_k,{\mathbf{\beta }}_k,{\mathbf{\gamma }}_k,{\mathbf{\theta }}_k,{\tilde{\mathbf{Z}}}_{-k},\mathbf{D},\mathbf{R})\propto \mathcal{N}({\tilde{z}}_{ik};{\tilde{\mu }}_{ik|-k},{\tilde{\sigma }}_{k|-k}^2)\mathbb{I}\{{\tilde{z}}_{ik}\in ({\ell }_{ik},u_{ik}]\},$$

$$p(\mathbf{G}|\tilde{\mathbf{W}})=\frac{p(\tilde{\mathbf{W}}|\mathbf{G})p(\mathbf{G})}{{\sum }_{\mathbf{G}}p(\tilde{\mathbf{W}}|\mathbf{G})p(\mathbf{G})},$$

Full text

Keywords: Gaussian copula; Mixed data; Multiple-response regression models; Sparse covariance matrix; Variable selection.

1 Introduction

The identification of important predictors in linear and non-linear regression models is one of the most frequently studied questions in statistical theory. In Bayesian statistics, this problem is known as Bayesian variable selection (BVS) and for Gaussian responses, there is an extensive literature for an efficient detection of important predictors in both single- and multi-response regression models. In the case of a single response, a list of relevant papers includes, but is not limited to, george1997approaches, liang2008mixtures, guan2011bayesian and Rockova2018, while brown1998multivariate, holmes2002accounting deal with the same problem in the multi-response case. BVS in single-response non-linear regression models has also received great attention. In dellaportas2002bayesian and forster2012reversible the corresponding methods are reviewed and advances are proposed.

More recently, there has been an increasing interest in the joint analysis of outcomes of diverse types, for instance, continuous, binary, categorical and count data, given their availability from studies involving multivariate data, see for example hoff2007extending, murray2013bayesian, zhang2015bayesian and bhadra2018inferring. In regression analysis, the most popular model used to account for the multivariate dependencies between any combination of discrete and/or continuous responses is the Gaussian copula regression (GCR) model (song2009joint) in which each response is associated with a (potentially different) set of predictors. When only Gaussian responses are considered, this is known as Seemingly Unrelated Regression (SUR) model (zellner1962efficient). Recent contributions for sparse Bayesian SUR models include, for instance, wang2010sparse and Deshpande2019. Bayesian methods for the estimation of the regression coefficients of the GCR model with a fixed set of predictors have also been proposed (pitt2006efficient). Despite the growing Bayesian literature regarding an efficient selection of important predictors, to the best of our knowledge, variable selection for the GCR model has not been attempted. We propose here the first fully Bayesian approach for model selection in regression models with multiple responses of diverse types.

The main obstacle of the application of BVS in single-response non-linear models, as well as in the SUR model and the GCR model with responses of diverse types, is the non-tractability of the marginal likelihood. To overcome this problem, Markov chain Monte Carlo (MCMC) algorithms are based either on Laplace approximation, see for example bove2011hyper, or on a Metropolis-Hastings (M-H) step in which the dimension of the proposal distribution is not fixed at each iteration (forster2012reversible). The latter is an application of the Reversible Jump algorithm (green1995reversible) which is known to experience a low acceptance rate when the trans-dimension proposal distribution is not devised carefully, resulting in MCMC samplers with poor mixing (brooks2003efficient). In addition, in current applications, the number of predictors is often very large and any MCMC algorithm for BVS in both linear and non-linear regression has to be designed carefully in order to explore successfully the ultra-high dimensional model space consisting of all the possible subsets of predictors (bottolo2010evolutionary; Lamnisos2009).

The main contribution of this paper is the development of a Bayesian approach for the joint update of the latent binary vector of important predictors and the corresponding vector of non-zero regression coefficients for each response of the GCR model. By using the proposed strategy, we perform BVS without any approximation. We also avoid the Reversible Jump algorithm by utilizing the Gaussian latent variables of the GCR model to construct a proposal distribution defined on a fixed-dimensional space. For Gaussian responses and for the outcomes that can be modeled as a dependent version of a Gaussian response, e.g., the Probit model for binary data, the designed proposal distribution allows the “implicit marginalisation” of the non-zero regression coefficients in the M-H acceptance probability (holmes2006bayesian), favoring the exploration of the predictors model space in a more efficient manner.

Modeling the dependence amongst the responses is another key aspect in GCR. Until recently devising an efficient MCMC algorithm for a structured (constrained) covariance matrix that includes the identifiability conditions for the non-linear responses has been a difficult task. Here, we follow the solution proposed by talhouk2012efficient and specify a conjugate prior on the correlation matrix. We utilize the parameter expansion for data augmentation (liu1999parameter; van2001art) to expand the correlation into a covariance matrix. We also adopt the idea of covariance selection to obtain a parsimonious representation of the dependence amongst the responses. Based on the theory of decomposable Gaussian graphical models (lauritzen1996graphical), we use the hyper-inverse Wishart distribution as the prior density for the covariance matrix. This prior allows some of the off-diagonal elements of the inverse covariance matrix to be identical to zero and to estimate the conditional dependence pattern of the observations (webb2008bayesian).

We tested the performance of our model, Bayesian Variable Selection for Gaussian Copula Regression (BVSGCR), in a comprehensive simulation study and compared the performance of our new approach with conventional Bayesian methods for the selection of important predictors in single-response (linear and non-linear) regression models, see holmes2006bayesian, fruhwirth2009improved and dvorzak2016sparse.

We also applied the proposed method to two real data sets. The first data set includes a combination of nine continuous, binary and ordered categorical responses. These responses are phenotypic traits of a rare disorder called Ataxia-Telangiectasia. We analyzed one of the largest cohorts of patients, consisting of $46$ individuals affected by the disease (schon2018genotype). Our model borrows information across the outcomes in order to identify important associations between the responses and a set of genetic and immunological predictors that have been collected in the same study. The second data set consists of four counts and one ordered categorical response which are measured in $122$ individuals suffering from Temporal Lobe Epilepsy (johnson2016systems). Our interest lies in the identification of associations between the responses and $162$ correlated genes that have been identified in a recent gene-network analysis related to cognition abilities and epilepsy (johnson2016systems). Finally, in both real data sets, we compared the predictive ability of the BVSGCR model with widely-used single-response linear and non-linear sparse regression models.

The rest of the paper is organized as follows. In Section 2 we provide a brief presentation of the GCR model and the prior distributions on the regression coefficients and the correlation structure. In Section 3 we describe the novel MCMC algorithm that we propose for BVS when a combination of discrete and/or continuous responses are considered. Section LABEL:sec:simstudy presents the results of the simulation study and in Section LABEL:sec:real we apply the proposed model on two real data sets with missing values in the outcome variables which led to a straightforward modification of the designed MCMC algorithm. Finally, in Section LABEL:sec:discussion we conclude with a short discussion.

2 Gaussian copula regression model

In the following, all vectors, in bold font, are understood as column vectors and the superscript “$T$”is used to denote the transpose of a vector or a matrix. Matrices are also indicated in bold font. The lower-case notation will be reserved for the observations with the corresponding random variables in capital letters.

2.1 Gaussian copulas

An $m$-variate function $C(u_{1},\ldots,u_{m})$, where $C:[0,1]^{m}\rightarrow[0,1],$ is called a copula if it is a continuous distribution function and each marginal is a uniform distribution function on $[0,1]$. sklar1959fonctions proved that any joint cumulative distribution function (cdf) of continuous random variables can be completely specified by its marginal distributions and a unique copula $C$. If $F_{1}(\cdot),\ldots,F_{m}(\cdot)$ are the marginal cdfs of a combination of $m$ continuous and discrete random variables $Y_{1},\ldots,Y_{m}$, their joint cdf can be specified through a specific copula function $C$ as

[TABLE]

A copula function which is commonly used for modeling the dependence structure of any combinations of continuous and discrete variables is the Gaussian copula, see for example hoff2007extending and murray2013bayesian. The Gaussian copula $C$ is specified through the function

[TABLE]

where $\Phi_{m}(\cdot;\bm{R})$ is the cdf of an $m$-variate Gaussian distribution with zero mean vector and correlation matrix $\bm{R}$ and $\Phi^{-1}(\cdot)$ is the inverse of the univariate standard normal cdf. Thus, taking in eq. (1) $u_{k}=F_{Y_{k}}(y_{k})$, for each $k=1,\ldots,m$, we specify the cdf of $\bm{Y}=(Y_{1},\ldots,Y_{m})$ to be the Gaussian copula function. xue2000multivariate proves that the density of the Gaussian copula is

[TABLE]

where $\tilde{\bm{z}}$ is an $m$-dimensional vector and $\tilde{z}_{k}=\Phi^{-1}(u_{k})$, $k=1,\ldots,m$, is known as the normal score which follows marginally the standard Gaussian distribution. Then, $\widetilde{\bm{Z}}$ is an $m$-variate Gaussian distribution with zero mean vector and correlation matrix $\bm{R}$, i.e., $\widetilde{\bm{Z}}\sim\operatorname{N}_{m}(\bm{0},\bm{R})$. If all the marginal distributions are continuous, the matrix $\bm{R}$ can be interpreted as the correlation matrix of the elements of $\bm{Y}$ and zeros in its inverse imply the conditional independence among the corresponding elements of $\bm{Y}$. However, in the presence of discrete random variables, the notion of conditional independence has to be interpreted with care; zeros in $\bm{R}^{-1}$ imply that the observed variables are independent conditionally only on the latent variables (webb2008bayesian). Note also that if $\bm{R}$ is the identity matrix the elements of $\bm{Y}$ can be considered to be independent despite the presence of discrete variables, see xue2000multivariate, song2009joint and talhouk2012efficient for a detailed discussion.

2.2 Regression model

Let $\bm{Y}=(\bm{y}_{1}^{T},\ldots,\bm{y}_{n}^{T})$ be the $(n\times m)$-dimensional matrix of observed data where each $\bm{y}_{i}=(y_{i1},\ldots,y_{im})$ consists of any combination of $m$ discrete and/or continuous responses. Moreover, let $\bm{x}_{ik}$ be the $p_{k}$-dimensional vector of predictors for the $i$th sample and the $k$th response. To model the marginal distribution of each $y_{ik}$, we specify its cdf $F_{k}(y_{ik};\bm{x}_{ik},\bm{\beta}_{k},\bm{\theta}_{k})$ to be the cdf of any parametric distribution. Our notation emphasizes the dependence of each $y_{ik}$ on the $p_{k}$-dimensional response-specific vector of predictors $\bm{x}_{ik}$, the associated regression coefficients $\bm{\beta}_{k}$ and the response-specific parameters $\bm{\theta}_{k}$.

The GCR model is described by the transformation

[TABLE]

where $\tilde{z}_{i1},\ldots,\tilde{z}_{im}$ are realizations from $\widetilde{\bm{Z}}_{i}$ and $h_{ik}(\cdot)=\Phi^{-1}\{F_{k}(\cdot;\bm{x}_{ik},\bm{\beta}_{k},\bm{\theta}_{k})\}$ for each $i=1,\ldots,n$ and $k=1,\ldots,m$.

By assuming that each $F_{k}(\cdot)$ is a member of the exponential family, we obtain the multivariate Generalized Linear Model presented in song2009joint as the multivariate extension of the well-known single-response Generalized Linear Model (mccullagh1989generalized) with $\bm{\theta}_{k}$ the specific vector of parameters for the $k$th response. The SUR model is a special case of eq. (3) when all margins are univariate Gaussian with mean $\bm{x}_{ik}^{T}\bm{\beta}_{k}$ and variance $\theta_{k}$ and $\bm{R}$ is the correlation matrix. The multi-response Probit regression model of chib1998analysis is obtained from eq. (3) by specifying each margin to be the cdf of a Bernoulli random variable with probability of success $\Phi(\bm{x}_{ik}^{T}\bm{\beta}_{k})$. Finally, setting $F_{k}(y_{ik};\bm{x}_{ik},\bm{\beta}_{k},\bm{\theta}_{k})=\Phi(\theta_{kc}-\bm{x}_{ik}^{T}\bm{\beta}_{k})$, $c=1,\ldots,C_{k}-1$, with $C_{k}$ the number of categories, eq. (3) becomes a regression model for a combination of binary and ordinal observations and $\bm{\theta}_{k}=(\theta_{k1},\ldots\theta_{k,C_{k}-1})$ consists of the cut-points for the $k$th ordinal observation (mccullagh1980regression). In both multi-response Probit and ordinal regression models the matrix $\bm{R}$ is in the correlation form for identifiability conditions (chib1998analysis).

Eq. (3) implies that the joint likelihood function of the observations $\bm{Y}$ conditional on the correlation matrix $\bm{R}$, the regression coefficients $\bm{B}=\operatorname{vec}(\bm{\beta}_{1},\ldots,\bm{\beta}_{m})$ and the parameters $\bm{\Theta}=\operatorname{vec}(\bm{\theta}_{1},\ldots,\bm{\theta}_{m})$ is an intractable function of non-computable high-dimensional integrals (song2009joint).

2.3 Prior distributions

In this section, we specify the prior distributions on the regression coefficients $\bm{B}$, the associated sparsity prior on the inclusion probability and the correlation matrix $\bm{R}$ of the model. By noting that for different choices of the marginal cdfs $F_{k}(\cdot)$, $k=1,\ldots,m$, we have different vectors of parameters $\bm{\theta}_{k}$, we will assign their prior distributions differently in each simulated and real data example. See Supp. Mat. Section S.1 for details.

2.3.1 Variable selection

We use a hierarchical non-conjugate model to assign a prior distribution on the regression coefficients of the GCR model defined in eq. (3). A point mass at zero is specified on the regression coefficients of the unimportant predictors, whereas a Gaussian distribution is assigned to the non-zero regression coefficients (george1997approaches). By utilizing the binary latent vector $\bm{\gamma}_{k}=(\gamma_{k1},\ldots,\gamma_{kp_{k}})$, $k=1,\ldots,m$, where $\gamma_{kj}$ is $1$ if, for the $j$th predictor and the $k$th response, the regression coefficient is different from zero and [math] otherwise, we assume that for each $k$

[TABLE]

[TABLE]

[TABLE]

where $\delta_{0}$ denotes a point mass at zero and $v$ is a fixed value. It is common practice to standardize the predictor variables, taking $v=1$ in order to place appropriate prior mass on reasonable values of the non-zero regression coefficients (hans2007shotgun). Integrating out $\pi_{k}$, it is readily shown that marginally

[TABLE]

where $|\bm{\gamma}_{k}|=\sum_{j=1}^{p_{k}}\gamma_{kj}$. The hyper-parameters $a_{k}$ and $b_{k}$ can be chosen using prior information about the number of important covariates associated with the $k$th response and its variance. See (kohn2001nonparametric) for further details on the elicitation of the hyper-parameters $a_{k}$ and $b_{k}$ in eq. (5).

The sparsity prior (5) is a common choice in single-response BVS as well as in the sparse SUR model (wang2010sparse). When the predictors are common to all responses, the primary inferential question may shift to the identification of key predictors that exert their influence to a large fraction of responses at the same time. When a large number of responses are regressed independently on the same set of predictors, Richardson2010 proposed to decompose the a priori “cell” inclusion probability into its marginal effects, i.e., $\pi_{kj}=\pi_{k}\times\rho_{j}$, $\pi_{kj}\in[0,1]$. The idea behind this decomposition is to control the level of sparsity for each response $k$ through a suitable choice of the hyperparameters $a_{k}$ and $b_{k}$ as in (5), while $\rho_{j}\sim\mathrm{Ga}(c,d)$ captures the relative “propensity” of predictor $j$ to influence several responses at the same time. For a detailed discussion regarding the effect of the prior $p(\rho_{j})$ on BVS in multiple-response regression models, see Ruffieux2020.

In the following, for ease of notation, we denote by $\bm{\beta}_{\bm{\gamma}_{k}}$ all non-zero elements of $\bm{\beta}_{k}$ and analogously, we indicate by $\bm{x}_{i,\bm{\gamma}_{k}}$, the elements of $\bm{x}_{ik}$ corresponding to those elements of ${\bm{\gamma}_{k}}$ equal to 1. Let $\bm{X}_{k}$ be the $n\times p_{k}$ matrix of all available predictors for the $k$th response. Similarly to the vector case, $\bm{X}_{\bm{\gamma}_{k}}$ denotes the $n\times|\bm{\gamma}_{k}|$ matrix in which the columns are selected according to the latent binary vector $\bm{\gamma}_{k}$. Finally, we refer to $\bm{\Gamma}=(\bm{\gamma}_{1},\ldots,\bm{\gamma}_{m})$ as the latent binary matrix.

2.3.2 Correlation matrix R and adjacency matrix G

To specify a prior distribution on the correlation matrix $\bm{R}$, we follow talhouk2012efficient and use the parameter expansion for data augmentation strategy (liu1999parameter) to expand the correlation matrix into a covariance matrix. This choice allows us to specify conjugate prior distributions on the resulting covariance matrix. In particular, we utilize the hyper-inverse Wishart (dawid1993hyper) to allow for zero entries in the inverse of the covariance matrix. Details can be summarized as follows.

First, to expand $\bm{R}$ into a covariance matrix, we define the transformation $\widetilde{\bm{W}}=\widetilde{\bm{Z}}\bm{D}$, where $\widetilde{\bm{Z}}$ is the $n\times m$ matrix of the Gaussian latent variables and $\bm{D}$ is an $m\times m$ diagonal matrix with elements $\delta_{k}$, $k=1,\ldots,m$. Then,

[TABLE]

where $\bm{\Sigma}=\bm{D}\bm{R}\bm{D}$ and $\bm{I}_{n}$ is a diagonal matrix of dimension $n$ that encodes the independence assumption amongst the observations. Then, a conjugate prior distribution can be assigned on $\bm{\Sigma}$ and updated at each iteration of the MCMC algorithm before it is projected back to $\bm{R}$ using the inverse transformation $\bm{R}=\bm{D}^{-1}\bm{\Sigma}\bm{D}^{-1}$. We also utilize the theory of decomposable models to perform a conjugate analysis of the covariance structure of the model since the hyper-inverse Wishart distribution is a conjugate prior distribution for the covariance matrix $\bm{\Sigma}$ with respect to the adjacency matrix $\bm{G}$ of a decomposable graph $\mathcal{G}$. The diagonal elements in the adjacency matrix are always restricted to be 1 to ensure the positive definiteness of $\bm{G}$.

Second, we assign the following prior structure on $\bm{G}$ and $\bm{D}$. Let $g_{\ell}$, $\ell=1,\ldots,m(m-1)/2$, be the binary indicator for the presence of the $\ell$th off-diagonal edge in the lower triangular part of the adjacency matrix $\bm{G}$ of the decomposable graph $\mathcal{G}$. We assume that

[TABLE]

For a detailed discussion regarding compatible priors on $\bm{G}$ for decomposable graphs, see Bornn2011.

We denote by $p(\bm{G})$ the induced marginal prior on the adjacency matrix and define the joint distribution of $\bm{D}$, $\bm{R}$ and $\bm{G}$ as

[TABLE]

where

[TABLE]

with $r^{kk}$ the $k$th diagonal element of $\bm{R}^{-1}$. Assuming a uniform prior for $\bm{R}|\bm{G}$, it can be shown that, for a decomposable graph $\mathcal{G}$, $\bm{\Sigma}|\bm{G}\sim\mathrm{H}\mathrm{I}\mathrm{W}_{\mathcal{G}}(2,\bm{I}_{m})$ (talhouk2012efficient).

3 MCMC sampling strategy

We are interested in sampling from the joint posterior distribution $p(\bm{B},\bm{\Gamma},\bm{\Theta},\widetilde{\bm{Z}},\bm{D},\bm{R},\bm{G}|\bm{Y})$. To draw samples from the specified model, we design a novel MCMC algorithm which proceeds as follows. We first update the regression coefficients $\bm{B}$, the latent binary matrix $\bm{\Gamma}$, the parameters $\bm{\Theta}$ and the Gaussian latent variables $\widetilde{\bm{Z}}$ in $m$ blocks. In particular, for each $k=1,\ldots,m$, we sample jointly $(\bm{\beta}_{k},\bm{\gamma}_{k})$, then we update $\bm{\theta}_{k}$ and, finally, we draw realizations from the Gaussian latent variables $\widetilde{\bm{Z}}_{k}=(\widetilde{Z}_{1k},\ldots,\widetilde{Z}_{nk})$. In the last step, we draw the correlation matrix $\bm{R}$ and the adjacency matrix $\bm{G}$ from their full conditional distributions.

Algorithm 1 presents the pseudo code of the designed MCMC algorithm where the subscript “$-k$” implies that the corresponding matrix (or vector) consists of all the elements except those that are related to the $k$th response.

For the large majority of responses, drawing $(\bm{\beta}_{k},\bm{\gamma}_{k})$ in step 4 is complicated due to the unavailability of the full conditional distributions in closed-form expression and a M-H step is therefore required. However, because of the uncertainty associated with the latent binary vector $\bm{\gamma}_{k}$, the proposal distribution requires a different dimension at each MCMC iteration. In this framework, a commonly used tool is the Reversible Jump algorithm (green1995reversible), although it may experience a low acceptance rate, resulting in a MCMC sampler with poor mixing (brooks2003efficient; Lamnisos2009).

To avoid the Reversible Jump algorithm, we design a proposal distribution defined on a fixed-dimensional space that can be used for the joint update of $\bm{\beta}_{k}$ and $\bm{\gamma}_{k}$. It is worth noticing that we conduct the update of $(\bm{\beta}_{k},\bm{\gamma}_{k})$ by first integrating out the latent variables $\widetilde{\bm{Z}}_{k}$. This accelerates the convergence of the proposed MCMC algorithm since conditioning on $\widetilde{\bm{Z}}_{k}$ induces many restrictions on the admissible values of $(\bm{\beta}_{k},\bm{\gamma}_{k})$, see pitt2006efficient. Drawing samples from the posterior distributions of the remaining parameters and the Gaussian latent variables of the model can be conducted by using standard MCMC algorithms which we also describe briefly in this section.

3.1 Proposal distribution for variable selection

To sample jointly the regression coefficients $\bm{\beta}_{k}$ and the latent binary vector $\bm{\gamma}_{k}$ for each $k=1,\ldots,m$ (step 4 Algorithm 1), we design a M-H step that targets the distribution with density

[TABLE]

where $p(\bm{\beta}_{k}|\bm{\gamma}_{k})$ and $p(\bm{\gamma}_{k})$ are the prior densities on $\bm{\beta}_{k}$ and $\bm{\gamma}_{k}$ as defined in eq. (4) and eq. (6), respectively.

In the case of a continuous response $\bm{y}_{k}$, we have from eq. (2) that

[TABLE]

where we condition on $\widetilde{\bm{Z}}_{-k}$ with an abuse of notation (pitt2006efficient) since $\tilde{z}_{ik}=\Phi^{-1}\{F_{k}(y_{ik};\bm{x}_{i,\bm{\gamma}_{k}},\bm{\beta}_{\bm{\gamma}_{k}},\bm{\theta}_{k})\}$ is a shorter expression of the transformation in eq. (3), $r^{k\ell}=(\bm{R}^{-1})_{k\ell}$ and $f_{k}(y_{ik};\bm{x}_{i,\bm{\gamma}_{k}},\bm{\beta}_{\bm{\gamma}_{k}},\bm{\theta}_{k})$ denotes the probability density function of $F_{k}(y_{ik};\bm{x}_{i,\bm{\gamma}_{k}},\bm{\beta}_{\bm{\gamma}_{k}},\bm{\theta}_{k})$.

If $\bm{y}_{k}$ is a discrete response, we have that

[TABLE]

where $\ell_{ik}=\Phi^{-1}\{F_{k}(y_{ik}-1;\bm{x}_{i,\bm{\gamma}_{k}},\bm{\beta}_{\bm{\gamma}_{k}},\bm{\theta}_{k})\}$ and $u_{ik}=\Phi^{-1}\{F_{k}(y_{ik};\bm{x}_{i,\bm{\gamma}_{k}},\bm{\beta}_{\bm{\gamma}_{k}},\bm{\theta}_{k})\}$ with $\tilde{\mu}_{ik|-k}=\bm{R}_{k,-k}\bm{R}_{-k,-k}^{-1}\tilde{\bm{z}}_{i,-k}$ and $\tilde{\sigma}_{k|-k}^{2}=1-\bm{R}_{k,-k}\bm{R}_{-k,-k}^{-1}\bm{R}_{-k,k}$. Thus, the conditional density in eq. (10) will be intractable for the majority of the continuous and for all the discrete distributions that can be used for the marginal modeling of the $k$th response.

Instead of relying on the Reversible Jump algorithm or the Laplace approximation to sample from eq. (10), we utilize a M-H step with a fixed-dimensional proposal distribution obtained by re-parameterizing the Gaussian latent variables $\widetilde{\bm{Z}}_{k}$. More precisely, for each $i=1,\ldots,n$ and $k=1,\ldots,m$, we set

[TABLE]

where $Z_{ik}$ is marginally normal distributed with mean $\mu_{ik}=\bm{x}_{i,\bm{\gamma}_{k}}^{T}\bm{\beta}_{\bm{\gamma}_{k}}$ and variance either $\sigma^{2}_{k}=\delta_{k}^{2}$ (if continuous) or $\sigma^{2}_{k}=1$ (if discrete). Thus, the vector of realizations $\bm{z}_{k}$ is a sufficient statistics for $(\bm{\beta}_{k},\bm{\gamma}_{k})$, whereas $\tilde{\bm{z}}_{k}$ is ancillary (Yu2011). Irrespectively of the type of the response, eq. (13) creates a link between the Gaussian Copula model, where the Gaussian latent variables are non-centered, and BVS that requires a centered parametrization. With this transformation, we also aim to keep the advantages of using (data augmentation) centered auxiliary variables when performing BVS (holmes2006bayesian). In particular, the fact that the vector $\bm{z}_{k}$ retains information about the likelihood is key since, regardless of the response’s type, it allows an effective update of $\bm{\beta}_{k}$ given a change in the predictors set.

To construct a proposal distribution for a M-H step that targets eq. (10), we replace the intractable density $p(\bm{y}_{k}|\bm{\beta}_{k},\bm{\gamma}_{k},\bm{\theta}_{k},\widetilde{\bm{Z}}_{-k},\bm{R})$ with the Gaussian density of the latent variable $\bm{Z}_{k}$ conditioned on $\widetilde{\bm{Z}}_{-k}$. Exploiting the fact that $\bm{\beta}_{k}$ is quadratic in $p(\bm{z}_{k}|\bm{\beta}_{k},\bm{\gamma}_{k},\widetilde{\bm{Z}}_{-k},\bm{D},\bm{R})$, given a candidate value of the latent binary vector $\bm{\gamma}_{k}^{\ast}$, we propose the non-zero regression coefficients $\bm{\beta}_{k}^{\ast}$ from the distribution with density

[TABLE]

where $\sigma_{k|-k}^{2}$ is the diagonal element of the conditional covariance matrix of $\bm{Z}_{k}|\widetilde{\bm{Z}}_{-k}$, $\bm{V}_{\bm{\gamma}_{k}^{\ast}}^{-1}=\sigma_{k|-k}^{-2}\bm{X}_{\bm{\gamma}^{\ast}_{k}}^{T}\bm{X}_{\bm{\gamma}^{\ast}_{k}}+v^{-1}\bm{I}_{\left|\bm{\gamma}^{\ast}_{k}\right|}$ and $\tilde{\bm{\mu}}_{k|-k}$ is the conditional mean of $\widetilde{\bm{Z}}_{k}|\widetilde{\bm{Z}}_{-k}$. The proposal distribution in eq. (14) shows two important features: (i) in the covariance matrix $\bm{V}_{\bm{\gamma}^{\ast}_{k}}$, the term $\bm{X}_{\bm{\gamma}^{\ast}_{k}}^{T}\bm{X}^{\ast}_{\bm{\gamma}_{k}}$ is rescaled by the conditional variance $\sigma_{k|-k}^{-2}$ and (ii) the proposal mean is shifted by the rescaled conditional mean $\sigma_{k}\tilde{\bm{\mu}}_{k|-k}$, where the factor $\sigma_{k}$ brings $\bm{z}_{k}$ and $\tilde{\bm{\mu}}_{k|-k}$ on the same scale. Thus, the realizations $\bm{z}_{k}$ of the Gaussian latent variables $\bm{Z}_{k}|\widetilde{\bm{Z}}_{-k}$ in eq. (14) are recentered to remove the effect of conditioning on $\widetilde{\bm{Z}}_{-k}$, i.e., $\operatorname{E}(\bm{Z}_{k}|\widetilde{\bm{Z}}_{-k}-\sigma_{k}\tilde{\bm{\mu}}_{k|-k})=\operatorname{E}(\bm{Z}_{k})=\bm{X}_{\bm{\gamma}_{k}}\bm{\beta}_{\bm{\gamma}_{k}}$, see Supp. Mat. Section S.3 for details.

Remark 1**.**

If $\bm{R}$ is the identity matrix, then the proposal distribution (14) becomes

[TABLE]

since $\bm{V}_{\bm{\gamma}^{\ast}_{k}}^{-1}=\sigma_{k}^{-2}\bm{X}_{\bm{\gamma}^{\ast}_{k}}^{T}\bm{X}_{\bm{\gamma}^{\ast}_{k}}+v^{-1}\bm{I}_{\bm{\gamma}^{\ast}_{k}}$ and the conditional mean $\tilde{\bm{\mu}}_{k|-k}$ is equal to zero. When no information is shared across responses, the proposal distribution in eq. (15) coincides with the full conditional distribution of the non-zero regression coefficients used in single-response Probit ($\sigma^{2}_{k}=1$) regression model (albert1993bayesian; holmes2006bayesian). It is worth noticing that, regardless of the similarities with the Probit model, this proposal density has been obtained as a special case of eq. (14).

Remark 2**.**

For some continuous and discrete data including the Gaussian, binary, ordered and nominal categorical responses, the proposal density in eq. (14) allows the “implicit marginalisation” of the regression coefficients when the joint update of $(\bm{\beta}_{k},\bm{\gamma}_{k})$ in the M-H step is performed. For this type of responses, the acceptance probability $\alpha$ of the joint move is

[TABLE]

where $\bm{m}_{\bm{\gamma}^{\ast}_{k}}=\bm{V}_{\bm{\gamma}^{\ast}_{k}}\bm{X}_{\bm{\gamma}^{\ast}_{k}}^{T}\sigma_{k|-k}^{-2}(\bm{z}_{k}-\sigma_{k}\tilde{\bm{\mu}}_{k|-k})$. Similarly to holmes2006bayesian, in eq. (16) the current and proposed value of the regression coefficients $(\bm{\beta}_{k},\bm{\beta}_{k}^{\ast})$ do not appear. For details, see Supp. Mat. Section S.4.

Finally, since the main focus of this paper is to provide a proposal distribution for the regression coefficients given the latent binary vector, any proposal distribution for $\bm{\gamma}_{k}^{\ast}$ can be used. Here, we use a modified version of the proposal distribution designed by guan2011bayesian. This proposal makes efficient use of the inexpensive evaluation of the marginal association between each response and the predictors in order to propose a new latent vector $\bm{\gamma}_{k}^{\ast}$. See Supp. Mat. Section S.2 for a detailed description.

3.2 Sampling the Gaussian latent variables

Next, we describe the sampling strategy for the Gaussian latent variables (step 6 Algorithm 1). If the $k$th response is continuous then the transformation (3) is a one-to-one transformation. In this case $\tilde{\bm{z}}_{k}$ is updated deterministically by setting $\tilde{z}_{ik}=\Phi^{-1}\{F_{k}(y_{ik};\bm{x}_{i,\bm{\gamma}_{k}},\bm{\beta}_{\bm{\gamma}_{k}},\bm{\theta}_{k})\}$ for each $i=1,\ldots,n$ and $k=1,\ldots,m$. For a discrete response $\bm{y}_{k}$, we have that

[TABLE]

where $\tilde{\mu}_{ik|-k}$, $\tilde{\sigma}_{k|-k}^{2}$, $\ell_{ik}$ and $u_{ik}$ are defined in eq. (12). Therefore, each $\tilde{z}_{ik}$ has to be sampled from the Gaussian distribution $\mathrm{N}(\tilde{\mu}_{ik|-k},\tilde{\sigma}_{k|-k}^{2})$ truncated on the interval $(\ell_{ik},u_{ik}]$.

3.3 Sampling the matrix D, the correlation matrix R and the adjacency matrix G

To update $(\bm{D},\bm{R})$ and $\bm{G}$ in the designed MCMC Algorithm 1 (step 8 and 9), we follow talhouk2012efficient and work in the space of the scaled Gaussian latent variables $\widetilde{\bm{W}}=\widetilde{\bm{Z}}\bm{D}$. In practice, we obtain samples from $p(\bm{\Sigma}|\widetilde{\bm{W}},\bm{G})$ and transform them using the inverse transformation $\bm{R}=\bm{D}^{-1}\bm{\Sigma}\bm{D}^{-1}$ where $\bm{D}$ is sampled from its prior distribution in eq. (9). To sample the adjacency matrix $\bm{G}$ from $p(\bm{G}|\bm{Y},\widetilde{\bm{Z}},\bm{D},\bm{R})$, we target the distribution with density

[TABLE]

where the summation in the denominator is over all the decomposable graphs $\mathcal{\bm{G}}$, $p(\bm{G})$ is defined through eq. (8) and $p(\widetilde{\bm{W}}|\bm{G})$ can be computed analytically due to the tractability of the hyper-inverse Wishart distribution, see Supp. Mat. Section S.5 for further details. We sample from eq. (18) by using a M-H step in which, conditionally on the current adjacency matrix $\bm{G}$, a new graph is proposed by adding or deleting an edge between two vertices whose index has been chosen randomly between the vertices that belong to a decomposable graph. The proposed graph is then accepted or rejected using the accept/reject mechanism of the M-H step which targets the density in eq. (18). Finally, conditionally on $\widetilde{\bm{W}}$ and the updated adjacency matrix $\bm{G}$, we sample $\bm{\Sigma}$ from its conditional distribution $\mathrm{H}\mathrm{I}\mathrm{W}_{\mathcal{G}}(2+n,\bm{I}_{m}+\widetilde{\bm{W}}^{T}\widetilde{\bm{W}})$.