Diagnosing Multiple Imputation and Model Congeniality

Why Multiple Imputation Needs Diagnostics

Multiple imputation produces several completed data sets, but producing them is not evidence that the imputation procedure is adequate. The imputed values must be plausible given the observed information, the computation must be stable enough for its output to be useful, and the imputation models must support the analysis that will be applied to the completed data. These are related questions, but they are not interchangeable.

A useful diagnosis therefore has several parts. Compare imputations with observed data in settings where that comparison is meaningful; inspect chains and summaries for computational problems; and check whether the imputation model preserves the relationships required by the substantive model. None of these checks can verify assumptions about values that are never observed. In particular, an attractive diagnostic display does not establish that missingness is ignorable.

Definition: An imputation-model diagnostic is a check of whether the fitted imputation procedure produces values and estimates that are plausible in light of observed data and its stated model. A computational diagnostic checks whether the algorithm has stabilized and explored its intended distribution. Model compatibility means that the conditional distributions used for imputation can arise from a common joint distribution. Congeniality is a stronger relationship between the imputation procedure and the substantive analysis: broadly, the imputation model must preserve the features needed by the analysis, and the analysis must not demand information or structure that the imputation procedure has excluded.

Compatibility and congeniality are not synonyms. Compatibility is a mathematical property of conditional models. Congeniality also concerns the relationship between the information and assumptions available to the imputer and those used by the analyst. In practice, a key concern is whether the imputation models retain the outcome-predictor relationships, nonlinear terms, interactions, and design features that the substantive analysis will use.

Observed-Data and Masked-Value Checks

Start by comparing imputed values with observed values for the same variable, but do not expect their overall distributions to be identical automatically. Observed values and missing values may belong to different covariate patterns. If missingness depends on observed predictors, the imputed values should be judged conditionally on those predictors, not only by an overall histogram.

For a variable $X$, useful comparisons may include its distribution among observed and imputed cases within categories of an observed predictor, or its relationship with an outcome among those cases. Check plausible ranges, means, quantiles, tail behavior, and associations. For binary or categorical variables, compare category frequencies within relevant strata. A close marginal match can conceal a poor conditional match; a marginal difference can also be appropriate if the cases receiving imputations have different observed covariates.

When the observed data are sufficiently informative, an especially direct check is to hide some observed values temporarily. Fit the imputation procedure using the remaining data, predict the held-out values, and compare imputations with the values that were deliberately hidden. Repeating this process across several splits gives a measure of predictive performance in observed cases. For continuous variables, assess bias, prediction error, and interval coverage; for categorical variables, assess calibration and classification probabilities as appropriate. This evaluates prediction where values are known, not the truth of predictions for values missing in the original data.

Worked Example: A Conditional Distribution Check

Suppose a continuous biomarker $X$ is observed for some patients and imputed for others. Patients are grouped by a fully observed disease-severity indicator $Z$. Among patients with $Z=0$, observed $X$ has mean 10 and standard deviation 2, while the imputed values have mean 10.4 and standard deviation 2.1. Among patients with $Z=1$, observed $X$ has mean 20 and standard deviation 3, while imputed values have mean 19.9 and standard deviation 3.1.

The pooled observed mean is not enough to assess fit: if the observed and imputed groups have different proportions with $Z=1$, their overall means may differ even when each within-group comparison is close. The relevant comparisons are $10.4-10=0.4$ for $Z=0$ and $19.9-20=-0.1$ for $Z=1$. These differences are small relative to the within-group standard deviations, so the summaries do not reveal a large location or scale discrepancy.

This is only a preliminary check. We would also inspect quantiles and tails, and check whether the imputation procedure reproduces the association of $X$ with other observed predictors. Similar means and standard deviations do not rule out skewness, multimodality, or a misspecified conditional relationship.

Result: Suppose a validation value $X$ is observed, and a validation procedure hides it using a mechanism that is independent of $X$ conditional on the predictors supplied to the imputation model. If the validation cases are representative of the cases to which that model will be applied, performance on the hidden values estimates predictive performance for that same observed-data population. It does not, without further assumptions, establish predictive performance for the originally missing values.

The distinction follows from what is being sampled. Artificially hidden values come from the distribution of observed $X$ values under the validation design. The original missing values may have a different distribution if missingness depends on $X$ even after conditioning on the predictors. Validation can therefore expose poor prediction and calibration, but it cannot identify a distribution for values that are systematically unobserved without assumptions about missingness.

Worked Example: Interpreting Masked-Value Validation

Suppose 200 observed biomarker values are randomly divided into training and validation sets. The fitted imputation model predicts the 50 held-out values with average prediction error $-0.2$ units and root mean squared error 3 units. Its nominal 90% predictive intervals contain 44 of the 50 held-out values.

