14 Survival analysis

ในงานวิจัยที่กระทำกับผู้ป่วย หรือแม้กระทั้งเซลล์นั้น บางครั้งจะมีความจำเป็นที่ต้องทำการวิเคราะห์ข้อมูลเพื่อเปรียบเทียบสร้างแบบจำลองที่สามารถทำนายเวลาที่ใช้ก่อนที่จะเกิดเหตุการณ์ที่ท่านสนใจ (Event) เช่น เวลาที่ผู้ป่วยจะเสียชีวิตจากโรคมะเร็งนับตั้งแต่วันวินิจฉัย เวลาของเซลล์ที่จะตายหลังจากใส่สารบางอย่างที่สนใจ เป็นต้น

ลักษณะพิเศษของการวิเคราะห์ Survival analysis คือ มีการวิเคราะห์โดยใช้ปัจจัยที่เรียกว่า Censoring ร่วม ซึ่งคือการที่ เหตุการณ์ที่คาดหวังว่าจะเกิดนั้นไม่มาถึงแม้ว่าจะครบตามเวลาที่ผู้วิจัยสังเกตการณ์แล้ว ซึ่งทำให้ไม่สามารถมั่นใจได้ว่าเหตุการณ์นั้นจะเกิดต่อไปหรือไม่ ณ เวลาหลังจากนี้

การ Censor โดยหลักมี 3 แบบ คือ

Right-censor (ดังรูป) คือ ไม่แน่ใจข้อมูลการเกิด event เวลาสุดท้ายที่พบ ซึ่งพบมากที่สุด
Left-censor คือ ไม่แน่ใจข้อมูลช่วงเวลาเริ่มต้น เช่น Diagnosis วันไหน
Interval-censor คือ เวลาช่วงใดช่วงหนึ่งหายไป

14.1 Kaplein-Meier estimate (KM)

คือ กราฟแสดงอัตราการเกิดของเหตุการณ์เมื่อเทียบกับเวลาที่ผ่านไป โดยหลักการคำนวณ Survival คือ

\[ \hat{S}(t) = \prod_{i: t_{i} \leq t}(1- \frac{d_{i}}{n_{i}})\]

\[ \text{Survival proability} = \frac{n.risk - n.event}{n.risk} \]

อธิบายหลักการอย่างง่ายของ KM นั่นคือ ทุกเคสทียังไม่เกิดเหตุการณ์นั้น จะเป็น เคสที่เสี่ยงต่อการเกิดเหตุการณ์ (At risk) ซึ่งจำนวนเคสแรกเริ่มที่เสี่ยง (Number at risk) จะเท่ากับจำนวนเคสทั้งหมด (Sample size) โดยจะนับการเกิด Event ตามปกติ เพียงแต่ถ้าเคสนั้นถูก Censor นั้น n.risk จะลดลงด้วย ทำให้อัตราการเกิด Event ยังไม่เปลี่ยนแปลง ดังตัวอย่างตามตาราง

time	n.risk	n.event	n.censor	calculation	estimate
5	10	1	0	9/10	0.9
11	9	3	0	0.9 x 6/9	0.6
12	6	1	0	0.6 x 5/6	0.5
92	5	0	1	0.5 x 5/5	0.5
105	4	0	1	0.5 x 4/4	0.5
173	3	0	1	0.5 x 3/3	0.5
174	2	0	1	0.5 x 2/2	0.5
175	1	0	1	0.5 x 0/1	0

เมื่อนำไปพล็อตกราฟแล้วจะได้ผลดังนี้

สังเกตุว่าส่วนที่ Censor (มีสัญลักษณ์ +) จะไม่มีการตกลงของกราฟ แต่เมื่อถึงเวลาที่มี Event เกิดขึ้น การตกลงของกราฟจะสูงกว่าเมื่อไม่มี Censor นำมาก่อน

14.1.1 การสร้าง KM ใน R

ตัวอย่างข้อมูลของผู้ป่วยมะเร็งรังไข่ที่ได้รับการรักษาโดยการผ่าตัด

library(dplyr) 
library(survival)
names(ovarian)

## [1] "futime"   "fustat"   "age"      "resid.ds" "rx"       "ecog.ps"

อธิบายตัวแปร:

age = อายุ
futime = ระยะเวลาติดตามตั้งแต่วินิจฉัยจนเสียชีวิต/มาพบแพทย์ครั้งสุดท้าย
fustat = 0 - censor, 1 - dead
resid.ds = มีชิ้นส่วนของมะเร็งหลงเหลือหลังจากการผ่าตัด (ผ่าตัดได้ไม่หมด)
rx = กลุ่มการรักษา
ecog.ps = ECOG performance status คะแนนน้อยแปลว่าผู้ป่วยมีสุขภาพโดยรวมดี

เมื่อใช้ Function Surv() จะทำการเปลี่ยน futime ให้รับรู้การ Censor สังเกตว่าผู้ป่วยที่ ไม่เกิดเหตุการณ์จะมีสัญลักษณ์ + อยู่ข้างหลัง บ่งบอกว่าข้อมูลนั้นถูก Censor นั่นหมายความว่า ผู้ป่วยจะเกิดเหตุการณ์หรือไม่ก็ได้หลังจากนี้ เพียงแต่ผู้วิจัยไม่สามารถทราบได้แล้ว

censored_df <- ovarian |>    
  mutate(censored_futime = Surv(ovarian$futime, ovarian$fustat)) 

head(select(censored_df, futime, fustat, censored_futime))

การพล็อต KM นั้นสามารถทำได้โดยใช้ Package survminer โดยเริ่มจากการสร้างตาราง Survival curve จากคำสั่ง survfit()

library(survminer)

ovarian_surv <- survfit(
  Surv(futime/365.25, fustat) ~ 1, data = ovarian # เปลี่ยนเป็นปี
  )

ovarian_surv |> tidy() |> head(10)

หลังจากนั้นใช้คำสั่ง ggsurvplot() เพื่อทำการสร้างกราฟ

ggsurvplot(ovarian_surv, risk.table = TRUE, break.time.by = 0.25, 
           surv.median.line = "hv")

จะเห็นว่าผู้ป่วยกลุ่มนี้มี Median survival อยู่ประมาณ 1.75 ปี

ovarian_surv_resid <- survfit(
  Surv(futime/365.25, fustat) ~ resid.ds, data = ovarian
  )

ggsurvplot(ovarian_surv_resid, risk.table = TRUE, break.time.by = 0.25,
           surv.median.line = "hv", pval = TRUE)

จะพบว่า ถ้าผ่าตัดแล้วไม่เหลือร่องรอยของโรค จะมีอัตราการรอดชีวิตที่ดีกว่า แต่ยังไม่ถึงระดับมีนัยสำคัญ

14.2 Log-rank test

Log-rank test คือ Non-parametric test สำหรับ Univariate analysis ที่เปรียบเทียบความแตกต่างของอัตราการเกิด Event ว่าแตกต่างอย่างมีนัยสำคัญหรือไม่

ovarian_surv_diff <- survdiff(
  Surv(futime/365.25, fustat) ~ resid.ds, data = ovarian
  )

ovarian_surv_diff

## Call:
## survdiff(formula = Surv(futime/365.25, fustat) ~ resid.ds, data = ovarian)
## 
##             N Observed Expected (O-E)^2/E (O-E)^2/V
## resid.ds=1 11        3     6.26      1.70      3.62
## resid.ds=2 15        9     5.74      1.85      3.62
## 
##  Chisq= 3.6  on 1 degrees of freedom, p= 0.06

Observed คือ จำนวน event ที่เกิดขึ้นในแต่ละกลุ่ม
Expected คือ จำนวน event ที่คาดว่าจะเกิดขึ้นในแต่ละกลุ่ม
(O-E)^2/E = Chi-square statistics ของค่า observed และ expected
Chisq = ผลสุดท้ายของ Chi-square statistics = sum(O-E)^2/E
p = p-value ของ Chi-square statistics

ค่า Log-rank นี้ สามารถแสดงใน KM ได้โดยใช้ pval = TRUE ตามหัวข้อเบื้องต้น

14.3 Cox-proportional hazard (CPH) model

คือ Semi-parametric model ซึ่งวัด Risk ของการเกิด Event นั้นๆ โดยมีสมการ คือ

\[ h(t) = h_{0}(t) \times exp(b_{1}X_{1} + b_{2}X_{2}+ ... + b_{s}X_{p}) \]

\[ ln(\frac{h(t)}{h_{0}(t)}) = b_{1}X_{1} + b_{2}X_{2}+ ... + b_{s}X_{p} \]

