Kaplan-Meier และ Cox สำหรับเหตุการณ์แรก: KM บรรยาย Cox สร้างแบบจำลอง

Read the English version

บทคัดย่อ

การวิเคราะห์การรอดชีพแบบเหตุการณ์แรก (first-event survival analysis) ถามว่านานเท่าไรจนเกิดเหตุการณ์ครั้งแรก ตัวอย่างเช่นการกลับเป็นซ้ำครั้งแรก การนอนโรงพยาบาลครั้งแรก และการเสียชีวิต ส่วนเหตุการณ์ที่เกิดซ้ำภายหลังไม่นับในการวิเคราะห์นี้ เครื่องมือสองอย่างแบ่งงานกัน Kaplan-Meier ใช้บรรยาย คือประมาณฟังก์ชันการรอดชีพ S(t) ซึ่งเป็นความน่าจะเป็นที่จะยังปลอดเหตุการณ์เกินเวลา t ใช้วาดกราฟการรอดชีพ ประมาณความน่าจะเป็นของการรอดชีพ ประมาณมัธยฐานการรอดชีพเมื่อไปถึง และเปรียบเทียบกลุ่มด้วยตา แบบจำลอง Cox ใช้สร้างแบบจำลอง คือประมาณฮาซาร์ดสัมพัทธ์ (relative hazard) ระหว่างกลุ่ม และปรับตัวแปรร่วม (covariates) ได้ ผลอย่าง HR = 0.70 หมายความว่าฮาซาร์ดที่ประมาณได้ของกลุ่มหนึ่งต่ำกว่า ภายใต้ข้อสมมติของแบบจำลอง ข้อสมมติที่สำคัญที่สุดคือ proportional hazards คือ hazard ratio คงที่โดยประมาณตลอดเวลา จากนั้นมีกับดักสามอย่าง คือ hazard ratio ไม่ใช่ risk ratio, Kaplan-Meier ไม่ได้ปรับอะไรเลย และมัธยฐานอาจยังไม่ถึง บทความนี้สอนหลักเดียว คือ เหตุการณ์แรก Kaplan-Meier บรรยาย Cox สร้างแบบจำลอง

🔍 คลิกเพื่อขยาย คลิกที่ใดก็ได้เพื่อปิด

รายงานการทดลองที่มีทั้งกราฟและอัตราส่วน

แพทย์โรคหัวใจคนหนึ่งอ่านรายงานการทดลองหัวใจล้มเหลว ผลลัพธ์คือเวลาตั้งแต่การสุ่มเข้ากลุ่มจนถึงการเสียชีวิตจากทุกสาเหตุ รูปหนึ่งแสดงกราฟ Kaplan-Meier ของแต่ละกลุ่ม อีกบรรทัดหนึ่งให้ hazard ratio จาก Cox เท่ากับ 0.70 และเพื่อนร่วมงานสรุปว่าเสียชีวิตน้อยลง 30%

การทดลองนี้เป็นของสมมติ ก่อนจะเชื่อตัวเลขใด เธอถามคำถามง่าย ๆ สองคำถาม: การวิเคราะห์นี้ตอบคำถามอะไร และเครื่องมือแต่ละอย่างบอกอะไรเพิ่ม

นาฬิกาหนึ่งเรือนต่อผู้ป่วยหนึ่งคน หยุดที่เหตุการณ์แรก

ผู้ป่วยแต่ละคนให้เวลาหนึ่งค่า: จากจุดเริ่มนับ เช่นการสุ่มเข้ากลุ่ม ไปจนถึงเหตุการณ์แรก ผู้ป่วยที่นอนโรงพยาบาลสามครั้งจะนับเฉพาะเวลาจนนอนโรงพยาบาลครั้งแรก ครั้งที่สองและครั้งที่สามไม่มีบทบาทในการวิเคราะห์นี้ ตอนที่ 6 ของชุดบทความนี้นับการนอนโรงพยาบาลทุกครั้งแทน

ผลลัพธ์แบบเหตุการณ์แรกที่พบบ่อย ได้แก่

