MCMC Method for Arbitrary Missing Data

The MI Procedure

MCMC Method for Arbitrary Missing Data

The Markov Chain Monte Carlo (MCMC) method originated in physics as a tool for exploring equilibrium distributions of interacting molecules. In statistical applications, it is used to generate pseudo-random draws from multidimensional and otherwise intractable probability distributions via Markov chains. A Markov chain is a sequence of random variables in which the distribution of each element depends only on the value of the previous one.

In MCMC simulation, one constructs a Markov chain long enough for the distribution of the elements to stabilize to a stationary distribution, which is the distribution of interest. By repeatedly simulating steps of the chain, the method simulates draws from the distribution of interest. Refer to Schafer (1997) for a detailed discussion of this method.

In Bayesian inference, information about unknown parameters is expressed in the form of a posterior probability distribution. This posterior distribution is computed using Bayes' theorem

$p({{\theta}} | y)=\frac {p( y | {{\theta}}) p( {{\theta}})} { \int p( y | {{\theta}}) p({{\theta}}) d{{\theta}} }$

MCMC has been applied as a method for exploring posterior distributions in Bayesian inference. That is, through MCMC, one can simulate the entire joint posterior distribution of the unknown quantities and obtain simulation-based estimates of posterior parameters that are of interest.

In many incomplete data problems, the observed-data posterior $p({{\theta}} | Y_{obs})$ is intractable and cannot easily be simulated. However, when Y_obs is augmented by an estimated/simulated value of the missing data Y_mis, the complete-data posterior $p({{\theta}} | Y_{obs}, Y_{mis})$ is much easier to simulate. Assuming that the data are from a multivariate normal distribution, data augmentation can be applied to Bayesian inference with missing data by repeating the following steps:

1. The imputation I-step: Given an estimated mean vector and covariance matrix, the I-step simulates the missing values for each observation independently. That is, if you denote the variables with missing values for observation i by Y_i(mis) and the variables with observed values by Y_i(obs), then the I-step draws values for Y_i(mis) from a conditional distribution for Y_i(mis) given Y_i(obs).

2. The posterior P-step: Given a complete sample, the P-step simulates the posterior population mean vector and covariance matrix. These new estimates are then used in the next I-step. Without prior information about the parameters, a noninformative prior is used. You can also use other informative priors. For example, a prior information about the covariance matrix can be helpful to stabilize the inference about the mean vector for a near singular covariance matrix.

The two steps are iterated long enough for the results to be reliable for a multiply imputed data set (Schafer 1997, p. 72). That is, with a current parameter estimate ${{\theta}}^{(t)}$ at the tth iteration, the I-step draws Y_mis^(t+1) from $p( Y_{mis} | Y_{obs}, {{\theta}}^{(t)} )$ and the P-step draws ${{\theta}}^{(t+1)}$ from $p( {{\theta}} | Y_{obs}, Y_{mis}^{(t+1)} )$ .

This creates a Markov chain

$(Y_{mis}^{(1)}, {{\theta}}^{(1)})$ , $(Y_{mis}^{(2)}, {{\theta}}^{(2)})$ , ... ,

which converges in distribution to $p( Y_{mis}, {{\theta}} | Y_{obs})$ .Assuming the iterates converge to a stationary distribution, the goal is to simulate an approximately independent draw of the missing values from this distribution.

To validate the imputation results, you should repeat the process with different random number generators and starting values based on different initial parameter estimates.

The next three sections provide details for the imputation step, Bayesian estimation of the mean vector and covariance matrix, and the posterior step.

Imputation Step

In each iteration, starting with a given mean vector

and covariance matrix ${\Sigma}$ , the imputation step draws values for the missing data from the conditional distribution Y_mis given Y_obs.

Suppose $mu=[ mu_{1}', mu_{2}' ] '$ is the partitioned mean vector of two sets of variables, Y_obs and Y_mis, where $mu_{1}$ is the mean vector for variables Y_obs and ${mu}_{2}$ is the mean vector for variables Y_mis.

Also suppose

