기본 콘텐츠로 건너뛰기

생존분석- 콕스 비례위험모형 (6)





 콕스 비례위험 모형 역시 여러 가지 모델을 지닐 수 있습니다. 이 가운데 가장 좋은 모형이 어떤 것인지를 알아내는 과정은 쉽지 않지만, 앞서 다른 분석 방법에서 본 것처럼 몇 가지 판단할 수 있는 방법이 있습니다. 콕스 비례위험 모형 역시 AIC (Akaike information criterion)를 기준으로 가장 적합한 모형을 찾는 방법이 있습니다. kidney 데이터에서 이를 찾아보겠습니다. 




library(survival)
library(MASS)

kd3<-coxph data="kidney)</span" disease="" factor="" sex="" status="=1)~age" ties="efron" time="" urv="">
summary(kd3)

model_selection<- kd3="" model="" selection="" span="" stepaic="" stepwise="" trace="TRUE)">
model_selection


 stepAIC를 통해 stepwise model selection이 가능한데, 당연히 AIC가 가장 작은 모델을 선택합니다. 만약 과정 없이 최종 모델만 보고 싶다면 trace + FALSE로 지정하면 됩니다. 


> kd3<-coxph data="kidney)</span" disease="" factor="" sex="" status="=1)~age" ties="efron" time="" urv="">
> summary(kd3)
Call:
coxph(formula = Surv(time, status == 1) ~ age + factor(sex) + 
    factor(disease), data = kidney, ties = "efron")

  n= 76, number of events= 58 

                        coef exp(coef)  se(coef)      z Pr(>|z|)    
age                 0.003181  1.003186  0.011146  0.285   0.7754    
factor(sex)2       -1.483137  0.226925  0.358230 -4.140 3.47e-05 ***
factor(disease)GN   0.087957  1.091941  0.406369  0.216   0.8286    
factor(disease)AN   0.350794  1.420195  0.399717  0.878   0.3802    
factor(disease)PKD -1.431108  0.239044  0.631109 -2.268   0.0234 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                   exp(coef) exp(-coef) lower .95 upper .95
age                   1.0032     0.9968   0.98151    1.0253
factor(sex)2          0.2269     4.4067   0.11245    0.4579
factor(disease)GN     1.0919     0.9158   0.49238    2.4216
factor(disease)AN     1.4202     0.7041   0.64880    3.1088
factor(disease)PKD    0.2390     4.1833   0.06939    0.8235

Concordance= 0.697  (se = 0.041 )
Rsquare= 0.207   (max possible= 0.993 )
Likelihood ratio test= 17.65  on 5 df,   p=0.003
Wald test            = 19.9  on 5 df,   p=0.001
Score (logrank) test = 20.13  on 5 df,   p=0.001

> model_selection<- kd3="" model="" selection="" span="" stepaic="" stepwise="" trace="TRUE)">
Start:  AIC=368.16
Surv(time, status == 1) ~ age + factor(sex) + factor(disease)

                  Df    AIC
- age              1 366.24
               368.16
- factor(disease)  3 372.69
- factor(sex)      1 381.12

Step:  AIC=366.24
Surv(time, status == 1) ~ factor(sex) + factor(disease)

                  Df    AIC
               366.24
- factor(disease)  3 370.74
- factor(sex)      1 379.15
> model_selection
Call:
coxph(formula = Surv(time, status == 1) ~ factor(sex) + factor(disease), 
    data = kidney, ties = "efron")

                      coef exp(coef) se(coef)      z        p
factor(sex)2       -1.4774    0.2282   0.3569 -4.140 3.48e-05
factor(disease)GN   0.1392    1.1494   0.3635  0.383   0.7017
factor(disease)AN   0.4132    1.5116   0.3360  1.230   0.2188
factor(disease)PKD -1.3671    0.2549   0.5889 -2.321   0.0203

