A Worked Spline Analysis of a Nonlinear Dose Response

From Dose Measurements to Fitted Curves

A dose-response relationship describes how the mean outcome changes with dose. In many applications, that change is not well represented by one straight line: the response may rise slowly at low doses and more quickly at higher doses, or its rate of change may vary smoothly. Piecewise regression and splines give two ways to represent such curvature. This worked analysis uses the same hypothetical dataset for both, so their fitted values and boundary behavior can be compared directly.

Suppose dose \(x\) is measured in milligrams and outcome \(y\) in response units. The five hypothetical observations below are illustrative, not measurements from a real study.

Dose \(x\) (mg)	Observed outcome \(y\) (response units)
0	2.00
2	4.44
4	6.72
6	10.60
8	19.12

For illustration, treat the observed doses as fixed and suppose outcomes at different doses come from independent experimental units. A regression model describes the conditional mean outcome at each dose; the observed outcome can differ from that mean because of an error. The usual working assumptions for the least-squares analyses here are that the errors have mean zero, are independent, and have a common variance. These assumptions are not established by this small table. In particular, the examples demonstrate fitting and interpretation, not reliable inferential conclusions.

A Continuous Piecewise Linear Fit

A continuous piecewise linear model with a knot at 4 mg has one straight segment on each side of the knot, joined without a jump. Let \(t=x-4\), so the knot is at \(t=0\), and write \(t_+=\max(t,0)\). The model and its conditional mean are

Definition: A continuous piecewise linear regression model with knot \(\kappa\) can be written as \(E(Y\mid x)=\beta_0+\beta_1(x-\kappa)+\beta_2(x-\kappa)_+\). The fitted mean is continuous at \(\kappa\). Its slope to the left is \(\beta_1\), and its slope to the right is \(\beta_1+\beta_2\).

The term \((x-\kappa)_+\) equals zero at and below the knot, and equals \(x-\kappa\) above it. Thus \(\beta_2\) is the change in slope at the knot, not the slope of the upper segment by itself. Here the model is \(\hat y=a+bt+ct_+\); its fitted slope is \(b\) below 4 mg and \(b+c\) above 4 mg.

Worked Example: Fitting the Continuous Piecewise Linear Model

For the five observations, the values of \(t=x-4\) are \(-4,-2,0,2,4\), and the corresponding values of \(t_+\) are \(0,0,0,2,4\). Fit \(\hat y=a+bt+ct_+\) by least squares.

Least squares chooses \(a,b,c\) to minimize the sum of squared residuals, \(\sum_i(y_i-a-bt_i-ct_{i+})^2\). Setting each partial derivative to zero gives the normal equations. The needed sums from the table are \(\sum t=0\), \(\sum t_+=6\), \(\sum t^2=40\), \(\sum tt_+=20\), \(\sum t_+^2=20\), \(\sum y=42.88\), \(\sum ty=80.80\), and \(\sum t_+y=97.68\).

\begin{aligned} 5a+6c&=42.88,\\ 40b+20c&=80.80,\\ 6a+20b+20c&=97.68. \end{aligned}

For example, the first equation follows by summing the fitted values and equating them to \(\sum y\); the second and third follow after multiplying by \(t\) and \(t_+\), respectively, and summing. From the second equation, \(b=2.02-0.5c\). From the first, \(a=8.576-1.2c\). Substituting both into the third gives \(6(8.576-1.2c)+20(2.02-0.5c)+20c=97.68\), so \(91.856+2.8c=97.68\), and \(c=2.08\). Therefore \(b=0.980\) and \(a=6.080\).

The fitted equation is \(\hat y=6.080+0.980t+2.080t_+\). Its left-segment slope is \(0.980\) response units per mg, and its right-segment slope is \(0.980+2.080=3.060\) response units per mg. The fitted mean at the knot is \(6.080\) response units from either segment, so there is no jump.

The fitted values at the observed doses are \(2.160,\ 4.120,\ 6.080,\ 12.200,\ 18.320\) response units. For example, at \(x=6\), \(t=2\) and \(t_+=2\), giving \(\hat y=6.080+0.980(2)+2.080(2)=12.200\). The residual there is \(10.60-12.20=-1.60\) response units. Across all five observations the residual sum of squares is \(3.7376\) response-units squared. The nonzero residuals show that the piecewise line does not pass through every observation.

A Cubic Regression Spline Fit

A cubic regression spline uses cubic polynomial pieces joined with smoothness constraints. For this example, use a cubic regression spline with one interior knot at 4 mg, expressed in the truncated-power basis. This is an ordinary cubic regression spline, not a natural cubic spline. The model has five coefficients:

E(Y\mid x)=\gamma_0+\gamma_1x+\gamma_2x^2+\gamma_3x^3+\gamma_4(x-4)_+^3.

The hinge-cubic term is zero for \(x\leq 4\) and equals \((x-4)^3\) for \(x>4\). It allows the cubic polynomial to change above the knot while keeping the fitted curve and its first two derivatives continuous there. The five-dose example has exactly five observations, so a full-rank five-parameter model can interpolate all five. An exact fit in this setting is a consequence of flexibility and the chosen data, not proof that the curve is appropriate for other dose-response data.

Worked Example: Fitting the Cubic Spline to the Same Data

Consider the coefficient values \(\gamma_0=2\), \(\gamma_1=1.5\), \(\gamma_2=-0.20\), \(\gamma_3=0.030\), and \(\gamma_4=0.040\). Substitution gives the fitted curve \(\hat y=2+1.5x-0.20x^2+0.030x^3+0.040(x-4)_+^3\).

