Survival Evaluation: Optimize the Partial Chance of the Cox Mannequin | by Nicolo Cosimo Albanese | Dec, 2022

Discovering the coefficients that maximize the log-partial probability in Python

1. Introduction
2. The Cox proportional hazard mannequin
3. Optimization drawback
4. Implementation
5. Conclusions
6. References

Survival evaluation encompasses a set of statistical strategies for describing time to occasion information.

On this submit, we introduce a well-liked survival evaluation algorithm, the Cox proportional hazards model¹. Then, we outline its log-partial probability and the gradient, and optimize it to search out one of the best set of mannequin parameters by a sensible Python instance.

We outline the survival fee as the share of sufferers who haven’t skilled the antagonistic occasion (e.g. demise) after a sure time period.

The Cox proportional hazards mannequin can assess the affiliation between variables and survival fee. Given a set of covariates x, it defines the hazard perform as:

From the method, we observe that the hazard perform h(t|x) is proportional to a baseline hazard perform h₀(t) and relative dangers exp(βx).

The underlying hazard perform h₀(t) doesn’t rely upon the covariates. Because the type of h₀(.) is unspecified, the mannequin is semi-parametric.

Allow us to clarify the which means of the mannequin coefficients by a simplified situation with one covariate solely. Allow us to think about a threat issue xᵢ, for instance smoke, as binary variable (0: non smoker vs. 1: smoker). The Cox mannequin might be expressed as h(t|xᵢ)= h₀(t)exp(βxᵢ), the place exp(β) signifies the relative threat of antagonistic occasion given by smoking over not smoking:

• Danger given by smoking:
  (xᵢ=1): h₀(t)exp(β⋅xᵢ) = h₀(t)exp(β⋅1) = h₀(t)exp(β)
• Danger given by not smoking:
  (xᵢ=0): h₀(t)exp(β⋅xᵢ) = h₀(t)exp(β⋅0) = h₀(t)
• Relative threat = threat given by smoking / threat given by not smoking:
  h₀(t)exp(β) / h₀(t) = exp(β)

The relative threat exp(β) — additionally referred to as hazard ratio — is fixed and doesn’t rely upon time.

In Information Science, the duty of “becoming” a mannequin to a dataset signifies the seek for the set of mannequin parameters that optimizes a sure goal perform. Some widespread examples are the minimization of a loss perform or the maximization of a log-likelihood.

In our case, we have to estimate β with out figuring out h₀(.). To this goal, Cox proposed to maximise the partial likelihood²:

Within the earlier equation:

• Okay is the set of chronologically ordered occasion (demise) instances: t₁ < t₂ < ... <tₖ.
• R(tⱼ) identifies the set of topics in danger at time tⱼ.

Intuitively, the partial chances are a product of the conditional chances of seeing the antagonistic occasions over the set of noticed occasion instances, given the set of sufferers in danger at these instances and beneath the belief of proportional hazards.

We are able to observe that:

• L(β) is impartial from ho(t), that may stay unspecified.
• L(β) doesn’t account for the precise occasion instances, however solely their order.

We might derive the log-partial probability as:

Within the earlier equation:

• N is the variety of topics.
• θ = exp(βx).
• δⱼ signifies the occasion (1: demise, 0: in any other case).

To suit the Cox mannequin, it’s needed to search out the β coefficients that reduce the detrimental log-partial probability.

We recall that the detrimental partial chances are, most often, a strictly convex³ perform. Thus, it has a singular world minimizer.

Allow us to import the wanted libraries:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import reduce
from sksurv.datasets import load_whas500

We load the Worcester Coronary heart Assault Examine dataset⁴ out there within the scikit-survival⁵ package deal. Particularly:

• We concentrate on two covariates:
  – afb: atrial fibrillation (0: no, 1: sure)
  – mitype: MI sort (0: non Q-wave, 1: Q-wave)
• We alter the information to account for ties, i.e. these affected person who skilled the antagonistic occasion on the similar time. Because of the assumption of a steady hazard, the Cox mannequin doesn’t admit ties. For the sake of simplicity, we add a small quantity of random noise to every occasion date to take away them.
• We order the dataset by date, because the partial probability requires ordered occasion instances.