Likelihood ratio test=17.56  on 4 df, p=0.001501
n= 76, number of events= 58 


 여기서는 sex와 disease만 넣은 변수가 가장 좋은 모델로 나왔습니다. 하지만 만약 우리가 진짜 알려는 독립 변수는 age이고 sex와 disease는 보정변수이기만 하다면 age를 뺄 수 없습니다. 그리고 아무래도 age가 중요한 영향인자라면 주된 독립 변수가 아니라도 뺄 수 없습니다. 선형 회귀나 로지스틱 회귀 모형과 마찬가지로 참고일 뿐이지 모델을 결정하는 절대적 방법은 아니지만, 여러 개의 변수가 존재할 때 참고할 수 있는 방법임은 분명합니다. 


 콕스 비례위험 모형은 이름처럼 각 변수의 위험도가 시간에 따라 비례해 증가하거나 감소한다고 가정하고 있습니다. 그러나 몇몇 변수들은 실제로는 그렇지 않을 가능성이 있습니다. 이럴 때 가장 간단한 해결책은 비례 가정을 만족시키지 않는 변수를 뺀 모형을 만드는 것이지만, 앞서 설명한 것처럼 그 변수가 빼기 어려운 중요한 보정 변수인 경우 고민에 빠지게 됩니다. 이 경우 층화(strata)가 한 가지 방법이 될 수 있습니다. 기본 개념은 비례 정도가 다른 집단에서 (예를 들어 남자와 여자, 노인과 청년 등)을 여러 개의 층으로 만들어 하나의 HR 값을 찾는 것입니다. 


 여기서는 age가 보려는 변수인데 주요 보정 변수인 sex가 비례 위험 가정을 만족시키지 않는 것처럼 보인다고 가정하겠습니다. 이 경우 sex를 모델에서 빼기 어렵기 때문에 성별에 따른 각각의 값을 구하거나 혹은 층화 콕스 비례 위험모형을 적용해야 합니다. strata()를 통해 적용할 수 있습니다. 


 #strata

kd3<-coxph age="" data="kidney)</span" factor="" sex="" status="=1)~" ties="efron" time="" urv="">
summary(kd3)

kd4<-coxph age="" data="kidney)</span" sex="" status="=1)~" strata="" ties="efron" time="" urv="">
summary(kd4)

> kd3<-coxph age="" data="kidney)</span" factor="" sex="" status="=1)~" ties="efron" time="" urv="">
> summary(kd3)
Call:
coxph(formula = Surv(time, status == 1) ~ age + factor(sex), 
    data = kidney, ties = "efron")

  n= 76, number of events= 58 

                  coef exp(coef)  se(coef)      z Pr(>|z|)   
age           0.002032  1.002034  0.009246  0.220  0.82607   
factor(sex)2 -0.829314  0.436349  0.298955 -2.774  0.00554 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

             exp(coef) exp(-coef) lower .95 upper .95
age             1.0020      0.998    0.9840     1.020
factor(sex)2    0.4363      2.292    0.2429     0.784

Concordance= 0.662  (se = 0.045 )
Rsquare= 0.089   (max possible= 0.993 )
Likelihood ratio test= 7.12  on 2 df,   p=0.03
Wald test            = 8.02  on 2 df,   p=0.02
Score (logrank) test = 8.45  on 2 df,   p=0.01

> kd4<-coxph age="" data="kidney)</span" sex="" status="=1)~" strata="" ties="efron" time="" urv="">
> summary(kd4)
Call:
coxph(formula = Surv(time, status == 1) ~ age + strata(sex), 
    data = kidney, ties = "efron")

  n= 76, number of events= 58 

        coef exp(coef) se(coef)     z Pr(>|z|)
age 0.008007  1.008039 0.009253 0.865    0.387

    exp(coef) exp(-coef) lower .95 upper .95
age     1.008      0.992    0.9899     1.026

Concordance= 0.547  (se = 0.048 )
Rsquare= 0.01   (max possible= 0.982 )
Likelihood ratio test= 0.77  on 1 df,   p=0.4
Wald test            = 0.75  on 1 df,   p=0.4
Score (logrank) test = 0.75  on 1 df,   p=0.4


 결과를 보면 층화를 한 경우와 아닌 경우 모두 통계적 유의성이 없지만, 비례 계수가 1.002에서 1.008로 바뀌는 것을 볼 수 있습니다. 만약 sex가 비례 위험 가정을 만족하지 않는다면 층화를 통해 이 문제를 극복하고 더 정확한 값을 구할 수 있습니다. strata()를 적용한 변수에 대해서는 HR 값이 얻어지지 않는다는 점도 주목해야 합니다. 당연히 값이 나올 수 없습니다. 다만 이 결과를 믿기 위해서는 sex와 age간의 중요한 교호작용(상호작용)이 있는지를 검증해야 합니다. 


