paranet：具有弹性网络正则化的参数生存模型

该paranet软件包允许拟合具有右删失的参数生存模型，该模型具有 L1 和 L2 正则化惩罚（弹性网络），python具有sklearn类似的语法。目前支持三种参数分布：

1. 指数的
2. 威布尔
3. 冈佩尔茨

之所以选择这些分布，是因为在实证研究中的使用频率以及它们可以建模的危害分布范围。此外，可以使用逆向方法廉价地提取样本，从而允许分位数和随机数据生成方法在很短的时间内运行。

X 和 Y 包目前支持参数模型，但这些包都不支持正则化。弹性净生存模型可以与 Z 拟合，但这仅适用于 Cox-PH 模型。虽然 Cox 模型是生存建模的一个非常重要的工具，但它对大规模数据集的主要限制是 i) 它无法轻松推断个体生存时间，并且 ii) 它的损失函数是非参数的，并且具有运行时间增长 O(n^2) 而不是 O(n)。

该paranet软件包允许用户在右删失数据上拟合高维线性模型，然后提供对事件发生时间结果的个性化或分组预测。在客户流失数据上拟合参数模型可以让数据科学回答有趣的问题，例如：“在这 100 个客户中，我们什么时候首先预计其中 10% 会流失？”或“对于这个新客户，他们在什么时间点离开我们的风险最高（即最大危险）？”，或“对于现有客户，他们在 10 个月后流失的概率是多少？”。

(0) 安装

该paranet软件包可以使用pip install paranet=0.1. 注意这个包已经用 python 3.9+ 测试过了。使用较早版本的 python 可能会导致错误。

(1) 基本语法

该类parametric是这个包的济贫院模型。初始化模型时，用户总是需要指定dist参数。这可以是有效分布类型的列表或字符串。这个列表的长度没有限制，但如果它是 length k，那么后续的时间测量将需要一个向量或一个带有k列的矩阵。

尽管每个输入参数都在文档字符串中定义，但有几个参数会在整个过程中频繁出现，并且为了方便起见在此处定义。

1. x：(n,p)协变量数组。用户可以手动添加截距和缩放数据，也可以设置add_int=True和scale_x=True.
2. t: 一(n,k)组时间测量值，应该是非零的。如果 $k\geq 0$ 则模型假定每一列对应于（可能）不同的分布。
3. d：一(n,k)组审查指标，其值应为 0 或 1。按照惯例，0 对应于经过审查的观测值，1 对应于未经审查的观测值。
4. gamma：正则化的强度（参见第（2）节进一步描述）。如果这个变量是一个向量或者一个矩阵，它必须对应于的列数x（如果给定的话包括截距）。
5. rho：相对 L1/L2 强度，其中 0 对应于 L2-only（即 Ridge 回归），1 对应于 L1-only（即 Lasso）。

作为一项规则 paranet，将尝试在可能的情况下进行广播。例如，如果时间测量数组是t.shape==(n,k)，dist=='weibull'那么它将假设 的每一列t都是 Weibull 分布。相反，如果t.shape==(n,)和，它会为每个分发dist=['weibull','gompertz']广播副本。t

该类parametric具有 X 键方法。如果x、t或d已初始化，则参数可以留空。

1. fit(x, t, d, gamma, rho)：将适合给定gamma/rho惩罚的弹性网络模型，并启用类似hazard执行的方法。
2. find_lambda_max(x, t, d, gamma)：使用次梯度的 KKT 条件来确定将gamma除尺度和形状参数之外的所有协变量清零所需的最大值。
3. {hazard,survival,pdf}(t, x)：与predictis类似sklearn，这些方法提供了对危险、生存和密度函数的个性化估计。
4. quantile(percentile, x)：提供个体化生存分布的分位数。
5. rvs(n_sim, censoring)：为审查目标生成一定数量的样本。

初始化parametric类时，用户可以包含数据矩阵，这些矩阵将被保存以供以后需要它们的方法使用。但是，在后面的方法中指定这些参数将始终覆盖（而不是替换）这些继承的属性。

1. dist：必需的参数，它是一个字符串或一个列表，其元素必须是以下之一：指数、weibull 或 gompertz。
2. alpha：可以提前手动定义形状参数（需要匹配 的维度dist）。
3. beta: 可以预先定义尺度参数（需要与 和 的维数相匹配dist）x。
4. scale_x: 将协变量标准化，使其均值为 0，方差为 1。在使用任何正则化时，强烈建议这样做。如果此参数设置为 True，请始终提供协变量的原始形式，因为它们将在推理期间进行缩放。
5. scale_t：将时间向量归一化为模型拟合的最大值，这有助于解决溢出问题。在推理过程中，输出将始终返回到原始比例。但是，系数会因此而改变。

