Restricted Cubic Splines for Covariate Relationships

Why Use a Restricted Cubic Spline?

A linear covariate effect may be too rigid when a model parameter changes at different rates across the covariate range. A polynomial can represent curvature, but its behavior may become extreme near the edges or outside the observed range. A restricted cubic spline offers a flexible middle ground: it joins smooth cubic pieces across specified knots while requiring the relationship to be linear beyond the boundary knots.

In an NLME model, the spline describes how a typical model parameter depends on a continuous covariate. For example, it could describe how typical clearance varies with body weight. The spline does not itself specify the random-effects distribution or the residual model; those remain parts of the broader model. Here the focus is how to build and evaluate the fixed covariate function.

Definition: A restricted cubic spline is a function formed from a linear term and a set of truncated-cubic basis functions. It is twice continuously differentiable, and it is constrained to be linear below the smallest knot and above the largest knot. The smallest and largest knots are called the boundary knots; knots between them are interior knots.

Constructing the Basis

Let \(x\) be a continuous covariate and let its \(K\) ordered knots be \(\xi_1<\xi_2<\cdots<\xi_K\), where \(K\geq 3\). A truncated power \(u_+^3\) means \(u^3\) when \(u>0\), and 0 otherwise. For each \(j=1,\ldots,K-2\), define the restricted cubic spline basis function

d_j(x)=(x-\xi_j)_+^3 -\frac{\xi_K-\xi_j}{\xi_K-\xi_{K-1}}(x-\xi_{K-1})_+^3 +\frac{\xi_{K-1}-\xi_j}{\xi_K-\xi_{K-1}}(x-\xi_K)_+^3.

A spline contribution to a model parameter's predictor is then

f(x)=\beta_0+\beta_1 x+\sum_{j=1}^{K-2}\gamma_j d_j(x).

There are \(K-2\) nonlinear basis functions, plus the intercept and linear term. Thus, with four knots, there are two nonlinear basis functions. The knots define where the basis can change curvature; the fitted coefficients \(\gamma_j\) determine the actual shape. A knot does not force a peak, a trough, or any particular change in direction.

In practice, it is often useful to rescale the covariate first. For a reference value \(x_{\mathrm{ref}}\) and scale \(s>0\), define \(z=(x-x_{\mathrm{ref}})/s\), place the knots on the \(z\) scale, and construct \(d_j(z)\). The choice of units changes the numerical basis values and coefficient sizes, but not the fitted relationship when the model is transformed consistently.

Worked Example: Building Two Basis Values

Suppose the four knots on a rescaled covariate \(z\) are \(0,1,2,3\). Compute the basis values at \(z=2.5\). These are hypothetical design values, not observations from a real study.

Here \(K=4\), so there are \(K-2=2\) nonlinear basis functions. The denominator in the basis formula is \(\xi_4-\xi_3=3-2=1\). For \(j=1\), the knot-dependent multipliers are \((3-0)/1=3\) and \((2-0)/1=2\). Since the positive-part terms at knot 3 are zero at \(z=2.5\),

$$ d_1(2.5)=(2.5-0)^3-3(2.5-2)_+^3+2(2.5-3)_+^3 =15.625-3(0.5)^3+0=15.250. $$

For \(j=2\), the multipliers are \((3-1)/1=2\) and \((2-1)/1=1\), giving

$$ d_2(2.5)=(2.5-1)^3-2(2.5-2)_+^3+(2.5-3)_+^3 =3.375-2(0.5)^3+0=3.125. $$

Thus the basis vector at \(z=2.5\) is \((15.250,\,3.125)\). Large-looking values are possible because the raw basis uses cubes; their size depends on the covariate scale. Coefficients are estimated for this chosen basis, so one should not interpret a basis value or coefficient in isolation as a covariate effect.

Why the Tails Are Linear

Below the smallest knot \(\xi_1\), all positive-part terms in every \(d_j(x)\) are zero. The spline contribution from the nonlinear basis is therefore zero there, and \(f(x)\) reduces to \(\beta_0+\beta_1x\).

