Нахождение коэффициентов, которые максимизируют частичное логарифмическое правдоподобие в Python
Оглавление
- "Введение"
- Модель пропорциональных рисков Кокса
- Задача оптимизации
- "Выполнение"
- Выводы
- "Использованная литература"
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
Примечание: набор данных whas500
находится в свободном доступе для использования из пакета scikit-survival
. Пакет scikit-survival
находится под лицензией GPL v3.