(2) 玩具示例

下面的代码块显示了如何将三个参数分布拟合到由协变量生成的单个数据数组中。有关与 jupyter notebook 一起使用的其他演示，paranet请参阅示例文件夹。

# Load modules
import numpy as np
import pandas as pd
import plotnine as pn
from scipy import stats
from paranet.models import parametric

# (i) Create a toy dataset
n, p, seed = 100, 5, 3
x = stats.norm().rvs([n,p],seed)
shape = 2
b0 = 0.25
beta = stats.norm(scale=0.5).rvs([p,1],seed)
eta = x.dot(beta).flatten() + b0
scale = np.exp(eta)
t = (-np.log(stats.uniform().rvs(n,seed))/scale)**(1/shape)
d = np.ones(n)

# (ii) Fit the (unregularized) model
mdl = parametric(dist=['exponential', 'weibull', 'gompertz'], x=x, t=t, d=d, scale_x=False, scale_t=False)
mdl.fit()

# (iii) Plot the individual survival, hazard, and density functions for five "new" observations
n_points = 100
n_new = 4
t_range = np.exp(np.linspace(np.log(0.25), np.log(t.max()), n_points))
x_new = stats.norm().rvs([n_new,p],seed)

# We can at look at the hazard for first out-of-sample individual
# Notice that for the exponential distribution (first column) the hazard is independent of time which is as expected
print(np.log(mdl.hazard(t_range, np.tile(x_new[[0]],[n_points,1]))).round(2))

# We can then comprehensively calculate this for each method
methods = ['hazard', 'survival', 'pdf']
holder = []
for j in range(n_new):
    x_j = np.tile(x_new[[j]],[n_points,1])
    for method in methods:
        res_j = getattr(mdl, method)(t_range, x_j)
        if method == 'hazard':
            res_j = np.log(res_j)
        res_j = pd.DataFrame(res_j, columns = mdl.dist).assign(time=t_range,method=method, sample=j+1)
        holder.append(res_j)

# Plot the results
res = pd.concat(holder).melt(['sample','time','method'],None,'dist')

gg_res = (pn.ggplot(res, pn.aes(x='time', y='value', color='dist')) + 
    pn.theme_bw() + pn.geom_line() + 
    pn.scale_color_discrete(name='Distribution') + 
    pn.facet_grid('method~sample',scales='free',labeller=pn.labeller(sample=pn.label_both)))

(3) 概率分布参数化

每个参数生存分布由尺度 $\lambda$ 和除指数分布外的形状 $\alpha$ 参数定义。每个分布都已经过参数化，因此规模参数的值越高，表示“风险”越高。密度函数如下所示。比例和形状参数也必须为正，除了形状参数可以为正或负的 Gompertz 分布的情况。

$$ \begin{align*} f(t;\lambda, \alpha) &= \begin{cases} \lambda \exp{ -\lambda t } & \text{ if Exponential} \ \alpha \lambda t^{ \alpha-1} \exp{ -\lambda t^{\alpha} } & \text{ if Weibull} \ \lambda \exp{ \alpha t } \exp{ -\frac{\lambda}{\alpha}( e^{\alpha t} - 1) } & \text{ if Gompertz} \ \end{cases} \end{align*} $$

当从单变量分布到多变量分布时，我们假设尺度参数采用的是参数（$\eta$）的线性组合的指数变换（以确保正性）。通过平衡数据似然性与系数的大小 ($R$) 来进行优化，如下所示。

$$ \begin{align*} \lambda_i &= \exp\Big( \beta_0 + \sum_{j=1}^p x_{ij}\beta_j \Big) \ R(\beta;\gamma,\rho) &= \gamma\big(\rho | \beta_{1:} | 1 + 0.5(1-\rho)|\beta {1:}| 2^2\big) \ \ell(\alpha,\beta, \gamma,\rho) &= \begin{cases} -n^{-1}\sum {i=1}^n\delta_i\log\lambda_i - \lambda_i t_i + R(\beta;\gamma,\rho ) & \text{ 如果指数} \ -n^{-1}\sum_{i=1}^n\delta_i[\log(\alpha\lambda_i)+(\alpha-1)\log t_i] - \lambda t_i^\alpha + R(\beta;\gamma,\rho) & \text{ if Weibull} \ -n^{-1}\sum_{i=1}^n\delta_i[\log\lambda + \alpha t] - \frac{\lambda}{\alpha}(\exp{\alpha t_i } -1) + R(\beta;\gamma,\rho) & \text{ if Gompertz} \ \end{cases} \结束{对齐*} $$

