기본 콘텐츠로 건너뛰기

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

댓글

이 블로그의 인기 게시물

중등도 코로나 19 환자에서 제한적인 효과만을 확인한 렘데시비르

   현재 렘데시비르는 인공 호흡기 치료가 필요한 환자 등 중증 코로나 19 환자의 치료제로 사용되고 있습니다. 다만 공급이 부족해서 갑자가 중증 코로나 19 환자가 늘어날 경우 수급 문제가 불거질 수밖에 없습니다. 따라서 이를 중등도 코로나 19 환자 (경증이나 무증상은 당연히 치료가 필요없거나 대증 요법 정도) 사용한다는 것은 지금 당장에 급한 연구는 아니지만, 나중에 수급량이 많아지고 가격이 낮아질 때를 생각하면 의미 있는 연구 결과입니다.   렘데시비르 제조사인 길리어드 사이언스의 다아애나 브레이너드 박사 ( Dr. Diana Brainard of Gilead Sciences )는 병원에 입원이 필요했던 중등도 코로나 19 환자 584명을 렘데시비르 5일, 10일 투여 실험군과 위약을 투여한 대조군으로 나눠 임상 경과, 사망 등 중요한 합병증 발생 등에서 차이를 보이는지 조사했습니다. 연구 대상자의 40% 이상은 당뇨 같은 만성 질환이 있었으며 평균 연령은 57세였습니다.   연구 결과 대조군의 2%, 5일 투여군의 2%, 10일 투여군의 1%가 사망했으나 사망자 수가 4명, 3명, 2명 정도로 통계적 차이를 판단하기에는 사망자 수가 적어 결론을 내리기 어려웠습니다. 연구팀은 5일 투여군에서 임상 경과가 다소 호전되는 것을 판단했으나 10일 투여군에서는 통계적 차이를 검증할 수 없었습니다. 결론적으로 말해 중등도 코로나 19 환자에서 렘데시비르 투여는 확실한 임상적 이득이 없었습니다.     5일 투여량만 해도 가격이 3120달러라는 점을 감안할 때 중등도 환자에서 렘데시비르 투여는 설령 공급이 원활하다고 해도 적절한 치료법이 아닌 것으로 판단됩니다. 더구나 지금처럼 공급이 불안해 중증 환자 투여 물량 확보가 아쉬운 상황에서는 말할 것도 없습니다.  이번 연구는 렘데시비르보다 효과가 뛰어나면서 저렴한 차세대 항바이러스 치료제의 필요성을 다시 보여준 것으로 생각됩니다....

178Tb/s - 광섬유 통신 신기록 수립

  ( Lidia Galdino, lead researcher on the study. Credit: James Tye / UCL )  현재 가정에 들어오는 초고속 광통신은 10Gb/s 정도가 가장 빠른 편입니다. 그러나 사실 광섬유를 통해서 보낼 수 있는 정보량은 이것보다 훨씬 많습니다. 현재 개발 중인 기술 가운데서는 이미 1Tb/s 이상이 속도를 지닌 것들이 많아 미래에는 테라 인터넷 시대가 가능할 것으로 예상됩니다. 그런 가운데 영국 유니버시티 칼리지 런던 (UCL)의 리디아 갈디노( Lidia Galdino )가 이끄는 연구팀이 178Tb/s라는 신기록을 세웠다는 소식입니다.   광섬유 통신을 위한 주파수는 4.5THz까지 높일 수 있습니다. 최근에 개발된 신기술은 9THz 벽도 돌파했습니다. 그러나 연구팀은 16.8THz라는 새로운 영역에 도달했습니다. 주파수가 높을 수록 더 많은 정보를 전송할 수 있습니다. 하지만 16.8THz의 주파수로도 178Tb/s의 데이터 전송 속도는 사실 불가능합니다. 연구팀은  Geometric Shaping (GS) 방식으로 레이저 신호의 페이즈, 밝기, 편광 수준을 조절해 여러 신호가 서로 간섭 없이 전달되게 만들어 이런 신기록을 수립했습니다.   물론 새로운 기술이 항상 그렇듯이 바로 상용화가 가능한 것은 아닙니다. 특히 광섬유는 장거리 정보 전송이 가능해야 할 뿐 아니라 오류 없이 장시간 작동이 가능해야 합니다. 자꾸 에러가 생기거나 고장난다면 아무리 속도가 빨라도 무용지물입니다. 이 단계를 극복하고 경제성까지 갖춰 실용화 되는 경우는 사실 소수에 지나지 않습니다.   하지만 경쟁적으로 Tb 급 대역폭을 지닌 광통신 기술이 개발되고 있어 테라바이트 급 정보도 순식간에 전송할 수 있는 시대가 올 것은 분명합니다. 과거 전화선을 이용한 56K 모뎀에서 지금의 광섬유 기반 초고속 인터넷이 온 것처럼 속도는 계속 빨라질 수밖에 없기 때문입니다...

코로나 19 감염 위험성이 큰 동물은?

  ( A new genomic study ranks the potential of the SARS-CoV-2 spike protein to bind to the ACE2 receptor site in 410 vertebrate animals. Old World primates and great apes, which have identical amino acids at the binding site as humans, are predicted to have a very high propensity for binding ACE2 and are likely susceptible to SARS-CoV-2 infection. Credit: Matt Verdolivo/UC Davis )    코로나 19를 일으키는  SARS-CoV-2는 사스 및 메르스와 마찬가지로 박쥐를 자연 숙주로 삼는 베타 코로나바이러스 계통입니다. 따라서 본래 박쥐에 살던 바이러스가 중간 단계를 거쳐 사람으로 넘어온 것으로 추정됩니다. 그런데 이 이야기는 SARS-CoV-2 역시 사람에만 감염되는 것이 아니라 종을 넘어 다른 포유류나 조류, 혹은 척추동물에 감염될 수 있음을 시사합니다. 코로나 바이러스가 침투하는 경로인 ACE2 수용체는 다른 척추동물에서도 볼 수 있기 때문입니다. 다만 SARS-CoV-2와 ACE2 수용체의 결합 수준은 동물마다 큰 차이가 있어 감염 위험도는 모두 다를 수밖에 없습니다.   캘리포니아 대학 데이비스 캠퍼스 (UC Davis)의 조안나 다마스 (Joana Damas)와 그 동료들은 410종의 척추동물에서 ACE2 수용체 유전자를 조사해 코로나 19 감염 위험도를 조사했습니다. 연구팀에 따르면 ACE2 수용체 전체가 아니라 25개의 아미노산 배열이 특히 결합력에 큰 영향을 미쳐 감염 가능성을 결정합니다.   연구 결과 쉽게 예상할 수 있듯이 사람과 가까운 침팬치, 고릴라, 보노보...