อัตรา ความน่าจะเป็น และความน่าจะเป็นของการเปลี่ยนสถานะ: คำนวณโดยไม่หลงหน่วยเวลา
On this page
Read the English version
บทคัดย่อ
คลินิกโรคไตแห่งหนึ่งนับการเริ่มฟอกไตครั้งแรกได้ 120 ครั้ง ในการติดตาม 1,000 คน-ปี ของผู้ใหญ่ที่เป็นโรคไตเรื้อรังระยะที่ 4 คิดเป็นอัตรา 0.12 ต่อคน-ปี แบบจำลองของคลินิกย้ายผู้ป่วยทุก 3 เดือน คือรอบยาว 0.25 ปี และต้องการโอกาสที่จะเริ่มฟอกไตภายในหนึ่งรอบ ทางลัดที่ชวนให้ใช้คือคูณสองค่านี้ 0.12 × 0.25 = 0.0300 แต่ทางลัดนี้ให้ค่าสูงเกินจริง เพราะผู้ป่วยที่เริ่มฟอกไตแล้วออกจากกลุ่มที่ยังเสี่ยงอยู่ โอกาสที่ถูกต้องคือหนึ่งลบสัดส่วนที่ยังไม่ได้ฟอกไต ได้ 0.0296 โอกาสที่จะย้ายจากสถานะหนึ่งไปอีกสถานะหนึ่งภายในหนึ่งรอบนี้ เรียกว่าความน่าจะเป็นของการเปลี่ยนสถานะ (transition probability) ในห้าปีช่องว่างจะกว้างขึ้น คือ 45.12% เทียบกับ 60.00% ความเสี่ยงห้าปีต้องแปลงกลับผ่านอัตรา ห้ามหารด้วยจำนวนปี และผลของการรักษาต้องคูณกับอัตรา ไม่ใช่คูณกับโอกาส เมื่อการเสียชีวิตแข่งกับการฟอกไต ทางออกทั้งสองใช้ผู้ป่วยกลุ่มเดียวกัน โอกาสที่จะอยู่ต่อ เริ่มฟอกไต และเสียชีวิตจึงรวมได้หนึ่งพอดี
คลินิกเดียว อัตราเดียว และแบบจำลองที่ขยับทุก 3 เดือน
คลินิกโรคไตในโรงพยาบาลจังหวัดแห่งหนึ่งติดตามผู้ใหญ่ที่เป็นโรคไตเรื้อรังระยะที่ 4 ซึ่งเป็นระยะก่อนไตวาย บันทึกของคลินิกรวมได้ 1,000 คน-ปี โดยผู้ป่วยหนึ่งคนที่ติดตามนานหนึ่งปีนับเป็นหนึ่งคน-ปี ในช่วงนั้นมีผู้ป่วย 120 คนเริ่มฟอกไตเป็นครั้งแรก เมื่อหารจำนวนเหตุการณ์ด้วยเวลาที่เสี่ยง จะได้ 120 / 1,000 = 0.12 ครั้งต่อคน-ปี
ทีมของคลินิกกำลังสร้างแบบจำลองที่แบ่งผู้ป่วยเป็นสถานะสุขภาพไม่กี่สถานะ คือยังอยู่ระยะที่ 4 ฟอกไตแล้ว และเสียชีวิต แบบจำลองย้ายผู้ป่วยระหว่างสถานะทุก 3 เดือน ตามช่วงนัดของคลินิก ช่วงเวลาคงที่นี้เรียกว่ารอบ (cycle) ในที่นี้ยาว 0.25 ปี ในแต่ละรอบ แบบจำลองต้องการตัวเลขหนึ่งตัว คือโอกาสที่ผู้ป่วยระยะที่ 4 จะเริ่มฟอกไตก่อนนัดครั้งถัดไป
คำตอบที่ชวนให้ใช้คือคูณอัตราด้วยความยาวรอบ 0.12 × 0.25 = 0.0300 หรือ 30 คนจากผู้ป่วยทุก 1,000 คน บทความนี้แก้นิสัยนั้น อัตรากับโอกาสเป็นปริมาณต่างชนิดกัน และเวลาต้องอยู่ในการแปลงด้วย ตัวเลขทั้งหมดเป็นตัวอย่างสมมติเพื่อการสอน
อัตราคือความเร็ว โอกาสคือสัดส่วน
อัตรานับว่าเหตุการณ์เกิดเร็วแค่ไหน คือจำนวนเหตุการณ์หารด้วยเวลาที่คนเสี่ยง ในที่นี้คือ 0.12 ครั้งต่อคน-ปี ส่วนโอกาส ซึ่งเรียกว่าความน่าจะเป็นหรือความเสี่ยงก็ได้ คือสัดส่วนของกลุ่ม คือส่วนของคนที่เกิดเหตุการณ์แล้วเมื่อสิ้นช่วงเวลาที่กำหนด สัดส่วนจะเกิน 1 ไม่ได้ เพราะกลุ่มหนึ่งจะเริ่มฟอกไตครั้งแรกมากกว่าจำนวนผู้ป่วยไม่ได้ ส่วนอัตราไม่มีเพดานแบบนั้น
การคูณอัตราด้วยเวลาจึงไม่ได้สัดส่วน ในช่วง 10 ปี ทางลัดให้ 0.12 × 10 = 1.20 ซึ่งมากกว่าจำนวนผู้ป่วยทั้งหมดของคลินิก อัตรายังเป็นความเร็วที่ดีอยู่ เพียงแต่ไม่ใช่สัดส่วน
อัตรา ณ ขณะหนึ่งในกลุ่มผู้ป่วยที่ยังไม่เกิดเหตุการณ์เรียกว่าฮาซาร์ด (hazard) บทความนี้สมมติว่าฮาซาร์ดคงที่ คือค่า 0.12 ของคลินิกใช้ได้ทุกปีของการติดตาม นี่เป็นข้อสมมติที่ต้องตรวจกับข้อมูล ไม่ใช่ข้อเท็จจริง (Fleurence and Hollenbeak)
ทำไมอัตรา × เวลาจึงสูงเกินจริง
ผู้ป่วยแต่ละคนเริ่มฟอกไตครั้งแรกได้เพียงครั้งเดียว ผู้ป่วยที่เริ่มแล้วจะออกจากกลุ่มที่ยังเสี่ยง อัตราเดิมจึงใช้กับกลุ่มที่เล็กลงเล็กน้อยทุกสัปดาห์ ทางลัดไม่สนใจเรื่องนี้และบวกผู้เริ่มฟอกไตใหม่ต่อไปเหมือนไม่มีใครออกเลย วิธีที่ถูกต้องจึงติดตามสัดส่วนที่ยังไม่ได้ฟอกไตแทน
สัดส่วนนี้ไม่ได้ลดลงเป็นเส้นตรง ในทุกช่วงเวลาที่ยาวเท่ากัน สัดส่วนนี้ลดลงเป็นสัดส่วนเท่าเดิมของคนที่ยังเหลืออยู่ แต่ละช่วงจึงหักออกน้อยกว่าช่วงก่อนเล็กน้อย เครื่องคิดเลขคำนวณเรื่องนี้ด้วยปุ่ม ex ซึ่งเป็นฟังก์ชันเอ็กซ์โพเนนเชียล (Briggs and colleagues) ตัวอย่างคำนวณด้านล่างแทนตัวเลขของคลินิกลงไป
ตัวอย่างคำนวณ: โอกาสต่อรอบ 3 เดือนของคลินิก
ใช้อัตราของคลินิก 0.12 ครั้งต่อคน-ปี และหนึ่งรอบยาว 0.25 ปี
-
คูณอัตราด้วยเวลา
\[ 0.12 \times 0.25 = 0.03 \]
ผลคูณนี้คือค่าอัตราที่สะสมครบหนึ่งรอบ ยังไม่ใช่โอกาส ทางลัดหยุดตรงนี้และเรียก 0.0300 ว่าโอกาส
-
หาสัดส่วนที่ยังไม่ได้ฟอกไต
\[ e^{-0.03} = 0.9704 \]
การหักสัดส่วนเท่าเดิมของคนที่เหลือทีละช่วงสั้น ๆ คูณกันแล้วได้ e ยกกำลังลบของค่าอัตราที่สะสมในรอบ คือกดปุ่ม e^x ด้วย −0.03 จากผู้ป่วยระยะที่ 4 จำนวน 1,000 คนเมื่อเริ่มรอบ มี 970.4 คนที่ยังไม่ได้ฟอกไตเมื่อถึงนัดครั้งถัดไป
-
หาส่วนที่เหลือ
\[ 1 - 0.9704 = 0.0296 \]
ดังนั้นผู้ป่วย 29.6 คนจากทุก 1,000 คนเริ่มฟอกไตในรอบนี้ เทียบกับ 30 คนของทางลัด
-
ตรวจโดยแปลงย้อนกลับ
\[ -\ln(0.9704) / 0.25 \approx 0.12 \]
ปุ่ม ln บนเครื่องคิดเลขย้อนกลับปุ่ม e^x ลอการิทึมธรรมชาติของสัดส่วนที่ยังไม่ได้ฟอกไต เมื่อกลับเครื่องหมายแล้วหารด้วยความยาวรอบ จะได้อัตราของคลินิกคืนมา คือ 0.12 เมื่อปัดสองตำแหน่ง ส่วน R ได้ 0.12 พอดีเพราะใช้ค่าที่ไม่ปัดเศษ ซึ่งยืนยันว่าคำนวณถูก
ผลลัพธ์: โอกาสต่อรอบของแบบจำลองคือ 0.0296 ต่ำกว่า 0.0300 ของทางลัดเล็กน้อย ในรอบสั้นหนึ่งรอบทางลัดยังใกล้เคียง หัวข้อถัดไปจะแสดงว่าเมื่อไรที่ไม่ใกล้
ชื่อเรียก: ความน่าจะเป็นของการเปลี่ยนสถานะ
ค่า 0.0296 ของคลินิกคือตัวเลขที่แบบจำลองต้องการ คือโอกาสที่จะย้ายจากสถานะหนึ่งคือระยะที่ 4 ไปอีกสถานะหนึ่งคือการฟอกไต ภายในหนึ่งรอบยาว 0.25 ปี โอกาสนี้เรียกว่าความน่าจะเป็นของการเปลี่ยนสถานะ (transition probability) ซึ่งแปลงมาจากอัตราที่นับเหตุการณ์ต่อหน่วยเวลา โดยให้เวลาอยู่ในการคำนวณด้วย (Siebert and colleagues)
เขียนเป็นสัญลักษณ์ได้ว่า ให้ $S$ แทนสัดส่วนที่ยังไม่ได้ฟอกไตเมื่อผ่านเวลา $t$ ที่อัตราคงที่ $r$ จะได้ $S = e^{-rt}$ อัตรา $r$ ในรอบยาว $t$ จึงให้ความน่าจะเป็นของการเปลี่ยนสถานะ $p = 1 - e^{-rt}$ และแปลงความน่าจะเป็นกลับเป็นอัตราได้ด้วย $r = -\ln(1 - p) / t$
ตรวจตัวเลขของหนึ่งรอบด้วย R
# Illustrative clinic counts: first starts of dialysis among people with stage 4 kidney disease.
# Writes: nothing to disk; the results print in the console.
starts <- 120
person_years <- 1000
cycle_years <- 0.25
rate <- starts / person_years
still_free <- exp(-rate * cycle_years)
per_thousand <- 1000 * c(rate * cycle_years, 1 - still_free, still_free)
cat("Rate per person-year:", rate, " shortcut over 10 years:", sprintf("%.2f", rate * 10), "\n")
cat("Shortcut for one 3-month cycle:", formatC(rate * cycle_years, format = "f", digits = 4), "\n")
cat("Share still free, exp(-0.03):", signif(still_free, 4), "\n")
cat("Cycle probability, 1 minus that share:", signif(1 - still_free, 3), "\n")
cat("Of 1,000 at stage 4, shortcut, starts, still free:", round(per_thousand, 1), "\n")
cat("Rate recovered from the probability:", -log(still_free) / cycle_years, "\n")Rate per person-year: 0.12 shortcut over 10 years: 1.20
Shortcut for one 3-month cycle: 0.0300
Share still free, exp(-0.03): 0.9704
Cycle probability, 1 minus that share: 0.0296
Of 1,000 at stage 4, shortcut, starts, still free: 30 29.6 970.4
Rate recovered from the probability: 0.12ต้องใช้ R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · ใช้ได้กับ R 4.6.1 · ดาวน์โหลด rates-probabilities-and-transition-probabilities.R
สคริปต์ฉบับเต็ม
# Rates, Probabilities and Transition Probabilities: Keep Time in the Calculation
# Follow one kidney clinic from a counted rate to the chances a 3-month model needs.
# ติดตามคลินิกโรคไตหนึ่งแห่ง จากอัตราที่นับได้ไปสู่ความน่าจะเป็นที่แบบจำลองรอบละ 3 เดือนต้องใช้
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the inputs and the exits as explicit tables.
library(dplyr) # Calculate the converted probabilities.
library(ggplot2) # Display the exact curve against the shortcut without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
starts <- 120
person_years <- 1000
cycle_years <- 0.25
rate <- 0.12
years <- 5
hazard_ratio <- 0.70
rate_dialysis <- 0.12
rate_death <- 0.08
stopifnot(abs(starts / person_years - rate) < 1e-12)
decimals <- function(value, digits) formatC(value, format = "f", digits = digits)
percent <- function(value, digits) paste0(formatC(100 * value, format = "f", digits = digits), "%")
#%%@note 2. Convert the rate over one cycle and over five years
# Survival falls as exp(-rate x time); the probability is one minus that survival.
# สัดส่วนที่ยังไม่เกิดเหตุการณ์ลดลงแบบ exp(-อัตรา x เวลา) ความน่าจะเป็นคือหนึ่งลบสัดส่วนนั้น
still_free <- exp(-rate * cycle_years)
per_thousand <- 1000 * c(rate * cycle_years, 1 - still_free, still_free)
risk_five <- 1 - exp(-rate * years)
rate_back <- -log(1 - risk_five) / years
stopifnot(abs(rate_back - rate) < 1e-12)
yearly <- c(risk_five / years, 1 - exp(-rate_back))
per_cycle <- c(risk_five / 20, 1 - exp(-rate_back * cycle_years))
#%%@note 3. Apply the treatment effect to the rate, then convert
# The trial's ratio multiplies the rate, never the probability; the foundation uses survival 0.60.
# อัตราส่วนจากการทดลองคูณกับอัตรา ไม่ได้คูณกับความน่าจะเป็น ตัวอย่างพื้นฐานใช้การรอดชีพ 0.60
risk_hr <- 1 - 0.60^c(1, hazard_ratio)
treated <- rate * hazard_ratio
five_year <- 1 - exp(-c(rate, treated) * years)
#%%@note 4. Share one survival between two exits
# Add the exit rates, find who stays, then split the leavers by each rate's share.
# รวมอัตราของทางออกทั้งสอง หาสัดส่วนที่อยู่ต่อ แล้วแบ่งผู้ที่ออกตามสัดส่วนของแต่ละอัตรา
exits <- tibble::tibble(exit = c("dialysis", "death"), rate = c(rate_dialysis, rate_death)) |>
dplyr::mutate(separate = 1 - exp(-rate * cycle_years))
rate_total <- sum(exits$rate)
stay <- exp(-rate_total * cycle_years)
row_joint <- c(stay, (1 - stay) * exits$rate / rate_total)
separate <- exits$separate
stopifnot(abs(sum(row_joint) - 1) < 1e-12)
#%%@note 5. Display the model
# Plot the exact risk against rate x time over ten years; no graphics file is saved.
# แสดงความเสี่ยงที่ถูกต้องเทียบกับอัตราคูณเวลาในสิบปี โดยไม่บันทึกไฟล์ภาพ
curve <- tibble::tibble(time = seq(0, 10, by = 0.25)) |>
dplyr::mutate(exact = 1 - exp(-rate * time), shortcut = rate * time)
print(ggplot2::ggplot(curve, ggplot2::aes(x = time)) +
ggplot2::geom_line(ggplot2::aes(y = exact)) +
ggplot2::geom_line(ggplot2::aes(y = shortcut), linetype = "dashed") +
ggplot2::labs(x = "Years", y = "Share who started dialysis"))
#%%@note 6. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Rate per person-year:", rate, " shortcut over 10 years:", sprintf("%.2f", rate * 10), "\n")
cat("Shortcut for one 3-month cycle:", formatC(rate * cycle_years, format = "f", digits = 4), "\n")
cat("Share still free, exp(-0.03):", signif(still_free, 4), "\n")
cat("Cycle probability, 1 minus that share:", signif(1 - still_free, 3), "\n")
cat("Of 1,000 at stage 4, shortcut, starts, still free:", round(per_thousand, 1), "\n")
cat("Rate recovered from the probability:", -log(still_free) / cycle_years, "\n")
cat("Five-year risk, exact and rate x time:", percent(c(risk_five, rate * years), 2), "\n")
cat("Five-year risk and share still free:", signif(c(risk_five, 1 - risk_five), 4), "\n")
cat("Rate recovered from the share:", rate_back, "\n")
cat("One year, divided by 5 and through the rate:", percent(yearly, 2), "\n")
cat("One cycle, divided by 20 and through the rate:", signif(per_cycle, 3), "\n")
cat("Survival 0.60, risks and ratio:", percent(risk_hr, 1), round(risk_hr[2] / risk_hr[1], 2), "\n")
cat("Treated rate, cycle probability:", treated, signif(1 - exp(-treated * cycle_years), 3), "\n")
cat("Five-year risks, control and treated:", percent(five_year, 2), "\n")
cat("Five-year risk ratio:", sprintf("%.3f", five_year[2] / five_year[1]), "\n")
cat("Total rate, each exit's share:", rate_total, c(rate_dialysis, rate_death) / rate_total, "\n")
cat("Stay, dialysis, death:", decimals(row_joint, 4), " row sum:", sum(row_joint), "\n")
cat("Leave in one cycle, joint and separate:", decimals(c(1 - stay, sum(separate)), 4), "\n")
cat("Separate dialysis and death:", decimals(separate, 4), "\n")
เปลี่ยนช่วงเวลาผ่านอัตรา ไม่ใช่ด้วยการหาร
ในหนึ่งรอบ ทางลัดยังใกล้เคียง คือ 30 เทียบกับ 29.6 ต่อผู้ป่วย 1,000 คน แต่เมื่อยืดเวลาออกไป ช่องว่างจะกว้างขึ้น ในช่วง 5 ปี สัดส่วนที่เริ่มฟอกไตจริงคือ $1 - e^{-(0.12 \times 5)}$ = 45.12% ขณะที่อัตรา × เวลาให้ 60.00%
ความเสี่ยงห้าปีที่มีผู้รายงาน แปลงกลับเป็นหนึ่งรอบ
ทีนี้ลองคิดย้อนทาง สมมติว่าคลินิกไม่มีข้อมูลนับของตัวเอง และทะเบียนผู้ป่วยระดับประเทศรายงานเพียงว่า ผู้ป่วยระยะที่ 4 ที่เริ่มฟอกไตภายใน 5 ปีของการติดตามคิดเป็น 45.12% เราเลือกตัวเลขนี้ให้ตรงกับอัตราของคลินิก เพื่อเทียบสองวิธีได้ ถ้าหารด้วย 20 รอบใน 5 ปี จะได้ 0.0226 ต่อรอบ
วิธีที่ถูกต้องคือย้อนผ่านอัตรา ขั้นแรกแปลงความเสี่ยงห้าปีเป็นอัตรา สัดส่วนที่ยังไม่ได้ฟอกไตเมื่อครบ 5 ปีคือ 1 − 0.4512 = 0.5488 หาลอการิทึมธรรมชาติของค่านี้ กลับเครื่องหมาย แล้วหารด้วย 5 คือ $-\ln(0.5488) / 5$ = 0.12 ต่อคน-ปี จากนั้นแปลงอัตรานั้นเป็นหนึ่งรอบเหมือนเดิม ได้ 0.0296
การหารจึงทำให้โอกาสต่ำไปเกือบหนึ่งในสี่ เพราะอะไร อัตรานี้ใช้กับเฉพาะผู้ป่วยที่ยังไม่ได้ฟอกไต และกลุ่มนี้หดลงเรื่อย ๆ ปีแรก ๆ จึงมีผู้เริ่มฟอกไตมากกว่าปีหลัง ๆ การเกลี่ย 45.12% เท่า ๆ กันทำให้ทุกปีได้ส่วนเท่ากัน และให้ปีแรกน้อยเกินไป คือ 9.02% เทียบกับค่าที่ถูกต้อง 11.31%
ตรวจตัวเลขห้าปีด้วย R
# The same illustrative rate read over five years, then brought back to one year and one cycle.
# A long-interval risk is converted back through the rate, never divided by the years.
rate <- 0.12
years <- 5
percent <- function(value, digits) paste0(formatC(100 * value, format = "f", digits = digits), "%")
risk_five <- 1 - exp(-rate * years)
rate_back <- -log(1 - risk_five) / years
cat("Five-year risk, exact and rate x time:", percent(c(risk_five, rate * years), 2), "\n")
cat("Five-year risk and share still free:", signif(c(risk_five, 1 - risk_five), 4), "\n")
cat("Rate recovered from the share:", rate_back, "\n")
yearly <- c(risk_five / years, 1 - exp(-rate_back))
cat("One year, divided by 5 and through the rate:", percent(yearly, 2), "\n")
per_cycle <- c(risk_five / 20, 1 - exp(-rate_back * 0.25))
cat("One cycle, divided by 20 and through the rate:", signif(per_cycle, 3), "\n")Five-year risk, exact and rate x time: 45.12% 60.00%
Five-year risk and share still free: 0.4512 0.5488
Rate recovered from the share: 0.12
One year, divided by 5 and through the rate: 9.02% 11.31%
One cycle, divided by 20 and through the rate: 0.0226 0.0296ต้องใช้ R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · ใช้ได้กับ R 4.6.1 · ดาวน์โหลด rates-probabilities-and-transition-probabilities.R
สคริปต์ฉบับเต็ม
# Rates, Probabilities and Transition Probabilities: Keep Time in the Calculation
# Follow one kidney clinic from a counted rate to the chances a 3-month model needs.
# ติดตามคลินิกโรคไตหนึ่งแห่ง จากอัตราที่นับได้ไปสู่ความน่าจะเป็นที่แบบจำลองรอบละ 3 เดือนต้องใช้
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the inputs and the exits as explicit tables.
library(dplyr) # Calculate the converted probabilities.
library(ggplot2) # Display the exact curve against the shortcut without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
starts <- 120
person_years <- 1000
cycle_years <- 0.25
rate <- 0.12
years <- 5
hazard_ratio <- 0.70
rate_dialysis <- 0.12
rate_death <- 0.08
stopifnot(abs(starts / person_years - rate) < 1e-12)
decimals <- function(value, digits) formatC(value, format = "f", digits = digits)
percent <- function(value, digits) paste0(formatC(100 * value, format = "f", digits = digits), "%")
#%%@note 2. Convert the rate over one cycle and over five years
# Survival falls as exp(-rate x time); the probability is one minus that survival.
# สัดส่วนที่ยังไม่เกิดเหตุการณ์ลดลงแบบ exp(-อัตรา x เวลา) ความน่าจะเป็นคือหนึ่งลบสัดส่วนนั้น
still_free <- exp(-rate * cycle_years)
per_thousand <- 1000 * c(rate * cycle_years, 1 - still_free, still_free)
risk_five <- 1 - exp(-rate * years)
rate_back <- -log(1 - risk_five) / years
stopifnot(abs(rate_back - rate) < 1e-12)
yearly <- c(risk_five / years, 1 - exp(-rate_back))
per_cycle <- c(risk_five / 20, 1 - exp(-rate_back * cycle_years))
#%%@note 3. Apply the treatment effect to the rate, then convert
# The trial's ratio multiplies the rate, never the probability; the foundation uses survival 0.60.
# อัตราส่วนจากการทดลองคูณกับอัตรา ไม่ได้คูณกับความน่าจะเป็น ตัวอย่างพื้นฐานใช้การรอดชีพ 0.60
risk_hr <- 1 - 0.60^c(1, hazard_ratio)
treated <- rate * hazard_ratio
five_year <- 1 - exp(-c(rate, treated) * years)
#%%@note 4. Share one survival between two exits
# Add the exit rates, find who stays, then split the leavers by each rate's share.
# รวมอัตราของทางออกทั้งสอง หาสัดส่วนที่อยู่ต่อ แล้วแบ่งผู้ที่ออกตามสัดส่วนของแต่ละอัตรา
exits <- tibble::tibble(exit = c("dialysis", "death"), rate = c(rate_dialysis, rate_death)) |>
dplyr::mutate(separate = 1 - exp(-rate * cycle_years))
rate_total <- sum(exits$rate)
stay <- exp(-rate_total * cycle_years)
row_joint <- c(stay, (1 - stay) * exits$rate / rate_total)
separate <- exits$separate
stopifnot(abs(sum(row_joint) - 1) < 1e-12)
#%%@note 5. Display the model
# Plot the exact risk against rate x time over ten years; no graphics file is saved.
# แสดงความเสี่ยงที่ถูกต้องเทียบกับอัตราคูณเวลาในสิบปี โดยไม่บันทึกไฟล์ภาพ
curve <- tibble::tibble(time = seq(0, 10, by = 0.25)) |>
dplyr::mutate(exact = 1 - exp(-rate * time), shortcut = rate * time)
print(ggplot2::ggplot(curve, ggplot2::aes(x = time)) +
ggplot2::geom_line(ggplot2::aes(y = exact)) +
ggplot2::geom_line(ggplot2::aes(y = shortcut), linetype = "dashed") +
ggplot2::labs(x = "Years", y = "Share who started dialysis"))
#%%@note 6. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Rate per person-year:", rate, " shortcut over 10 years:", sprintf("%.2f", rate * 10), "\n")
cat("Shortcut for one 3-month cycle:", formatC(rate * cycle_years, format = "f", digits = 4), "\n")
cat("Share still free, exp(-0.03):", signif(still_free, 4), "\n")
cat("Cycle probability, 1 minus that share:", signif(1 - still_free, 3), "\n")
cat("Of 1,000 at stage 4, shortcut, starts, still free:", round(per_thousand, 1), "\n")
cat("Rate recovered from the probability:", -log(still_free) / cycle_years, "\n")
cat("Five-year risk, exact and rate x time:", percent(c(risk_five, rate * years), 2), "\n")
cat("Five-year risk and share still free:", signif(c(risk_five, 1 - risk_five), 4), "\n")
cat("Rate recovered from the share:", rate_back, "\n")
cat("One year, divided by 5 and through the rate:", percent(yearly, 2), "\n")
cat("One cycle, divided by 20 and through the rate:", signif(per_cycle, 3), "\n")
cat("Survival 0.60, risks and ratio:", percent(risk_hr, 1), round(risk_hr[2] / risk_hr[1], 2), "\n")
cat("Treated rate, cycle probability:", treated, signif(1 - exp(-treated * cycle_years), 3), "\n")
cat("Five-year risks, control and treated:", percent(five_year, 2), "\n")
cat("Five-year risk ratio:", sprintf("%.3f", five_year[2] / five_year[1]), "\n")
cat("Total rate, each exit's share:", rate_total, c(rate_dialysis, rate_death) / rate_total, "\n")
cat("Stay, dialysis, death:", decimals(row_joint, 4), " row sum:", sum(row_joint), "\n")
cat("Leave in one cycle, joint and separate:", decimals(c(1 - stay, sum(separate)), 4), "\n")
cat("Separate dialysis and death:", decimals(separate, 4), "\n")
การรักษาชะลออัตรา จึงต้องใช้ผลการรักษากับอัตรา
อัตราส่วนฮาซาร์ด (hazard ratio) คือการเปรียบเทียบความเร็วที่เหตุการณ์เกิดในคนที่ยังไม่เกิดเหตุการณ์ทีละขณะ ไม่ใช่อัตราส่วนของความเสี่ยงตลอดช่วงเวลา ถ้ากลุ่มควบคุมรอดชีพ 0.60 คือหกในสิบคนยังมีชีวิตเมื่อครบห้าปี และอัตราส่วนฮาซาร์ดคือ 0.70 ความเสี่ยงห้าปีที่จะเสียชีวิตของกลุ่มควบคุมและกลุ่มที่ได้ยาจะเป็น 40.0% และ 30.1% ตามลำดับ หรือ risk ratio 0.75 อ่านบทเรียนเต็มได้ที่ อัตราส่วนฮาซาร์ดไม่ใช่อัตราส่วนความเสี่ยง: ข้อสมมติความเป็นสัดส่วนของฮาซาร์ด เวลารอดชีพเฉลี่ยแบบจำกัดช่วง และข้อมูลนำเข้าแบบจำลอง ตัวอย่างนี้เป็นคนละกลุ่มและคนละเหตุการณ์กับคลินิก ค่า 0.75 จึงไม่ใช่ตัวเลขของคลินิก ตัวเลขของคลินิกอยู่ถัดไป
การทดลองที่วัดเวลาจนเกิดเหตุการณ์มักรายงานผลแบบนี้ สมมติว่าการทดลองยาชะลอไตเสื่อมตัวหนึ่งรายงานอัตราส่วนฮาซาร์ด 0.70 อัตราส่วนนี้ปรับความเร็ว จึงต้องคูณกับอัตรา 0.12 × 0.70 = 0.084 ครั้งต่อคน-ปีสำหรับผู้ป่วยที่ได้ยา จากนั้นจึงแปลง ในหนึ่งรอบ โอกาสของผู้ที่ได้ยาคือ $1 - e^{-(0.084 \times 0.25)}$ = 0.0208
อย่าเอา 0.70 ไปคูณกับโอกาสหรือความเสี่ยงช่วงยาว ที่อัตรา 0.12 ของคลินิก ความเสี่ยงห้าปีที่ถูกต้องคือ 45.12% เมื่อไม่ได้ยา และ 34.30% เมื่อได้ยา คิดเป็น risk ratio 0.760 ไม่ใช่ 0.70 สองอัตราส่วนนี้ห่างกันเพราะผู้ป่วยที่ยังเสี่ยงอยู่ในกลุ่มที่ไม่ได้ยาลดลงเร็วกว่า ในรอบสั้นหนึ่งรอบทั้งสองแทบเท่ากัน การใช้ค่าหนึ่งแทนอีกค่าจึงพลาดได้ง่ายโดยไม่รู้ตัว
ตรวจตัวเลขของการรักษาด้วย R
# Apply an illustrative hazard ratio to the rate first, then convert; synthetic numbers.
# Writes: nothing to disk; the results print in the console.
rate <- 0.12
hazard_ratio <- 0.70
cycle_years <- 0.25
years <- 5
percent <- function(value, digits) paste0(formatC(100 * value, format = "f", digits = digits), "%")
risk_hr <- 1 - 0.60^c(1, hazard_ratio)
treated <- rate * hazard_ratio
five_year <- 1 - exp(-c(rate, treated) * years)
cat("Survival 0.60, risks and ratio:", percent(risk_hr, 1), round(risk_hr[2] / risk_hr[1], 2), "\n")
cat("Treated rate, cycle probability:", treated, signif(1 - exp(-treated * cycle_years), 3), "\n")
cat("Five-year risks, control and treated:", percent(five_year, 2), "\n")
cat("Five-year risk ratio:", sprintf("%.3f", five_year[2] / five_year[1]), "\n")Survival 0.60, risks and ratio: 40.0% 30.1% 0.75
Treated rate, cycle probability: 0.084 0.0208
Five-year risks, control and treated: 45.12% 34.30%
Five-year risk ratio: 0.760ต้องใช้ R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · ใช้ได้กับ R 4.6.1 · ดาวน์โหลด rates-probabilities-and-transition-probabilities.R
สคริปต์ฉบับเต็ม
# Rates, Probabilities and Transition Probabilities: Keep Time in the Calculation
# Follow one kidney clinic from a counted rate to the chances a 3-month model needs.
# ติดตามคลินิกโรคไตหนึ่งแห่ง จากอัตราที่นับได้ไปสู่ความน่าจะเป็นที่แบบจำลองรอบละ 3 เดือนต้องใช้
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the inputs and the exits as explicit tables.
library(dplyr) # Calculate the converted probabilities.
library(ggplot2) # Display the exact curve against the shortcut without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
starts <- 120
person_years <- 1000
cycle_years <- 0.25
rate <- 0.12
years <- 5
hazard_ratio <- 0.70
rate_dialysis <- 0.12
rate_death <- 0.08
stopifnot(abs(starts / person_years - rate) < 1e-12)
decimals <- function(value, digits) formatC(value, format = "f", digits = digits)
percent <- function(value, digits) paste0(formatC(100 * value, format = "f", digits = digits), "%")
#%%@note 2. Convert the rate over one cycle and over five years
# Survival falls as exp(-rate x time); the probability is one minus that survival.
# สัดส่วนที่ยังไม่เกิดเหตุการณ์ลดลงแบบ exp(-อัตรา x เวลา) ความน่าจะเป็นคือหนึ่งลบสัดส่วนนั้น
still_free <- exp(-rate * cycle_years)
per_thousand <- 1000 * c(rate * cycle_years, 1 - still_free, still_free)
risk_five <- 1 - exp(-rate * years)
rate_back <- -log(1 - risk_five) / years
stopifnot(abs(rate_back - rate) < 1e-12)
yearly <- c(risk_five / years, 1 - exp(-rate_back))
per_cycle <- c(risk_five / 20, 1 - exp(-rate_back * cycle_years))
#%%@note 3. Apply the treatment effect to the rate, then convert
# The trial's ratio multiplies the rate, never the probability; the foundation uses survival 0.60.
# อัตราส่วนจากการทดลองคูณกับอัตรา ไม่ได้คูณกับความน่าจะเป็น ตัวอย่างพื้นฐานใช้การรอดชีพ 0.60
risk_hr <- 1 - 0.60^c(1, hazard_ratio)
treated <- rate * hazard_ratio
five_year <- 1 - exp(-c(rate, treated) * years)
#%%@note 4. Share one survival between two exits
# Add the exit rates, find who stays, then split the leavers by each rate's share.
# รวมอัตราของทางออกทั้งสอง หาสัดส่วนที่อยู่ต่อ แล้วแบ่งผู้ที่ออกตามสัดส่วนของแต่ละอัตรา
exits <- tibble::tibble(exit = c("dialysis", "death"), rate = c(rate_dialysis, rate_death)) |>
dplyr::mutate(separate = 1 - exp(-rate * cycle_years))
rate_total <- sum(exits$rate)
stay <- exp(-rate_total * cycle_years)
row_joint <- c(stay, (1 - stay) * exits$rate / rate_total)
separate <- exits$separate
stopifnot(abs(sum(row_joint) - 1) < 1e-12)
#%%@note 5. Display the model
# Plot the exact risk against rate x time over ten years; no graphics file is saved.
# แสดงความเสี่ยงที่ถูกต้องเทียบกับอัตราคูณเวลาในสิบปี โดยไม่บันทึกไฟล์ภาพ
curve <- tibble::tibble(time = seq(0, 10, by = 0.25)) |>
dplyr::mutate(exact = 1 - exp(-rate * time), shortcut = rate * time)
print(ggplot2::ggplot(curve, ggplot2::aes(x = time)) +
ggplot2::geom_line(ggplot2::aes(y = exact)) +
ggplot2::geom_line(ggplot2::aes(y = shortcut), linetype = "dashed") +
ggplot2::labs(x = "Years", y = "Share who started dialysis"))
#%%@note 6. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Rate per person-year:", rate, " shortcut over 10 years:", sprintf("%.2f", rate * 10), "\n")
cat("Shortcut for one 3-month cycle:", formatC(rate * cycle_years, format = "f", digits = 4), "\n")
cat("Share still free, exp(-0.03):", signif(still_free, 4), "\n")
cat("Cycle probability, 1 minus that share:", signif(1 - still_free, 3), "\n")
cat("Of 1,000 at stage 4, shortcut, starts, still free:", round(per_thousand, 1), "\n")
cat("Rate recovered from the probability:", -log(still_free) / cycle_years, "\n")
cat("Five-year risk, exact and rate x time:", percent(c(risk_five, rate * years), 2), "\n")
cat("Five-year risk and share still free:", signif(c(risk_five, 1 - risk_five), 4), "\n")
cat("Rate recovered from the share:", rate_back, "\n")
cat("One year, divided by 5 and through the rate:", percent(yearly, 2), "\n")
cat("One cycle, divided by 20 and through the rate:", signif(per_cycle, 3), "\n")
cat("Survival 0.60, risks and ratio:", percent(risk_hr, 1), round(risk_hr[2] / risk_hr[1], 2), "\n")
cat("Treated rate, cycle probability:", treated, signif(1 - exp(-treated * cycle_years), 3), "\n")
cat("Five-year risks, control and treated:", percent(five_year, 2), "\n")
cat("Five-year risk ratio:", sprintf("%.3f", five_year[2] / five_year[1]), "\n")
cat("Total rate, each exit's share:", rate_total, c(rate_dialysis, rate_death) / rate_total, "\n")
cat("Stay, dialysis, death:", decimals(row_joint, 4), " row sum:", sum(row_joint), "\n")
cat("Leave in one cycle, joint and separate:", decimals(c(1 - stay, sum(separate)), 4), "\n")
cat("Separate dialysis and death:", decimals(separate, 4), "\n")
สองทางออกใช้ผู้ป่วยกลุ่มเดียวกัน
ผู้ป่วยระยะที่ 4 บางคนเสียชีวิตก่อนที่จะต้องฟอกไต การเสียชีวิตจึงเป็นทางออกที่แข่งขันกัน (competing exit) เมื่อผู้ป่วยออกทางนี้ไปแล้ว การฟอกไตจะมาก่อนไม่ได้อีก ลองนึกถึงห้องรอที่มีสองประตู คือฟอกไตกับเสียชีวิต ทั้งสองประตูรับผู้ป่วยจากกลุ่มเดียวกันในช่วง 3 เดือนเดียวกัน
แบบจำลองต้องบอกได้ว่าผู้ป่วยระยะที่ 4 ทุกคนอยู่ที่ไหนเมื่อถึงนัดครั้งถัดไป คือยังอยู่ระยะที่ 4 ฟอกไตแล้ว หรือเสียชีวิต โอกาสทั้งสามนี้อยู่ด้วยกันในแถวเดียวของตารางในแบบจำลอง และเพราะผู้ป่วยทุกคนไปอยู่ที่เดียวเท่านั้น แถวนี้จึงต้องรวมได้ 1 สำหรับทุกสถานะเริ่มต้น
สมมติว่าอัตราการเสียชีวิตก่อนฟอกไตของคลินิกคือ 0.08 ต่อคน-ปี ถ้าแปลงแต่ละอัตราแยกกัน จะได้ 0.0296 สำหรับการฟอกไต และ 0.0198 สำหรับการเสียชีวิต รวมออก 0.0494 การคำนวณแต่ละครั้งถือว่าประตูของตัวเองเป็นประตูเดียว จึงยังนับผู้ป่วยที่ออกทางประตูอีกบานไปแล้วด้วย
วิธีร่วมจะบวกอัตราทั้งสองก่อน 0.12 + 0.08 = 0.20 ต่อคน-ปี แล้วหาสัดส่วนที่ยังไม่ฟอกไตและยังไม่เสียชีวิต (สัดส่วนที่อยู่ต่อ) $e^{-(0.20 \times 0.25)}$ = 0.9512 จากนั้นแบ่งผู้ที่ออก 0.0488 ตามส่วนของแต่ละอัตราในอัตรารวม คือ 0.12 / 0.20 = 0.60 สำหรับการฟอกไต และ 0.40 สำหรับการเสียชีวิต ได้ 0.0488 × 0.60 = 0.0293 และ 0.0488 × 0.40 = 0.0195 แถวนี้รวมได้พอดี 0.9512 + 0.0293 + 0.0195 = 1
ส่วนการแปลงแยกกันรวมออกได้ 0.0494 ซึ่งมากกว่า 0.0488 ที่ออกจริง
ดังนั้นเมื่อแบบจำลองมีการเสียชีวิตด้วย โอกาสฟอกไตต่อรอบคือ 0.0293 ไม่ใช่ 0.0296 ค่า 0.0296 ใช้ได้เฉพาะเมื่อไม่นับการเสียชีวิต
การแบ่งแบบนี้ใช้ได้เมื่ออัตราทั้งสองคงที่และไม่มีใครย้ายอีกภายในรอบ ถ้าผู้ป่วยออกแล้วกลับมาได้ หรือผ่านสถานะอื่นภายในหนึ่งรอบ ให้ใช้แบบจำลองที่ยอมให้ย้ายแบบนั้น แล้วตรวจว่าทุกแถวยังรวมได้ 1 (Siebert and colleagues; Briggs and colleagues)
ตรวจทางออกทั้งสองด้วย R
# Two exits share the same people: dialysis first or death first; synthetic rates per year.
# Writes: nothing to disk; the results print in the console.
rate_dialysis <- 0.12
rate_death <- 0.08
cycle_years <- 0.25
decimals <- function(value, digits) formatC(value, format = "f", digits = digits)
rate_total <- rate_dialysis + rate_death
stay <- exp(-rate_total * cycle_years)
row_joint <- c(stay, (1 - stay) * c(rate_dialysis, rate_death) / rate_total)
separate <- 1 - exp(-c(rate_dialysis, rate_death) * cycle_years)
cat("Total rate, each exit's share:", rate_total, c(rate_dialysis, rate_death) / rate_total, "\n")
cat("Stay, dialysis, death:", decimals(row_joint, 4), " row sum:", sum(row_joint), "\n")
cat("Leave in one cycle, joint and separate:", decimals(c(1 - stay, sum(separate)), 4), "\n")
cat("Separate dialysis and death:", decimals(separate, 4), "\n")Total rate, each exit's share: 0.2 0.6 0.4
Stay, dialysis, death: 0.9512 0.0293 0.0195 row sum: 1
Leave in one cycle, joint and separate: 0.0488 0.0494
Separate dialysis and death: 0.0296 0.0198ต้องใช้ R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · ใช้ได้กับ R 4.6.1 · ดาวน์โหลด rates-probabilities-and-transition-probabilities.R
สคริปต์ฉบับเต็ม
# Rates, Probabilities and Transition Probabilities: Keep Time in the Calculation
# Follow one kidney clinic from a counted rate to the chances a 3-month model needs.
# ติดตามคลินิกโรคไตหนึ่งแห่ง จากอัตราที่นับได้ไปสู่ความน่าจะเป็นที่แบบจำลองรอบละ 3 เดือนต้องใช้
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the inputs and the exits as explicit tables.
library(dplyr) # Calculate the converted probabilities.
library(ggplot2) # Display the exact curve against the shortcut without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
starts <- 120
person_years <- 1000
cycle_years <- 0.25
rate <- 0.12
years <- 5
hazard_ratio <- 0.70
rate_dialysis <- 0.12
rate_death <- 0.08
stopifnot(abs(starts / person_years - rate) < 1e-12)
decimals <- function(value, digits) formatC(value, format = "f", digits = digits)
percent <- function(value, digits) paste0(formatC(100 * value, format = "f", digits = digits), "%")
#%%@note 2. Convert the rate over one cycle and over five years
# Survival falls as exp(-rate x time); the probability is one minus that survival.
# สัดส่วนที่ยังไม่เกิดเหตุการณ์ลดลงแบบ exp(-อัตรา x เวลา) ความน่าจะเป็นคือหนึ่งลบสัดส่วนนั้น
still_free <- exp(-rate * cycle_years)
per_thousand <- 1000 * c(rate * cycle_years, 1 - still_free, still_free)
risk_five <- 1 - exp(-rate * years)
rate_back <- -log(1 - risk_five) / years
stopifnot(abs(rate_back - rate) < 1e-12)
yearly <- c(risk_five / years, 1 - exp(-rate_back))
per_cycle <- c(risk_five / 20, 1 - exp(-rate_back * cycle_years))
#%%@note 3. Apply the treatment effect to the rate, then convert
# The trial's ratio multiplies the rate, never the probability; the foundation uses survival 0.60.
# อัตราส่วนจากการทดลองคูณกับอัตรา ไม่ได้คูณกับความน่าจะเป็น ตัวอย่างพื้นฐานใช้การรอดชีพ 0.60
risk_hr <- 1 - 0.60^c(1, hazard_ratio)
treated <- rate * hazard_ratio
five_year <- 1 - exp(-c(rate, treated) * years)
#%%@note 4. Share one survival between two exits
# Add the exit rates, find who stays, then split the leavers by each rate's share.
# รวมอัตราของทางออกทั้งสอง หาสัดส่วนที่อยู่ต่อ แล้วแบ่งผู้ที่ออกตามสัดส่วนของแต่ละอัตรา
exits <- tibble::tibble(exit = c("dialysis", "death"), rate = c(rate_dialysis, rate_death)) |>
dplyr::mutate(separate = 1 - exp(-rate * cycle_years))
rate_total <- sum(exits$rate)
stay <- exp(-rate_total * cycle_years)
row_joint <- c(stay, (1 - stay) * exits$rate / rate_total)
separate <- exits$separate
stopifnot(abs(sum(row_joint) - 1) < 1e-12)
#%%@note 5. Display the model
# Plot the exact risk against rate x time over ten years; no graphics file is saved.
# แสดงความเสี่ยงที่ถูกต้องเทียบกับอัตราคูณเวลาในสิบปี โดยไม่บันทึกไฟล์ภาพ
curve <- tibble::tibble(time = seq(0, 10, by = 0.25)) |>
dplyr::mutate(exact = 1 - exp(-rate * time), shortcut = rate * time)
print(ggplot2::ggplot(curve, ggplot2::aes(x = time)) +
ggplot2::geom_line(ggplot2::aes(y = exact)) +
ggplot2::geom_line(ggplot2::aes(y = shortcut), linetype = "dashed") +
ggplot2::labs(x = "Years", y = "Share who started dialysis"))
#%%@note 6. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Rate per person-year:", rate, " shortcut over 10 years:", sprintf("%.2f", rate * 10), "\n")
cat("Shortcut for one 3-month cycle:", formatC(rate * cycle_years, format = "f", digits = 4), "\n")
cat("Share still free, exp(-0.03):", signif(still_free, 4), "\n")
cat("Cycle probability, 1 minus that share:", signif(1 - still_free, 3), "\n")
cat("Of 1,000 at stage 4, shortcut, starts, still free:", round(per_thousand, 1), "\n")
cat("Rate recovered from the probability:", -log(still_free) / cycle_years, "\n")
cat("Five-year risk, exact and rate x time:", percent(c(risk_five, rate * years), 2), "\n")
cat("Five-year risk and share still free:", signif(c(risk_five, 1 - risk_five), 4), "\n")
cat("Rate recovered from the share:", rate_back, "\n")
cat("One year, divided by 5 and through the rate:", percent(yearly, 2), "\n")
cat("One cycle, divided by 20 and through the rate:", signif(per_cycle, 3), "\n")
cat("Survival 0.60, risks and ratio:", percent(risk_hr, 1), round(risk_hr[2] / risk_hr[1], 2), "\n")
cat("Treated rate, cycle probability:", treated, signif(1 - exp(-treated * cycle_years), 3), "\n")
cat("Five-year risks, control and treated:", percent(five_year, 2), "\n")
cat("Five-year risk ratio:", sprintf("%.3f", five_year[2] / five_year[1]), "\n")
cat("Total rate, each exit's share:", rate_total, c(rate_dialysis, rate_death) / rate_total, "\n")
cat("Stay, dialysis, death:", decimals(row_joint, 4), " row sum:", sum(row_joint), "\n")
cat("Leave in one cycle, joint and separate:", decimals(c(1 - stay, sum(separate)), 4), "\n")
cat("Separate dialysis and death:", decimals(separate, 4), "\n")
ตรวจทั้งกระบวนการ ไม่ใช่แค่สูตร
ก่อนนำตัวเลขที่แปลงแล้วเข้าแบบจำลอง ให้เขียนสถานะเริ่มต้น สถานะปลายทาง ความยาวรอบ และทางออกอื่นไว้ แปลงย้อนกลับแล้วตรวจว่าได้อัตราเดิมคืนมา เหมือนที่ 0.12 ของคลินิกกลับมา จากนั้นตรวจว่าทุกแถวของโอกาสรวมได้ 1 ก่อนเริ่มรันแบบจำลอง
การตรวจเหล่านี้จับความผิดพลาดของเลขคณิตได้ แต่จับข้อสมมติที่ผิดไม่ได้ อัตราคงที่เป็นคำกล่าวเกี่ยวกับตัวโรค ถ้าการเริ่มฟอกไตเร็วขึ้นในปีหลัง ๆ อัตรา 0.12 ค่าเดียวก็ผิด ไม่ว่าจะแปลงอย่างระวังเพียงใด อัตราส่วนฮาซาร์ดที่ไม่คงที่ที่ 0.70 ตลอดการติดตามก็เช่นกัน รอบที่สั้นลงลดความคลาดเคลื่อนของทางลัดได้ แต่แก้นิยามเหตุการณ์ที่ผิดหรือประชากรที่ไม่ตรงไม่ได้
คำศัพท์ที่ควรแยกให้ชัด
- Rate (อัตรา)
- จำนวนเหตุการณ์ต่อหน่วยเวลาที่เสี่ยง เช่น เริ่มฟอกไต 0.12 ครั้งต่อคน-ปี และไม่มีเพดาน
- Hazard (ฮาซาร์ด)
- อัตรา ณ ขณะหนึ่งในกลุ่มคนที่ยังไม่เกิดเหตุการณ์
- Probability (risk) (ความน่าจะเป็น (ความเสี่ยง))
- สัดส่วนของกลุ่มที่เกิดเหตุการณ์แล้วเมื่อสิ้นช่วงเวลาที่กำหนด และไม่เกิน 1 สำหรับทุกกลุ่ม
- Cycle (รอบ)
- ช่วงเวลาคงที่ที่แบบจำลองย้ายผู้ป่วยระหว่างสถานะ ในที่นี้คือ 0.25 ปี
- Transition probability (ความน่าจะเป็นของการเปลี่ยนสถานะ)
- โอกาสที่จะย้ายจากสถานะหนึ่งไปอีกสถานะหนึ่งภายในหนึ่งรอบ ซึ่งแปลงมาจากอัตรา ในที่นี้คือ 0.0296
- Hazard ratio (อัตราส่วนฮาซาร์ด)
- อัตราของกลุ่มที่ได้รักษาหารด้วยอัตราของกลุ่มที่ไม่ได้รักษา ให้คูณกับอัตราก่อน แล้วจึงแปลง
- Competing exit (ทางออกที่แข่งขันกัน)
- เหตุการณ์อื่น เช่น การเสียชีวิต ที่พาผู้ป่วยออกจากสถานะเริ่มต้นก่อน
เชื่อมกับบทความพื้นฐาน
บทความพื้นฐานเรื่องอัตราส่วนฮาซาร์ดและ proportional hazards แสดงว่าอัตราส่วน 0.70 กลายเป็น risk ratio 0.75 ในห้าปีได้อย่างไร ส่วนบทความเรื่อง risk ratio และผลต่างความเสี่ยง อธิบายฝั่งความเสี่ยง บทความตอนหลังของชุดนี้นำโอกาสเหล่านี้ไปใส่แบบจำลองกลุ่มประชากร ทีละรอบ
ลองใช้ในแอป
เปิดแอปแบบจำลองกลุ่มประชากรของเราเพื่อลองเปลี่ยนค่าข้อมูลนำเข้าและดูผลลัพธ์เปลี่ยนตาม ต้องเข้าสู่ระบบและใช้แผน Pro หรือ Max
แอปเริ่มจากค่าตัวอย่างของแอปเอง ตัวเลขในแอปจึงไม่ตรงกับบทความนี้
เอกสารอ้างอิง
- Fleurence Rachael L, Hollenbeak Christopher S. Rates and Probabilities in Economic Modelling. PharmacoEconomics. 2007;25(1):3-6. doi:10.2165/00019053-200725010-00002. https://doi.org/10.2165/00019053-200725010-00002
- Siebert Uwe, Alagoz Oguzhan, Bayoumi Ahmed M., Jahn Beate, Owens Douglas K., Cohen David J., Kuntz Karen M.. State-Transition Modeling: A Report of the ISPOR-SMDM Modeling Good Research Practices Task Force-3. Value in Health. 2012;15(6):812-820. doi:10.1016/j.jval.2012.06.014. https://doi.org/10.1016/j.jval.2012.06.014
- Andrew Briggs, Karl Claxton, Mark Sculpher. Decision Modelling for Health Economic Evaluation (Handbooks for Health Economic Evaluation). Oxford University Press; 2006. ISBN 9780198526629. https://openlibrary.org/isbn/9780198526629
ลำดับการเรียนในชุดนี้
- คู่มือชุดบทความการสร้างแบบจำลองเศรษฐศาสตร์สุขภาพ: จากหลักฐานสู่การตัดสินใจ
- การเลือกแบบจำลองเพื่อการตัดสินใจ: ให้เส้นทางทางคลินิกกำหนดโครงสร้าง
- ข้อมูลนำเข้าของแบบจำลองที่ตรวจสอบได้: สร้างตารางพารามิเตอร์ที่สืบย้อนถึงที่มา
- การค้นหาข้อมูลนำเข้าของแบบจำลองจากวรรณกรรม: เลือกค่าประมาณให้ตรงกับคำถาม
- อัตรา ความน่าจะเป็น และความน่าจะเป็นของการเปลี่ยนสถานะ: คำนวณโดยไม่หลงหน่วยเวลา (บทความนี้)
- สร้างแผนภูมิต้นไม้การตัดสินใจทีละขั้น: ตามรอยผลลัพธ์ ต้นทุน และสุขภาพ
- สร้างแบบจำลองมาร์คอฟแบบกลุ่มประชากรทีละขั้น: ติดตามประชากรผ่านเวลา
- การวิเคราะห์ผลกระทบด้านงบประมาณในทางปฏิบัติ: นับคนให้ชัดก่อนคำนวณเงิน
- ICER และผลประโยชน์สุทธิในรูปตัวเงินเมื่อมีหลายทางเลือก: หาขอบเขตทางเลือกที่มีประสิทธิภาพ
- การวิเคราะห์ความไวแบบอาศัยความน่าจะเป็น: จากการแจกแจงพารามิเตอร์สู่ CEAC
- การประเมินความสมเหตุสมผลและตรวจสอบแบบจำลอง: ทดสอบคำถาม สมการ และโค้ด
- มูลค่าคาดหมายของข้อมูลสมบูรณ์: หลักฐานที่ดีขึ้นอาจมีมูลค่าเพียงใด
ประเด็นสำคัญ
- อัตราคือความเร็ว 0.12 ครั้งต่อคน-ปี ส่วนความน่าจะเป็นคือสัดส่วน ซึ่งไม่เกิน 1 สำหรับทุกกลุ่ม
- ให้เวลาอยู่ในการแปลง: $1 - e^{-(0.12 \times 0.25)}$ = 0.0296 ต่อรอบ 3 เดือนเมื่อการฟอกไตเป็นทางออกเดียว ไม่ใช่ทางลัด 0.0300
- แปลงความเสี่ยงช่วงยาวเป็นช่วงเวลาใหม่โดยผ่านอัตรา: 45.12% ใน 5 ปี ให้ 0.0296 ต่อรอบ ไม่ใช่ 0.0226 จากการหาร
- ใช้อัตราส่วนฮาซาร์ดกับอัตราก่อนแล้วจึงแปลง: 0.12 × 0.70 = 0.084 ให้ 0.0208 ต่อรอบ และในช่วง 5 ปี risk ratio คือ 0.760 ไม่ใช่ 0.70
- ให้ทางออกที่แข่งขันกันใช้ผู้ป่วยกลุ่มเดียวกัน แถวจึงรวมได้ 1: อยู่ต่อ 0.9512 เริ่มฟอกไต 0.0293 และเสียชีวิต 0.0195
อ่านต่อในวิกิ: อัตราส่วนฮาซาร์ดไม่ใช่อัตราส่วนความเสี่ยง: ข้อสมมติความเป็นสัดส่วนของฮาซาร์ด เวลารอดชีพเฉลี่ยแบบจำกัดช่วง และข้อมูลนำเข้าแบบจำลอง · อัตราส่วนความเสี่ยง อัตราส่วนออดส์ ผลต่างความเสี่ยง และจำนวนที่ต้องรักษา: ตารางสองคูณสองเดียว คำตอบสี่แบบ