Nonlinear Covariate Effects in Cause-Specific Cox Models

Why Allow a Nonlinear Covariate Effect?

A cause-specific Cox model can represent how a continuous covariate $x$ relates to the instantaneous rate of a particular event cause without assuming that the relationship is linear. A linear term assumes that a one-unit increase in $x$ has the same multiplicative effect on the cause-specific hazard at every value of $x$. That may be too restrictive: the hazard might rise rapidly at low values and then level off, or be U-shaped. A smooth function can represent such patterns without dividing the covariate into arbitrary groups.

The goal is not to make a curve flexible for its own sake. The fitted shape should be supported by the observed data and should answer a meaningful question about the cause-specific hazard. A smooth covariate effect does not, by itself, describe absolute event probabilities or cumulative incidence.

Definition: For event cause $j$, let $T$ be the time to the first event of any cause, and let $X$ denote the cause. The cause-specific hazard at time $t$, conditional on covariates, is the instantaneous event rate for cause $j$ among individuals who have had no event by $t$. A cause-specific Cox model with a smooth effect of continuous covariate $x$ is $$ \lambda_j(t\mid x,\mathbf z) = \lambda_{j0}(t)\exp\{f(x)+\mathbf z^\mathsf{T}\boldsymbol\beta\}, $$ where $\lambda_{j0}(t)$ is an unspecified baseline hazard for cause $j$, $f(x)$ is the covariate effect, and $\mathbf z$ contains any other covariates with coefficients $\boldsymbol\beta$. A linear effect sets $f(x)=\beta_x x$; a nonlinear model allows a more general smooth $f$.

When fitting the model for cause $j$, an event from another cause removes that individual from the risk set for cause $j$ at that time. This is the cause-specific hazard formulation, not a claim that competing event times are independent. Inference requires an appropriate cause-specific model and non-informative censoring, conditional on the covariates included in the analysis. The proportional-hazards assumption requires the covariate effect to be constant over time: $f(x)$ can be nonlinear in $x$, but it is not a function of $t$ in this model.

Restricted Cubic Splines and the Model

A restricted cubic spline represents $f(x)$ using joined cubic polynomials. The pieces meet smoothly at the knots: the function and its first two derivatives are continuous there. The restriction is that the function is linear below the smallest knot and above the largest knot. That avoids increasingly erratic cubic behaviour at the extremes, although linear tails do not guarantee that estimates in sparsely observed tails are precise.

Definition: A restricted cubic spline with $K$ knots represents a covariate effect as a linear component plus $K-2$ nonlinear components, for $K-1$ covariate degrees of freedom (not counting the model intercept, which is absorbed into the baseline hazard). Its knots include the boundary knots. A common specification uses quantiles of the observed covariate distribution for knot locations; the exact locations and number of knots should be chosen before interpreting the resulting curve.

Writing the spline in basis form makes its role in the Cox model explicit:

f(x)=\beta_1x+\sum_{m=2}^{K-1}\beta_m B_m(x), \qquad \lambda_j(t\mid x,\mathbf z) = \lambda_{j0}(t)\exp\left\{\beta_1x+\sum_{m=2}^{K-1}\beta_mB_m(x)+\mathbf z^\mathsf{T}\boldsymbol\beta\right\}.

The functions $B_m(x)$ are spline basis functions determined by the knot locations and a specified spline parameterization. Their coefficients should not usually be interpreted one at a time: the basis columns are mathematical components, and the scientific interpretation concerns the combined fitted function $f(x)$. Software may scale or parameterize the basis differently, so coefficient values are not comparable across parameterizations unless the basis is the same.

For a comparison at the same time $t$, and with the other covariates held fixed, the cause-specific hazard ratio for $x_1$ versus $x_0$ is

\operatorname{HR}_j(x_1,x_0) = \exp\{f(x_1)-f(x_0)\}.

The difference in fitted effects, not a single spline coefficient, determines this hazard ratio. Its confidence interval should use the estimated covariance among all spline coefficients. A test of the nonlinear terms asks whether the data support curvature beyond a linear effect; it is distinct from a test of whether the covariate has any association with the hazard.

Fitting and Testing Nonlinearity

