데이터마이너를 꿈꾸며

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을 이용한 누구나 하는 통계분석