kd5<-coxph age="" data="kidney)</span" sex="" status="=1)~" strata="" ties="efron" time="" urv="">
summary(kd5)

> kd5<-coxph age="" data="kidney)</span" sex="" status="=1)~" strata="" ties="efron" time="" urv="">
> summary(kd5)
Call:
coxph(formula = Surv(time, status == 1) ~ age * strata(sex), 
    data = kidney, ties = "efron")

  n= 76, number of events= 58 

                         coef exp(coef) se(coef)      z Pr(>|z|)
age                  -0.01339   0.98670  0.01692 -0.791    0.429
age:strata(sex)sex=2  0.02893   1.02935  0.02009  1.440    0.150

                     exp(coef) exp(-coef) lower .95 upper .95
age                     0.9867     1.0135    0.9545     1.020
age:strata(sex)sex=2    1.0294     0.9715    0.9896     1.071

Concordance= 0.559  (se = 0.047 )
Rsquare= 0.035   (max possible= 0.982 )
Likelihood ratio test= 2.74  on 2 df,   p=0.3
Wald test            = 2.68  on 2 df,   p=0.3
Score (logrank) test = 2.71  on 2 df,   p=0.3


 P = 0.15로 유의한 교호작용은 없으며 층화 모형이 의미 있다고 할 수 있습니다. 만약 다른 변수도 층화를 통해 구해야 한다면 같은 과정을 반복하면 됩니다. 


kd3<-coxph age="" data="kidney)</span" disease="" factor="" sex="" status="=1)~" strata="" ties="efron" time="" urv="">
summary(kd3)

kd4<-coxph age="" data="kidney)</span" disease="" sex="" status="=1)~" strata="" ties="efron" time="" urv="">
summary(kd4)

kd5<-coxph age="" data="kidney)</span" disease="" sex="" status="=1)~" strata="" ties="efron" time="" urv="">
summary(kd5)


 그리고 만약 연속변수가 비례 위험 모형을 만족시키지 않는 경우 변수를 몇 개의 그룹으로 쪼개 문제를 피할 수도 있습니다. 예를 들어 나이 60세를 기준으로 나누는 등입니다. 그래도 비례 위험 가정이 만족되지 않으면 역시 층화를 시도하거나 그룹을 나눠 분석할 수 있습니다. 

댓글

이 블로그의 인기 게시물

