เลือก vce() แบบไหนเมื่อใด: standard error แบบอิงแบบจำลอง แบบ robust หรือแบบ cluster-robust

Read the English version

บทคัดย่อ

ค่าสัมประสิทธิ์การถดถอยกับค่าคลาดเคลื่อนมาตรฐานของมันมาจากการคำนวณคนละส่วน และออปชัน vce() ของ Stata เป็นตัวกำหนดส่วนหลัง ในข้อมูลจำลองของผู้ป่วย 600 คนจากโรงพยาบาล 40 แห่ง โรงพยาบาลครึ่งหนึ่งมีโปรแกรมติดตามผู้ป่วยหลังจำหน่าย อัตราส่วนอัตราอุบัติการณ์ (incidence rate ratio) ของการกลับมานอนโรงพยาบาลซ้ำที่สัมพันธ์กับโปรแกรมเท่ากับ 0.72 ไม่ว่าจะเลือกวิธีใด ค่าคลาดเคลื่อนมาตรฐานของลอการิทึมของอัตราส่วนนี้เท่ากับ 0.0686 เมื่อเชื่อความแปรปรวนแบบ Poisson, 0.0816 เมื่อใช้ vce(robust) และ 0.1721 เมื่อใช้ vce(cluster hospital) วิธี cluster-robust รวมส่วนที่ผู้ป่วยทุกคนของโรงพยาบาลเดียวกันมีต่อค่าประมาณเข้าด้วยกันก่อนยกกำลังสอง จึงยอมให้ผู้ป่วยในโรงพยาบาลเดียวกันสัมพันธ์กันได้ ช่วงความเชื่อมั่น 95% ของวิธีนี้ คือ 0.51 ถึง 1.01 เป็นช่วงที่สะท้อนโปรแกรมซึ่งจัดสรรตามโรงพยาบาล แบบจำลอง negative binomial ให้อัตราส่วนเกือบเท่ากัน แต่ยังถือว่าผู้ป่วยเป็นอิสระต่อกัน ตัวประมาณความแปรปรวนไม่มีตัวใดเปลี่ยนสิ่งที่ประมาณ หรือขจัดตัวกวนระดับโรงพยาบาลได้ บทความนี้ให้ตารางตัดสินใจสำหรับ vce() คำสั่ง R ที่เทียบเท่ากัน และประโยคสำหรับส่วนวิธีการและส่วนผลลัพธ์ของผลลัพธ์แบบนับและแบบสองค่า

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

ค่าคลาดเคลื่อนมาตรฐานสามค่าสำหรับโปรแกรมเดียว

ทีมคุณภาพของเครือข่ายโรงพยาบาลแห่งหนึ่งศึกษาการกลับมานอนโรงพยาบาลซ้ำ (readmission) คือจำนวนครั้งที่ผู้ป่วยแต่ละคนต้องนอนโรงพยาบาลในช่วงหนึ่งปีหลังจำหน่าย โรงพยาบาลครึ่งหนึ่งของเครือข่ายมีโปรแกรมติดตามผู้ป่วยหลังจำหน่าย ทีมต้องการทราบว่าผู้ป่วยของโรงพยาบาลเหล่านั้นกลับมานอนซ้ำน้อยกว่าหรือไม่ ข้อมูลในที่นี้เป็นข้อมูลจำลอง: ผู้ป่วย 600 คน จากโรงพยาบาล 40 แห่ง แห่งละ 15 คน

แผนการวิเคราะห์ระบุ vce(robust) ผู้ทบทวนร่างรายงานขอ vce(cluster hospital) และนักสถิติเสนอให้ใช้แบบจำลอง negative binomial แทน ซึ่งเป็นแบบจำลองข้อมูลนับที่ยอมให้มีความกระจายเกินกว่า Poisson อัตราส่วนอัตราอุบัติการณ์ (incidence rate ratio) ของโปรแกรม คืออัตราส่วนของจำนวนครั้งที่นอนโรงพยาบาลเฉลี่ยต่อผู้ป่วยในโรงพยาบาลที่มีโปรแกรมต่อค่านั้นในโรงพยาบาลอื่น แทบไม่ขยับระหว่างสามวิธี ส่วนช่วงความเชื่อมั่นขยับ: สองวิธีไม่ครอบคลุม 1 และอีกหนึ่งวิธีครอบคลุมค่านั้น

บทความนี้อธิบายว่า vce() ควบคุมอะไร แต่ละทางเลือกสมมติอะไร และเลือกจากการออกแบบการศึกษาอย่างไรก่อนเห็นผลใด ๆ จากนั้นจึงอธิบายวิธีรายงาน

vce() ควบคุมอะไร

คำสั่งถดถอยประมาณสองสิ่ง คือสัมประสิทธิ์ และความแปรปรวนของสัมประสิทธิ์เหล่านั้นข้ามการศึกษาซ้ำ สิ่งที่สองคือเมทริกซ์ความแปรปรวนและความแปรปรวนร่วม (variance-covariance matrix, VCE) ความแปรปรวนที่ประมาณได้ของสัมประสิทธิ์แต่ละตัวอยู่บนเส้นทแยงมุม และความแปรปรวนร่วมของสัมประสิทธิ์แต่ละคู่อยู่นอกเส้นทแยงมุม ค่าคลาดเคลื่อนมาตรฐานคือรากที่สองของเส้นทแยงมุม

ออปชัน vce() เลือกวิธีประมาณเมทริกซ์นี้ [1] มันไม่แตะต้องสัมประสิทธิ์ ค่าประมาณจุดจึงอยู่ที่เดิม สิ่งที่เปลี่ยนคือค่าคลาดเคลื่อนมาตรฐาน และเปลี่ยนช่วงความเชื่อมั่นกับ P value ทุกตัวที่สร้างจากมัน

ความแปรปรวนร่วมสำคัญเมื่อผลลัพธ์รวมสัมประสิทธิ์หลายตัว ความแปรปรวนของผลรวมของสัมประสิทธิ์สองตัว ตามที่ lincom คำนวณ รวมสองเท่าของความแปรปรวนร่วมของมันด้วย ส่วน nlcom ใช้เมทริกซ์ทั้งหมดสำหรับการรวมแบบไม่เป็นเชิงเส้น ทั้งสองจึงรับผลของ vce() ไปเต็ม ๆ

Stata: การประมาณแบบอิงแบบจำลองและ VCE ของมัน

