tabwhiteก่อนเริ่มจะขอสรุปก่อนนะครับ เป้าหมายของบทความนี้คือการอธิบายที่มาที่ไปของ survival analysis และวิธีการสร้าง Kaplan-Meier curveรวมถึงความหมายของhazard, hazard ratio, right censoring ก็เลยจัดเป็นหมวด stat & R
tabwhiteไว้บทความต่อไปค่อยพูดถึงการ appraise paper ที่มีการใช้ survival analysis, การจับผิดกราฟ, การแปลผล HR vs Median survival time และ selection bias จากสิ่งที่เรียกว่า Left truncation อีกที


เวลานั้นสำคัญไฉน? [Time-to-Event]
tabwhiteวันนี้เราจะมาพูดกันถึงรูปแบบการวิเคราะห์ข้อมูลอีกประเภทหนึ่งที่มีเรื่องของเวลามาเกี่ยวข้อง เพื่อให้เห็นความสำคัญของ“เวลา”เราลองคิดตามตัวอย่างตลกๆนะครับ
tabwhiteสมมติผมต้องการศึกษาว่าคนอายุซัก40ปี ที่มีLDLสูงๆเนี่ยกินstatinเทียบกับplaceboแล้วจะตายลดลงหรือไม่? แต่ปรากฎว่าผมต้องส่งผลงานวิจัยในอีก2เดือน ก็เลยให้ยาแล้วติดตามไปแค่2เดือน ผลก็คือ…2เดือนยังไม่ผ่านIRBครับ…ล้อเล่นครับ(แต่จริงT-T) ผลก็คือทั้งสองกลุ่มมีอัตราตายที่2เดือนไม่ต่างกัน … คือไม่มีใครตายเลย!
tabwhiteในทางกลับกันครับสมมติว่าอาจารย์ท่านหนึ่งทำการศึกษานี้ค้างไว้แล้วผมมาสานต่อ ติดตามไปรวมกัน60ปีนับจากวันแรกที่กินยา ผลก็คือทั้งสองกลุ่มอัตราตายที่60ปีหลังกินยาไม่ต่างกัน … คือไม่มีใครเหลือให้ศึกษาแล้ว!
tabwhiteนี่แปลว่าstatinไม่มีประโยชน์หรือเปล่า? คงไม่ใช่ จะเห็นได้ว่าคนเราทุกคนเกิดมาก็มีmortality rate=100%ทั้งนั้นแหละ ดังนั้นถ้าอยากจะเปรียบเทียบว่าคนที่กินstatinมีผลลดอัตราตายได้ไหม วิธีหนึ่งก็คือการกำหนดว่าเราจะดูเป็น 5-year mortality rate, 10-y mortality rate ฯลฯ
//ผมอยากเน้นว่า”อัตราการเกิดevent”หรือ”Risk”นั้น จะต้องมีเวลากำกับเสมอ โดยทั่วไปเราจะละไว้ในฐานที่เข้าใจว่าเวลานั้นก็คือช่วงเวลาที่ทำการศึกษานั่นแหละ สั้นยาวแล้วแต่eventที่สนใจ — ดังนั้นRisk(มีอีกชื่อที่ชัดเจนกว่าว่าCumulative incidence)จึงเป็นค่าที่ขึ้นกับระยะเวลาในการศึกษา

