R -기초 생존 분석, PBC (primary biliary cirrhosis) 데이터

PBC 데이터를 통한 생존분석

PBC 데이터 다운받을 수 있는 사이트

http://www4.stat.ncsu.edu/~boos/var.select/pbc.html

http://www4.stat.ncsu.edu/~boos/var.select/pbc.dat.txt

이 데이터셋은 mayo clinic에서 수집된 데이터로 일차성 담즙성 간경화증(primary biliary cirrhosis) 환자들의 생존에 대한 데이터입니다. 

데이터 핸들링

library(survival)

data <- read.table("pbc.dat")

# Column rename
colnames(data) <- c('id', 'futime', 'status', 'drug', 'age', 'sex', 'ascites', 'hepato', 'spiders', 'edema', 'bili', 'chol', 'albumin', 'copper', 'alk_phos', 'sgot', 'trig', 
                   'platelet', 'protime', 'stage')

head(data)

id	futime	status	drug	age	sex	ascites	hepato	spiders	edema	bili	chol	albumin	copper	alk_phos	sgot	trig	platelet	protime	stage
1	400	2	1	21464	1	1	1	1	1.0	14.5	261	2.60	156	1718.0	137.95	172	190	12.2	4
2	4500	0	1	20617	1	0	1	1	0.0	1.1	302	4.14	54	7394.8	113.52	88	221	10.6	3
3	1012	2	1	25594	0	0	0	0	0.5	1.4	176	3.48	210	516.0	96.10	55	151	12.0	4
4	1925	2	1	19994	1	0	1	1	0.5	1.8	244	2.54	64	6121.8	60.63	92	183	10.3	4
5	1504	1	2	13918	1	0	1	1	0.0	3.4	279	3.53	143	671.0	113.15	72	136	10.9	3
6	2503	2	2	24201	1	0	1	0	0.0	0.8	248	3.98	50	944.0	93.00	63	.	11.0	3

생존분석 자료의 기본 자료 구성

1. Survival 데이터에서는 time을 나타내는 변수와 event를 나타내는 변수가 있다. 

2. time은 event 까지의 시간을 나타내며, 

3. event 변수는 우리의 관심 event (예를 들어, 사망) 혹은 censoring을 나타낸다. (이는 주로 0, 1로 나타내어진다.)

PBC 데이터에서는  status = 0=alive, 1=liver transplant, 2=dead 인데,  0,1을 중도 절단된 censoring, 2 를 관심 event로 코딩을 해봅시다.

data$status[(data$status == 1)] = 0
data$status[(data$status == 2)] = 1

drug = 1 : 페니실린

drug = 2 : 플라시보이므로 , 해석상의 이점을 얻기 위해

drug = 1 페니실린, drug = 0은 플라시보로 재코딩해줍니다.

data <- data[data$drug != '.', ]
data$drug <- as.character(data$drug)
data$drug[(data$drug == 2)] = 0
data$drug <- factor(data$drug)

1. KM estimation

우선 Survival function을 추정하는 비모수적인 방법으로 많이 쓰이는 Kaplan-Meier analysis 를 해보겠습니다.

# event, censoring 구분
Y = Surv(data$futime, data$status)

# KM estimation
fit = survfit(Y~data$drug)
summary(fit)

우선 Surv라는 함수를 통해 무엇이 time과 event인지를 알려줍니다. 그리고 drug에 따라 생존을 분석할 것이므로 ~data$drug를 입력하면 아래처럼 lifetable을 만들어줍니다.

