สัมประสิทธิ์เท่าเดิม แต่ standard error ต่างกัน: คำนวณ sandwich estimator ด้วยมือ

Read the English version

บทคัดย่อ

การถดถอยแบบ Poisson ที่ประมาณสองครั้ง ครั้งแรกแบบปกติและครั้งที่สองด้วย robust standard error พิมพ์สัมประสิทธิ์หนึ่งค่าคู่กับ standard error สองค่า สัมประสิทธิ์มาจากแบบจำลองของค่าเฉลี่ย ส่วน standard error มาจากการคำนวณความแปรปรวนแยกต่างหาก ชุดข้อมูลตัวอย่างสองชุดที่คำนวณด้วยมือ ชุดละผู้ป่วยสี่คน ซึ่งน้อยเกินกว่าจะใช้อนุมานได้ มีค่าเฉลี่ยเท่ากันที่ 5 ครั้งของการนอนโรงพยาบาล แบบจำลอง Poisson ที่มีเฉพาะจุดตัดแกน (intercept-only) ให้สัมประสิทธิ์ ln 5 = 1.609 และ standard error แบบอิงแบบจำลอง 0.224 เท่ากันทั้งสองชุด ตัวประมาณแบบแซนด์วิช (sandwich estimator) อ่านความกระจายจากข้อมูลแทน ส่วน bread ของมันคือ J = 20 เท่ากันทั้งสองชุด ส่วน meat คือ K เท่ากับ 2 ในชุดหนึ่งและ 64 ในอีกชุดหนึ่ง standard error แบบ sandwich HC0 จึงเป็น 0.071 และ 0.400 ตัวประมาณแบบแซนด์วิชทำให้ standard error เล็กลงเมื่อจำนวนนับแปรปรวนน้อยกว่าที่สมมติไว้ และทำให้กว้างขึ้นเมื่อแปรปรวนมากกว่า vce(robust) ของ Stata ยังคูณความแปรปรวนด้วย n/(n - 1) จึงพิมพ์ค่าที่ใหญ่กว่า บทความนี้คำนวณทีละขั้นด้วยมือ ตรวจซ้ำด้วย Stata และ R และต่อไปยังผลลัพธ์แบบสองค่า

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

แบบจำลองเดียว standard error สองค่า

ทีมเภสัชกรรมของโรงพยาบาลแห่งหนึ่งตรวจสอบจำนวนผู้ป่วยรับใหม่ทั่วเครือข่ายโรงพยาบาลในภูมิภาคเพื่อวางแผนคลังยา แฟ้มข้อมูลของทีมมีจำนวนผู้ป่วยรับใหม่หนึ่งวันของโรงพยาบาล 100 แห่ง จำนวนนับแกว่งตั้งแต่ 0 ถึง 10 รอบค่าเฉลี่ย 5 รายต่อวัน โรงพยาบาลเหล่านี้เป็นของสมมติและข้อมูลเป็นข้อมูลจำลอง

นักวิเคราะห์ประมาณการถดถอยแบบ Poisson (Poisson regression) ซึ่งเป็นแบบจำลองมาตรฐานสำหรับข้อมูลจำนวนนับ และสร้างแบบจำลองให้ลอการิทึมของค่าเฉลี่ยของจำนวนนับ เมื่อไม่มีตัวแปรร่วม สัมประสิทธิ์ตัวเดียวของมันคือลอการิทึมของค่าเฉลี่ย นักวิเคราะห์รันใน Stata สองครั้ง ครั้งแรกแบบปกติ และอีกครั้งด้วย option vce(robust) ซึ่งสั่งให้ Stata ประมาณความแปรปรวนแบบ robust (sandwich)

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

สัมประสิทธิ์เดียวกันกับ standard error สองค่า

ความกระจายสองชนิด

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

ลองนึกภาพการศึกษาเดียวกันที่ทำซ้ำ 1,000 ครั้ง ทุกครั้งใช้ตัวอย่างใหม่จากประชากรเดียวกัน standard error ก็คือส่วนเบี่ยงเบนมาตรฐานของสัมประสิทธิ์ 1,000 ค่านั้นโดยแท้ คือรากที่สองของความแปรปรวนชนิดที่สอง การศึกษาเพียงครั้งเดียวไม่เคยเห็นการรันอื่น จึงต้องประมาณความกระจายนั้นจากข้อมูลของตัวเอง

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

