ทำความเข้าใจ VCE: Robust Standard Error, Sandwich Estimator และ vce(cluster id)

Read the English version

บทคัดย่อ

Regression coefficient กับความไม่แน่นอนของมันมาจากการคำนวณสองส่วนที่แยกกัน และออปชัน vce() ของ Stata เปลี่ยนเฉพาะส่วนที่สอง Standard Error ไม่ใช่ความกระจายของข้อมูล: เมื่อ Variance เท่ากับ 18 จากโรงพยาบาล 100 แห่ง Standard Error ของค่าเฉลี่ยเท่ากับ 0.4243 Stata ประมาณ Variance-Covariance Matrix หรือ VCE โดย Standard Error คือรากที่สองของค่าบนแนวทแยง และต้องใช้ Covariance ในนั้นทุกครั้งที่รวม coefficients เข้าด้วยกัน Model-based VCE เชื่อ Variance ของ Model ส่วน Robust Sandwich วาง Meat ที่สร้างจาก Squared residuals ไว้ระหว่าง Bread ที่กลับค่าแล้วสองชุด เมื่อคำนวณด้วยมือสำหรับโรงพยาบาล 100 แห่ง Intercept-only Poisson Model ให้ Model-based SE เท่ากับ 0.0447 และ Robust SE เท่ากับ 0.0849 โดย coefficient เท่าเดิม ใน Modified Poisson Regression สำหรับ Binary outcome Robust SE อาจเป็นค่าที่เล็กกว่า บทความนี้ยังเปรียบเทียบ vce(robust) กับ vce(cluster id) แสดงสิ่งที่ Robust SE แก้ไม่ได้ และให้ประโยคสำหรับส่วน Methods และ Results

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

Coefficient เท่าเดิม แต่ความไม่แน่นอนต่างกัน

ทีมวิจัยทีมหนึ่งนับจำนวนผู้ป่วยที่ Admit ต่อวันในโรงพยาบาล 100 แห่ง แล้ว Fit Poisson Model ใน Stata พร้อม vce(robust) เพราะงานครั้งก่อนก็ใช้แบบนี้ อีก Dataset หนึ่งติดตามผู้ป่วยกลุ่มเดิมผ่านการ Admit หลายครั้ง และโค้ดของมันเขียนว่า vce(cluster patient_id) หลายคนเข้าใจทั้งสองคำสั่งเพียงว่าเป็นคำสั่งสำหรับปรับ Standard Error โดยไม่ทราบว่าทำไมต้องปรับ หรือเกิดอะไรขึ้นภายใน Model

บทความนี้จะอธิบายตั้งแต่ Variance, Standard Error, Variance-Covariance Estimator (VCE), Model-based SE, Robust SE และ Sandwich Estimator ไปจนถึงการประยุกต์กับ Poisson Regression, Modified Poisson Regression และ Clustered Data

ทำไมเราต้องมี Standard Error

สมมติว่าเราต้องการศึกษาจำนวนผู้ป่วยที่ Admit ต่อวันในโรงพยาบาล 100 แห่ง เมื่อเก็บข้อมูลครบแล้ว พบว่าจำนวน Admit เฉลี่ยเท่ากับ 5 คนต่อวัน

แต่ค่าเฉลี่ยเพียงค่าเดียวไม่ได้บอกว่าข้อมูลของแต่ละโรงพยาบาลแตกต่างกันมากแค่ไหน ลองพิจารณาข้อมูลสองชุดที่มี Mean เท่ากัน Dataset A กระจายตัวน้อย: Mean 5, Variance 1 ส่วน Dataset B กระจายตัวมาก: Mean 5, Variance 16 รูปที่ 1 ซึ่งแสดงเป็นตัวเลขในตารางด้านล่าง แสดงโรงพยาบาลชุดละแปดแห่ง เพื่อให้เห็นรูปแบบการกระจาย

ถ้าเราเก็บข้อมูลใหม่จากโรงพยาบาลอีก 100 แห่ง ค่าเฉลี่ยที่ประมาณได้จากประชากรแบบ B มีแนวโน้มเปลี่ยนแปลงระหว่างตัวอย่างมากกว่าประชากรแบบ A ภายใต้สมมติฐานว่าหน่วยข้อมูลเป็นอิสระและมีการแจกแจงเหมือนกันภายในแต่ละประชากร นี่คือเหตุผลที่ต้องมี Standard Error

รูปที่ 1: Mean เท่ากัน แต่ความกระจายต่างกัน

Variance ต่างจาก Standard Error อย่างไร

