Comparison of Models for Competing Risk Analysis

[Theoretical background]

Competing Risk 분석에서 모델의 예측 성능을 평가하기 위해 일반적으로 사용되는 지표는 다음과 같다.

1. Harrell’s C-index

개념:
관찰된 사건 발생 순서와 모델의 예측 위험 순서가 얼마나 잘 일치하는지를 나타낸다.
범위는 0부터 1이고, 0.5는 무작위 예측과 동일한 성능이며 1에 가까울수록 예측 성능이 좋다.

계산 방법:
Fine-Gray 모델을 이용하여 Competing Risk 데이터를 일반적인 생존 분석 형태로 변환한 후
Cox 비례위험 모델로 분석한 결과를 바탕으로 계산한다.
R에서는 survival::finegray()로 데이터를 변환하고, survival::coxph()와 survival::concordance()로 C-index를 구한다.

2. Wolbers’ C-index

개념:
Harrell’s C-index 개념을 확장하여 Competing Risk 상황에서 특정 원인의 사건 발생을 예측할 때 사용한다.

계산 방법:
Fine-Gray 모델(FGR)을 직접 적용하여 사건 특화된 C-index를 구한다.
Bootstrapping을 이용하여 반복 표본추출로 신뢰구간을 얻는다.
R에서는 riskRegression::FGR()과 pec::cindex()를 사용한다.

3. AUC (Area Under Curve)

개념:
특정 시점에서 사건 발생 여부를 이진 분류로 간주하고, 예측의 민감도와 특이도를 종합적으로 평가한다.
1에 가까울수록 예측 성능이 뛰어나고, 0.5에 가까울수록 무작위와 유사하다.

계산 방법:
riskRegression::Score() 함수를 이용해 특정 시점(예: 5년, 10년)의 AUC와 신뢰구간을 계산한다.

4. Brier Score

개념:
예측 확률과 실제 관찰된 사건의 발생 여부 간의 차이를 측정하는 지표이다.
낮을수록 성능이 우수하고, 0에 가까울수록 완벽한 예측이다.

계산 방법:
AUC와 같은 방식으로 riskRegression::Score() 함수로 특정 시점의 Brier Score와 신뢰구간을 구한다.

[Preprocess]

melanoma 데이터에서 필요한 변수(sex, age, thickness, ulcer)를 추출하고, status 변수는 censor=0인 Competing Risk 형태로 재구성한다. 시간을 연 단위로 변환하여 분석을 위한 데이터(melanoma_dt)를 준비한다.
status: 1=melanoma 사망, 2=생존, 3=melanoma 외 사망
status_competing: 0=생존, 1=melanoma 사망, 2=melanoma 외 사망

suppressPackageStartupMessages({
  library(data.table);library(survival);library(boot);library(magrittr);library(riskRegression);library(prodlim);library(pec);library(rmarkdown)
})

melanoma_dt <- as.data.table(melanoma) %>% .[,c("sex","ulcer"):=lapply(.SD,as.factor),.SDcols=c("sex","ulcer")] %>% 
  .[,status_competing:=factor(ifelse(status==1,1,ifelse(status==2,0,2)), levels=0:2)] %>% 
  .[,time_Y:=time/365.25] %>% .[,-c("year","status","time")]

model_var_list <- list(c("sex","age"), c("sex","age","thickness"), c("sex","age","ulcer"), c("sex","age","thickness","ulcer"))

Hist_formula_list <- list( Hist(time_Y, status_competing) ~ sex+age,
                           Hist(time_Y, status_competing) ~ sex+age+thickness,
                           Hist(time_Y, status_competing) ~ sex+age+ulcer,
                           Hist(time_Y, status_competing) ~ sex+age+thickness+ulcer )

[Each model]

4가지 변수 조합의 모델을 설정하고 각각의 Harrell C-index, Wolbers C-index, AUC, Brier score를 계산한다. Wolbers’ C-index의 신뢰구간은 bootstrap 방법을 이용해 반복 표본추출로 계산한다.

