Immortal Time Bias

Epilogue

วันก่อนจะอ่านpaperเรื่อง Statin & new onset DM in elderly women แล้วพบกับสิ่งที่เรียกว่า Immortal time bias ครับ Bias ชนิดนี้ผมเคยอ่านมาแล้วครั้งหนึ่ง พอเข้าใจ concept อยู่ แต่ยอมรับว่าสังเกตได้ยากเหมือนกัน ซึ่งไม่แปลกเพราะมีคนพบ bias ชนิดนี้ซ่อนอยู่ในหลายๆงานวิจัยที่เป็น observation โดยเฉพาะที่ใช้ survival analysis
จุดมุ่งหมายของบทความนี้คือความพยายามในการ“แสดง”ให้เห็นที่มาของมันอย่างชัดเจน รวมถึงวิธีการแก้ไข โดยทั้งการสร้างบทความนี้และcodeทั้งหมด ทำใน R ด้วยเครื่องมือที่เรียกว่า R Markdown นะครับ มาลองติดตามกัน

Classic Case

ตัวอย่างที่คลาสสิกสุดๆ คงเป็น cohort ที่เปรียบเทียบคนไข้ที่เข้าคิวทำ heart transplant แล้วพบว่า คนที่ได้ทำมีระยะเวลาการรอดชีวิตสูงกว่ากลุ่มที่ไม่ได้ทำอย่างมีนัยสำคัญ โดยนับ“เวลาที่มีชีวิต” ตั้งแต่วันที่เข้าคิวจะ transplant จนถึงวันตาย … ฟังดูก็ตรงไปตรงมาและน่าจะเชื่อได้ใช่ไหมครับ? ลองมาดูข้อมูลนี้ที่ผมทำเทียมขึ้นมานะครับ
คำอธิบาย: step1: เป็นการจำลองข้อมูลของคน 50 คนที่มี hazard = 5% แปลว่าในทุกๆ 1 หน่วยเวลา จะมีโอกาสเกิด event คือ ตาย ได้ 5%(เช่น วัน เดือน หรือปี ในกรณีนี้คงนับเป็นวันก็ได้

  library(survival)
  set.seed(2010)
  hazard=0.05
  subject <- data.frame(id=1:50,time=0,event=1,transplant=0)
  for(i in 1:50){
    while(runif(1,0,1)>=hazard)
      subject$time[i]<-subject$time[i]+1
  }

step2: ดังนั้นทั้งกลุ่มนี้ถ้าเราสุ่มๆ แบ่งคนมา 25:25 ก็น่าจะมี survival พอๆกัน มีระยะเวลามีชีวิตรอดพอๆกัน … จากข้อมูลแรก ผมจะลองสุ่มคนมา 25 คน ให้ได้ทำ heart transplant แล้วแสดงกราฟข้อมูลของคนทั้ง 50 คน แกนXคือระยะเวลาตั้งแต่เข้าคิวจนถึงวันตาย แต่ละเส้น = คน 1 คน นะครับ โดย เส้นสีน้ำเงิน = ได้ transplant, สีแดง = ไม่ได้ทำนะครับ

  lucky.draw <- sample(1:50,25,replace=F)
  subject$transplant[lucky.draw]<- 1
  subject<-subject[order(subject$transplant,subject$time),]
  subject$id <- 1:50
  barplot(subject$time,horiz=TRUE,main="Plot#1 True Random Transplant",
          col=ifelse(subject$transplant==1,"blue","red"))

  plot(survfit(Surv(subject$time,subject$event)~subject$transplant),
       col=c("red","blue"))

step2X: แต่ในความเป็นจริง การได้ทำ transplant ไม่ใช่ว่าจะทำก็ได้ทำเลย สมมติว่าต้องมีเวลาเตรียมตัวประมาณ 10 วันตั้งแต่เข้าคิว ผมจะลองเลือกคนที่ได้ทำ transplant ใหม่ โดยโชคดี จะตกอยู่กับคนที่รอคิวมาเกิน 10 วันเท่านั้นนะครับ

  subject$transplant <- rep(0,50)
  list10 <- subject[subject$time>=10,"id"]
  lucky.draw <- sample(list10,25,replace=F)
  subject$transplant[lucky.draw]<- 1
  subject<-subject[order(subject$transplant,subject$time),]
  barplot(subject$time,horiz=TRUE,main="Plot#2 Real World Transplant",
          col=ifelse(subject$transplant==1,"blue","red"))
  abline(v=10,lty=2)

  plot(survfit(Surv(subject$time,subject$event)~subject$transplant),
       col=c("red","blue"))

ถ้าคิดตามมาเรื่อยๆอย่างนี้ ทุกคนคงจะเห็นได้ว่า ที่กลุ่มสีน้ำเงินที่ได้ transplant ดูมีชีวิตยืนยาวกว่านั้นก็เพราะ พวกเขาต้องผ่าน“การทดสอบ”มาแล้ว 10 วัน จึงจะมีโอกาสได้ทำ transplant จริงๆ … แปลว่าถ้าใครตายก่อน 10 วันก็จะอดได้ทำ transplant, พูดกลับกันได้ว่า คนที่ได้ทำ transplant ไม่มีวันจะตายก่อน 10 วัน … ช่วงเวลา 10 วันนี้นี่แหละครับ ที่เราเรียกว่า “Immortal Time”

Immortal time bias ทำให้เกิดอะไร?

ปัจจัยที่สำคัญก็คือ index date หรือวันที่เริ่มนับว่าแต่ละ subject เข้าในการศึกษา กับ exposure determination date หรือวันที่แต่ละ subject เข้านิยามที่กำหนดว่า exposure status ของ subject นั้นเป็นอย่างไร จากตัวอย่างเดิม แท่งสีเหลืองคือ immortal time ที่เกิดขึ้น

  subject$txdate <- ifelse(subject$transplant==1,10,0)
  barplot(subject$time,horiz=TRUE,main="Plot#3 Immortal time",
          col=ifelse(subject$transplant==1,"blue","red"))
  abline(v=10,lty=2)
  barplot(height=subject$txdate,horiz=TRUE,col="yellow",add=TRUE)

ทีนี้ปัญหาในการแปลผลเกิดจากขั้นตอนการวิเคราะห์ ซึ่งมี 2 แบบหลักๆ
1. นับ immortal time ไปด้วย แบบนี้จะทำให้กลุ่มที่มี exposure มีโอกาสเกิด outcome ช้ากว่ากลุ่มที่ไม่มี exposure
2. ตัด immortal time ทิ้งไปเลย กล่าวคือเลื่อน index date ของกลุ่ม exposure ไปที่วันที่ได้รับ exposure ซึ่งส่งผลได้หลายแบบ
…2.1 ในกลุ่ม exposure มีแนวโน้มจะทำให้เกิด outcome เร็วขึ้นเมื่อเทียบกับกลุ่ม control ที่มี index date ตามเดิม
…2.2 จริงๆแล้วช่วงก่อนที่จะมี exposure เกิดขึ้น ก็คือช่วงที่เหมือน control นั่นเอง แถมเป็นการรอดชีวิตของกลุ่ม control ด้วยนะ แต่กลับไม่เอามาคิดเป็นคะแนนบวกให้กลุ่ม control (แต่ถ้าตายในช่วงนี้ จะจัดเป็น control โดยอัตโนมัติ แถมคิดแต้มลบ) date ตามเดิม

  subject$time <- subject$time-subject$txdate
  barplot(subject$time,horiz=TRUE,main="Plot#4 Remove immortal time",
          col=ifelse(subject$transplant==1,"blue","red"))

  plot(survfit(Surv(subject$time,subject$event)~subject$transplant),
       col=c("red","blue"))

ซึ่งในตัวอย่างนี้จะเห็นว่าเมื่อตัด immortal time ทิ้งไป กลุ่มtransplantก็มีsurvivalลดลง แต่ไม่ถึงกับห่วยกว่าอีกกลุ่มนะครับ ความรุนแรงของผลกระทบนี้ก็แปรไปตามแต่ละการศึกษา

What are conditions leading to this bias ?

มีหลาย design ที่ทำให้เกิด bias แบบนี้ได้
1. index date ยึดตามเวลาที่ fix และ exposure เกิดตามมาทีหลัง
2. index date ยึดตามเหตุการณ์บางอย่าง เช่น admit, discharge, diagnosis และ exposure เกิดตามมาทีหลัง
3. index date ยึดตาม definition ของ exposure ซึ่งมีความเป็นลำดับชั้นกัน เช่น ได้กินยา6เดือนเทียบกับกินยา3เดือน (ส่วนต่าง 3 เดือนคือ immortal time) ใช้ยา1ตัวเทียบกับยา2ตัว (ถ้ายาตัวที่ 2 มักเป็น add-on therapy แปลว่าช่วงระหว่างรอดูผลของยาตัวแรกคือ immortal time) หรือนิยามการ expose ว่าต้องได้รับการสั่งยา 2 ครั้ง (ระยะห่างระหว่าง2ครั้งก็เป็น immortal time)
4. อันนี้ชวนงงที่สุด คือกรณีที่ expose มีการถอยเข้าๆออกๆ ระหว่างการติดตามได้ แต่เราคิดโดยอาศัยเกณฑ์ตอนแรกเข้าเท่านั้น เช่น index date เป็นวันวินิจฉัยRheumatoid แล้วเราเทียบตับอักเสบในคนที่ได้MTXตั้งแต่แรก กับไม่ได้ตั้งแต่แรก จะเห็นได้ว่า กลุ่มที่ไม่ได้ตั้งแต่แรกนั้นในอนาคตอาจจะได้MTXก็เป็นได้ และถ้ากลุ่มcontrolนี้ตับอักเสบหลังจากได้MTX ก็จะจัดว่าเป็นตับอักเสบในกลุ่ม control ถือได้ว่าช่วงเวลาหลังได้ delay MTX นี้เป็น immortal time ของกลุ่มMTXยังไงยังงั้น

How to eliminate this bias?

จริงๆแล้วมีหลายวิธี วิธีที่เข้าใจง่ายที่สุดน่าจะเป็นการ design ที่ดีไม่ให้เกิด immortal time ขึ้น เช่น ถ้าเราจะนับว่าexposeเมื่อได้ยา6เดือน ก็ต้องนับ index date ของ exposure ในวันที่ได้ยาครบ 6 เดือน และกลุ่ม control ที่ไม่ได้ exposure นั้น index date ต้องหาทางทำยังไงก็ได้ ให้ match กับกลุ่ม exposure เช่น ต้องเป็นคนที่ไม่ได้ยาและมีระยะของโรคเท่ากัน มีprognosisเท่าๆกัน (เหมือนตัด immortal time ของกลุ่ม exposure ออกไปพร้อมๆกับ mortal time ของกลุ่ม control ด้วย) ข้อเสียคือพูดง่ายแต่ทำยาก และ N จะลดลงอย่างมาก

วิธีการที่ดีกว่าจึงเป็นการวิเคราะห์ทางสถิติแบบที่เรียกว่ามี time-dependent covariate โดยหลักการคือ ถ้า exposure status (รวมถึงปัจจัยอื่นๆที่เราคิดว่ามีผลต่อ outcome) มีการเปลี่ยนแปลงไม่คงที่ เช่น กรณีของ transplant ถ้าเรายังคงนับ index date เป็นวันที่มาเข้าคิว ช่วงก่อนที่คนคนหนึ่งจะได้ทำ transplant ถ้าเขารอดชีวิต (หรือตาย) ก็ควรนับเป็นสิ่งที่เกิดในกลุ่มที่ยังไม่ได้ transplant
ตัวอย่างลองคิดถึงคนคนนึงมีชีวิตอยู่หลังจากเข้าคิว transplant มาได้ 10 ปี (เว่อมาก) แล้วในที่สุดเขาก็ได้ทำ transplant แล้วตายวันนั้นเลย … สามัญสำนึกเราคงบอกได้ว่าระยะเวลามีชีวิตรอด 10 ปีนี้ ควรนำมาคิดในกลุ่มไม่ได้ทำ transplant มากกว่า … หลังจากได้ทำ transplant แล้วนั่นแหละจึงจะนำมาคิดว่าเป็นเหตุการณ์ในกลุ่ม transplant
ข้างล่างนี้เป็นตัวอย่างตารางเก็บข้อมูลแบบที่มี time-dependent covariate สังเกต subject id#4 จะมีการ split ช่วงเวลาออกเป็น 2 ช่วงคือก่อนและหลัง transplant ครับ

jasa1[1:10,1:5]

##     id start stop event transplant
## 1    1     0   49     1          0
## 2    2     0    5     1          0
## 102  3     0   15     1          1
## 3    4     0   35     0          0
## 103  4    35   38     1          1
## 4    5     0   17     1          0
## 5    6     0    2     1          0
## 6    7     0   50     0          0
## 104  7    50  674     1          1
## 7    8     0   39     1          0

ในอีกหน้าหนึ่งเป็น code R ที่ผมทำไว้ให้ตัวเองอ่านว่าจะทำตารางเก็บข้อมูลแบบนี้ได้อย่างไร พร้อมคำบรรยายภาษาอังกฤษ ท่านใดสนใจก็ลองไปติดตามดูได้ครับ ข้อมูลทั้งหมดใช้ dataset ตัวอย่างที่อยู่ใน R package survival อยู่แล้ว