Variance ใช้อธิบายความกระจายของข้อมูล ส่วน Standard Error (SE) ใช้อธิบายความไม่แน่นอนของค่าประมาณจากการสุ่มตัวอย่างซ้ำ ตัวอย่างเช่น หากมี $n$ observations ที่เป็นอิสระกัน และแต่ละตัวมี Variance เท่ากับ $\sigma^2$ Variance ของค่าเฉลี่ยคือ

ดังนั้น

สมมติว่าจำนวน Admit ของโรงพยาบาล 100 แห่งมี Mean เท่ากับ 5 และ Variance เท่ากับ 18

ความหมายคือ ค่าเฉลี่ยที่ประมาณได้มี Sampling uncertainty ประมาณ 0.4243 ซึ่งวัดบน scale เดียวกับจำนวน Admit ค่า 0.4243 ไม่ใช่ Standard Deviation ของโรงพยาบาลแต่ละแห่ง แต่เป็น Standard Error ของ estimated mean

ทีนี้ให้ Variance เท่าเดิมที่ 18 แต่เพิ่มจำนวนโรงพยาบาลจาก 100 เป็น 400

SE ลดลงประมาณครึ่งหนึ่ง เพราะจำนวน observations เพิ่มขึ้นสี่เท่า จึงต้องแยกให้ได้ว่า จำนวนข้อมูล (sample size) กับความกระจายของข้อมูล (variance) ต่างก็มีผลต่อ Standard Error

VCE คืออะไร

VCE ย่อมาจาก Variance-Covariance Estimator คือวิธีการที่ใช้ประมาณ Variance และ Covariance ของค่าประมาณทางสถิติ ใน Stata เมื่อเราเขียน

poisson admissions treatment age, vce(robust)

คำว่า vce(robust) หมายถึง ให้ใช้ Robust Variance-Covariance Estimator ในการประมาณความไม่แน่นอนของ Regression coefficients [1]

สิ่งที่ Stata คำนวณไม่ใช่ SE เพียงค่าเดียว แต่เป็น Variance-Covariance Matrix ของ Regression coefficients ทั้งหมด

ทำไมต้องเป็น Matrix

สมมติเรามี Regression Model

โดย $\beta_0$ คือ Intercept และ $\beta_1$ คือ Regression coefficient ของ $X$ จากการวิเคราะห์ เราได้

ค่าประมาณทั้งสองมีความไม่แน่นอนของตัวเอง และยังอาจเปลี่ยนแปลงร่วมกันเมื่อสุ่มข้อมูลใหม่ด้วย VCE จึงต้องประมาณ Matrix ที่มีหน้าตาแบบนี้

VCE อยู่ตรงไหนใน Stata

หลังจาก Fit Regression Model แล้ว Stata จะเก็บผลสำคัญไว้ใน object ตามตารางด้านล่าง สามารถเรียกดูได้ด้วยคำสั่ง

matrix list e(b)
matrix list e(V)
estat vce

ถ้า Model มี 4 coefficients รวม Intercept แล้ว e(b) จะมีขนาด 1 × 4 ขณะที่ e(V) มีขนาด 4 × 4 รายละเอียดนี้เป็นโครงสร้างมาตรฐานของ Stata estimation results [1]

VCE เกี่ยวข้องกับ Confidence Interval และ P-value อย่างไร

ให้สมมติว่าการทำ Regression ได้ $\hat\beta_1 = 0.40$ และ Standard Error เท่ากับ 0.10 เราสามารถคำนวณ Wald 95% Confidence Interval บน coefficient scale ด้วย Normal approximation

ถ้าปรับมาใช้ Robust VCE แล้ว SE เปลี่ยนเป็น 0.19 โดย coefficient ยังคงเป็น 0.40

รูปที่ 2: coefficient เดียวกัน สองช่วง

Model-based VCE คืออะไร

ใน Regression Model ทั่วไป Software ต้องมีวิธีประมาณ Variance ของ coefficients วิธีมาตรฐานที่พบบ่อยใน Maximum Likelihood Estimation คือ Model-based VCE ซึ่งใช้ข้อมูลจาก Likelihood และสมมติฐานของ Model รากที่สองของค่าบนแนวทแยงของมันคือ Model-based SE ในหลาย Model เราสามารถประมาณจาก Inverse ของ Observed Information Matrix หรือ Negative Hessian

โดย

นี่เป็นสูตรของ observed-information VCE ไม่ใช่สูตรที่ใช้แทน VCE ทุกชนิดได้โดยไม่มีเงื่อนไข