tabwhiteอีกวิธีหนึ่งคือการดูว่ากลุ่มไหนตายก่อนกัน ถ้ากินstatinแล้วตายช้ากว่า-time to event นานกว่า-ก็น่าจะแปลว่าดีกว่ากลุ่มที่ไม่กินstatinใช่ไหมครับ? วิธีการวิเคราะห์ตัวแปรตามที่เป็น”time to event”เช่นนี้ เราเรียกว่า”Survival analysis”ตามที่มาว่าตอนแรกใช้ศึกษาเรื่องการรอดชีวิต แต่จริงๆแล้ว การวิเคราะห์แบบนี้ยังใช้ได้ในกรณีอื่นๆ เช่น ถ้าเป็นมะเร็งแล้วเราให้ยาไปอยากรู้ว่าจะลดการกลับเป็นซ้ำได้ไหมก็อาจจะดูว่ายาสูตรไหนมีแนวโน้มจะกลับเป็นซ้ำเร็วช้ากว่ากัน(time to recurrence) ถ้าเราจะคิดค้นยาแก้หวัดเราก็คงอยากจะวัดผลว่ายานี้ทำให้หายหวัดเร็วขึ้นกี่วัน(time to cure-เพราะทุกคนก็ต้องหายอยู่แล้ว) หรือจะเป็นระยะเวลาที่คนไข้aphasiaกลับมาพูดก็ได้ ระยะเวลาที่คนไข้dementiaต้องเข้าnursing homeก็ได้ ฯลฯ


แล้วSurvival analysisมันมีอะไรพิเศษ? [Censoring]
tabwhiteผมคิดว่าถ้าคนอ่านไม่เคยรู้จักSurvival analysisมาก่อน(ซึ่งคงมีน้อย)แล้วผมบอกว่า ไหนลองเอาตัวแปรtime to eventมาวิเคราะห์ซิ สิ่งที่จะเกิดขึ้นก็น่าจะเป็นการแสดงผลด้วยmeanหรือmedianใช่ไหมครับ ถ้าเป็นการเปรียบเทียบก็ใช้ t test หรือ Mann-Whitney U Test ซึ่งผมคิดว่าไม่ผิดเลยนะครับ แต่ทีนี้มันมีแง่มุมพิเศษอย่างหนึ่งก็คือ ในการศึกษาเมื่อเราตรวจติดตามผู้ป่วยไปนานๆ บางคนก็ไม่มีeventเกิดขึ้นซะที หรือถ้าติดตามเป็นสิบๆปีบางคนก็หายไปจากการศึกษา หรือแม้กระทั่งตายไปก่อนจะเกิดeventขึ้น (เช่น โดนรถชนตายตอนที่มะเร็งยังไม่กลับมา เศร้า T-T) … ลองคิดตามนะครับว่า ถ้าเราติดตามคนคนนึงมา5ปีแล้วเค้ายังไม่ตาย แต่ดันย้ายบ้านหนีไป แบบนี้เราจะทำอะไรกับข้อมูลของคนคนนี้ดี? ตัดออกจากการศึกษาเลยแล้วสูญเสียความจริงที่ว่า”อย่างน้อย5ปีนี้ไม่เกิดeventขึ้น”ทิ้งไป หรือจะสมมติว่าตายที่5ปีพอดี? คำตอบก็คือการลงข้อมูลว่าเขาคนนี้”เกิดeventที่เวลา>5ปี”หรือนิยมเขียนว่า”5+” เราเรียกการทำแบบนี้กับข้อมูลว่า“Right Censoring” (แน่นอนว่าเมื่อมีrightก็ต้องมีleftนะครับ และยังมีintervalอีกคำนึง แต่ในงานวิจัยทางคลินิกเราไม่ค่อยได้ใช้ เลยขอละไว้) ซึ่งการนำค่าที่censorไปใช้ก็คือ “ถือว่าก่อนเวลา5ปีคนนี้ไม่มีevent” และ “หลังเวลา5ปีไม่มีคนนี้อยู่ในการศึกษา” ฟังดูป่วงๆใช่ไหมครับ เรามาดูตอนต่อไปเลยดีกว่าว่ามันแปลว่าอย่างไร
//Right censoring เกิดได้จากสาเหตุใหญ่ๆ2แบบก็คือ 1.)ผู้วิจัยจงใจยุติการติดตาม เช่น ตั้งเป้าไว้ว่าจะติดตามไป 10ปี, ติดตามถึงวันที่ 31ธค65, วิจัยจนกว่าจะเกิดeventครบ100ครั้ง 2.)ไม่สามารถได้ข้อมูลeventของsubjectนั้นแล้ว เช่น loss-follow-up, ตายก่อน(ในกรณีที่eventที่สนใจไม่ใช่การตาย)


