기본 콘텐츠로 건너뛰기

생존분석- 콕스 비례위험모형 (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세를 기준으로 나누는 등입니다. 그래도 비례 위험 가정이 만족되지 않으면 역시 층화를 시도하거나 그룹을 나눠 분석할 수 있습니다. 

댓글

이 블로그의 인기 게시물

벨 V-280 Valor 시험 비행 성공

( The V-280 Valor flew for the first time at Bell Helicopter's Amarillo Assembly Center in Texas(Credit: Bell Helicopter/YouTube) )  앞서 소개드린 V-280 발러가 첫 번째 비행 테스트에 성공했다는 소식입니다. V-22 오스프리의 소형화 버전이라고 할 수 있는 V-280 발러는  미 육군의 차세대 헬기 사업인 Future Vertical Lift (FVL)에 입찰을 시도하는 틸트로터기로 현재 미 육군이 주력으로 사용하는 블랙호크 헬기와 비슷한 체급입니다. 다만 틸트로터기인 만큼 최고 속도나 항속 거리면에서 더 유리합니다. 스펙은 이전 포스트를 참조해 주시기   이전 포스트:  https://blog.naver.com/jjy0501/221115245986  (동영상)   V-280 발러는 틸트로터기의 더 대중화 될 수 있을지를 검증하는 중요한 무대가 될 것입니다. V-22 오스프리의 경우 복잡한 구조로 인해 가격이 너무 비싸져서 사실 미국은 몰라도 그 동맹국에 널리 도입되기는 어려운 부분이 있습니다. V-280 역시 가격이 아주 저렴할 것 같지는 않지만, 좀 더 합리적인 대안은 될 수 있을 것 같습니다. 만약 성공적인 결과가 나오면 한국을 포함한 미국의 동맹국에서 도입을 검토할 수 있을지 모르겠다는 생각입니다.   참고  https://newatlas.com/bell-v-280-valor-maiden-flight/52663/

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 ...