• เวลาจนกลับเป็นซ้ำครั้งแรก (time to first recurrence)
• เวลาจนเสียชีวิต (time to death)
• เวลาจนนอนโรงพยาบาลครั้งแรก (time to first admission)
• เวลาจนเกิดกล้ามเนื้อหัวใจตายครั้งแรก (time to first myocardial infarction)
• เวลาจนการรักษาล้มเหลว (time to treatment failure)

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

Kaplan-Meier บรรยาย: ฟังก์ชันการรอดชีพ

ให้ $T$ เป็นเวลาจนเกิดเหตุการณ์แรก Kaplan-Meier ประมาณฟังก์ชันการรอดชีพ (survival function)

อ่านว่า ความน่าจะเป็นที่จะยังปลอดเหตุการณ์เกินเวลา $t$ Kaplan-Meier ประมาณค่านี้จากการติดตามที่สิ้นสุดก่อนเวลาในผู้ป่วยบางคน โดยไม่มีแบบจำลองสำหรับรูปร่างของเส้น [1]

Kaplan-Meier ใช้บรรยายเป็นหลัก (primarily descriptive) มีประโยชน์สำหรับ

• วาดกราฟการรอดชีพ (survival curves)
• ประมาณความน่าจะเป็นของการรอดชีพ (survival probabilities)
• ประมาณมัธยฐานการรอดชีพ (median survival) เมื่อไปถึง
• เปรียบเทียบกลุ่มด้วยตา

การใช้งานทั้งสี่อย่างในผู้ป่วยสิบคน

ผู้ป่วยหัวใจล้มเหลวสิบคนถูกติดตามนับจากการสุ่มนานถึง 24 เดือน และเหตุการณ์คือการเสียชีวิตจากทุกสาเหตุ ห้าคนเสียชีวิตเมื่อเดือนที่ 6, 10, 14, 18 และ 20 หนึ่งคนขาดการติดตามเมื่อเดือนที่ 8 และอีกสี่คนยังมีชีวิตเมื่อการติดตามสิ้นสุดที่เดือนที่ 24

ณ เวลาที่เสียชีวิตแต่ละครั้ง $t_j$ เส้นจะถูกคูณด้วยสัดส่วนของชุดเสี่ยง (at-risk set) ที่รอดผ่านเวลานั้น ชุดเสี่ยงคือทุกคนที่ยังปลอดเหตุการณ์และยังอยู่ในการติดตามก่อนถึง $t_j$ เล็กน้อย:

โดย $n_j$ คือจำนวนในชุดเสี่ยง และ $d_j$ คือจำนวนการเสียชีวิตที่ $t_j$ ผู้ป่วยที่ถูกเซ็นเซอร์ออกจากชุดเสี่ยงไปโดยไม่ทำให้เส้นลดเป็นขั้น

อ่านการใช้งานทั้งสี่อย่างจากกราฟ

วาดกราฟ คอลัมน์การรอดชีพเมื่อวาดเป็นขั้นบันได คือกราฟ Kaplan-Meier: ราบระหว่างการเสียชีวิต และลดลงหนึ่งขั้นที่การเสียชีวิตแต่ละครั้ง

ความน่าจะเป็นของการรอดชีพ ค่าประมาณ ณ เดือนที่กำหนดคือขั้นสุดท้ายที่เดือนนั้นหรือก่อนหน้า: 0.7875 ที่เดือนที่ 12 และ 0.45 ที่เดือนที่ 24

มัธยฐานการรอดชีพ เมื่อไปถึง มัธยฐานการรอดชีพคือเวลาแรกที่เส้นลดถึงครึ่งหนึ่งหรือต่ำกว่า ผู้ป่วยห้าในสิบคนเสียชีวิต เส้นลดจาก 0.5625 เป็น 0.45 ที่เดือนที่ 20 มัธยฐานการรอดชีพจึงเท่ากับ 20 เดือน

เปรียบเทียบกลุ่มด้วยตา เมื่อวาดเส้นหนึ่งเส้นต่อกลุ่มบนแกนเดียวกัน ตาจะเปรียบเทียบกลุ่มได้ ส่วนถัดไปทำเช่นนี้กับสองกลุ่มของการทดลองจำลอง

เปรียบเทียบสองกลุ่ม: การทดลองจำลอง