การแสดงผลSurvival curve [Kaplan-Meier Plot]
tabwhiteถ้าจะพูดให้ยาวก็ต้องเริ่มเกริ่นตั้งแต่การทำlife-table, Product-limit estimate ซึ่งเพื่อไม่ให้ยาวเกินไป(ยาวมากแล้ว) จะขอแสดงวิธีการทำ Kaplan-Meier plot เลยละกันนะครับ
tabwhiteใช้ข้อมูลจาก rats ซึ่งเป็นข้อมูลตัวอย่างใน package survival ของ R ซึ่งเป็นข้อมูลหนูที่โดนจับกรอกยาจนเกิดtumorขึ้นมา สำหรับการสาธิตขั้นแรกผมจะสมมติให้ไม่มีการcensorก่อนนะครับ ทางซ้ายสุดคือข้อมูลดิบซึ่งขั้นต่ำควรจะมีตัวแปรtime-to-eventและตัวแปรว่าcensorหรือไม่(1=event,0=censored)
surv1.png
tabwhite ต่อไปดูตารางที่อยู่เหนือกราฟ คอลัมน์ hazard นั้นคือโอกาสที่หนูที่อยู่มาถึงเวลา[t]จะเกิดtumorขึ้น=event/N at risk อาจกล่าวได้ว่าriskปกติที่พวกเรารู้จักกันคือ”cumulative incidence” ส่วนhazardคือ”instantaneous incidence”นะครับ ส่วนsurvivalก็คือโอกาสที่หนูที่อยู่มาถึงเวลานี้จะอยู่รอดต่อไปซึ่งก็คือ1-hazardนั่นแหละ ค่าsurvivalนี้นำไปคูณกับCumulative survivalในเวลาก่อนหน้านี้ = cumulative survivalอันใหม่ (เช่น ช่วงที่2 มีโอกาสรอดชีวิต 0.979 แต่ก่อนจะมาถึงช่วงที่2ได้ต้องรอดชีวิตจากช่วงที่1มาก่อน ซึ่งมีโอกาส 0.98 นำมาคูณกันก็ได้เป็น 0.96 พอดี) และก็คำนวณอย่างนี้เรื่อยมา เสร็จแล้วนำcumulative survivalไปplotเทียบกับเวลา ก็จะได้ Kaplan-Meier curveออกมา
tabwhiteทีนี้ลองเป็นข้อมูลที่มีการcensorข้อมูลบ้าง(ในที่นี้คือหนูอาจจะตายไปเองก่อนเกิดเนื้องอกหรือหลุดจากกรงก็ไม่ทราบได้)
surv2.png
tabwhitetabwhiteดูตารางเหนือกราฟ เมื่อeventบางส่วนกลายเป็นการcensor hazardกับsurvivalนั้นคำนวณแบบเดิมโดยไม่นำcensorมาคิดในเวลานั้น แต่พอถึงช่วงเวลาถัดไปจะเห็นว่าได้หักคนที่censorแล้วออกไปจาก N at risk ทีนี้อยากให้สังเกตว่าในภาพแรก เวลามีevent1ครั้งcumulative survivalจะลดลง2%เสมอ(1ใน50ตัว) แต่ในภาพที่สองจะพบว่า หลังจากเริ่มมีการcensorแล้ว1eventดึงcumulative survivalลงมามากกว่า2%! กลับขึ้นไปดูวิธีการนำcensorไปใช้ข้างบน ผมบอกไว้ว่า

การนำค่าที่censorไปใช้ก็คือ “ถือว่าก่อนเวลา5ปีคนนี้ไม่มีevent” และ “หลังเวลา5ปีไม่มีคนนี้อยู่ในการศึกษา”

