Multiple imputation & Bayesian analysis

Let \(Y\) denote a complete observation with density \(p_Y(y\mid\theta)\), and let \(D=(Y_1,\ldots,Y_n)\) be an independent sample. Both frequentist and Bayesian analyses begin with the same sampling model and likelihood:

The difference lies in how uncertainty about \(\theta\) is represented. A frequentist analysis regards the true value \(\theta_0\) as fixed and studies the behavior of an estimator across repeated samples. A Bayesian analysis specifies a prior density \(q(\theta)\) and updates it after observing \(D\). The prior expresses uncertainty about the parameter before these observations are incorporated.

Bayes’ theorem combines the likelihood and prior to give the posterior density:

The denominator normalizes the product so that the posterior integrates to one. It is also the marginal density of the observed dataset,

This integral often has no convenient closed form. Nevertheless, the unnormalized product \(L(D\mid\theta)q(\theta)\) is enough to identify many familiar posterior families and to construct simulation algorithms. Once the posterior has been obtained, inference proceeds by summarizing its location, spread, and shape.

Different posterior summaries answer different questions. A mean or median provides a measure of location, while the variance and quantiles describe uncertainty. For example, the posterior mean is

The posterior mode is the value at which the posterior density is largest:

Computing these summaries requires a specified prior. A common construction selects a parametric family,

in which \(\gamma\) is a hyperparameter controlling features such as the prior mean or spread. Choices of \(\gamma\) therefore belong to the model specification and should be distinguished from the unknown parameter \(\theta\) about which inference is sought.

Bayesian estimators can also be studied from a repeated-sampling perspective. Consider the posterior mode. Because the marginal density of the data does not depend on the candidate value of \(\theta\), its defining optimization can be written as

In a regular model, the curvature of the log likelihood increases with the sample size because it is a sum of independent contributions. A fixed, sufficiently smooth prior then has diminishing influence near the likelihood maximum, and the posterior mode approaches the MLE. Consistency and efficiency follow under the same kind of identification and regularity conditions needed for the likelihood analysis, together with suitable prior support. The local expansion below makes those qualifications explicit.

The normal approximation to a posterior is a local approximation, not a consequence of Bayes’ formula alone. In a regular, identifiable finite-dimensional model, with the true parameter in the interior and a prior density positive and continuous near it, a second-order expansion gives

on \(n^{-1/2}\) neighborhoods of the likelihood maximum. The resulting covariance is approximately \(\{nI(\widehat\theta)\}^{-1}\). At boundaries, with nonidentifiability, or under substantial prior restrictions, this approximation need not hold. Bayesian credible intervals describe posterior probability; repeated-sampling coverage is an additional property to establish under assumptions.

Suppose \(Y_1,\ldots,Y_n\) are conditionally independent Bernoulli observations with success probability \(\theta\), and assign the prior \(\theta\sim\operatorname{Beta}(\alpha,\lambda)\). The positive hyperparameters \(\alpha\) and \(\lambda\) control the prior distribution. Write \(s=\sum_{i=1}^nY_i\) for the observed number of successes. The likelihood and prior density are, respectively,

Multiplying them adds the exponents of \(\theta\) and \(1-\theta\). Terms independent of \(\theta\) enter only the normalizing constant, leaving the posterior kernel

The exponents identify this kernel as a beta density, so the posterior belongs to the same family as the prior:

Its normalized density follows directly from the beta normalizing constant:

This calculation illustrates why recognizing a kernel is useful. Instead of evaluating the marginal likelihood by integration, we identify a probability family whose integral is already known. The result supplies both the posterior distribution and its summaries, while making the contribution of the observed successes and failures explicit.

The beta–Bernoulli example can also be read as a decomposition of information. With \(a=\alpha\), \(b=\lambda\), and \(s=\sum_iY_i\), the posterior mean is

It averages a prior mean and an empirical mean with weights determined by their effective sample sizes. The posterior variance is

The next Bernoulli observation has success probability \((a+s)/(a+b+n)\) after integrating over the posterior. This predictive averaging will later be performed for missing observations rather than future observations.

Suppose \(Y_1,\cdots, Y_n\) are a random sample from \(N(\mu,\tau^{-1})\). First assume that \(\tau\) is known and \(\mu\) is unknown and that, as the prior for \(\mu\),