For a cause-specific model, the spline terms enter the linear predictor just as other covariates do. At distinct event times for cause $j$, the partial likelihood contribution for an event at time $t_i$ is proportional to the event individual's relative hazard divided by the sum of relative hazards among those still at risk for cause $j$. For the model above, that relative hazard is $\exp\{f(x_i)+\mathbf z_i^\mathsf{T}\boldsymbol\beta\}$. Thus, the spline changes how individuals' covariates compare within risk sets; the baseline hazard remains unspecified.

Suppose a linear model and a spline model use the same individuals, outcome definition, other covariates, and risk-set construction. The linear model is nested within the spline model: it sets all $K-2$ nonlinear coefficients to zero. A likelihood-ratio test can compare their partial log likelihoods. Under the null hypothesis of linearity and the usual large-sample regularity conditions, twice the log-likelihood difference is approximately chi-square with $K-2$ degrees of freedom. This test addresses nonlinearity, not whether the overall covariate effect is zero. A separate joint test of all $K-1$ spline coefficients addresses that broader null.

Worked Example: Comparing a Linear Effect With a Spline

In a hypothetical analysis of cause $j$, a linear cause-specific Cox model has partial log likelihood $-218.6$. A model using a four-knot restricted cubic spline for the same continuous covariate, with all other terms unchanged, has partial log likelihood $-214.9$. Assume distinct event times for the stated calculation and that the spline fit is nested in the usual way.

With four knots, the spline has $4-1=3$ covariate degrees of freedom: one linear and two nonlinear. The null hypothesis for the comparison is that the two nonlinear coefficients are both zero. The likelihood-ratio statistic is

2\{\ell_{\text{spline}}-\ell_{\text{linear}}\} = 2\{-214.9-(-218.6)\} = 2(3.7) = 7.4.

The reference distribution is approximately $\chi^2$ with $2$ degrees of freedom, giving $p\approx0.0247$. At the 0.05 significance level, this is evidence against a purely linear effect in this model. It does not establish that the curve is the true functional form, nor does it show that the overall covariate effect is nonzero; that requires a separate joint test. The result also depends on the model and its assumptions being adequate.

Worked Interpretation of a Fitted Curve

The practical quantity of interest is often a contrast between two covariate values. The model's proportional-hazards structure means that this contrast does not vary with time, provided the model assumptions hold. The hazard ratio is conditional on being event-free at the time of comparison and on the other covariates in the model.

Worked Example: Contrasting Two Covariate Values

Suppose a hypothetical fitted cause-specific Cox model uses $x$ measured in years and the effect $f(x)=0.025x-0.40B_2(x)+0.18B_3(x)$. For illustration, the fitted basis values are $B_2(40)=0$, $B_3(40)=0$, $B_2(60)=0.35$, and $B_3(60)=-0.10$. The estimated standard error for the complete contrast $f(60)-f(40)$, calculated using the coefficient covariance matrix, is $0.160$.

First calculate the fitted effect at each value:

$$ f(40)=0.025(40)-0.40(0)+0.18(0)=1.000, $$

$$ f(60)=0.025(60)-0.40(0.35)+0.18(-0.10) =1.500-0.140-0.018 =1.342. $$

The contrast is $1.342-1.000=0.342$, so the estimated cause-specific hazard ratio comparing 60 with 40 years is

\widehat{\operatorname{HR}}_j(60,40) = \exp(0.342) = 1.408.

A 95% Wald confidence interval for the log hazard ratio is $0.342\pm1.96(0.160)$, or $(0.0284,0.6556)$. Exponentiating the endpoints gives an approximate hazard-ratio interval of $(1.029,1.926)$. In this hypothetical model, at a given time and for otherwise identical covariates, the estimated cause-specific hazard for an individual with $x=60$ is 1.408 times that for one with $x=40$. The interval describes uncertainty in this estimated contrast; it is not a probability statement about the fixed parameter.

Choosing Complexity and Checking the Fit

Spline flexibility is a modelling choice, not a guarantee of improved prediction or a correct relationship. More knots allow more curvature but also require more information. Knot locations should lie within the observed covariate range, and the data need adequate support across that range, especially near the boundaries. Report the number and locations of knots, how they were selected, and the scale and range of $x$.