tabwhiteใช่แล้วครับ นั่นก็คือตัวหารน้อยลงนั่นเอง ผลของแต่ละคนที่เกิดeventจะมีน้ำหนักมากขึ้นหลังจากมีการcensor ถ้าดูกราฟจะเห็นว่ามีขีดตั้งๆขึ้นมา ขีดพวกนี้เป็นตัวแทนของการcensorครับ ซึ่งถ้าข้อมูลเยอะมากๆขีดจะยุ่บยั่บไปหมดทำให้บางทีหลายงานวิจัยไม่แสดงขีดนี้ให้เห็น

Screenshot from 2017-03-06 21:55:38.png
tabwhiteอีกตัวอย่างนึงข้อมูลการตายในรพ.ของคนไข้heart transplant อันนี้จะให้เห็นว่าสามารถคำนวณและplot95%CIของcumulative survivalได้ด้วย (ซึ่งมีวิธีแสดงผลหลายแบบจะไม่กล่าวถึงในที่นี้) จากกราฟนี้เราสามารถดูด้วยตาและบอกข้อมูลได้หลายอย่าง ดังนี้ครับ
สีเขียว) แม้กราฟนี้จะไม่แสดงการcensor แต่เราก็เดาได้ว่ามีเกิดขึ้นแน่ๆ เพราะกราฟตกเยอะในช่วงท้ายๆ ถ้าไม่ใช่เพราะcensorอีกสาเหตุหนึ่งก็คือในเวลานั้นมีคนตายพร้อมกันหลายคน! รพ.ไฟดับเหรอ ก็ไม่น่านะครับ
สีน้ำเงิน) เราสามารถบอก median survival time ได้จากการการเล็งตำแหน่งที่cumulative survivalเหลือ0.5ว่าเป็นเวลาเท่าไร
สีแดง) ถ้าเราอยากรู้ 2,000-day survival rate ก็บอกได้ด้วยวิธีคล้ายๆกัน


เทียบ Survival curve ทำอย่างไร? [Log-rank test & Hazard ratio]
tabwhiteถ้ามีข้อมูลsurvival≥2กลุ่มที่เราอยากจะเปรียบเทียบกัน สถิติที่ใช้ในงานวิจัยส่วนใหญ่ก็คือLog-rank testครับ ซึ่งในรายละเอียดก็คือการทำchi-sqaure testในทุกๆจุดเวลาที่มีevent เกิดขึ้นแล้วนำผลมารวมกัน ผลลัพธ์จะออกมาเป็นp-valueซึ่งบอกแค่ว่าsurvival curve 2 เส้นต่างกันหรือไม่ แต่ไม่ได้บอกขนาดและทิศทางของความต่าง ซึ่งต้องใช้วิธีอื่นเพื่อคำนวณ hazard ratio ซึ่งมีวิธีการแปลผลคล้ายๆกับ RR ครับ แต่วิธีที่ใช้คำนวณนั้นมีตัวเลือกและมีความซับซ้อนอยู่เหมือนกันจะไม่ขอลงลึก ทีนี้ลองมาดูข้อมูลrecurrence free survivalในคนไข้AMLยุคโบราณ(ต้องออกตัวก่อนเดี๋ยวโดนHematoเขม่น สมัยนี้การรักษาดีมากแล้วครับ)

Rplot.png

tabwhiteระหว่างกลุ่มที่ได้กับไม่ได้ยาถามว่าการเกิดrecurrenceต่างกันหรือไม่? ถ้ามองด้วยตาก็จะเห็นว่าเส้นแดงดูจะดีกว่าเส้นน้ำเงินนิดหน่อยแต่ก็ไม่ทิ้งกันมาก แต่ว่ามันsignificantหรือเปล่า? ถ้าทำLog-rank testจะได้ p 0.147 แปลว่าไม่พบความต่างกันอย่างมีนัยสำคัญ ถ้าใช้Cox Propotional Hazard Regression Model จะได้ HR ของกลุ่มที่ไม่ได้ยาเทียบกับกลุ่มได้ยา = 2.07(95%CI 0.76-5.62) ก็แปลว่ามีแนวโน้มว่าถ้าไม่ได้ยาจะกลับเป็นซ้ำเร็วกว่าแต่ก็ไม่sigเช่นกันเพราะคร่อม 1

