数据科学与机器学习案例之信用卡欺诈识别(严重类失衡数据建模)
数据科学与机器学习案例之Stacking集成方法对鸢尾花进行分类
数据科学案例之生存分析与二手车定价
生存分析
生存分析的基本概念
生存分析是一种将事件的结果与出现此类结果所经历的时间结合在一起的统计分析方法。生存分析包含了:事件、生存时间、删失数据、生存函数、风险函数、半衰期等概念。
事件:
起始事件与终点事件。比如我们在研究二手车分析定价
中起始事件是二手汽车平台发布车辆出售信息,终点事件就是用户下单购买二手车。
生存时间:
生存时间不关注事件发生的瞬时时间,只关注事件的起始事件开始到终点事件发生之间的时间间隔。生存时间服从偏态分布:指数分布、 威布尔分布,对数logistic分布。
删失数据:
是指在观察时间窗内无法获取事件发生时间,其中包括右删失、左删失与区间删失。右删失是指知道观察的起始时间,但无法获取终点事件发生的具体时间,主要由以下的原因导致:研究对象在观察时间窗内未发生有效事件、研究对象在观察时间窗内由于某些原因被丢失、研究对象在事件发生之前由于非事件原因脱离有效观测。
生存函数:
对于观察事件窗内的任意时刻,生存函数是到该时刻仍未发生该事件的概率。生存函数是每个时刻生存概率的乘积, S ( t ) = p ( T > t ) = 1 − F ( t ) S(t)=p(T>t)=1-F(t) S(t)=p(T>t)=1−F(t),其中 F ( t ) F(t) F(t)代表生存时间的累积分布函数,表示事件发生时间未超过时刻 t t t的概率。生存函数具有以下三个性质: S ( t ) ∈ [ 0 , 1 ] S(t)\in[0,1] S(t)∈[0,1],且 S ( t ) S(t) S(t)单调递减;在起始时刻 t = 0 t=0 t=0时,所有对象均处于存活状态,此时 S ( t ) = 1 S(t)=1 S(t)=1;当时间 t t t趋于无穷大时,生存概率趋近于0, S ( t ) = 0 S(t)=0 S(t)=0.
风险函数与累积风险函数:
风险函数又称条件死亡率,是指在时刻 t t t之前未发生任何事件而恰好在时刻 t t t发生事件的概率。风险函数的公式为: h ( t ) = lim h → 0 p ( T < t + h ∣ T ≥ t ) h h(t)=\lim_{h \to 0} \frac{p(T<t+h | T\geq t)}{h} h(t)=limh→0hp(T<t+h∣T≥t),当生存时间 t t t为连续型随机变量时,风险函数表示为: h ( t ) = f ( t ) S ( t ) = − d l n [ S ( t ) ] d t h(t)=\frac{f(t)}{S(t)}=\frac{-dln[S(t)]}{dt} h(t)=S(t)f(t)=dt−dln[S(t)],累积风险函数表示为: H ( t ) = ∫ 0 t h ( u ) d u = − l n [ S ( t ) ] H(t)=\int_0^t {h(u)}du=-ln[S(t)] H(t)=∫0th(u)du=−ln[S(t)],生存函数与风险函数的关系表达为: S ( t ) = e − H ( t ) = e − ∫ 0 t h ( u ) d u S(t)=e^{-H(t)}=e^{-\int_0^t {h(u)} du} S(t)=e−H(t)=e−∫0th(u)du.
中位生存时间:
指一半的个体未发生终点事件的时间。生存分析中由于删失数据的存在,无法直接以平均时间反应时间发生的时间水平,因此生存分析中常以存活概率 50 % 50\% 50%时对应消耗的时间描述。
生存分析的统计分析方法
对生存函数的刻画:
KM曲线法
,假设 t 1 、 t 2 、 t 3 、,,, t k t_{1}、t_{2}、t_{3}、,,,t_{k} t1、t2、t3、,,,tk共有 k k k个观测时刻及 N N N个样本, a j a_{j} aj代表在时刻 t j t_{j} tj发生事件的人数, b j b_{j} bj是在 ( t j , t j + 1 ) (t_{j},t_{j+1}) (tj,tj+1)时间段内删失的人数,则在 t j t_{j} tj时刻处于风险中的人数表示为: c j = ( a j + b j ) + ( a j + 1 + b j + 1 ) + . . . + ( a k + b k ) c_{j}=(a_{j}+b_{j})+(a_{j+1}+b_{j+1})+...+(a_{k}+b_{k}) cj=(aj+bj)+(aj+1+bj+1)+...+(ak+bk),
则 t j t_{j} tj时刻的风险率为 d j n j \frac{d_{j}}{n_{j}} njdj,KM曲线
对应生存函数的估计为: S ( t ) = ∏ j : t j ≤ t n j − d j n j S(t)=\displaystyle\prod_{j:t_{j}\leq t}^{} \frac{n_{j}-d_{j}}{n_{j}} S(t)=j:tj≤t∏njnj−dj.
通过直接观察我们无法确定生存曲线之间是否具有显著性差异,引入了log-rank与Breslow检验方法
。
Nelson-Aalen累计风险曲线:
H ( t ) = ∏ j : t j ≤ t d j n j H(t)=\displaystyle\prod_{j:t_{j}\leq t}^{} \frac{d_{j}}{n_{j}} H(t)=j:tj≤t∏njdj.
生存函数的参数估计:
参数模型是假定生存时间是符合某种分布,通常情形下生存时间的分布主要有指数分布、威布尔分布、对数正态分布和对数logistic分布。下面只是对指数分布做简要介绍:指数分布的生存函数为: S ( t ) = e − λ t S(t)=e^{-\lambda t} S(t)=e−λt,风险函数表示为: h ( t ) = λ h(t)=\lambda h(t)=λ.
半参COX等比例风险回归模型
解决问题:
KM曲线
只能在单一的分类变量间进行比较,也就是说KM曲线
只描述了该单变量与生存时间之间的关系而忽略了其他变量的影响,由此产生了半参COX等比例风险回归模型
.
半参COX等比例风险回归模型公式为:
h ( t ; x ) = h 0 ( t ) e β x h(t;x)=h_{0}(t)e^{\beta x} h(t;x)=h0(t)eβx,模型中的 h 0 ( t ) h_0(t) h0(t)为非参数部分; ( β 1 , β 2 , β 3 , . . . β p ) (\beta_{1},\beta_{2},\beta_{3},... \beta_{p}) (β1,β2,β3,...βp)为模型的参数。 h 0 ( t ) h_0(t) h0(t)是所有自变量取值为 0 0 0的风险率。每个自变量对生存率的影响不随时间的改变而改变。
对于模型的解释:当自变量 x p x_{p} xp为连续性变量时,回归系数 β p \beta_{p} βp是除了 x p x_{p} xp之外其他变量均不变的情况下, x p x_{p} xp每增加一个单位,风险率变化 e β p e^{\beta_{p}} eβp倍;当自变量为分类型变量时,假设 x p x_{p} xp的类别有两类: 0 0 0与 1 1 1,回归系数 β j \beta_{j} βj表示除了 x p x_{p} xp之外其他变量均不变的情况下,类别 1 1 1相对于类别 0 0 0的风险率是 e β p e^{\beta_{p}} eβp倍.
二手车分析定价
生存曲线的绘制以及差异对比
绘制二手车销售情况随时间变化的整体生存曲线与不同颜色下的生存曲线
生存数据的分析主要是比较不同的处理方法以及各种因素(变化)对生存函数的影响,而不是单纯的寻找拟合的模型。
> library(survival)
> df = read.csv('data.csv')
> fit = survfit(Surv(在售时长,是否卖出) ~ 颜色,data = df)
> summary(fit)
Call: survfit(formula = Surv(在售时长, 是否卖出) ~ 颜色, data = df)
颜色=Black
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 1755 46 0.9738 0.00381 0.9663 0.981
1 1665 43 0.9486 0.00530 0.9383 0.959
2 1578 38 0.9258 0.00634 0.9135 0.938
3 1507 36 0.9037 0.00718 0.8897 0.918
4 1436 28 0.8861 0.00777 0.8710 0.901
5 1366 27 0.8685 0.00832 0.8524 0.885
6 1317 29 0.8494 0.00886 0.8322 0.867
7 1271 23 0.8340 0.00926 0.8161 0.852
8 1222 23 0.8184 0.00965 0.7997 0.837
9 1180 29 0.7982 0.01011 0.7787 0.818
10 1122 13 0.7890 0.01031 0.7690 0.809
11 1096 17 0.7768 0.01057 0.7563 0.798
12 1061 12 0.7680 0.01075 0.7472 0.789
13 1031 10 0.7605 0.01090 0.7394 0.782
14 1010 14 0.7500 0.01111 0.7285 0.772
15 987 12 0.7409 0.01128 0.7191 0.763
16 969 10 0.7332 0.01142 0.7112 0.756
17 950 4 0.7301 0.01148 0.7080 0.753
18 944 7 0.7247 0.01157 0.7024 0.748
19 931 5 0.7208 0.01164 0.6984 0.744
20 919 4 0.7177 0.01170 0.6951 0.741
21 913 4 0.7145 0.01175 0.6919 0.738
22 906 4 0.7114 0.01180 0.6886 0.735
23 897 4 0.7082 0.01186 0.6853 0.732
24 890 4 0.7050 0.01191 0.6821 0.729
25 883 3 0.7026 0.01195 0.6796 0.726
26 877 2 0.7010 0.01198 0.6779 0.725
27 874 7 0.6954 0.01207 0.6722 0.719
28 860 2 0.6938 0.01209 0.6705 0.718
29 853 4 0.6905 0.01214 0.6671 0.715
30 846 36 0.6612 0.01258 0.6370 0.686
31 770 35 0.6311 0.01299 0.6062 0.657
32 699 28 0.6058 0.01332 0.5803 0.633
33 638 21 0.5859 0.01357 0.5599 0.613
34 580 37 0.5485 0.01403 0.5217 0.577
35 511 29 0.5174 0.01437 0.4900 0.546
36 444 34 0.4778 0.01479 0.4496 0.508
37 382 26 0.4452 0.01510 0.4166 0.476
38 333 26 0.4105 0.01538 0.3814 0.442
39 287 32 0.3647 0.01565 0.3353 0.397
40 240 16 0.3404 0.01574 0.3109 0.373
41 205 14 0.3172 0.01585 0.2876 0.350
42 172 11 0.2969 0.01597 0.2672 0.330
43 150 15 0.2672 0.01611 0.2374 0.301
44 122 6 0.2540 0.01618 0.2242 0.288
45 107 8 0.2350 0.01631 0.2052 0.269
46 89 9 0.2113 0.01647 0.1813 0.246
47 65 16 0.1593 0.01678 0.1296 0.196
48 41 8 0.1282 0.01672 0.0993 0.166
49 17 7 0.0754 0.01819 0.0470 0.121
颜色=Brown
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 200 6 0.970 0.0121 0.9466 0.994
1 190 7 0.934 0.0176 0.9003 0.969
2 178 4 0.913 0.0201 0.8747 0.954
3 170 4 0.892 0.0223 0.8491 0.937
4 165 3 0.876 0.0238 0.8301 0.923
5 160 3 0.859 0.0252 0.8112 0.910
6 155 7 0.820 0.0280 0.7673 0.877
7 146 5 0.792 0.0297 0.7361 0.853
8 140 3 0.775 0.0306 0.7175 0.838
9 135 2 0.764 0.0313 0.7049 0.828
10 133 1 0.758 0.0315 0.6987 0.822
11 131 1 0.752 0.0318 0.6924 0.817
12 126 1 0.746 0.0321 0.6859 0.812
13 122 1 0.740 0.0324 0.6792 0.807
14 119 2 0.728 0.0331 0.6657 0.796
17 115 1 0.721 0.0334 0.6589 0.790
19 111 2 0.708 0.0340 0.6448 0.778
20 108 3 0.689 0.0349 0.6236 0.761
33 100 1 0.682 0.0352 0.6161 0.755
50 99 2 0.668 0.0359 0.6014 0.742
51 91 6 0.624 0.0377 0.5543 0.703
52 77 2 0.608 0.0385 0.5369 0.688
53 67 2 0.590 0.0394 0.5173 0.672
54 63 3 0.562 0.0407 0.4872 0.647
55 54 2 0.541 0.0418 0.4648 0.629
56 51 2 0.520 0.0427 0.4422 0.610
57 47 2 0.497 0.0437 0.4188 0.591
58 41 5 0.437 0.0460 0.3553 0.537
59 32 2 0.410 0.0470 0.3270 0.513
60 26 1 0.394 0.0478 0.3104 0.499
61 24 1 0.377 0.0485 0.2933 0.486
65 18 1 0.356 0.0502 0.2705 0.470
66 17 1 0.335 0.0514 0.2484 0.453
67 15 2 0.291 0.0534 0.2028 0.417
69 13 1 0.268 0.0538 0.1812 0.397
71 12 2 0.224 0.0533 0.1402 0.357
74 9 1 0.199 0.0529 0.1180 0.335
75 8 1 0.174 0.0518 0.0971 0.312
80 6 1 0.145 0.0506 0.0731 0.287
81 5 1 0.116 0.0481 0.0514 0.261
100 2 1 0.058 0.0475 0.0116 0.289
111 1 1 0.000 NaN NA NA
> survminer::ggsurvplot(fit,
+ data = df,
+ conf.int = TRUE,
+ pval = TRUE,
+ surv.median.line = "hv",
+ risk.table = TRUE,
+ surv.plot.height = .6,
+ legend.position = c(75,.7),
+ risk.table.height = .3,
+ add.all = TRUE)
绘制不同访问次数的生存曲线
> df$class = ifelse(df$访问次数 > mean(df$访问次数),'较多的访问次数','较少的访问次数')
> table(df$class)
较多的访问次数 较少的访问次数
3501 7695
> fit1 = survfit(Surv(在售时长,是否卖出) ~ class,data = df)
> survminer::ggsurvplot(fit1,
+ data = df,
+ conf.int = TRUE,
+ pval = TRUE,
+ surv.median.line = "hv",
+ legend.position = c(75,.7))
>