การอ่านเชิงบรรยายแบบเดียวกันใช้กับข้อมูลที่ใหญ่ขึ้นได้ การทดลองจำลองหนึ่งชุดรับผู้ใหญ่ 1,500 คนหลังเข้ารักษาตัวด้วยหัวใจล้มเหลว กลุ่มละ 750 คน ติดตามตั้งแต่การสุ่มจนถึงการเสียชีวิตจากทุกสาเหตุเป็นเวลา 24 ถึง 36 เดือน และมีผู้เสียชีวิต 825 คน

การรอดชีพ Kaplan-Meier ที่ 12 เดือนเท่ากับ 0.688 ในกลุ่มควบคุม และ 0.743 ในกลุ่มที่ได้รับการรักษา ที่ 24 เดือนเท่ากับ 0.473 และ 0.569 เส้นทั้งสองลดถึงครึ่งหนึ่งภายในช่วงติดตาม มัธยฐานทั้งสองจึงไปถึง: 22.6 เดือนในกลุ่มควบคุม และ 28.3 เดือนในกลุ่มที่ได้รับการรักษา

ตัวเลขเหล่านี้บรรยายว่าแต่ละกลุ่มเป็นอย่างไร การสุ่มทำให้สองกลุ่มเปรียบเทียบกันได้โดยเฉลี่ย แต่ตัวเส้นเองไม่ได้ปรับอะไรเลย: เป็นค่าประมาณหนึ่งค่าต่อกลุ่ม ไม่ใช่แบบจำลอง

R: การรอดชีพ Kaplan-Meier ที่ 12 และ 24 เดือน

# ---------------------------------------------------------------- all-cause death: Kaplan-Meier
# time origin = randomisation; censoring = end of study or loss to follow-up
km <- survfit(Surv(fu_months, death) ~ arm, data = w3)
s <- summary(km, times = c(12, 24))
print(s)

> km <- survfit(Surv(fu_months, death) ~ arm, data = w3)

> s <- summary(km, times = c(12, 24))

> print(s)
Call: survfit(formula = Surv(fu_months, death) ~ arm, data = w3)

                arm=0
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12    498     230    0.688  0.0171        0.656        0.722
   24    332     153    0.473  0.0186        0.438        0.511

                arm=1
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12    542     190    0.743  0.0161        0.713        0.776
   24    400     125    0.569  0.0184        0.534        0.606

Cox สร้างแบบจำลอง: ฮาซาร์ดสัมพัทธ์

ฮาซาร์ด (hazard) คืออัตราการเกิดเหตุการณ์ ณ เวลา $t$ ในผู้ป่วยที่ยังปลอดเหตุการณ์ เป็นอัตราต่อหนึ่งหน่วยเวลา ไม่ใช่ความน่าจะเป็นที่จะเกิดเหตุการณ์ไปแล้ว อัตราส่วนฮาซาร์ด (hazard ratio, HR) คือฮาซาร์ดของกลุ่มหนึ่งหารด้วยฮาซาร์ดของอีกกลุ่มหนึ่ง

แบบจำลอง Cox ประมาณฮาซาร์ดสัมพัทธ์ (relative hazard) ระหว่างกลุ่ม โดยปรับตัวแปรร่วม (covariates) ได้ [2] ตัวแปรร่วมที่ใส่ในแบบจำลอง เช่นอายุหรือความรุนแรงของหัวใจล้มเหลว จะถูกปรับ ดังนั้น adjusted hazard ratio มาจากแบบจำลอง Cox ไม่ใช่จากกราฟ Kaplan-Meier แบบจำลองไม่ระบุรูปร่างของฮาซาร์ดตามเวลา

ผลลัพธ์ที่พบบ่อย:

หมายความว่าฮาซาร์ดที่ประมาณได้ของกลุ่มหนึ่งต่ำกว่าอีกกลุ่มหนึ่ง ภายใต้ข้อสมมติของแบบจำลอง

การทดลองเดียวกัน สรุปด้วยแบบจำลอง

ในการทดลองจำลอง stcox arm ให้ HR เท่ากับ 0.78 ของกลุ่มที่ได้รับการรักษาเทียบกับกลุ่มควบคุม โดยมีช่วงเชื่อมั่น (confidence interval, CI) 95% คือ 0.68 ถึง 0.90 ฮาซาร์ดของการเสียชีวิตที่ประมาณได้ในกลุ่มที่ได้รับการรักษาต่ำกว่าในกลุ่มควบคุม ภายใต้ข้อสมมติของแบบจำลอง