Above the largest knot \(\xi_K\), all the positive-part terms are active. For a given \(j\), write \(A=(\xi_K-\xi_j)/(\xi_K-\xi_{K-1})\) and \(B=(\xi_{K-1}-\xi_j)/(\xi_K-\xi_{K-1})\). The coefficient of \(x^3\) in \(d_j(x)\) is \(1-A+B=0\). The coefficient of \(x^2\) is \(-3\xi_j+3A\xi_{K-1}-3B\xi_K=0\). Hence \(d_j(x)\) is at most linear above \(\xi_K\). This cancellation is the restriction that prevents cubic growth in the tails.

At each knot, the truncated-cubic components and their first two derivatives join continuously. The resulting function is smooth through the knots, with curvature changing gradually rather than through abrupt corners. Smoothness does not mean curvature is constant, nor does it mean that the function must be monotone.

Key result: With \(K\) ordered knots, the restricted cubic spline has \(K-2\) nonlinear basis functions. It is linear outside the boundary knots and smooth at the knots. Its interior curvature is determined jointly by the knot locations and the fitted coefficients.

Calculating a Covariate-Dependent Parameter

A spline can be placed directly on a parameter scale when that scale is appropriate. If the parameter must remain positive, a common choice is to place the spline on its log scale:

\log P_{\mathrm{typ}}(z)=\beta_0+\beta_1z+\sum_{j=1}^{K-2}\gamma_jd_j(z), \qquad P_{\mathrm{typ}}(z)=\exp\left(\beta_0+\beta_1z+\sum_{j=1}^{K-2}\gamma_jd_j(z)\right).

Here \(P_{\mathrm{typ}}(z)\) is the typical value of the parameter at covariate value \(z\). The coefficients are treated as fitted values when calculating a prediction. A log-scale spline guarantees a positive predicted parameter, but it does not by itself ensure that the covariate function is biologically plausible or supported by the data.

Worked Example: Predicting Typical Clearance

Consider a hypothetical fitted model for typical clearance \(CL_{\mathrm{typ}}\), measured in L/h. Let \(z=(\mathrm{WT}-70\ \mathrm{kg})/(10\ \mathrm{kg})\), with knots \(0,1,2,3\), corresponding to weights 70, 80, 90, and 100 kg. Suppose the fitted coefficients are \(\beta_0=\log(4.00)\), \(\beta_1=0.100\), \(\gamma_1=-0.00400\), and \(\gamma_2=0.0120\). At 95 kg, \(z=2.5\), and the basis values from the previous calculation are \(d_1=15.250\) and \(d_2=3.125\).

Substitute these values into the fitted log-parameter equation:

\log CL_{\mathrm{typ}}(2.5) =\log(4.00)+0.100(2.5)-0.00400(15.250)+0.0120(3.125) =\log(4.00)+0.2265.

Exponentiating gives

CL_{\mathrm{typ}}(2.5)=4.00\exp(0.2265) \approx 4.00(1.2542)=5.017\ \mathrm{L/h}.

Under this hypothetical fitted relationship, typical clearance at 95 kg is approximately 5.017 L/h. This is a model-based typical parameter value, not an individual prediction incorporating that person's random effect.

Worked Example: Comparing an Interior Value With a Tail Value

Use the same hypothetical model and coefficients to calculate \(CL_{\mathrm{typ}}\) at 75 kg and 105 kg. These values illustrate the spline's behavior; the 105 kg value is outside the upper boundary knot and should be treated as extrapolation.

At 75 kg, \(z=0.5\), which is below the first knot. Both basis functions are zero. Therefore

\log CL_{\mathrm{typ}}(0.5)=\log(4.00)+0.100(0.5)=\log(4.00)+0.0500, \qquad CL_{\mathrm{typ}}(0.5)=4.00e^{0.0500}\approx4.205\ \mathrm{L/h}.

At 105 kg, \(z=3.5\), above the last knot. Direct substitution gives \(d_1(3.5)=3.5^3-3(1.5)^3+2(0.5)^3=33.000\) and \(d_2(3.5)=2.5^3-2(1.5)^3+(0.5)^3=9.000\). Thus