Call: survfit(formula = Y ~ data$drug)

                data$drug=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   51    154       1    0.994 0.00647        0.981        1.000
   77    153       1    0.987 0.00912        0.969        1.000
  110    152       1    0.981 0.01114        0.959        1.000
  130    151       1    0.974 0.01282        0.949        0.999
  186    150       1    0.968 0.01428        0.940        0.996
  191    149       1    0.961 0.01559        0.931        0.992
  207    148       1    0.955 0.01679        0.922        0.988
  216    147       1    0.948 0.01788        0.914        0.984
  264    146       2    0.935 0.01986        0.897        0.975
  304    144       1    0.929 0.02075        0.889        0.970
  321    143       1    0.922 0.02160        0.881        0.965
  326    142       1    0.916 0.02240        0.873        0.961
  460    141       1    0.909 0.02317        0.865        0.956
  549    140       1    0.903 0.02389        0.857        0.951
  552    139       1    0.896 0.02459        0.849        0.946
  597    138       1    0.890 0.02525        0.841        0.941
  611    137       1    0.883 0.02589        0.834        0.935
  708    136       1    0.877 0.02650        0.826        0.930
  733    135       1    0.870 0.02709        0.819        0.925
  769    134       1    0.864 0.02765        0.811        0.920
  786    133       1    0.857 0.02820        0.804        0.914
  790    131       1    0.851 0.02873        0.796        0.909
  797    130       1    0.844 0.02925        0.789        0.903
  850    128       1    0.837 0.02975        0.781        0.898
  853    127       1    0.831 0.03024        0.774        0.892
  859    126       1    0.824 0.03071        0.766        0.887
  890    125       1    0.818 0.03116        0.759        0.881
  930    124       1    0.811 0.03160        0.751        0.875
  943    123       1    0.804 0.03203        0.744        0.870
  974    122       1    0.798 0.03244        0.737        0.864
 1080    118       1    0.791 0.03286        0.729        0.858
 1165    115       1    0.784 0.03328        0.722        0.852
 1212    114       1    0.777 0.03370        0.714        0.846
 1217    111       1    0.770 0.03411        0.706        0.840
 1356    103       1    0.763 0.03459        0.698        0.834
 1413    101       1    0.755 0.03506        0.690        0.827
 1427     98       1    0.748 0.03554        0.681        0.821
 1444     95       1    0.740 0.03603        0.672        0.814
 1487     93       1    0.732 0.03651        0.664        0.807
 1536     91       1    0.724 0.03698        0.655        0.800
 1786     79       1    0.715 0.03763        0.645        0.792
 1847     76       1    0.705 0.03829        0.634        0.784
 2090     69       1    0.695 0.03908        0.622        0.776
 2419     56       1    0.683 0.04030        0.608        0.766
 2466     53       1    0.670 0.04155        0.593        0.756
 2503     51       1    0.657 0.04276        0.578        0.746
 2769     40       1    0.640 0.04473        0.558        0.734
 2796     38       1    0.623 0.04662        0.538        0.722
 2847     35       1    0.605 0.04857        0.517        0.709
 3090     32       1    0.587 0.05060        0.495        0.695
 3170     29       1    0.566 0.05275        0.472        0.680
 3244     28       1    0.546 0.05460        0.449        0.664
 3358     26       1    0.525 0.05640        0.425        0.648
 3395     24       1    0.503 0.05814        0.401        0.631
 3428     22       1    0.480 0.05983        0.376        0.613
 3445     21       1    0.457 0.06119        0.352        0.595
 3762     15       1    0.427 0.06427        0.318        0.573
 3839     13       1    0.394 0.06719        0.282        0.551
 3853     12       1    0.361 0.06916        0.248        0.526

                data$drug=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   41    158       1    0.994 0.00631        0.981        1.000
   71    157       1    0.987 0.00889        0.970        1.000
  131    156       1    0.981 0.01086        0.960        1.000
  140    155       1    0.975 0.01250        0.950        0.999
  179    154       1    0.968 0.01393        0.941        0.996
  198    153       1    0.962 0.01521        0.933        0.992
  223    152       1    0.956 0.01637        0.924        0.988
  334    151       1    0.949 0.01744        0.916        0.984
  348    150       1    0.943 0.01844        0.908        0.980
  388    149       1    0.937 0.01937        0.900        0.975
  400    148       1    0.930 0.02025        0.892        0.971
  515    147       1    0.924 0.02108        0.884        0.966
  673    145       1    0.918 0.02187        0.876        0.962
  694    144       1    0.911 0.02263        0.868        0.957
  750    141       1    0.905 0.02337        0.860        0.952
  762    140       1    0.898 0.02408        0.852        0.947
  799    139       1    0.892 0.02476        0.845        0.942
  824    138       1    0.885 0.02541        0.837        0.937
  904    134       1    0.879 0.02607        0.829        0.931
  971    132       1    0.872 0.02671        0.821        0.926
  980    131       1    0.866 0.02732        0.814        0.921
  999    130       1    0.859 0.02791        0.806        0.915
 1000    129       1    0.852 0.02848        0.798        0.910
 1012    128       1    0.846 0.02902        0.791        0.904
 1037    127       1    0.839 0.02955        0.783        0.899
 1077    126       1    0.832 0.03005        0.775        0.893
 1083    125       1    0.826 0.03054        0.768        0.888
 1152    124       1    0.819 0.03101        0.760        0.882
 1170    122       1    0.812 0.03148        0.753        0.876
 1191    121       2    0.799 0.03236        0.738        0.865
 1235    117       1    0.792 0.03279        0.730        0.859
 1297    114       1    0.785 0.03323        0.723        0.853
 1350    111       1    0.778 0.03368        0.715        0.847
 1360    110       1    0.771 0.03411        0.707        0.841
 1434    105       1    0.764 0.03456        0.699        0.834
 1492    100       1    0.756 0.03505        0.690        0.828
 1576     97       1    0.748 0.03554        0.682        0.821
 1657     93       1    0.740 0.03606        0.673        0.814
 1682     92       1    0.732 0.03655        0.664        0.807
 1690     91       2    0.716 0.03748        0.646        0.793
 1741     87       1    0.708 0.03794        0.637        0.786
 1827     82       1    0.699 0.03845        0.628        0.779
 1925     78       1    0.690 0.03899        0.618        0.771
 2055     72       1    0.681 0.03960        0.607        0.763
 2081     71       1    0.671 0.04019        0.597        0.755
 2105     70       1    0.661 0.04074        0.586        0.746
 2224     65       1    0.651 0.04137        0.575        0.738
 2256     63       1    0.641 0.04198        0.564        0.729
 2288     61       1    0.630 0.04259        0.552        0.720
 2297     60       1    0.620 0.04315        0.541        0.710
 2386     54       1    0.608 0.04385        0.528        0.701
 2400     53       1    0.597 0.04450        0.516        0.691
 2540     47       1    0.584 0.04533        0.502        0.680
 2583     42       1    0.570 0.04634        0.486        0.669
 2598     41       1    0.556 0.04725        0.471        0.657
 2689     38       1    0.542 0.04822        0.455        0.645
 3086     28       1    0.522 0.05023        0.433        0.631
 3222     24       1    0.501 0.05264        0.407        0.615
 3282     22       1    0.478 0.05495        0.381        0.599
 3574     18       1    0.451 0.05795        0.351        0.580
 3584     17       1    0.425 0.06032        0.322        0.561
 4079      8       1    0.372 0.07247        0.254        0.545
 4191      7       1    0.319 0.07922        0.196        0.519

