Spline Terms and Imputation Compatibility

Why a Spline Changes the Imputation Problem

A spline lets an analysis model represent a relationship whose slope changes across the range of a predictor. If that predictor is incomplete, imputing it with a conditional model that does not reflect the spline relationship can distort the association the analysis is intended to estimate. The issue is not resolved simply by adding the observed outcome to the predictor’s imputation equation: the form of the conditional distribution matters too.

Consider a continuous outcome \(Y\), an incomplete continuous predictor \(X\), and fully observed covariates \(Z\). A one-knot piecewise-linear analysis model can be written using the hinge function \((X-k)_+=\max(0,X-k)\), where \(k\) is a specified knot. The spline is linear on each side of the knot, continuous at the knot, and generally has a kink there.

Definition: A one-knot piecewise-linear spline with knot \(k\) uses the basis terms \(X\) and \((X-k)_+\). In a model with conditional mean \(\alpha+\gamma Z+\beta_1X+\beta_2(X-k)_+\), the slope is \(\beta_1\) for \(X<k\) and \(\beta_1+\beta_2\) for \(X>k\). At \(X=k\), the fitted mean is continuous; if \(\beta_2\ne 0\), it generally has no unique derivative.

For an incomplete \(X\), the relevant imputation target is \(X\mid Y,Z\), not just the analysis mean \(Y\mid X,Z\). A compatible conditional must reflect the likelihood contribution of the spline analysis model as well as a model for the distribution of \(X\) given \(Z\). This is a stronger requirement than including \(Y\) as a single linear predictor in a regression for \(X\).

Deriving the Conditional for a Missing Predictor

For a transparent derivation, assume the following model. Conditional on \(Z\), the predictor follows a normal distribution with mean \(\mu(Z)\) and variance \(\tau^2\). Given \(X,Z\), the outcome is normal with constant variance \(\sigma^2\) and a piecewise-linear mean. Assume the parameters and knot are fixed for this derivation. The joint density factors as \(p(X\mid Z)p(Y\mid X,Z)\); therefore the conditional density of \(X\) given \(Y,Z\) is proportional to their product.

X\mid Z\sim N(\mu(Z),\tau^2),\qquad Y\mid X,Z\sim N\!\left(\alpha+\gamma Z+\beta_1X+\beta_2(X-k)_+,\sigma^2\right)

On either side of \(k\), the outcome mean is an ordinary line in \(X\). Write it as \(a_j+b_jX\), incorporating the fixed value of \(Z\) into \(a_j\). Below the knot, \(a_1=\alpha+\gamma Z\) and \(b_1=\beta_1\). Above the knot, \(a_2=\alpha+\gamma Z-\beta_2k\) and \(b_2=\beta_1+\beta_2\). Multiplying the normal prior density and normal outcome likelihood gives a quadratic exponent in \(X\) within each region.

p(x\mid y,z)\ \propto\ \exp\!\left\{-\frac{(x-\mu)^2}{2\tau^2} -\frac{(y-a_j-b_jx)^2}{2\sigma^2}\right\}, \quad x\in R_j

Completing the square shows that the density within region \(R_j\) is proportional to a normal density with variance \(v_j\) and mean \(m_j\), restricted to that region. The two regions have different means and variances in general. Their normal kernels also have different constants after completing the square, so their relative posterior probabilities are not obtained merely by comparing \(m_1\) and \(m_2\).

v_j=\left(\frac{1}{\tau^2}+\frac{b_j^2}{\sigma^2}\right)^{-1}, \qquad m_j=v_j\left\{\frac{\mu}{\tau^2} +\frac{b_j(y-a_j)}{\sigma^2}\right\}

To draw from the full conditional, calculate each region’s integrated density, including its completing-the-square constant and the probability that the corresponding normal kernel falls in that region. Normalize those two integrals to obtain the probabilities of being below or above the knot. First select a region using those probabilities; then draw from its normal kernel truncated to that region. This procedure follows directly from the product \(p(x\mid z)p(y\mid x,z)\).

