Нахождение коэффициентов, которые максимизируют частичное логарифмическое правдоподобие в Python

Оглавление

  1. "Введение"
  2. Модель пропорциональных рисков Кокса
  3. Задача оптимизации
  4. "Выполнение"
  5. Выводы
  6. "Использованная литература"

1. Введение

Анализ выживания включает в себя набор статистических методов для описания времени до данных о событии.

В этом посте мы представляем популярный алгоритм анализа выживания — модель пропорциональных рисков Кокса¹. Затем мы определяем его логарифмически-частную вероятность и градиент и оптимизируем его, чтобы найти наилучший набор параметров модели на практическом примере Python.

2. Модель пропорциональных рисков Кокса

Мы определяем коэффициент выживаемости как процент пациентов, у которых не возникло неблагоприятного события (например, смерти) через определенный период времени.

Модель пропорциональных рисков Кокса может оценить связь между переменными и выживаемостью. Учитывая набор ковариат x, он определяет функцию опасности как:

Из формулы видно, что функция риска h(t|x) пропорциональна базовой функции риска h₀(t) и относительным рискам exp(βx).

Основная функция опасности h₀(t) не зависит от ковариат. Поскольку форма h₀(.) не указана, модель является полупараметрической.

Поясним значение коэффициентов модели с помощью упрощенного сценария только с одной ковариантой. Рассмотрим фактор риска xᵢ, например курение, как двоичную переменную (0: некурящий против 1: курящий). Модель Кокса может быть выражена как h(t|xᵢ)= h₀(t)exp(βxᵢ), где exp(β) указывает на относительный риск неблагоприятного события, связанного с курением, по сравнению с отказом от курения:

  • Риск, связанный с курением:
    (xᵢ=1): h₀(t)exp(β⋅xᵢ) = h₀(t)exp(β⋅1) = h₀(t)exp(β)
  • Риск, связанный с отказом от курения:
    (xᵢ=0): h₀(t)exp(β⋅xᵢ) = h₀(t)exp(β⋅0) = h₀(t)
  • Относительный риск = риск, связанный с курением / риск, связанный с отказом от курения:
    h₀(t)exp(β) / h₀(t) = exp(β)

Относительный риск exp(β) — также называемый коэффициентом риска — постоянен и не зависит от времени.

3. Проблема оптимизации

В науке о данных задача «подгонки» модели под набор данных означает поиск набора параметров модели, оптимизирующего определенную целевую функцию. Некоторыми распространенными примерами являются минимизация функции потерь или максимизация логарифмического правдоподобия.

В нашем случае нам нужно оценить β, не зная h₀(.). С этой целью Кокс предложил максимизировать частичное правдоподобие²:

В предыдущем уравнении:

  • K представляет собой набор хронологически упорядоченных моментов времени события (смерти): t₁ < t₂ < ... <tₖ.
  • R(tⱼ) определяет группу субъектов, подверженных риску в момент времени tⱼ.

Интуитивно, частичная вероятность является произведением условных вероятностей наблюдения неблагоприятных событий в течение набора времени наблюдаемых событий, учитывая набор пациентов, подверженных риску в эти моменты времени, и при допущении пропорциональных опасностей.

Мы можем заметить, что:

  • L(β) не зависит от ho(t), который может оставаться неуказанным.
  • L(β) учитывает не фактическое время событий, а только их порядок.

Мы могли бы получить частную логарифмическую вероятность как:

В предыдущем уравнении:

  • N - количество предметов.
  • θ = exp(βx).
  • δⱼ указывает событие (1: смерть, 0: иначе).

Чтобы соответствовать модели Кокса, необходимо найти коэффициенты β, которые минимизируют отрицательное логарифмическое частичное правдоподобие.

Напомним, что отрицательное частичное правдоподобие в большинстве случаев является строго выпуклой³ функцией. Таким образом, он имеет уникальный глобальный минимизатор.

4. Реализация

Давайте импортируем необходимые библиотеки:

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