invisible(capture.output({
  table_1 <- lapply (1:length(model_var_list), function(i) { 
  # Harrell_C_index_info
  fgr_data_for_Harrell_C_index <- survival::finegray(as.formula(paste0("Surv(time_Y, status_competing) ~ ",paste(model_var_list[[i]],collapse="+"))), data=melanoma_dt, etype=1)
  fgr_model_for_Harrell_C_index <- coxph(as.formula(paste0("Surv(fgstart, fgstop, fgstatus) ~ ",paste(model_var_list[[i]],collapse="+"))), data=fgr_data_for_Harrell_C_index)
  Harrell_C_index_info <- survival::concordance(fgr_model_for_Harrell_C_index)
  
  # Wolbers_C_index_info
  fgr_model_for_Wolbers_C_index <- riskRegression::FGR(Hist_formula_list[[i]], data=melanoma_dt, cause=1)
  Wolbers_C_index_info <- pec::cindex(object=fgr_model_for_Wolbers_C_index, formula=Hist_formula_list[[i]], data=melanoma_dt, cause=1, confInt=T, verbose=F)
  
  cindex_values <- numeric()
  for(b in 1:20){
    indices <- sample(seq_len(nrow(melanoma_dt)), size=nrow(melanoma_dt), replace=T)
    boot_data <- melanoma_dt[indices]
    fgr_model_for_Wolbers_C_index_boot <- riskRegression::FGR(Hist_formula_list[[i]], data=boot_data, cause=1)
    Wolbers_C_index_info_boot <- pec::cindex(object=fgr_model_for_Wolbers_C_index_boot, formula=Hist_formula_list[[i]], data=boot_data, cause=1, verbose=F)
    cindex_values[[b]] <- Wolbers_C_index_info_boot$AppCindex$FGR
  }
  Wolbers_C_index_lower <- quantile(cindex_values, probs=0.05/2)
  Wolbers_C_index_upper <- quantile(cindex_values, probs=(1-(0.05/2)))

  # AUC_and_Brier_info
  fgr_model_for_AUC_and_Brier <- riskRegression::FGR(Hist_formula_list[[i]], data=melanoma_dt, cause=1)
  AUC_and_Brier_info <- riskRegression::Score(list(fgr_model_for_AUC_and_Brier), formula=Hist_formula_list[[i]] , data=melanoma_dt, times=c(5,10), null.model=F)

  # return
  data.table( "Model"=LETTERS[i],
              "Harrell_C_index"=paste0(sprintf("%.04f",Harrell_C_index_info$concordance),"(",
                                       sprintf("%.04f",Harrell_C_index_info$concordance-qnorm(0.975)*sqrt(Harrell_C_index_info$var)),"-",
                                       sprintf("%.04f",Harrell_C_index_info$concordance+qnorm(0.975)*sqrt(Harrell_C_index_info$var)),")"),
              "Wolbers_C_index"= paste0(sprintf("%.04f",Wolbers_C_index_info$AppCindex$FGR),"(",
                                        sprintf("%.04f",Wolbers_C_index_lower),"-",
                                        sprintf("%.04f",Wolbers_C_index_upper),")"),
              "AUC_t=5"=paste0(sprintf("%.04f",AUC_and_Brier_info$AUC$score[1]$AUC)," (",
                               sprintf("%.04f",AUC_and_Brier_info$AUC$score[1]$lower),"-",
                               sprintf("%.04f",AUC_and_Brier_info$AUC$score[1]$upper),")"),
              "AUC_t=10"=paste0(sprintf("%.04f",AUC_and_Brier_info$AUC$score[2]$AUC)," (",
                               sprintf("%.04f",AUC_and_Brier_info$AUC$score[2]$lower),"-",
                               sprintf("%.04f",AUC_and_Brier_info$AUC$score[2]$upper),")"),
              "Brier_t=5"=paste0(sprintf("%.04f",AUC_and_Brier_info$Brier$score[1]$Brier)," (",
                               sprintf("%.04f",AUC_and_Brier_info$Brier$score[1]$lower),"-",
                               sprintf("%.04f",AUC_and_Brier_info$Brier$score[1]$upper),")"),
              "Brier_t=10"=paste0(sprintf("%.04f",AUC_and_Brier_info$Brier$score[2]$Brier)," (",
                               sprintf("%.04f",AUC_and_Brier_info$Brier$score[2]$lower),"-",
                               sprintf("%.04f",AUC_and_Brier_info$Brier$score[2]$upper),")")
  ) }) %>% 
    do.call(rbind,.)
}))