* ---- 1. Model-based standard errors: vce(oim), the default -------------------------
poisson admissions i.program i.older, irr vce(oim) nolog
estat vce

. * ---- 1. Model-based standard errors: vce(oim), the default -------------------------
. poisson admissions i.program i.older, irr vce(oim) nolog

Poisson regression                                      Number of obs =    600
                                                        LR chi2(2)    =  68.92
                                                        Prob > chi2   = 0.0000
Log likelihood = -957.42519                             Pseudo R2     = 0.0347

------------------------------------------------------------------------------
  admissions |        IRR   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
   1.program |   .7189098   .0493347    -4.81   0.000     .6284364    .8224084
     1.older |    1.59639   .1112163     6.71   0.000     1.392637    1.829953
       _cons |   1.301232    .080132     4.28   0.000     1.153285    1.468159
------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.

. estat vce

Covariance matrix of coefficients of poisson model

             | admissions
             |          1.          1.
        e(V) |    program       older       _cons
-------------+-----------------------------------
admissions   |
   1.program |  .00470929
     1.older | -.00006111   .00485354
       _cons |  -.0019389  -.00296885   .00379229

R: การประมาณเดียวกันและเมทริกซ์เดียวกัน

# ---- 1. Model-based standard errors (the inverse information, as Stata's vce(oim)) ----
fit <- glm(admissions ~ program + older, family = poisson, data = d)
print(summary(fit)$coefficients)
b <- coef(fit)
V_oim <- vcov(fit)
print(V_oim)

> fit <- glm(admissions ~ program + older, family = poisson,
+     data = d)

> print(summary(fit)$coefficients)
              Estimate Std. Error   z value     Pr(>|z|)
(Intercept)  0.2633118 0.06158143  4.275832 1.904251e-05
program     -0.3300193 0.06862416 -4.809084 1.516236e-06
older        0.4677449 0.06966719  6.713992 1.893708e-11

> b <- coef(fit)

> V_oim <- vcov(fit)

> print(V_oim)
             (Intercept)       program         older
(Intercept)  0.003792272 -1.938898e-03 -2.968830e-03
program     -0.001938898  4.709275e-03 -6.111248e-05
older       -0.002968830 -6.111248e-05  4.853517e-03

แบบอิงแบบจำลอง: เชื่อความแปรปรวนของแบบจำลอง

การถดถอยแบบ Poisson (Poisson regression) สร้างแบบจำลองให้ข้อมูลนับ: ลอการิทึมของค่าเฉลี่ยจำนวนนับของผู้ป่วยเป็นฟังก์ชันเชิงเส้นของตัวแปรร่วม และสมมติว่าความแปรปรวนของจำนวนนับเท่ากับค่าเฉลี่ย ค่าเริ่มต้นของ Stata สำหรับแบบจำลองนี้คือ vce(oim) ซึ่งกลับค่าเมทริกซ์สารสนเทศที่สังเกตได้ (observed information matrix) อันวัดว่าลอการิทึมของความควรจะเป็นโค้งแหลมเพียงใดรอบค่าประมาณ ความโค้งนั้นคำนวณราวกับว่าความแปรปรวนแบบ Poisson ถูกต้อง ผลที่ได้คือค่าคลาดเคลื่อนมาตรฐานแบบอิงแบบจำลอง: 0.0686 สำหรับสัมประสิทธิ์ของโปรแกรม ซึ่งคือลอการิทึมของอัตราส่วนอัตราอุบัติการณ์ของมัน

ข้อมูลนับทางคลินิกมักฝ่าฝืนข้อสมมติเรื่องความแปรปรวน เมื่อจำนวนนับแปรปรวนมากกว่าค่าเฉลี่ย เรียกว่าการกระจายเกินกว่าที่แบบจำลองกำหนด (overdispersion) ค่าคลาดเคลื่อนมาตรฐานแบบอิงแบบจำลองมักเล็กเกินไป เมื่อแปรปรวนน้อยกว่า เรียกว่าการกระจายน้อยกว่าที่แบบจำลองกำหนด (underdispersion) มักใหญ่เกินไป ตอนที่ 1 แสดงทั้งสองทิศทางด้วยตัวเลข

การตรวจอย่างเร็วคือสถิติ dispersion แบบ Pearson ซึ่งเทียบกำลังสองของส่วนเหลือกับความแปรปรวนแบบ Poisson และอยู่ใกล้ 1 เมื่อความแปรปรวนนั้นเหมาะสม ในที่นี้ค่าเท่ากับ 1.41 ดังที่ผลลัพธ์ Stata ด้านล่างแสดง ความกระจายส่วนเกินมีที่มาที่ทราบในการจำลองนี้: ผู้ป่วยของโรงพยาบาลเดียวกันร่วมรับอิทธิพลของโรงพยาบาลนั้นต่อจำนวนครั้งที่นอนโรงพยาบาล

Stata: สถิติ dispersion แบบ Pearson

* Pearson dispersion: sum of (y - mu)^2/mu over n - 3 (1 when the Poisson variance is right)
predict double mu, n
generate double pearson2 = (admissions - mu)^2/mu
quietly summarize pearson2
canon pearson_dispersion 4 r(sum)/(e(N) - 3)

. * Pearson dispersion: sum of (y - mu)^2/mu over n - 3 (1 when the Poisson variance is right)
. predict double mu, n

. generate double pearson2 = (admissions - mu)^2/mu

. quietly summarize pearson2

. canon pearson_dispersion 4 r(sum)/(e(N) - 3)
CANON g2.pearson_dispersion 1.4052

แบบ robust: ผู้ป่วยเป็นอิสระต่อกัน แต่ความแปรปรวนยังเป็นที่สงสัย

vce(robust) เก็บ bread จากแบบจำลองไว้ แต่เอา meat มาจากข้อมูล score contribution ของผู้ป่วยแต่ละคนเข้าสู่ meat เป็นกำลังสอง $K = \sum_i s_i s_i^{\top}$ และความแปรปรวนคือ sandwich $J^{-1} K J^{-1}$ ในตัวอย่างขนาดใหญ่ วิธีนี้ยังใช้ได้เมื่อข้อสมมติเรื่องความแปรปรวนผิด ตราบใดที่แบบจำลองค่าเฉลี่ยถูกต้องและผู้ป่วยเป็นอิสระต่อกัน [2]