บนข้อมูลชุดเดียวกัน เครื่องมือทั้งสองทำงานต่างกัน Kaplan-Meier บรรยายแต่ละกลุ่มด้วยเส้น การรอดชีพ ณ เดือนที่ระบุ และมัธยฐาน ส่วน Cox สรุปการเปรียบเทียบเป็นอัตราส่วนค่าเดียวที่ใส่ตัวแปรร่วมได้ด้วย อัตราส่วนค่าเดียวนั้นตั้งอยู่บนข้อสมมติหนึ่ง ซึ่งเป็นเรื่องของส่วนถัดไป

Stata: แบบจำลอง Cox สำหรับการเสียชีวิตจากทุกสาเหตุ

* ---------------------------------------------------------------- all-cause death: Cox model
* hazard ratio for treatment versus control (Breslow ties; there are no tied times)
stcox arm

. * ---------------------------------------------------------------- all-cause death: Cox model
. * hazard ratio for treatment versus control (Breslow ties; there are no tied times)
. stcox arm

        Failure _d: death==1
  Analysis time _t: fu_months

Iteration 0:  Log likelihood = -5648.1152
Iteration 1:  Log likelihood = -5641.8514
Iteration 2:  Log likelihood = -5641.8514
Refining estimates:
Iteration 0:  Log likelihood = -5641.8514

Cox regression with no ties

No. of subjects =       1,500                           Number of obs =  1,500
No. of failures =         825
Time at risk    = 29,562.8078
                                                        LR chi2(1)    =  12.53
Log likelihood = -5641.8514                             Prob > chi2   = 0.0004

------------------------------------------------------------------------------
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         arm |   .7812631   .0545882    -3.53   0.000      .681275    .8959261
------------------------------------------------------------------------------

ข้อสมมติสำคัญ: proportional hazards

แบบจำลอง Cox แบบปกติสมมติว่า hazard ratio คงที่โดยประมาณตลอดเวลา HR ค่าเดียวจึงบรรยายการติดตามทั้งช่วงได้

ถ้ากราฟการรอดชีพบ่งว่าผลของการรักษา

• เปลี่ยนขนาดไปมากตามเวลา
• ตัดกัน (crosses over) หรือ
• กลับทิศทาง (reverses direction)

HR คงที่ค่าเดียวอาจทำให้เข้าใจผิดได้ [3, 4] เมื่ออัตราส่วนเปลี่ยนไปจริง ตัวเลขค่าเดียวที่แบบจำลอง Cox คืนให้คือค่าเฉลี่ยถ่วงน้ำหนักของอัตราส่วนที่เปลี่ยนไป โดยน้ำหนักขึ้นกับว่าเหตุการณ์เกิดเมื่อใดและติดตามผู้ป่วยนานเท่าไร [4]

การทดลองที่ผลกลับทิศ

รุ่นดัดแปลงของการทดลองจำลองถูกสร้างให้ HR ของการเสียชีวิตในกลุ่มที่ได้รับการรักษาเท่ากับ 0.60 ในเดือนที่ 0 ถึง 12 และ 1.10 หลังจากนั้น ผลจึงกลับทิศเมื่อครบหนึ่งปี แบบจำลอง Cox ค่าเดียวที่ fit กับรุ่นดัดแปลงนี้ให้ 0.84 (CI 95% 0.74 ถึง 0.96) ตัวเลขนี้ไม่ได้บรรยายช่วงใดเลย: สูงเกินไปสำหรับปีแรก และต่ำเกินไปสำหรับช่วงหลังจากนั้น

Stata: hazard ratio ค่าเดียวของการทดลองรุ่นดัดแปลง

* one hazard ratio averaged over follow-up, then the PH check
stcox arm

. * one hazard ratio averaged over follow-up, then the PH check
. stcox arm

        Failure _d: death==1
  Analysis time _t: fu_months
       ID variable: id