ตัวอย่างคำนวณด้วยมือ: ข้อมูลสองชุดที่มีค่าเฉลี่ยเท่ากัน

ตัวอย่างคำนวณด้วยมือ ข้อมูลสองชุดบันทึกจำนวนครั้งการนอนโรงพยาบาลของผู้ป่วยสี่คนในช่วงเวลาเดียวกัน ชุด A มีจำนวนนับ 4, 5, 5 และ 6 ส่วนชุด B มี 1, 1, 9 และ 9 ทั้งสองชุดรวมกันได้ 20 ครั้ง ค่าเฉลี่ยจึงเป็น 5 เท่ากัน

ชุด A กระจุกตัวใกล้กัน ส่วนชุด B แยกเป็นสองคู่ที่ต่างกันมากความแปรปรวนเชิงประจักษ์ (empirical variance) วัดความกระจายนั้น คือค่าเฉลี่ยของระยะห่างกำลังสองจากค่าเฉลี่ย โดยหารด้วย $n$ ซึ่งเป็นจำนวนผู้ป่วย การหารด้วย $n$ แทน $n - 1$ ตรงกับสูตร sandwich ที่ใช้ด้านล่าง ในชุด A ระยะห่างกำลังสองจากค่าเฉลี่ยคือ 1, 0, 0 และ 1

ความแปรปรวนเชิงประจักษ์ของชุด A จึงเป็น 0.5 ในชุด B ทุกจำนวนนับอยู่ห่างจากค่าเฉลี่ย 4 ครั้ง ระยะห่างกำลังสองแต่ละค่าจึงเป็น 16

ความแปรปรวนเชิงประจักษ์ของชุด B จึงเป็น 16 สูงกว่า 0.5 ของชุด A มาก

สิ่งที่แบบจำลอง Poisson สมมติ

การถดถอยแบบ Poisson ที่มีเฉพาะจุดตัดแกน (intercept-only) คือไม่มีตัวแปรร่วม จะประมาณลอการิทึมของค่าเฉลี่ยของจำนวนนับ

ในที่นี้ $\bar y$ คือค่าเฉลี่ยของจำนวนนับ และ $\hat\beta_0$ คือจุดตัดแกนที่ประมาณได้ เครื่องหมายหมวก (hat) แสดงค่าประมาณ ข้อมูลทั้งสองชุดมีค่าเฉลี่ย 5 จึงได้สัมประสิทธิ์เดียวกัน ค่าเฉลี่ยที่แบบจำลองคาด $\hat\mu = 5$ เท่ากันสำหรับผู้ป่วยทุกคนในชุดข้อมูลหนึ่ง

แบบจำลอง Poisson ยังสมมติว่าความแปรปรวนของจำนวนนับเท่ากับค่าเฉลี่ยของมัน ในที่นี้คือ 5 ชุด A ซึ่งมีความแปรปรวนเชิงประจักษ์ 0.5 แสดงการกระจายน้อยกว่าที่แบบจำลองกำหนด (underdispersion) คือกระจายน้อยกว่าที่แบบจำลองสมมติ ชุด B ที่ 16 แสดงการกระจายเกินกว่าที่แบบจำลองกำหนด (overdispersion) คือกระจายมากกว่าที่แบบจำลองสมมติ

standard error แบบอิงแบบจำลอง

standard error แบบอิงแบบจำลอง (model-based standard error) คือค่าที่ได้จากความแปรปรวนซึ่งแบบจำลองสมมติไว้ สำหรับแบบจำลอง Poisson ที่มีเฉพาะจุดตัดแกน เมื่อมีผู้ป่วย $n$ คนและค่าเฉลี่ยที่แบบจำลองคาดคือ $\hat\mu$ ค่านี้คือ