par(mai=c(1,1,1,1))
plot(fit,main="KM Curve", xlab="Time(Week)", ylab=expression(paste(hat(S),"(t)")))

# Log rank test (difference between groups)
log_rank = survdiff(Surv(data$futime, data$status) ~ data$drug, data)
log_rank

결과를 보시면 플라시보와 페니실린 처방군간에 생존률에 큰 차이가 없는 것을 확인할 수 있습니다.

두 군간에 차이가 있는지 검정하는 방법인 log_rank 테스트에서도 p value > 0.05로 두 군간에 차이가 유의하지 않습니다.

Call:
survdiff(formula = Surv(data$futime, data$status) ~ data$drug, 
    data = data)

              N Observed Expected (O-E)^2/E (O-E)^2/V
data$drug=0 154       60     61.8    0.0513     0.102
data$drug=1 158       65     63.2    0.0502     0.102
 

Chisq= 0.1 on 1 degrees of freedom, p= 0.7

2. Cox regression

생존분석을 다변수로 할 때, Cox regression을 주로 활용합니다. 이는 Cox proportional hazard model로 불리기도 합니다.

Cox regression은 다음과 같은 이유로 많이 쓰입니다.

1. 다양한 변수들의 계수와 hazard ratio를 추정하는데 좋다. -> 이를 통해 변수별 생존에 대한 중요도를 비교하기 쉽다.