100 테슬라급 자기장 도달

 미국의 로스 알라모스 국립 연구소 (Los Alamos National Laboratory) 에서 과학자들이 지금까지 인간이 개발한 가장 강력한 자기장을 발생시키는 장치 개발에 도전하고 있습니다. 자기장의 세기를 나타내는 방법으로 자기력선의 밀도를 나타내기 위해 단위 면적당 자기력선의 수를 표시하는 단위인 테슬라 (T) 가 있습니다. (1T = 1Wb/㎡  웨버 (Wb) 는 자속의 단위)   의료용으로 사용되는 초전도체를 이용한 MRI 의 경우 1.5 - 3 테슬라급의 강력한 자기장으로 인체 내부를 볼 수 있게 만들지만 과학 연구용으로 이보다 더 강력한 자기장이 필요할 수 있습니다. 최근에 등장한 90 테슬라급 자기장에 이어 이번에 로스 알라모스 국립 연구소에서는 100 테슬라급인 100.75 T 를 실현 했다고 합니다.   이를 구현한 것은 18000 파운드 (8.16 톤 정도) 의 코일과 여기에 에너지를 공급할 1200 메가줄 (Megajoule) 급 모터 제네레이터등의 설비입니다.  (   The 1,200-megajoule motor generator that powers the magnetic pulse.  )  이와 같은 연구를 통해 알아내고자 하는 것은  Quantum Phase transitions and new ultra high field magnetic states Electronic Structure determination Topologically protected states of matter  로 요약할 수 있다고 합니다.   아무튼 수 T 급 MRI 만 해도 자기장의 힘이 엄청난데 100 T 라니 엄청난 자기장이네요. 이는 지구 자기장 세기보다 200만배 강력한 (물론 좁은 범위에서 작용하는 자기장이라 지구 전체...

고대 양서류 이야기 (2) - 악어를 닮은 거대 양서류들

  페름기에는 다양한 양막류가 진화해서 앞서 소개한 육상형 템노스폰딜리는 점차 설 자리를 잃게 됩니다. 하지만 양서류는 본래 자신의 서식지인 물과 습지에서 여전히 번성을 누렸습니다. 당시에는 악어류 같은 대형 양서형 파충류도 없던 시절이었기 때문에 이와 비슷한 생태학적 지위는 여전히 양서류의 몫이었습니다. 이들에 대한 이야기는 제 책인 포식자에서 비교적 간단히 다뤘는데, 오늘은 여기에 대한 보충 설명입니다.  책 정보:  http://book.naver.com/bookdb/book_detail.nhn?bid=13347200 Yes 24:  http://www.yes24.com/24/goods/58772859 11번가:  http://books.11st.co.kr/product/SellerProductDetail.tmall?method=getSellerProductDetail&prdNo=1977867160 알라딘:  http://www.aladin.co.kr/shop/wproduct.aspx?ItemId=134877825 교보문고:  http://www.kyobobook.co.kr/product/detailViewKor.laf?ejkGb=KOR&mallGb=KOR&barcode=9788970447988&orderClick=LAG&Kc= 인터파크 :  http://book.interpark.com/product/BookDisplay.do?_method=detail&sc.prdNo=279593764&sc.saNo=003002003&bid1=search_auto&bid2=detail&bid3=prd_img&bid4=001 영풍문고:  http://www.ypbooks.co.kr/book.yp?bookcd=100843205&gubun=NV ...

이빨이 다시 진화한 개구리

  ( CT scans of Gastrotheca guentheri skulls revealed what appeared to be identical rows of teeth on both the upper and lower jaws, which researchers later confirmed through dissection. Credit: Florida Museum/Daniel Paluh )  개구리는 2억년 전 진화 과정에서 이빨을 잃어버리고 큰 턱과 혀를 이용해 곤충 같은 작은 먹이를 잡아 먹는 방향으로 진화했습니다. 파충류나 포유류 같은 다른 사지류와의 경쟁에서 밀려 양서류가 쇠퇴하고 멸종하던 시기에도 개구리는 여전히 생존할 수 있었던 비결입니다.   이빨이 없는 덕분에 큰 혀를 발사하기 편해졌을지는 모르지만, 이빨이 없으면 종종 불편할 때가 있습니다. 씹는 대신 삼키기 때문에 음식을 씹을 수 없다는 점은 문제되지 않지만, 필사적으로 달아나려는 먹이를 잡기 힘들기 때문입니다. 이런 이유 때문에 일부 개구리는 이빨 같이 보이는 엄니 (fang)을 지니고 있으나 이는 사라진 이빨이 다시 난 것이 아니라 다른 부분이 진화한 것입니다.   이는 돌로의 법칙 ( Dollo's Law )으로 알려져 있습니다. 진화 과정에서 퇴화한 부분이 다시 생겨나지 않으며 대신 필요하면 다른 부분이 진화해 그 역할을 대신한다는 것입니다. 예를 들어 아가미가 사라진 사지 동물은 다시 물에 들어온다고 해도 아가미가 다시 생기진 않습니다. 대신 고래처럼 폐가 커져서 그 기능을 대신하게 됩니다.   하지만 모든 법칙엔 예외가 있기 마련입니다. 남미에서 발견된 멸종 위기 개구리 중 하나인 구엔터 유대류 개구리 ( Gastrotheca guentheri, Guenther's marsupial frog, dentate marsupial frog)는 완전한 형태의 이빨을 지니고 있습니다. 참고로 유대류 개구리라는 명칭은...