ซึ่ง \(h(t)/h_{0}(t)\) นั้นคือ Hazard ratio (HR) หรือ ความเสี่ยงของการเกิด event นั้นๆ

การวิเคราะห์ CPH นั้นมีข้อดีกว่า Log-rank คือสามารถประมาณการเชิงปริมาณ (Quantitative measurement) ผ่าน HR และสามารถวิเคราะห์สมการแบบ Multivariate analysis ได้

ovarian_cox <- coxph(Surv(futime/365.25, fustat) ~ 
                       resid.ds + age + factor(rx), data = ovarian) 

ovarian_cox

## Call:
## coxph(formula = Surv(futime/365.25, fustat) ~ resid.ds + age + 
##     factor(rx), data = ovarian)
## 
##                coef exp(coef) se(coef)      z       p
## resid.ds     0.6964    2.0065   0.7585  0.918 0.35858
## age          0.1285    1.1372   0.0473  2.718 0.00657
## factor(rx)2 -0.8489    0.4279   0.6392 -1.328 0.18416
## 
## Likelihood ratio test=16.77  on 3 df, p=0.0007889
## n= 26, number of events= 12

exp(coef) = HR = \(h(t)/h_{0}(t)\) ในที่นี้ท่านสามารถอภิปรายได้ว่า อายุที่เพิ่มขึ้น 1 ปีนั้น ส่งผลให้เกิดอัตราการเสียชีวิตในผู้ป่วยมะเร็งรังไข่เพิ่มขึ้น 1.13 เท่า (13%) และมีนัยสำคัญทางสถิติ
p ในตาราง คือ ค่าคำนวณ \(p\)-value จาก Wald’s test ของแต่ละตัวแปรว่ามีผลต่ออัตราการรอดชีวิตหรือไม่
p ข้างล่าง คือ overall \(p\) จาก Likelihood ratio test ว่าจากทั้งหมด มีตัวแปรใดตัวแปรหนึ่งส่งผลให้อัตรากการรอดชีวิตเปลี่ยนไปอย่างมีนัยสำคัญทางสถิติหรือไม่

14.3.1 การตรวจสอบ Assumption validity ของ CPH

CPH นั้นมี Assumption ดังนี้:

ตัวแปรแต่ละกลุ่มมีอัตราการเกิด Event ที่แตกต่างกัน
HR เท่ากันทุกช่วงเวลา เช่น ที่ 1, 2, 5 ปี อัตราส่วนการเสียชีวิตระหว่างตัวแปรเท่ากันหมด
ตัวแปรมีความสัมพันธ์แบบ Linear continuous variable
ไม่จำเป็นต้องทราบลักษณะการกระจายตัวของข้อมูลก่อน (จึงเป็น Semi-parametric model)

สามารถตรวจสอบ HR ได้โดยใช้ Proportionality assumption test จาก Schoenfeld residuals โดย cox.zph()

ovarian_coxzph <- cox.zph(ovarian_cox) 
ggcoxzph(ovarian_coxzph)

โดยการทดสอบนี้ จะทำการเปรียบเทียบ Residuals ระหว่าง Risk-weight average กับ ตัวแปรนั้นๆ ว่ามีการเปลี่ยนแปลงไปในทิศทางใดทิศทางหนึ่งหรือไม่ ถ้ามี (\(p\) < 0.05) หมายความว่า เวลาที่ผ่านไปอาจจะส่งผลให้ HR นั้นมีความแตกต่างกัน ซึ่งจะต้องทำ Time-varying CPH เพิ่มเติม

โดยในข้อมูล ovarian นี้ ไม่มีตัวใดที่ \(p\)-value < 0.05 จึงถือได้ว่า อัตราส่วนนั้นคงที่ และทำให้ CPH นั้น Valid

ในส่วนของ Linearity สามารถตวจสอบโดยใช้ ggcoxfunctional()

ovarian_linear_age <- ggcoxfunctional(Surv(futime/365.25, fustat) ~ age +                                 + I(log(age)) + I(sqrt(age)), data = ovarian) 
ovarian_linear_age

จะเห็นว่า age นั้นการเพิ่มขึ้นแบบ Linearity โดยมี Deviation เล็กน้อย

13 Regression model

15 Principal component analysis