\log CL_{\mathrm{typ}}(3.5) =\log(4.00)+0.100(3.5)-0.00400(33.000)+0.0120(9.000) =\log(4.00)+0.326,

CL_{\mathrm{typ}}(3.5)=4.00e^{0.326}\approx5.542\ \mathrm{L/h}.

The fitted log-parameter predictor is linear beyond the boundary knot, but its slope need not equal the slope below the lower boundary. Also, the spline's linear-tail constraint does not make extrapolation reliable: the covariate relationship outside the observed range may differ from that imposed by the model.

Assumptions, Diagnostics, and Alternatives

A restricted cubic spline is a representation of a covariate relationship, not a guarantee that the relationship has been correctly specified. Its use assumes that a continuous, smooth function is a reasonable description over the modeled range and that the chosen boundary-tail behavior is acceptable. Knot locations should lie in a region with meaningful covariate support, and there must be enough information across the range to estimate the additional coefficients. With sparse data, apparent curvature can be driven by a few observations or become weakly identified.

Assess the fitted relationship on a plot over the observed covariate range, including the boundary knots and the data distribution. Look for abrupt or implausible changes in slope, excessive curvature in regions with little data, and predictions that conflict with the parameter's scientific meaning. In an NLME analysis, also check whether the covariate function is supported by the model fit and whether conclusions are sensitive to reasonable knot choices. The function's shape and uncertainty matter more than a test of one spline coefficient at a time.

The earlier material on piecewise linear covariate effects uses linear segments meeting at specified knots; it can be easier to interpret but generally permits changes in slope at the knots. A single linear effect is simpler but cannot represent curvature. A polynomial is compact but may behave poorly near the covariate extremes. A restricted cubic spline trades additional coefficients for smooth, flexible curvature while preserving linear tails. The choice should reflect scientific expectations, data support, and the intended use of predictions.

Common pitfall: It is incorrect to assume that each knot is a peak or that adding knots automatically improves the model. Knots set locations where the basis can contribute to curvature; the fitted coefficients determine the shape. Adding knots also adds nonlinear basis functions and can increase estimation uncertainty, especially when data are sparse.

Choosing and Reporting a Spline

The number and locations of knots are modeling choices. Knots are often placed at selected quantiles of the covariate distribution so that the fitted curvature is supported across the data range. Whatever rule is used, report the knot locations and the covariate scale. If the covariate is centered and rescaled, state the reference value and scale as well; otherwise, basis values and coefficients cannot be reproduced or interpreted correctly.

Avoid relying on the spline to establish behavior far beyond the observed data. Its linear tails are a mathematical constraint, not evidence that the true relationship is linear there. Likewise, a visually smooth fitted curve can still be uncertain. Interpret the estimated relationship in regions where the data provide support, and consider uncertainty in the fitted curve when making substantive conclusions.

Check Your Understanding

Use the definitions and calculations above to answer these questions. The solutions follow each question.

How many nonlinear basis functions are used with six knots? Solution: The number is \(K-2=6-2=4\).
With knots \(0,1,2,3\), what are \(d_1(z)\) and \(d_2(z)\) at \(z=0.5\)? Solution: This is below the first knot, so every positive-part term is zero. Hence \(d_1(0.5)=d_2(0.5)=0\).
For the hypothetical clearance model above, why does the predicted parameter at \(z=2.5\) equal about 5.017 L/h rather than the log predictor value? Solution: The model specifies \(\log CL_{\mathrm{typ}}\), so the parameter is obtained by exponentiating: \(\log CL_{\mathrm{typ}}=\log(4.00)+0.2265\), and \(CL_{\mathrm{typ}}=4.00e^{0.2265}\approx5.017\) L/h.
What does the restricted-tail condition imply beyond the largest knot, and what does it not imply? Solution: It implies that the spline function is linear beyond that knot. It does not imply that the extrapolated relationship is scientifically correct or supported by observations.