생존분석 - Competing Risk를 고려한 생존 분석

Competing Risk를 고려한 생존 분석

Competing risk 분석은 생존분석에서 관심 event 외의 다른 event가 있을 때, 분석하는 방법입니다. 예를 들어, 특정 약이 폐암으로 인한 사망을 줄이는가에 대해 관심을 갖고 연구하려고 하는데 심혈관 질환, 뇌혈관 질환 등으로 사망한 경우 Competing risk 분석을 통해 이 약이 폐암으로 인한 사망은 줄여도, 다른 사망에 대해서는 반대의 효과를 줄 수 있기 때문에 이를 확인할 필요가 있습니다.

Competing risk는 기본적인 생존분석 모형의 확장된 버전입니다. 데이터의 형태는 기본 생존 분석의 [time to event, event( censoring, event 를 나타내는 변수)]를 따르지만, vent가 censoring, event 0/1 로 나누어진 것이 아니라, 3개 이상의 범주형 변수를 갖게 됩니다. 

R의 경우 timereg, cmprsk 를 이용하여 competing risk 분석을 수행할 수 있습니다. 

먼저 데이터를 읽어옵니다. 

fol <- read.csv("data.csv", header=TRUE)

age	path1	hgb	ldh	clinstg	blktxcat	relsite	ch	rt	survtime	stat	dftime	dfcens	resp	stnum
56	NH04	140	NA	2	1	B		Y	0.6981520	1	0.238193018	1	CR	1
36	NH02	130	NA	2	1	D		Y	14.5023956	1	12.418891170	1	CR	2
39	NH02	140	NA	2	3		Y	Y	4.9144422	1	0.002737851	1	NR	3
37	NH03	140	NA	1	1			Y	15.6851472	1	15.685147159	1	CR	4
61	NH04	110	NA	2	2			Y	0.2354552	1	0.002737851	1	NR	5
69	NH02	120	NA	1	1			Y	8.4188912	1	8.418891170	1	CR	6

이 데이터에는 1,2 두 가지의 사망원인이 있습니다. 

cause
  0   1   2 
193 272  76 

1. 관심 사건 외 다른 이벤트들을 censoring으로 취급하여 분석

library(data.table)
fol <- setDT(fol)

fol$y1 = fol$cause
fol$y2 = fol$cause
fol[cause==2, y1:=0]
fol[cause==1, y2:=0]
fol[y2==2 , y2:=1]

원래는 0이 censoring 인데, 관심요인이 1일 때는 2를 censoring, 관심요인이 2일 때는 1을 censoring으로 바꾸어주는 코드입니다. 이 때, data.table으로 dataframe으로 변환한 후, := 키워드를 이용하였습니다.

① Cause = 1에 대한 Model (Cause=2 를 Censoring 취급)

=> age, hgb, clinstg를 보정해서 blktxcat의 effect 파악

library(survival)

cause1 <- Surv(fol$dftime, fol$y1)
cox1 <- coxph(cause1 ~ age+hgb+blktxcat+clinstg, data=fol)

Call:
coxph(formula = y0 ~ age + hgb + blktxcat + clinstg, data = fol)

  n= 541, number of events= 272 

             coef exp(coef) se(coef)     z Pr(>|z|)    
age      0.023784  1.024069 0.004733 5.025 5.02e-07 ***
hgb      0.003078  1.003083 0.004075 0.755 0.449974    
blktxcat 0.207383  1.230454 0.057329 3.617 0.000298 ***
clinstg  0.329831  1.390733 0.138171 2.387 0.016981 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
age          1.024     0.9765    1.0146     1.034
hgb          1.003     0.9969    0.9951     1.011
blktxcat     1.230     0.8127    1.0997     1.377
clinstg      1.391     0.7190    1.0608     1.823

Concordance= 0.623  (se = 0.019 )
Rsquare= 0.082   (max possible= 0.997 )
Likelihood ratio test= 46.23  on 4 df,   p=2e-09
Wald test            = 45.16  on 4 df,   p=4e-09
Score (logrank) test = 45.84  on 4 df,   p=3e-09