# load the whas500 dataset
X, goal = load_whas500()
# allow us to think about two covariates
cols = ["afb", "mitype"]
df = X[cols].
rename(columns={cols[0]: "v1", 
cols[1]: "v2"}).
astype(float)
# extract occasions and respective instances
df["event"], df["time"] = [list(i) for i in zip(*target)]
# add random noise to the occasion time to keep away from ties
df.time = df.time.apply(lambda x : x + np.random.regular(2, 1))
# type observations by descending occasion time
df = df.sort_values("time", ascending=False).reset_index(drop=True)
# examine first rows
df.head(10)

v = df[["v1", "v2"]].to_numpy()
time, occasion = df.time.to_numpy(), df.occasion.to_numpy().astype(int)

Now we have to outline the target perform for the optimization activity, i.e. the detrimental log-partial probability, which we’re going to reduce:

Observe: in commonplace Machine Studying issues, X typically describes the enter options. In our case, as an alternative, the unknown variable is β, for which we’re looking for one of the best set of values.

def get_theta(x):
'''
Return theta as per detrimental log-partial probability
of the Cox mannequin and its gradient.
It assumes enter options v to be ordered by occasion time.
Args:
- x: beta coefficients 
Output:
- theta_l: cumulative theta <numpy.ndarray>
- theta_l_v: cumulative theta by options <numpy.ndarray>
'''
theta = np.exp(np.dot(v, x))
theta_l = np.cumsum(theta)
theta_l_v = np.cumsum(v * theta.reshape(-1,1), axis=0)
return theta_l, theta_l_v
def objective_function(x):
'''
Return the detrimental log-partial probability 
of the Cox mannequin
Args:
- x: beta coefficients <numpy.ndarray>
Output:
- Unfavorable log-partial probability of the Cox mannequin
'''
theta_l, _ = get_theta(x)
return -np.sum(occasion * (np.dot(v, x) - np.log(theta_l)))

We derive the gradient of the target perform, i.e., its spinoff with respect to β, as follows:

def gradient(x):
'''
Return the gradient of the detrimental log-partial
probability of the Cox mannequin
Args:
- x: beta coefficients <numpy.ndarray>
Output:
- Gradient of the target perform
'''
theta_l, theta_l_v = get_theta(x)
return -np.sum(occasion.reshape(-1,1) * (v-(theta_l_v/theta_l.reshape(-1,1))), axis=0)

We are able to now initialize β and carry out the minimization activity. We use the Newton-CG⁶ algorithm and the scipy package deal:

# beginning values for beta
beta = np.array([1, 1])
opt_result = reduce(
objective_function,
beta, 
methodology = "Newton-CG",
jac = gradient
)
opt_result

The outcomes are:

• β₁ = 0.5293
• β₂ = -0.6541

We are able to match the Cox mannequin on the identical enter information utilizing a library and confirm that we might receive the identical set of values for β:

from sksurv.linear_model import CoxPHSurvivalAnalysis
mannequin = CoxPHSurvivalAnalysis()
model_fit = mannequin.match(
df[["v1", "v2"]], 
np.array(checklist(zip(df.occasion, df.time)), dtype=np.dtype("bool, float")))
model_fit.coef_

Certainly, the β coefficients are virtually an identical, with small variations after the seventh decimal place.

Allow us to plot the estimated optimum and the target perform:

def objective_function_in_x(x0, x1):
'''
Return the detrimental log-partial probability 
of the Cox mannequin evaluated within the given beta
Args:
- x0: enter beta_0 <numpy.ndarray>
- x1: enter beta_1 <numpy.ndarray>
Output:
- goal perform in beta_0, beta_1 <numpy.ndarray>
'''
v0, v1, l = v[:,0], v[:,1], v.form[0]
theta = np.exp(x0*v0 + x1*v1)
return -np.sum(occasion * ((x0*v0 + x1*v1) - np.log(np.cumsum(theta).reshape(-1, l))))
def get_plot_data(inf=-5, sup=5, measurement=10):
'''
Return a three-dim sq. field with the analysis
of the detrimental log-partial probability of the Cox mannequin
Args:
- inf: min worth of the field, default: -5 <int>
- sup: min worth of the field, default: 5 <int>
- measurement: measurement of the output coordinates arrays, default: 10 <int>
Output:
- x0: enter beta_0 <numpy.ndarray>
- x1: enter beta_1 <numpy.ndarray>
- z: goal perform in beta_0, beta_1 <numpy.ndarray>
'''
x0, x1 = np.linspace(inf, sup, measurement), np.linspace(inf, sup, measurement)
x0, x1 = np.meshgrid(x0, x1)
z = np.zeros((measurement, measurement))
for i in vary(0, x0.form[0]):
for j in vary(0, x0.form[1]):
z[i][j] = objective_function_in_x(x0[i][j], x1[i][j])
return x0, x1, z
def get_min_obj_function(mannequin):
'''
Return coordinates of native optimum discovered with optimization
Args:
- mannequin: occasion of <scipy.optimize._optimize.OptimizeResult>
Output:
- x0: optimum for beta_0 <numpy.ndarray>
- x1: optimum for beta_1 <numpy.ndarray>
- z: goal perform within the optimum <numpy.ndarray>
'''
x0, x1 = np.array(mannequin.x[0]), np.array(mannequin.x[1])
z = objective_function_in_x(x0, x1)
return x0, x1, z
# generate information
x0, x1, z = get_plot_data(-5, 5, 10)
x0_min, x1_min, z_min = get_min_obj_function(opt_result)
# plot the target perform and the native optimum
fig = plt.determine(figsize=plt.figaspect(0.4))
# left subplot
ax = fig.add_subplot(1, 2, 1, projection="3d")
ax.contour3D(x0, x1, z, 100, cmap="viridis")
ax.scatter(x0_min, x1_min, z_min, marker="o", colour="crimson", linewidth=50000)
ax.set_xlabel("$β_1$")
ax.set_ylabel("$β_2$")
ax.set_zlabel("- Log-Partial Chance")
# proper subplot
ax = fig.add_subplot(1, 2, 2, projection="3d")
ax.contour3D(x0, x1, z, 100, cmap="viridis")
ax.scatter(x0_min, x1_min, z_min, marker="o", colour="crimson", linewidth=50000)
ax.set_xlabel("$β_1$")
ax.set_ylabel("$β_2$")
ax.set_zlabel("- Log-Partial Chance")
ax.view_init(10, 30)
fig.suptitle("Unfavorable log-partial probability of the Cox mannequin with native optimum", fontsize=10);

Observe: the optimization drawback with the beforehand outlined capabilities might be solved with any variety of enter variables v. Nonetheless, the 3-D plot requires to think about solely two. In reality, a 3-D plot can solely show one β coefficient for every axis.

From the plot, it’s seen that the detrimental log-partial chances are a convex loss perform.

Within the context of survival evaluation, we launched the Cox proportional hazard mannequin and match it on enter information. Particularly, we wrote the detrimental log-partial probability and its gradient in Python. Then, we minimized it to search out one of the best set of mannequin parameters.

[1] D. R. Cox, Regression Fashions and Life-Tables, Journal of the Royal Statistical Society. Collection B (Methodological), Vol. 34, n°2., pp. 187–220, 1972.

[2] https://en.wikipedia.org/wiki/Likelihood_function

[3] C. M. Winson et al., Fenchel duality of Cox partial probability with an utility in survival kernel studying, Synthetic Intelligence in Medication,
vol. 116, 102077, 2021.

[4] S. Pölsterl, scikit-survival: A Library for Time-to-Occasion Evaluation Constructed on High of scikit-learn, Journal of Machine Studying Analysis, vol. 21, no. 212, pp. 1–6, 2020 (Package deal web site).

[5] https://scikit-survival.readthedocs.io/en/steady/api/generated/sksurv.datasets.load_whas500.html

[6] https://docs.scipy.org/doc/scipy/reference/optimize.minimize-newtoncg.html#optimize-minimize-newtoncg

Observe: the dataset whas500 is freely out there to be used from the scikit-survival package deal. The scikit-survivalpackage deal is licensed beneath the GPL v3.