Nonlinear Omics Effects On Survival

Why Model a Curved Omics Association?

A continuous omics measurement need not have a constant effect on log hazard across its range. A small increase may matter little at low abundance but more at high abundance, or the association may flatten after a threshold-like region. A linear term in a Cox model cannot represent such curvature: it imposes a constant change in log hazard for each one-unit increase in the measurement.

Restricted cubic splines provide a smooth, flexible representation while avoiding the unstable behavior that unrestricted high-degree polynomials can have at the extremes. The goal is not to discover a biological threshold automatically. It is to estimate a continuous exposure-response curve, quantify its uncertainty, and interpret comparisons against a clearly stated reference measurement.

Model and Spline Definition

Let $T$ denote event time, let $\delta$ indicate whether the event was observed, let $x$ be a continuous omics measurement recorded at baseline, and let $z$ be a vector of adjustment covariates. A Cox model with a spline effect is

h(t\mid x,z)=h_0(t)\exp\{g(x)+z^\mathsf{T}\gamma\}, \qquad g(x)=\beta_x x+\sum_{j=1}^{K-2}\theta_j B_j(x).

Here $h_0(t)$ is an unspecified baseline hazard, $K$ is the number of ordered knots $\kappa_1<\cdots<\kappa_K$, and $B_j(x)$ are nonlinear spline basis functions. The linear term and $K-2$ nonlinear terms together give a restricted cubic spline with $K$ knots. The basis is chosen so the fitted function is linear below the first knot and above the last knot, while it can curve between them.

Definition: For $a\in\mathbb{R}$, define $(x-a)_+=\max(x-a,0)$. One standard restricted cubic spline basis is $$ B_j(x)=(x-\kappa_j)_+^3 -\frac{\kappa_K-\kappa_j}{\kappa_K-\kappa_{K-1}}(x-\kappa_{K-1})_+^3 +\frac{\kappa_{K-1}-\kappa_j}{\kappa_K-\kappa_{K-1}}(x-\kappa_K)_+^3, \quad j=1,\ldots,K-2. $$ Together with $x$, these basis functions span the restricted cubic spline. Rescaling the basis changes coefficient magnitudes but not the fitted curve when the model is transformed consistently.

The model assumes proportional hazards: for fixed $x$ and $z$, the hazard ratio comparing two individuals does not vary with time. It also assumes the specified covariate effects are adequately represented, observations are independent across individuals, and censoring is non-informative conditional on the modeled covariates for the likelihood analysis. Event times and censoring must be recorded in a way appropriate for the chosen Cox-model implementation. A spline changes the functional form of $x$'s effect; it does not remove these assumptions or address time-varying effects.

Contrasts Relative to a Reference

A spline coefficient by itself is generally not a useful effect summary. Its value depends on the basis and knot locations. Instead, choose a reference measurement $x_0$, such as the median or a clinically interpretable value, and compare the fitted log hazards at $x$ and $x_0$, holding the adjustment covariates fixed.

\log\{\mathrm{HR}(x,x_0)\}=g(x)-g(x_0), \qquad \mathrm{HR}(x,x_0)=\exp\{g(x)-g(x_0)\}.

This is a conditional hazard ratio, not a ratio of survival probabilities. For example, an estimated hazard ratio of 1.30 means that at a given time, among individuals with the same modeled adjustment covariates, the fitted instantaneous event hazard at $x$ is 30% higher than at $x_0$, under the proportional-hazards model. It does not say that the probability of the event by a specified time is 30% higher.

Worked Example: Calculating a Hazard Ratio From a Spline Curve

Suppose a hypothetical analysis uses an omics measurement $x$ in assay units, with four knots at 0, 1, 2, and 3. For this basis, the two nonlinear functions are $B_1(x)=(x-0)_+^3-3(x-2)_+^3+2(x-3)_+^3$ and $B_2(x)=(x-1)_+^3-2(x-2)_+^3+(x-3)_+^3$. Assume the fitted curve is $g(x)=0.20x-0.015B_1(x)+0.040B_2(x)$. Compare $x=2.5$ with the reference $x_0=1$, at the same values of all other covariates.

At $x=2.5$, the positive-part terms give $B_1(2.5)=2.5^3-3(0.5^3)=15.625-0.375=15.25$ and $B_2(2.5)=1.5^3-2(0.5^3)=3.375-0.250=3.125$. At the reference, $B_1(1)=1$ and $B_2(1)=0$. Therefore,