Assumptions for this derivation: The predictor model is normal given \(Z\); the outcome model is normal with constant variance and a one-knot piecewise-linear mean; and the parameters and knot are treated as fixed. This gives the stated conditional under this model. It does not establish compatibility with arbitrary conditional models for other variables in a chained-equations system.

Worked Examples

Worked Example: Reading the Spline Basis Correctly

Suppose the fitted mean is \(12+1.2X+0.8(X-5)_+\), with outcome measured in points and \(X\) measured in years. The knot is \(k=5\).

For \(X<5\), the hinge is zero, so the fitted mean is \(12+1.2X\), with slope \(1.2\) points per year. For \(X>5\), the mean is \(12+1.2X+0.8(X-5)=8+2.0X\), with slope \(2.0\) points per year. At \(X=5\), both expressions give \(18\) points, so the fitted mean is continuous. The left and right slopes differ, so there is no unique derivative at the knot.

At \(X=6.2\), the hinge value is \(6.2-5=1.2\). The fitted mean is \(12+1.2(6.2)+0.8(1.2)=12+7.44+0.96=20.4\) points. Omitting the hinge would give \(12+1.2(6.2)=19.44\) points, a difference of \(20.4-19.44=0.96\) points. That difference illustrates why the hinge is part of the analysis relationship that imputation must preserve.

Worked Example: Calculating the Two Conditional Kernels

Suppose, hypothetically, that \(X\mid Z\sim N(4,4)\), the outcome variance is \(\sigma^2=4\), and the outcome mean is \(10+X+2(X-5)_+\). For a case with \(Y=18\), calculate the region-specific conditional normal parameters. The first argument of \(N(\cdot,\cdot)\) is the mean and the second is the variance.

Below the knot, \(a_1=10\) and \(b_1=1\). Thus \(v_1=(1/4+1/4)^{-1}=2\), and \(m_1=2\{4/4+1(18-10)/4\}=2(1+2)=6\). The below-knot contribution is proportional to a \(N(6,2)\) density restricted to \(X\leq 5\).

Above the knot, \(10+X+2(X-5)=3X\), so \(a_2=0\) and \(b_2=3\). Thus \(v_2=(1/4+9/4)^{-1}=0.4\), and \(m_2=0.4\{4/4+3(18-0)/4\}=0.4(1+13.5)=5.8\). The above-knot contribution is proportional to a \(N(5.8,0.4)\) density restricted to \(X>5\).

These are truncated kernels, not two alternative full conditional distributions. To generate a draw, their integrated contributions must first be normalized to determine which region applies. In particular, the unrestricted kernel means alone do not determine the region probabilities.

Worked Example: Why a Linear Imputation Can Miss the Spline

Use the preceding analysis mean, \(12+1.2X+0.8(X-5)_+\). Suppose an imputation equation for missing \(X\) uses \(Y\) only as a linear predictor, with no mechanism for the change in the outcome relationship across the knot. At \(X=6.2\), the spline analysis mean is \(20.4\) points, while the mean from the analysis relationship with the hinge omitted is \(19.44\) points. The discrepancy is \(0.96\) points.

This calculation does not prove that every linear regression of \(X\) on \(Y\) is biased by exactly \(0.96\), or even in a particular direction. It shows that the analysis likelihood treats values on the two sides of the knot differently. A simple linear conditional for \(X\mid Y,Z\) does not in general reproduce the piecewise-normal conditional derived from that likelihood and a predictor model. A model-based draw using the two region kernels instead carries the specified change in slope into the imputation.

Practical Implementation and Checks

A direct way to preserve the spline relationship is to use a substantive-model-compatible conditional for the missing predictor. Under the assumptions above, this means sampling from the density proportional to the analysis likelihood times the predictor model. In a chained-equations procedure, such a conditional can be used for the update of \(X\), while other incomplete variables are updated using their own conditional models. The resulting set of conditionals is not automatically guaranteed to correspond to one coherent joint distribution; compatibility must be considered for the system as a whole.

