【入門】個人生存數據分析

Primary Survival analysis

一、數據說明

數據是5.3號爬取的選手「17_shou」在亞服的30場比賽數據。 數據裡面包含了: -index（序號）； -time（生存時間）； -status（生存狀態，死亡=1，生存=0）； -heals（回復用品使用的次數）（離散型）； -boosts（能量用品使用的次數）（離散型）； 利用這些數據可以做出最基本的生存分析。

二、數據打開

library(readxl)

pubgsr<-read_xlsx("C:/Users/Fan/Desktop/pubgsr.xlsx")

str(pubgsr)#查看各變數類型

## Classes tbl_df, tbl and data.frame: 30 obs. of 8 variables:

## $ 序號 : num 1 2 3 4 5 6 7 8 9 10 ...

## $ 回復用品的使用數量: num 12 10 3 0 7 1 7 0 0 3 ...

## $ 能量用品的使用數量: num 7 8 7 0 1 0 13 1 0 4 ...

## $ time : num 1850 1675 1314 163 1492 ...

## $ 是否存活 : chr "byplayer" "byplayer" "byplayer" "byplayer" ...

## $ status : chr "1" "1" "1" "1" ...

## $ heals : num 1 1 1 0 1 1 1 0 0 1 ...

## $ boosts : num 1 1 1 0 1 0 1 1 0 1 ...

pubgsr$status<-as.numeric(pubgsr$status)

pubgsr$boosts<-as.numeric(pubgsr$boosts)

pubgsr$heals<-as.numeric(pubgsr$heals)#將需要的變數轉化為數值型

str(pubgsr)#檢查一遍

## Classes tbl_df, tbl and data.frame: 30 obs. of 8 variables:

## $ 序號 : num 1 2 3 4 5 6 7 8 9 10 ...

## $ 回復用品的使用數量: num 12 10 3 0 7 1 7 0 0 3 ...

## $ 能量用品的使用數量: num 7 8 7 0 1 0 13 1 0 4 ...

## $ time : num 1850 1675 1314 163 1492 ...

## $ 是否存活 : chr "byplayer" "byplayer" "byplayer" "byplayer" ...

## $ status : num 1 1 1 1 1 1 0 1 1 1 ...

## $ heals : num 1 1 1 0 1 1 1 0 0 1 ...

## $ boosts : num 1 1 1 0 1 0 1 1 0 1 ...

Surv：用於創建生存數據對象 survfit：創建KM生存曲線或是Cox調整生存曲線 survdiff：用於不同組的統計檢驗

三、構築生存數據對象

library(survival)

library(survminer)

## Loading required package: ggplot2

## Loading required package: ggpubr

## Loading required package: magrittr

attach(pubgsr)

surv<-Surv(pubgsr$time,pubgsr$status)#格式是Surv(時間，事件)

summary(surv)

## time status

## Min. : 87.0 Min. :0.0

## 1st Qu.: 234.2 1st Qu.:1.0

## Median : 846.5 Median :1.0

## Mean : 868.5 Mean :0.9

## 3rd Qu.:1451.5 3rd Qu.:1.0

## Max. :1910.0 Max. :1.0

四、非參數法——Kaplan-Meier乘積極限法

pubgfit1<-survfit(surv~1)#表示不分組

summary(pubgfit1)

## Call: survfit(formula = surv ~ 1)

##

## time n.risk n.event survival std.err lower 95% CI upper 95% CI

## 87 30 1 0.967 0.0328 0.9045 1.000

## 88 29 1 0.933 0.0455 0.8482 1.000

## 93 28 1 0.900 0.0548 0.7988 1.000

## 94 27 1 0.867 0.0621 0.7532 0.997

## 163 26 1 0.833 0.0680 0.7101 0.978

## 164 25 1 0.800 0.0730 0.6689 0.957

## 206 24 1 0.767 0.0772 0.6293 0.934

## 232 23 1 0.733 0.0807 0.5910 0.910

## 241 22 1 0.700 0.0837 0.5538 0.885

## 266 21 1 0.667 0.0861 0.5176 0.859

## 314 20 1 0.633 0.0880 0.4824 0.832

## 343 19 1 0.600 0.0894 0.4480 0.804

## 452 18 1 0.567 0.0905 0.4144 0.775

## 496 17 1 0.533 0.0911 0.3816 0.745

## 725 16 1 0.500 0.0913 0.3496 0.715

## 968 15 1 0.467 0.0911 0.3183 0.684

## 1083 14 1 0.433 0.0905 0.2878 0.652

## 1150 13 1 0.400 0.0894 0.2581 0.620

## 1174 12 1 0.367 0.0880 0.2291 0.587

## 1286 11 1 0.333 0.0861 0.2010 0.553

## 1314 10 1 0.300 0.0837 0.1737 0.518

## 1330 9 1 0.267 0.0807 0.1473 0.483

## 1492 8 1 0.233 0.0772 0.1220 0.446

## 1543 7 1 0.200 0.0730 0.0978 0.409

## 1580 6 1 0.167 0.0680 0.0749 0.371

## 1675 5 1 0.133 0.0621 0.0535 0.332

