데이터마이너를 꿈꾸며

6.5 준모수적 방법

Cox의 비례 모형

모수적 모형은 가정이 타당할 때는 상당히 효율적이지만 어떤 모형이 적당한가에 대한 지식이 없다면 함부로 사용하기가 곤란하다. 이에 반해 Cox(1972)의 비례 위험 모형(proportional hazards model)은 준모수적(semi-parametric)방법으로서 생존 시간의 분포에 대한 가정을 필요로 하지 않는다. 또한 시간에 따라 바뀌는 공변량(time-dependent variable)의 경우에도 분석할 수 있다는 장점이 있어, 생존자료의 분석에 매우 자주 사용된다.

[예] Prenctice(1973)에 소개된 것으로 40명의 폐암 환자의 생존시간을 조사한 것이다. 40명의 환자 중 21명은 기존 치료방법인 처리1에, 나머지 19명은 새로운 치료방법은 처리2에 할당되었으며, 생존시간에 영향을 미칠 것으로 생각되는 공변량은 다음과 같다.

X1: 진단시의 환자상태(Perfermance Status : 0~100 점)

X2: 환자의 나이(단위 : 년)

X3: 진단 후 연구 참여시까지의 시간(단위 : 월)

Trt: 치료 방법 (1 – 기존 치료, 2 – 새로운 치료)

TYPE: 종양의 유형 (squamous, small, adeno, large)

1. 데이터 입력

> setwd('c:/Rwork/')

> lung<-read.table('lung.txt',header=T)

> colnames(lung)<-c('time','status','x1','x2','x3','trt','type')

> head(lung)

  time status x1 x2 x3 trt type

1  411      1 70 64  5   1    1

2  126      1 60 63  9   1    1

3  118      1 70 65 11   1    1

4   82      1 40 69 10   1    1

5    8      1 40 63 58   1    1

6   25      0 70 48  9   1    1

2. Cox비례위험모형에 근거한 생존시간 분석

> library(survival)

> coxfit1 = coxph(Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung)

> summary(coxfit1)

Call:

coxph(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +

    factor(type), data = lung)

  n= 40, number of events= 37

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

x1            -0.060281  0.941500  0.013777 -4.375 1.21e-05 ***

x2            -0.015086  0.985027  0.022340 -0.675   0.4995   

x3             0.001201  1.001201  0.011886  0.101   0.9195   

factor(trt)2  -0.448171  0.638795  0.431302 -1.039   0.2988   

factor(type)2  0.279682  1.322709  0.547259  0.511   0.6093   

factor(type)3  1.418190  4.129638  0.625283  2.268   0.0233 * 

factor(type)4  0.361145  1.434971  0.479210  0.754   0.4511   

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

x1               0.9415     1.0621    0.9164    0.9673

x2               0.9850     1.0152    0.9428    1.0291

x3               1.0012     0.9988    0.9781    1.0248

factor(trt)2     0.6388     1.5654    0.2743    1.4876

factor(type)2    1.3227     0.7560    0.4525    3.8663

factor(type)3    4.1296     0.2422    1.2125   14.0655

factor(type)4    1.4350     0.6969    0.5610    3.6707

Concordance= 0.764  (se = 0.058 )

Rsquare= 0.524   (max possible= 0.994 )

Likelihood ratio test= 29.66  on 7 df,   p=0.0001097

Wald test            = 26.29  on 7 df,   p=0.0004479

Score (logrank) test = 30.67  on 7 df,   p=7.138e-05

- coxph: cox의 비례위험 모형을 ㅈ거합함.

-> 결과 해석:

- 진단시의 환자상태(x1)에 대한 회귀계수가 -0.060281로 유의하다(p-value < 0.001). 위험률(exp(coef))를 보면 진단시의 환자상태에 대한 점수가 1점 더 높은 환자는 다른 조건이 동일한 환자와 비교할 때 0.9415(exp(-0.060281))배 위험하다는 것을 알 수 있다.

exp(-coef) = 1.0621은 exp(coef)의 역수인 1/0.9415이다.

- 나이(x2), 진단 후 연구 참여시까지의 시간(x3), 처리(trt)는 모두 p-value가 0.05보다 커서 유의한 영향을 미치지 않는다.

- TYPE: 종양의 유형 (squamous, small, adeno, large)을 보면 factor(type)3인 adeno의 경우 squamous와 차이를 보인다.(p-value: 0.0233)