Мы загружаем набор данных Worcester Heart Attack Study⁴, доступный в пакете scikit-survival⁵. В частности:

  • Мы фокусируемся на двух ковариатах:
    - afb: фибрилляция предсердий (0: нет, 1: да)
    - mitype: тип ИМ (0: отсутствие зубца Q, 1: зубец Q)
  • Мы корректируем данные для учета связей, т. е. тех пациентов, которые испытали нежелательное явление в одно и то же время. Из-за предположения о непрерывной опасности модель Кокса не допускает ничьих. Для простоты мы добавляем небольшое количество случайного шума к каждой дате события, чтобы удалить их.
  • Мы упорядочиваем набор данных по дате, так как частичная вероятность требует упорядоченного времени события.
# load the whas500 dataset
X, target = load_whas500()

# let us consider two covariates
cols = ["afb", "mitype"]

df = X[cols].\
        rename(columns={cols[0]: "v1", 
                        cols[1]: "v2"}).\
        astype(float)

# extract events and respective times
df["event"], df["time"] = [list(i) for i in zip(*target)]

# add random noise to the event time to avoid ties
df.time = df.time.apply(lambda x : x + np.random.normal(2, 1))

# sort observations by descending event time
df = df.sort_values("time", ascending=False).reset_index(drop=True)

# inspect first rows
df.head(10)

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

Теперь нам нужно определить целевую функцию для задачи оптимизации, то есть отрицательноелогарифмическое частичное правдоподобие, которое мы собираемся минимизировать:

Примечание. В стандартных задачах машинного обучения X обычно описывает входные функции. В нашем случае вместо этого неизвестная переменная — β, для которой мы пытаемся найти наилучший набор значений.

def get_theta(x):
    '''
    Return theta as per negative log-partial likelihood
    of the Cox model and its gradient.
    It assumes input features v to be ordered by event time.
    
    Args:
        - x: beta coefficients 
    
    Output:
        - theta_l: cumulative theta <numpy.ndarray>
        - theta_l_v: cumulative theta by features <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 negative log-partial likelihood 
    of the Cox model
    
    Args:
        - x: beta coefficients <numpy.ndarray>
    Output:
        - Negative log-partial likelihood of the Cox model
    '''
    theta_l, _ = get_theta(x)
    return -np.sum(event * (np.dot(v, x) - np.log(theta_l)))

Мы получаем градиент целевой функции, т. е. ее производную по β, следующим образом:

def gradient(x):
    '''
    Return the gradient of the negative log-partial
    likelihood of the Cox model
    
    Args:
        - x: beta coefficients <numpy.ndarray>
    Output:
        - Gradient of the objective function
    '''
    theta_l, theta_l_v = get_theta(x)
    return -np.sum(event.reshape(-1,1) * (v-(theta_l_v/theta_l.reshape(-1,1))), axis=0)

Теперь мы можем инициализировать β и выполнить задачу минимизации. Мы используем алгоритм Newton-CG⁶ и пакет scipy:

# starting values for beta
beta = np.array([1, 1])

opt_result = minimize(
    objective_function,
    beta, 
    method = "Newton-CG",
    jac = gradient
)

opt_result

Результаты:

  • β₁ = 0.5293
  • β₂ = -0.6541

Мы можем подогнать модель Кокса к тем же входным данным, используя библиотеку, и убедиться, что мы получим тот же набор значений для β:

from sksurv.linear_model import CoxPHSurvivalAnalysis

model = CoxPHSurvivalAnalysis()
model_fit = model.fit(
    df[["v1", "v2"]], 
    np.array(list(zip(df.event, df.time)), dtype=np.dtype("bool, float")))

model_fit.coef_

Действительно, коэффициенты β почти идентичны, с небольшими отличиями после седьмого знака после запятой.

Построим расчетный оптимум и целевую функцию:

def objective_function_in_x(x0, x1):
    '''
    Return the negative log-partial likelihood 
    of the Cox model evaluated in the given beta
    
    Args:
        - x0: input beta_0 <numpy.ndarray>
        - x1: input beta_1 <numpy.ndarray>
    Output:
        - objective function in beta_0, beta_1 <numpy.ndarray>
    '''
    v0, v1, l = v[:,0], v[:,1], v.shape[0]
    theta = np.exp(x0*v0 + x1*v1)
    return -np.sum(event * ((x0*v0 + x1*v1) - np.log(np.cumsum(theta).reshape(-1, l))))

def get_plot_data(inf=-5, sup=5, size=10):
    '''
    Return a three-dim square box with the evaluation
    of the negative log-partial likelihood of the Cox model
    
    Args:
      - inf: min value of the box, default: -5 <int>
      - sup: min value of the box, default: 5 <int>
      - size: size of the output coordinates arrays, default: 10 <int>
    Output:
      - x0: input beta_0 <numpy.ndarray>
      - x1: input beta_1 <numpy.ndarray>
      - z: objective function in beta_0, beta_1 <numpy.ndarray>
    '''
    x0, x1 = np.linspace(inf, sup, size), np.linspace(inf, sup, size)
    x0, x1 = np.meshgrid(x0, x1)
    z = np.zeros((size, size))
    for i in range(0, x0.shape[0]):
        for j in range(0, x0.shape[1]):
            z[i][j] = objective_function_in_x(x0[i][j], x1[i][j])
    return x0, x1, z

def get_min_obj_function(model):
    '''
    Return coordinates of local optimum found with optimization
    
    Args:
      - model: instance of <scipy.optimize._optimize.OptimizeResult>
    Output:
      - x0: optimum for beta_0 <numpy.ndarray>
      - x1: optimum for beta_1 <numpy.ndarray>
      - z: objective function in the optimum <numpy.ndarray>
    '''
    x0, x1 = np.array(model.x[0]), np.array(model.x[1])
    z = objective_function_in_x(x0, x1)
    return x0, x1, z

# generate data
x0, x1, z = get_plot_data(-5, 5, 10)
x0_min, x1_min, z_min = get_min_obj_function(opt_result)

# plot the objective function and the local optimum
fig = plt.figure(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", color="red", linewidth=50000)
ax.set_xlabel("$β_1$")
ax.set_ylabel("$β_2$")
ax.set_zlabel("- Log-Partial Likelihood")

# right 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", color="red", linewidth=50000)
ax.set_xlabel("$β_1$")
ax.set_ylabel("$β_2$")
ax.set_zlabel("- Log-Partial Likelihood")
ax.view_init(10, 30)
fig.suptitle("Negative log-partial likelihood of the Cox model with local optimum", fontsize=10);

Примечание: задача оптимизации с ранее определенными функциями может быть решена с любым количеством входных переменных v. Тем не менее, трехмерный график требует рассмотрения только двух. Фактически, трехмерный график может отображать только один коэффициент β для каждой оси.

Из графика видно, что отрицательное логарифмическое частичное правдоподобие представляет собой выпуклую функцию потерь.

5. Выводы

В контексте анализа выживаемости мы представили модель пропорциональных рисков Кокса и подогнали ее к входным данным. В частности, мы написали отрицательное логарифмическое частичное правдоподобие и его градиент на Python. Затем мы минимизировали его, чтобы найти наилучший набор параметров модели.

6. Ссылки

[1] Д. Р. Кокс, Регрессионные модели и таблицы продолжительности жизни, Журнал Королевского статистического общества. Серия Б (Методическая), Том. 34, № 2., стр. 187–220, 1972.

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

[3] К. М. Винсон и др., Двойственность Фенхеля частичного правдоподобия Кокса с применением в обучении ядра выживания, Искусственный интеллект в медицине,
vol. 116, 102077, 2021.

[4] С. Пёльстерл, scikit-survival: библиотека для анализа времени до события, построенная на основе scikit-learn, Journal of Machine Learning Research, vol. 21, нет. 212, стр. 1–6, 2020 (Сайт пакета).

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

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

Примечание: набор данных whas500 находится в свободном доступе для использования из пакета scikit-survival. Пакет scikit-survival находится под лицензией GPL v3.