ในที่นี้ $n\hat\mu = 4 \times 5 = 20$ เป็นจำนวนนับรวมที่คาดไว้ standard error จึงเป็น $\sqrt{1/20} = \sqrt{0.05} = 0.224$ บนมาตรวัดลอการิทึม เท่ากันทั้งสองชุดข้อมูล สูตรนี้เห็นเพียง $n$ และ $\hat\mu$ จึงแยกชุด A ออกจากชุด B ไม่ได้

แซนด์วิช: bread และ meat

ตัวประมาณแบบแซนด์วิช (sandwich estimator) คงสัมประสิทธิ์ไว้ แต่คำนวณความแปรปรวนของมันใหม่จากข้อมูล [1, 2] standard error ของมันมักเรียกว่า robust standard error รูปแบบพื้นฐานที่ไม่มีตัวคูณปรับสำหรับตัวอย่างเล็กเรียกว่า HC0 สำหรับแบบจำลอง Poisson ที่มีเฉพาะจุดตัดแกน ทุกส่วนของแซนด์วิชเป็นตัวเลขตัวเดียว

score contribution ของแต่ละหน่วยสังเกต (score contribution) ของผู้ป่วยคนที่ $i$ คือ $s_i = y_i - \hat\mu$ จำนวนนับ $y_i$ ของผู้ป่วยลบด้วยค่าเฉลี่ยที่แบบจำลองคาด เมื่อมีตัวแปรร่วม ค่านี้คือส่วนตกค้างนั้นคูณค่าตัวแปรร่วมของผู้ป่วย ส่วนส่วน meat ของ sandwich estimator (meat) คือ $K$ ซึ่งรวมค่าเหล่านี้หลังยกกำลังสอง

ส่วน meat จึงอ่านจากข้อมูล ผู้ป่วยแต่ละคนเพิ่มระยะห่างกำลังสองจากค่าเฉลี่ยที่แบบจำลองคาด ส่วนส่วน bread ของ sandwich estimator (bread) คือ $J$ วัดว่าสมการประมาณค่า (estimating equation) เปลี่ยนแรงเพียงใดเมื่อสัมประสิทธิ์เปลี่ยน สมการนั้นคือสมการที่โปรแกรมแก้เพื่อหาสัมประสิทธิ์ สำหรับแบบจำลองนี้ bread คือ

ค่านี้คือจำนวนนับรวมที่คาดไว้ ซึ่งเป็นผลคูณที่อยู่ในสูตรแบบอิงแบบจำลอง ความแปรปรวนแบบแซนด์วิชวางส่วนผกผันของ bread ไว้ทั้งสองข้างของ meat

รากที่สองของค่านี้คือ standard error แบบ HC0 ความแปรปรวนแบบอิงแบบจำลองคือส่วนผกผันของ bread เพียงอย่างเดียว $1/(n\hat\mu)$ หากจำนวนนับแปรปรวนตามที่แบบจำลอง Poisson สมมติ $K$ จะใกล้ $J$ และ standard error ทั้งสองแบบจะเกือบเท่ากัน เมื่อมีตัวแปรร่วม $J$ และ $K$ กลายเป็นเมทริกซ์ ดังที่บทความเรื่อง robust standard errorแสดงไว้

ทางลัดที่ใช้ได้กับแบบจำลองนี้เท่านั้น

สำหรับแบบจำลอง Poisson ที่มีเฉพาะจุดตัดแกน ขั้นตอนทั้งหมดยุบเหลือบรรทัดเดียว

สำหรับชุด B ได้ $\sqrt{64}/20 = 8/20 = 0.400$ ทางลัดนี้ใช้ได้กับแบบจำลองที่มีเฉพาะจุดตัดแกนเท่านั้น เพราะมีเพียงที่นั่นที่ bread เป็นตัวเลขตัวเดียว $n\hat\mu$ เมื่อมีตัวแปรร่วม standard error แต่ละค่ามาจากเส้นทแยงมุมของเมทริกซ์ $J^{-1} K J^{-1}$

แซนด์วิชขยับไปทางไหน

ในแบบจำลองที่มีเฉพาะจุดตัดแกนนี้ meat หารด้วย $n$ คือความแปรปรวนเชิงประจักษ์ ส่วน bread หารด้วย $n$ คือความแปรปรวนที่แบบจำลองสมมติ ดังนั้นที่นี่ standard error แบบ HC0 จะเกิน standard error แบบอิงแบบจำลองพอดีเมื่อ $K$ เกิน $J$