회귀계수가 0이라는 가설에 대한 검정통계량. 차례로 우도비 검정(-2 LOG L), Wald 검정통계량, Score 검정이 있다. 세 통계랑 모두 가설을 기각하게 되므로 (p-value < 0.05), 적어도 하나의 공변량은 의미가 있다는 것(= 적어도 모든 변수가 0이라는 귀무가설을 기각)을 알 수 있다.

Likelihood ratio test= 29.66  on 7 df,   p=0.0001097

Wald test            = 26.29  on 7 df,   p=0.0004479

Score (logrank) test = 30.67  on 7 df,   p=7.138e-05

3. 생존함수추정치

> fit4 = survfit(coxfit1)

> summary(fit4)

Call: survfit(formula = coxfit1)

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

    1     40       1 9.90e-01 1.10e-02     9.68e-01        1.000

    2     39       1 9.79e-01 1.67e-02     9.47e-01        1.000

    8     38       2 9.54e-01 2.71e-02     9.03e-01        1.000

   10     36       1 9.37e-01 3.32e-02     8.74e-01        1.000

   11     35       1 9.20e-01 3.87e-02     8.47e-01        0.999

   12     34       2 8.83e-01 4.84e-02     7.93e-01        0.984

   15     32       1 8.64e-01 5.33e-02     7.65e-01        0.975

   16     31       1 8.44e-01 5.79e-02     7.38e-01        0.965

   18     30       1 8.22e-01 6.27e-02     7.07e-01        0.954

   19     29       1 7.96e-01 6.77e-02     6.74e-01        0.941

   20     28       1 7.67e-01 7.27e-02     6.37e-01        0.924

   21     27       1 7.34e-01 7.73e-02     5.97e-01        0.903

   43     25       1 6.97e-01 8.22e-02     5.54e-01        0.879

   44     24       1 6.61e-01 8.62e-02     5.12e-01        0.854

   51     23       1 6.26e-01 8.94e-02     4.73e-01        0.828

   54     22       1 5.84e-01 9.23e-02     4.28e-01        0.796

   56     21       1 5.44e-01 9.39e-02     3.87e-01        0.763

   82     20       1 5.06e-01 9.46e-02     3.50e-01        0.730

   84     19       1 4.64e-01 9.54e-02     3.11e-01        0.694

   90     18       1 4.25e-01 9.51e-02     2.74e-01        0.659

  100     17       1 3.81e-01 9.46e-02     2.34e-01        0.620

  118     15       1 3.35e-01 9.36e-02     1.94e-01        0.580

  126     14       1 2.93e-01 9.12e-02     1.59e-01        0.540

  153     13       1 2.53e-01 8.79e-02     1.28e-01        0.500

  164     12       1 2.14e-01 8.37e-02     9.97e-02        0.461

  177     11       1 1.80e-01 7.81e-02     7.66e-02        0.421

  200     10       1 1.40e-01 7.08e-02     5.20e-02        0.377

  201      9       1 1.07e-01 6.19e-02     3.41e-02        0.333

  231      8       1 8.00e-02 5.25e-02     2.21e-02        0.289

  250      6       1 5.18e-02 4.19e-02     1.06e-02        0.253

  287      5       1 2.88e-02 2.98e-02     3.81e-03        0.218

  340      4       1 9.18e-03 1.50e-02     3.74e-04        0.225

  411      3       1 2.19e-03 5.10e-03     2.30e-05        0.210

  991      2       1 1.28e-04 5.20e-04     4.54e-08        0.362

  999      1       1 3.36e-10 5.39e-09     7.43e-24        1.000

 > names(fit4)

 [1] "n"         "time"      "n.risk"    "n.event"   "n.censor"  "surv"     

 [7] "type"      "cumhaz"    "std.err"   "upper"     "lower"     "conf.type"

[13] "conf.int"  "call"     

> fit4$surv

 [1] 9.895862e-01 9.787704e-01 9.543895e-01 9.372414e-01 9.195316e-01

 [6] 8.834541e-01 8.637107e-01 8.439345e-01 8.215162e-01 7.962148e-01

[11] 7.670584e-01 7.344065e-01 7.344065e-01 6.974181e-01 6.609451e-01

[16] 6.256915e-01 5.835925e-01 5.436422e-01 5.055202e-01 4.643798e-01

[21] 4.249339e-01 3.809762e-01 3.809762e-01 3.352450e-01 2.932721e-01

