High-Dimensional Cox Regression With Multi-Omics

Why High-Dimensional Cox Regression Is Different

A multi-omics study may measure thousands of features on a few hundred patients, while only a subset experience the event during follow-up. In that setting, fitting an ordinary Cox model with one coefficient per feature is not simply a matter of asking software to handle a large table. The data may not contain enough independent information to identify those coefficients, and small changes in the observations can produce large changes in the fitted estimates.

The central issue is the information available in the observed event-risk comparisons. The number of patients matters, but so do the number of events, the composition of the risk sets, and the degree of redundancy among the omics features. A model can therefore be unstable even when the number of predictors is smaller than the number of patients.

Definition: In a p-greater-than-n setting, the number \(p\) of candidate predictors exceeds the number \(n\) of observed individuals. In a Cox model, the number of observed events and the variation in predictors across risk sets also constrain how much information is available about the coefficient vector.

The Cox Model and Its Information

Let \(Y_i\) be the observed follow-up time for individual \(i\), let \(\delta_i=1\) indicate an observed event and \(\delta_i=0\) indicate censoring, and let \(x_i\) be a vector of \(p\) covariates, including any omics features in the model. The proportional hazards model specifies

h_i(t\mid x_i)=h_0(t)\exp(x_i^\mathsf{T}\beta),

where \(h_0(t)\) is an unspecified baseline hazard and \(\beta\) is the coefficient vector. For the standard partial-likelihood expression below, assume independent individuals, proportional hazards, covariates measured as specified by the model, and censoring independent of event time conditional on the modeled covariates. Also assume distinct event times for this expression; with tied event times, a specified tie-handling approximation or exact method is needed.

At event time \(t_j\), let \(i(j)\) denote the individual who experiences the event, and let \(R_j\) be the risk set: those still under observation and event-free just before \(t_j\). The partial likelihood compares the event individual with everyone in that risk set:

L(\beta)=\prod_{j=1}^{E} \frac{\exp(x_{i(j)}^\mathsf{T}\beta)} {\sum_{k\in R_j}\exp(x_k^\mathsf{T}\beta)}, \qquad \ell(\beta)=\log L(\beta),

Here \(E=\sum_i\delta_i\) is the number of observed events. The exponential quantity \(\exp(x_k^\mathsf{T}\beta)\) is the individual’s relative hazard under the model. Multiplying every such quantity in a risk set by the same positive constant leaves that risk-set factor unchanged. Equivalently, adding a common constant to every linear predictor \(x_k^\mathsf{T}\beta\) in that risk set leaves its factor unchanged.

Key result: At any finite \(\beta\), the contribution of event \(j\) to the score is the event individual’s covariate vector minus the risk-set weighted mean. The observed information is a sum of weighted risk-set covariances. A coefficient direction with no variation in these risk-set comparisons has no information in the partial likelihood.

To see this, define the risk-set weights

w_{kj}(\beta)= \frac{\exp(x_k^\mathsf{T}\beta)} {\sum_{l\in R_j}\exp(x_l^\mathsf{T}\beta)}.

These weights are nonnegative and sum to one over \(R_j\). Differentiating the event contribution \(x_{i(j)}^\mathsf{T}\beta-\log\{\sum_{k\in R_j}\exp(x_k^\mathsf{T}\beta)\}\) gives the score contribution

U_j(\beta)=x_{i(j)}-\sum_{k\in R_j}w_{kj}(\beta)x_k.

Differentiating again gives the negative weighted covariance of the covariates in that risk set. Thus the observed information, the negative Hessian of \(\ell(\beta)\), is

I(\beta)=\sum_{j=1}^{E}\left[ \sum_{k\in R_j}w_{kj}x_kx_k^\mathsf{T} - \left(\sum_{k\in R_j}w_{kj}x_k\right) \left(\sum_{k\in R_j}w_{kj}x_k\right)^\mathsf{T} \right].

Each bracketed term is a covariance matrix and is positive semidefinite. Its rank is at most \(|R_j|-1\), because a weighted covariance cannot contain more independent directions than the number of risk-set members minus one. Consequently, the rank of \(I(\beta)\) is at most \(\min\left(p,\sum_j\min(p,|R_j|-1)\right)\). The event count alone does not give a rank-at-most-\(E\) bound, since a single event contribution can have rank greater than one.