2. Cox regression은 robust 하다. (참에 매우 근접한 값을 추정해준다.)

3. 계산된 hazard가 양수이기 때문에 타당한 값이다.

4. baseline hazard를 몰라도 회귀계수, hazard ratio를 추정할 수 있다.

Cox regression의 이론적 배경에 대해서는 다음에 기회가 있을 때 다루어보도록 하겠습니다. 

R에서 Cox regression을 수행하는 것은 매우 간단한데, 앞서 정의한 event, time 변수인 Y를 종속변수로 놓고 분석대상인 변수들을 ~ +로 연결하면 됩니다.  

# Cox regression
cox = coxph(Y~data$drug + data$age + data$sex + data$stage + data$edema + data$bili + data$albumin)
summary(cox) 

Call:
coxph(formula = Y ~ data$drug + data$age + data$sex + data$stage + 
    data$edema + data$bili + data$albumin)

  n= 312, number of events= 125 

                   coef  exp(coef)   se(coef)      z Pr(>|z|)    
data$drug1    2.569e-02  1.026e+00  1.855e-01  0.138 0.889882    
data$age      6.736e-05  1.000e+00  2.651e-05  2.541 0.011057 *  
data$sex     -5.940e-01  5.521e-01  2.550e-01 -2.330 0.019832 *  
data$stage1  -2.248e+00  1.056e-01  1.023e+00 -2.197 0.027992 *  
data$stage2  -8.573e-01  4.243e-01  2.950e-01 -2.906 0.003659 ** 
data$stage3  -4.884e-01  6.136e-01  2.186e-01 -2.234 0.025487 *  
data$stage4          NA         NA  0.000e+00     NA       NA    
data$edema    9.740e-01  2.649e+00  3.148e-01  3.095 0.001971 ** 
data$bili     1.263e-01  1.135e+00  1.529e-02  8.263  < 2e-16 ***
data$albumin -9.288e-01  3.950e-01  2.556e-01 -3.634 0.000279 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

             exp(coef) exp(-coef) lower .95 upper .95
data$drug1      1.0260     0.9746   0.71322    1.4760
data$age        1.0001     0.9999   1.00002    1.0001
data$sex        0.5521     1.8112   0.33497    0.9101
data$stage1     0.1056     9.4656   0.01423    0.7844
data$stage2     0.4243     2.3568   0.23801    0.7564
data$stage3     0.6136     1.6296   0.39979    0.9419
data$stage4         NA         NA        NA        NA
data$edema      2.6486     0.3776   1.42921    4.9083
data$bili       1.1347     0.8813   1.10117    1.1692
data$albumin    0.3950     2.5315   0.23936    0.6519

Concordance= 0.839  (se = 0.029 )
Rsquare= 0.427   (max possible= 0.983 )
Likelihood ratio test= 174  on 9 df,   p=<2e-16
Wald test            = 186.2  on 9 df,   p=<2e-16
Score (logrank) test = 287  on 9 df,   p=<2e-16

그러면 이렇게 회귀계수와 회귀 계수의 exponential을 돌려주게 됩니다. Cox regression이 좋은 점은 회귀계수에 exponential을 취했을 때, 해당 변수가 1단위 증가했을 때의 hazard ratio가 됩니다. 예를 들어, age의 경우, 다른 변수들을 보정하고, 1살 증가했을 때, hazard가 1.0001배 증가한다. 라고 볼 수 있습니다. sex의 경우, 다른 모든 변수를 보정했을 때, 남자에 비해 여자가 hazard가 0.55 배라고 볼 수 있습니다 (이 데이터에서 0=male, 1=female로 코딩이 되어있기 때문에). p-value를 보면 어떠한 변수가 생존에 영향을 크게 끼치는 지를 비교해볼 수 있습니다.