Why several incomplete predictors require a joint audit
When two predictors are incomplete, chained equations alternate between them: a model imputes one predictor using the other and the outcome, then another model imputes the second predictor. It is tempting to check each imputation equation separately for familiar terms. But an equation can look reasonable on its own while the set of equations does not correspond to any joint distribution that preserves the planned analysis model.
This issue is especially consequential when the analysis regression contains both an interaction and a nonlinear effect. Here we use two incomplete quantitative predictors, \(X\) and \(Z\), and an outcome \(Y\). The substantive regression includes \(XZ\) and \(X^2\). The key question is not simply whether these terms appear somewhere in the imputation equations. It is whether the conditional distributions used for imputation can be the conditionals of a joint model whose outcome regression is the intended one.
Compatibility is a mathematical property of the models, not a guarantee supplied by chained-equation software. It also differs from the missingness assumption: a missing-at-random assumption, when appropriate, concerns the mechanism producing missing values; it does not make incompatible conditional models compatible.
A joint model that preserves the analysis regression
Consider this analysis model:
To derive full conditionals for imputing \(X\) and \(Z\), the outcome regression alone is not enough. We also need a model for the predictors. Suppose \((X,Z)\) has a bivariate normal distribution. Then \(X\mid Z\) and \(Z\mid X\) are normal, with linear conditional means. Combining this predictor model with the displayed outcome model defines a joint distribution for \((Y,X,Z)\). Its conditionals provide a precise reference for auditing the chained equations.
Write the predictor conditionals as \(X\mid Z=z\sim N(m_X(z),\tau_X^2)\) and \(Z\mid X=x\sim N(m_Z(x),\tau_Z^2)\). These two normal conditionals are not arbitrary: they must come from the same bivariate normal distribution. The variances and slopes therefore obey restrictions implied by a shared covariance matrix.
For fixed \(x\), collect the terms in the outcome mean that involve \(z\):
The likelihood contribution for \(z\), given \(x\) and \(y\), is normal in \(z\). Multiplying it by the normal predictor conditional \(Z\mid X=x\) gives another normal kernel. Completing the square yields:
Thus, even though this full conditional is normal in \(z\), its mean and variance generally vary with \(x\) through nonlinear expressions. In particular, the denominator contains \(\{\beta_2+\beta_3x\}^2\), and the numerator involves \(x^2\) through \(A(x)\). A simple linear regression of \(Z\) on \(Y\) and \(X\) is not generally this conditional distribution.
For \(X\), fix \(z\) and define \(C(z)=\beta_0+\beta_2z\) and \(D(z)=\beta_1+\beta_3z\). The joint-model full conditional has kernel:
If \(\beta_4\ne0\), squaring the quadratic expression in \(x\) produces a fourth-degree term in the log density. That conditional is generally not normal. Consequently, a normal linear regression for \(X\) may fail to reproduce the full conditional even if it includes \(Y\), \(Z\), and their interaction.
Worked audits
Worked Example: Deriving the Conditional for \(Z\)
Suppose the hypothetical analysis model has \(\beta_0=0.5\), \(\beta_1=0.8\), \(\beta_2=-0.4\), \(\beta_3=0.3\), \(\beta_4=0.2\), and \(\sigma^2=1.44\). Let the hypothetical predictor distribution be bivariate normal with \(\operatorname{Var}(X)=4\), \(\operatorname{Var}(Z)=1\), and \(\operatorname{Cov}(X,Z)=1\), both means zero. Its conditional distributions are \(X\mid Z\sim N(Z,3)\) and \(Z\mid X\sim N(0.25X,0.75)\). For \(x=1\) and \(y=3\), calculate the full conditional for \(Z\).
First, \(A(1)=0.5+0.8(1)+0.2(1)^2=1.5\), and \(B(1)=-0.4+0.3(1)=-0.1\). The prior conditional mean is \(m_Z(1)=0.25\). The conditional precision is:
The numerator for the mean is:
Therefore the conditional mean is \(0.229167/1.340278\approx0.1710\), and the conditional variance is \(1/1.340278\approx0.7461\). Under this specified joint model, \(Z\mid X=1,Y=3\) is normal with mean about 0.1710 and variance about 0.7461. The result depends on the interaction and outcome residual as well as the predictor distribution.
Worked Example: Checking the Shape of the Conditional for \(X\)
Use the same hypothetical model and predictor distribution, and consider \(z=0\), \(y=2\). Since \(X\mid Z=0\sim N(0,3)\), its prior conditional mean is zero and its variance is 3. The outcome mean at these fixed values is \(\eta=0.5+0.8x+0.2x^2\). The likelihood residual is \(2-\eta=1.5-0.8x-0.2x^2\), so the log-kernel, apart from a constant, is:
Expanding the square gives \(2.25-2.4x+0.04x^2+0.32x^3+0.04x^4\). The log-kernel therefore contains cubic and quartic terms; it is not the quadratic log-kernel of a normal distribution. For a numerical check, at \(x=0\) the log-kernel is \(-2.25/2.88=-0.78125\). At \(x=1\), the residual is \(0.5\), so the log-kernel is \(-0.25/2.88-1/6\approx-0.2535\). These calculations illustrate the shape implied by the model; they do not turn the conditional into a normal distribution. A normal linear imputation for \(X\) is therefore not the exact conditional under this joint model.
Worked Example: Auditing a Pair of Proposed Linear Imputations
Suppose an analyst proposes \(X\mid Y,Z\) as a normal linear regression with predictors \(Y\), \(Z\), and \(YZ\), and proposes \(Z\mid Y,X\) as a normal linear regression with predictors \(Y\), \(X\), \(X^2\), and \(YX\). These equations include recognizable terms related to the analysis model. Are they, by that fact alone, compatible with the joint model just specified?
No. The exact \(X\mid Y,Z\) kernel in the preceding example has cubic and quartic powers of \(X\) in its log density, whereas a normal linear regression conditional for \(X\) has a quadratic log density in \(X\). The proposed equation cannot equal that exact conditional when \(\beta_4=0.2\).
The proposed equation for \(Z\) also does not establish compatibility. Under the specified joint model, its exact conditional mean has denominator \(1/\tau_Z^2+\{\beta_2+\beta_3X\}^2/\sigma^2\), and its variance is the reciprocal of that expression. A normal regression with a linear predictor generally has constant residual variance and a mean linear in the listed predictors, so it will not generally match this conditional either. Including \(X^2\) and \(YX\) may enrich the mean but does not supply the required \(X\)-dependent variance and denominator. The audit must compare full conditional distributions, not just lists of terms.
What a compatibility audit can establish
The derivation gives a practical sequence for auditing multiple incomplete predictors. First state the intended outcome model. Next specify, at least conceptually, a joint predictor distribution. Then derive the full conditional for each incomplete predictor, conditioning on the outcome and all other predictors. Finally compare the proposed chained-equation model with those derived conditionals, checking both distributional shape and parameter dependence.
Make the interaction, nonlinear terms, and their coefficients explicit before evaluating any imputation equation.
The outcome regression alone does not determine the distributions of incomplete predictors given the outcome.
Check its shape, conditional mean, conditional variance, and dependence on the other predictor and outcome.
Matching some terms in a regression formula is not sufficient if the implied conditional variance or functional form differs.
If the proposed conditionals equal the conditionals derived from a single joint distribution, they are compatible with that joint model. If an exact equality fails, the proposed set is not compatible with that particular joint model. This is a model-specific conclusion: different predictor distributions could imply different full conditionals, and compatibility with some other joint distribution is a separate question.
In practice, exact full conditionals may be inconvenient. A richer approximate imputation model can preserve important features of the analysis relationship, but term inclusion alone does not prove exact compatibility. Joint-model methods or substantive-model-compatible methods offer alternatives when the target regression structure makes ordinary conditional regressions inadequate. These alternatives also require defensible distributional assumptions and appropriate implementation; they do not eliminate the need to inspect the model.
Assumptions, diagnostics, and limitations
The calculations above rely on a correctly specified Gaussian outcome model with constant variance and a bivariate normal predictor distribution. They also assume the displayed coefficients and variance describe the target analysis relationship. If the predictors are categorical, heavily skewed, bounded, or have a different joint structure, these specific conditional formulas need not apply. A compatibility audit should use a predictor model appropriate to the variables at hand.
A useful diagnostic is to look for structural mismatches: a normal conditional where the derived log density is quartic; a constant residual variance where the exact conditional variance depends on another predictor; or a linear conditional mean where the derived mean is nonlinear or rational. These are warning signs of incompatibility with the proposed joint model. Numerical instability, implausible imputed values, or marked changes in important fitted relationships across imputed data sets can also prompt review, but none alone proves incompatibility. They are not substitutes for deriving or otherwise assessing the conditional models.
- Do not infer compatibility from a shared term list. Interactions and squared terms can induce nonstandard conditional distributions and variances.
- Do not infer incompatibility with every possible joint model from one failed comparison. A failed equality rules out the specified joint model, not all possible joint distributions.
- Do not confuse compatibility with correctness. A set of conditionals can be compatible but still describe the data poorly or fail to represent the intended analysis relationship.
No single simple joint model is guaranteed to represent every realistic predictor distribution. With several incomplete covariates, the number and complexity of implied conditionals increase, and different modeling choices can produce different conclusions. Compatibility is therefore a valuable audit of internal coherence, not a universal guarantee of unbiased estimates, adequate uncertainty, or good imputations.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Why is the outcome regression alone insufficient to derive \(X\mid Y,Z\) and \(Z\mid Y,X\)?
- For the model in this tutorial, what feature makes the conditional log density for \(X\) generally non-quadratic when \(\beta_4\ne0\)?
- Does including \(X^2\) in a linear regression for \(Z\) establish compatibility? Explain.
- What does a failed comparison between proposed conditionals and a specified joint model establish?
Solutions. (1) The outcome regression specifies \(Y\mid X,Z\), not the joint distribution of \(X,Z\); a predictor distribution is also needed to derive the other full conditionals. (2) Squaring the quadratic expression in \(X\) creates a fourth-degree term in the log density, and generally also a cubic term. (3) No. The exact conditional mean and variance can have nonlinear, rational dependence on \(X\), and compatibility concerns the full distribution, not just one term. (4) It rules out compatibility with that specified joint model; it does not rule out compatibility with every possible joint distribution.