如何用pymc3（威布尔分布回归）进行简单的生存分析？

Question

我是 pymc3 的新手，我读过《黑客贝叶斯方法》，并尽最大努力完成 pymc3 中现有的生存分析教程。但是，我不明白如何编写/解释“生存函数”。

对于这个问题，我从 NIST 定义的威布尔分布中生成了一些虚拟数据：

n = 1000
alpha = 1 
gam = 0.5
mu = 0
noise = np.random.normal(0, 0.025, [n, 1])

x = np.random.rand(n, 1)*10
f_x = (gam/alpha)*(((gam-mu)/alpha)**(gam-1))*np.exp(-((x-mu)/alpha)**gam)
y = f_x + noise

由于我想创建一个包含审查数据和未经审查数据的模型，例如贝叶斯参数化的 pymc3教程，因此我实现了截止并审查了这些数据点：

cens = np.array([1 if k < 7.5 else 0 for k in x])

然后我开始根据先验建立我的模型：

with pm.Model() as survival_model:

     alpha0 = pm.Normal('alpha0', mu=1, sigma = 1)
     gam0 = pm.Normal('gam0', mu=0.5, sigma = 1)
     mu0 = pm.Normal('mu0', mu=0.0, sigma = 1)
     noise0 = pm.Normal('noise0', mu=0.0, sigma = 0.05)

现在麻烦开始了，我知道我需要定义一个考虑删失值的似然函数，并将所有参数作为输入来输出似然度。我认为对于截尾值，我必须找到描述 P(Y > y) 的方程。通常，我可以使用 CDF，但在本例中，我发现使用 Matlab 和 Mathematica 不定积分不存在。我应该怎么办？

Answer 1

您可以在 pymc3 中使用自定义分布创建似然函数。特别是pm.DensityDist班级。我将使用 surpyval 预先存在的函数方法。

这要求输入是对数似然。对于观测值，对数似然是密度函数的对数。对于审查观测值，对数似然实际上只是分布的累积风险函数的负数。将所有这些值相加并返回。

完整示例如下：

import pymc3 as pm
from surpyval import Weibull

# Create 100 random variables with alpha=50 and beta=2
rvs = Weibull.random(100, 50, 2)

# Set all values above 60 to be 60.. 
# i.e all above 60 are censored
rvs[rvs >= 60] = 60

# separate censored and observed
censored = rvs[rvs == 60]
failures = rvs[rvs < 60]

with pm.Model() as survival_model:
     alpha = pm.Normal('alpha', mu=100, sigma = 10)
     beta = pm.Gamma('beta', alpha=5, beta = 2)
     
     # Where the magic happens
     def logp(failures, censored):
          # likelihood for observed is the log of the density function
          failed_log_like = Weibull.log_df(failures, alpha, beta)
          # The log-likelihood for censored observations is the negative
          # of the cumulative hazard function.
          censored_log_like = -Weibull.Hf(censored, alpha, beta)
          # Return the sum of all of it
          return censored_log_like.sum() + failed_log_like.sum()

     weibull_neg_ll = pm.DensityDist('weibull_neg_ll', 
                                     logp,
                                     observed={'failures' : failures,
                                               'censored' : censored})
     start = pm.find_MAP()
     step = pm.NUTS()
     trace = pm.sample(10000, step, start, random_seed=123, progressbar=True)
     pm.traceplot(trace)