ตอนที่ 1 คำนวณ bread และ meat เป็นตัวเลข รวมถึงตัวคูณ $n/(n - 1)$ ที่ Stata ใช้ โดย $n$ คือจำนวนหน่วยสังเกต รูปแบบเมทริกซ์และตัวแปรผัน HC0 ถึง HC3 อยู่ในบทความเรื่องrobust standard error

ในที่นี้ robust standard error ของลอการิทึมของอัตราส่วนอัตราอุบัติการณ์คือ 0.0816 เป็น 1.19 เท่าของค่าแบบอิงแบบจำลอง ช่วงความเชื่อมั่น 95% ของอัตราส่วนอัตราอุบัติการณ์กว้างขึ้นเล็กน้อย เป็น 0.61 ถึง 0.84 จาก 0.63 ถึง 0.82 วิธีนี้ยังถือว่าผู้ป่วย 600 คนเป็นอิสระต่อกัน ซึ่งข้อมูลเหล่านี้ไม่สนับสนุน

แบบ cluster-robust: โรงพยาบาลเป็นอิสระต่อกัน แต่ผู้ป่วยสัมพันธ์กัน

ผู้ป่วยของโรงพยาบาลเดียวกันใช้บุคลากร แนวปฏิบัติตอนจำหน่าย และพื้นที่รับผู้ป่วยร่วมกัน จำนวนครั้งที่นอนโรงพยาบาลของพวกเขาจึงมีแนวโน้มขึ้นลงไปด้วยกัน สหสัมพันธ์ภายในกลุ่ม (intraclass correlation, ICC) วัดสิ่งนี้ คือสหสัมพันธ์ระหว่างผลลัพธ์ของผู้ป่วยสองคนในกลุ่ม (cluster) เดียวกัน บทความเรื่องการทดลองแบบสุ่มเป็นกลุ่ม (cluster randomized trials) แสดงว่ามันกำหนดการออกแบบการศึกษาอย่างไร cluster-robust standard error ยอมให้มีสหสัมพันธ์แบบใดก็ได้เช่นนี้ โดยเปลี่ยนเฉพาะ meat

ขั้นแรก score contribution ของผู้ป่วยทุกคนในกลุ่ม $g$ ซึ่งในที่นี้คือโรงพยาบาลหนึ่งแห่ง ถูกบวกเป็นคะแนนของกลุ่มหนึ่งค่า $S_g$

จากนั้นสร้าง meat จากคะแนนของกลุ่ม $G$ ค่า หนึ่งค่าต่อโรงพยาบาล โดย $G$ คือจำนวนกลุ่ม ความแปรปรวนกลายเป็น $J^{-1} K_{\text{cl}} J^{-1}$ [3] การยกกำลังสองผลรวมเก็บผลคูณของคะแนนทุกคู่ที่มาจากโรงพยาบาลเดียวกันไว้ เมื่อโรงพยาบาลหนึ่งดันจำนวนนับของผู้ป่วยทุกคนขึ้นหรือลงพร้อมกัน ส่วนเหลือของพวกเขามีเครื่องหมายเดียวกัน และผลคูณเหล่านี้ทำให้ meat ใหญ่ขึ้น

วิธีนี้มีเงื่อนไขสองอย่าง กลุ่มต้องเป็นอิสระต่อกัน และต้องมีกลุ่มจำนวนมากพอ เพราะ meat ตอนนี้อาศัยคะแนนของกลุ่ม $G$ ค่า ไม่ใช่ผู้ป่วย $n$ คน [3] คำสั่ง poisson ของ Stata คูณความแปรปรวนนี้ด้วย $G/(G - 1)$ และโค้ด R ขอให้ vcovCL() ใช้ตัวคูณเดียวกัน ผลของทั้งสองภาษาด้านล่างจึงตรงกัน

สมการประมาณค่าวางนัยทั่วไป (generalized estimating equations, GEE) ใช้ meat ที่รวมตามกลุ่มแบบเดียวกัน รอบแบบจำลองที่สมมติรูปแบบสหสัมพันธ์อย่างง่ายระหว่างผู้ป่วยในกลุ่มเดียวกัน (เรียกว่าสหสัมพันธ์สำหรับทำงาน หรือ working correlation เป็นการคาดเดารูปแบบนั้นที่ไม่จำเป็นต้องถูกต้องทุกประการ) [3] คู่มือการวัดซ้ำอธิบายเรื่องนี้

ตัวอย่างที่ลงมือคำนวณ: หนึ่งโปรแกรม ค่าคลาดเคลื่อนมาตรฐานสามค่า

ข้อมูลจำลองนี้สร้างขึ้นโดยรู้ค่าจริง มีผู้ป่วย 600 คนใน 40 โรงพยาบาล โรงพยาบาลละ 15 คน โรงพยาบาลหมายเลขคู่มีโปรแกรม ผู้ป่วย 300 คนจึงอยู่ในโรงพยาบาลที่มีโปรแกรม ตัวแปรแบบใช่หรือไม่ใช่ระบุผู้ป่วย 302 คนที่อยู่ในกลุ่มอายุมาก

จำนวนครั้งที่ผู้ป่วยแต่ละคนนอนโรงพยาบาลตลอดหนึ่งปีเป็นไปตามการแจกแจงแบบ Poisson โรงพยาบาลแต่ละแห่งมีตัวคูณของตัวเอง $e^{u}$ โดย $u$ สุ่มจากการแจกแจงปกติที่มีค่าเฉลี่ย 0 และส่วนเบี่ยงเบนมาตรฐาน 0.4 ตัวคูณนี้คูณค่าเฉลี่ยของผู้ป่วยทุกคนในโรงพยาบาลนั้น ผู้ป่วยของโรงพยาบาลเดียวกันจึงคล้ายกัน

ในโรงพยาบาลที่ไม่มีโปรแกรมและมีตัวคูณเท่ากับ 1 ค่าเฉลี่ยของผู้ป่วยอายุน้อยคือ 1.2 ค่านี้ถูกคูณด้วย 0.75 เมื่อมีโปรแกรม และคูณด้วย 1.6 สำหรับผู้ป่วยอายุมาก เมื่อเฉลี่ยข้ามโรงพยาบาล ค่าเฉลี่ยทุกค่าสูงขึ้นเล็กน้อย เพราะตัวคูณมีค่าเฉลี่ยสูงกว่า 1 เล็กน้อย