$$ g(2.5)-g(1) =[0.20(2.5)-0.015(15.25)+0.040(3.125)] -[0.20(1)-0.015(1)+0.040(0)] =0.39625-0.185=0.21125. $$

The fitted hazard ratio is $\exp(0.21125)=1.235$. In this hypothetical model, the fitted hazard at 2.5 assay units is about 23.5% higher than at 1 assay unit, conditional on the adjustment covariates. This is an estimated association, not evidence that changing the measurement would cause a change in hazard.

Uncertainty for the Fitted Contrast

Let $\eta$ contain the estimated spline and linear coefficients, and let $c(x,x_0)$ be the vector of differences between the corresponding basis values at $x$ and $x_0$. Then the estimated log hazard ratio is $c(x,x_0)^\mathsf{T}\hat{\eta}$. If $\widehat{\mathrm{Var}}(\hat{\eta})$ is the estimated coefficient covariance matrix, the model-based standard error is

SE\{\widehat{\log\mathrm{HR}}(x,x_0)\} =\sqrt{c(x,x_0)^\mathsf{T}\widehat{\mathrm{Var}}(\hat{\eta})c(x,x_0)}.

A large-sample 95% confidence interval on the log-hazard-ratio scale is the estimate plus or minus $1.96$ standard errors. Exponentiating both endpoints gives an interval for the hazard ratio. This interval describes the uncertainty of the estimated contrast under the model; it is not a 95% probability that the fixed true hazard ratio lies inside the particular interval.

Worked Example: Confidence Interval for a Reference Contrast

Continue with the hypothetical model and contrast above. The coefficient contrast is $c=(2.5-1,\;B_1(2.5)-B_1(1),\;B_2(2.5)-B_2(1))=(1.5,14.25,3.125)$. Suppose the estimated covariance matrix of the three coefficients is diagonal with variances $0.0100,\;0.00000400,\;0.0000160$, respectively. For this illustration, assume the off-diagonal covariances are zero.

The estimated variance of the log hazard ratio is

$$ 1.5^2(0.0100)+14.25^2(0.00000400)+3.125^2(0.0000160) =0.022500+0.00081225+0.00015625 =0.0234685. $$

Thus the standard error is $\sqrt{0.0234685}=0.1532$. The 95% interval for the log hazard ratio is $0.21125\pm1.96(0.1532)=[-0.0890,\;0.5115]$. Exponentiating gives a hazard-ratio interval of $[\exp(-0.0890),\exp(0.5115)]=[0.915,\;1.668]$. Because the interval includes 1, these data do not establish a difference in hazard at the 5% two-sided level for this particular contrast. The point estimate remains 1.235; the interval indicates substantial uncertainty about its magnitude and direction.

Why the Restricted Form Works

The restriction to linear tails follows from how the basis is constructed. Below the first knot, all positive-part terms are zero, so each $B_j(x)=0$. Above the last knot, all positive-part terms are active. The coefficients multiplying the cubic terms then sum to zero:

1-\frac{\kappa_K-\kappa_j}{\kappa_K-\kappa_{K-1}} +\frac{\kappa_{K-1}-\kappa_j}{\kappa_K-\kappa_{K-1}}=0.

The quadratic terms also cancel after expanding each cube. Consequently, beyond the last knot each $B_j(x)$ is at most linear in $x$, rather than growing cubically. Since the model also includes $x$, the fitted log hazard is linear in either tail. This restriction reduces the risk of extreme extrapolated curvature while leaving the curve flexible within the knot range; it does not guarantee that tail estimates are precise.

Worked Example: Testing for Added Nonlinearity

Suppose a hypothetical study fits two Cox models to the same individuals, with the same adjustment covariates and the same fixed four-knot placement. The reduced model has a single linear term for the omics measurement. The larger model uses the restricted cubic spline, adding two nonlinear coefficients. Assume both models are unpenalized and satisfy the regularity conditions for a likelihood-ratio comparison. Their partial log-likelihoods are $-120.0$ for the reduced model and $-117.0$ for the spline model.

