Imputing Predictors for a Quadratic Regression

Why a quadratic analysis model changes imputation

Suppose an analysis models a continuous outcome as a quadratic function of a predictor. The predictor's linear and squared terms together describe a curved relationship: the outcome may increase and then decrease, or decrease and then increase. If the predictor is incomplete, its imputation needs to preserve the information in both terms. A routine linear regression of the predictor on the outcome does not generally do that.

The key issue is that the analysis model describes the outcome given the predictor, while imputation needs the predictor given the observed outcome and other variables. Reversing a conditional relationship is not a matter of simply swapping the regression variables. The conditional distribution must follow from a model for the predictor as well as the stated outcome model.

Definition: Let \(X\) be an incomplete continuous predictor, \(Y\) a continuous outcome, and \(Z\) a vector of other observed predictors. A quadratic analysis model has conditional mean \(E(Y\mid X,Z)=\beta_0+\beta_1X+\beta_2X^2+\boldsymbol{\gamma}^{T}Z\). An imputation model for \(X\) is compatible with this analysis model when its conditional distribution for \(X\mid Y,Z\) follows from a joint model that includes this outcome relationship.

For a worked derivation, assume the following two models:

X\mid Z \sim N\bigl(\mu(Z),\tau^2\bigr), \qquad Y\mid X,Z \sim N\bigl(m(X,Z),\sigma^2\bigr),

where \(\mu(Z)=\alpha_0+\boldsymbol{\alpha}^{T}Z\), \(m(X,Z)=\beta_0+\beta_1X+\beta_2X^2+\boldsymbol{\gamma}^{T}Z\), and both variances are positive. These are modeling assumptions, not consequences of the quadratic regression alone. In particular, the normal model for \(X\mid Z\) must be reasonable enough for the application; a different predictor distribution would produce a different conditional distribution for imputation.

Deriving the predictor's conditional distribution

For given parameter values, the joint density factors as \(f(x,y\mid z)=f(y\mid x,z)f(x\mid z)\). Dividing by the marginal density of \(Y\mid Z\), which is constant with respect to \(x\), gives the conditional density needed for imputation:

f(x\mid y,z) \propto \exp\left\{ -\frac{\left[y-m(x,z)\right]^2}{2\sigma^2} -\frac{\left[x-\mu(z)\right]^2}{2\tau^2} \right\}.

The first term rewards predictor values whose quadratic-model predictions are close to the observed outcome. The second term weights those values according to the predictor distribution given \(Z\). Both terms matter: the analysis model alone does not specify how plausible different values of \(X\) are before observing \(Y\).

This also explains why the resulting distribution is generally not normal. When \(\beta_2\ne 0\), squaring the quadratic expression \(y-m(x,z)\) creates a fourth-degree polynomial in \(x\) in the log density. Its leading term is negative, \(-\beta_2^2x^4/(2\sigma^2)\), so the density is integrable under the stated model. But its shape need not be a single bell curve: there can be more than one region of high density. If \(\beta_2=0\), the exponent is quadratic in \(x\), and with the stated normal models the conditional distribution is normal. That special case does not justify using a normal conditional model when the quadratic term is present.

Key result: Under a normal model for \(X\mid Z\) and a normal-error quadratic model for \(Y\mid X,Z\), the exact conditional density of \(X\mid Y,Z\) is proportional to the product of the outcome likelihood and the predictor density shown above. It is generally nonnormal when \(\beta_2\ne 0\), and may be multimodal.

Worked Example: Evaluating the Conditional Density at Candidate Values

Suppose a hypothetical study uses \(X\mid Z\sim N(1.5+0.4Z,1)\) and \(Y\mid X,Z\sim N(10+2X-0.5X^2+Z,4)\). For a participant with \(Z=1\) and observed \(Y=12\), compare the conditional density at \(x=1\) and \(x=2\).

The predictor mean is \(\mu(1)=1.5+0.4(1)=1.9\). At \(x=1\), the predicted outcome is \(10+2(1)-0.5(1)^2+1=12.5\). Ignoring the same normalizing constant at both candidate values, the log density kernel is