## 1850 4 1 0.100 0.0548 0.0342 0.293

plot(pubgfit1,main="Survivor",xlim=c(0,2165),col = "red",xlab = "time",ylab="rate")

函數圖像的縱軸是生存的概率，橫軸是生存的時間，兩條虛線代表著置信區間根據@特玩網提供的毒圈數據(2017.11.16),比賽最長的持續時間為36分05秒，即2165秒。

我們從函數圖像可以看出，17_shou在500秒之前是生存率下降的最快的階段。換算成毒圈，就是在第一個圈還沒縮完之前的生存率下降的比較快，這個從17_shou的直播里也能夠驗證這一點，主播在天梯當中一般都是跳機場、P城；沙漠圖則是Pecado和皇宮這些鋼槍的地方，極其影響生存率。

從生存數據統計表中可以看到，17_shou這30長比賽的中位生存時間是725秒，剛好是第一個毒圈縮完。可以看到17_shou鋼槍是很兇了。

五、生存分組並進行生存比較

上面我們僅僅查看了選手在30場比賽的生存率分析，但是還沒有探究有什麼因素影響，現在我們構建以回復用品和能量用品為分組的依據並用以這30場比賽的分組。以「有無使用該物品」作為分類的依據，對於heals變數和boosts變數：1=有，0=無。下面我們構築分組

pubgfit2<-survfit(Surv(time,status)~heals+boosts,data = pubgsr)

summary(pubgfit2)

## Call: survfit(formula = Surv(time, status) ~ heals + boosts, data = pubgsr)

##

## heals=0, boosts=0

## time n.risk n.event survival std.err lower 95% CI upper 95% CI

## 87 7 1 0.857 0.132 0.6334 1.000

## 88 6 1 0.714 0.171 0.4471 1.000

## 93 5 1 0.571 0.187 0.3008 1.000

## 94 4 1 0.429 0.187 0.1822 1.000

## 163 3 1 0.286 0.171 0.0886 0.922

## 164 2 1 0.143 0.132 0.0233 0.877

## 241 1 1 0.000 NaN NA NA

##

## heals=0, boosts=1

## time n.risk n.event survival std.err lower 95% CI upper 95% CI

## 266 2 1 0.5 0.354 0.125 1

## 452 1 1 0.0 NaN NA NA

##

## heals=1, boosts=0

## time n.risk n.event survival std.err lower 95% CI upper 95% CI

## 206 2 1 0.5 0.354 0.125 1

## 314 1 1 0.0 NaN NA NA

##

## heals=1, boosts=1

## time n.risk n.event survival std.err lower 95% CI upper 95% CI

## 232 19 1 0.947 0.0512 0.8521 1.000

## 343 18 1 0.895 0.0704 0.7669 1.000

## 496 17 1 0.842 0.0837 0.6931 1.000

## 725 16 1 0.789 0.0935 0.6259 0.996

## 968 15 1 0.737 0.1010 0.5632 0.964

## 1083 14 1 0.684 0.1066 0.5041 0.929

## 1150 13 1 0.632 0.1107 0.4480 0.890

## 1174 12 1 0.579 0.1133 0.3946 0.850

## 1286 11 1 0.526 0.1145 0.3435 0.806

## 1314 10 1 0.474 0.1145 0.2949 0.761

## 1330 9 1 0.421 0.1133 0.2485 0.713

## 1492 8 1 0.368 0.1107 0.2045 0.664

## 1543 7 1 0.316 0.1066 0.1629 0.612

## 1580 6 1 0.263 0.1010 0.1240 0.558

## 1675 5 1 0.211 0.0935 0.0881 0.503

## 1850 4 1 0.158 0.0837 0.0559 0.446

text<-c("heals=0,boosts=0","heals=0,boosts=1","heals=1,boosts=0","heals=1,boosts=1")

plot(pubgfit2,main="Survivor2",xlim=c(0,2165),xlab = "time",ylab="rate",col = c(2,3,4,5))#或直接採用ggsurvplot()

legend(1400,1,lty=1:2,legend=text,col=c(2,3,4,5))

#diff()差分函數

diff1<-survdiff(surv~heals,data = pubgsr,rho=1)#rho=1是Wilcon法檢驗，更容易發現早期生存率差異,生存時間呈對數正太分布時效率較高

diff2<-survdiff(surv~boosts,data = pubgsr)#rho=0,是log-rank時序檢驗，呈weibull分布或屬於比例風險模型時效率較高

#似然比檢驗則是用於生存時間呈指數分布

diff3<-survdiff(surv~heals+boosts,data = pubgsr,rho = 1)

diff1

## Call:

## survdiff(formula = surv ~ heals, data = pubgsr, rho = 1)

##

## N Observed Expected (O-E)^2/E (O-E)^2/V

## heals=0 9 7.53 1.77 18.82 27.4

## heals=1 21 7.77 13.53 2.46 27.4

##

## Chisq= 27.4 on 1 degrees of freedom, p= 1.63e-07

diff2

## Call:

## survdiff(formula = surv ~ boosts, data = pubgsr)