[Compare two model]

1. Compare Harrell_C_index

Harrell’s C-index는 각 모델 간의 concordance 함수를 통해 비교하며, 차이의 유의성을 Z 검정으로 평가한다

table_2_Harrell_C_index <- lapply(1:length(model_var_list), function(i) {
  fgr_data_for_Harrell_C_index <- survival::finegray(as.formula(paste0("Surv(time_Y, status_competing) ~ ",paste(model_var_list[[i]],collapse="+"))), data=melanoma_dt, etype=1)
  fgr_model_for_Harrell_C_index <- coxph(as.formula(paste0("Surv(fgstart, fgstop, fgstatus) ~ ",paste(model_var_list[[i]],collapse="+"))), data=fgr_data_for_Harrell_C_index)
 
  lapply(1:length(model_var_list), function(j) {
    fgr_data_for_Harrell_C_index_new <- finegray(as.formula(paste0("Surv(time_Y, status_competing) ~ ",paste(model_var_list[[j]],collapse="+"))), data=melanoma_dt, etype=1)
    fgr_model_for_Harrell_C_index_new <- coxph(as.formula(paste0("Surv(fgstart, fgstop, fgstatus) ~ ",paste(model_var_list[[j]],collapse="+"))), data=fgr_data_for_Harrell_C_index_new)
  
    ctest <- concordance(fgr_model_for_Harrell_C_index, fgr_model_for_Harrell_C_index_new)
    contr <- c(-1, 1)
    dtest <- contr %*% coef(ctest)
    dvar <- contr %*% vcov(ctest) %*% contr
    z <- dtest/sqrt(dvar)
    p_value_temp <- 2 * (1 - pnorm(abs(z))) 
    p_value <- ifelse(p_value_temp<0.001,"<0.001",sprintf("%.03f",p_value_temp))
    return(p_value)
  }) %>% unlist()
}) %>% do.call(rbind,.) %>% as.data.table() %>% setnames(c("A","B","C","D")) %>% cbind(Harrell_C_index=c("A","B","C","D"),.) %>% replace(is.na(.),"-")
rmarkdown::paged_table(table_2_Harrell_C_index)

2. Compare Wolbers_C_index

Wolbers’ C-index는 bootstrap을 이용한 paired-sample t-test로 두 모델 간 차이를 평가한다. bootstrapping이므로 pair(model1, model2)에 대해 1번만 코드 실행한다.

set.seed(10)
invisible(capture.output({
   table_2_Wolbers_C_index <- lapply(1:(length(model_var_list)-1), function(i) {
  cindex_values <- numeric()
  for(b in 1:20){
    indices <- sample(seq_len(nrow(melanoma_dt)), size=nrow(melanoma_dt), replace=T)
    boot_data <- melanoma_dt[indices]
    fgr_model_for_Wolbers_C_index_boot <- riskRegression::FGR(Hist_formula_list[[i]], data=boot_data, cause=1)
    Wolbers_C_index_info_boot <- pec::cindex(object=fgr_model_for_Wolbers_C_index_boot, formula=Hist_formula_list[[i]], data=boot_data, cause=1, verbose=F)
    cindex_values[[b]] <- Wolbers_C_index_info_boot$AppCindex$FGR
  }

  lapply ((i+1):length(model_var_list), function(j) { # A와 B,C,D 비교, B와 C,D 비교, C와 D 비교 
    cindex_values_new <- numeric()
    for(b in 1:20){
      indices_new <- sample(seq_len(nrow(melanoma_dt)), size=nrow(melanoma_dt), replace=T)
      boot_data_new <- melanoma_dt[indices_new]
      fgr_model_for_Wolbers_C_index_boot_new <- riskRegression::FGR(Hist_formula_list[[j]], data=boot_data_new, cause=1)
      Wolbers_C_index_info_boot_new <- pec::cindex(object=fgr_model_for_Wolbers_C_index_boot_new, formula=Hist_formula_list[[j]], data=boot_data_new, cause=1, verbose=F)
      cindex_values_new[[b]] <- Wolbers_C_index_info_boot_new$AppCindex$FGR
    }
      
    p_value_temp <- t.test(cindex_values-cindex_values_new, mu=0)$p.value
    return( ifelse(p_value_temp<0.001,"<0.001",sprintf("%.03f",p_value_temp)) )
  }) %>% unlist() %>% c(rep("-",(i-1)),.)
}) %>% do.call(rbind,.) %>% as.data.table() %>% setnames(c("B","C","D")) %>% 
  cbind("A"=rep("-",3),.) %>% rbind("D"=data.table("-","-","-","-"), use.names=F) %>% cbind(Wolbers_C_index=c("A","B","C","D"),.) %>% replace(is.na(.),"-")
}))
rmarkdown::paged_table(table_2_Wolbers_C_index)

