Link และ Family: ตัวหนึ่งกำหนดมาตรวัดผล อีกตัวกำหนดค่าคลาดเคลื่อนมาตรฐาน
On this page
Read the English version
บทคัดย่อ
แบบจำลองเชิงเส้นวางนัยทั่วไป (generalized linear model, GLM) มีสามส่วน ตัวทำนายเชิงเส้น (linear predictor) หรือ eta รวมตัวแปรร่วมเข้าด้วยกัน ฟังก์ชันเชื่อมโยง (link function) เชื่อมค่าเฉลี่ยของผลลัพธ์กับ eta และ family ให้ฟังก์ชันความแปรปรวน (variance function) link กำหนดมาตรวัดผล คือ identity ให้ผลต่างความเสี่ยง (risk difference) log ให้ risk ratio และ logit ให้ odds ratio ส่วน family กำหนดค่าคลาดเคลื่อนมาตรฐาน และเมื่อมีตัวแปรร่วม กำหนดว่าผู้ป่วยแต่ละคนมีน้ำหนักเท่าใด ภายใต้ log link ความเสี่ยงเท่ากับเลขชี้กำลังของ eta จึงต้องให้ eta ไม่เกินศูนย์ และใกล้กำแพงนั้นการประมาณแบบ log-binomial จะล้มเหลว Poisson ดัดแปลง (modified Poisson) คง log link ไว้และไม่ชนกำแพง ในทะเบียนจำลองของผู้ป่วยกระดูกสะโพกหัก 20,000 คน แบบจำลองทั้งสองให้ผลตรงกันทุกประการเมื่อมีการผ่าตัดเร็วเพียงตัวเดียว ต่างกันเล็กน้อยหลังปรับด้วยอายุและเพศ และ log-binomial ล้มเหลวเมื่อใส่อายุและความเปราะบาง โดยรวม link ควรตามมาตรวัดผล family ที่ยืมมาจาก Poisson ต้องใช้ค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช ซึ่งประมาณจากส่วนเหลือที่สังเกตได้ และความต่างเล็กน้อยหลังปรับสะท้อนการถ่วงน้ำหนัก ไม่ใช่ความผิดพลาด
นักวิจัยสองคน ทะเบียนเดียว มาตรวัดผลสองแบบ
ทีมทะเบียนผู้ป่วยประชุมทบทวนการวิเคราะห์ภาวะสับสนเฉียบพลัน (delirium) หลังผ่าตัดกระดูกสะโพกหัก ภาวะนี้คือสภาวะสับสนที่เกิดขึ้นเฉียบพลัน และเกิดใน 31.0% ของผู้สูงอายุ 20,000 คนในทะเบียนจำลองนี้ ตัวแปรสัมผัสคือการผ่าตัดเร็ว ภายในวันแรกหลังรับไว้ในโรงพยาบาล
นักวิจัยคนแรกใช้ logistic regression และรายงาน odds ratio (อัตราส่วนออดส์) นักวิจัยคนที่สองรายงาน risk ratio (อัตราส่วนความเสี่ยง) จากแบบจำลองที่ใช้ log link คือแบบจำลองบนลอการิทึมของความเสี่ยง เมื่อใช้ family (การแจกแจงที่สมมติให้ผลลัพธ์) แบบ binomial แล้วเปลี่ยนเป็นแบบ Poisson ค่า risk ratio ของเธอตรงกันถึงสี่ตำแหน่งทศนิยม แต่ช่วงเชื่อมั่นไม่ตรงกัน
ผู้เชี่ยวชาญด้านระเบียบวิธีที่เป็นประธานการประชุมถามว่าการเลือกแต่ละอย่างทำหน้าที่อะไร สรุปสั้น ๆ คือ link เลือกมาตรวัดผล ส่วน family เลือกค่าคลาดเคลื่อนมาตรฐาน
สามส่วนของแบบจำลองเชิงเส้นวางนัยทั่วไป
แบบจำลองเชิงเส้นวางนัยทั่วไป (GLM) มีสามส่วน [1, 2] ส่วนแรกคือตัวทำนายเชิงเส้น (linear predictor) $\eta$ (อักษรกรีก eta) ซึ่งเป็นผลรวมถ่วงน้ำหนักของตัวแปรร่วมของผู้ป่วย:
$$\eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots$$แต่ละ $x$ คือค่าของตัวแปรร่วม เช่นตัวบ่งชี้การผ่าตัดเร็ว และแต่ละ $\beta$ คือสัมประสิทธิ์ที่ต้องประมาณ ผลรวมนี้เขียนย่อเป็น $X\beta$ และมีค่าได้ทั้งลบและบวก
ส่วนที่สองคือฟังก์ชันเชื่อมโยง (link function) $g(\cdot)$ ซึ่งเชื่อมค่าเฉลี่ยของผลลัพธ์กับตัวทำนายเชิงเส้น:
$$g(\mu) = \eta$$ในที่นี้ $\mu$ คือค่าเฉลี่ยของผลลัพธ์ $Y$ ในผู้ป่วยที่มีค่าตัวแปรร่วมเหมือนกัน สำหรับผลลัพธ์ที่ลงรหัส 1 เมื่อเกิดภาวะสับสนเฉียบพลันและ 0 เมื่อไม่เกิด $\mu$ คือความเสี่ยง ฟังก์ชันผกผันของ link เปลี่ยนตัวทำนายเชิงเส้นกลับเป็นค่าเฉลี่ย
ส่วนที่สามคือ family ซึ่งเป็นการแจกแจงที่สมมติให้ผลลัพธ์กระจายรอบค่าเฉลี่ย ในการประมาณ หน้าที่หลักของมันคือให้ฟังก์ชันความแปรปรวน (variance function) $V(\mu)$ คือความแปรปรวนของผลลัพธ์ขึ้นกับค่าเฉลี่ยอย่างไร
ผลลัพธ์ family และ link: การจับคู่หลักของ GLM
| ผลลัพธ์ | Family, $V(\mu)$ | Link, $g(\mu)$ | สัมประสิทธิ์ของตัวแปรสัมผัส |
|---|---|---|---|
| ต่อเนื่อง เช่นความดันโลหิต | Gaussian (แจกแจงปกติ) ค่าคงที่ | Identity, $\mu$ | ผลต่างของค่าเฉลี่ย |
| สองค่า (0 หรือ 1) | Binomial, $\mu(1-\mu)$ | Identity, $\mu$ | ผลต่างความเสี่ยง |
| สองค่า (0 หรือ 1) | Binomial, $\mu(1-\mu)$ | Log, $\log \mu$ | Log risk ratio |
| สองค่า (0 หรือ 1) | Binomial, $\mu(1-\mu)$ | Logit, $\log\{\mu/(1-\mu)\}$ | Log odds ratio |
| สองค่า (0 หรือ 1) แบบ modified Poisson | Poisson เป็น working family, $\mu$, พร้อมค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช | Log, $\log \mu$ | Log risk ratio |
| จำนวนนับของเหตุการณ์ | Poisson, $\mu$ | Log, $\log \mu$ | Log rate ratio โดยใช้ลอการิทึมของเวลาติดตามเป็นออฟเซต (offset) คือพจน์ที่สัมประสิทธิ์ถูกตรึงไว้ที่ 1 |
link กำหนดมาตรวัดผลอย่างไร
ให้ตัวแปรสัมผัสแบบสองค่า $A$ ลงรหัส 1 สำหรับการผ่าตัดเร็วและ 0 สำหรับกรณีอื่น และให้ $\mu_1$ กับ $\mu_0$ เป็นความเสี่ยงเมื่อมีและไม่มีการผ่าตัดเร็ว เมื่อใช้ identity link $g(\mu) = \mu$ แบบจำลองคือ $\mu = \beta_0 + \beta_1 A$ ดังนั้น $\beta_1$ คือผลต่างความเสี่ยง (risk difference) $\mu_1 - \mu_0$ เมื่อใช้ binomial family link นี้ต้องรักษาความเสี่ยงที่แบบจำลองพอดีทุกค่าให้อยู่ระหว่าง 0 กับ 1 [3]
ใน Stata ตั้ง link ภายใน glm เช่น glm delirium surg24, family(binomial) link(identity) ส่วนใน R คือ glm(delirium ~ surg24, family = binomial(link = "identity"), data = d) ในที่นี้ surg24 คือตัวแปรสัมผัส $A$ และ d คือข้อมูลทะเบียนที่ใช้ต่อไปนี้
Log link (ln): อัตราส่วนจากสมบัติของลอการิทึม
เมื่อใช้ log link แบบจำลองคือ $\log \mu = \beta_0 + \beta_1 A$ โดย $\log$ คือลอการิทึมธรรมชาติ (ln) เมื่อนำกลุ่มที่ไม่สัมผัสไปลบจากกลุ่มที่สัมผัสจะได้
$$\log \mu_1 - \log \mu_0 = \log \frac{\mu_1}{\mu_0} = \beta_1$$ผลต่างของลอการิทึมคือลอการิทึมของอัตราส่วน ดังนั้น $\beta_1$ คือ log risk ratio และเอกซ์โปเนนเชียลของมัน $e^{\beta_1}$ คือ risk ratio เอง ส่วนข้อมูลนับ พีชคณิตเดียวกันให้ rate ratio [2]
Logit link: log odds
logit link ใช้ลอการิทึมของออดส์ (odds) คือความเสี่ยงหารด้วยหนึ่งลบความเสี่ยง:
$$\operatorname{logit}(\mu) = \log \frac{\mu}{1-\mu} = \beta_0 + \beta_1 A$$การลบแบบเดียวกันทำให้ $e^{\beta_1}$ เป็น odds ratio ฟังก์ชันผกผัน $\mu = e^{\eta}/(1 + e^{\eta})$ อยู่ระหว่าง 0 กับ 1 อย่างเคร่งครัดสำหรับทุกค่าของ $\eta$ logit link จึงไม่เคยเรียกร้องความเสี่ยงที่สูงกว่า 1 และสำหรับความเสี่ยงสองค่าเดียวกัน odds ratio อยู่ห่างจาก 1 มากกว่า risk ratio เสมอเมื่อความเสี่ยงสองค่าต่างกัน และห่างมากกว่ามากเมื่อผลลัพธ์พบบ่อย
ตัวอย่างคำนวณด้วยมือ: ความเสี่ยงหนึ่งคู่ สาม link
ตัวอย่างคำนวณด้วยมือ สมมติว่าผู้ป่วย 40 คนจาก 100 คนที่สัมผัส และ 20 คนจาก 100 คนที่ไม่สัมผัส เกิดภาวะสับสนเฉียบพลัน
-
ความเสี่ยงสองค่า
\[ \mu_1 = \frac{40}{100} = 0.400, \quad \mu_0 = \frac{20}{100} = 0.200 \]
$\mu_1$ คือความเสี่ยงในกลุ่มที่สัมผัส และ $\mu_0$ คือความเสี่ยงในกลุ่มที่ไม่สัมผัส
-
Identity link: ผลต่างความเสี่ยง
\[ \beta_1 = \mu_1 - \mu_0 = 0.400 - 0.200 = 0.200 \]
สัมประสิทธิ์คือผลต่างนั้นเอง
-
Log link: risk ratio
\[ e^{\beta_1} = \frac{0.400}{0.200} = 2.00, \quad \beta_1 = \log 2.00 = 0.69 \]
สัมประสิทธิ์คือ 0.69 และเอกซ์โปเนนเชียลของมันคือ risk ratio
-
ออดส์ในแต่ละกลุ่ม
\[ \frac{0.400}{0.600} = 0.667, \quad \frac{0.200}{0.800} = 0.250 \]
ออดส์แต่ละค่าคือความเสี่ยงหารด้วยหนึ่งลบความเสี่ยง
-
Logit link: odds ratio
\[ e^{\beta_1} = \frac{0.667}{0.250} = 2.67, \quad \beta_1 = \log 2.67 = 0.98 \]
สัมประสิทธิ์คือ 0.98 และเอกซ์โปเนนเชียลของมันคือ odds ratio
ผลลัพธ์: ความเสี่ยงหนึ่งคู่ให้ผลต่างความเสี่ยง 0.200 risk ratio 2.00 และ odds ratio 2.67 โดย link เป็นตัวตัดสินว่าสัมประสิทธิ์รายงานค่าใด
family: จุดที่การแจกแจงเข้ามามีบทบาท
link และตัวทำนายเชิงเส้นกำหนดค่าเฉลี่ย ส่วน family อธิบายว่าผลลัพธ์ที่สังเกตได้กระจายรอบค่าเฉลี่ยนั้นอย่างไร เพราะผู้ป่วยสองคนที่มีตัวแปรร่วมเหมือนกันทุกประการก็ยังต่างกันได้ family เข้าสู่การประมาณผ่านฟังก์ชันความแปรปรวน [2]:
- binomial: $V(\mu) = \mu(1-\mu)$ ซึ่งเป็นความแปรปรวนของผลลัพธ์ที่ลงรหัส 0 หรือ 1
- Poisson: $V(\mu) = \mu$ ซึ่งเป็นความแปรปรวนของจำนวนนับ
- Gaussian (แจกแจงปกติ): $V(\mu)$ คงที่ คือความกระจายเท่ากันที่ทุกค่าเฉลี่ย
สำหรับผลลัพธ์แบบสองค่า มีเพียงความแปรปรวนแบบ binomial ที่ถูกต้อง ความแปรปรวนแบบ Poisson $\mu$ มากกว่า $\mu(1-\mu)$ ที่ทุกระดับความเสี่ยง และช่องว่างยิ่งกว้างเมื่อความเสี่ยงสูงขึ้น ค่าคลาดเคลื่อนมาตรฐานแบบอิงแบบจำลอง (model-based) สร้างจาก $V(\mu)$ ดังนั้น family ที่ไม่ตรงกับผลลัพธ์จะพิมพ์ค่าคลาดเคลื่อนมาตรฐานที่ผิด แม้สัมประสิทธิ์จะถูกต้อง
ทำไมต้องมี family: น้ำหนักและค่าคลาดเคลื่อนมาตรฐาน
GLM หาสัมประสิทธิ์โดยแก้สมการสกอร์ (score equations) หนึ่งสมการต่อหนึ่งสัมประสิทธิ์ แต่ละสมการตั้งผลรวมถ่วงน้ำหนักของส่วนเหลือ (residual) ให้เท่ากับศูนย์ [2] เมื่อใช้ log link สมการเหล่านี้คือ
$$\sum_{i=1}^{n} x_i \, \frac{\mu_i}{V(\mu_i)} \, (y_i - \mu_i) = 0$$ในที่นี้ $y_i$ คือผลลัพธ์ของผู้ป่วยคนที่ $i$ และ $\mu_i$ คือความเสี่ยงที่แบบจำลองพอดี $x_i$ เก็บค่าตัวแปรร่วมของผู้ป่วย โดยมี 1 สำหรับจุดตัดแกน และ $n$ คือจำนวนผู้ป่วย ตัวคูณ $\mu_i/V(\mu_i)$ คือน้ำหนักของส่วนเหลือแต่ละตัว
เมื่อใช้ Poisson family น้ำหนักนั้นเท่ากับ 1 ส่วนเหลือทุกตัวจึงนับเท่ากัน เมื่อใช้ binomial family น้ำหนักคือ $1/(1-\mu_i)$ ผู้ป่วยที่มีความเสี่ยงที่แบบจำลองพอดีสูงจึงดึงแรงกว่า
สมการสกอร์ของ Poisson มีค่าคาดหมายเป็นศูนย์ที่สัมประสิทธิ์จริง ตราบที่แบบจำลองค่าเฉลี่ย $\log \mu_i = X_i\beta$ ถูกต้อง ไม่ว่าความแปรปรวนจริงจะเป็นอย่างไร ค่าที่ได้จากการแก้สมการเหล่านี้จึงมีความคงเส้นคงวา (consistent) สำหรับ log risk ratio คือเข้าใกล้ค่าจริงเมื่อตัวอย่างใหญ่ขึ้น [4] มีเพียงค่าคลาดเคลื่อนมาตรฐานของ Poisson ซึ่งสร้างจาก $V(\mu) = \mu$ เท่านั้นที่ผิด
ค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช (sandwich standard error) หรือเรียกว่า robust standard error มาแทนที่ค่านั้นโดยประมาณความแปรปรวนจากส่วนเหลือที่สังเกตได้ แทนที่จะใช้ $V(\mu)$ แซนด์วิชซ่อมอะไรได้และซ่อมอะไรไม่ได้ อธิบายในตอนที่ 4 ของชุดบทความนี้
กำแพง: ทำไม log link จึงเรียกร้องความเสี่ยงที่สูงกว่า 1 ได้
ภายใต้ log link ฟังก์ชันผกผันคือ $\mu = e^{\eta}$ ซึ่งเป็นบวกสำหรับทุกค่าของ $\eta$ แต่ไม่มีเพดาน เมื่อ $\eta = 0$ ค่านี้เท่ากับ 1 และทุก $\eta$ ที่เป็นบวกให้ค่าเฉลี่ยที่สูงกว่า 1 ความเสี่ยงสูงกว่า 1 ไม่ได้ ดังนั้นแบบจำลอง log-binomial ซึ่งคือ binomial family ที่ใช้ log link ต้องรักษาตัวทำนายเชิงเส้นของผู้ป่วยทุกคนให้ไม่เกิน 0:
$$\eta_i = X_i\beta \le 0 \quad \text{for every patient } i$$อสมการนี้คือกำแพง ซึ่งเป็นขอบเขตจำกัด (boundary) ของสัมประสิทธิ์ โดยมีหนึ่งข้อบังคับต่อผู้ป่วยหนึ่งคน [3, 5]
ปัญหาเริ่มเมื่อค่าที่พอดีที่สุดอยู่บนหรือใกล้กำแพง เรื่องนี้เกิดเมื่อรูปแบบตัวแปรร่วมบางแบบ (covariate pattern) คือชุดค่าตัวแปรร่วมหนึ่งชุด มีความเสี่ยงใกล้ 1 หรือเมื่อแบบจำลองดันความเสี่ยงให้เกิน 1 ขั้นตอนคำนวณเสนอก้าวที่ข้ามกำแพงซ้ำแล้วซ้ำอีก แล้วต้องตัดก้าวให้สั้นลง ผลคือไม่ถึงการลู่เข้า (convergence) ซึ่งเป็นจุดที่ค่าประมาณต่อเนื่องหยุดเปลี่ยนแปลง หรือหยุดบนขอบเขตพร้อมคำเตือน [5]
ตัวอย่างคำนวณด้วยมือ: รูปแบบตัวแปรร่วมที่อยู่เลยกำแพง
ตัวอย่างคำนวณด้วยมือ แบบจำลอง log-binomial สมมติ risk ratio ค่าเดียวสำหรับทุกรูปแบบตัวแปรร่วม สมมติว่ารูปแบบหนึ่งมีความเสี่ยงพื้นฐาน 0.60 และ risk ratio ของตัวแปรสัมผัสในแบบจำลองคือ 2.00
-
ความเสี่ยงที่แบบจำลองต้องการ
\[ \mu = 0.60 \times 2.00 = 1.20 \]
ผู้ป่วยที่สัมผัสในรูปแบบนั้นจะต้องมีความเสี่ยง 1.20
-
ตัวทำนายเชิงเส้นภายใต้ log link
\[ \eta = \log 1.20 = 0.18 \]
ค่าเฉลี่ยนั้นสอดคล้องกับตัวทำนายเชิงเส้น 0.18
-
อยู่ฝั่งไหนของกำแพง
\[ \eta = 0.18 > 0 \;\Rightarrow\; \mu = e^{\eta} > 1 \]
ตัวทำนายเชิงเส้นที่เป็นบวกให้ค่าเฉลี่ยสูงกว่า 1 ซึ่งไม่ใช่ความเสี่ยง
ผลลัพธ์: แบบจำลอง log-binomial แทนรูปแบบนี้ไม่ได้ การประมาณจึงติดค้างอยู่ที่กำแพงหรือหยุดบนขอบเขต
modified Poisson ผ่านกำแพงได้อย่างไร
Poisson ดัดแปลง (modified Poisson regression) คง log link ไว้ แต่เปลี่ยน binomial family เป็น Poisson family โดยใช้ Poisson เป็น working family คือให้สมการสกอร์โดยไม่ถือว่าเป็นการแจกแจงของผลลัพธ์ [4] ค่าเฉลี่ยของ Poisson คือจำนวนนับที่คาดหมาย ค่า $e^{\eta}$ ที่สูงกว่า 1 จึงยอมให้มีได้ และการประมาณไม่มีกำแพง สมการค่าเฉลี่ยไม่เปลี่ยน ดังนั้น $e^{\beta_1}$ ยังคงเป็น risk ratio และรายงานพร้อมค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช
ค่าที่แบบจำลองพอดีซึ่งสูงกว่า 1 ไม่ใช่ความเสี่ยง มันแสดงว่าแบบจำลองค่าเฉลี่ยเชิงเส้นบนสเกลลอการิทึมไม่อาจถูกต้องสมบูรณ์สำหรับผู้ป่วยเหล่านั้น Poisson ดัดแปลงเป็นวิธีที่ใช้ได้จริงในการประมาณ risk ratio ไม่ใช่แบบจำลองที่ใกล้ความจริงกว่าสำหรับผลลัพธ์แบบสองค่า [6]
ทะเบียนจำลอง: ตัวแปรร่วมตัวเดียว risk ratio เท่ากัน
ในทะเบียนจำลองที่ใช้ต่อจากนี้ ข้อมูลถูกสร้างขึ้นให้ผู้ป่วยที่เปราะบาง และผู้ที่มีภาวะสมองเสื่อมหรือใช้ยาต้านการแข็งตัวของเลือด รอผ่าตัดนานกว่า ทุกการเปรียบเทียบด้านล่างจึงมีตัวกวน (confounding) การประมาณเหล่านี้แสดงว่า link และ family ทำงานอย่างไร ไม่ได้ประมาณผลของการผ่าตัดเร็ว
การประมาณชุดแรกมีเพียงการผ่าตัดเร็ว แบบจำลองจึงอิ่มตัว (saturated) คือมีหนึ่งพารามิเตอร์ต่อหนึ่งกลุ่ม และทำซ้ำความเสี่ยงที่สังเกตได้ทั้งสองค่า ผู้ป่วยทุกคนในกลุ่มเดียวกันมีความเสี่ยงที่แบบจำลองพอดีและน้ำหนักเท่ากัน สมการสกอร์จึงตั้งความเสี่ยงที่แบบจำลองพอดีของแต่ละกลุ่มให้เท่ากับความเสี่ยงที่สังเกตได้ ไม่ว่าน้ำหนักจะเป็นเท่าใด และ family ทั้งสองให้ risk ratio เท่ากัน
การผ่าตัดเร็วเพียงตัวเดียว: สอง family หนึ่ง risk ratio
| แบบจำลอง | Risk ratio | ช่วงเชื่อมั่น 95% | ค่าคลาดเคลื่อนมาตรฐานของ log risk ratio |
|---|---|---|---|
| Log-binomial ค่าคลาดเคลื่อนมาตรฐานแบบอิงแบบจำลอง | 0.3909 | 0.3709 ถึง 0.4120 | 0.0268 |
| Poisson ดัดแปลง ค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช | 0.3909 | 0.3709 ถึง 0.4120 | 0.0268 |
| Poisson ค่าคลาดเคลื่อนมาตรฐานแบบอิงแบบจำลอง | 0.3909 | 0.3683 ถึง 0.4149 | 0.0304 |
ทำไมช่วงเชื่อมั่นแบบอิงแบบจำลองของ Poisson จึงกว้างเกินไป
ทั้งสามแถวมี risk ratio เดียวกันคือ 0.3909 ค่าคลาดเคลื่อนมาตรฐานของ log risk ratio แบบอิงแบบจำลองของ Poisson คือ 0.0304 มากกว่าค่าแบบแซนด์วิช 0.0268 เพราะ $V(\mu) = \mu$ ประเมินความแปรปรวนของผลลัพธ์แบบสองค่าสูงเกินจริง ช่วงเชื่อมั่นของมัน 0.3683 ถึง 0.4149 จึงกว้างเกินไป ค่าแบบแซนด์วิชเท่ากับค่าคลาดเคลื่อนมาตรฐานของ log-binomial ถึงทศนิยมสี่ตำแหน่ง ซึ่งฟังก์ชันความแปรปรวนของ log-binomial ถูกต้องในกรณีนี้
Stata: การผ่าตัดเร็วเพียงตัวเดียวภายใต้สอง family
glm delirium surg24, family(binomial) link(log)
* model-based CI (identical bread in both languages for a saturated model)
ciz c2.rr_crude_logbin surg24 exp
* log-RR SE, log-binomial, model-based
canon c2.se_crude_logbin _se[surg24]
quietly glm delirium surg24, family(poisson) link(log)
scalar sc_pn = _se[surg24]
glm delirium surg24, family(poisson) link(log) vce(robust)
matrix bpc = e(b)
* modified Poisson, robust CI
ciz c2.rr_crude_poisson surg24 exp
* Poisson model-based log-RR SE (too large for a binary outcome) versus the sandwich log-RR SE
canon c2.se_crude_poisson_naive scalar(sc_pn)
canon c2.se_crude_poisson_robust _se[surg24]
* naive Poisson 95% CI
canon c2.rr_crude_poisson_naive.lo exp(_b[surg24] - invnormal(0.975)*scalar(sc_pn))
canon c2.rr_crude_poisson_naive.hi exp(_b[surg24] + invnormal(0.975)*scalar(sc_pn)). glm delirium surg24, family(binomial) link(log)
Iteration 0: Log likelihood = -15980.666
Iteration 1: Log likelihood = -11656.381
Iteration 2: Log likelihood = -11603.266
Iteration 3: Log likelihood = -11601.752
Iteration 4: Log likelihood = -11601.751
Iteration 5: Log likelihood = -11601.751
Generalized linear models Number of obs = 20,000
Optimization : ML Residual df = 19,998
Scale parameter = 1
Deviance = 23203.50245 (1/df) Deviance = 1.160291
Pearson = 20000 (1/df) Pearson = 1.0001
Variance function: V(u) = u*(1-u) [Bernoulli]
Link function : g(u) = ln(u) [Log]
AIC = 1.160375
Log likelihood = -11601.75122 BIC = -174846.4
------------------------------------------------------------------------------
| OIM
delirium | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
surg24 | -.9392798 .0267983 -35.05 0.000 -.9918036 -.886756
_cons | -.8693602 .0109964 -79.06 0.000 -.8909127 -.8478077
------------------------------------------------------------------------------
. * model-based CI (identical bread in both languages for a saturated model)
. ciz c2.rr_crude_logbin surg24 exp
CANON w1.c2.rr_crude_logbin 0.3909
CANON w1.c2.rr_crude_logbin.lo 0.3709
CANON w1.c2.rr_crude_logbin.hi 0.4120
. * log-RR SE, log-binomial, model-based
. canon c2.se_crude_logbin _se[surg24]
CANON w1.c2.se_crude_logbin 0.0268
. quietly glm delirium surg24, family(poisson) link(log)
. scalar sc_pn = _se[surg24]
. glm delirium surg24, family(poisson) link(log) vce(robust)
Iteration 0: Log pseudolikelihood = -13125.243
Iteration 1: Log pseudolikelihood = -12910.664
Iteration 2: Log pseudolikelihood = -12910.633
Iteration 3: Log pseudolikelihood = -12910.633
Generalized linear models Number of obs = 20,000
Optimization : ML Residual df = 19,998
Scale parameter = 1
Deviance = 13415.26583 (1/df) Deviance = .6708304
Pearson = 13797 (1/df) Pearson = .689919
Variance function: V(u) = u [Poisson]
Link function : g(u) = ln(u) [Log]
AIC = 1.291263
Log pseudolikelihood = -12910.63291 BIC = -184634.7
------------------------------------------------------------------------------
| Robust
delirium | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
surg24 | -.9392798 .026799 -35.05 0.000 -.9918049 -.8867547
_cons | -.8693602 .0109967 -79.06 0.000 -.8909132 -.8478071
------------------------------------------------------------------------------
. matrix bpc = e(b)
. * modified Poisson, robust CI
. ciz c2.rr_crude_poisson surg24 exp
CANON w1.c2.rr_crude_poisson 0.3909
CANON w1.c2.rr_crude_poisson.lo 0.3709
CANON w1.c2.rr_crude_poisson.hi 0.4120
. * Poisson model-based log-RR SE (too large for a binary outcome) versus the sandwich log-RR SE
. canon c2.se_crude_poisson_naive scalar(sc_pn)
CANON w1.c2.se_crude_poisson_naive 0.0304
. canon c2.se_crude_poisson_robust _se[surg24]
CANON w1.c2.se_crude_poisson_robust 0.0268
. * naive Poisson 95% CI
. canon c2.rr_crude_poisson_naive.lo exp(_b[surg24] - invnormal(0.975)*scalar(sc_pn))
CANON w1.c2.rr_crude_poisson_naive.lo 0.3683
. canon c2.rr_crude_poisson_naive.hi exp(_b[surg24] + invnormal(0.975)*scalar(sc_pn))
CANON w1.c2.rr_crude_poisson_naive.hi 0.4149R: การประมาณชุดเดียวกัน ด้วยแพ็กเกจ sandwich
cb <- rr_glm(delirium ~ surg24, binomial(link = "log"), d)
canon_ci("c2.rr_crude_logbin", cb$est, cb$lo, cb$hi) # model-based SE (identical bread in both languages here)
canon("c2.se_crude_logbin", cb$se) # log-RR SE, log-binomial, model-based
cp <- rr_glm(delirium ~ surg24, poisson(link = "log"), d, TRUE)
canon_ci("c2.rr_crude_poisson", cp$est, cp$lo, cp$hi) # modified Poisson, robust CI
canon("c2.se_crude_poisson_naive", cp$se_naive) # Poisson model-based log-RR SE (too large for a binary outcome)
canon("c2.se_crude_poisson_robust", cp$se) # sandwich log-RR SE
canon("c2.rr_crude_poisson_naive.lo", exp(log(cp$est) - z * cp$se_naive)) # naive Poisson 95% CI, lower
canon("c2.rr_crude_poisson_naive.hi", exp(log(cp$est) + z * cp$se_naive)) # naive Poisson 95% CI, upper> cb <- rr_glm(delirium ~ surg24, binomial(link = "log"),
+ d)
> canon_ci("c2.rr_crude_logbin", cb$est, cb$lo, cb$hi)
CANON w1.c2.rr_crude_logbin 0.3909
CANON w1.c2.rr_crude_logbin.lo 0.3709
CANON w1.c2.rr_crude_logbin.hi 0.4120
> canon("c2.se_crude_logbin", cb$se)
CANON w1.c2.se_crude_logbin 0.0268
> cp <- rr_glm(delirium ~ surg24, poisson(link = "log"),
+ d, TRUE)
> canon_ci("c2.rr_crude_poisson", cp$est, cp$lo, cp$hi)
CANON w1.c2.rr_crude_poisson 0.3909
CANON w1.c2.rr_crude_poisson.lo 0.3709
CANON w1.c2.rr_crude_poisson.hi 0.4120
> canon("c2.se_crude_poisson_naive", cp$se_naive)
CANON w1.c2.se_crude_poisson_naive 0.0304
> canon("c2.se_crude_poisson_robust", cp$se)
CANON w1.c2.se_crude_poisson_robust 0.0268
> canon("c2.rr_crude_poisson_naive.lo", exp(log(cp$est) -
+ z * cp$se_naive))
CANON w1.c2.rr_crude_poisson_naive.lo 0.3683
> canon("c2.rr_crude_poisson_naive.hi", exp(log(cp$est) +
+ z * cp$se_naive))
CANON w1.c2.rr_crude_poisson_naive.hi 0.4149ปรับด้วยอายุและเพศ: ความต่างเล็กน้อย ไม่ใช่ความผิดพลาด
เมื่อเพิ่มอายุและเพศ ผู้ป่วยแต่ละคนจะมีความเสี่ยงที่แบบจำลองพอดีเป็นของตัวเอง น้ำหนักจึงต่างกันภายในแต่ละกลุ่ม น้ำหนักของ log-binomial คือ $1/(1-\mu_i)$ ให้ความสำคัญกับผู้ป่วยความเสี่ยงสูงมากกว่าน้ำหนัก 1 ของ Poisson ดัดแปลง หากแบบจำลองค่าเฉลี่ยถูกต้องสมบูรณ์ ค่าประมาณทั้งสองจะมีความคงเส้นคงวา แต่ก็ยังต่างกันเล็กน้อยในตัวอย่างใดตัวอย่างหนึ่ง
ในทะเบียน ค่าทั้งสองคือ 0.5196 และ 0.5179 ข้อมูลจำลองชุดนี้สร้างขึ้นให้การผ่าตัดเร็วมีผลต่างกันระหว่างผู้ป่วยที่เปราะบางกับไม่เปราะบาง จึงไม่มี risk ratio ค่าเดียวที่อธิบายผู้ป่วยทุกคนได้ การถ่วงน้ำหนักสองแบบจึงเฉลี่ยผลที่แตกต่างกันออกมาเป็นค่าสรุปที่ต่างกันเล็กน้อย ความต่างเล็กน้อยแบบนี้เป็นสิ่งที่คาดไว้ และไม่ใช่ความผิดพลาด
ความต่างเล็กน้อยไม่ได้แสดงว่าแบบจำลองค่าเฉลี่ยถูกต้อง เพราะทั้งสองการประมาณใช้แบบจำลองนั้นร่วมกัน ส่วนความต่างมากเป็นเหตุผลให้ตรวจแบบจำลองค่าเฉลี่ยว่าขาดปฏิกิริยาสัมพันธ์ (interaction) ผลของอายุที่โค้ง หรือค่าที่แบบจำลองพอดีสุดโต่งหรือไม่
ปรับด้วยอายุและเพศ: ใกล้กันแต่ไม่เท่ากัน
| แบบจำลอง | Risk ratio | ช่วงเชื่อมั่น 95% | ค่าคลาดเคลื่อนมาตรฐานของ log risk ratio |
|---|---|---|---|
| Log-binomial, Stata (observed information) | 0.5196 | 0.4924 ถึง 0.5484 | 0.0275 |
| Log-binomial, R (expected information) | 0.5196 | 0.4928 ถึง 0.5479 | 0.0270 |
| Poisson ดัดแปลง ค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช | 0.5179 | 0.4911 ถึง 0.5463 | 0.0272 |
| Poisson ค่าคลาดเคลื่อนมาตรฐานแบบอิงแบบจำลอง | 0.5179 | ไม่ได้คำนวณ | 0.0315 |
ค่าคลาดเคลื่อนมาตรฐานสองค่าสำหรับแบบจำลอง log-binomial เดียว
ซอฟต์แวร์สร้างค่าคลาดเคลื่อนมาตรฐานแบบอิงแบบจำลองจากinformation (สารสนเทศ) คือความโค้งของลอการิทึมของภาวะน่าจะเป็น (log-likelihood) ซึ่งเป็นปริมาณที่การประมาณทำให้สูงสุด ณ จุดยอด ค่าเริ่มต้นของ Stata ใช้ความโค้งที่สังเกตได้ และให้ค่าคลาดเคลื่อนมาตรฐาน 0.0275 ส่วน glm ของ R ใช้ค่าคาดหมายของมัน และให้ 0.0270
ทั้งสองแบบตรงกันสำหรับ canonical link ของแต่ละ family คือ link ที่สมการสกอร์ให้ส่วนเหลือทุกตัวมีน้ำหนัก 1 ได้แก่ logit สำหรับ binomial และ log สำหรับ Poisson สำหรับ log link กับ binomial family ทั้งสองแบบโดยทั่วไปต่างกัน เมื่อตัวแปรร่วมทำให้ผู้ป่วยมีความเสี่ยงที่แบบจำลองพอดีต่างกัน ในการประมาณแบบอิ่มตัวข้างต้นทั้งสองแบบตรงกัน [2]
Stata: ปรับด้วยอายุและเพศ
* adjusted for age and sex
quietly glm delirium surg24 age_c female, family(poisson) link(log)
scalar sc_p2n = _se[surg24]
glm delirium surg24 age_c female, family(poisson) link(log) vce(robust)
matrix bp2 = e(b)
scalar sc_p2 = exp(_b[surg24])
scalar sc_p2lo = exp(_b[surg24] - invnormal(0.975)*_se[surg24])
scalar sc_p2hi = exp(_b[surg24] + invnormal(0.975)*_se[surg24])
scalar sc_p2r = _se[surg24]
* adjusted fits: the log-binomial starts from the Poisson fit
glm delirium surg24 age_c female, family(binomial) link(log) from(bp2). * adjusted for age and sex
. quietly glm delirium surg24 age_c female, family(poisson) link(log)
. scalar sc_p2n = _se[surg24]
. glm delirium surg24 age_c female, family(poisson) link(log) vce(robust)
Iteration 0: Log pseudolikelihood = -12338.777
Iteration 1: Log pseudolikelihood = -12134.247
Iteration 2: Log pseudolikelihood = -12134.177
Iteration 3: Log pseudolikelihood = -12134.177
Generalized linear models Number of obs = 20,000
Optimization : ML Residual df = 19,996
Scale parameter = 1
Deviance = 11862.35488 (1/df) Deviance = .5932364
Pearson = 13271.16135 (1/df) Pearson = .6636908
Variance function: V(u) = u [Poisson]
Link function : g(u) = ln(u) [Log]
AIC = 1.213818
Log pseudolikelihood = -12134.17744 BIC = -186167.8
------------------------------------------------------------------------------
| Robust
delirium | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
surg24 | -.6578832 .0271581 -24.22 0.000 -.7111121 -.6046544
age_c | .0450466 .0008786 51.27 0.000 .0433246 .0467686
female | -.0201594 .020605 -0.98 0.328 -.0605445 .0202257
_cons | -1.100182 .0202124 -54.43 0.000 -1.139798 -1.060567
------------------------------------------------------------------------------
. matrix bp2 = e(b)
. scalar sc_p2 = exp(_b[surg24])
. scalar sc_p2lo = exp(_b[surg24] - invnormal(0.975)*_se[surg24])
. scalar sc_p2hi = exp(_b[surg24] + invnormal(0.975)*_se[surg24])
. scalar sc_p2r = _se[surg24]
. * adjusted fits: the log-binomial starts from the Poisson fit
. glm delirium surg24 age_c female, family(binomial) link(log) from(bp2)
Iteration 0: Log likelihood = -10500.15
Iteration 1: Log likelihood = -10413.165
Iteration 2: Log likelihood = -10400.121
Iteration 3: Log likelihood = -10399.792
Iteration 4: Log likelihood = -10399.791
Generalized linear models Number of obs = 20,000
Optimization : ML Residual df = 19,996
Scale parameter = 1
Deviance = 20799.58236 (1/df) Deviance = 1.040187
Pearson = 19229.2475 (1/df) Pearson = .9616547
Variance function: V(u) = u*(1-u) [Bernoulli]
Link function : g(u) = ln(u) [Log]
AIC = 1.040379
Log likelihood = -10399.79118 BIC = -177230.6
------------------------------------------------------------------------------
| OIM
delirium | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
surg24 | -.654659 .0274789 -23.82 0.000 -.7085167 -.6008014
age_c | .039937 .000757 52.75 0.000 .0384532 .0414208
female | .0025709 .0169232 0.15 0.879 -.030598 .0357398
_cons | -1.098012 .0185469 -59.20 0.000 -1.134363 -1.06166
------------------------------------------------------------------------------R: ปรับด้วยอายุและเพศ แล้วปรับด้วยอายุและความเปราะบาง
# adjusted fits: the log-binomial needs starting values (taken from the Poisson fit) to converge
p2 <- rr_glm(delirium ~ surg24 + age_c + female, poisson(link = "log"), d, TRUE)
l2 <- rr_glm(delirium ~ surg24 + age_c + female, binomial(link = "log"), d, start = coef(p2$fit))
canon("c2.rr_adj2_logbin", l2$est)
# R's glm uses the expected-information SE for the non-canonical log link (Stata's ML glm: observed)
canon("c2.rr_adj2_logbin.lo.r", l2$lo); canon("c2.rr_adj2_logbin.hi.r", l2$hi)
canon("c2.se_adj2_logbin.r", l2$se) # log-RR SE, adjusted log-binomial (expected information)
canon("c2.rr_adj2_poisson", p2$est)
canon("c2.rr_adj2_poisson.lo", p2$lo); canon("c2.rr_adj2_poisson.hi", p2$hi) # modified Poisson, robust CI
canon("c2.se_adj2_poisson_naive", p2$se_naive) # Poisson model-based log-RR SE, adjusted
canon("c2.se_adj2_poisson_robust", p2$se) # sandwich log-RR SE, adjusted
g2 <- rr_glm(delirium ~ surg24 + age_c + female, gaussian(link = "log"), d, TRUE, start = coef(p2$fit))
canon("c2.rr_adj2_gaussian", g2$est) # Gaussian log link, adjusted for age and sex
canon("c2.rr_adj2_gaussian.lo.r", g2$lo); canon("c2.rr_adj2_gaussian.hi.r", g2$hi) # expected-information bread
# adjusted for age and frailty: the Poisson fit gives fitted risks above 1, so the log-binomial cannot start
paf <- rr_glm(delirium ~ surg24 + age_c + frail, poisson(link = "log"), d, TRUE)
canon("c2.rr_adjaf_poisson", paf$est)
canon("c2.rr_adjaf_poisson.lo", paf$lo); canon("c2.rr_adjaf_poisson.hi", paf$hi)
canon("c2.se_adjaf_poisson_naive", paf$se_naive) # Poisson model-based log-RR SE, age and frailty
canon("c2.se_adjaf_poisson_robust", paf$se) # sandwich log-RR SE, age and frailty
canon("c2.adjaf_poisson_max_fitted", max(fitted(paf$fit))) # largest fitted risk: above 1 means the log-binomial wall
lbaf <- tryCatch(glm(delirium ~ surg24 + age_c + frail, family = binomial(link = "log"), data = d, start = coef(paf$fit)),
error = function(e) { cat("NOTE log-binomial (age, frailty) stopped:", conditionMessage(e), "\n"); NULL })
canon_n("c2.adjaf_logbin_converged.r", !is.null(lbaf) && isTRUE(lbaf$converged)) # 1 = converged> p2 <- rr_glm(delirium ~ surg24 + age_c + female, poisson(link = "log"),
+ d, TRUE)
> l2 <- rr_glm(delirium ~ surg24 + age_c + female, binomial(link = "log"),
+ d, start = coef(p2$fit))
> canon("c2.rr_adj2_logbin", l2$est)
CANON w1.c2.rr_adj2_logbin 0.5196
> canon("c2.rr_adj2_logbin.lo.r", l2$lo)
CANON w1.c2.rr_adj2_logbin.lo.r 0.4928
> canon("c2.rr_adj2_logbin.hi.r", l2$hi)
CANON w1.c2.rr_adj2_logbin.hi.r 0.5479
> canon("c2.se_adj2_logbin.r", l2$se)
CANON w1.c2.se_adj2_logbin.r 0.0270
> canon("c2.rr_adj2_poisson", p2$est)
CANON w1.c2.rr_adj2_poisson 0.5179
> canon("c2.rr_adj2_poisson.lo", p2$lo)
CANON w1.c2.rr_adj2_poisson.lo 0.4911
> canon("c2.rr_adj2_poisson.hi", p2$hi)
CANON w1.c2.rr_adj2_poisson.hi 0.5463
> canon("c2.se_adj2_poisson_naive", p2$se_naive)
CANON w1.c2.se_adj2_poisson_naive 0.0315
> canon("c2.se_adj2_poisson_robust", p2$se)
CANON w1.c2.se_adj2_poisson_robust 0.0272
> g2 <- rr_glm(delirium ~ surg24 + age_c + female, gaussian(link = "log"),
+ d, TRUE, start = coef(p2$fit))
> canon("c2.rr_adj2_gaussian", g2$est)
CANON w1.c2.rr_adj2_gaussian 0.5294
> canon("c2.rr_adj2_gaussian.lo.r", g2$lo)
CANON w1.c2.rr_adj2_gaussian.lo.r 0.5021
> canon("c2.rr_adj2_gaussian.hi.r", g2$hi)
CANON w1.c2.rr_adj2_gaussian.hi.r 0.5582
> paf <- rr_glm(delirium ~ surg24 + age_c + frail, poisson(link = "log"),
+ d, TRUE)
> canon("c2.rr_adjaf_poisson", paf$est)
CANON w1.c2.rr_adjaf_poisson 0.6413
> canon("c2.rr_adjaf_poisson.lo", paf$lo)
CANON w1.c2.rr_adjaf_poisson.lo 0.6094
> canon("c2.rr_adjaf_poisson.hi", paf$hi)
CANON w1.c2.rr_adjaf_poisson.hi 0.6749
> canon("c2.se_adjaf_poisson_naive", paf$se_naive)
CANON w1.c2.se_adjaf_poisson_naive 0.0320
> canon("c2.se_adjaf_poisson_robust", paf$se)
CANON w1.c2.se_adjaf_poisson_robust 0.0260
> canon("c2.adjaf_poisson_max_fitted", max(fitted(paf$fit)))
CANON w1.c2.adjaf_poisson_max_fitted 1.1066
> lbaf <- tryCatch(glm(delirium ~ surg24 + age_c + frail,
+ family = binomial(link = "log"), data = d, start = coef(paf$fit)),
+ error = function(e) {
+ cat("NOTE log-binomial (age, frailty) stopped:", conditionMessage(e),
+ "\n")
+ NULL
+ })
NOTE log-binomial (age, frailty) stopped: cannot find valid starting values: please specify some
> canon_n("c2.adjaf_logbin_converged.r", !is.null(lbaf) &&
+ isTRUE(lbaf$converged))
CANON w1.c2.adjaf_logbin_converged.r 0อายุและความเปราะบาง: ทะเบียนชนกำแพง
เมื่อใส่อายุและความเปราะบางในแบบจำลอง Poisson ดัดแปลงให้ risk ratio 0.6413 (ช่วงเชื่อมั่น 95% คือ 0.6094 ถึง 0.6749) ค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิชของ log risk ratio คือ 0.0260 เทียบกับ 0.0320 แบบอิงแบบจำลอง ค่าที่แบบจำลองพอดีสูงสุดคือ 1.1066 ในผู้ป่วยเปราะบางที่อายุมากที่สุดและไม่ได้ผ่าตัดเร็ว
การจำลองข้อมูลกำหนดให้ความเสี่ยงจริงของภาวะสับสนเฉียบพลันในทะเบียนนี้สูงถึง 0.991 แต่ไม่เกิน 1 ผลของอายุแบบเส้นตรงบนสเกลลอการิทึมไม่อาจแบนราบลงในแบบเดียวกัน ค่าที่แบบจำลองพอดีของผู้ป่วยเปราะบางที่อายุมากที่สุดจึงเกิน 1
แบบจำลอง log-binomial ที่เริ่มจากสัมประสิทธิ์ Poisson เหล่านี้ไม่ลู่เข้า Stata วนซ้ำจนถึงขีดจำกัดและเตือนว่ามีค่าเฉลี่ยที่เป็นไปไม่ได้ (inadmissible) คือค่าที่ความเสี่ยงเป็นไม่ได้ ส่วน glm ของ R หยุดด้วยข้อผิดพลาดก่อนการวนซ้ำครั้งแรก เพราะค่าเริ่มต้นทำให้ค่าเฉลี่ยที่แบบจำลองพอดีบางค่าเกิน 1 อยู่แล้ว
ขั้นต่อไปที่กำหนดไว้ล่วงหน้า คือการถอยจากแบบจำลอง log-binomial ไปยัง Poisson ดัดแปลงและวิธีอื่นต่อไป เป็นเรื่องของตอนที่ 3 ของชุดบทความนี้
Stata: อายุและความเปราะบาง จุดที่การประมาณ log-binomial ชนกำแพง
quietly glm delirium surg24 age_c frail, family(poisson) link(log)
scalar sc_pafn = _se[surg24]
glm delirium surg24 age_c frail, family(poisson) link(log) vce(robust)
matrix bpaf = e(b)
ciz c2.rr_adjaf_poisson surg24 exp
* model-based versus sandwich log-RR SE, age and frailty
canon c2.se_adjaf_poisson_naive scalar(sc_pafn)
canon c2.se_adjaf_poisson_robust _se[surg24]
predict double mu_af, mu
summarize mu_af
* largest fitted risk: above 1 means the log-binomial wall
canon c2.adjaf_poisson_max_fitted r(max)
capture noisily glm delirium surg24 age_c frail, family(binomial) link(log) from(bpaf) iterate(100). quietly glm delirium surg24 age_c frail, family(poisson) link(log)
. scalar sc_pafn = _se[surg24]
. glm delirium surg24 age_c frail, family(poisson) link(log) vce(robust)
Iteration 0: Log pseudolikelihood = -11571.438
Iteration 1: Log pseudolikelihood = -11364.541
Iteration 2: Log pseudolikelihood = -11364.358
Iteration 3: Log pseudolikelihood = -11364.358
Generalized linear models Number of obs = 20,000
Optimization : ML Residual df = 19,996
Scale parameter = 1
Deviance = 10322.71624 (1/df) Deviance = .5162391
Pearson = 12654.60943 (1/df) Pearson = .632857
Variance function: V(u) = u [Poisson]
Link function : g(u) = ln(u) [Log]
AIC = 1.136836
Log pseudolikelihood = -11364.35812 BIC = -187707.4
------------------------------------------------------------------------------
| Robust
delirium | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
surg24 | -.4442503 .0260218 -17.07 0.000 -.4952521 -.3932485
age_c | .0323723 .0009141 35.41 0.000 .0305807 .034164
frail | 1.152159 .0275096 41.88 0.000 1.098241 1.206077
_cons | -1.827827 .0254834 -71.73 0.000 -1.877774 -1.777881
------------------------------------------------------------------------------
. matrix bpaf = e(b)
. ciz c2.rr_adjaf_poisson surg24 exp
CANON w1.c2.rr_adjaf_poisson 0.6413
CANON w1.c2.rr_adjaf_poisson.lo 0.6094
CANON w1.c2.rr_adjaf_poisson.hi 0.6749
. * model-based versus sandwich log-RR SE, age and frailty
. canon c2.se_adjaf_poisson_naive scalar(sc_pafn)
CANON w1.c2.se_adjaf_poisson_naive 0.0320
. canon c2.se_adjaf_poisson_robust _se[surg24]
CANON w1.c2.se_adjaf_poisson_robust 0.0260
. predict double mu_af, mu
. summarize mu_af
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
mu_af | 20,000 .31015 .2711498 .0390369 1.106572
. * largest fitted risk: above 1 means the log-binomial wall
. canon c2.adjaf_poisson_max_fitted r(max)
CANON w1.c2.adjaf_poisson_max_fitted 1.1066
. capture noisily glm delirium surg24 age_c frail, family(binomial) link(log) from(bpaf) iterate(100)
Iteration 0: Log likelihood = -10777.424 (not concave)
Iteration 1: Log likelihood = -10776.706 (not concave)
Iteration 2: Log likelihood = -10776.648 (not concave)
Iteration 3: Log likelihood = -10776.642 (not concave)
Iteration 4: Log likelihood = -10776.642 (not concave)
Iteration 5: Log likelihood = -10776.642 (not concave)
Iteration 6: Log likelihood = -10776.642 (not concave)
Iteration 7: Log likelihood = -10776.642 (not concave)
Iteration 8: Log likelihood = -10776.642 (not concave)
Iteration 9: Log likelihood = -10776.642 (not concave)
Iteration 10: Log likelihood = -10776.642 (not concave)
Iteration 11: Log likelihood = -10776.642 (not concave)
Iteration 12: Log likelihood = -10776.642 (not concave)
Iteration 13: Log likelihood = -10776.642 (not concave)
Iteration 14: Log likelihood = -10776.642 (not concave)
Iteration 15: Log likelihood = -10776.642 (not concave)
Iteration 16: Log likelihood = -10776.642 (not concave)
Iteration 17: Log likelihood = -10776.642 (not concave)
Iteration 18: Log likelihood = -10776.642 (not concave)
Iteration 19: Log likelihood = -10776.642 (not concave)
Iteration 20: Log likelihood = -10776.642 (not concave)
Iteration 21: Log likelihood = -10776.642 (not concave)
Iteration 22: Log likelihood = -10776.642 (not concave)
Iteration 23: Log likelihood = -10776.641 (not concave)
Iteration 24: Log likelihood = -10776.641 (not concave)
Iteration 25: Log likelihood = -10776.641 (not concave)
Iteration 26: Log likelihood = -10776.641 (not concave)
Iteration 27: Log likelihood = -10776.641 (not concave)
Iteration 28: Log likelihood = -10776.641 (not concave)
Iteration 29: Log likelihood = -10776.641 (not concave)
Iteration 30: Log likelihood = -10776.641 (not concave)
Iteration 31: Log likelihood = -10776.641 (not concave)
Iteration 32: Log likelihood = -10776.641 (not concave)
Iteration 33: Log likelihood = -10776.641 (not concave)
Iteration 34: Log likelihood = -10776.641 (not concave)
Iteration 35: Log likelihood = -10776.641 (not concave)
Iteration 36: Log likelihood = -10776.641 (not concave)
Iteration 37: Log likelihood = -10776.641 (not concave)
Iteration 38: Log likelihood = -10776.641 (not concave)
Iteration 39: Log likelihood = -10776.641 (not concave)
Iteration 40: Log likelihood = -10776.641 (not concave)
Iteration 41: Log likelihood = -10776.641 (not concave)
Iteration 42: Log likelihood = -10776.641 (not concave)
Iteration 43: Log likelihood = -10776.641 (not concave)
Iteration 44: Log likelihood = -10776.641 (not concave)
Iteration 45: Log likelihood = -10776.641 (not concave)
Iteration 46: Log likelihood = -10776.641 (not concave)
Iteration 47: Log likelihood = -10776.641 (not concave)
Iteration 48: Log likelihood = -10776.641 (not concave)
Iteration 49: Log likelihood = -10776.641 (not concave)
Iteration 50: Log likelihood = -10776.641 (not concave)
Iteration 51: Log likelihood = -10776.641 (not concave)
Iteration 52: Log likelihood = -10776.641 (not concave)
Iteration 53: Log likelihood = -10776.64 (not concave)
Iteration 54: Log likelihood = -10776.64 (not concave)
Iteration 55: Log likelihood = -10776.64 (not concave)
Iteration 56: Log likelihood = -10776.64 (not concave)
Iteration 57: Log likelihood = -10776.64 (not concave)
Iteration 58: Log likelihood = -10776.64 (not concave)
Iteration 59: Log likelihood = -10776.64 (not concave)
Iteration 60: Log likelihood = -10776.64 (not concave)
Iteration 61: Log likelihood = -10776.64 (not concave)
Iteration 62: Log likelihood = -10776.64 (not concave)
Iteration 63: Log likelihood = -10776.64 (not concave)
Iteration 64: Log likelihood = -10776.64 (not concave)
Iteration 65: Log likelihood = -10776.64 (not concave)
Iteration 66: Log likelihood = -10776.64 (not concave)
Iteration 67: Log likelihood = -10776.64 (not concave)
Iteration 68: Log likelihood = -10776.64 (not concave)
Iteration 69: Log likelihood = -10776.64 (not concave)
Iteration 70: Log likelihood = -10776.64 (not concave)
Iteration 71: Log likelihood = -10776.64 (not concave)
Iteration 72: Log likelihood = -10776.64 (not concave)
Iteration 73: Log likelihood = -10776.64 (not concave)
Iteration 74: Log likelihood = -10776.64 (not concave)
Iteration 75: Log likelihood = -10776.64 (not concave)
Iteration 76: Log likelihood = -10776.64 (not concave)
Iteration 77: Log likelihood = -10776.64 (not concave)
Iteration 78: Log likelihood = -10776.64 (not concave)
Iteration 79: Log likelihood = -10776.64 (not concave)
Iteration 80: Log likelihood = -10776.64 (not concave)
Iteration 81: Log likelihood = -10776.64 (not concave)
Iteration 82: Log likelihood = -10776.639 (not concave)
Iteration 83: Log likelihood = -10776.639 (not concave)
Iteration 84: Log likelihood = -10776.639 (not concave)
Iteration 85: Log likelihood = -10776.639 (not concave)
Iteration 86: Log likelihood = -10776.639 (not concave)
Iteration 87: Log likelihood = -10776.639 (not concave)
Iteration 88: Log likelihood = -10776.639 (not concave)
Iteration 89: Log likelihood = -10776.639 (not concave)
Iteration 90: Log likelihood = -10776.639 (not concave)
Iteration 91: Log likelihood = -10776.639 (not concave)
Iteration 92: Log likelihood = -10776.639 (not concave)
Iteration 93: Log likelihood = -10776.639 (not concave)
Iteration 94: Log likelihood = -10776.639 (not concave)
Iteration 95: Log likelihood = -10776.639 (not concave)
Iteration 96: Log likelihood = -10776.639 (not concave)
Iteration 97: Log likelihood = -10776.639 (not concave)
Iteration 98: Log likelihood = -10776.639 (not concave)
Iteration 99: Log likelihood = -10776.639 (not concave)
Iteration 100: Log likelihood = -10776.639 (not concave)
convergence not achieved
Generalized linear models Number of obs = 20,000
Optimization : ML Residual df = 19,998
Scale parameter = 1
Deviance = 21553.27781 (1/df) Deviance = 1.077772
Pearson = 1.00000e+10 (1/df) Pearson = 500050.9
Variance function: V(u) = u*(1-u) [Bernoulli]
Link function : g(u) = ln(u) [Log]
AIC = 1.077864
Log likelihood = -10776.6389 BIC = -176496.7
------------------------------------------------------------------------------
| OIM
delirium | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
surg24 | -.4160061 .0221925 -18.75 0.000 -.4595027 -.3725096
age_c | .0323723 7.76e-10 4.2e+07 0.000 .0323723 .0323723
frail | 1.152157 1.78e-08 6.5e+07 0.000 1.152157 1.152157
_cons | -1.827829 . . . . .
------------------------------------------------------------------------------
Warning: Parameter estimates produce inadmissible mean estimates in one or
more observations.
Warning: Convergence not achieved.ความเข้าใจผิดที่พบบ่อยและวิธีแก้
-
"family เพียงบอกชนิดของผลลัพธ์ ดังนั้นจะเลือกแบบใดก็ไม่สำคัญ"
ในทะเบียน risk ratio แบบ log link ค่าเดียวกันคือ 0.3909 มีค่าคลาดเคลื่อนมาตรฐานของ log risk ratio เท่ากับ 0.0268 เมื่อใช้ binomial family และ 0.0304 เมื่อใช้ความแปรปรวนแบบอิงแบบจำลองของ Poisson
วิธีแก้: family กำหนดฟังก์ชันความแปรปรวน ซึ่งเป็นตัวขับค่าคลาดเคลื่อนมาตรฐาน และในแบบจำลองที่ปรับตัวแปรร่วมก็ขับน้ำหนักของข้อสังเกต Poisson family กับผลลัพธ์แบบสองค่าต้องใช้ค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช
-
"robust standard error ช่วยแก้แบบจำลองที่ระบุผิด"
Poisson ดัดแปลงต้องมีแบบจำลองค่าเฉลี่ยที่ถูกต้อง แซนด์วิชเพียงแทนที่ความแปรปรวนแบบ Poisson ซึ่งผิดสำหรับผลลัพธ์แบบสองค่า
วิธีแก้: robust standard error ช่วยแก้ค่าคลาดเคลื่อนมาตรฐานเมื่อข้อสมมติเรื่องความแปรปรวนผิด แต่ไม่ได้แก้แบบจำลองค่าเฉลี่ยที่ผิด และสัมประสิทธิ์ที่ลำเอียงก็ยังคงลำเอียงอยู่
-
"risk ratio ของ log-binomial กับ Poisson ดัดแปลงที่ปรับแล้วต่างกัน แสดงว่าอันใดอันหนึ่งผิด"
สองแบบถ่วงน้ำหนักผู้ป่วยด้วยฟังก์ชันความแปรปรวนต่างกัน ได้ 0.5196 และ 0.5179 หลังปรับด้วยอายุและเพศในทะเบียน
วิธีแก้: ควรคาดว่ามีช่องว่างเล็กน้อย แต่อย่าอ่านว่าเป็นหลักฐานว่าแบบจำลองค่าเฉลี่ยถูกต้อง เพราะทั้งสองแบบใช้แบบจำลองนั้นร่วมกัน เมื่อช่องว่างมาก ให้ตรวจแบบจำลองค่าเฉลี่ย เช่นปฏิกิริยาสัมพันธ์ รูปร่างของตัวแปรร่วมต่อเนื่อง และค่าที่แบบจำลองพอดีสุดโต่ง ก่อนโทษวิธีใดวิธีหนึ่ง
-
"Poisson ดัดแปลงใช้ได้เฉพาะกับผลลัพธ์ที่พบน้อย"
มันประมาณ risk ratio ได้กับผลลัพธ์ที่พบบ่อยด้วย เช่นภาวะสับสนเฉียบพลันในทะเบียน
วิธีแก้: ความพบน้อยสำคัญกับคำถามอีกข้อหนึ่ง คือ odds ratio ใกล้เคียง risk ratio หรือไม่
สิ่งที่ควรทำในการวิเคราะห์ของคุณเอง
- เขียนมาตรวัดผลลงในโปรโตคอล แล้วเลือก link ที่ให้มาตรวัดนั้น คือ identity, log หรือ logit
- เมื่อ family ไม่ตรงกับผลลัพธ์ เช่นใช้ Poisson family กับผลลัพธ์แบบสองค่า ให้รายงานค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช
- เมื่อเปรียบเทียบค่าคลาดเคลื่อนมาตรฐานของ log-binomial ระหว่างโปรแกรม ให้ตรวจว่าแต่ละโปรแกรมใช้ observed information หรือ expected information
อภิธานศัพท์
- linear predictor (ตัวทำนายเชิงเส้น (linear predictor))
- ผลรวมถ่วงน้ำหนักของตัวแปรร่วมของผู้ป่วย คือ $\eta$
- link function (ฟังก์ชันเชื่อมโยง (link function))
- ฟังก์ชัน $g(\cdot)$ ที่เชื่อมค่าเฉลี่ยของผลลัพธ์กับ $\eta$ และกำหนดมาตรวัดผล
- family (family (การแจกแจงของผลลัพธ์))
- การแจกแจงที่สมมติให้ผลลัพธ์ ซึ่งให้ฟังก์ชันความแปรปรวน
- variance function (ฟังก์ชันความแปรปรวน (variance function))
- $V(\mu)$ คือความแปรปรวนขึ้นกับค่าเฉลี่ยอย่างไร: $\mu(1-\mu)$ สำหรับ binomial และ $\mu$ สำหรับ Poisson
- log-binomial model (แบบจำลอง log-binomial)
- GLM แบบ binomial ที่ใช้ log link ให้ risk ratio
- modified Poisson regression (Poisson ดัดแปลง (modified Poisson regression))
- GLM แบบ Poisson ที่ใช้ log link และค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช ประมาณกับผลลัพธ์แบบสองค่า
- boundary
- ขอบของสัมประสิทธิ์ที่ยอมให้มีได้ สำหรับแบบจำลอง log-binomial คือความเสี่ยงที่แบบจำลองพอดีเท่ากับ 1
- convergence
- การลู่เข้า (convergence): จุดที่ค่าประมาณต่อเนื่องของขั้นตอนคำนวณหยุดเปลี่ยนแปลง
- sandwich standard error (ค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช (robust standard error))
- ค่าคลาดเคลื่อนมาตรฐานจากส่วนเหลือที่สังเกตได้ ใช้ได้เมื่อฟังก์ชันความแปรปรวนผิดแต่แบบจำลองค่าเฉลี่ยถูกต้อง
เอกสารอ้างอิง
- Nelder JA, Wedderburn RWM. Generalized linear models. J R Stat Soc Ser A. 1972;135(3):370-384. doi:10.2307/2344614 https://doi.org/10.2307/2344614
- McCullagh P, Nelder JA. Generalized linear models. 2nd ed. London: Chapman and Hall; 1989. doi:10.1007/978-1-4899-3242-6 https://doi.org/10.1007/978-1-4899-3242-6
- Wacholder S. Binomial regression in GLIM: estimating risk ratios and risk differences. Am J Epidemiol. 1986;123(1):174-184. doi:10.1093/oxfordjournals.aje.a114212 https://doi.org/10.1093/oxfordjournals.aje.a114212
- Zou G. A modified Poisson regression approach to prospective studies with binary data. Am J Epidemiol. 2004;159(7):702-706. doi:10.1093/aje/kwh090 https://doi.org/10.1093/aje/kwh090
- Williamson T, Eliasziw M, Fick GH. Log-binomial models: exploring failed convergence. Emerg Themes Epidemiol. 2013;10:14. doi:10.1186/1742-7622-10-14 https://doi.org/10.1186/1742-7622-10-14
- Naimi AI, Whitcomb BW. Estimating risk ratios and risk differences using regression. Am J Epidemiol. 2020;189(6):508-510. doi:10.1093/aje/kwaa044 https://doi.org/10.1093/aje/kwaa044
ประเด็นสำคัญ
- link กำหนดมาตรวัดผล: identity ให้ผลต่างความเสี่ยง log ให้ risk ratio และ logit ให้ odds ratio
- family กำหนดฟังก์ชันความแปรปรวน ซึ่งขับค่าคลาดเคลื่อนมาตรฐาน และเมื่อมีตัวแปรร่วม ขับน้ำหนักที่ผู้ป่วยแต่ละคนได้รับ
- ในแบบจำลอง log-binomial ความเสี่ยงที่แบบจำลองพอดีจะไม่เกิน 1 ก็ต่อเมื่อ eta ไม่เกิน 0 ดังนั้น eta ของผู้ป่วยทุกคนต้องอยู่ในช่วงนั้น ใกล้กำแพงนั้นการประมาณ log-binomial จะไม่ลู่เข้าหรือหยุดบนขอบเขต
- Poisson ดัดแปลงคง log link โดยไม่มีกำแพง และให้ risk ratio ที่คงเส้นคงวาเมื่อแบบจำลองค่าเฉลี่ยถูกต้อง พร้อมค่าคลาดเคลื่อนมาตรฐานแบบแซนด์วิช
- เมื่อมีตัวแปรสัมผัสแบบสองค่าตัวเดียว แบบจำลองทั้งสองให้ risk ratio เท่ากัน เมื่อมีตัวแปรร่วมทั้งสองถ่วงน้ำหนักผู้ป่วยต่างกันและให้ค่าต่างกันเล็กน้อย ซึ่งไม่ใช่ความผิดพลาด
อ่านต่อในวิกิ: [[risk-regression-models-epidemiology-th]] [[glm-families-effect-measures]]