ชุด A มี $K = 2$ เทียบกับ $J = 20$ แซนด์วิชจึงทำให้ standard error ของมันเล็กลง ชุด B มี $K = 64$ เทียบกับ $J = 20$ แซนด์วิชจึงทำให้ standard error กว้างขึ้น เมื่อมีตัวแปรร่วม $J$ และ $K$ เป็นเมทริกซ์ และการเปรียบเทียบเพียงครั้งเดียวระหว่างสองค่าตัดสินทิศทางไม่ได้ ทิศทางขึ้นกับว่าส่วนตกค้างขนาดใหญ่อยู่ที่ใด ดังที่บทความเรื่อง robust standard errorอธิบายไว้

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

การคำนวณเดียวกันใน Stata และ R

สคริปต์จำลองของชุดบทความนี้พิมพ์ข้อมูลทั้งสองชุดเป็นจำนวนนับคงที่ แล้วทำซ้ำแต่ละขั้นด้วยคำสั่งเมทริกซ์ ใน Stata คำสั่ง matrix accum สร้าง $J$ โดยถ่วงผู้ป่วยแต่ละคนด้วย $\hat\mu$ และสร้าง $K$ โดยถ่วงแต่ละคนด้วยส่วนตกค้างกำลังสอง แผงด้านล่างแสดงชุด B สำหรับชุด A บรรทัดเดียวกันพิมพ์ $J = 20$ $K = 2$ และความแปรปรวนแบบแซนด์วิช 0.005