Three Ways Estimates Become Unreliable

Worked Example: More Coefficients Than Event Comparisons

Suppose a hypothetical study has \(n=6\) patients, \(E=4\) observed events, and \(p=8\) omics predictors. Even if the predictors vary, what can be said about the information matrix?

For each event, its risk-set covariance has rank no greater than \(|R_j|-1\), and no greater than \(p\). Adding four such matrices gives

\operatorname{rank}\{I(\beta)\} \leq \min\left(p,\sum_{j=1}^{4}\min(p,|R_j|-1)\right).

Because there are six patients, all risk-set covariate contrasts lie in a subspace of dimension at most five, so the information matrix has rank at most 5. Thus the \(8\times8\) matrix is singular, with at least three coefficient directions having zero curvature at that \(\beta\). Therefore, the usual inverse-information calculation cannot provide a finite standard error for all eight coefficients. This calculation does not say which particular features are uninformative; that depends on their values in the risk sets.

Exact non-identifiability occurs when different coefficient vectors produce the same risk-set comparisons and hence the same partial likelihood. One source is an exact linear dependence among predictor columns. Another is a predictor that is constant across all individuals: its contribution adds the same amount to every linear predictor in each risk set and cancels from the partial likelihood. These facts mean that \(p>n\) is a strong warning, but the risk-set contrasts—not merely the raw predictor count—determine the identifiable directions.

Worked Example: Duplicate Omics Features

Suppose two measured features are exact duplicates across all patients. To illustrate the information calculation, take three risk-set covariate contrasts to be \(z_1=(1,1)\), \(z_2=(2,2)\), and \(z_3=(-1,-1)\). For a two-person risk set at \(\beta=0\), the two members have equal weights, so its information contribution is one quarter of the contrast outer product. The resulting information is

I(0)=\frac{1}{4}\sum_{j=1}^{3}z_jz_j^\mathsf{T} =\frac{1}{4}(1+4+1) \begin{pmatrix}1&1\\1&1\end{pmatrix} =1.5\begin{pmatrix}1&1\\1&1\end{pmatrix}.

This matrix has eigenvalue \(3\) in the direction \((1,1)\) and eigenvalue \(0\) in the direction \((1,-1)\). The likelihood can learn about the sum \(\beta_1+\beta_2\), but not how that sum is divided between the two duplicate features. For example, \((\beta_1,\beta_2)=(2,0)\) and \((0,2)\) give the same linear predictors and partial likelihood. A numerical optimizer may return one solution, but that does not make the individual coefficients identifiable.

Near-duplicate features create a related but less absolute problem. The information matrix may be invertible yet have a very small eigenvalue. In that weakly informed direction, modest changes in the sample, censoring pattern, or preprocessing can cause large coefficient changes. This is why a model can return finite estimates and still be unreliable.

Worked Example: Monotone Likelihood and a Coefficient That Diverges

Consider a one-feature Cox partial likelihood with three event contributions, each comparing an event individual with a single other individual. Suppose the event individual has feature value \(1\), and the other risk-set member has value \(0\), at each of the three event times. Each factor is \(\exp(\beta)/(1+\exp(\beta))\), so

L(\beta)=\left\{\frac{\exp(\beta)}{1+\exp(\beta)}\right\}^{3}.

At \(\beta=0\), \(L(0)=(1/2)^3=0.125\). At \(\beta=\log(3)\), each factor is \(3/(1+3)=0.75\), so \(L(\log 3)=0.75^3=0.421875\). As \(\beta\) tends to infinity, every factor tends to 1 and \(L(\beta)\) tends to 1, but no finite \(\beta\) attains that value. The likelihood increases without reaching a finite maximum. This is called monotone likelihood; the ordinary finite maximum-likelihood estimate does not exist. In practice, a fitting algorithm may report a very large coefficient, a large standard error, or failure to converge.

Assumptions, Diagnostics, and Alternatives