This gives

Hence,

Next, assume that \(\mu\) is known and \(\tau\) is unknown and that

i.e., \(q(\tau)\propto \tau^{\delta_0/2-1}\exp(-\gamma_0\tau/2)\). The posterior distribution of \(\tau\) is

It follows that

Finally, assume that both \(\mu\) and \(\tau\) are unknown. Specify the following prior \(q(\mu,\tau)= q(\mu\mid \tau)q(\tau)\), where

and

It can be shown that

Write \(\kappa_0=\tau_0\) for the prior precision multiplier, reserving \(\tau\) for the sampling precision. Define

To obtain these quantities, expand the two sums of squares and complete the square in \(\mu\):

Consequently \(\tau\mid D\sim\operatorname{Gamma}(a_n,b_n)\) and \(\mu\mid\tau,D\sim N\{\mu_n,(\kappa_n\tau)^{-1}\}\). Integrating out \(\tau\) gives a Student \(t\) distribution for \(\mu\) with \(2a_n\) degrees of freedom, location \(\mu_n\), and scale squared \(b_n/(a_n\kappa_n)\). The predictive distribution for a new \(Y\) is also Student \(t\), with scale squared \(b_n(\kappa_n+1)/(a_n\kappa_n)\). The extra term in the predictive scale is the future observation’s residual variation.

When little substantive prior information is available, one may seek a prior that changes slowly over the region favored by the likelihood. Such choices are often described as vague, flat, or noninformative; reference-prior constructions provide a more formal approach to related questions. These terms do not make a prior neutral in every parameterization. The simplest example is a constant density,

whose behavior depends on the parameter space and on the coordinate in which the density is declared constant. On an unbounded parameter space it need not integrate to a finite value, so posterior propriety must be checked separately.

This is an improper prior in the sense that

Improper priors are sometimes useful because they may lead to proper posterior distributions, i.e.,

However, there is another drawback of using the “uniform” prior of \(q(\theta)\propto 1\), which is that it depends on the parameterization. Suppose \(\psi\) is another parameterization of \(\theta\). Under the “uniform” prior for \(\theta\), the prior for \(\psi\) is

which is different from the “uniform” prior for \(\psi\).

An alternative choice for non-informative priors that overcomes the non-uniqueness issue is to use

where \(\mathcal I(\theta)\) is the information for \(\theta\) based on one observation. In the case that \(\theta\) is a vector

For example, suppose that, given \(\theta\), \(Y_1,\cdots, Y_n\) are a random sample from \(N(\theta,1)\). Then \(\mathcal I(\theta)=1\). So the information-based non-informative prior coincides with the “uniform” prior, and

So the posterior is proper and

As another example, suppose that, given \(\theta\), \(Y_1,\cdots, Y_n\) are a random sample from \(\mbox{Poisson}(\theta)\). Because \(\mathcal I(\theta)=\theta^{-1}\), we take the prior to be

This gives

So

The posterior is also the basis for prediction. Suppose the observations \(D=(Y_1,\ldots,Y_n)\) arise from \(p_Y(y\mid\theta)\) and the prior is \(q(\theta)\). Bayes’ theorem gives

Now consider a future observation \(Y\) that is independent of the recorded sample conditional on \(\theta\). Its distribution is obtained by averaging the sampling density over parameter uncertainty:

This is the posterior predictive distribution. It includes both variation in a new outcome at a fixed parameter and uncertainty about that parameter. The same averaging will be central to multiple imputation, where the outcome to be predicted is an unrecorded part of an existing observation.

In the example of \(N(\theta, 1)\) and \(q(\theta)\propto 1\), we had that

Since for another \(Y\) (independent of \(D\)) we have that \(Y\mid \theta\sim N(\theta, 1)\), after some calculus, we find that

In general, the posterior distribution does not have an analytic form. So the following two steps are needed in order to draw Y from its predictive distribution. For \(j=1,2,\cdots\),

The same reasoning applies to missing values, except that their conditional distribution must respect the observed part of the same record. Under MAR, distinct data and response parameters, and a factorized prior, the joint posterior can be written

Under nonignorability the response-mechanism factor remains in this expression. Missing observations are therefore latent variables with a distribution, not unknown constants to be replaced once and then treated as data.