Stata: แบบจำลองแบบอิงแบบจำลอง แล้วจึงถึง bread และ meat

    * model-based standard error: the Poisson model assumes variance = mean
    poisson admissions if dataset == `d', nolog
    canon `k'.b0 4 _b[_cons]
    canon `k'.se_model 4 _se[_cons]
    * bread and meat by hand (one column, the constant)
    predict double mu if dataset == `d', n
    generate double r2 = (admissions - mu)^2 if dataset == `d'
    matrix accum J = one if dataset == `d' [iweight = mu], noconstant
    matrix accum K = one if dataset == `d' [iweight = r2], noconstant
    matrix Ji = invsym(J)
    matrix V = Ji*K*Ji
    matrix list J
    matrix list K
    matrix list V

Poisson regression                                      Number of obs =      4
                                                        LR chi2(0)    =   0.00
                                                        Prob > chi2   =      .
Log likelihood = -13.414897                             Pseudo R2     = 0.0000

------------------------------------------------------------------------------
  admissions | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   1.609438   .2236068     7.20   0.000     1.171177    2.047699
------------------------------------------------------------------------------
CANON g1.b.b0 1.6094
CANON g1.b.se_model 0.2236
(4 missing values generated)
(4 missing values generated)
(obs=20)
(obs=64)

symmetric J[1,1]
     one
one   20

symmetric K[1,1]
     one
one   64

symmetric V[1,1]
     one
one  .16

R: การประมาณเดียวกัน bread และ meat และ vcovHC()

  # model-based standard error: the Poisson model assumes variance = mean
  fit <- glm(admissions ~ 1, family = poisson, data = s)
  print(summary(fit)$coefficients)
  canon(paste0(k, ".b0"), coef(fit)[[1]])
  canon(paste0(k, ".se_model"), sqrt(vcov(fit)[1, 1]))
  # bread and meat by hand (one column, the constant)
  bm <- bread_meat(matrix(1, n, 1), y, fitted(fit))
  print(bm)
  canon(paste0(k, ".bread"), bm$J[1, 1])
  canon(paste0(k, ".meat"), bm$K[1, 1])
  canon(paste0(k, ".v_model"), bm$Ji[1, 1], 6)
  canon(paste0(k, ".v_hc0"), bm$V[1, 1], 6)
  canon(paste0(k, ".se_hc0"), sqrt(bm$V[1, 1]))
  # Stata's robust standard error: the sandwich package's HC0 times n/(n - 1)
  canon(paste0(k, ".se_robust"), sqrt(sandwich(fit)[1, 1] * n / (n - 1)))
  canon(paste0(k, ".robust_factor"), n / (n - 1))
  # the same HC0 variance from sandwich::vcovHC() (for a GLM it equals sandwich(fit) and the hand value above)
  V_hc0 <- sandwich::vcovHC(fit, type = "HC0")
  print(V_hc0)
  stopifnot(isTRUE(all.equal(V_hc0[1, 1], bm$V[1, 1])))

            Estimate Std. Error  z value     Pr(>|z|)
(Intercept) 1.609438  0.2236067 7.197629 6.126851e-13
CANON g1.b.b0 1.6094
CANON g1.b.se_model 0.2236
$J
     [,1]
[1,]   20

$K
     [,1]
[1,]   64

$Ji
     [,1]
[1,] 0.05

$V
     [,1]
[1,] 0.16

CANON g1.b.bread 20.0000
CANON g1.b.meat 64.0000
CANON g1.b.v_model 0.050000
CANON g1.b.v_hc0 0.160000
CANON g1.b.se_hc0 0.4000
CANON g1.b.se_robust 0.4619
CANON g1.b.robust_factor 1.3333
            (Intercept)
(Intercept)        0.16

ทำไม Stata พิมพ์ robust standard error ที่ใหญ่กว่า

ประมาณชุด B ด้วย vce(robust) แล้ว Stata พิมพ์ standard error 0.4619 ไม่ใช่ 0.400 option นี้สั่งให้ Stata ประมาณความแปรปรวนแบบ robust (sandwich) [3] สำหรับคำสั่งที่ประมาณด้วยวิธีภาวะน่าจะเป็นสูงสุด (maximum likelihood) เช่น poisson หรือ logistic Stata คูณความแปรปรวนนั้นด้วย $n/(n - 1)$ ค่าที่คำนวณด้วยมือข้างต้นจึงเรียกว่า HC0 คือแซนด์วิชที่ไม่มีตัวคูณนั้น

เมื่อมีผู้ป่วยสี่คน ตัวคูณนี้ใหญ่

ความแปรปรวนเพิ่มขึ้นด้วยตัวคูณ 1.333 standard error เพิ่มขึ้นด้วยรากที่สองของตัวคูณนั้น

standard error จึงเพิ่มขึ้นด้วยตัวคูณ 1.155 ชุด A กลายเป็น $0.071 \times 1.155 = 0.082$ และชุด B กลายเป็น $0.400 \times 1.155 = 0.462$ ค่าเหล่านี้คือค่าที่ Stata พิมพ์ ปัดเป็นทศนิยมสามตำแหน่ง การรันจึงยืนยันตัวคูณนี้ เมื่อมีโรงพยาบาล 100 แห่ง ตัวคูณนี้สูงกว่า 1 เพียงเล็กน้อย

Stata: การประมาณเดียวกันด้วย vce(robust)

    * Stata's robust standard error: the same sandwich times n/(n - 1)
    poisson admissions if dataset == `d', vce(robust) nolog

Poisson regression                                      Number of obs =      4
                                                        Wald chi2(0)  =      .
                                                        Prob > chi2   =      .
Log pseudolikelihood = -13.414897                       Pseudo R2     = 0.0000

------------------------------------------------------------------------------
             |               Robust
  admissions | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   1.609438   .4618802     3.48   0.000     .7041693    2.514706
------------------------------------------------------------------------------

ช่วงเชื่อมั่น: ศูนย์กลางเดียวกัน ความกว้างต่างกัน

กลับไปที่โรงพยาบาล 100 แห่งในตารางตอนต้น จำนวนนับของพวกเขาแปรปรวนตั้งแต่ 0 ถึง 10 โดยมีความแปรปรวนเชิงประจักษ์ 18 เทียบกับค่าเฉลี่ย 5 ส่วน meat $K$ คือ 1,800 เทียบกับ bread $J$ ที่ 500 แซนด์วิชจึงทำให้ standard error กว้างขึ้น standard error แบบ HC0 คือ $\sqrt{1800}/500 = 0.0849$