[26] 2.533265e-01 2.142532e-01 1.795911e-01 1.401508e-01 1.065182e-01

[31] 8.001270e-02 5.176743e-02 2.881620e-02 9.176925e-03 2.194374e-03

[36] 1.281990e-04 3.359491e-10

4. 생존함수그래프

> plot(survfit(coxfit1),xlab='time',ylab='Survival function',xlim=c(0,998.9))

> legend(500,1.0,c('누적한계추정치','95%신뢰구간'),lty=c(1,2))

5. 누적함수그래프

누적위험함수의 추정치 그래프를 통해 위험함수의 형태를 짐작할 수 있다. 예를 들어 단조 증가하는 직선형태는 위험함수가 시간에 대해 일정하다는 것을 의미하며, 위쪽으로 휘는 모양이면 시간이 지남에 따라 일정하다는 것을 의미하며, 위쪽으로 휘는 모양이면 시간이 지남에 따라 위험함수가 증가하고, 아래 방향으로 휘면 감소한다는 것을 의미한다.

> H.hat = -log(fit4$surv)

> H.hat = c(H.hat,tail(H.hat,1)) 

> plot(c(fit4$time,1100),H.hat,xlab='time',ylab='comulative hazard function',type='s')

 - tail: 벡터, 매트릭스 데이터에서 마지막 n개의 행들을 선택함. 예를 들어 tail(H.hat,1)라고 입력하면 H.hat 벡터의 마지막 1개의 성분을 선택. 해당 처리를 하는 이유는 위험함수 곡선은 시간에 대한 상승곡선이므로, 마지막 함수값을 추가함을써 곡선의 모양을 자연스럽게 하기 위한 처리임.

6. 비례성 검토를 위한 로그-로그 그림

비례성의 가정이 타당한 것인지 검토하는 방법에 대해 알아보자. 폐암자료에서 두 처리그룹별로 시간에 따른 log(-log S^ (t)) 의 그래프를 그렸을 때 평행하게 되는 가를 볼 수 있다.

> coxfit2 = coxph(Surv(time,status)~x1+x2+x3+strata(trt)+factor(type),data=lung) 

> plot(survfit(coxfit2),fun='cloglog',lty=1:2,col=c('red','blue'))#fun='cloglog', 두 그래프가 대체적으로 평행하므로 비례성이 타당하다고 할 수 있다.

> legend('topleft',c('처리1','처리2'),lty=1:2,col=c('red','blue'))

- strata(trt): 처리를 층으로 입력해주면, 처리를 층으로 입력하여, 처리1과 처리2로 나누어 Cox 비례 위험 모형을 추정함.

- fun='cloglog'를 입력하면 로그-로그 그림을 그릴 수 있다.

- 비례성검토를 위한 로그-로그 그림(log-log plot)은 그래프가 평행하므로 비례성의 가정이 타당하다 할 수 있다.

7. 로그 – 랭크 테스트 log – rank test

1) 종양 유형 별 생존함수의 차이

> survdiff(Surv(time,status)~factor(type),data=lung)

Call:

survdiff(formula = Surv(time, status) ~ factor(type), data = lung)

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

factor(type)=1 14       12    16.72     1.330     2.825

factor(type)=2 11       10     6.93     1.356     1.771

factor(type)=3  5        5     2.18     3.651     4.099

factor(type)=4 10       10    11.17     0.123     0.187

 Chisq= 7.4  on 3 degrees of freedom, p= 0.0614

X²(df=3)=7.4이고, p-value가 0.0614로 0.05보다는 약간 크지만 상당히 의미 있음을 알 수 있다.

2) 치료법 생존함수의 차이

> survdiff(Surv(time,status)~x1,data=lung)

Call:

survdiff(formula = Surv(time, status) ~ x1, data = lung)

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

x1=20 2        2    0.128   27.3119   28.6351

x1=30 4        4    1.565    3.7893    4.2568

x1=40 7        7    1.932   13.2948   15.1966

x1=50 4        3    3.298    0.0270    0.0305

x1=60 7        7    6.642    0.0193    0.0247

x1=70 9        7   12.420    2.3650    3.7914

x1=80 6        6    6.982    0.1382    0.1761

x1=90 1        1    4.033    2.2811    3.7772

 Chisq= 58.8  on 7 degrees of freedom, p= 2.65e-10

X²(df=7)=58.8이고, p-value가 0.05보다 매우 작기 때문에 의미 있음을 알 수 있다.

8. anova를 이용하여 모형 비교