เพื่อเข้าใจง่ายขึ้น ให้มองว่า Hessian เป็นข้อมูลที่บอกว่า Likelihood เปลี่ยนเร็วเพียงใดเมื่อเราเปลี่ยนค่า Regression coefficient ถ้า Likelihood มีความโค้งมากในบริเวณค่าที่ประมาณได้ ก็มักหมายความว่า coefficient ถูกระบุได้ค่อนข้างแม่นยำ แต่ถ้า Likelihood ค่อนข้างแบน ค่าประมาณ coefficient จะมีความไม่แน่นอนมากขึ้น

ตัวอย่างใน Poisson Regression

สมมติเราสร้าง Poisson Regression Model

โดยกำหนด

หรือเทียบเท่ากับ $\mu_i = \exp(X_i^{\top}\boldsymbol\beta)$ Poisson distribution มีคุณสมบัติสำคัญคือ

แปลว่า Conditional variance เท่ากับ Conditional mean ถ้า Model คาดว่าโรงพยาบาลหนึ่งแห่งมีจำนวน Admit เฉลี่ย 5 คน $E(Y_i \mid X_i) = 5$ Poisson Model ก็จะใช้ $Var(Y_i \mid X_i) = 5$ สำหรับการคำนวณ Model-based uncertainty ภายใต้สมมติฐาน Poisson

แล้วถ้าข้อมูลจริงมี Variance = 18?

กรณีนี้ Conditional variance สูงกว่าค่าเฉลี่ยตาม Poisson assumption

เรียกว่า Overdispersion เมื่อ Poisson Model ใช้ Variance = 5 แต่ความผันผวนจริงสูงกว่า การคำนวณ Model-based SE อาจประเมินความไม่แน่นอนต่ำเกินไป

ต้องระวังว่า การเปรียบเทียบ Mean และ Variance ของข้อมูลรวมทั้ง Dataset ยังไม่ใช่การพิสูจน์ Conditional Overdispersion เมื่อ Model มี predictors หลายตัว เพราะความแตกต่างของค่าเฉลี่ยระหว่างกลุ่มก็ทำให้ Overall variance สูงขึ้นได้ สิ่งที่เราต้องประเมินคือ Variance ภายใต้เงื่อนไขของ predictors ที่ระบุไว้ใน Model

Robust VCE คืออะไร

Robust VCE เป็นวิธีประมาณ Variance-Covariance Matrix ที่ลดการพึ่งพาสมมติฐานบางส่วนของ Model โดยเฉพาะสมมติฐานเกี่ยวกับ Variance วิธีที่ใช้กันแพร่หลายคือ Huber-White Sandwich Variance Estimator และรากที่สองของค่าบนแนวทแยงของมันคือ Robust (Sandwich) SE [2]

ใน Poisson Regression หาก Conditional mean model ถูกต้อง แต่ Conditional variance ไม่เป็นไปตาม Poisson distribution เรายังสามารถประมาณ Regression coefficients ด้วย Poisson estimating equations และใช้ Robust VCE สำหรับ Statistical inference ได้ แนวคิดนี้เป็นพื้นฐานของ Poisson Pseudo-Maximum Likelihood และ Modified Poisson Regression

เปรียบเทียบ Model-based กับ Robust VCE

Sandwich Estimator คำนวณอย่างไร

เมื่อเข้าใจแล้วว่าเราต้องประมาณ Variance ของ coefficients ขั้นต่อไปคือทำความเข้าใจวิธีคำนวณ Sandwich Estimator มีสูตรพื้นฐานว่า

สูตรนี้เป็นรูปแบบพื้นฐานสำหรับ independent observations และยังไม่รวม finite-sample correction มันถูกเรียกว่า Sandwich เพราะมีองค์ประกอบสามส่วน

• Bread $A^{-1}$ ด้านซ้าย
• Meat $B$ ตรงกลาง
• Bread $A^{-1}$ ด้านขวา

Bread คืออะไร

สำหรับ Poisson Regression ที่มี Log link

โดย $X_i$ คือ Vector ของ predictors รวม Intercept Bread ทำหน้าที่เชื่อมความผันผวนของ Estimating equations เข้ากับความผันผวนของ Estimated coefficients มันสะท้อนความไวของสมการที่ใช้ประมาณ coefficient ต่อการเปลี่ยนแปลงของค่า coefficient นั้น

Meat คืออะไร

ก่อนอื่นต้องรู้จัก Score contribution ของ Observation แต่ละตัว สำหรับ Poisson Regression

