Statistics & Gene Expression Data Analysis Note 10: Gibbs Sampling in Binary Regression

A second specific example of Gibbs sampling - canonical Markov chain Monte Carlo method for calculations in Bayesian binary regressions, and a trivial extension of Gibbs in linear regressions with multiple shrinkage.

♠ Probit Model

• Expression level vector x_{j} on array (tumor sample, etc) j = 1, …, n

• Binary outcome: z_{j} = 0 or 1

• Probit probability model: π_{j} = P r(z_{j} = 1) (conditional on chosen predictor variables and model parameters: in full, π_{j} = P r(z_{j} = 1|x_{j}, β))

• Linear regression model based on regression function μ_{j} = x′_{j}β (perhaps x_{j} is extended with a leading 1 to include an intercept term β_{0})

• Probit regression:

π_{j} = Φ(μ_{j})

where Φ is standard normal cumulative distribution function

♠ Prior and Posterior

• Example: β ∼ N(0, C^{−1}) and we’ll take the precision as diagonal

• Posterior for regression parameters

p(β|z) ∝ exp(−β′Cβ) ∏_{j=1}^{n} Φ(μ_{j})^{z_{j}}(1 − Φ(μ_{j}))^{1−z_{j}}

which can be evaluated, numerically optimized (it is unimodal) using standard NR routines,

• Exercise: Write code (in C/C++ or other) to implement a Newton-Raphson search for the mode of the posterior density. As ALWAYS in such problems, work with the log-posterior, i.e., maximise the log of the target function for numerical stability

• Under vague reference prior (C → 0) the posterior is the normalised likelihood function. The MLE and related information can be computer using R/Splus glm functions, Matlab glmfit function, or (quite easily) by user-written code as a special case

♠ Latent Variables and Data Augmentation

• For each sample, recognise an underlying latent variable y_{j}

• π_{j} = P r(y_{j} > 0) when y_{j} ∼ N(μ_{j}, 1) = N(x′_{j}β, 1)

• z_{j} = 1 if and only if y_{j} > 0 - e.g., latent variable is positive for ER+ cases, negative for ER− cases

- could precisely classify cases if we could observe the latent y_{j}, but we do not; result is the binary probability model

- MCMC calculations impute these “missing” values along with values of the parameters β

♠ Conditional Posteriors

• p(β|y, z)

If y were known, we have a linear regression of the y_{j} with regression variables x_{j}, parameter β and error variance 1. Notice that this simplifies (it is a special case) of the linear regression with unknown error variance: we simply set the precision parameter ϕ to be 1, and then the analysis is all the same otherwise. In the MCMC, we just skip the simulation of ϕ values at each iteration, since we know ϕ = 1 here.

The actual value of z is at this point irrelevant – the information they contain is already there in the (current, imputed or candidate) values of y : Formally, β is conditionally independent of z given y. The posterior (see earlier notes) is multivariate normal

β|y, z ∼ 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

• p(y|β, z)

If β were known, then the y_{j} are independent normals but subject to the information provided by the z_{j} - in each case, we just condition the initial normal on the information that y_{j} must be positive (if z_{j} = 1) or negative (if z_{j} = 0). The result may be written simply in terms of the posterior cumulative distribution function, as

P (y_{j}|z_{j} = 1) = [Φ(y_{j} − μ_{j}) − (1 − π_{j})]/π_{j}, for y_{j} > 0,

and

P (y_{j}|z_{j} = 0) = Φ(y_{j} − μ_{j})/(1 − π_{j}), for y_{j} < 0

(check that you can derive this). Or, for any value of z_{j} = 0, 1,

P (y_{j}|z_{j}) = [Φ(y_{j} − μ_{j}) − z_{j}(1 − π_{j})]/[z_{j}π_{j} + (1 − z_{j})(1 − π_{j})].

Simulated values of the y_{j} are then drawn, independently, via the inverse CDF approach: generate u_{j} ∼ U(0, 1) and solve for y_{j} in P (y_{j}|z_{j}) = u_{j}. It can be written trivially as

y_{j} = μ_{j} + Φ^{−1}{z_{j}(1 − π_{j}) + u_{j}(z_{j} + (1 − π_{j})(1 − 2z_{j}))}

with Φ^{−1} being the normal quantile function (inverse CDF - qnorm function).

♠ Including Additional Parameters

MCMC neatly extends to include other parameters. A key example is the multiple shrinkage prior model that is of general interest, but of particular relevance in factor regression models (see earlier notes on this).

• Suppose C = diag(γ_{0}, γ_{1}, …, γ_{p}) with prior variances τ_{i} = γ_{i}^{−1} defining element-wise shrinkage parameters for the individual predictor variables. The above discussion all applied now explicitly conditional on values of C, so that the simulation iterations can run with simulations of C coupled in too. That requires priors on elements of C; if these are independent gamma priors, γ_{j} ∼ Ga(k/2, h/2) for each j, say, then the relevant conditional posteriors are also independent gammas, namely

γ_{j} ∼ 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.