(4) 审查是如何计算的？

调用该parametric.rvs方法时，用户可以指定删失值。在paranet中，审查是由指数分布产生的，其值小于实际值。正式地：

$$ \begin{align*} T^{\text{obs}} &= \begin{cases} T^{\text{act}} & \text{ if } T^{\text{act}} < C \ C & \text{ 否则} \end{cases} \ C &\sim \text{Exp}(\lambda_C) \ \end{align*} $$

当然还有其他可能产生审查的过程（例如EXAMPLE）。在删失过程中使用指数分布的原因是为了解决一个（相对）简单的优化问题，即找到单个尺度参数 ($\lambda_C$)，从而获得 $\phi$ 的（渐近）删失概率：

$$ \begin{align*} \phi(\lambda_C) &= P(C \leq T_i) = \int_0^\infty f_T(u) F_C(u) du, \ \lambda_C^* &= \arg\min_ λ | \phi(\lambda) - \phi^* |_2^2, \end{align*} $$

其中$F_C(u)$ 是以$\lambda_C$ 作为尺度参数的指数分布的CDF，$f_T(u)$ 是目标分布的密度（例如 Weibull-pdf）。找到尺度参数相当于一个可以用 执行的求根问题scipy。对于多变量情况，找到单个尺度参数更加复杂，因为需要对 $\lambda_i$ 本身的分布做出假设，这是随机的。虽然很容易生成审查特定的分布（即$C_i$），但这将打破非信息审查假设，因为审查随机变量现在是已实现风险评分的函数。这paranet包假设协变量来自标准正态分布：$x_{ij} \sim N(0,1)$ 使得 $\eta_i \sim N(0, |\beta|^2_2)$ 和 $\lambda_i \sim \text{对数正态}(0, |\beta|^2_2)$。重要的是，至少要对数据进行标准化，以使这一假设可信。

$$ \begin{align*} P(C \leq T) &= \int_0^\infty \Bigg( \int_0^\infty P(C \leq T_i) di \Bigg) F_C(u) du \ &= \ int_0^\infty\int_0^\infty F_C(u)f_{i}(u) f_\lambda(i) du di , \end{align*} $$

其中 $f_{i}(u)$ 是在 $u$ 处评估的目标分布的密度，而 $f_\lambda(i)$ 是在 $i$ 处评估的对数正态分布的 pdf。这是一个要解决的更复杂的积分，paranet目前在对值网格进行积分时使用蛮力方法，而不是使用双求积，因为后一种方法在运行时间方面被证明是非常昂贵的。

(5) 优化是如何发生的？

与 不同glmnet的是，这些paranet包不使用坐标下降 (CD)。相反，这个包使用 L1 范数的平滑近似来允许直接优化，scipy如下所示。glmnet由于存在形状参数，参数生存模型不容易适应 使用的迭代重加权最小二乘 (IRLS) 方法。总之，指数模型可以很容易地利用现有的基于 CD 的弹性网络求解器来拟合。转向近端梯度下降将能够直接优化 L1 范数损失，并代表一个可能的未来版本。

$$ \begin{align*} R(\beta;\gamma,\rho) &= \gamma\big(\rho | \beta_{1:} | 1 + 0.5(1-\rho)|\beta {1 :}|_2^2\big) \ |\beta| &\approx \sqrt{\beta^2 + \epsilon} \ \frac{\partial R}{\partial \beta} &\approx \gamma\Bigg(\rho \frac{\beta}{\sqrt{\beta ^2+\epsilon}} + (1-\rho)\beta\Bigg) \end{align*} $$

(6) 作出贡献

如果你有兴趣为这个包做出贡献，请随时给我发电子邮件或提出拉取请求。该包的主要类和功能已经过重要的单元测试，为确保更改不会破坏包，建议source tests/run_pytest_{univariate,multivariate,elnet}.sh在进行任何最终合并之前运行。这个包是用特定的 conda 环境构建的，开发人员可以使用conda env create -f paranet.yml，然后conda activate paranet.