② Cause = 2에 대한 Model (Cause=1 를 Censoring 취급)

=> age, hgb, clinstg를 보정해서 blktxcat의 effect 파악

cause2 <- Surv(fol$dftime, fol$y2)

cox2 <- coxph(cause2 ~ age+hgb+blktxcat+clinstg, data=fol)

summary(cox2)

Call:
coxph(formula = cause2 ~ age + hgb + blktxcat + clinstg, data = fol)

  n= 541, number of events= 76 

             coef exp(coef) se(coef)     z Pr(>|z|)    
age      0.087125  1.091033 0.011440 7.616 2.62e-14 ***
hgb      0.005945  1.005962 0.008452 0.703    0.482    
blktxcat 0.164626  1.178952 0.109133 1.508    0.131    
clinstg  0.353410  1.423915 0.285801 1.237    0.216    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
age          1.091     0.9166    1.0668     1.116
hgb          1.006     0.9941    0.9894     1.023
blktxcat     1.179     0.8482    0.9519     1.460
clinstg      1.424     0.7023    0.8132     2.493

Concordance= 0.764  (se = 0.039 )
Rsquare= 0.131   (max possible= 0.759 )
Likelihood ratio test= 75.68  on 4 df,   p=1e-15
Wald test            = 62.84  on 4 df,   p=7e-13
Score (logrank) test = 67.68  on 4 df,   p=7e-14

2. Cumulative incidence curve를 통한 분석

Cumulative incidence curve의 정의 

개념 : Cumulative incidence curve란 competing risk를 고려한 1-생존함수 의 추정값 이다. 

$$ CIC_c(t) = P(T_c \lt t) $$

위 식에서 CIC란 어떤 사람의 c라는 event의 발생 시각인 Tc가 t보다 작을 확률입니다. 만약 competing risk가 없다면 이것은  간단하게 계산됩니다. 일반적으로 competring risk 가 없는 상황에서는 이것은 1 - (Kaplen-Meier 생존함수)로 추정하게 됩니다. 하지만 competing risk가 있는 경우, 다른 원인으로도 사람이 죽을 수 있기 때문에 Kaplen-Meier 방법을 그대로 사용하면 특정 시점 t 이전에 event 가 발생했을 확률이 과대추정됩니다. 

CIC 방법에서는 proportion of dying을 아래와 같은 방식으로 추정하게 됩니다. 

$$ CIC_c(t) = \int_{0}^{t} h_c(u)S(u) du $$

이 방식을 직관적으로 이해하면 Incidence를 모든 t에 대해서 더해준다고 이해할 수 있는데 incidence는 특정 시점에서의 cause-specific hazard와 생존함수를 곱해준 값입니다. 이 때, 생존함수는 인구집단에서의 생존함수입니다. 만약 이 생존 함수를 해당 cause에 대한 생존함수를 사용한다면 1-KM과 같은 값이 나오게 됩니다. 하지만 인구집단에서의 생존함수를 사용하여 competing risk를 고려하는 것으로 이해할 수 있습니다. 

Cumulative incidence curve는 cmprsk 라는 라이브러리를 통해 분석할 수 있습니다. Cumulative incidence는 cause-specific hazard function으로 계산되며, proportion of dying을 계산해줍니다. 이 때, 관심 요인 blktxcat을 지정해줍니다. blktxcat을 약의 종류라고 생각을하고 0,1,2,3 총 4개의 약의 종류가 있는 것입니다. 이 4개의 약이 1,2 두 가지 사망 원인에 어떻게 작용하는지를 아래 표를 통해 확인할 수 있습니다=> 각각의 strata 별로 proportion of dying의 예측값을 보여주고 있습니다. 

예를 들어, 25달이 지난 후, 사망원인 1을 보면, 0의 약을 쓴 사람은 52%가 죽은 반면, 3 약을 쓴 사람은 66%가 죽었습니다.  

library(cmprsk)
CIC <- cuminc(fol$dftime, fol$cause, fol$blktxcat)
CIC

Tests:
       stat           pv df