-\frac{(12-12.5)^2}{2(4)}-\frac{(1-1.9)^2}{2(1)} =-0.03125-0.405=-0.43625.

At \(x=2\), the prediction is \(10+2(2)-0.5(2)^2+1=13\), giving

-\frac{(12-13)^2}{2(4)}-\frac{(2-1.9)^2}{2(1)} =-0.125-0.005=-0.130.

The density kernel at \(x=1\), divided by the kernel at \(x=2\), is \(\exp(-0.43625+0.130)=\exp(-0.30625)\approx0.7362\). Thus the density at \(x=1\) is about 0.736 times the density at \(x=2\), with the model parameters fixed. These are density values, not probabilities that \(X\) equals either exact value; a continuous variable has probability zero at any single point.

Why the conditional distribution can have two modes

A quadratic mean can assign similar outcome means to predictor values on opposite sides of its turning point. If the predictor distribution does not strongly favor one side, observing the outcome may leave two separated plausible regions for \(X\). A single normal imputation model may then place too much probability in the less plausible values between those regions.

Worked Example: Locating Two Modes

Consider the hypothetical models \(X\sim N(0,4)\) and \(Y\mid X\sim N(5+X^2,1)\), with observed \(Y=9\). Apart from a constant, the log conditional density is

\ell(x)=-\frac{(9-5-x^2)^2}{2}-\frac{x^2}{8} =-\frac{(4-x^2)^2}{2}-\frac{x^2}{8}.

Differentiate to find stationary points:

\ell'(x)=2x(4-x^2)-\frac{x}{4} =x(7.75-2x^2).

Thus \(x=0\) or \(x^2=3.875\), giving \(x\approx-1.969\) and \(x\approx1.969\). At zero, \(\ell''(0)=7.75>0\), so zero is a local minimum of the log density. At either nonzero stationary point, \(\ell''(x)=7.75-6x^2=7.75-23.25=-15.5<0\), so each is a local maximum. The conditional density therefore has two modes, near \(-1.969\) and \(1.969\). This is a direct consequence of the quadratic outcome relationship and the assumed predictor distribution.

A sampler that explores only one mode can produce an imputation chain that looks stable locally but fails to represent the other plausible region. This is a sampling failure, not evidence that one mode can be discarded. Starting chains in different regions and examining trace plots and between-chain behavior can help reveal the problem. A sampler must target the full conditional density, including its separated regions.

Implementing an imputation that respects the quadratic term

For each missing value of \(X\), a compatible update targets \(f(x\mid y,z)\), rather than an unrelated regression chosen only because \(X\) is the variable being imputed. The density can be sampled with a suitable numerical method, such as slice sampling, Metropolis–Hastings, or a carefully checked grid-based method. The method must explore the entire relevant support and, especially for a multimodal density, move adequately between modes.

Specify the two model components.
State the predictor distribution given observed covariates and the quadratic outcome model, including their error distributions and variance assumptions.

Construct the conditional target.
Multiply the outcome density at the observed \(Y\) by the predictor density at each candidate \(x\); normalize if the chosen sampler requires a normalized density.

Sample missing predictor values.
Use a sampler that explores the target's full shape. Where feasible, draw model parameters as well as predictor values so the imputations reflect parameter uncertainty.

Recompute the analysis term.
After drawing \(X\), calculate \(X^2\) from that same draw. Do not impute an independently varying squared term.

In a chained-equations workflow, compatibility is about the conditional distribution being used for the incomplete predictor, not merely whether a quadratic term appears somewhere in the analysis formula. If an imputation procedure uses a convenient regression approximation instead of the derived conditional, it should not be assumed to reproduce the target automatically. Parameter uncertainty, convergence, and the distributional assumptions remain important.

Worked Example: Why an Ordinary Linear Imputation Can Mislead

Use the hypothetical models \(X\sim N(0,1)\) and \(Y\mid X\sim N(5+X^2,1)\), and suppose the observed outcome is \(Y=9\). Compare the conditional density kernels at \(x=0\) and \(x=2\).

