Seeing this post by Tim Morris on the difference between iterative imputation and Gibbs sampling reminded me of some articles that my colleagues and I have written on the topic:
[2001] Using conditional distributions for missing-data imputation. Discussion of “Conditionally specified distributions” by B. Arnold et al. Statistical Science 3, 268-269. (Andrew Gelman and T. E. Raghunathan).
[2014] Multiple imputation for continuous and categorical data: Comparing joint and conditional approaches. Political Analysis 22, 497-519. (Jonathan Kropko, Ben Goodrich, Andrew Gelman, and Jennifer Hill).
[2014] On the stationary distribution of iterative imputations. Biometrika 1, 155-173. (Jingchen Liu, Andrew Gelman, Jennifer Hill, Yu-Sung Su, and Jonathan Kropko).
It’s something that’s been discussed in the Bayesian statistics community for a long time—I remember talking about it with Jun Liu around 1992 or so.
For a very simple example, consider these two incoherent conditional specifications:
x|y ~ normal(0, 1)
y|x ~ normal(x, 1).
These are obviously incoherent: in the specification for x|y, x and y are independent; in the specification for y|x, they are dependent.
What happens if you do Gibbs updating (following Morris, I won’t call it “Gibbs sampling,” as there is no underlying joint distribution)?
Iteration 1: Draw x from normal(0, 1)
Iteration 2: Draw y from normal(x, 1)
Iteration 3: Draw x from normal(0, 1)
Iteration 4: Draw y from normal(x, 1)
Etc.
What is the joint distribution of (x,y)? It depends on when you stop. If you stop at an odd iteration, the joint distribution is (x,y) ~ MVN((0,0), ((1,0),(0,2))). If you stop at an even iteration, the joint distribution is (x,y) ~ MVN((0,0), ((1,1),(1,2))). If you want a single joint distribution, you could stop at a random iteration, in which case the joint distribution will depend on your random process, but for any reasonable choice of random stopping point will be approximately (x,y) ~ 0.5*MVN((0,0), ((1,0),(0,2))) + 0.5*MVN((0,0), ((1,1),(1,2))).
Things get more complicated when you have more than two variables, because then there are other possible updating schedules. You could update x,y,z,x,y,z,…, or x,y,x,z,y,z,…, or do a random update—and these could lead to different joint distributions, even after averaging over a random stopping point. There have to be some conditions under which you could bound the variation around these possible distributions—at the extreme case where the conditional specifications are coherent and the process is ergodic, there’s no variation at all—but I don’t know what they are.
Iterative imputation for missing data is kinda like the above, except the conditional distributions are themselves estimated from available data. In some ways this makes the problem harder to analyze (which is why the math in the third article linked above gets so complicated), but in other ways it makes things better that the distributions are estimated from a common dataset, as this should enforce some approximate coherence among the conditional specifications.
Lots to chew about, from practical, theoretical, and computational perspectives.
In that case, is it correct to say “incoherent” actually means sampling from the average of distinct “coherent” distributions?
Eg, same as the average person has about one ovary?
In general, taking some simple example and treating it as a special case seems it is a fertile strategy for finding new insights (eg, Euclids 5th postulate).
Were these methods introduced because it used to be too hard to do joint Bayesian inference for the missing data and model parameters? Nowadays, I’d prefer to actually model the missing data as additional parameters and turn the Bayesian crank. To simulate GLMs, you can take a multivariate normal model of the data with link functions. Of course, we can’t do that in Stan if there is discrete data that can’t be marginalized (of course, if it can’t be marginalized, it’s probably combinatorial and probably leads to NP-hard inference, like say, variable selection with spike-and-slab priors).
Bob:
The reason for these methods is that the only easily available methods for modeling multivariate distributions are the multivariate normal and mixtures of multivariate normals. Usual state of the art is to just use multivariate normal, but that can have problems, as discussed in my paper with Kobi from a few years back. In practice, I think iterative regression imputation can give better results than using the multivariate normal (that’s the topic of the Kropko et al. paper cited in the above post).
Fundamentally I don’t think the big problem is the discreteness of the data. Yes, in political and social surveys almost all our variables are discrete, but I think it’s usually fine to model using continuous variables and then round at the end, or model using latent continuous variables. The real issue is that there’s structure in the data that should be easy enough to pick up in a regression but would be hard to put in a big dumb multivariate normal model.
I know that some approaches have been proposed to fit large multivariate nonparametric models to do missing-data imputation, and I could imagine that working well, if used as part of a larger workflow. If there’s black-box imputation, I’d like general methods of diagnosing problems and improving the imputation models. Again, that comes up in my paper with Kobi, and we had big plans to continue in that direction–we put some of this in our “mi” package in R–but we never completely followed through with it. I think this is an important and fertile area of research . . . OK, once I’m finished typing this comment, I’ll add it to my list of unwritten research articles.
I really really don’t think the incoherence of iterative imputation is any sort of practical problem. It’s more of a theoretical oddity than anything else. On the flip side, this suggests the possibility of a theoretical result (generalizing from the Liu et al. paper cited in the above post) demonstrating that, indeed, the incoherence is something you don’t really have to worry about, that in practice you can just use a random updating scheme or whatever.
I understand that the difficulty of modeling non-parametric multivariate distributions is why people turn to conditional modeling-based imputation. But it’s not clear to me why, for conditional modeling-based imputation, people turn to iterative regression (ie MICE). As far as I can tell, the original reason for MICE is that, back then, people knew how to do regression (ie predicting the conditional mean), but didn’t know how to do conditional modeling (ie predicting the full conditional distribution of even a single feature). But if one knows how to do nonparametric conditional modeling, one can simply autoregressively sample from the conditional distributions. Thus one need not iteratively revisit previously-sampled features, as in MICE. And we now do know how to do nonparametric conditional modeling, with either diffusion modeling or hierarchical partitioning (see my recent preprint: https://arxiv.org/abs/2407.05593). So it’s not clear to me that MICE, with its attendant incoherence problem, is really needed.
Calvin:
“MICE” is a computer program that stands for Multiple Imputation using Chained Equations. I don’t like this term because the imputation doesn’t use any equations. It uses regression models.
As discussed in the above-linked paper by Kropko et al., there are lots of different ways to impute missing data, and different methods correspond to different models. Iterative imputation implicitly corresponds to a model, too–it’s just that the implicit joint distribution of the imputations is complicated and not consistent with the conditionals being fit.
I find this normal example interesting. When teaching Gibbs, people often justify why it is enough to specify the joint purely by the conditionals. There is a Hammersley–Clifford theorem seemingly suggests that you can get the joint from the conditional:
f(x,y) = \frac{f_{Y|X}(y \mid x)}{\int \left[ f_{Y|X}(y \mid x) / f_{X|Y}(x \mid y) \right] dy}.
Using this formula, the incoherent normal example leads to f(x, y)= normal (x | 0,1) normal(y| x, 1), which surprising is the even iteration result. Why does it only matches the even iteration? I don’t know.
(obviously you can change x and y in this formula, but the integral diverges. Is at least one integral diverges = incoherence? I do not know.)
So I guess the Hammersley–Clifford theorem can only be applied to coherent conditionals, which (to me) is not well-known.
In additional to MI, in likelihood-free computation, sometimes we need to estimate p(y|theta) and p(\theta|y) together. In this case both theta and y can be very high dimensional objects. Is the incoherence prevalent in this regressions? If there is incoherence, is there some way to post-hoc adjust the estimated conditionals? … many questions
Yuling:
Yes, the Hammersley-Clifford theorem assumes as a premise that the conditional distributions are compatible with each other.