โดย $Y_i - \hat\mu_i$ คือ Raw residual หรือความแตกต่างระหว่างค่าที่สังเกตได้กับค่าที่ Model ทำนาย จากนั้นนำ Score contribution ของแต่ละ Observation มาสร้าง Meat

แทนสูตร Score ของ Poisson ลงไป

สังเกตว่าสูตรนี้ใช้ Squared residual ของแต่ละ observation ดังนั้นหากข้อมูลมีความผันผวนมาก Score contributions ก็สามารถสะท้อนความผันผวนนั้นใน Meat ได้

ทำไม Model-based กับ Robust VCE จึงต่างกัน

ถ้า Poisson variance assumption ถูกต้อง

ดังนั้นในเชิงค่าคาดหมาย Meat จะสอดคล้องกับ Bread

Sandwich VCE จึงใกล้เคียงกับ

ซึ่งก็คือ Model-based VCE แต่ถ้า Conditional variance ไม่เป็นไปตาม Poisson assumption

Meat กับ Bread จะไม่จำเป็นต้องเท่ากันอีกต่อไป จึงทำให้

นี่คือเหตุผลทางคณิตศาสตร์ว่าทำไม Robust SE จึงแตกต่างจาก Model-based SE Stata ใช้หลักการ Huber-White Sandwich ใน vce(robust) และขยายแนวคิดนี้ไปสู่ Cluster-robust VCE ได้เช่นกัน

Worked Example: คำนวณ Robust SE ด้วยมือทีละขั้น

ทดลองคำนวณด้วย Stata

เราสามารถสร้าง Dataset เดียวกันและทดลองใน Stata ได้โดยไม่ต้องใช้ข้อมูลผู้ป่วยจริง คำสั่งของตัวอย่างคำนวณด้วยมือมีดังนี้

* Create synthetic admission count data
clear

input admissions frequency
0  20
1  25
5  10
9  25
10 20
end

* Expand to 100 independent observations
expand frequency
drop frequency

* Examine the data
count
summarize admissions, detail

* Model 1: Conventional Poisson
poisson admissions
estimates store poisson_model

display "Coefficient = " _b[_cons]
display "Model-based SE = " _se[_cons]

* Model 2: Robust Poisson
poisson admissions, vce(robust)
estimates store poisson_robust

display "Coefficient = " _b[_cons]
display "Robust SE = " _se[_cons]

* Inspect coefficients and VCE
matrix list e(b)
matrix list e(V)
estat vce

* Compare the two models
estimates table poisson_model poisson_robust, b se

เมื่อรันแล้วให้สังเกตผลสำคัญตามตารางด้านล่าง

Stata: โรงพยาบาล 100 แห่ง แบบ Model-based และ 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
------------------------------------------------------------------------------

Modified Poisson Regression: ทำไมจึงต้องใช้ Robust SE

นี่คืออีกหนึ่งการประยุกต์ที่สำคัญมากใน Clinical Epidemiology คือการประมาณ Risk Ratio (RR) จาก Binary outcome ด้วย Modified Poisson Regression ซึ่งเป็น Log-link Poisson Model ที่ Fit กับ outcome ที่มีค่า 0 หรือ 1 สมมติเราศึกษาความสัมพันธ์ระหว่าง Inadequate Antenatal Care (ANC) กับ Neonatal Morbidity ภายใน 28 วัน Outcome เป็น Binary variable

Binary outcome มี Variance เท่าไร

สำหรับ Binary outcome ที่มีโอกาสเกิดเหตุการณ์เท่ากับ $p$

แต่ Variance ของมัน คือ Bernoulli variance

ต่างจาก Poisson ที่กำหนด $Var(Y) = p$ ตัวอย่างเช่น หาก Risk เท่ากับ 0.40

ตัวอย่างคำนวณ RR จากข้อมูล Binary

คำสั่ง Stata สำหรับ Modified Poisson

* Binary outcome: 0 = No morbidity, 1 = Morbidity
* Independent observations

poisson morbidity28 i.inadequate_anc ///
    i.age_group i.low_ses i.coverage_gap, ///
    irr vce(robust)

แม้ Stata จะแสดงหัวตารางว่า IRR ซึ่งย่อมาจาก Incidence Rate Ratio เพราะเป็นคำสั่ง poisson, irr แต่ถ้า Outcome เป็น Binary และ Model กำลังประมาณ Risk ภายใน Follow-up horizon เดียวกัน โดยไม่มี Person-time offset ค่า $\exp(\hat\beta)$ จะถูกตีความเป็น RR ภายใต้ Log-risk mean model ชื่อ IRR ที่ Stata แสดงไม่ได้ทำให้ค่าประมาณกลายเป็น Incidence Rate Ratio โดยอัตโนมัติ