With full data D, Bayesian inference is based on the posterior

if the posterior density \(p(\theta\mid D)\) does not have a closed-form, the Gibbs sampler can be used to draw from (\(\ref{eq:full_bayes}\)). When missing data are present, denote \(D=(D_{\mathrm{obs}}, D_{\mathrm{mis}})\). Then the Bayesian analysis is focused on \([\theta\mid D_{\mathrm{obs}}]\), the kernel of which usually does not have a closed form. Assuming data are MAR,

The first option \(q(\theta)=q(\beta\mid \alpha)q(\alpha)\) is appealing when using the “information principle” to construct the prior. Specifically, denote \(\mathcal I(\beta\mid Z)\) as the information of \(\beta\) resulting from the conditional likelihood \(p(Y\mid Z;\beta)\), that is,

Then, the marginal information of \(\beta\) is

Then, using the information principle, we may choose

Sometimes calculating \(\mathcal I(\beta\mid \alpha)\) through (\(\ref{eq:beta_info}\)) is difficult because the integral may not be analytic. In such cases, one has to construct non-informative priors by other means. Write \(D_Y=(Y_1,\cdots, Y_n)\) and \(D_Z=(Z_1,\cdots, Z_n)\). Denote

and

This gives

Thus, the joint distribution of \((D,\theta)\) is

We describe the Gibbs sampling steps for \(D_{\mathrm{mis}}, \beta,\) and \(\alpha\) and the associated kernel functions for drawing from the conditional distributions. For simplicity, the iteration index is omitted. Write \(Z_i=(Z_{i,\mathrm{obs}},Z_{i,\mathrm{mis}})\) and \(Z_{i,\mathrm{mis}}\) is missing at random.

The full conditional calculations use the specified working joint prior

This improper prior must be accompanied by a posterior propriety check. It should not be described as the general multivariate Jeffreys prior. For \((\mu,V)\) in a \(p\)-variate normal model, using the unique entries of symmetric \(V\) as coordinates, the information determinant is proportional to \(|V|^{-p}\) and the joint Jeffreys density is proportional to \(|V|^{-p/2}\). The \(|V|^{-1/2}\) expression in the lecture calculation holds for that block only when \(p=1\). Moreover, a factor depending on conditioning parameters cannot be discarded when defining a joint prior from improper conditional pieces. Specifying the working joint prior directly removes that ambiguity while preserving the following Gibbs example.

Let

where \(Z\) consists of 1 and a \(p\)-dimensional vector of continuous covariates and \(\epsilon\sim N(0,\tau^{-1})\). Assume that the last \(p\) components of \(Z\) follow multivariate normal with mean \(\mu\) and variance \(V^{-1}\). Write \(\theta=(\beta,\tau,\mu,V)\). We have seen how to compute the MLE of \(\theta\) using EM algorithm (with explicit E and M steps). Now we conduct Bayesian inference. For prior construction, we first use the information principle to construct non-informative priors for \((\beta,\tau)\) given \((\mu, V)\).

By straightforward calculation,

so that

Hence,

where we recall that \(E[Z^{\otimes 2}\mid \mu, V]\) is a \((p+1)\)-dimensional square matrix.

It remains to calculate \(\det E[Z^{\otimes 2}\mid \mu, V]\). Because we can write

It is easy to see that

To calculate its determinant, we use the general rule that

Therefore we have

where \(|A|=\det A\) for any matrix \(A\). Thus,

For the multivariate normal block, \(\det\mathcal I(\mu,V)\propto|V|^{-p}\) and its joint Jeffreys density is \(|V|^{-p/2}\). Rather than multiplying improperly normalized conditional pieces, we specify the following working joint prior for this example:

The normal model permits explicit updates of the parameter blocks. Given the current parameter values, missing covariates are drawn using the conditional multivariate normal distribution already derived for the E step. With a completed dataset in hand, the remaining updates concern \(\beta\), \(\tau\), \(\mu\), and \(V\). To derive the regression update, let \(\widehat\beta\) be the complete-data least-squares estimate and complete the square in the residual sum of squares:

The first exponential factor does not depend on the candidate regression coefficient. The second identifies the conditional normal kernel in \(\beta\) and shows how the design matrix determines its posterior precision.

Likewise,

It can be seen that