vce(robust) ของ Stata เพิ่มตัวคูณ $n/(n - 1)$ และพิมพ์ 0.0853 ซึ่งเป็น 1.91 เท่าของค่าแบบอิงแบบจำลอง 0.0447 ช่วงเชื่อมั่น 95% ทั้งสองช่วงมีศูนย์กลางที่สัมประสิทธิ์เดียวกัน 1.6094 และช่วงแบบ robust กว้างเป็น 1.91 เท่า บนมาตรวัดลอการิทึม ช่วงนี้ไปตั้งแต่ 1.4423 ถึง 1.7766 เทียบกับ 1.5218 ถึง 1.6971 การยกกำลังด้วยเลขชี้กำลังกับขีดจำกัดทั้งสองให้ขีดจำกัดของจำนวนนับเฉลี่ย ซึ่งยังคงเป็น 5 ครั้งต่อวัน

กลุ่มโรงพยาบาลที่คงที่อีก 100 แห่งมีค่าเฉลี่ยเท่ากัน แต่จำนวนนับมีเพียง 4, 5 หรือ 6 ความแปรปรวนเชิงประจักษ์ของกลุ่มนี้คือ 0.5 ช่วงเชื่อมั่นแบบ robust ของมันคือ 1.5816 ถึง 1.6373 จึงแคบกว่าช่วงแบบอิงแบบจำลองของมันเอง

แบบจำลองเดียวสามารถรวมทั้งสองกลุ่มได้ โดยมีตัวบ่งชี้ (indicator) สำหรับกลุ่มที่แปรปรวน สัมประสิทธิ์ของตัวบ่งชี้คือลอการิทึมของอัตราส่วนอัตราอุบัติการณ์ (rate ratio) คืออัตราส่วนของจำนวนนับเฉลี่ยต่อวันของสองกลุ่ม เมื่อค่าเฉลี่ยเท่ากัน rate ratio เป็น 1.00 และมีเพียงความกว้างของช่วงที่ขึ้นกับ standard error ช่วงแบบ robust ของมันกว้างกว่า เพราะความกระจายส่วนเกินของกลุ่มที่แปรปรวนมากกว่าส่วนที่ขาดของกลุ่มที่คงที่

ช่วงเชื่อมั่นแบบอิงแบบจำลองและแบบ robust ในโรงพยาบาลสองกลุ่ม

Stata: กลุ่มที่แปรปรวน แบบอิงแบบจำลองและแบบ robust

    poisson admissions if volatile == `g', nolog
    scalar sc_b = _b[_cons]
    scalar sc_sm = _se[_cons]
    canon `k'.b0 4 _b[_cons]
    canon `k'.se_model 4 _se[_cons]
    predict double mu if volatile == `g', n
    generate double r2 = (admissions - mu)^2 if volatile == `g'
    matrix accum J = one if volatile == `g' [iweight = mu], noconstant
    matrix accum K = one if volatile == `g' [iweight = r2], noconstant
    matrix Ji = invsym(J)
    matrix V = Ji*K*Ji
    canon `k'.bread 4 el(J, 1, 1)
    canon `k'.meat 4 el(K, 1, 1)
    canon `k'.v_model 6 el(Ji, 1, 1)
    canon `k'.v_hc0 6 el(V, 1, 1)
    canon `k'.se_hc0 4 sqrt(el(V, 1, 1))
    canon `k'.ratio_hc0_model 4 sqrt(el(V, 1, 1))/scalar(sc_sm)
    poisson admissions if volatile == `g', vce(robust) nolog

Poisson regression                                      Number of obs =    100
                                                        LR chi2(0)    =   0.00
                                                        Prob > chi2   =      .
Log likelihood = -365.2899                              Pseudo R2     = 0.0000

------------------------------------------------------------------------------
  admissions | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   1.609438   .0447214    35.99   0.000     1.521786     1.69709
------------------------------------------------------------------------------
CANON g1.volatile.b0 1.6094
CANON g1.volatile.se_model 0.0447
(100 missing values generated)
(100 missing values generated)
(obs=500)
(obs=1800)
CANON g1.volatile.bread 500.0000
CANON g1.volatile.meat 1800.0000
CANON g1.volatile.v_model 0.002000
CANON g1.volatile.v_hc0 0.007200
CANON g1.volatile.se_hc0 0.0849
CANON g1.volatile.ratio_hc0_model 1.8974