Useful checks include plotting the estimated hazard ratio curve against a clearly stated reference value, with a confidence band, and displaying the covariate distribution or event counts across the same range. Look for abrupt changes unsupported by the data, very wide uncertainty at the extremes, and a curve driven by a few observations. Compare prespecified plausible knot choices or modestly different degrees of freedom to assess sensitivity. A nonlinear covariate term does not fix a violation of proportional hazards: assess proportionality separately, using the cause-specific Cox diagnostics appropriate to the model. Likewise, inspect whether the functional form of other continuous covariates is plausible.

Worked Example: Reading Tail Uncertainty Before Reporting a Curve

A hypothetical analysis uses $x$, measured in milligrams per litre, with five restricted-cubic-spline knots at 12, 24, 38, 55, and 72. The fitted curve is nearly flat from 24 to 55, but only four individuals have $x>65$, and the 95% confidence band widens substantially above 60. The reference value is $x=38$.

Five knots give $5-1=4$ covariate degrees of freedom, including three nonlinear degrees of freedom. Because the boundary knots are 12 and 72, the fitted effect is constrained to be linear below 12 and above 72. That mathematical restriction does not provide strong empirical information for $x>65$: with only four observations in that region, estimates and contrasts there may remain imprecise. A wide band is consistent with that lack of support.

A defensible report would show the curve relative to $x=38$, display the data support, and describe the tail estimates as uncertain rather than interpreting their apparent direction as a stable pattern. As a sensitivity check, the analyst could compare a prespecified lower-complexity spline or a scientifically justified restriction of the reporting range. The model should not be selected solely because one knot choice produces a more striking curve.

Alternatives, Limits, and Common Pitfalls

A linear term is simpler and easier to summarize, but it can miss curvature. Categorizing a continuous variable can make group comparisons easy to present, but loses information, depends on cut points, and imposes discontinuous changes between categories. Restricted cubic splines give a smooth curve with linear tails and are often a useful compromise. Other smooth approaches, such as penalized splines, control flexibility through a penalty rather than a fixed set of knots; their smoothing choices and inference require their own care. These approaches all model the cause-specific hazard and do not turn it into an absolute risk model.

The cause-specific Cox coefficients describe associations in the cause-specific hazard among those still event-free. They are not automatically causal effects, and they do not alone determine the probability of experiencing cause $j$ by a given time. Interpretation must also respect the covariates adjusted for, the event definition, and the model assumptions.

Common mistake: Treating a spline coefficient as the hazard ratio for a one-unit increase in $x$, or reading a significant test of the nonlinear terms as proof of a particular curve. Correction: The effect is the combined function $f(x)$; compute contrasts as $\exp\{f(x_1)-f(x_0)\}$ using all relevant coefficients. A nonlinearity test provides evidence against the specified linear model, not proof that the fitted shape is correct.

Other warning signs include too many spline degrees of freedom for the available event information, sparse observations near the boundary knots, and conclusions based on a curve without uncertainty intervals. A spline can reduce functional-form bias, but it cannot repair incorrect risk sets, unaddressed informative censoring, or a time-varying covariate effect that has been forced to be proportional.

Check Your Understanding

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

With six knots, how many covariate degrees of freedom does a restricted cubic spline have, and how many are nonlinear?
A fitted spline gives $f(x_1)=0.7$ and $f(x_0)=0.2$. What is the estimated cause-specific hazard ratio for $x_1$ versus $x_0$, holding other covariates fixed?
A likelihood-ratio comparison of a linear model with a five-knot spline gives statistic $9.1$. What degrees of freedom are used to test nonlinearity?
Why should a curve with a wide confidence band in a sparsely observed tail be interpreted cautiously, even when the spline is constrained to be linear beyond its boundary knot?

Solutions. (1) Six knots give $6-1=5$ covariate degrees of freedom, of which $6-2=4$ are nonlinear. (2) The contrast is $0.7-0.2=0.5$, so the hazard ratio is $\exp(0.5)=1.649$. (3) A five-knot spline has $5-2=3$ nonlinear degrees of freedom, so the reference distribution is approximately $\chi^2_3$. (4) The linear-tail restriction constrains the shape but does not create data support. With few observations, the estimated tail contrast can be imprecise, as reflected by the wide interval.