Also,

Further,

Then, we have that

All four conditional distributions of the components of the parameter have closed forms and can be directly sampled from.

Let \(A=\sum_i Z_iZ_i^T\), \(\widehat\beta=A^{-1}\sum_iZ_iY_i\), and \(\operatorname{RSS}(\beta)=\sum_i(Y_i-\beta^TZ_i)^2\). Since

the conditional distribution of \(\beta\) is \(N\{\widehat\beta,(\tau A)^{-1}\}\). Multiplying the likelihood power \(\tau^{n/2}\) by the prior power \(\tau^{(p-1)/2}\) gives the gamma shape \((n+p+1)/2\) for \(\tau\mid\beta,D\), with rate \(\operatorname{RSS}(\beta)/2\).

Similarly, completing the square in \(\sum_i(W_i-\mu)^TV(W_i-\mu)\) gives \(\mu\mid V,D\sim N_p\{\overline W,(nV)^{-1}\}\). The conditional density of \(V\) is proportional to

In the chapter’s rate-matrix Wishart convention, this is \(\operatorname{Wishart}(n+p-1,S_\mu)\). A software routine that expects a scale matrix instead requires \(S_\mu^{-1}\). Its mean is \((n+p-1)S_\mu^{-1}\). The degrees of freedom must exceed \(p-1\) and the rate matrix must be positive definite. Full-rank response designs, nonzero residual sums of squares, and sufficient covariate information are essential; the mere existence of formal conditional expressions does not establish propriety of the joint posterior with missing data.

The preceding conjugate calculations rely on the quadratic structure of the normal likelihood. In a generalized linear model, the corresponding parameter conditionals usually have no familiar closed form. For a canonical link and \(\phi=1\), for example, the response likelihood is

The nonlinear cumulant term makes it difficult to choose a prior \(q(\beta\mid\alpha)\) that yields a recognizable posterior family. An information-based prior can also be difficult to compute because it requires an expectation over the covariate distribution. These difficulties affect computation, rather than the underlying posterior construction: the likelihood and prior still determine a conditional kernel that a suitable simulation method can use.

This is because in this case

and for nonlinear link \(\mu(\cdot)\), the expectation \(E[\mathcal I(\beta|Z)|\alpha]\) typically does not have a closed form. In such cases, it is convenient to take \(q(\beta\mid \alpha)\propto 1\) if this results in proper posteriors. Then, Gibbs sampler generally needs to sweep through each univariate component of \(\beta\) and \(\alpha\). In doing so, notice that

and

Nonignorable missingness introduces another component into the model. Let \(R\) identify the observation pattern, and specify a response-mechanism model with parameter \(\psi\):

The full record is now \((Y,R)\), and its likelihood includes the probability of the recorded pattern as well as the data density:

A joint prior \(q(\theta,\psi)\) completes the Bayesian specification. Unlike the ignorable case, the response-mechanism factor generally affects the distribution used to fill in missing values. Consequently, changing that model can change the posterior for the scientific parameter, even when the observed values themselves are held fixed.

The resulting Gibbs sampler alternates among missing values, data-model parameters, and response-mechanism parameters. Their conditional kernels are

The first kernel is also the one used to simulate missing values in a nonignorable Monte Carlo E step. The other two hold the completed data fixed while updating their respective parameter blocks. At each step, only the variable being sampled is allowed to vary; the remaining arguments specify the conditioning information.

Nonignorability changes the missing-value update as well as the missingness-parameter update. A draw of \(Z_{\mathrm{mis}}\) must be weighted by how compatible the completed record is with its observed response pattern. The complete-data kernel is

For the \(\psi\) update, hold \(\theta\) and completed data fixed but treat \(\psi\) as the varying argument. A prior can regularize a poorly identified direction, but posterior concentration in such a direction may reflect prior information rather than information in the observations. Sensitivity to those assumptions is part of the analysis.

Let \(\widehat\theta^{(k)}\) be the estimate obtained from imputed dataset \(k\), and let \(\widehat V^{(k)}\) be its estimated sampling covariance, for example the inverse observed-information matrix. The pooled point estimate is the average of the completed-data estimates:

Two covariance components enter the uncertainty calculation. The average complete-data covariance measures within-imputation variation,