หากใส่ exposure(tmorbidity28) เพื่อใช้ Person-time เป็น denominator โครงสร้าง Mean model จะเปลี่ยนไปเป็น Rate model ซึ่งต้องตีความ effect measure ใหม่

อีกเรื่องที่ต้องระวังคือ Log-link Risk Model อาจทำนายค่าความน่าจะเป็นสูงกว่า 1 ในบาง Covariate profiles ได้ Robust SE ไม่ได้แก้ปัญหานี้

vce(robust) หรือ vce(cluster patient_id)?

จนถึงตอนนี้เรากำลังสมมติว่า Observations ทุกตัวเป็นอิสระจากกัน แต่ข้อมูลทางคลินิกหลายประเภทไม่ได้เป็นอิสระจริง ตัวอย่างเช่น เราศึกษาจำนวน Hospital Admissions ของผู้ป่วย 100 คน โดยผู้ป่วยแต่ละคนมีข้อมูลติดตามหลายครั้ง Observations จากผู้ป่วยคนเดียวกันอาจมีความสัมพันธ์กัน

ถ้าใช้ vce(robust) ธรรมดาโดยปฏิบัติต่อทุกแถวเสมือนเป็น Independent observation เราอาจประเมินความไม่แน่นอนผิด จึงต้องพิจารณา Cluster-robust VCE และ Cluster-robust SE ที่ได้จากมัน

Cluster-robust VCE ทำงานอย่างไร

สำหรับ Robust VCE ปกติ Meat ใช้ Score contribution ของ observation แต่ละตัว

แต่ถ้ามีผู้ป่วยหลายคน และแต่ละคนมีข้อมูลหลายแถว เราจะรวม Score contributions ภายในผู้ป่วยคนเดียวกันก่อน กำหนดให้ $g$ เป็น Cluster ในที่นี้คือ Patient

จากนั้นสร้าง Cluster-level Meat

จึงได้

โดยยังไม่ได้แสดง finite-sample corrections แนวคิดนี้ช่วยรองรับความสัมพันธ์ของ Observations ภายใน Cluster โดยมีเงื่อนไขสำคัญสองอย่าง คือ Clusters ต้องมีความเป็นอิสระระหว่างกัน และโดยทั่วไปต้องมีจำนวน Independent clusters มากเพียงพอ [5] วิธีนี้ได้รับการประยุกต์กับ Modified Poisson Regression สำหรับ correlated binary outcomes ด้วย [4]

คำสั่ง Stata

* Independent observations
poisson admissions i.treatment c.age, vce(robust)

* Repeated observations within each patient
poisson admissions i.treatment c.age, ///
    vce(cluster patient_id)

vce(cluster patient_id) ไม่ได้หมายความว่า Software จะสร้าง Matched pairs, แก้ Confounding หรือเพิ่ม Patient fixed effects มันเปลี่ยนวิธีประมาณ Variance โดยใช้ Cluster เป็นหน่วยอิสระในการคำนวณ Score variability

ถ้าจำนวน Clusters น้อย Cluster-robust SE แบบมาตรฐานอาจให้ Confidence Interval และ P-value ที่ไม่น่าเชื่อถือ ต้องพิจารณา Small-sample adjustments หรือวิธีการอนุมานที่เหมาะสมกับจำนวน Clusters

Stata มี VCE กี่ประเภท

Stata รองรับ VCE หลายรูปแบบ โดยแบบที่ใช้ได้จริงขึ้นอยู่กับ Estimation command [1]

OPG, bootstrap และ jackknife

vce(opg) ไม่ได้เท่ากับ Robust Sandwich VCE แม้จะใช้ Score information เช่นกัน เพราะ Robust Sandwich ใช้ทั้ง Bread และ Meat

ส่วน vce(bootstrap) และ vce(jackknife) ไม่ได้เหนือกว่า Robust SE ในทุกกรณี หากใช้กับ Clustered data ก็ต้องกำหนด Resampling unit ให้สอดคล้องกับการออกแบบการศึกษา

ตารางตัดสินใจสำหรับ Clinical Research

Robust SE ทำอะไรไม่ได้

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

Robust SE ไม่ได้ทำให้ Model ถูกต้องทุกอย่าง

