Statistics & Gene Expression Data Analysis Note 9: Gibbs Sampling in Linear Regression

A specific example of Gibbs sampling - canonical Markov chain Monte Carlo method for calculations in Bayesian statistical models, as well as a key applied model class.

♠ Simulation of Posteriors

• For any statistical model, simulated values of posterior distributions are nowadays the standard in statistical computation: summarise large samples of parameter values from a posterior distribution to easily understand the information contained in that posterior about the parameters

• MCMC methods, such as Gibbs sampling, generate simulations sequentially using Markov chains

• Gibbs sampling represents an approach in which sets of coupled conditional posteriors derived from the model are used for these simulations

• Here’s a key example

♠ Prior and Posterior in the Linear Regression Model: Multiple Shrinkage Example

• Example: Now take β ∼ N(0, C^{−1}) and we’ll take the precision as diagonal, with diagonal elements (multiple shrinkage factors) now to be estimated along with the regression parameter and precision ϕ. Note that this is a bit different to the conditionally conjugate prior of the standard regression earlier – the error precision is no longer included in the prior for β. This is a modification of the analysis that relaxes that strict adoption of a traditional conjugate prior – we are freed from the need to impose strict conjugacy now that we’re using MCMC methods.

• Write C^{−1} = diag(τ_{1}, τ_{1}, …, τ_{p}) with prior variances τ_{i} defining element-wise shrinkage parameters for the individual predictor variables.

• Now allow for learning on the shrinkage factors using priors that are independent gamma priors, τ_{i}^{−1} ∼ Ga(k/2, h/2) for each j, say, then the relevant conditional posteriors are all available in conjugate forms.

♠ Conditional Posteriors

• p(β|y, ϕ, C)

We have a linear regression so that the standard results follow: The posterior (see earlier notes) is multivariate normal

β|y, ϕ, C ∼ N(b, B^{−1})

with B = C + H′Hϕ and b = B^{−1}ϕH′y where H = X′ is the n × p design matrix. For multivariate normals, use programmed functions (multivariate normal simulation is very standard - e.g.. rMNorm.m or similar) or use direct Cholesky decomposition: e.g., by hand in matlab

b + chol(inv(B)) ∗ randn(p, 1)

for a single draw

• If we know y, β and C, then the conditional posterior for the error precision ϕ depends only on the data and β through the residual ϵ = y − Hβ, and the relevant conditional posterior for sampling values from is just the updated gamma distribution ϕ|y, β, C ∼ Ga((a + n)/2, (b + ϵ′ϵ)/2).

• The relevant conditional posteriors are also independent gammas, namely

τ_{j}^{−1} ∼ Ga((k + 1)/2, (h + β_{j}^{2})/2).

This allows for learning on differential shrinkage parameters across variables.

• Exercise: Develop the MCMC with these priors.

• Think about choices of prior parameters: One general way to think about ranges of relevant parameter values in binary regression models is to consider how variation in μ_{j} translates through to the probability scale π_{j}. Absolute values of μ_{j} bigger than 2 or so lead to probabilities that are already very extreme.

Mike West  June 2, 2004  Page 0