Poisson regression                                      Number of obs =    100
                                                        Wald chi2(0)  =      .
                                                        Prob > chi2   =      .
Log pseudolikelihood = -365.2899                        Pseudo R2     = 0.0000

------------------------------------------------------------------------------
             |               Robust
  admissions | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   1.609438   .0852803    18.87   0.000     1.442292    1.776584
------------------------------------------------------------------------------

ผลลัพธ์แบบใช่หรือไม่ใช่: Poisson ดัดแปลง

การคำนวณแบบเดียวกันอธิบายทางเลือกที่พบบ่อยสำหรับผลลัพธ์แบบสองค่า ผลลัพธ์แบบใช่หรือไม่ใช่ที่มีความเสี่ยง $p$ มีความแปรปรวนแบบ Bernoulli (Bernoulli variance) เท่ากับ $p(1 - p)$ ในขณะที่แบบจำลอง Poisson สมมติว่าเท่ากับ $p$ ซึ่งคือค่าเฉลี่ย ที่ $p = 0.40$ ความแปรปรวนแบบ Bernoulli คือ $0.40 \times 0.60 = 0.24$ เทียบกับ Poisson ที่ 0.40 สำหรับความเสี่ยงใดก็ตามระหว่าง 0 ถึง 1 ค่า $p(1 - p)$ ต่ำกว่า $p$

แบบจำลอง Poisson จึงประเมินความแปรปรวนของผลลัพธ์แบบสองค่าสูงเกินไป standard error แบบอิงแบบจำลองของมันมักใหญ่เกินไป และแซนด์วิชอ่านความกระจายที่เล็กกว่านั้นจากข้อมูลPoisson ดัดแปลง (modified Poisson regression) ใช้คู่นี้พอดี [4] วิธีนี้ประมาณแบบจำลอง Poisson ที่มี log link สำหรับ risk ratio และเอา standard error จากแซนด์วิช

ชุดข้อมูลจำลองมีทารกแรกเกิด 200 คน ครึ่งหนึ่งได้รับการฝากครรภ์ (antenatal care, ANC) ไม่เพียงพอ ภาวะเจ็บป่วยของทารกแรกเกิดภายในสี่สัปดาห์หลังคลอดมีความเสี่ยง 0.40 หลังได้รับ ANC ไม่เพียงพอ และ 0.20 หลังได้รับ ANC เพียงพอ การประมาณแบบ Poisson ทั้งสองแบบให้ risk ratio 2.00 คือสัมประสิทธิ์ 0.6931 บนมาตรวัดลอการิทึม standard error แบบอิงแบบจำลองคือ 0.2739 ส่วนแซนด์วิชให้ 0.2345 แบบ HC0 และ 0.2351 จาก vce(robust)

ในกรณีนี้ robust standard error เป็นค่าที่เล็กกว่า ตรงตามที่การคำนวณแบบ Bernoulli ทำนายไว้ แบบจำลอง log-binomial (log-binomial model) สร้างแบบจำลองให้ลอการิทึมของความเสี่ยงและมีความแปรปรวนแบบ Bernoulli อยู่ในตัว มันให้ risk ratio เดียวกันคือ 2.00 และ standard error 0.2345 ซึ่งเป็นค่า HC0 บนมาตรวัดลอการิทึม ช่วงเชื่อมั่น 95% ของสัมประสิทธิ์ Poisson แคบลงจาก 0.1564 ถึง 1.2299 (แบบอิงแบบจำลอง) เป็น 0.2323 ถึง 1.1540 (แบบ robust)

การถดถอยโลจิสติกให้ odds ratio 2.67 แทน เพราะสร้างแบบจำลองให้ log odds ทั้งสองอัตราส่วนเป็นการเปรียบเทียบที่ไม่ได้ปรับของทารกแรกเกิด 200 คนกลุ่มเดียวกัน odds ratio อยู่ห่างจาก 1 มากกว่าเพราะผลลัพธ์พบบ่อย ไม่ใช่เพราะตัวกวนบทความเรื่อง risk ratio เมื่อผลลัพธ์พบบ่อยอธิบายว่าเมื่อใดควรเลือก Poisson ดัดแปลง และควรลองอะไรก่อน