ผู้ป่วยนอนโรงพยาบาลรวม 872 ครั้ง เฉลี่ย 1.45 ครั้งต่อคน โดยมีความแปรปรวน 2.34 ค่าเฉลี่ยคือ 1.22 ในโรงพยาบาลที่มีโปรแกรม และ 1.69 ในโรงพยาบาลอื่น ทุกแถวของตารางด้านล่างประมาณแบบจำลองเดียวกัน คือการถดถอยแบบ Poisson ของจำนวนครั้งที่นอนโรงพยาบาลบนโปรแกรมและกลุ่มอายุ มีเพียง vce() ที่เปลี่ยน

โปรแกรมภายใต้ตัวประมาณความแปรปรวนสามแบบ

อ่านตาราง

อัตราส่วนอัตราอุบัติการณ์เท่ากับ 0.72 ในทุกแถว เพราะ vce() ไม่เคยขยับสัมประสิทธิ์ cluster-robust standard error เป็น 2.11 เท่าของ robust และ 2.51 เท่าของแบบอิงแบบจำลอง ช่วงความเชื่อมั่นของมัน 0.51 ถึง 1.01 ครอบคลุม 1 ขณะที่อีกสองวิธีไม่ครอบคลุม

เหตุผลอยู่ที่วิธีจัดสรรโปรแกรม ผู้ป่วยทุกคนในโรงพยาบาลหนึ่งมีสถานะโปรแกรมและตัวคูณของโรงพยาบาลเดียวกัน การเปรียบเทียบจึงอาศัยโรงพยาบาล 40 แห่ง ไม่ใช่ผู้ป่วยที่เป็นอิสระต่อกัน 600 คน มีเพียง cluster-robust standard error ที่นับข้อมูลแบบนี้

กลุ่มอายุให้ผลต่างออกไป: แต่ละโรงพยาบาลมีผู้ป่วยอายุมากและอายุน้อยในสัดส่วนใกล้เคียงกัน ตัวคูณร่วมของโรงพยาบาลยกหรือลดทั้งสองกลุ่มพร้อมกัน จึงตัดกันไปเป็นส่วนใหญ่ในการเปรียบเทียบระหว่างกลุ่มอายุ cluster-robust standard error ของกลุ่มอายุ 0.0697 จึงไม่ใหญ่กว่า robust standard error 0.0818 ในตัวอย่างนี้เป็น 0.85 เท่าของมัน การจัดกลุ่มมีราคาสูงที่สุดสำหรับตัวแปรร่วมที่เหมือนกันสำหรับทุกคนในกลุ่ม

เนื่องจากแบบจำลองไม่มีพจน์ของโรงพยาบาล 0.72 จึงประมาณอัตราส่วนที่เฉลี่ยข้ามโรงพยาบาล เมื่อใช้ log link และตัวคูณของโรงพยาบาลไม่เกี่ยวกับโปรแกรม อัตราส่วนที่เฉลี่ยแล้วนี้เท่ากับอัตราส่วนภายในโรงพยาบาล ข้อมูลถูกสร้างโดยกำหนดอัตราส่วนนั้นที่ 0.75 มีเพียงค่าคงที่ที่ดูดซับตัวคูณของโรงพยาบาล ซึ่งเป็นเหตุให้ผลลัพธ์ Stata ข้างต้นแสดงค่าคงที่สูงกว่า 1.2

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

ในตัวอย่างอื่น ค่าประมาณขยับไปราวกับที่ cluster-robust standard error บอกไว้ ช่วงที่แคบกว่าสองช่วงจะพลาด 0.75 บ่อยกว่าระดับ 95% ที่ระบุไว้มาก อัตราครอบคลุม (coverage) คือความถี่ที่ช่วงเช่นนี้ครอบคลุมค่าจริง ซึ่งตัดสินจากตัวอย่างจำนวนมาก ดังที่บทความเรื่อง robust standard error ทำด้วยการจำลอง

Stata: vce(robust)

* ---- 3. vce(robust): sandwich, each patient independent ----------------------------
poisson admissions i.program i.older, irr vce(robust) nolog

. * ---- 3. vce(robust): sandwich, each patient independent ----------------------------
. poisson admissions i.program i.older, irr vce(robust) nolog

Poisson regression                                      Number of obs =    600
                                                        Wald chi2(2)  =  41.34
                                                        Prob > chi2   = 0.0000
Log pseudolikelihood = -957.42519                       Pseudo R2     = 0.0347

------------------------------------------------------------------------------
             |               Robust
  admissions |        IRR   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
   1.program |   .7189098   .0586611    -4.04   0.000     .6126588    .8435876
     1.older |    1.59639   .1306192     5.72   0.000     1.359854     1.87407
       _cons |   1.301232   .0873629     3.92   0.000     1.140792    1.484237
------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.

Stata: vce(cluster hospital)

* ---- 4. vce(cluster hospital): sandwich, hospitals independent ----------------------
poisson admissions i.program i.older, irr vce(cluster hospital) nolog

. * ---- 4. vce(cluster hospital): sandwich, hospitals independent ----------------------
. poisson admissions i.program i.older, irr vce(cluster hospital) nolog

Poisson regression                                      Number of obs =    600
                                                        Wald chi2(2)  =  45.15
                                                        Prob > chi2   = 0.0000
Log pseudolikelihood = -957.42519                       Pseudo R2     = 0.0347

                              (Std. err. adjusted for 40 clusters in hospital)
------------------------------------------------------------------------------
             |               Robust
  admissions |        IRR   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
   1.program |   .7189098   .1237477    -1.92   0.055     .5130444    1.007381
     1.older |    1.59639   .1112673     6.71   0.000      1.39255    1.830068
       _cons |   1.301232    .156567     2.19   0.029     1.027866    1.647302
------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.

R: robust และ cluster-robust standard error

# ---- 3. Robust: sandwich, each patient independent (HC0 times n/(n - 1), as Stata's vce(robust)) ----
V_rob <- sandwich(fit) * n / (n - 1)
se_rob <- sqrt(diag(V_rob))
canon("robust.se_program", se_rob[["program"]])
irr_ci("robust.irr_program", b[["program"]], se_rob[["program"]])
canon("robust.se_older", se_rob[["older"]])

