Как получить функцию плотности вероятности из CoxPHSurvivalAnalysis в scikit-survival?

Функцию плотности вероятности (PDF) можно получить из кумулятивной функции распределения (CDF) как:

f(t) = dF(t)/dt

Теперь в анализе выживания (SA) PDF (f(t)) можно выразить через функцию выживания S(t) и функцию риска h(t), которая определяется следующим образом:

f(t) = h(t) x S(t)

где S(t) = 1 - F(t) и h(t) = -dS(t)/dt x S(t) = dH(t)/dt

Итак, PDF f(t) можно выразить как: f(t) = dH(t)/dt x S(t)

Теперь, чтобы вычислить функцию риска f(t), нам нужна производная кумулятивной функции риска (CHF) H(t). Поскольку все CHF являются дискретными точками данных, нам нужен InterpolatedUnivariateSpline из библиотеки scipy, чтобы дифференцировать их. Он создает плавную сплайн-интерполяцию CHF, которую затем можно дифференцировать, чтобы получить h(t). Вот небольшая модификация вставленного кода:

# Import the necessary libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sksurv.linear_model import CoxPHSurvivalAnalysis
from scipy.interpolate import InterpolatedUnivariateSpline

# Define a function to compute the probability density function (pdf) 
# from the cumulative hazard function (chf) and survival function (sf).
def compute_pdf_from_chf_and_sf(chf, sf):
    # The hazard function is the derivative of the cumulative hazard function.
    # We use InterpolatedUnivariateSpline for spline interpolation to create a smooth 
    # function approximation of the CHF. This provides us with a smooth curve that 
    # passes through each data point, allowing us to differentiate the function and obtain 
    # the hazard function.
    chf_spline = InterpolatedUnivariateSpline(chf.x, chf(chf.x))
    hazard_function = chf_spline.derivative()(chf.x)
    
    # The pdf can be computed using the formula: pdf(t) = hazard(t) * survival(t)
    pdf = hazard_function * sf(chf.x)
    return chf.x, pdf

# Generate random data for demonstration purposes
# Here, we create a random dataset with one covariate and survival times.

np.random.seed(42)  # Setting a fixed seed.
data = np.random.randint(5, 30, size=10)
X_train = pd.DataFrame(data, columns=['covariate'])
y_train = np.array(np.random.randint(0, 100, size=10)/100, dtype=[('status', bool), ('target', float)])

# Initialize and fit the Cox Proportional Hazards model
estimator = CoxPHSurvivalAnalysis()
estimator.fit(X_train, y_train)

# Predict for new data points
X_test = pd.DataFrame({'covariate': [12, 2]})
cumulative_hazard_functions = estimator.predict_cumulative_hazard_function(X_test)
survival_functions = estimator.predict_survival_function(X_test)

# Plot the Cumulative Hazard, Survival, and PDF side by side in a single row
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for chf, sf in zip(cumulative_hazard_functions, survival_functions):
    # Compute the pdf using our defined function
    times, pdf_values = compute_pdf_from_chf_and_sf(chf, sf)
    
    # Plotting the cumulative hazard function
    axes[0].step(chf.x, chf(chf.x), where='post')
    
    # Plotting the survival function
    axes[1].step(sf.x, sf(sf.x), where='post')
    
    # Plotting the probability density function
    axes[2].step(times, pdf_values, where='post')

# Setting titles for each subplot
axes[0].set_title('Cumulative Hazard Functions')
axes[1].set_title('Survival Functions')
axes[2].set_title('Probability Density Functions')

# Display the plots
plt.tight_layout()
plt.show()

что приводит к

Ссылки : Машинное обучение для анализа выживания: опрос