3. Compare AUC_and_Brier

AUC와 Brier Score는 riskRegression::Score 함수를 이용하여 각 모델 간의 성능 차이를 비교한다.

invisible(capture.output({
  table_2_AUC_and_Brier <- lapply(1:length(model_var_list), function(i) {
  fgr_model_AUC_and_Brier <- riskRegression::FGR(Hist_formula_list[[i]], data=melanoma_dt, cause=1)
  
  lapply (1:length(model_var_list), function(j) {
    fgr_model_AUC_and_Brier_new <- riskRegression::FGR(Hist_formula_list[[j]], data=melanoma_dt, cause=1)
    
    AUC_and_Brier <- riskRegression::Score(list(fgr_model_AUC_and_Brier,fgr_model_AUC_and_Brier_new),
                                           formula=Hist(time_Y, status_competing) ~ 1, data=melanoma_dt, times=c(5,10), null.model=F)
    
    data.table("V1"=c(sprintf("%.03f",AUC_and_Brier$AUC$contrasts$p), sprintf("%.03f",AUC_and_Brier$Brier$contrasts$p))) %>% setnames(as.character(j))
  }) %>% do.call(cbind,.)
  })
}))

table_2_AUC_t_5 <- lapply(1:length(model_var_list), function(k){
  table_2_AUC_and_Brier[[k]][1]
}) %>% do.call(rbind,.) %>% setnames(c("A","B","C","D")) %>% cbind(AUC_t_5=c("A","B","C","D"),.)
for (i in seq(length(model_var_list))) {
  table_2_AUC_t_5[i,(i+1):="-"]
}
rmarkdown::paged_table(table_2_AUC_t_5)

table_2_AUC_t_10 <- lapply(1:length(model_var_list), function(k){
  table_2_AUC_and_Brier[[k]][2]
}) %>% do.call(rbind,.) %>% setnames(c("A","B","C","D")) %>% cbind(AUC_t_10=c("A","B","C","D"),.)
for (i in seq(length(model_var_list))) {
  table_2_AUC_t_10[i,(i+1):="-"]
}
rmarkdown::paged_table(table_2_AUC_t_10)

table_2_Brier_t_5 <- lapply(1:length(model_var_list), function(k){
  table_2_AUC_and_Brier[[k]][3]
}) %>% do.call(rbind,.) %>% setnames(c("A","B","C","D")) %>% cbind(Brier_t_5=c("A","B","C","D"),.)
for (i in seq(length(model_var_list))) {
  table_2_Brier_t_5[i,(i+1):="-"]
}
rmarkdown::paged_table(table_2_Brier_t_5)

table_2_Brier_t_10 <- lapply(1:length(model_var_list), function(k){
  table_2_AUC_and_Brier[[k]][4]
}) %>% do.call(rbind,.) %>% setnames(c("A","B","C","D")) %>% cbind(Brier_t_10=c("A","B","C","D"),.) %>% replace(is.na(.),"-")
for (i in seq(length(model_var_list))) {
  table_2_Brier_t_10[i,(i+1):="-"]
}
rmarkdown::paged_table(table_2_Brier_t_10)

각 지표는 성능이 좋을수록 C-index와 AUC는 높게, Brier Score는 낮게 나타난다.
모델 간 비교 시 p-value가 작을수록 두 모델 간 성능의 유의한 차이가 있음을 의미한다.p-value가 0.05 미만이면 통계적으로 유의하게 두 모델의 성능이 다름을 나타낸다.
이러한 과정을 통해 Competing Risk 모델의 성능을 객관적으로 평가하고 비교할 수 있으며, 실무적 활용도가 높은 모델을 선택하는 데 중요한 기준을 제공한다.