สมมติ Conditional mean ที่แท้จริงคือ

แต่เรากลับ Fit Model ที่ไม่ใส่ $X^2$

แม้จะใช้

poisson outcome x, vce(robust)

เราก็ยังไม่สามารถสรุปได้ว่า coefficient ของ $X$ เป็นค่าที่ถูกต้องสำหรับความสัมพันธ์ตาม Model ที่แท้จริง Robust VCE อาจให้ความไม่แน่นอนของ Pseudo-true coefficient ซึ่งคือค่าที่ coefficient ของ Model ที่ระบุผิดลู่เข้าหา แต่ไม่ได้ทำให้ Coefficient นั้นกลายเป็น Parameter ทางคลินิกหรือ Causal effect ที่เราต้องการโดยอัตโนมัติ

เช่นเดียวกัน หากเรานำ Case-control sample ที่สุ่ม Cases และ Controls ในสัดส่วนที่ไม่สะท้อนประชากรต้นทางไป Fit Modified Poisson แล้วรายงาน RR โดยไม่คำนึงถึง Sampling design การเพิ่ม vce(robust) ก็ไม่ได้แก้ปัญหานี้ บทความเรื่องการสุ่มกลุ่มควบคุมแสดงว่าการออกแบบเป็นตัวกำหนด effect measure อย่างไร

Robust Poisson หรือ Negative Binomial?

หากศึกษา Count outcome และพบ Overdispersion เราอาจใช้ได้ทั้ง Robust Poisson และ Negative Binomial Regression แต่ทั้งสองวิธีไม่ได้มีเป้าหมายเหมือนกันทุกกรณี

Robust Poisson กำหนด Conditional mean

และใช้ Sandwich VCE โดยไม่จำเป็นต้องกำหนด Poisson variance ให้ถูกต้อง หาก Conditional mean ถูกต้องและ Independence assumptions เหมาะสม Poisson Pseudo-Maximum Likelihood สามารถให้ Consistent coefficient estimates ได้ แม้ Conditional variance จะต่างจาก Conditional mean

Negative Binomial Regression เป็นอีก Model ที่กำหนด Mean-Variance relationship อย่างชัดเจน เช่นรูปแบบ NB2

โดย $\alpha$ เป็น Overdispersion parameter [6] ถ้า Distribution และ Variance structure ของ Negative Binomial เหมาะสม Model อาจมีประสิทธิภาพสูงกว่าในการประมาณ Parameters บางชนิด และมีประโยชน์ในการจำลองหรือพยากรณ์ Count distribution แต่ไม่ได้หมายความว่า Negative Binomial จะดีกว่า Robust Poisson เสมอไป

ควรเลือกจากวัตถุประสงค์ของการศึกษา ความเหมาะสมของ Mean model และ Distribution assumptions ไม่ใช่เลือกจากการพบ Overdispersion เพียงอย่างเดียว บทความเรื่องอัตราแบบ Poisson (บทความภาษาอังกฤษ) อธิบาย Count model เหล่านี้โดยตรง

ตัวอย่างการรายงานผลในงานวิจัย

Methods สำหรับ Modified Poisson Regression:

Adjusted risk ratios and 95% confidence intervals were estimated using modified Poisson regression with a log link and robust sandwich variance estimation. The analysis assumed independent observations and a correctly specified conditional mean model.

ถ้าข้อมูลมี Clustering ภายในโรงพยาบาล และใช้ Cluster-robust VCE จริง ควรเปลี่ยนเป็น

Adjusted risk ratios and 95% confidence intervals were estimated using modified Poisson regression with a log link and cluster-robust sandwich variance estimation, accounting for within-hospital correlation.

สองข้อความนี้มีความหมายต่างกัน และไม่ควรนำไปใช้แทนกันโดยไม่ตรวจสอบ Study design

Results: สมมติผลวิเคราะห์ได้ $aRR = 2.00$, 95% CI 1.26 ถึง 3.17 สามารถรายงานว่า

Inadequate antenatal care was associated with a higher risk of neonatal morbidity within 28 days (adjusted RR 2.00, 95% CI 1.26 to 3.17).

ตัวเลขนี้เป็นตัวอย่างสำหรับการเขียนเท่านั้น ไม่ใช่ผลการวิเคราะห์จาก Dataset จริง และหากนำไปใส่ Manuscript จะต้องใช้ตัวเลขจาก Adjusted model ที่รันจริง

สรุปสิ่งที่ต้องเข้าใจ

อ่านต่อในวิกิ: [[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]]