ผลลัพธ์แบบสองค่าหนึ่งตัว แบบจำลองสามแบบ

Stata: Poisson ดัดแปลงด้วย vce(robust)

* modified Poisson: the same coefficient with Stata's robust (sandwich) standard error
poisson morbidity inadequate_anc, irr vce(robust) nolog

. * modified Poisson: the same coefficient with Stata's robust (sandwich) standard error
. poisson morbidity inadequate_anc, irr vce(robust) nolog

Poisson regression                                      Number of obs =    200
                                                        Wald chi2(1)  =   8.69
                                                        Prob > chi2   = 0.0032
Log pseudolikelihood = -128.84039                       Pseudo R2     = 0.0257

--------------------------------------------------------------------------------
               |               Robust
     morbidity |        IRR   std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc |          2   .4702186     2.95   0.003      1.26155    3.170702
         _cons |         .2   .0401004    -8.03   0.000     .1350089    .2962767
--------------------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.

Stata: การถดถอยโลจิสติกเพื่อเปรียบเทียบ

* logistic regression for contrast: an odds ratio, not a risk ratio
logistic morbidity inadequate_anc, nolog

. * logistic regression for contrast: an odds ratio, not a risk ratio
. logistic morbidity inadequate_anc, nolog

Logistic regression                                     Number of obs =    200
                                                        LR chi2(1)    =   9.66
                                                        Prob > chi2   = 0.0019
Log likelihood = -117.34141                             Pseudo R2     = 0.0395

--------------------------------------------------------------------------------
     morbidity | Odds ratio   Std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc |   2.666666   .8606626     3.04   0.002     1.416591    5.019872
         _cons |   .2500001      .0625    -5.55   0.000      .153158    .4080755
--------------------------------------------------------------------------------
Note: _cons estimates baseline odds.

สิ่งที่แซนด์วิชแก้ได้ และสิ่งที่จะตามมา

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

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

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

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

• เปรียบเทียบความกระจายที่แบบจำลองสมมติกับความกระจายในข้อมูล เมื่อไม่มีตัวแปรร่วม หมายถึงความแปรปรวนของผลลัพธ์เทียบกับค่าเฉลี่ยของมัน เมื่อมีตัวแปรร่วม ให้ใช้ส่วนตกค้างหลังประมาณแบบจำลอง เพราะตัวแปรร่วมเพียงอย่างเดียวก็ทำให้ความแปรปรวนดิบเกินค่าเฉลี่ยได้ ช่องว่างที่ใหญ่อาจเป็นสัญญาณว่าข้อสมมติเรื่องความแปรปรวนผิด
• รายงานว่าใช้ standard error แบบใด และตัวคูณปรับสำหรับตัวอย่างเล็กใดที่โปรแกรมของคุณใช้ เพื่อให้ผู้อื่นทำช่วงเชื่อมั่นซ้ำได้
• อ่านช่องว่างที่ใหญ่ระหว่าง standard error แบบอิงแบบจำลองกับแบบ robust เป็นสัญญาณให้ทบทวนแบบจำลองอีกครั้ง
• เมื่อมีหน่วยสังเกตน้อย ให้ระวังเมื่อ robust standard error เล็กกว่าแบบอิงแบบจำลอง
• สำหรับผลลัพธ์แบบสองค่าที่ต้องการ risk ratio ให้พิจารณาแบบจำลอง log-binomial ก่อน หากไม่ลู่เข้า ให้ถอยลงมาใช้ Poisson ดัดแปลงที่มี robust standard errorบทความเรื่อง risk ratio เมื่อผลลัพธ์พบบ่อยอธิบายลำดับนั้น
• ตรวจสคริปต์ยาว ๆ เทียบกับการคำนวณด้วยมือบนชุดข้อมูลเล็กจิ๋ว เช่นผู้ป่วยสี่คนในที่นี้

อ่านต่อในวิกิ: [[heteroskedasticity-robust-standard-errors-th]] [[risk-regression-models-epidemiology-th]] [[choosing-vce-robust-cluster-standard-errors-th]]