Survival analysis - Log Rank, Coxph

 

/*

 * http://sosal.kr/
 * made by so_Sal
 */

 

 

 

1. 파일의 구조 및 데이터 정보

2. 생존과 관련된 환자정보

3. LogRank test로 위의 두 그룹이 유의하게 생존에서 차이를 보이는지 확인하기

4. Cox-proprtional harzard regression analysis: 회귀를 이용한 생존분석

 

 

 

 

 

Log rank test

독립변수를 통해 group화 된 샘플이, 그룹간에 생존분포 차이가 있는지 확인할 수 있는 가설검정 (hypothesis test)

유의한 결과가 나온다면 즉 해당 독립변수는 prognostic factor로써 좋은 기능을 한다는 것.

 

 

Cox-proprtional harzard regression analysis

콕스 비례위험모형, 줄여서 coxph라고 한다. 회귀 기법으로 survival data를 분석할 수 있기 때문에, n개의 그룹으로 나누지 못하는 독립변수에 대해서도 prognostic factor로써 좋은 기능을 할 수 있는지 평가할 수 있다.

 

 

 

 

이 글에서는 자세한 통계적 계산 원리나 개념들은 다루지 않고, 오로지 R 프로그래밍의 survival 패키지 안에 구현되어 있는 survdiff, coxph 등의 함수를 사용해보고, 어떠한 정보들을 얻을 수 있으며 만들어진 모델로부터 어떤 정보를 가져올 수 있는지 등, 유저 측면에서 간단하게 설명하고자 한다.

 

 

 

coxph 모델을 돌리기 위한 데이터는 TCGA (The Cancer Genome Atlas) 에서 폐선암 - LUAD (lung adenocarcinoma) 환자의 clinical data를 다운받아서 가공한 것이다.

 

 

---

Data download:patients.txt

---

 

 

 

 

 

  

1. 파일의 구조 및 데이터 정보

 

 

> clinical <- read.table("patients.txt", sep="\t", header=TRUE, stringsAsFactors=FALSE)

> head(clinical) 

 

dimension 함수를 통해 읽어들인 데이터의 row와 column의 수가 어떻게 되는지 쉽게 알 수 있다.

> dim(clinical)
[1] 476   7

 

총 476명의 환자와 7종류의 데이터 column이 존재한다는 것을 알 수 있다.

 

 

column은 총 7개로 구성되어 있으며, Barcode는 환자를 인식하는 번호정도로 생각하면 된다.

Age는 day이며, Gender, Status, Stage, tobacco는 factor 데이터로 사용하면 더욱 편하기 때문에, factor로 변환해준다.

 

factor로 변환된 데이터는 summary 함수를 통해 각 level에 몇명의 환자(혹은 데이터)가 존재하는 지 쉽게 알 수 있다.

 

 

> clinical$Gender <- factor(clinical$Gender)
> summary(clinical$Gender)
FEMALE   MALE
   257    219

 

> clinical$Status <- factor(clinical$Status)

> summary(clinical$Status)
Alive  Dead
  358   118

 

> summary(clinical$Stage)
Stage I  Stage II Stage III  Stage IV      NA's 
     255        114          81          25           1

 

> clinical$Tobacco <- factor(clinical$Tobacco)

 

 

 

str (structure) 함수를 사용하면 clinical의 데이터의 구조를 쉽게 알 수 있다.

> str(clinical)

 

어떤 변수들이 존재하는지, 각 변수들의 데이터 타입은 어떻게 되는지 알 수 있다.

친절하게 Age, Time 등은 자동으로 int로 받아들여 졌으며, Barcode는 문자열로, 나머지는 위에서 factor로 형변환 되었기 때문에 factor로 표시되어졌다.

 

clinical 안에 존재하는 다양한 변수들에 접근하려면 clinical$variable 형식으로 접근해야 하지만, attach() 함수를 통해 frame name을 쓰지 않더라도 접근할 수 있다.

 

> attach(clinical) # clinical$Age를 Age라는 이름으로 접근이 가능하다.

 

 

 

 

 

2. 생존과 관련된 환자정보

 

2.1 본격적으로 생존분석에 들어가기에 앞서 survival이라는 library를 불러와야한다.

 

> library(survival)

survival이라는 library에는 여기서 다루고자 하는 coxph함수와 더불어, Surv, Survfit, Survreg 등 생존분석에 활용되는 다양한 함수 및 회귀기법 함수들이 내장되어 있다. 추가적인 내용 및 예제, 함수 사용법들은 아래 cran에서 제공하는 메뉴얼을 참조하면 된다.

https://cran.r-project.org/web/packages/survival/survival.pdf

 

 

2.2 Surv() - Survival object 생성 함수. 생존시간 및 생존상태를 입력으로 받아 object를 생성한다.

 

> Surv(Time, Status=="Dead")

 

두번째인자는 True (혹은 1)일 경우 죽은 환자로 인식하기 때문에, Status=="Alive"가 아닌 Status=="Dead"로 해주어야 한다.

Alive 상태인 환자들은 +라고 표시되며, 이는 right-censored data로 인식한다.

 

 

 

 