Iteration 0:  Log likelihood = -5860.3193
Iteration 1:  Log likelihood = -5857.1468
Iteration 2:  Log likelihood = -5857.1468
Refining estimates:
Iteration 0:  Log likelihood = -5857.1468

Cox regression with no ties

No. of subjects =       1,500                           Number of obs =  1,500
No. of failures =         861
Time at risk    = 29,758.7835
                                                        LR chi2(1)    =   6.35
Log likelihood = -5857.1468                             Prob > chi2   = 0.0118

------------------------------------------------------------------------------
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         arm |   .8421581   .0574503    -2.52   0.012     .7367606    .9626332
------------------------------------------------------------------------------

ใช้อะไรเมื่ออัตราส่วนเปลี่ยนไป

ทางเลือกที่เป็นไปได้ ได้แก่

• ผลที่เปลี่ยนตามเวลา (time-varying effects) ยอมให้ HR ต่างกันในแต่ละช่วงของการติดตาม การแบ่งการติดตามของรุ่นดัดแปลงที่เดือนที่ 12 ให้ 0.62 (CI 95% 0.51 ถึง 0.77) ในเดือนที่ 0 ถึง 12 และ 1.07 (CI 95% 0.89 ถึง 1.28) หลังเดือนที่ 12 ใกล้กับ 0.60 และ 1.10 ที่ฝังไว้ในการจำลอง โดยจุดตัดถูกกำหนดไว้ก่อนเห็นข้อมูล
• แบบจำลองการรอดชีพแบบยืดหยุ่น (flexible survival models) สร้างแบบจำลองฮาซาร์ดตามเวลาเป็นเส้นโค้งเรียบ ผลของการรักษาจึงเปลี่ยนไปตามเวลาอย่างราบรื่นได้ แทนที่จะถูกบังคับให้เป็นค่าคงที่ค่าเดียว
• ระยะเวลารอดชีพเฉลี่ยแบบจำกัดช่วง (restricted mean survival time, RMST) ที่ขอบเขตเวลาที่กำหนดไว้ล่วงหน้า พื้นที่ใต้กราฟการรอดชีพจนถึงขอบเขตเวลา (horizon) $\tau$ ที่เลือกไว้ก่อน คือเวลาปลอดเหตุการณ์เฉลี่ยภายในช่วงนั้น [5] ในรุ่นดัดแปลง จนถึง 24 เดือน ประมาณได้ 16.9 เดือนในกลุ่มควบคุม และ 18.7 เดือนในกลุ่มที่ได้รับการรักษา ตอนที่ 3 อธิบายเรื่องนี้โดยละเอียด

Stata: hazard ratio หนึ่งค่าต่อช่วงเวลาของการทดลองรุ่นดัดแปลง

* split follow-up at 12 months and give each period its own treatment hazard ratio
stsplit late, at(12)
generate byte arm_0_12 = arm * (late == 0)
generate byte arm_12plus = arm * (late == 12)
stcox arm_0_12 arm_12plus

. * split follow-up at 12 months and give each period its own treatment hazard ratio
. stsplit late, at(12)
(1,075 observations (episodes) created)

. generate byte arm_0_12 = arm * (late == 0)

. generate byte arm_12plus = arm * (late == 12)

. stcox arm_0_12 arm_12plus

        Failure _d: death==1
  Analysis time _t: fu_months
       ID variable: id

Iteration 0:  Log likelihood = -5860.3193
Iteration 1:  Log likelihood = -5849.6257
Iteration 2:  Log likelihood = -5849.6247
Refining estimates:
Iteration 0:  Log likelihood = -5849.6247

Cox regression with no ties

No. of subjects =       1,500                           Number of obs =  2,575
No. of failures =         861
Time at risk    = 29,758.7835
                                                        LR chi2(2)    =  21.39
Log likelihood = -5849.6247                             Prob > chi2   = 0.0000

------------------------------------------------------------------------------
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    arm_0_12 |   .6245692   .0650446    -4.52   0.000     .5092532    .7659976
  arm_12plus |   1.067557   .0982992     0.71   0.478     .8912784    1.278699
------------------------------------------------------------------------------

จากถ้อยคำในรายงานสู่วิธีการ

อ่านต่อในวิกิ: [[time-to-event-survival-analysis-th]]