At \(x=0\), the predicted outcome is 5, so the outcome residual is \(9-5=4\). The log kernel is \(-4^2/2-0^2/2=-8\). At \(x=2\), the prediction is 9 and the predictor log-density contribution is \(-2^2/2=-2\), so the log kernel is \(-0^2/2-2=-2\). The ratio of the kernel at zero to the kernel at two is

\frac{\exp(-8)}{\exp(-2)}=\exp(-6)\approx0.002479.

The model therefore gives a very low conditional density at zero compared with the density at two; by symmetry, the density at \(-2\) is the same as at \(2\). A simple linear regression of \(X\) on \(Y\) cannot represent this two-sided conditional shape with a single normal distribution. In this example, imputation concentrated around zero would place substantial weight in a region the specified joint model makes relatively implausible. The compatible conditional distribution instead retains both signs as plausible.

Diagnostics, alternatives, and limitations

Several checks are useful after specifying or implementing the conditional model:

Check the model assumptions. Examine whether a normal predictor model given \(Z\) and a normal-error quadratic outcome model are defensible. Skewness, heavy tails, or influential observations can make these assumptions inadequate.
Check the conditional shape and sampling. Plot conditional density curves for representative observed \(Y,Z\) values. For simulation-based draws, inspect trace plots and compare chains for mixing, including movement between modes.
Check the analysis relationship after imputation. Compare the observed and imputed predictor values across outcome ranges, and verify that \(X^2\) in each completed data set is calculated from its corresponding \(X\).
Do not mistake fit diagnostics for a missingness test. Imputed-versus-observed comparisons can identify conspicuous model problems, but they cannot establish that a missingness assumption holds.

The main alternative is a convenient regression imputation, such as a normal linear regression of \(X\) on \(Y\) and \(Z\). It is easier to fit and may be an approximation in some settings, but it does not generally equal the conditional density derived here. Including \(Y^2\) as a predictor in that regression does not, by itself, make it the exact conditional model. The exact conditional is determined by the predictor distribution and the quadratic outcome model together.

A joint model that specifies a distribution for \(X\) and the quadratic regression for \(Y\) offers a direct route to the needed conditional distribution. A flexible approximation may be more practical in a given application, but its adequacy should be assessed rather than presumed. Neither approach removes the need to model the data appropriately. The result also concerns a single incomplete continuous predictor under particular distributional assumptions; other predictor types, multiple incomplete predictors, or other outcome models require their own conditional derivations.

As with other multiple-imputation procedures, valid inference requires appropriate missing-data and model assumptions, and imputations should reflect uncertainty in unknown parameters. Treating fitted parameters as known can understate uncertainty. A technically correct conditional density also cannot compensate for a misspecified outcome relationship or an unsuitable predictor model. The aim is not to force every application into this normal-quadratic specification, but to make the imputation model respect the nonlinear term the analysis intends to estimate.

Common pitfall: “The analysis includes \(X^2\), so I will add \(X^2\) as a predictor in the regression used to impute \(X\).” This is not a valid general recipe: \(X^2\) is unknown when \(X\) is missing and must equal the square of the imputed \(X\). Derive or otherwise justify the conditional distribution of \(X\mid Y,Z\), then calculate \(X^2\) from each imputed value.

Check Your Understanding

Answer each question, then check the solutions below.

Under the normal predictor and normal-error quadratic outcome models, what two density components determine \(f(x\mid y,z)\)?
Why is the conditional distribution generally not normal when \(\beta_2\ne0\)?
For \(X\sim N(0,1)\), \(Y\mid X\sim N(5+X^2,1)\), and \(Y=9\), what is the ratio of the conditional density kernels at \(x=0\) and \(x=2\)?
After imputing \(X\), how should the analysis term \(X^2\) be obtained?

Solutions. (1) It is proportional to \(f(y\mid x,z)f(x\mid z)\): the outcome density and the predictor density. (2) The squared residual from a quadratic mean produces a fourth-degree term in the log density, so its shape is not generally the quadratic log density of a normal distribution. (3) At zero the log kernel is \(-8\); at two it is \(-2\). The ratio is \(\exp(-8)/\exp(-2)=\exp(-6)\approx0.002479\). (4) Square each imputed value of \(X\); do not impute \(X^2\) independently.