# ---- 4. Cluster-robust: hospitals independent (HC0 with the G/(G - 1) factor, as Stata's vce(cluster)) ----
V_cl <- vcovCL(fit, cluster = ~hospital, type = "HC0", cadjust = TRUE)
se_cl <- sqrt(diag(V_cl))
canon("cluster.n_clusters", length(unique(d$hospital)), 0)
canon("cluster.irr_program", exp(b[["program"]]))
canon("cluster.se_program", se_cl[["program"]])
irr_ci("cluster.irr_program", b[["program"]], se_cl[["program"]])
canon("cluster.se_older", se_cl[["older"]])
irr_ci("cluster.irr_older", b[["older"]], se_cl[["older"]])
canon("ratio.cluster_robust_program", se_cl[["program"]] / se_rob[["program"]])
canon("ratio.cluster_oim_program", se_cl[["program"]] / se_oim[["program"]])
canon("ratio.robust_oim_program", se_rob[["program"]] / se_oim[["program"]])
canon("ratio.cluster_robust_older", se_cl[["older"]] / se_rob[["older"]])

> V_rob <- sandwich(fit) * n/(n - 1)

> se_rob <- sqrt(diag(V_rob))

> canon("robust.se_program", se_rob[["program"]])
CANON g2.robust.se_program 0.0816

> irr_ci("robust.irr_program", b[["program"]], se_rob[["program"]])
CANON g2.robust.irr_program_lb 0.6127
CANON g2.robust.irr_program_ub 0.8436

> canon("robust.se_older", se_rob[["older"]])
CANON g2.robust.se_older 0.0818

> V_cl <- vcovCL(fit, cluster = ~hospital, type = "HC0",
+     cadjust = TRUE)

> se_cl <- sqrt(diag(V_cl))

> canon("cluster.n_clusters", length(unique(d$hospital)),
+     0)
CANON g2.cluster.n_clusters 40

> canon("cluster.irr_program", exp(b[["program"]]))
CANON g2.cluster.irr_program 0.7189

> canon("cluster.se_program", se_cl[["program"]])
CANON g2.cluster.se_program 0.1721

> irr_ci("cluster.irr_program", b[["program"]], se_cl[["program"]])
CANON g2.cluster.irr_program_lb 0.5130
CANON g2.cluster.irr_program_ub 1.0074

> canon("cluster.se_older", se_cl[["older"]])
CANON g2.cluster.se_older 0.0697

> irr_ci("cluster.irr_older", b[["older"]], se_cl[["older"]])
CANON g2.cluster.irr_older_lb 1.3926
CANON g2.cluster.irr_older_ub 1.8301

> canon("ratio.cluster_robust_program", se_cl[["program"]]/se_rob[["program"]])
CANON g2.ratio.cluster_robust_program 2.1095

> canon("ratio.cluster_oim_program", se_cl[["program"]]/se_oim[["program"]])
CANON g2.ratio.cluster_oim_program 2.5083

> canon("ratio.robust_oim_program", se_rob[["program"]]/se_oim[["program"]])
CANON g2.ratio.robust_oim_program 1.1890

> canon("ratio.cluster_robust_older", se_cl[["older"]]/se_rob[["older"]])
CANON g2.ratio.cluster_robust_older 0.8518

Poisson แบบ robust หรือ negative binomial: เป้าหมายสองแบบ

negative binomial regression เก็บแบบจำลองค่าเฉลี่ยของ Poisson ไว้ แล้วเพิ่มพารามิเตอร์ความกระจาย $\alpha$ จำนวนนับที่มีค่าเฉลี่ย $\mu$ จึงมีความแปรปรวน $\mu + \alpha\mu^2$ [4] มันสร้างแบบจำลองให้การแจกแจงทั้งหมดของจำนวนนับ ส่วน Poisson แบบ robust อ้างน้อยกว่า คือสร้างแบบจำลองเฉพาะค่าเฉลี่ย แล้วปล่อยให้ sandwich จัดการความแปรปรวน บทความเรื่องอัตราแบบ Poissonอธิบายแบบจำลองนี้เอง

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

ทั้งสองยังถ่วงน้ำหนักผู้ป่วยต่างกัน negative binomial ถ่วงส่วนเหลือของผู้ป่วยแต่ละคนด้วย $1/(1 + \alpha\mu)$ ผู้ป่วยที่คาดว่ามีจำนวนนับสูงจึงมีน้ำหนักน้อยลง ส่วน Poisson ถ่วงส่วนเหลือทุกตัวเท่ากัน

หากแบบจำลองค่าเฉลี่ยผิด เช่นโปรแกรมได้ผลดีกว่าในโรงพยาบาลที่ยุ่งกว่า ทั้งสองอาจรายงานอัตราส่วนที่ต่างกัน แต่ละค่าเป็นค่าเฉลี่ยถ่วงน้ำหนักคนละแบบของผลที่แปรผัน และไม่มีการเลือก vce() ใดทำให้ทั้งสองตรงกัน หัวข้อสิ่งที่ vce() เปลี่ยนไม่ได้ด้านล่างจะกลับมาที่ประเด็นนี้อีกครั้ง

ในข้อมูลจำลอง nbreg ให้อัตราส่วนอัตราอุบัติการณ์ 0.72 โดยมีช่วงความเชื่อมั่น 95% ที่ 0.62 ถึง 0.84 ต่างจากอัตราส่วนแบบ Poisson เพียงทศนิยมตำแหน่งที่สาม เพราะน้ำหนักต่างกัน $\alpha$ ที่ประมาณได้คือ 0.23 และสถิติอัตราส่วนความควรจะเป็นสำหรับ $\alpha = 0$ คือ 34.41 จำนวนนับจึงกระจายเกินกว่าที่ Poisson กำหนด ค่าคลาดเคลื่อนมาตรฐานของลอการิทึมของอัตราส่วนอัตราอุบัติการณ์ 0.0795 ใกล้เคียงค่าของ Poisson แบบ robust ที่ 0.0816

ทั้งสองถือว่าผู้ป่วยเป็นอิสระต่อกัน ในที่นี้ความกระจายเกินและการจัดกลุ่มมีที่มาเดียวกัน คือตัวคูณของโรงพยาบาล และ negative binomial ดูดซับความแปรปรวนส่วนเกินโดยไม่รู้ว่าใครใช้ร่วมกัน เมื่อเพิ่ม vce(cluster hospital) ใน Stata ค่าคลาดเคลื่อนมาตรฐานของมันสูงขึ้นเป็น 0.1702 และช่วงความเชื่อมั่นเป็น 0.52 ถึง 1.01 ผลลัพธ์ Stata ชุดสุดท้ายในหัวข้อนี้แสดงการประมาณนั้น