The likelihood-ratio statistic is twice the difference in partial log-likelihoods: $2[-117.0-(-120.0)]=6.00$. Under the null hypothesis that both added nonlinear coefficients are zero, this is approximately distributed as a chi-square random variable with 2 degrees of freedom. For 2 degrees of freedom, the upper-tail probability at 6.00 is $\exp(-6.00/2)=0.0498$. Thus the comparison gives borderline evidence, at the 5% level, that the spline improves fit over a linear effect in this specified model. It does not identify a threshold, establish causality, or test the proportional-hazards assumption.

Fitting Choices, Diagnostics, and Alternatives

Knot placement and number control flexibility. Knots are often placed at selected quantiles of the observed measurement distribution, with boundary knots near the extremes and interior knots spread through the data. Three to five knots are common starting points, not universal rules. More knots add flexibility and parameters, which can make estimates unstable when the number of events is small or measurements are concentrated in a narrow range. The measurement scale must be stated: a log-transformed abundance has a different interpretation and curve than untransformed abundance.

Inspect the fitted curve and its uncertainty across the observed measurement range, and note how many observations and events support each region. A sharp bend driven by a few extreme measurements, very wide intervals in a tail, or large changes after modest knot adjustments are warnings of weak support. Also check whether the proportional-hazards assumption is plausible; a curved effect across $x$ and a time-varying effect across $t$ are different issues. If the curve changes over time, a time-varying-effect model is needed rather than merely adding more spline knots for $x$.

The main alternatives answer different modeling needs. A linear Cox term is simpler and more precise when linearity is adequate, but it can miss curvature. Categorizing the measurement can make group summaries easy to present, but discards within-group information and makes results depend on cut points. Fractional polynomials offer a different low-dimensional smooth functional form, but do not have the same restricted-tail construction. These alternatives do not remove the need to check model assumptions or support in the data.

Approach	Useful feature	Main caution
Linear Cox effect	Simple, compact interpretation	Imposes a constant log-hazard change per unit
Restricted cubic spline	Smooth curve with linear tails	Needs enough events and careful knot choices
Categorized measurement	Can communicate group comparisons	Information loss and cut-point sensitivity
Fractional polynomial	Flexible parametric curve	Different functional-form and extrapolation behavior

In multi-omics work, a spline for one measurement is not a substitute for addressing correlated features, feature selection, or a clear adjustment strategy. If the feature or knot locations were selected after examining outcomes, ordinary confidence intervals and tests may understate the uncertainty induced by that selection. Likewise, the usual large-sample inference can be unreliable with few events, separation-like patterns, or sparse observations near the tails.

Common Mistake

A frequent error is to interpret a spline coefficient as the hazard ratio for a one-unit increase in the omics measurement. That is generally wrong: the nonlinear basis coefficients are tied to the chosen basis, and the effect depends on the starting and ending values of $x$. The correction is to compute $g(x)-g(x_0)$, state the reference and units, and exponentiate the full contrast. A second mistake is to call a curved fitted association a time-varying effect. Curvature means the association varies with the omics value; a time-varying effect means it varies over follow-up time.

Key takeaway: A restricted cubic spline models a potentially nonlinear association between a continuous omics measurement and log hazard, with linear behavior beyond its boundary knots. Interpret the fitted function through hazard-ratio contrasts against an explicit reference, and assess uncertainty, event support, and proportional hazards separately.

Check Your Understanding

Work through the questions, then compare your reasoning with the solutions below.

What does a restricted cubic spline assume about the fitted log-hazard curve beyond its boundary knots?
If $\hat g(x)-\hat g(x_0)=\log(0.80)$, what is the fitted hazard ratio, and how should it be described?
Why does testing that all nonlinear spline coefficients equal zero compare a spline effect with a linear effect?
Name one diagnostic warning that the spline curve is poorly supported by the data, and explain why adding knots may not fix it.

Solutions. (1) The curve is linear in each tail beyond the boundary knots; it is not forced to be flat. (2) The fitted hazard ratio is $\exp\{\log(0.80)\}=0.80$, meaning the fitted hazard at $x$ is 20% lower than at $x_0$, conditional on the modeled covariates and under proportional hazards. (3) The spline model contains a linear term plus nonlinear basis terms; setting all nonlinear coefficients to zero leaves only the linear effect. (4) A sharp tail bend supported by very few events is a warning, because the tail shape is weakly identified; adding knots can increase instability rather than create information in that region.