1 17.189363 0.0006461085  3
2  4.596607 0.2038333549  3
Estimates and Variances:
$est
             5         10         15        20        25        30
0 1 0.28456914 0.40301775 0.47654623 0.5235252 0.5235252 0.5235252
1 1 0.36620667 0.44808409 0.47719815 0.5096313 0.5541112        NA
2 1 0.42528724 0.54412907 0.59678622 0.5967862 0.5967862        NA
3 1 0.50604626 0.63733974 0.66670849 0.6667085 0.6667085        NA
0 2 0.02540834 0.08048914 0.15465458 0.1546546 0.1546546 0.2619280
1 2 0.11604642 0.16329519 0.18920323 0.2309031 0.2679696        NA
2 2 0.05245607 0.07941283 0.13288099 0.2003667 0.2409361        NA
3 2 0.03960396 0.06897272 0.09834147 0.1374998 0.2353957        NA

$var
               5           10          15          20          25          30
0 1 0.0010438691 0.0013684404 0.001665324 0.002070841 0.002070841 0.002070841
1 1 0.0022693095 0.0025757378 0.002723652 0.003459139 0.004840443          NA
2 1 0.0018165618 0.0021166396 0.002288543 0.002288543 0.002288543          NA
3 1 0.0025253539 0.0027149325 0.002701947 0.002701947 0.002701947          NA
0 2 0.0001266488 0.0004372677 0.001007219 0.001007219 0.001007219 0.006742183
1 2 0.0010075372 0.0014431850 0.001683726 0.002346285 0.003481490          NA
2 2 0.0003774854 0.0005948500 0.001226766 0.002385591 0.003610060          NA
3 2 0.0003819491 0.0007880136 0.001162941 0.001791598 0.011245328          NA

plot(CIC, xlab="month", ylab="proportion dying")

위와 같이 그림도 그릴 수 있습니다. 위쪽에 위치한 4개의 그래프가 사망원인 1에 대한 4가지 약의 효과 그래프이고, 아래쪽에 있는 4개의 그래프가 사망원인 2에 대한 것입니다. 전반적으로 사망원인 1로 죽는 사람이 많다는 것을 확인할 수 있습니다. 또한 약의 경우 사망원인 1에 대해 어느 정도 연관성이 있는 것으로 보입니다. 각 4개의 약에 대해 proportion of dying에 차이를 보여줍니다. 하지만 사망원인 2에 대해서는 약과 별로 상관 없는 것을 확인할 수 있습니다. 

3. Conditional probability function

library(Cprob)

CPC <- cpf(Hist(dftime, cause)~blktxcat, data=fol, failcode=1)
CPC
plot(CPC, xlab="month", ylab="proportion dying")

Conditional probability function도 마찬가지로 그릴 수 있습니다. 이 때, failcode로 관심 event를 지정해줍니다. 

4. Proportional-odds Model for the Conditional Probability Function

CPC_reg=cpfpo(Hist(dftime, cause)~blktxcat+age, data=fol, failcode=1, tis=seq(0, 70, 5))
CPC_reg

Warning message:

“glm.fit: algorithm did not converge”

Test for non-significant effect 

                  stat       pvalue

(Intercept) -26.769122 0.000000e+00

blktxcat      3.952040 7.748788e-05

age           5.041975 4.607511e-07

 Test for constant fit 

                   coef     exp.coef     SE.coef       stat     pvalue

(Intercept) -9.51827834 7.349609e-05 0.355569317 -24.964365 0.00000000

blktxcat     0.27674405 1.318829e+00 0.070025626   1.790873 0.07331363

age          0.02991248 1.030364e+00 0.005932691  10.309364 0.00000000

Conditional probability function의 proportional odds 모델을 이용하면 regression을 통해 각 factor들의 odds ratio를 구할 수 있습니다. 

data.csv

데이터 출처

Thomas H. Scheike, Mei-Jie Zhang

Title: Analyzing Competing Risk Data Using the R timereg Package

Journal of statistical software

https://www.jstatsoft.org/article/view/v038i02