Stata: negative binomial regression

* ---- 5. Negative binomial: another model for the variance, still one patient at a time ----
nbreg admissions i.program i.older, irr nolog

. * ---- 5. Negative binomial: another model for the variance, still one patient at a time ----
. nbreg admissions i.program i.older, irr nolog

Negative binomial regression                            Number of obs =    600
                                                        LR chi2(2)    =  48.68
Dispersion: mean                                        Prob > chi2   = 0.0000
Log likelihood = -940.21983                             Pseudo R2     = 0.0252

------------------------------------------------------------------------------
  admissions |        IRR   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
   1.program |   .7226321   .0574734    -4.08   0.000     .6183271    .8445322
     1.older |   1.590514   .1276263     5.78   0.000      1.35905      1.8614
       _cons |    1.30083   .0918756     3.72   0.000     1.132665    1.493962
-------------+----------------------------------------------------------------
    /lnalpha |  -1.455862   .2286047                     -1.903919   -1.007805
-------------+----------------------------------------------------------------
       alpha |   .2331992   .0533104                      .1489836    .3650192
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation to incidence-rate ratios.
Note: _cons estimates baseline incidence rate.
LR test of alpha=0: chibar2(01) = 34.41                Prob >= chibar2 = 0.000

R: glm.nb จากแพ็กเกจ MASS

# ---- 5. Negative binomial: another model for the variance, still one patient at a time ----
nb <- glm.nb(admissions ~ program + older, data = d)
print(summary(nb))

> nb <- glm.nb(admissions ~ program + older, data = d)

> print(summary(nb))

Call:
glm.nb(formula = admissions ~ program + older, data = d, init.theta = 4.288179278,
    link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.26300    0.07037   3.737 0.000186 ***
program     -0.32486    0.07953  -4.085 4.41e-05 ***
older        0.46406    0.08024   5.783 7.32e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(4.2882) family taken to be 1)

    Null deviance: 703.84  on 599  degrees of freedom
Residual deviance: 652.94  on 597  degrees of freedom
AIC: 1888.4

Number of Fisher Scoring iterations: 1


              Theta:  4.288
          Std. Err.:  0.980

 2 x log-likelihood:  -1880.440

Stata: negative binomial ร่วมกับ vce(cluster hospital)

* negative binomial with a cluster-robust standard error (Stata only)
nbreg admissions i.program i.older, irr vce(cluster hospital) nolog
canon nb_cluster.se_program.stata 4 _se[1.program]

. * negative binomial with a cluster-robust standard error (Stata only)
. nbreg admissions i.program i.older, irr vce(cluster hospital) nolog

Negative binomial regression                            Number of obs =    600
                                                        Wald chi2(2)  =  44.64
Dispersion: mean                                        Prob > chi2   = 0.0000
Log pseudolikelihood = -940.21983                       Pseudo R2     = 0.0252

                              (Std. err. adjusted for 40 clusters in hospital)
------------------------------------------------------------------------------
             |               Robust
  admissions |        IRR   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
   1.program |   .7226321   .1229716    -1.91   0.056     .5176864    1.008713
     1.older |   1.590514   .1109611     6.65   0.000     1.387248    1.823564
       _cons |    1.30083   .1573644     2.17   0.030     1.026239    1.648892
-------------+----------------------------------------------------------------
    /lnalpha |  -1.455862   .3898306                     -2.219916   -.6918084
-------------+----------------------------------------------------------------
       alpha |   .2331992   .0909082                      .1086182    .5006699
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation to incidence-rate ratios.
Note: _cons estimates baseline incidence rate.

. canon nb_cluster.se_program.stata 4 _se[1.program]
CANON g2.nb_cluster.se_program.stata 0.1702

การเลือก vce(): เริ่มจากการออกแบบ

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

OPG การสุ่มซ้ำ และกรณีที่มีกลุ่มน้อย

มีออปชันหนึ่งที่ดูเหมือน sandwich แต่ไม่ใช่ vce(opg) หรือ outer product of the gradients กลับค่า meat อย่างเดียว เหมือน vce(oim) จึงใช้ได้เฉพาะเมื่อความแปรปรวนของแบบจำลองถูกต้อง ในที่นี้ให้ 0.0595 ซึ่งเล็กกว่าค่าแบบอิงแบบจำลอง 0.0686 เสียอีก ผลลัพธ์ Stata ท้ายหัวข้อนี้แสดงการประมาณนั้น

bootstrap และ jackknife ประมาณแบบจำลองใหม่บนข้อมูลที่สุ่มซ้ำ [1] ทั้งสองคำนึงถึงการจัดกลุ่มก็ต่อเมื่อสุ่มโรงพยาบาลซ้ำทั้งแห่ง ซึ่งใน Stata เกิดขึ้นเมื่อให้ cluster(hospital) ภายใน vce() เท่านั้น [1]

cluster-robust standard error พึ่งพากลุ่มจำนวนมาก เมื่อมีกลุ่มน้อย ค่านี้มักเล็กเกินไป และการแก้สำหรับตัวอย่างขนาดเล็ก โดยตัดสินนัยสำคัญกับการแจกแจง t อาจเหมาะกว่า คู่มือการวัดซ้ำอธิบายไว้ ควรรายงานจำนวนกลุ่มคู่กับผลแบบ cluster-robust ทุกครั้ง

โรงพยาบาล 40 แห่งในที่นี้ ซึ่งครึ่งหนึ่งมีโปรแกรม มักถือว่าเพียงพอสำหรับ cluster-robust standard error แบบธรรมดา การแก้สำหรับตัวอย่างขนาดเล็กจะทำให้ช่วงกว้างขึ้นบ้าง และอาจรายงานเป็นการวิเคราะห์ความไว (sensitivity analysis)

Stata: vce(opg)

* ---- 2. vce(opg): the outer product of the scores alone ---------------------------
poisson admissions i.program i.older, irr vce(opg) nolog
canon opg.se_program 4 _se[1.program]

. * ---- 2. vce(opg): the outer product of the scores alone ---------------------------
. poisson admissions i.program i.older, irr vce(opg) nolog

Poisson regression                                      Number of obs =    600
                                                        LR chi2(2)    =  68.92
                                                        Prob > chi2   = 0.0000
Log likelihood = -957.42519                             Pseudo R2     = 0.0347