> coxfit1 = coxph(Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung)

> coxfit2 = coxph(Surv(time,status)~x1+factor(type),data=lung)

> anova(coxfit2,coxfit1)

Analysis of Deviance Table

 Cox model: response is  Surv(time, status)

 Model 1: ~ x1 + factor(type)

 Model 2: ~ x1 + x2 + x3 + factor(trt) + factor(type)

   loglik Chisq Df P(>|Chi|)

1 -88.143                  

2 -87.516 1.254  3    0.7401

anova()는 Cox Regression의 모형을 비교할 때 LRT(Likelihood ratio test)를 사용한다. 환자의 나이(x2), 진단 후 연구 참여시까지의 시간( x3)의 p-value는 0.7401으로 유의하지 않다.

출처: 보건정보데이터 분석 (이태림, 이재원, 김주한, 장대흥 공저), R을 이용한 누구나 하는 통계분석

비모수적 방법 non-parametric method

6.2.1생명표 방법 life table method

생명표 방법 개념 추가 설명

예] 협심증(angina pectoris)이 있는 2,418명의 남성들에 대한 생존자료이다. 이 생명표에서 경과기관에 따른 생존확률을 구해보자

1. 데이터 입력

> library(KMsurv)

> setwd('C:/Rwork')

> 협심증환자자료 = read.csv('협심증환자자료.csv')

> head(협심증환자자료)

  time censor freq

1  0.5      1  456

2  0.5      0    0

3  1.5      1  226

4  1.5      0   39

5  2.5      1  152

6  2.5      0   22

> attach(협심증환자자료)

2. lifetab 함수 연구

> lifetab

function (tis, ninit, nlost, nevent)

(이하 생략)

-> tis: 시간, 길이는 17, ninit = 전체 데이터, 2418명, nlost = 중도절단, nevent 사건발생 

3. nevent 자료 만들기

> which(censor==1) #사망 데이터

 [1]  1  3  5  7  9 11 13 15 17 19 21 23 25 27 29 31

> 집단.사망=협심증환자자료 [which(censor==1),]

> head(집단.사망)

   time censor freq

1   0.5      1  456

3   1.5      1  226

5   2.5      1  152

7   3.5      1  171

9   4.5      1  135

11  5.5      1  125

 4. nlost 자료 만들기

> which(censor==0) #절단 데이터

 [1]  2  4  6  8 10 12 14 16 18 20 22 24 26 28 30 32

> 집단.절단=협심증환자자료[which(censor==0),]

> head(집단.절단)

   time censor freq

2   0.5      0    0

4   1.5      0   39

6   2.5      0   22

8   3.5      0   23

10  4.5      0   24

12  5.5      0  107

5. ninit전체데이터 자료 만들기

> 사망자수=집단.사망[,3]

> 사망자수

 [1] 456 226 152 171 135 125  83  74  51  42  43  34  18   9   6   0

> 절단자수=집단.절단[,3]

> 절단자수

 [1]   0  39  22  23  24 107 133 102  68  64  45  53  33  27  23  30

> 합=사망자수+절단자수

> 합

 [1] 456 265 174 194 159 232 216 176 119 106  88  87  51  36  29  30

> sum(합) #2418

[1] 2418

6. tis 시간 변수 자료 만들기

> 년 = floor(집단.사망$time) #floor 내림값

> 년

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

> lt=length(년)

> lt

[1] 16

> length(년)

[1] 17

> 년[lt+1] = NA

> 년[lt+1]

[1] NA

7. 생명표 만들기 

> 생명표=lifetab(년,sum(합),절단자수,사망자수)

> 생명표

      nsubs nlost  nrisk nevent      surv        pdf     hazard     se.surv

0-1    2418     0 2418.0    456 1.0000000 0.18858561 0.20821918 0.000000000

1-2    1962    39 1942.5    226 0.8114144 0.09440394 0.12353102 0.007955134

2-3    1697    22 1686.0    152 0.7170105 0.06464151 0.09440994 0.009179397

3-4    1523    23 1511.5    171 0.6523689 0.07380423 0.11991585 0.009734736

4-5    1329    24 1317.0    135 0.5785647 0.05930618 0.10804322 0.010138361

5-6    1170   107 1116.5    125 0.5192585 0.05813463 0.11859583 0.010304216

6-7     938   133  871.5     83 0.4611239 0.04391656 0.10000000 0.010379949

