数据科学案例之生存分析与二手车定价

数据科学与机器学习案例之客户的信用风险与预测

数据科学与机器学习案例之信用卡欺诈识别（严重类失衡数据建模)

数据科学与机器学习案例之汽车目标客户销售策略研究

数据科学与机器学习案例之WiFi定位系统的位置预测

数据科学与机器学习案例之Stacking集成方法对鸢尾花进行分类

生存分析

生存分析的基本概念

生存分析是一种将事件的结果与出现此类结果所经历的时间结合在一起的统计分析方法。生存分析包含了：事件、生存时间、删失数据、生存函数、风险函数、半衰期等概念。
事件：起始事件与终点事件。比如我们在研究二手车分析定价中起始事件是二手汽车平台发布车辆出售信息，终点事件就是用户下单购买二手车。
生存时间：生存时间不关注事件发生的瞬时时间，只关注事件的起始事件开始到终点事件发生之间的时间间隔。生存时间服从偏态分布：指数分布、 威布尔分布，对数logistic分布。
删失数据：是指在观察时间窗内无法获取事件发生时间，其中包括右删失、左删失与区间删失。右删失是指知道观察的起始时间，但无法获取终点事件发生的具体时间，主要由以下的原因导致：研究对象在观察时间窗内未发生有效事件、研究对象在观察时间窗内由于某些原因被丢失、研究对象在事件发生之前由于非事件原因脱离有效观测。

生存函数：对于观察事件窗内的任意时刻，生存函数是到该时刻仍未发生该事件的概率。生存函数是每个时刻生存概率的乘积，$S(t)=p(T>t)=1-F(t)$，其中$F(t)$代表生存时间的累积分布函数，表示事件发生时间未超过时刻$t$的概率。生存函数具有以下三个性质：$S(t)\in [0,1]$,且$S(t)$单调递减；在起始时刻$t=0$时，所有对象均处于存活状态，此时$S(t)=1$；当时间$t$趋于无穷大时，生存概率趋近于0，$S(t)=0$.
风险函数与累积风险函数：风险函数又称条件死亡率，是指在时刻$t$之前未发生任何事件而恰好在时刻$t$发生事件的概率。风险函数的公式为：$h(t)={\lim }_{h\to 0}\frac{p(T<t+h|T\ge t)}{h}$，当生存时间$t$为连续型随机变量时，风险函数表示为：$h(t)=\frac{f(t)}{S(t)}=\frac{-dln[S(t)]}{dt}$，累积风险函数表示为：$H(t)={\int }_0^{t}h(u)du=-ln[S(t)]$，生存函数与风险函数的关系表达为：$S(t)=e^{-H(t)}=e^{-{\int }_0^{t}h(u)du}$.
中位生存时间：指一半的个体未发生终点事件的时间。生存分析中由于删失数据的存在，无法直接以平均时间反应时间发生的时间水平，因此生存分析中常以存活概率$50\%$时对应消耗的时间描述。

生存分析的统计分析方法

对生存函数的刻画：KM曲线法，假设$t_1、t_2、t_3、，，，t_k$共有$k$个观测时刻及$N$个样本，$a_j$代表在时刻$t_j$发生事件的人数，$b_j$是在$(t_j,t_{j+1})$时间段内删失的人数，则在$t_j$时刻处于风险中的人数表示为：$c_j=(a_j+b_j)+(a_{j+1}+b_{j+1})+...+(a_k+b_k)$,
则$t_j$时刻的风险率为$\frac{d_j}{n_j}$，KM曲线对应生存函数的估计为：$S(t)=\prod _{j:t_j\le t}\frac{n_j-d_j}{n_j}$.
通过直接观察我们无法确定生存曲线之间是否具有显著性差异，引入了log-rank与Breslow检验方法。
Nelson-Aalen累计风险曲线：$H(t)=\prod _{j:t_j\le t}\frac{d_j}{n_j}$.
生存函数的参数估计：参数模型是假定生存时间是符合某种分布，通常情形下生存时间的分布主要有指数分布、威布尔分布、对数正态分布和对数logistic分布。下面只是对指数分布做简要介绍：指数分布的生存函数为：$S(t)=e^{-\lambda t}$，风险函数表示为：$h(t)=\lambda$.

半参COX等比例风险回归模型

解决问题：KM曲线只能在单一的分类变量间进行比较，也就是说KM曲线只描述了该单变量与生存时间之间的关系而忽略了其他变量的影响，由此产生了半参COX等比例风险回归模型.
半参COX等比例风险回归模型公式为：$h(t;x)=h_0(t)e^{\beta x}$，模型中的$h_0(t)$为非参数部分；$({\beta }_1,{\beta }_2,{\beta }_3,...{\beta }_p)$为模型的参数。$h_0(t)$是所有自变量取值为$0$的风险率。每个自变量对生存率的影响不随时间的改变而改变。
对于模型的解释：当自变量$x_p$为连续性变量时，回归系数${\beta }_p$是除了$x_p$之外其他变量均不变的情况下，$x_p$每增加一个单位，风险率变化$e^{{\beta }_p}$倍；当自变量为分类型变量时，假设$x_p$的类别有两类：$0$与$1$，回归系数${\beta }_j$表示除了$x_p$之外其他变量均不变的情况下，类别$1$相对于类别$0$的风险率是$e^{{\beta }_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) 


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

绘制不同访问次数的生存曲线

> 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)) 
> 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

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 

1
2
3
4
5
6
7
8
9
10
11

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

对于COX半参数回归模型的解释：
车辆行驶路程、照片数量、价格、访问次数对二手车的销售速度有着负向影响，而在分类变量中模型以黑色车为标准，棕色车相对于黑色车销售速度更慢，白色、银色相较于黑色车相对较好。
以较多的访问次数为标准，较少的访问次数确实提高了二手车的销售速度。

这里给大家遗留了一个问题R语言COX半参数回归模型如何输出每个样本的生存概率以及如何绘制每个样本的生存曲线。