------------------------------------------------------------------------------
             |                 OPG
  admissions |        IRR   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
   1.program |   .7189098   .0427671    -5.55   0.000     .6397901     .807814
     1.older |    1.59639   .0975423     7.66   0.000     1.416215    1.799488
       _cons |   1.301232   .0741003     4.62   0.000      1.16381    1.454881
------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.

. canon opg.se_program 4 _se[1.program]
CANON g2.opg.se_program 0.0595

สิ่งที่ vce() เปลี่ยนไม่ได้

ทุกวิธีข้างต้นเปลี่ยนว่าค่าประมาณดูไม่แน่นอนเพียงใด ไม่มีวิธีใดเปลี่ยนสิ่งที่ถูกประมาณ ปริมาณเป้าหมายของการประมาณ (estimand) คือปริมาณที่การวิเคราะห์ตั้งใจประมาณ ถูกกำหนดโดยคำถาม ได้แก่ประชากร ผลลัพธ์ การเปรียบเทียบ และตัวสรุป ในที่นี้คืออัตราส่วนของจำนวนครั้งที่นอนโรงพยาบาลเฉลี่ยตลอดหนึ่งปีในโรงพยาบาลที่มีโปรแกรมต่อค่านั้นในโรงพยาบาลอื่น เปรียบเทียบผู้ป่วยกลุ่มอายุเดียวกัน สมมติอัตราส่วนเดียวสำหรับทั้งสองกลุ่มอายุ และเฉลี่ยข้ามโรงพยาบาล: อัตราส่วนอัตราอุบัติการณ์ระดับประชากร (marginal) หรือเฉลี่ยทั้งประชากร (population-averaged)

robust standard error ช่วยแก้ค่าคลาดเคลื่อนมาตรฐานเมื่อข้อสมมติเรื่องความแปรปรวนผิด แต่ไม่ได้แก้แบบจำลองค่าเฉลี่ยที่ผิด และสัมประสิทธิ์ที่ลำเอียงก็ยังคงลำเอียงอยู่

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

ตัวประมาณความแปรปรวนก็ไม่เปลี่ยนสิ่งที่ตัวอย่างหนึ่งประมาณได้ การเพิ่ม vce(robust) ให้แบบจำลอง Poisson ดัดแปลงที่ประมาณกับตัวอย่างแบบผู้ป่วย-กลุ่มควบคุม (case-control) ไม่ได้เปลี่ยนค่าประมาณให้เป็นอัตราส่วนความเสี่ยง วิธีสุ่มกลุ่มควบคุมเป็นตัวกำหนดมาตรวัด ดังที่บทความเรื่องการสุ่มกลุ่มควบคุมแสดง

การรายงาน: ประโยคสำหรับส่วนวิธีการและส่วนผลลัพธ์

ส่วนวิธีการควรระบุตัวประมาณความแปรปรวน เหตุผลที่ใช้ และจำนวนกลุ่ม ส่วนผลลัพธ์ควรบอกว่ารายงานช่วงความเชื่อมั่นแบบใด

สำหรับผลลัพธ์แบบใช่หรือไม่ใช่ คือการกลับมานอนโรงพยาบาลซ้ำภายในหนึ่งปี อัตราส่วนความเสี่ยงมาจาก Poisson ดัดแปลง (modified Poisson regression) คือแบบจำลอง Poisson ที่ใช้ log link ประมาณกับผลลัพธ์ 0 หรือ 1 ร่วมกับความแปรปรวนแบบ sandwich [5] ต้องใช้ sandwich เพราะความแปรปรวนแบบ Poisson ซึ่งคือความเสี่ยง $p$ เอง มากกว่าความแปรปรวนแบบ Bernoulli (Bernoulli variance) ของผลลัพธ์แบบใช่หรือไม่ใช่ ซึ่งคือ $p(1 - p)$ Zou และ Donner ขยายวิธีนี้ไปยังข้อมูลสองค่าที่จับกลุ่มกัน โดยใช้ความแปรปรวนแบบ cluster-robust [6] ขั้นตอนที่นำไปสู่แบบจำลองนี้อยู่ในบทความเรื่องอัตราส่วนความเสี่ยงสำหรับผลลัพธ์ที่พบบ่อย

ในข้อมูลจำลอง ความเสี่ยงของการกลับมานอนโรงพยาบาลซ้ำ ซึ่งพิมพ์ไว้ในผลลัพธ์ Stata ชุดแรกด้านล่าง คือ 0.667 ในโรงพยาบาลที่มีโปรแกรม และ 0.763 ในโรงพยาบาลอื่น หลังปรับกลุ่มอายุ อัตราส่วนความเสี่ยงคือ 0.87 โดยมีช่วงความเชื่อมั่นแบบ robust ที่ 0.79 ถึง 0.96 ช่วงความเชื่อมั่นแบบ cluster-robust ซึ่งเป็นช่วงที่ควรรายงานเพราะโปรแกรมถูกจัดสรรตามโรงพยาบาล คือ 0.75 ถึง 1.01 อัตราส่วนความเสี่ยงนี้เฉลี่ยข้ามโรงพยาบาลด้วย และสมมติอัตราส่วนเดียวสำหรับทั้งสองกลุ่มอายุ การจำลองไม่ได้กำหนดอัตราส่วนความเสี่ยงที่แท้จริงไว้

Stata: ความเสี่ยงของการกลับมานอนโรงพยาบาลซ้ำในแต่ละกลุ่ม

quietly summarize any_admission if program == 1
canon any.risk_program 4 r(mean)
quietly summarize any_admission if program == 0
canon any.risk_control 4 r(mean)

. quietly summarize any_admission if program == 1

. canon any.risk_program 4 r(mean)
CANON g2.any.risk_program 0.6667

. quietly summarize any_admission if program == 0

. canon any.risk_control 4 r(mean)
CANON g2.any.risk_control 0.7633

Stata: Poisson ดัดแปลง แบบ robust และ cluster-robust

* ---- 6. Binary outcome (any admission): modified Poisson for the risk ratio ------------
poisson any_admission i.program i.older, irr vce(robust) nolog
canon any.rr_program 4 exp(_b[1.program])
canon any.se_robust 4 _se[1.program]
canon any.rr_lb_robust 4 exp(_b[1.program] - invnormal(0.975)*_se[1.program])
canon any.rr_ub_robust 4 exp(_b[1.program] + invnormal(0.975)*_se[1.program])
poisson any_admission i.program i.older, irr vce(cluster hospital) nolog