7-8     722   102  671.0     74 0.4172073 0.04601094 0.11671924 0.010450930

8-9     546    68  512.0     51 0.3711964 0.03697464 0.10483042 0.010578887

9-10    427    64  395.0     42 0.3342218 0.03553750 0.11229947 0.010717477

10-11   321    45  298.5     43 0.2986843 0.04302654 0.15523466 0.010890741

11-12   233    53  206.5     34 0.2556577 0.04209376 0.17941953 0.011124244

12-13   146    33  129.5     18 0.2135639 0.02968456 0.14937759 0.011396799

13-14    95    27   81.5      9 0.1838794 0.02030570 0.11688312 0.011765989

14-15    59    23   47.5      6 0.1635737 0.02066194 0.13483146 0.012259921

15-NA    30    30   15.0      0 0.1429117         NA         NA 0.013300258

           se.pdf   se.hazard

0-1   0.007955134 0.009697769

1-2   0.005975178 0.008201472

2-3   0.005069200 0.007649121

3-4   0.005428013 0.009153696

4-5   0.004945997 0.009285301

5-6   0.005033980 0.010588867

6-7   0.004690538 0.010962697

7-8   0.005175094 0.013545211

8-9   0.005024599 0.014659017

9-10  0.005307615 0.017300846

10-11 0.006269963 0.023601647

11-12 0.006847514 0.030646128

12-13 0.006682743 0.035110295

13-14 0.006514794 0.038894448

14-15 0.008035120 0.054919485

15-NA          NA          NA

-> 5번째 열 surv는 생존 함수, 7번째 열 hazard는 위험 함수. 협심증 환자의 약 19%가 1년 이내, 약 28%가 2년 이내에 사망하는 것으로 추정된다. 따라서 협심증 환자가 2년 이상 생존할 확률은 약 72%가 된다는 것을 알 수 있다. 또한 5년 이상 생존할 확률은 약 52%이다.

#전체 환자 중 생존확율이 70% 되는 지점

> names(생명표)

 [1] "nsubs"     "nlost"     "nrisk"     "nevent"    "surv"      "pdf"     

 [7] "hazard"    "se.surv"   "se.pdf"    "se.hazard"

> which.min(abs(생명표$surv-0.7))

[1] 3

> 생명표[3,]

    nsubs nlost nrisk nevent      surv        pdf     hazard     se.surv

2-3  1697    22  1686    152 0.7170105 0.06464151 0.09440994 0.009179397

       se.pdf   se.hazard

2-3 0.0050692 0.007649121

-> abs:절대값, which.min 최소값, 또는 생명표에서 3번째 행을 보면 survival이 70%에 가장 가까이 접근했음을 알 수 있다. 2년 이상 생존할 확률은 72%

#전체 환자 중 생존확율이 50% 되는 지점

> which.min(abs(생명표$surv-0.5))

[1] 6

> 생명표[6,]

    nsubs nlost  nrisk nevent      surv        pdf    hazard    se.surv

5-6  1170   107 1116.5    125 0.5192585 0.05813463 0.1185958 0.01030422

        se.pdf  se.hazard

5-6 0.00503398 0.01058887

-> 5년이 지나면 surv 즉 생존할 확률이 약 52%가 됨을 알 수 있다.

8. 데이터 시각화 – 생존함수 추정치 그래프

 > plot(년[1:lt], 생명표[,5],type='s',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)

> plot(년[1:lt], 생명표[,5],type='o',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)         

9. 데이터 시각화 – 위험함수 추정치 그래프

> names(생명표)

 [1] "nsubs"     "nlost"     "nrisk"     "nevent"    "surv"      "pdf"     

 [7] "hazard"    "se.surv"   "se.pdf"    "se.hazard"

> mean(생명표$hazard)

[1] NA

> 생명표$hazard

 [1] 0.20821918 0.12353102 0.09440994 0.11991585 0.10804322 0.11859583

 [7] 0.10000000 0.11671924 0.10483042 0.11229947 0.15523466 0.17941953

[13] 0.14937759 0.11688312 0.13483146         NA

> mean(생명표$hazard,na.rm=T) #NA 값 제외한 평균

[1] 0.1294874

> plot(년[1:lt], 생명표[,7],type='s',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)

> abline(h=0.1294874)

> plot(년[1:lt], 생명표[,7],type='o',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)>

> abline(h=0.1294874)

 출처: 보건정보 데이터분석 (이태림, 이재원, 김주한, 장대흥 공저)