while the sample covariance of the estimates measures between-imputation variation:

For \(K\geq2\), Rubin’s pooled covariance estimator combines them as

The factor \(1+K^{-1}\) includes the Monte Carlo uncertainty from using a finite number of imputations. Each covariance is calculated on the same parameter scale; averaging standard errors or separate significance tests would give a different procedure.

For a scalar target with large complete-data degrees of freedom, define the relative increase in variance by

An approximate reference degrees of freedom and the resulting confidence interval are

When the complete-data sample is small, a further degrees-of-freedom adjustment is needed. For vector targets, the covariance matrices must be pooled on a common parameter scale. These distinctions explain why point estimates, standard errors, and confidence limits can require different numbers of imputations to become numerically stable.

Holding \(\widehat\theta^{init}\) fixed across imputations omits its uncertainty. As a result, the usual pooling formula (\(\ref{eq:mi_var}\)) need not estimate the variance of the resulting MI estimator correctly. There are also settings in which an initially efficient estimator produces a less efficient MI estimator after this procedure (Tsiatis, 2006; Davidian, 2017). An approximate way to restore parameter variation is to draw a new value at the start of imputation \(k\):

Here \(\widehat\Sigma^{init}\) estimates the covariance of \(\widehat\theta^{init}\). Missing values are then drawn from \([D_{\mathrm{mis}}\mid D_{\mathrm{obs}};\theta_k]\). This construction uses a large-sample normal approximation for parameter uncertainty, so its justification depends on the quality of both that approximation and the initial estimator.

Drawing a fresh parameter value addresses the missing variability in the plug-in construction, but it cannot correct a distribution centered at an inconsistent estimate. A Bayesian predictive construction instead averages over the parameter posterior:

This distribution propagates parameter uncertainty together with the conditional variability of the missing values. Under a correctly specified identifiable model, suitable prior support, and appropriate regularity, it provides a basis for proper imputation. Valid frequentist pooling also requires a suitable complete-data analysis and an appropriate relationship between the imputation and analysis models. Predictive simulation alone does not remove those requirements.

The predictive distribution contains two distinct sources of uncertainty. For a scalar missing quantity \(W\), the law of total variance gives

The first term averages the residual uncertainty at a fixed parameter; the second measures how the conditional prediction changes across the posterior. Fixing \(\theta\) at a fitted value generally loses the second component, even when random residual errors are still added. Deterministic mean imputation can lose both components and alter associations.

Initialize the parameter and alternate two conditional draws. At iteration \(j+1\), first draw a new set of missing values,

and then update the parameter using the newly completed data:

These are the same data-augmentation updates used for a fully Bayesian analysis. Their interpretation depends on which output is retained. After a burn-in of \(M\) iterations and suitable convergence checks, the parameter draws \(\theta^{(M+1)},\theta^{(M+2)},\ldots\) describe the parameter posterior. For multiple imputation, retain completed datasets instead:

The observed entries remain fixed in every dataset. The completed-data analyses and their pooling then supply the inferential output, with Monte Carlo dependence between retained draws taken into account when choosing and assessing the imputation run.

So far, the same probability model has supplied both the imputations and the completed-data analysis. To distinguish their roles, write \(L(D\mid\theta)\) for the analyst’s likelihood and \(L^{imp}(D\mid\theta^*)\) for the imputer’s likelihood. The latter determines both \([\theta^*\mid D]\) and \([D_{\mathrm{mis}}\mid D_{\mathrm{obs}};\theta^*]\). Using the same model gives

Meng’s (1994) notion of congeniality concerns whether imputation and analysis can be justified together under a compatible inferential construction. Identical likelihoods provide one natural way to arrange this compatibility, but they are not its general definition. What matters is whether the imputation preserves the distributional features and uncertainty relevant to the analysis.

Multivariate normal imputation is especially simple under a monotone missingness pattern. Let \(Y=(Y_1,\ldots,Y_p)\) and factor its distribution into successive regressions on the preceding variables. For \(j=1,\ldots,p\), take

The coefficient \(\beta_{j0}\) is an intercept, and \(\tau_j\) is the conditional precision. This representation connects a joint normal distribution to a sequence of ordinary normal regressions. The original course discussed software examples including SAS PROC MI, NORM, AMELIA, and MICE; the statistical distinction is between a joint-model construction such as this one and separately specified conditional models, which are considered below.