At \(x=0\), the prediction is \(2\). At \(x=2\), it is \(2+1.5(2)-0.20(2^2)+0.030(2^3)=4.44\). At \(x=4\), the hinge term is zero, giving \(2+6-3.2+1.92=6.72\). At \(x=6\), it is \(2+9-7.2+6.48+0.040(2^3)=10.60\). At \(x=8\), it is \(2+12-12.8+15.36+0.040(4^3)=19.12\). These are exactly the five observed outcomes.

The design matrix has full rank for these five distinct dose values: its columns are \(1,x,x^2,x^3,(x-4)_+^3\). The first four observations at \(x=0,2,4,6\) give a full-rank cubic-polynomial block for the first four columns, while the last column is zero at those doses and positive at \(x=8\). Thus the coefficients are uniquely determined by these exact interpolating values. The residual sum of squares is zero, but there are no residual degrees of freedom to assess how well the model generalizes.

Predictions and Boundary Behavior

A fitted curve is often used to predict the mean outcome at doses between observed values. Such predictions are interpolations. Predictions outside the observed dose range are extrapolations and depend strongly on how the model extends beyond the data. The following calculations compare the two fitted curves at three doses; their values are predictions of mean response, not guaranteed outcomes for individual experimental units.

Worked Example: Comparing Predictions at 3, 5, and 7 mg

At 3 mg, the piecewise fit has \(t=-1\) and \(t_+=0\), so \(\hat y=6.080+0.980(-1)=5.100\). The cubic spline gives \(2+1.5(3)-0.20(3^2)+0.030(3^3)=5.510\).

At 5 mg, the piecewise fit has \(t=t_+=1\), so \(\hat y=6.080+0.980+2.080=9.140\). For the cubic spline, the hinge is \(1^3\), giving \(2+7.5-5+3.75+0.040=8.290\).

At 7 mg, the piecewise fit gives \(6.080+0.980(3)+2.080(3)=15.260\). The cubic spline gives \(2+10.5-0.20(49)+0.030(343)+0.040(27)=14.070\). The predictions differ because the models distribute curvature differently between the observed doses. With only five observations, this comparison alone cannot establish which curve is preferable.

Boundary behavior: Below the smallest observed dose, this piecewise fit continues its leftmost straight line, with slope \(0.980\) response units per mg. Above the largest observed dose, it continues its rightmost straight line, with slope \(3.060\) response units per mg. This cubic regression spline instead continues its outer cubic polynomial pieces; it is not constrained to become linear at the boundaries. Neither extension is validated merely because the curve fits the observed range.

For example, at 10 mg the piecewise prediction is \(6.080+3.060(10-4)=24.440\) response units. The cubic spline prediction is \(2+15-20+30+0.040(6^3)=35.640\) response units. The widening difference illustrates why an extrapolated value can be highly sensitive to the chosen form. It is not evidence that either prediction is accurate. Subject-matter constraints and data collected near the target dose would be important to its interpretation.

Why the Comparison Matters

A continuous piecewise linear fit makes its segment slopes directly readable and changes slope abruptly at its knot. A cubic regression spline can represent smoother changes in slope and curvature, but its coefficients are not generally interpretable as simple segment slopes. The earlier material on interpreting spline curves emphasizes reading changes from the fitted curve rather than treating every basis coefficient as a slope.

Both fits also inherit the limitations of the data and assumptions. Here the five-parameter cubic spline uses all five observations to determine its coefficients, leaving no residual degrees of freedom. Its zero residual sum of squares is therefore not a fair demonstration that it captures the underlying dose-response process better than the three-parameter piecewise linear model. The piecewise model has nonzero residuals, but that fact alone does not show it is unsuitable. More observations, especially repeated measurements or measurements near the region of interest, would be needed for a meaningful assessment.

The least-squares calculations above assume independent errors with mean zero and common variance for the stated regression mean. If those assumptions are implausible, or if the dose measurements do not support the selected form, the numerical fit does not repair the problem. Likewise, a fitted association between dose and outcome does not by itself establish a causal effect.

Common Mistakes

Calling the hinge coefficient the upper-segment slope. In the piecewise model, \(2.080\) is the change in slope. The upper-segment slope is \(0.980+2.080=3.060\) response units per mg.
Assuming an exact spline fit is automatically best. Here the spline has as many coefficients as observations. Zero residual error is guaranteed by interpolation and does not demonstrate reliable predictions.
Treating a boundary prediction as if it were observed support. The models extend differently beyond 0–8 mg. A numerical prediction at 10 mg is an extrapolation, not a measured response.
Interpreting a fitted mean as an individual outcome. A prediction of 9.140 response units at 5 mg is the model's estimated mean response at that dose, not a promise that each individual response will equal 9.140.

Key takeaway: Fit each curve using its stated model, interpret its coefficients according to its basis, and calculate predictions from the fitted equation. Compare boundary extensions cautiously: a curve's behavior outside the observed dose range is determined by its form, not by evidence in the data.

Check Your Understanding

Use the equations and reasoning in this analysis to answer the following questions.

For the piecewise model, what are the fitted slopes below and above the 4 mg knot, and what does the hinge coefficient represent?
Calculate the piecewise-model prediction at 4 mg. Why does this model have the same fitted value from both sides of the knot?
Why does a zero residual sum of squares for the cubic spline not establish that its dose-response predictions are reliable?
Which fitted curve continues linearly beyond 8 mg, and which continues as a cubic polynomial?
Is the piecewise prediction at 10 mg an interpolation or an extrapolation? What caution should accompany it?