สรุปสิ่งที่อยากให้เข้าใจคือ
1. การวิเคราะห์ที่ตัวแปรตามเป็น time-to-event และมีการ censoring มีแนวคิดอย่างไร
2. Hazard vs Risk, Hazard ratio vs Relative risk
3. ที่มาของ Kaplan-Meier curve และองค์ประกอบในกราฟ

ตอนหน้าค่อยมาต่อกันที่การappriase paperที่มีการใช้ survival analysis นะครับ ขอย้ำว่า survivial analysisเป็นวิธีวิเคราะห์ข้อมูลเท่านั้น สามารถปรากฏอยู่ใน cohort, case-control, RCT ได้หมด เป็นการดู mortality, harm, benefit ได้หมด แค่ชื่อ Survival analysis ไปงั้นเอง


#R code for survival analysis
library("survival")
#extract data from dataset "rats" for simple demonstration
#select only female rat because males have very few tumor
data1 <- rats[rats$sex=="f"&rats$rx==1,c("time","status")]
#status 0 = censored, 1 = event occured
head(data1)
#sort by time
data1 <- data1[order(data1$time),]
#Surv function add censor mark (+) to time data
data1$surv <- Surv(data1$time,data1$status)
#Add FALSE survival data, assume no censor
data1$Xsurv <- Surv(data1$time,data1$status>-1)
#Prepare data for survival plot
Splot1 <- survfit(surv~1,data=data1)
SplotX <- survfit(Xsurv~1,data=data1)

#create new function to demonstrate life table
LifeTable <- function(S){
Life <- data.frame(
time=S$time,
N=S$n.risk,
event=S$n.event,
censor=S$n.censor
)
Life$hazard <- Life$event/Life$N
Life$surv <- 1-Life$hazard
Life$cumSurv<-S$surv
Life$nextN <- Life$N-Life$event-Life$censor
return(round(Life,4))
}
#display "life table" of no-censored data
LifeTable(SplotX)
plot(SplotX,conf.int = FALSE,xlab="time",ylab="cumulative survival",main="DATA WITH NO CENSOR")
#display "life table" of censored data
LifeTable(Splot1)
plot(Splot1,conf.int = FALSE,xlab="time",ylab="cumulative survival",main="DATA WITH CENSOR",mark.time = TRUE)

#quick code for plot KM of heart transplant patients
plot(survfit(Surv(stanford2$time,stanford2$status)~1))
#now for comparison
data2 <- rats[rats$sex=="f",]
#survival analysis, split data by $rx value
Splot2 <- survfit(Surv(data2$time,data2$status)~data2$rx)
summary(Splot2)
plot(Splot2,conf.int = FALSE,xlab="time",ylab="cumulative survival",main="DATA WITH COMPARE")
#Log-rank & CoxPH
LRnk <- survdiff(Surv(data2$time,data2$status)~data2$rx)
summary(coxph(Surv(data2$time,data2$status)~data2$rx))

#for AML data
data3 <- aml[aml$time!=161,]
data3
plot(survfit(Surv(data3$time,data3$status)~data3$x),col=c("red","blue"),
main="AML recurrence",xlab="Months",ylab="Recurrence free survival")
survdiff(Surv(data3$time,data3$status)~data3$x)
summary(coxph(Surv(data3$time,data3$status)~data3$x))
text(x=35,y=0.8,labels="Median survival = 31 months",col="red")
text(x=35,y=0.85,labels="With chemotherapy",col="red")
text(x=15,y=0.2,labels="Median survival = 23 months",col="blue")
text(x=15,y=0.25,labels="Without chemotherapy",col="blue")

Advertisements