The empirical coverage is $44/50=0.88$, close to the nominal 0.90. The average error is small compared with the root mean squared error, so these results do not suggest a large average bias in this validation split. They do not prove calibration: with only 50 cases, coverage is imprecisely estimated, and the validation cases may not represent the cases originally missing.

A useful next check is to repeat the split, examine error and coverage across important predictor strata, and check whether the model predicts the tails adequately. Reporting only the overall root mean squared error could conceal poor prediction in a clinically important subgroup.

Trace and Distribution Checks

When imputation is iterative, a trace plot can reveal behavior that a final completed data set cannot. Plot chain-specific summaries across iterations: for example, the mean and standard deviation of imputed values, selected regression coefficients, or a proportion of imputations in a scientifically important range. Chains that drift, jump between distinct regimes, or remain separated may indicate incomplete burn-in, poor mixing, unstable fitting, or a problematic imputation model.

No single trace proves convergence. A stable-looking plot can miss slow movement, and a flat trace does not show that the chain has explored all relevant regions. Examine multiple summaries and, when available, multiple independently initialized chains. Autocorrelation matters because a long sequence of highly dependent draws contains less information than the same number of nearly independent draws.

Theorem: Let $T_1,\ldots,T_n$ be a segment of a weakly stationary chain with finite variance $\gamma_0=\operatorname{Var}(T_t)$, and let $\rho_k=\operatorname{Corr}(T_t,T_{t+k})$. For the chain average $\bar T=n^{-1}\sum_{t=1}^n T_t$, $$ \operatorname{Var}(\bar T)=\frac{\gamma_0}{n}\left[1+2\sum_{k=1}^{n-1}\left(1-\frac{k}{n}\right)\rho_k\right]. $$ Thus, positive serial correlation can make the Monte Carlo variation of the average larger than $\gamma_0/n$.

Proof. Stationarity gives $\operatorname{Cov}(T_i,T_j)=\gamma_{|i-j|}$, where $\gamma_k=\gamma_0\rho_k$. Expanding the variance of a sum,

\operatorname{Var}(\bar T) =\frac{1}{n^2}\left(n\gamma_0+2\sum_{k=1}^{n-1}(n-k)\gamma_k\right) =\frac{\gamma_0}{n}\left[1+2\sum_{k=1}^{n-1}\left(1-\frac{k}{n}\right)\rho_k\right].

For each lag $k$, there are $n-k$ pairs $(i,j)$ with $j-i=k$, which accounts for the factor $n-k$. If the relevant autocorrelations are positive, the terms add to the independent-draw variance rather than canceling. This identity describes Monte Carlo precision under stationarity; it does not establish that an imputation algorithm has reached stationarity. $\square$

Worked Example: Autocorrelation and Monte Carlo Precision

Suppose a chain summary has marginal variance $\gamma_0=4$, and, for illustration, the lag-one correlation is $0.5$ while higher-lag correlations are zero. With $n=100$, the variance formula gives

\operatorname{Var}(\bar T) =\frac{4}{100}\left[1+2\left(1-\frac{1}{100}\right)(0.5)\right] =0.04(1.99)=0.0796.

The corresponding standard error is $\sqrt{0.0796}\approx 0.282$. If the 100 values had been independent with variance 4, the variance of their average would be $4/100=0.04$, with standard error $0.2$. The dependent chain therefore gives a less precise estimate of its mean, despite having the same number of recorded iterations. This calculation is illustrative: actual chains can have several nonzero autocorrelations.

Compatibility With the Substantive Model

A substantive analysis may rely on a particular conditional relationship that an imputation model can weaken or remove. For example, if the analysis includes an interaction or a quadratic predictor, the imputation model should preserve the information needed for that term. For an outcome model with a nonlinear predictor effect, imputing the predictor from a model that assumes only a linear relationship with the outcome can distort the fitted substantive relationship.

Compatibility can be made precise for simple conditional models. If two specified conditionals cannot be generated by any joint distribution, then repeatedly drawing from them does not amount to sampling from the conditionals of a single coherent joint model. The following Gaussian case gives a direct check.

Theorem: Suppose $X\mid Y=y\sim N(a+by,s_x^2)$ and $Y\mid X=x\sim N(c+dx,s_y^2)$, where $s_x^2>0$ and $s_y^2>0$. These conditionals arise from a nonsingular bivariate normal joint distribution if and only if $$ \frac{b}{s_x^2}=\frac{d}{s_y^2} \quad\text{and}\quad \frac{1}{s_x^2s_y^2}-\left(\frac{b}{s_x^2}\right)^2>0. $$

