发帖
雷达卡
另一种方法？
目前流行各种各样的预测模型和方法，为什么我们还需要生存分析这样的方法呢？首先，作为统计学流行的一个分支，生存分析可以在一段连续的时间进行预测。 当其他各类分析预测的方法和模型用来预测事件是否会出现时，生存分析可以用来预测在某一个时间点上事件是否会出现。所以，它需要时间的维度来完成对事件在某一时间发生概率的预测。 有的人觉得生存分析适用的领域在生物医疗科学，比如预测细菌或者细胞的生长死亡期望时间等，但是其实很多其他领域也用生存分析来预测，比如机械的寿命或者维护期限等。
How hard does it get …
直接套用生存分析的方法并不可取也并不会带来好处。 在使用生存分析之前，我们需要明白怎么使用它，为什么使用它。 比如有 Kaplan-Meier Curves, 通过 survival trees 或者 suvival forests 来建立Survival Function。
接下来我们从 R 中一步一步学习一下，如何建立一个有效的生存分析。使用的生存分析工具包为 survival
- # install.packages("survival")
- # Loading the package
- library("survival")
工具包中包括了一个用来使用的数据集。数据集中记录了10年中424位在 Mayo clinic 治疗的患有原发性胆汁性胆管炎的患者的数据。值得注意的是，424位患者中，312人参加了药物D-青霉胺的试验，其余112人希望保守治疗并且进行记录并追踪，但未参加药物试验。 这112例中有6例失踪。我们对数据集中的’时间’和’状态’特征特别感兴趣。 时间表示注册后的天数和最终状态（可以被删减censored，肝脏移植或死亡）间的时间。 有关数据集的更多详细信息可以从命令中读取。
- #Dataset description
- ?pbc
- #Fitting the survival model
- survival_func=survfit(Surv(pbc$time,pbc$status == 2)~1)
- survival_func
- #Plot the survival model
- plot(survival_func)
正如预期的那样，该图向我们展示了随着时间的流逝，生存率的下降的概率。 虚线是上下置信区间。 在的survfit（）函数中，我们将公式表示为～1，这表明我们要求函数仅根据生存对象来拟合模型。 输出和置信区间实际上是Kaplan-Meier估计值。 这一估计在医学研究生存分析中很突出。 Kaplan – Meier估计值是根据治疗后存活一段时间的总人数的患者人数（每位患者为一组数据）得出的。 我们可以用公式表示Kaplan – Meier函数：
- Ŝ(t)=∏(1-di/ni) for all i where ti≤t
- Here, di the number of events and ni is the total number of people at risk at time ti
What to make of the graph?
与其他使用测试样本并对其进行预测的机器学习技术不同，生存分析曲线是一个自我解释的曲线。从曲线上看，处理后约1000天存活的可能性大约为0.8或80％。我们可以类似地定义治疗后不同天数的生存概率。同时，我们也有置信区间范围，表示预期误差的幅度。例如，在存活1000天的情况下，上置信区间达到约0.85或85％，并下降到约0.75或75％。数据范围大约为10年或大约3500天，概率计算非常模糊，所以并不应该被接受。例如，如果想知道治疗后存活4500天的可能性，那么尽管上面的Kaplan – Meier图显示的范围在0.25到0.55之间，这本身就是一个很大的值以适应缺乏数据，但数据仍然不是足够充分，应该使用更好的数据来做出这样的估计。
替代模型：Cox比例危险模型
工具包中还包含了cox比例风险函数coxph（），并使用数据中的其他特征来创建更好的生存模型。 虽然数据没有处理缺失值，但我这里将跳过数据处理环节并直接拟合模型。 但是在具体的实践中，研究者需要研究数据并选择如何恰当地处理数据，以便最佳模型得到拟合。 由于本文的目的是让读者熟悉函数而不是处理函数，所以我跳过处理环节来应用函数。
- Fit Cox Model
- Cox_model = coxph(Surv(pbc$time,pbc$status==2) ~.,data=pbc)
- summary(Cox_model)
coxph(formula = Surv(pbc$time, pbc$status == 2) ~ ., data = pbc)
n= 276, number of events= 111
(142 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
id -2.729e-03 9.973e-01 1.462e-03 -1.866 0.06203 .
trt -1.116e-01 8.944e-01 2.156e-01 -0.518 0.60476
age 3.191e-02 1.032e+00 1.200e-02 2.659 0.00784 **
sexf -3.822e-01 6.824e-01 3.074e-01 -1.243 0.21378
ascites 6.321e-02 1.065e+00 3.874e-01 0.163 0.87038
hepato 6.257e-02 1.065e+00 2.521e-01 0.248 0.80397
spiders 7.594e-02 1.079e+00 2.448e-01 0.310 0.75635
edema 8.860e-01 2.425e+00 4.078e-01 2.173 0.02980 *
bili 8.038e-02 1.084e+00 2.539e-02 3.166 0.00155 **
chol 5.151e-04 1.001e+00 4.409e-04 1.168 0.24272
albumin -8.511e-01 4.270e-01 3.114e-01 -2.733 0.00627 **
copper 2.612e-03 1.003e+00 1.148e-03 2.275 0.02290 *
alk.phos -2.623e-05 1.000e+00 4.206e-05 -0.624 0.53288
ast 4.239e-03 1.004e+00 1.941e-03 2.184 0.02894 *
trig -1.228e-03 9.988e-01 1.334e-03 -0.920 0.35741
platelet 7.272e-04 1.001e+00 1.177e-03 0.618 0.53660
protime 1.895e-01 1.209e+00 1.128e-01 1.680 0.09289 .
stage 4.468e-01 1.563e+00 1.784e-01 2.504 0.01226 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
id 0.9973 1.0027 0.9944 1.000
trt 0.8944 1.1181 0.5862 1.365
age 1.0324 0.9686 1.0084 1.057
sexf 0.6824 1.4655 0.3736 1.246
ascites 1.0653 0.9387 0.4985 2.276
hepato 1.0646 0.9393 0.6495 1.745
spiders 1.0789 0.9269 0.6678 1.743
edema 2.4253 0.4123 1.0907 5.393
bili 1.0837 0.9228 1.0311 1.139
chol 1.0005 0.9995 0.9997 1.001
albumin 0.4270 2.3422 0.2319 0.786
copper 1.0026 0.9974 1.0004 1.005
alk.phos 1.0000 1.0000 0.9999 1.000
ast 1.0042 0.9958 1.0004 1.008
trig 0.9988 1.0012 0.9962 1.001
platelet 1.0007 0.9993 0.9984 1.003
protime 1.2086 0.8274 0.9690 1.508
stage 1.5634 0.6397 1.1020 2.218
Concordance= 0.849 (se = 0.031 )
Rsquare= 0.462 (max possible= 0.981 )
Likelihood ratio test= 171.3 on 18 df, p=0
Wald test = 172.5 on 18 df, p=0
Score (logrank) test = 286.1 on 18 df, p=0
Cox模型输出与线性回归输出相似。 R2只有46％不高，我们没有任何高度重要的特征。 最重要的特征似乎是年龄，胆红素（白细胞）和白蛋白。 让我们看看结果是怎样的。
- #Create a survival curve from the cox model
- Cox_curve <- survfit(Cox_model)
- plot(Cox_curve)
预测：使用生存森林
随机森林也可用于生存分析，R中的ranger包提供了功能。 然而，ranger功能不能处理缺失的值，所以我将使用一个较小的数据，所有行的NA值被删除。 这会将我的数据减少到只有276个观测值。
- #Using the Ranger package for survival analysis
- Install.packages("ranger")
- library(ranger)
让我们看一下随机森林重要性的结果：
- #Get the variable importance
- data.frame(sort(ranger_model$variable.importance,decreasing = TRUE))
- sort.ranger_model.variable.importance..decreasing...TRUE.
Lessons learned: Conclusion
虽然Survival软件包的Kaplan-Meier估计，Cox模型和ranger模型的输入数据都不同，但我们会使用ggplot将它们绘制在同一图上来比较这些方法。
- #Comparing models
- library(ggplot2)
- #Kaplan-Meier curve dataframe
- #Add a row of model name
- km <- rep("Kaplan Meier", length(survival_func$time))
- #Create a dataframe
- km_df <- data.frame(survival_func$time,survival_func$surv,km)
- #Rename the columns so they are same for all dataframes
- names(km_df) <- c("Time","Surv","Model")
- #Cox model curve dataframe
- #Add a row of model name
- cox <- rep("Cox",length(Cox_curve$time))
- #Create a dataframe
- cox_df <- data.frame(Cox_curve$time,Cox_curve$surv,cox)
- #Rename the columns so they are same for all dataframes
- names(cox_df) <- c("Time","Surv","Model")
我们在这里看到，考克斯模型是最具数据量和最大特征波动性的模型。处理过的数据点和处理缺失值的更好数据可能会使我们获得更好的R2和更稳定的曲线。生存与幸存对象有关，因此与机器学习完全不同，它与事件发生有关。了解这种技术对使用更多方法帮助我们解决问题非常重要，因为在这种情况下涉及时间的维度。希望本文能够为您提供生存分析的一瞥以及R中可用的，功能丰富的软件包。
——————————————————————————————————
Author Bio:
This article was contributed by Perceptive Analytics. Madhur Modi, Chaitanya Sagar, Vishnu Reddy and Saneesh Veetil contributed to this article.
Perceptive Analytics provides data analytics, data visualization, business intelligence and reporting services to e-commerce, retail, healthcare and pharmaceutical industries. Our client roster includes Fortune 500 and NYSE listed companies in the USA and India.
京公网安备 11010802022788号