If we choose the improper prior \(q(\beta_1,\tau_1,\cdots, \beta_p,\tau_p)\propto \prod_{j=1}^p\tau_j^{-1}\), then, one can show

where \(n_j\) is the number of complete cases for the \(j\)th regression model, \(\widehat\beta_j\) and \(\widehat\tau_j\) are the CC (thus observed-data) MLEs, and \(S_j\) is the design matrix for the regression model (\(\ref{eq:mono_mvn}\)). Moreover, because of the factorization in the likelihood, the \((\beta_j, \tau_j)\) are independent given \(D_{\mathrm{obs}}\). So, to draw \((\beta_j, \tau_j)\) from the posterior, one first draws \(\tau_j\) from \([\tau_j\mid D_{\mathrm{obs}}]\) and then \(\beta_j\) from \([\beta_j\mid D_{\mathrm{obs}}; \tau_j]\).

Now, to complete the imputation, it remains to draw \(D_{\mathrm{mis}}\) from \([D_{\mathrm{mis}}\mid D_{\mathrm{obs}};\theta^*]\), where \(\theta^*=(\beta_1,\tau_1,\cdots, \beta_p,\tau_p)\) are freshly drawn from their posterior. Because of the special structure, this can be done progressively in \(j=1,\cdots, p\). For each subject, assume the missing values in \(Y_{1},\cdots, Y_{j-1}\) have been imputed. Then, because of the imputation model (\(\ref{eq:mono_mvn}\)), \(Y_j\) can be imputed by

where \(\epsilon_j\) is simulated from \(N(0,1)\).

Under a monotone pattern, order variables so that subjects observed on variable \(j\) also have its predecessors. Let \(n_j\) be their number, let \(X_j\) be the \(n_j\times j\) design including an intercept, and let \(y_j\) be their responses. For the factorized regression prior \(q(\beta_j,\tau_j)\propto\tau_j^{-1}\), define

Integration over the \(j\) regression coefficients gives

Equivalently, \(\tau_j=\chi^2_{n_j-j}/\operatorname{RSS}_j\). If \(\widehat\tau_j=n_j/\operatorname{RSS}_j\), the multiplier is \(\widehat\tau_j/n_j\), not \(\widehat\tau_j^2/n_j\). These draws require \(n_j>j\), full column rank, and positive residual sum of squares. Use the same sampled coefficient and precision within one imputed dataset, but independent residual draws for different subjects. Progress from earlier to later variables, using previously imputed predecessors where necessary. Draw new parameter values for the next imputed dataset.

Now we look at what precisely is done in the imputation for the missing \(Y_j\). If \(Y_j\) is continuous, then the regression model can be chosen to be a linear regression, and the imputation is similar to what’s been described in the imputation step for the monotone missing case with multivariate normal model. If \(Y_j\) is categorical, then typically a GLM is used to model it, say, logistic regression. Denote the MLE of the regression parameter as \(\widehat\beta_j\). Then, drawing from the “posterior predictive distribution” is similar to what we have mentioned as the “sort-of” proper imputation. Specifically, we first draw \(\beta_j\) from

where \(\widehat\Sigma_j\) is an estimate for the variance of \(\widehat\beta_j\).

Then, \(Y_j\) is drawn from

where \(Z_j=(Y_1,\cdots, Y_{j-1}, Y_{j+1}, Y_p)^{\mathrm{T}}\). The MICE algorithm is implemented in, for example, SAS PROC MI using the FCS option and the R package MICE. Some caveats: The MICE/FCS approach is predicated on specification of full conditional models; that is, in our example for each \(j\), models for \(Y_j\) as a function of all other variables.

For nonignorable missingness, the imputation model must incorporate additional identifying assumptions. A selection-model construction specifies \(\pi(R\mid Y;\psi)\) and combines it with the imputer’s full-data likelihood:

A joint prior for \((\theta^*,\psi)\) then supplies the posterior predictive distribution used for completion. The analyst may apply the same complete-data calculation to datasets generated under several alternative assumptions, which makes MI useful for sensitivity analysis. The resulting scientific conclusions still depend on those assumptions, however; separating imputation from analysis does not make nonignorability irrelevant to the person interpreting the results.