Gibbs Sampling for Censored Regression

Gibbs sampling is very convenient for many "missing data" problems. To illustrate this situation, suppose we have the simple regression model

,

where the errors are iid . The problem is that some of the response variables are right-censored and we actually observe , where is a known censoring value. We know which observations are uncensored () and which observations are censored ().

We illustrate this situation by use of a picture inspired by a similar one in Martin Tanner's book. We show a scatterplot of the covariate against the observed . The points highlighted in red correspond to censored observations -- the actual observations (the ) exceed the censored values, which we indicate by arrows.


To apply Gibbs sampling to this situation, we imagine a complete data set where all of the 's are known. The unknowns in this problem are the regression coefficients , the error variance , and the 's corresponding to the censored observations.

The joint posterior of all unobservables (assuming the usual vague prior for regression) has the form

where is equal to 1 if the observation is not censored, and if the observation is censored at .

With the introduction of the complete data set, Gibbs sampling is straightforward.

Suppose one has initial estimates at the regression coefficients and the error variance. Then

(1) One simulates from the distribution of the missing data (the for the censored observations) given the parameters. Specifically is simulated from a normal() distribution censored below by .

(2) Given the complete data, one simulates values of the parameters using the usual approach for a normal linear regression model.

To do this on R, I have a couple of useful tools. The function rtruncated.R will simulate draws from an arbitrary truncated distribution

rtruncated=function(n,lo,hi,pf,qf,...)
qf(pf(lo,...)+runif(n)*(pf(hi,...)-pf(lo,...)),...)

For example, if one wishes to simulate 20 draws of a normal(mean = 10, sd = 2) distribution that is truncated below by 4, one writes

rtruncated(20, 4, Inf, pnorm, qnorm, 10, 2)

Also the function blinreg.R in the LearnBayes package will simulate draws of a regression coefficient and the error variance for a normal linear model with a noninformative prior.