Proof. Write the precision matrix of a nonsingular bivariate normal joint distribution as $Q$, and its linear natural-parameter vector as $h$. From the conditional distribution of $X\mid Y$, the coefficient of $x^2$ in the joint log density is $-1/(2s_x^2)$, and the coefficient of $xy$ is $b/s_x^2$. Therefore $Q_{11}=1/s_x^2$ and $Q_{12}=-b/s_x^2$. The conditional of $Y\mid X$ similarly requires $Q_{22}=1/s_y^2$ and $Q_{12}=-d/s_y^2$. A single joint precision matrix exists only if $b/s_x^2=d/s_y^2$.

A symmetric precision matrix is positive definite exactly when its first diagonal entry and determinant are positive. Here the first entry is $1/s_x^2>0$, and its determinant is $1/(s_x^2s_y^2)-(b/s_x^2)^2$. This proves necessity. Conversely, if both stated conditions hold, form $Q$ with these diagonal entries and common off-diagonal entry $-b/s_x^2=-d/s_y^2$. It is positive definite by the determinant condition, so it defines a nonsingular bivariate normal joint distribution. Choosing $h_1=a/s_x^2$ and $h_2=c/s_y^2$ gives conditional means $a+by$ and $c+dx$, respectively, with the stated conditional variances. Thus the conditionals are those of that joint distribution. $\square$

Worked Example: Checking Two Gaussian Conditionals

Consider $X\mid Y=y\sim N(2+0.5y,4)$ and $Y\mid X=x\sim N(1+0.25x,1)$. The variance convention is that the second argument of $N$ is the variance. The cross-precision terms agree:

\frac{b}{s_x^2}=\frac{0.5}{4}=0.125, \qquad \frac{d}{s_y^2}=\frac{0.25}{1}=0.25.

They do not agree, so these two conditionals cannot both be conditionals of the same joint Gaussian distribution. The issue is not that either conditional is individually an invalid normal distribution; each is well-defined on its own. The problem is their mutual incompatibility.

For comparison, change the second slope to $d=0.125$. The cross-precision terms then both equal $0.125$, and the determinant is

\frac{1}{4\cdot1}-(0.125)^2=0.25-0.015625=0.234375>0.

The conditionals are now compatible with a nonsingular joint Gaussian distribution. This check is specific to the stated linear Gaussian conditionals; it is not a general diagnostic for every chained-equations imputation procedure.

Practical Checks and Common Pitfalls

Compatibility alone does not guarantee congeniality with the analysis. In a practical review, write down the substantive model first, including its outcome, predictors, transformations, interactions, and important design variables. Then check whether each imputation model retains information needed for those features. For example, if the analysis includes $X^2$, an imputation model that only represents a linear association with $X$ may be inadequate. If the analysis estimates an interaction, imputation that omits the interaction can attenuate or otherwise distort it.

A common mistake is to treat a close match between the overall observed and imputed distributions as proof of a good imputation model. That reasoning ignores conditional structure: two groups can have identical pooled means but very different relationships with the outcome. The correction is to compare distributions and associations within meaningful observed-data strata, and to use masked-value validation where feasible.

Another mistake is to treat a stable trace as proof of valid imputations. Stability concerns computation, not whether the assumed conditional model is appropriate or whether missingness assumptions hold. Conversely, an unusual imputed distribution does not by itself prove model failure: the imputed cases may differ from observed cases because their observed predictors differ. Interpret diagnostics against the model and the missingness structure, not as isolated pass/fail tests.

Finally, no diagnostic based only on observed data can determine the distribution of values that are missing under a process depending on those unobserved values. Observed-data comparisons, masked validation, trace plots, and compatibility checks can reveal shortcomings and support transparent justification. They cannot establish untestable assumptions about the missing values.

Key takeaway: Diagnose multiple imputation on three distinct levels: predictive fit against observed information, computational stability of the imputation procedure, and compatibility with the substantive analysis. Favorable checks are evidence, not proof, and their limits should be stated.

Check Your Understanding

Use these questions to test whether you can distinguish the main diagnostic targets.

Why might the overall mean of imputed values differ from the overall mean of observed values even when within-stratum comparisons look reasonable?
What does randomly masking observed values assess, and why does it not generally verify predictions for values missing in the original data?
For a stationary chain, how does positive autocorrelation affect the variance of the chain average? What does that calculation not establish?
In the Gaussian conditional example, which equality checks whether the two conditional slopes imply the same cross-precision term?
Name two substantive-model features that an imputation procedure may need to preserve, and explain why omitting them can matter.
Why can neither a stable trace nor a close marginal distribution match establish that assumptions about the unobserved values are correct?