${\Sigma} &=& [ {\Sigma}_{11} & {\Sigma}_{12} \ {\Sigma}_{12}' & {\Sigma}_{22} \ ]$

is the partitioned covariance matrix for these variables, where ${\Sigma}_{11}$ is the covariance matrix for variables Y_obs, ${\Sigma}_{22}$ is the covariance matrix for variables Y_mis, and ${\Sigma}_{12}$ is the covariance matrix between variables Y_obs and variables Y_mis.

By using the sweep operator (Goodnight 1979) on the pivots of the ${\Sigma}_{11}$ submatrix, the matrix becomes

$[ {\Sigma}_{11}^{-1} & {\Sigma}_{11}^{-1} {\Sigma}_{12} \ -{\Sigma}_{12}' {\Sigma}_{11}^{-1} & {\Sigma}_{22.1} \ ]$

where ${\Sigma}_{22.1}={\Sigma}_{22} - {\Sigma}_{12}' {\Sigma}_{11}^{-1} {\Sigma}_{12}$ can be used to compute the conditional covariance matrix of Y_mis after controlling for Y_obs.

For an observation with the preceding missing pattern, the conditional distribution of Y_mis given Y_obs = y₁ is a multivariate normal distribution with the mean vector

${mu}_{2.1}={mu}_{2} + {\Sigma}_{12}' {\Sigma}_{11}^{-1} ( y_{1} - {mu}_{1} )$

and the conditional covariance matrix

${\Sigma}_{22.1}={\Sigma}_{22} - {\Sigma}_{12}' {\Sigma}_{11}^{-1} {\Sigma}_{12}$

Bayesian Estimation of the Mean Vector and Covariance Matrix

Suppose that Y = ( y₁', y₂', ..., y_n' )' is an (n×p) matrix made up of n (p×1) independent vectors y_i, each of which has a multivariate normal distribution with mean zero and covariance matrix ${\Lambda}$ .Then the SSCP matrix

$A=Y'Y=\sum_{i} { y_{i} y_{i}' }$

has a Wishart distribution $W( n, {\Lambda})$ .

When each observation y_i is distributed with a multivariate normal distribution with an unknown mean ,then the CSSCP matrix

$A=\sum_{i} { ( y_{i} - {\overline y} ) ( y_{i} - {\overline y} )' }$

has a Wishart distribution $W( n-1, {\Lambda})$ .

If A has a Wishart distribution $W( n, {\Lambda})$ ,then B = A^-1 has an inverted Wishart distribution $W^{-1}( n, {\Psi})$ , where n is the degrees of freedom and ${\Psi}={\Lambda}^{-1}$ is the precision matrix (Anderson 1984).

Note that, instead of using the parameter ${\Psi}={\Lambda}^{-1}$ for the inverted Wishart distribution, Schafer (1997) uses the parameter ${\Lambda}$ .

Suppose that each observation in the data matrix Y has a multivariate normal distribution with mean and covariance matrix ${\Sigma}$ .Then with a prior inverted Wishart distribution for ${\Sigma}$ and a prior normal distribution for

${\Sigma} \sim & & W^{-1} ( \, m, \, {\Psi} ) \ {mu} | {\Sigma} \sim & & N ( \, {mu}_{0}, \,\, \frac 1{\tau} {\Sigma})$

where $\tau \gt 0$ is a fixed number. The posterior distribution (Anderson 1984, p. 270; Schafer 1997, p. 152) is

${\Sigma} | Y \sim & & W^{-1} ( \, n+m, \,\, (n-1) S + {\Psi} + \frac{n \tau}{... ...rac{1}{n+\tau} \, (n \bar{y} + \tau {mu}_{0}), \,\, \frac 1{n+\tau} {\Sigma} )$

where (n-1) S is the CSSCP matrix.

Posterior Step

In each iteration, the posterior step simulates the posterior population mean vector

and covariance matrix ${\Sigma}$ from prior information for

and ${\Sigma}$ , and the complete sample estimates.

You can specify the prior parameter information using one of the following methods:

PRIOR=JEFFREYS, which uses a noninformative prior.
PRIOR=INPUT=, which provides a prior information for ${\Sigma}$ in the data set. Optionally, it also provides a prior information for in the data set.
PRIOR=RIDGE=, which uses a ridge prior.

The next four subsections provide details of the posterior step for different prior distributions.

1. A Noninformative Prior

Without prior information about the mean and covariance estimates, a noninformative prior can be used by specifying the PRIOR=JEFFREYS option. The posterior distributions (Schafer 1997, p. 154) are

${\Sigma}^{(t+1)} | Y \sim & & W^{-1} ( \, n-1, \, (n-1)S ) \ {mu}^{(t+1)} | ( ... ...{(t+1)}, Y) \sim & & N ( \, \bar{y}, \,\, \frac 1{\, n \,} {\Sigma}^{(t+1)} )$

2. An Informative Prior for and ${\Sigma}$

When prior information is available for the parameters

and ${\Sigma}$ , you can provide it with a SAS data set that you specify with the PRIOR=INPUT= option.

${\Sigma} \sim & & W^{-1} ( \, d^{*}, \, d^{*}S^{*} ) \ {mu} | {\Sigma} \sim & & N ( \, {mu}_{0}, \,\, \frac 1{n_{0}} {\Sigma} )$

To obtain the prior distribution for ${\Sigma}$ ,PROC MI reads the matrix S^* from observations in the data set with _TYPE_=`COV', and it reads n^*=d^*+1 from observations with _TYPE_=`N'.

To obtain the prior distribution for ,PROC MI reads the mean vector $mu_{0}$ from observations with _TYPE_=`MEAN', and it reads n₀ from observations with _TYPE_=`N_MEAN'. When there are no observations with _TYPE_=`N_MEAN', PROC MI reads n₀ from observations with _TYPE_=`N'.

The resulting posterior distribution, as described in the "Bayesian Estimation of the Mean Vector and Covariance Matrix" section, is given by

${\Sigma}^{(t+1)} | Y \sim & & W^{-1} ( \, n+d^{*}, \,\, (n-1)S + d^{*}S^{*} ... ...{0}} \, (n \bar{y} + n_{0} {mu}_{0}), \,\, \frac 1{n+n_{0}} {\Sigma}^{(t+1)} )$

where

$S_{m}=\frac{n n_{0}}{n + n_{0}} (\bar{y} - {mu}_{0}) (\bar{y} - {mu}_{0})'$

3. An Informative Prior for ${\Sigma}$

When the sample covariance matrix S is singular or near singular, prior information about ${\Sigma}$ can also be used without prior information about

to stabilize the inference about

.You can provide it with a SAS data set that you specify with the PRIOR=INPUT= option.

To obtain the prior distribution for ${\Sigma}$ ,PROC MI reads the matrix S^* from observations in the data set with _TYPE_=`COV', and it reads n^* from observations with _TYPE_=`N'.

Note that if the PRIOR=INPUT= data set also contains observations with _TYPE_=`MEAN', then a complete informative prior for both and ${\Sigma}$ will be used.

Corresponding to the prior for ${\Sigma}$

${\Sigma} \sim & & W^{-1} ( \, d^{*}, \, d^{*}S^{*} ) \$

the posterior distribution for ${\Sigma}$ (Anderson 1984, p. 269) is

${\Sigma}^{(t+1)} | Y \sim & & W^{-1} ( \, (n-1)+d^{*}, \,\, (n-1)S + d^{*}S^{*} ) \$

Thus, an estimate of ${\Sigma}$ is given by the weighted average

[1/((n-1) + d^*)] ( (n-1) S + d^* S^* )

and the posterior distribution for $( {mu}, {\Sigma})$ becomes

${\Sigma}^{(t+1)} | Y \sim & & W^{-1} ( \, (n-1)+d^*, \,\, (n-1) S + d^* S^{*}... ... Y ) \sim & & N ( \, \bar{y}, \,\, \frac 1{\,\, n \,\,} \, {\Sigma}^{(t+1)} )$

4. A Ridge Prior

A special case of the preceding adjustment is a ridge prior with S^* = Diag S (Schafer 1997, p. 156). That is, S^* is a diagonal matrix with diagonal elements equal to the corresponding elements in S.

You can request a ridge prior by using the PRIOR=RIDGE= option. You can explicitly specify the number $d^{*} \geq 1$ in the PRIOR=RIDGE=d^* option. Or you can implicitly specify the number by specifying the proportion p in the PRIOR=RIDGE=p option to request d^*= (n-1) p.

The posterior is then given by

${\Sigma}^{(t+1)} | Y \sim & & W^{-1} ( \, (n-1) + d^{*}, \,\, (n-1) S + d^{*}... ...}, Y ) \sim & & N ( \, \bar{y}, \, \frac 1{\,\, n \,\,} \, {\Sigma}^{(t+1)} )$

Chapter Contents
Previous
Next
Top