When a specialized sampler is unavailable, an approximate imputation model may be considered, but its adequacy should be assessed rather than assumed. A useful check is whether imputed values preserve the observed relationship between \(X\) and \(Y\) on both sides of \(k\), including the change in slope. Compare observed and imputed distributions of \(X\), the proportions on each side of the knot, and outcome-versus-predictor patterns across missingness groups. These checks can reveal implausible concentrations, an incorrect knot-region frequency, or a relationship that appears linear after imputation despite a spline in the analysis.

Specify the analysis relationship.
Write the hinge basis and knot explicitly, and state which predictor is incomplete.

Build the conditional target.
Combine the model for \(X\mid Z\) with the outcome likelihood \(Y\mid X,Z\); do not substitute the analysis mean alone for \(X\mid Y,Z\).

Check the two regions.
Verify the slopes, continuity at the knot, region-specific conditional kernels, and the way their probabilities are normalized.

Inspect completed data.
Compare imputed and observed predictor distributions and check that the spline-shaped outcome association is not lost.

Alternatives, Limitations, and a Common Pitfall

The direct model-based approach derives \(X\mid Y,Z\) from a specified model for \(X\mid Z\) and the spline outcome likelihood. It can represent the knot-induced change exactly under those assumptions, but requires suitable sampling and estimation. An ordinary linear regression imputation for \(X\) is simpler, but generally imposes a different conditional shape. A more flexible approximate conditional may capture curvature or changing association, but flexibility alone does not guarantee compatibility with the substantive spline model. The choice should be judged by the target analysis relationship and by diagnostics, not merely by which model is easiest to fit.

The derivation also has limits. A normal predictor distribution may be implausible, the outcome variance may not be constant, and a one-knot linear spline may not describe the scientific relationship. Parameters are treated as fixed in the algebra; a full imputation procedure must account for parameter uncertainty rather than treating estimated parameters as known. If the knot is data-selected, uncertainty from selecting it is not handled by this conditional derivation. Missingness assumptions and adequate model specification remain necessary; matching a spline term does not by itself solve nonignorable missingness or poor overlap.

Common pitfall: “The analysis model contains \(X\) and its hinge, so I will just include the outcome linearly when imputing \(X\).” This reasoning confuses \(Y\mid X,Z\) with \(X\mid Y,Z\). The latter depends on the predictor distribution as well as the outcome likelihood, and under the stated model has different truncated normal kernels on the two sides of the knot. Including the outcome is important, but a linear outcome term alone does not generally reproduce those kernels.

Key takeaway: A spline in an incomplete predictor creates a piecewise-shaped likelihood contribution to its imputation conditional. To preserve the analysis relationship, account for the hinge and the predictor distribution when constructing \(X\mid Y,Z\); then assess whether the imputed data retain the knot-specific pattern.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

For the mean \(3+2X-0.5(X-4)_+\), state the slope below and above the knot and the fitted mean at the knot.
Why is the conditional for an incomplete \(X\) not determined by the analysis mean \(Y\mid X,Z\) alone?
Under the assumptions in the derivation, what two quantities determine the region-specific variance \(v_j\)?
What must be done with the two region-specific kernels before drawing an imputed \(X\)?

Solutions. (1) The slopes are \(2\) below \(4\) and \(2-0.5=1.5\) above \(4\). At the knot the hinge is zero, so the fitted mean is \(3+2(4)=11\); the mean is continuous there. (2) The conditional \(X\mid Y,Z\) is proportional to \(p(X\mid Z)p(Y\mid X,Z)\), so it depends on both the predictor distribution and the outcome likelihood. (3) The region-specific slope \(b_j\), the predictor variance \(\tau^2\), and the outcome variance \(\sigma^2\) enter \(v_j=(1/\tau^2+b_j^2/\sigma^2)^{-1}\). (4) Integrate each kernel over its own region, normalize those integrals to get the region probabilities, select a region, and draw from its normal distribution truncated to that region.