2.3 kaplan-meier 카플란 메이어 graph 그리기

 

성별을 통해 환자들을 두 group으로 나누고 kaplan-meier를 통해 두 그룹이 차이가 있는지 생존상태를 확인해 볼 수 있다.

> fit = survfit(Surv(Time, Status=="Dead")~Gender)

> plot(fit, col=1:2)

> legend("topright", legend=levels(Gender), col=1:2, lty=1)

 

 

 

 

3. LogRank test로 위의 두 그룹이 유의하게 생존에서 차이를 보이는지 확인하기

 

 

> survdiff( Surv(Time, Status=="Dead") ~ Gender )

 

p-value: 0.993

성별은 폐 선암의 예후와 상관없다는 것을 알 수 있다.

 

 

> survdiff(Surv(Time,Status=="Dead")~(Age>median(Age,na.rm = TRUE)))

 

p-value: 0.0167

나이를 median 값을 기준으로 2그룹으로 나누었더니, 굉장히 유의하게 나왔다.

확실히 나이가 많은 그룹과 나이가 적은 그룹이 차이가 난다는 건 알 수 있지만,

상황에 따라서 이렇게 median으로 나누는 것은 좋지 못하다.

 

 

 

> survdiff(Surv(Time,Status=="Dead")~Stage)

 

 

p-value: 0.00000000934

암의 Stage로 나눈 환자들의 생존분포가 매우 유의하다는 것을 알 수 있다.

 

 

Kaplan-meier 그래프를 이용하여 그룹별로 생존분포가 어떤지 한눈에 확인할 수 있다. 

> fit = survfit(Surv(Time, Status=="Dead")~Stage)
> plot(fit, col=1:4)
> legend("topright", legend=levels(Stage), col=1:4, lty=1)

 

 

 

 

 

 

4. Cox-proprtional harzard regression analysis: 회귀를 이용한 생존분석

 

Log-rank test에서 유의하게 나온 Stage를 coxph함수를 통해 모델링 해보면 다음과 같다.

 

> coxph(formula = Surv(Time, Status == "Dead") ~ Stage)

 

Stage II, III, IV 독립변수 모두 p-value 0.05 이하로 유의하게 나왔다.

모델을 평가해주는 LRT (Likelihood ratio test) 역시 p-value가 굉장히 유의하게 나왔다.

Stage라는 독립변수는 I, II, III, IV로 총 4개의 level로 구성되어 있다.

factor변수는 회귀모델에서 각각의 level이 생존을 예측하는 독립변수로써 얼마나 큰 역할을 했는지를 서로간에 평가하게 된다.

일반적으로 Factor의 순서상 가장 앞에 있는 level이 reference가 되어, 나머지 level과 비교하게 되는데, 이 경우에는 Stage I이 reference가 되었다. 따라서 Stage I (Reference)는 Stage II, III, IV와 비교하여 얼마나 생존이 다른지 평가되어졌고, 그 결과 p-value가 모두 유의하게 나왔다.

 

또한 Stage라는 독립변수만으로 모델을 fitting한 결과, 모델의 p-value가 굉장히 작게 나온것으로 보아 Stage가 생존을 예측하는데 매우 좋은 prognostic factor라는 것을 알 수 있다.

 

 

 

- Coef, exp(Coef)

 

각 독립변수의 coef는 선형회귀처럼 종속변수 (y axis - survival)에 대한 기울기를 뜻한다. 따라서 +positive coef인 경우에 생존에 worse prognosis한 독립변수라고 할 수 있으며, -negative coef인 경우 protective effect를 가진 독립변수라고 해석할 수 있다. 또한 생존분석이기 때문에, 일반적인 선형회귀식에서 y = β0 + β1*x 에서 식 자체가 음수가 불가능 하다 (생존이 -일 수 없다.). 따라서 콕스비례위험 모델에서는 각 독립변수들이 그 자체로 exp() 되어져 지수간의 합으로 존재하기 때문에, 실제 Hazard ratio를 측정하기 위해선 exp(coef) (=HR)값을 사용해야 한다. reference로부터 나눠진 값이므로, exp(coef) 값이 1보다 크다면 reference보다 예후가 더 좋지 않음을 뜻하며, 1보다 작다면 예후가 더 좋은 prognostic factor로써 작용한다고 이해할 수 있다.

 

 

 

 

> coxph(formula = Surv(Time, Status == "Dead") ~ Age)

 

 

Age라는 독립변수 하나만으로 회귀모델을 만들었더니, 독립변수 Age는 p-value 0.1,

Age 하나만으로 생존을 예측하는 모델 자체의 p-value 0.101 이었다.

 

Age 하나만으로 예후를 판단하기에는 쉽지 않다.

 

 

> fit = coxph(formula = Surv(Time, Status == "Dead") ~ Age + Gender + Tobacco + Stage, data=clinical)

 

유의하게 나온 독립변수들을 summary 함수로 쉽게 확인할 수 있다.

 

HRplot을 이용하면 각 독립변수들의 Hazard ratio를 직관적으로 알 수 있다.

 

> library(moonBook)

> HRplot(fit, type=2, show.CI=TRUE, main="Hazard ratios of all individual variables")