##

## N Observed Expected (O-E)^2/E (O-E)^2/V

## boosts=0 9 9 1.8 28.78 36.4

## boosts=1 21 18 25.2 2.06 36.4

##

## Chisq= 36.4 on 1 degrees of freedom, p= 1.61e-09

diff3

## Call:

## survdiff(formula = surv ~ heals + boosts, data = pubgsr, rho = 1)

##

## N Observed Expected (O-E)^2/E (O-E)^2/V

## heals=0, boosts=0 7 6.23 1.000 27.388 35.879

## heals=0, boosts=1 2 1.30 0.767 0.371 0.489

## heals=1, boosts=0 2 1.47 0.600 1.252 1.565

## heals=1, boosts=1 19 6.30 12.933 3.402 29.619

##

## Chisq= 41.7 on 3 degrees of freedom, p= 4.75e-09

從函數圖像可以看出，17_shou在比賽當中使用了回復和能量兩種物品時生存率要顯著的增加，而使用單一種物品或者不使用的時候，生存率明顯低下。比較有趣的是，17_shou僅使用能量物品的生存時間要比僅使用回復藥品的生存時間要多出一點點。 查看差分函數的結果後可以發現，無論怎麼樣分組，分組的P值都是小於0.01,組間差異性顯著。

六、參數模型構建

test1<-survreg(surv~heals+boosts,data = pubgsr,dist="exponential")#dist是選擇分布的模式

summary(test1)

##

## Call:

## survreg(formula = surv ~ heals + boosts, data = pubgsr, dist = "exponential")

## Value Std. Error z p

## (Intercept) 4.81 0.348 13.82 1.89e-43

## heals 1.07 0.626 1.71 8.78e-02

## boosts 1.39 0.626 2.22 2.63e-02

##

## Scale fixed at 1

##

## Exponential distribution

## Loglik(model)= -201.3 Loglik(intercept only)= -212.5

## Chisq= 22.56 on 2 degrees of freedom, p= 1.3e-05

## Number of Newton-Raphson Iterations: 3

## n= 30

生存分析完整實例：http://blog.sciencenet.cn/blog-252888-719677.html

七、半參數法——Cox比例風險回歸模型

這個模型不需要對生存時間作出假定，卻可以通過模型來分析生存時間的分布規律和危險因素對生存時間的影響。但是此模型的前提是危險因素在每個生存時間分布中呈現同樣的影響。從圖像上看即各個樣本的生存曲線是不能存在交叉的，若有交叉代表著需要分層分類。此模型可以處理連續型變數

pubgfit3<-coxph(Surv(pubgsr$time,pubgsr$status)~heals+boosts,data = pubgsr)

plot(survfit(pubgfit3))

summary(pubgfit3)

## Call:

## coxph(formula = Surv(pubgsr$time, pubgsr$status) ~ heals + boosts,

## data = pubgsr)

##

## n= 30, number of events= 27

##

## coef exp(coef) se(coef) z Pr(>|z|)

## heals -2.09516 0.12305 0.76257 -2.747 0.00601 **

## boosts -2.96782 0.05142 0.92756 -3.200 0.00138 **

## ---

## Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

##

## exp(coef) exp(-coef) lower .95 upper .95

## heals 0.12305 8.127 0.027604 0.5485

## boosts 0.05142 19.449 0.008347 0.3167

##

## Concordance= 0.763 (se = 0.043 )

## Rsquare= 0.671 (max possible= 0.992 )

## Likelihood ratio test= 33.33 on 2 df, p=5.775e-08

## Wald test = 19.27 on 2 df, p=6.547e-05

## Score (logrank) test = 45.3 on 2 df, p=1.454e-10

test.ph<-cox.zph(pubgfit3,transform=rank)#PH假設的統計檢驗

test.ph

## rho chisq p

## heals -0.03200 0.03774 0.846

## boosts 0.00707 0.00139 0.970

## GLOBAL NA 0.03824 0.981

ggcoxzph(test.ph)#虛線表示擬合曲線上下2個單位的標準差。如果曲線偏離2個單位的標準差則表示不滿足比例風險假定

-函數結果分析 heals和boosts兩個變數的coef<0,均為保護性因素。p值均<0.01,變數顯著。 R方=0.672,說明函數解釋了67.2%的變數 似然比檢驗p值、Wald和Score檢驗的P值均<0.01，函數擬合不錯

-圖像分析 從圖像上看，能夠比較吻合亞服天梯中的存活情況，第一個圈的時候基本上就已經掉了一半的人，生存率驟降，但是到後面第二三個圈的時候存活率就慢慢穩定，掉人的速度就沒那麼快了。

-意義分析（參考https://wenku.baidu.com/view/8c289a7e910ef12d2bf9e72e.html?pn=50） 回歸方程:h(t)=h0(t)exp(-2.09516X1 + 0.05142X2) 不使用回復用品比使用回復用品對於死亡的危險度為0.12305，不使用能量藥品比使用能量藥品對於死亡的危險度為0.05142。可以看出不用回復藥品要比不用能量藥品更加的容易死亡。

推薦閱讀：