High dimension does not replace the Cox model’s assumptions. The proportional-hazards structure must be an adequate description for the modeled effects; censoring must be handled under an appropriate conditional independence assumption; and the event and risk-set data must be correctly recorded. The earlier material on proportional-hazards diagnostics concerns time patterns in effects: a time-varying effect is a model misspecification, not a cure for rank deficiency. Likewise, adding more omics measurements does not create more independent event information.

Useful warning signs include a singular or nearly singular information matrix, unusually large coefficient standard errors, coefficients that change greatly when a few patients are removed, extreme fitted linear predictors, and convergence warnings or estimates that keep growing as the iteration limit is increased. Check for duplicated or nearly duplicated features, features with little variation, and predictor patterns that nearly separate event individuals from other members of their risk sets. Compare the number of events with the number of candidate coefficients, but do not treat any single events-per-variable threshold as a universal guarantee of stability.

The approaches to high dimensionality differ in what they do to the model:

Approach	What it changes	Important caution
Unpenalized Cox regression	Estimates all specified coefficients by maximizing the partial likelihood.	Requires enough information for identifiable, stable estimates; a finite optimizer output is not proof of either.
Prespecified predictor reduction	Fits a smaller model based on scientific considerations or design.	Data-driven reduction can introduce selection bias and must be performed within the evaluation procedure when prediction is assessed.
Dimension reduction	Represents many features using fewer derived predictors.	The resulting effects may be harder to interpret at the individual-feature level; any outcome-guided construction must be accounted for in evaluation.
Penalized Cox regression	Constrains or shrinks coefficient estimates to stabilize fitting.	The penalty changes the estimation target and does not remove the need for appropriate validation and tuning.

The last three are broad strategies, not interchangeable fixes. For prediction, all outcome-dependent preprocessing and tuning must be conducted without using held-out evaluation outcomes; the earlier material on leakage audits and nested evaluation is relevant. For scientific inference, a stable predictive score does not by itself establish precise, unbiased effects for individual omics features.

Common Pitfalls and Practical Interpretation

“There are more patients than predictors, so the fit must be stable.” The event count and the independent variation within risk sets also matter; weak directions can remain even if \(p<n\).
“The software returned coefficients, so they are uniquely estimated.” An optimizer can return one point on a flat ridge or stop at a large finite value when the likelihood has no finite maximum.
“A small p-value for one feature settles the issue.” Standard errors and tests rely on the fitted model and its information. Rank deficiency, near-collinearity, or data-driven selection can make ordinary unpenalized inference misleading.
“A constant can be added to all exponential-scale risk scores without changing a factor.” The invariant operation is multiplication of those scores by a common positive constant. Equivalently, one may add a common constant to all linear predictors in the risk set.

When the information is weak or singular, individual unpenalized coefficients should not be interpreted as precise feature effects. The model may identify some combinations of features more clearly than others; that does not make each coefficient in the combination separately meaningful. Report the fitting problem and the diagnostics, and choose a modeling strategy that matches whether the objective is prediction, estimation, or feature-level interpretation.

Check Your Understanding

Answer each question, then compare with the solutions below.

Why is the observed information matrix in a Cox partial likelihood built from risk-set covariances?
A model has \(p=10\) coefficients and \(E=6\) observed events. What rank bound or conclusion about singularity follows from the event count alone?
Two features are exact duplicates. Which coefficient combination can be informed by the partial likelihood, and which direction is not identifiable?
In the monotone-likelihood example, what happens to the likelihood as \(\beta\) tends to infinity, and why is there no finite maximum?

Solutions. (1) Differentiating each event’s log partial-likelihood contribution gives the event covariate minus its weighted risk-set mean; differentiating again gives the negative weighted covariance. Summing the negative Hessians gives the information matrix. (2) The event count alone gives no rank-at-most-six bound, so it does not establish that the \(10\times10\) information matrix is singular. The risk-set covariance ranks and covariate variation are needed to determine whether ten independent coefficient directions can be identified. (3) The sum \(\beta_1+\beta_2\) can be informed; the difference direction \((1,-1)\) is not identifiable because changing the coefficients in opposite directions leaves their combined predictor unchanged. (4) Each factor tends to 1, so the product tends to 1, but it is less than 1 for every finite \(\beta\); the likelihood has a supremum rather than a finite attained maximum.