. * ---- 6. Binary outcome (any admission): modified Poisson for the risk ratio ------------
. poisson any_admission i.program i.older, irr vce(robust) nolog

Poisson regression                                      Number of obs =    600
                                                        Wald chi2(2)  =  23.03
                                                        Prob > chi2   = 0.0000
Log pseudolikelihood = -569.62054                       Pseudo R2     = 0.0058

-------------------------------------------------------------------------------
              |               Robust
any_admission |        IRR   std. err.      z    P>|z|     [95% conf. interval]
--------------+----------------------------------------------------------------
    1.program |   .8709482   .0446817    -2.69   0.007     .7876326     .963077
      1.older |   1.231833   .0645076     3.98   0.000     1.111672    1.364982
        _cons |   .6845157   .0322619    -8.04   0.000     .6241162    .7507604
-------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.

. canon any.rr_program 4 exp(_b[1.program])
CANON g2.any.rr_program 0.8709

. canon any.se_robust 4 _se[1.program]
CANON g2.any.se_robust 0.0513

. canon any.rr_lb_robust 4 exp(_b[1.program] - invnormal(0.975)*_se[1.program])
CANON g2.any.rr_lb_robust 0.7876

. canon any.rr_ub_robust 4 exp(_b[1.program] + invnormal(0.975)*_se[1.program])
CANON g2.any.rr_ub_robust 0.9631

. poisson any_admission i.program i.older, irr vce(cluster hospital) nolog

Poisson regression                                      Number of obs =    600
                                                        Wald chi2(2)  =  25.91
                                                        Prob > chi2   = 0.0000
Log pseudolikelihood = -569.62054                       Pseudo R2     = 0.0058

                               (Std. err. adjusted for 40 clusters in hospital)
-------------------------------------------------------------------------------
              |               Robust
any_admission |        IRR   std. err.      z    P>|z|     [95% conf. interval]
--------------+----------------------------------------------------------------
    1.program |   .8709482   .0657274    -1.83   0.067     .7511992    1.009787
      1.older |   1.231833   .0558421     4.60   0.000     1.127106    1.346291
        _cons |   .6845157   .0400294    -6.48   0.000     .6103887    .7676449
-------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.

ประโยคสำหรับนำไปปรับใช้

ประโยคในส่วนวิธีการแต่ละประโยคระบุแบบจำลอง มาตรวัด ตัวประมาณความแปรปรวน เหตุผล และหน่วยของกลุ่ม ประโยคในส่วนผลลัพธ์แต่ละประโยคระบุช่วงความเชื่อมั่นที่รายงาน

ส่วนวิธีการ ผลลัพธ์แบบใช่หรือไม่ใช่: ประมาณอัตราส่วนความเสี่ยงของการกลับมานอนโรงพยาบาลซ้ำภายในหนึ่งปีด้วย Poisson ดัดแปลงที่ปรับกลุ่มอายุ โดยใช้ cluster-robust standard error และถือโรงพยาบาล 40 แห่งเป็นกลุ่ม เพราะโปรแกรมถูกจัดสรรตามโรงพยาบาล

ส่วนผลลัพธ์ ผลลัพธ์แบบใช่หรือไม่ใช่: พบการกลับมานอนโรงพยาบาลซ้ำในผู้ป่วย 66.7% ของโรงพยาบาลที่มีโปรแกรม และ 76.3% ของโรงพยาบาลอื่น (อัตราส่วนความเสี่ยงที่ปรับแล้ว 0.87, ช่วงความเชื่อมั่น 95% 0.75 ถึง 1.01 จาก cluster-robust standard error)

ส่วนวิธีการ ผลลัพธ์แบบนับ: ประมาณอัตราส่วนอัตราอุบัติการณ์ของจำนวนครั้งที่นอนโรงพยาบาลในหนึ่งปีหลังจำหน่ายด้วยการถดถอยแบบ Poisson ที่ปรับกลุ่มอายุ โดยใช้ cluster-robust standard error และถือโรงพยาบาล 40 แห่งเป็นกลุ่ม เพราะโปรแกรมถูกจัดสรรตามโรงพยาบาล negative binomial regression ที่จัดกลุ่มแบบเดียวกันเป็นการวิเคราะห์ความไว

ส่วนผลลัพธ์ ผลลัพธ์แบบนับ: ผู้ป่วยในโรงพยาบาลที่มีโปรแกรมนอนโรงพยาบาลเฉลี่ย 1.22 ครั้ง เทียบกับ 1.69 ครั้งในโรงพยาบาลอื่น (อัตราส่วนอัตราอุบัติการณ์ที่ปรับแล้ว 0.72, ช่วงความเชื่อมั่น 95% 0.51 ถึง 1.01 จาก cluster-robust standard error; negative binomial ที่จัดกลุ่มแบบเดียวกัน 0.72, ช่วงความเชื่อมั่น 95% 0.52 ถึง 1.01)

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

สิ่งที่ควรทำในการวิเคราะห์ของคุณเอง

• ระบุหน่วยที่ถูกสุ่มเลือกหรือถูกจัดสรร และเขียนการเลือก vce() พร้อมเหตุผลลงในแผนการวิเคราะห์ก่อนประมาณแบบจำลอง
• ตรวจความกระจายของข้อมูลนับ แต่ให้การออกแบบการศึกษา ไม่ใช่การตรวจนั้น เป็นตัวตัดสินว่าจะจัดกลุ่มหรือไม่
• รายงานจำนวนกลุ่มทุกครั้งที่รายงานผลแบบ cluster-robust และพิจารณาการแก้สำหรับตัวอย่างขนาดเล็กเมื่อกลุ่มมีน้อย
• ระบุว่าเป้าหมายคือความสัมพันธ์ของค่าเฉลี่ยหรือการแจกแจงของจำนวนนับ ก่อนเลือก Poisson แบบ robust หรือ negative binomial
• รัน lincom หรือ nlcom หลังประมาณแบบจำลอง สำหรับทุกผลที่รวมสัมประสิทธิ์หลายตัว เพื่อให้ VCE ที่เลือกถูกนำไปใช้ต่อ

อ่านต่อในวิกิ: [[sandwich-standard-error-worked-by-hand-th]] [[heteroskedasticity-robust-standard-errors-th]] [[repeated-measures-modeling-guide-th]] [[poisson-regression-rates-lambda]] [[cluster-randomized-trials-design-ethics]]