log rank 检验
验证不同访问次数下的二手车销售生存曲线是否存在显著性差异
> surv_diff <- survdiff(Surv(在售时长, 是否卖出) ~ class, data = df)
> print(surv_diff)
Call:
survdiff(formula = Surv(在售时长, 是否卖出) ~ class, data = df)
N Observed Expected (O-E)^2/E (O-E)^2/V
class=较多的访问次数 3501 1771 2331 134.4 244
class=较少的访问次数 7695 4053 3493 89.6 244
Chisq= 244 on 1 degrees of freedom, p= <2e-16
log rank 检验所得到的结果显示p value
小于.05,我们有倾向认为两条生存曲线存在显著性差异。
半参数COX等比例风险模型预测
> fit1 = survival::coxph(Surv(在售时长,是否卖出) ~.,data = df[,-1])
> plot(survfit(fit1))
> summary(fit1)
Call:
survival::coxph(formula = Surv(在售时长, 是否卖出) ~ ., data = df[,
-1])
n= 11196, number of events= 5824
coef exp(coef) se(coef) z Pr(>|z|)
车辆行驶路程 -1.015e-05 1.000e+00 1.029e-06 -9.861 < 2e-16 ***
颜色Brown -1.023e+00 3.594e-01 1.480e-01 -6.916 4.63e-12 ***
颜色Silver 1.097e+00 2.995e+00 5.323e-02 20.606 < 2e-16 ***
颜色White 2.093e+00 8.111e+00 5.497e-02 38.076 < 2e-16 ***
照片数量 -2.791e-02 9.725e-01 5.157e-03 -5.412 6.22e-08 ***
价格 -2.508e-05 1.000e+00 5.409e-07 -46.358 < 2e-16 ***
访问次数 -3.770e-03 9.962e-01 1.529e-03 -2.465 0.0137 *
class较少的访问次数 2.027e-01 1.225e+00 3.929e-02 5.158 2.49e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
车辆行驶路程 1.0000 1.0000 1.0000 1.0000
颜色Brown 0.3594 2.7828 0.2689 0.4803
颜色Silver 2.9945 0.3339 2.6979 3.3238
颜色White 8.1111 0.1233 7.2826 9.0339
照片数量 0.9725 1.0283 0.9627 0.9824
价格 1.0000 1.0000 1.0000 1.0000
访问次数 0.9962 1.0038 0.9933 0.9992
class较少的访问次数 1.2247 0.8166 1.1339 1.3227
Concordance= 0.687 (se = 0.004 )
Likelihood ratio test= 4351 on 8 df, p=<2e-16
Wald test = 3510 on 8 df, p=<2e-16
Score (logrank) test = 3632 on 8 df, p=<2e-16
对于COX半参数回归模型的解释:
车辆行驶路程、照片数量、价格、访问次数对二手车的销售速度有着负向影响,而在分类变量中模型以黑色车为标准,棕色车相对于黑色车销售速度更慢,白色、银色相较于黑色车相对较好。
以较多的访问次数为标准,较少的访问次数确实提高了二手车的销售速度。
这里给大家遗留了一个问题R语言COX半参数回归模型如何输出每个样本的生存概率以及如何绘制每个样本的生存曲线
。