Best way to combine MCMC inference with multiple imputation? I can derive an MCMC algorithm for sampling from the posterior distribution of a parameter vector of interest, but only starting with a dataset that has no missing values. The actual dataset that I want to use for inference has substantial missingness in its covariates.
One approach would be to build a more complex MCMC algorithm that, for example, first fills in the missing data with draws from the missing values' posterior predictive distribution. However, this feels intractably hard.
What I'd rather do is use an off-the-shelf method to generate multiple imputations of the dataset (such as a MICE package), then run my existing MCMC algorithm on each imputed, complete dataset and then recombine into final estimates of (for example) a posterior expectation or a posterior interval for a parameter of interest.
Is there a body of literature that attempts to solve problems in this way? Or is there a much better way to do this? Or is this approach wrong-headed or infeasible? Any pointers would be helpful.
 A: One well-known approach is exactly what you describe (if I understood correctly): i.e. combine the inferences from the analyses of a large number of imputed datasets (each analyzed separately) by just using the equally weighted mixture distribution of the posterior distributions as the combined posterior distribution. In practice, this means that we just throw together all the (or a random subset) of the MCMC samples that we generated by analyzing each of the multiply imputed datasets on their own. We use this combined set of MCMC samples from all the models, as if they were a giant single set of MCMC samples. This is described in Bayesian Data Analysis (3rd edition by Gelman et al. 2014 on page 452). It's also investigated and reported to perform well with a sufficiently large number of imputations (e.g. 100) by Zhou and Reiter 2010 (Zhou, X. and Reiter, J.P. 2010. A note on Bayesian inference after multiple imputation. The American Statistician, 64(2), pp. 159-163.).
Note: checks for non-convergence need to be done within each imputed datasets. E.g. if you treat chains from separate imputations as if they were separate chains run on the same datasets, checks like Rhat etc. will tell you that you have a problem, because - unsurprisingly - differently imputed datasets do result in (slightly) different posterior distributions. With that caveat, I've essentially treated them like separate chains (e.g. if I had 4 chains per imputation and 100 imputation, I thereafter proceeded as if I had 400 chains).
There's alternatives to this. E.g. it is also clear that you could also run a single MCMC sampler across all imputations and at each MCMC iteration use the "mixture-likelihood" for all the imputations, but I'm not aware that any software does this easily (you can of course hand-code this in something like Stan). The alternative to have a single model that imputes and analyzes at the same time is of course also a serious option, but again harder to implement. Keeping imputation and analysis as separate steps also has advantages in terms of addressing different questions/estimands (e.g. jump-to-reference-group imputation becomes easier to implement).
A: (Thanks for the reminder! The answer has been revised to focus more on the originally posted question)
I think the paper by Zhou, as recommended above by @Björn, specifically discussed about Bayesian inference after multiple imputation.
Recently, I am also applying Bayesian conditional logistic regression after multiple imputation (m = 100) in my PhD project. I mixed the MCMC draws from each imputed dataset, and used them to approximate posterior distribution and calculate point estimate and credible interval. I also combined the MCMC draws within each imputed dataset and plotted the posterior distribution across all the datasets in the same graph, as a way to present the variation in posterior approximation across imputed datasets.
@Björn suggested above that treating chains from imputed datasets as separate chains and combine them across datasets. It seems reasonable but also I am a bit unsure. I am wondering if there is any reference supporting this. With this in doubt, as for checking MCMC sampling quality, I checked only for the first imputed dataset (or you could check any other one or multiple datasets) – calculating Rhat and effective sample size, and plotting a trace plot (and there are many other useful metrics and plots). Sorry I am new to the forum so I couldn't reply directly above to @Björn's response.
A: Are you familiar with censorship in survival analysis modeling? Censorship means that a patient was removed from the study without explicitly dying (at least in the context of survival analysis.) There are appropriate architectures to handle censorship; however, it's more a question of model architecture than the underlying MCMC sampling algorithm. For exploration of survival analysis via PyMC3 (uses HMC), see here.
In your write-up, it's not clear whether you have a survival analysis problem or something different. Nonetheless, I'd recommend approaching your problem from a model architecture perspective not an MCMC-variant-specific perspective. If you could provide some more details, we (the community) could provide some recommendations on how to handle observed and latent values (censoring.)
You didn't explicitly say that you're approaching this from a Bayesian perspective, however, MCMC is very common in the Bayesian paradigm, so I'll assume you prefer to be as Bayesian as possible. With that said, imputing values is not very Bayesian at all, as you'd be adding biases to your dataset without properly defining them via prior distributions. Censorship is the way to go so you have a "fork" so to speak that distinguishes the observed from the latent, and this will be captured nicely by your posterior distribution.
A: I wrote a paper arxiv.org/abs/1907.09090 that describes how the pseudo-marginal approach can impute missing data. 400 covariates sounds tough, though, to be completely honest. Depends on what kind of distributions you want to put on the columns, the number of rows, how you program everything. Intractable in your case? Probably, yes.
In section 3.3, we describe the approach (not ours) that sounds like what you want. Maybe some references in there will give you some formulas. Here’s a quote:

Multiple Imputation (MI) generates multiple complete data sets by
sampling several sets of plausible values for each missing data point
by sampling from the posterior predictive distribution [19], [20],
[7]. The same analysis is performed separately on each data set, and
the results are then combined. For example, in the context of
regression analysis, the model parameters derived from each imputed
dataset [sic] are combined by a simple average. The parameter
variances are calculated by averaging the individual variances from
each imputation, and the formula includes an additional term to
capture the between-imputation variance.

