logit, ologit และ clogit: link function เดียว แต่ตอบคำถามต่างกันสามแบบ

Clinical Epidemiology ResearchMethodology and Research Design THUniqcret doctor knowledges TH
logit, ologit และ clogit: link function เดียว แต่ตอบคำถามต่างกันสามแบบ
On this page

Read the English version

บทคัดย่อ

logistic regression, ordinal logistic regression และ conditional logistic regression ใช้ link เดียวกันคือ logit สัมประสิทธิ์ที่ยกกำลังด้วยเลขชี้กำลังจึงเป็น odds ratio ชนิดของผลลัพธ์และแบบแผนการสุ่มตัวอย่างเป็นตัวกำหนดแบบจำลอง ในทารกแรกเกิดจำลอง 500 คน การฝากครรภ์ไม่เพียงพอมี odds ratio เท่ากับ 2.01 สำหรับระดับภาวะแทรกซ้อนที่สูงขึ้น ค่านี้ปรับด้วยการคลอดก่อนกำหนด ซึ่งอาจเป็นตัวแปรสื่อกลาง (mediator) ดังนั้น 2.01 จึงไม่รวมผลใด ๆ ที่ส่งผ่านการคลอดก่อนกำหนด แบบจำลองแบบลำดับสมมติว่ามี odds ratio ค่าเดียวที่ทุกจุดตัดระหว่างระดับ Brant test พบหลักฐานของการเบี่ยงเบนจากข้อสมมตินี้สำหรับการคลอดก่อนกำหนด (P = 0.003) แต่ไม่พบสำหรับการฝากครรภ์ ในคู่ case-control ที่จับคู่ 20 คู่ logistic regression ทั่วไปให้ odds ratio 1.83 ส่วน conditional logistic regression ให้ 4.00 ทั้งสองค่าตอบคำถามต่างกัน คือเปรียบเทียบแบบรวมกลุ่มกับเปรียบเทียบภายในคู่ และแบบแผนการศึกษาเป็นตัวกำหนดว่าควรรายงานค่าใด มีเพียง 5 คู่ที่ต่างกันในสิ่งสัมผัส ค่า 4.00 จึงไม่แม่นยำ (ช่วงเชื่อมั่น 95% 0.45 ถึง 35.79) บทความนี้คำนวณค่าประมาณทั้งสองด้วยมือ แล้วแสดงโค้ด Stata และ R และตารางสำหรับตัดสินใจ


ภาพสรุป ตัวอย่างคำนวณด้วยมือ

สองแบบแผน สามคำสั่ง และ odds ratio สองค่า

หน่วยทารกแรกเกิดแห่งหนึ่งต้องการทราบว่าการฝากครรภ์ไม่เพียงพอ (inadequate antenatal care, ANC) สัมพันธ์กับภาวะแทรกซ้อนในทารกแรกเกิดมากขึ้นหรือไม่ คำว่าไม่เพียงพอในที่นี้หมายถึงการตรวจระหว่างตั้งครรภ์ที่น้อยเกินไปหรือเริ่มช้าเกินไป ทีมศึกษาคำถามนี้ด้วยแบบแผนสองแบบพร้อมกัน หน่วยนี้ ทารก และตัวเลขทุกตัวด้านล่างเป็นข้อมูลจำลอง

ใน birth cohort ทารกแรกเกิด 500 คนแต่ละคนมีระดับภาวะแทรกซ้อนภายใน 28 วัน ได้แก่ ไม่มี เล็กน้อย ปานกลาง หรือรุนแรง ระดับเหล่านี้มีลำดับ แต่ช่วงห่างระหว่างระดับไม่จำเป็นต้องเท่ากัน อีกส่วนหนึ่งเป็นการศึกษาแบบ case-control ซึ่งจับคู่ทารกที่มีภาวะแทรกซ้อน 20 คน (case) กับกลุ่มควบคุม (control) คนละหนึ่งคน control แต่ละคนเป็นทารกที่ยังไม่มีภาวะแทรกซ้อนในวันที่ case ของมันมี

นักวิเคราะห์รันคำสั่ง Stata สามคำสั่ง คือ logistic, ologit และ clogit ในคู่ที่จับคู่ logistic ให้ odds ratio 1.83 และ clogit ให้ 4.00 ทารก สิ่งสัมผัส และผลลัพธ์เป็นชุดเดียวกัน ผู้ทบทวนจึงถามว่าตัวเลขไหนถูก

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

แบบจำลองทุกแบบในบทความนี้อธิบายความน่าจะเป็น $p$ เช่น ความน่าจะเป็นที่ทารกแรกเกิดมีภาวะแทรกซ้อน ความน่าจะเป็นถูกจำกัดอยู่ระหว่าง 0 ถึง 1 ขณะที่สมการถดถอยรับค่าใดก็ได้ ฟังก์ชันเชื่อมโยง (link function) เชื่อมสองสิ่งนี้เข้าด้วยกัน โดยแปลง $p$ ไปอยู่บนมาตรวัดที่ปรับแบบจำลองเชิงเส้นได้

logit link คือ log ของ odds:

$$\operatorname{logit}(p) = \ln\left(\frac{p}{1 - p}\right)$$

ในที่นี้ $p/(1 - p)$ คือ odds และ $\ln$ คือลอการิทึมธรรมชาติ logit มีค่าตั้งแต่ลบอนันต์ถึงบวกอนันต์ ผลรวมของสัมประสิทธิ์คูณตัวแปรร่วม ซึ่งคือ ตัวทำนายเชิงเส้น (linear predictor) จึงรับค่าใดก็ได้ และยังแปลงกลับเป็นความน่าจะเป็นที่ถูกต้องได้

สำหรับสิ่งสัมผัสแบบสองค่า $x$ ที่รหัส 1 คือการฝากครรภ์ไม่เพียงพอ และ 0 คือเพียงพอ logistic regression ทั่วไปประมาณแบบจำลอง

$$\operatorname{logit} P(Y = 1 \mid x) = \beta_0 + \beta_1 x$$

ในที่นี้ $Y = 1$ หมายถึงทารกที่มีผลลัพธ์ และ $\beta_0$ คือ log odds เมื่อ $x = 0$ ความชัน $\beta_1$ คือการเปลี่ยนแปลงของ log odds เมื่อ $x$ เปลี่ยนจาก 0 เป็น 1 ผลต่างของ log odds สองค่าคือ log ของอัตราส่วนของมัน ดังนั้น $\beta_1$ คือ log odds ratio และ $e^{\beta_1}$ คือ odds ratio [1]

คำสั่ง logit ของ Stata พิมพ์ $\beta_1$ ส่วนคำสั่ง logistic ประมาณแบบจำลองเดียวกันและพิมพ์ $e^{\beta_1}$ ฟังก์ชัน glm() ของ R ที่ใช้ family = binomial ใช้ logit link เป็นค่าเริ่มต้นและพิมพ์ $\beta_1$

logistic regression ทั่วไปสมมติสองอย่างเกี่ยวกับข้อมูล คือผลลัพธ์เป็นแบบสองค่า และหน่วยสังเกตเป็นอิสระต่อกัน หมายถึงผลลัพธ์ของคนหนึ่งไม่บอกอะไรเกี่ยวกับอีกคนเมื่อกำหนดตัวแปรร่วมแล้ว ระดับภาวะแทรกซ้อนที่มีลำดับฝืนข้อสมมติแรก และคู่ที่จับคู่ฝืนข้อสมมติหลัง การจับคู่ยังผูก control แต่ละคนไว้กับ case ของมัน จึงต้องเข้าไปอยู่ในการเปรียบเทียบเอง ไม่ใช่แค่ในค่าคลาดเคลื่อนมาตรฐาน

ระบุชนิดของผลลัพธ์และแบบแผนการสุ่มตัวอย่างก่อน

การวิเคราะห์สามแบบในบทความนี้ใช้ logit link ร่วมกัน แต่ต้องใช้แบบจำลองต่างกัน การเลือกขึ้นกับสองข้อเท็จจริงที่ทราบก่อนวิเคราะห์ข้อมูลใด ๆ สิ่งแรกคือชนิดของผลลัพธ์ (outcome type) มีสองหมวด หรือหลายหมวดที่มีลำดับ สิ่งที่สองคือแบบแผนการสุ่มตัวอย่าง (sampling design) เป็นคนอิสระต่อกัน หรือเป็นชุดที่จับคู่ซึ่งประกอบด้วย case หนึ่งคนกับ control ของมัน

logit ปรากฏในตำแหน่งต่างกันเล็กน้อยในแต่ละแบบจำลอง logistic regression สร้างแบบจำลอง logit ของความน่าจะเป็นหนึ่งค่า ordinal logistic regression สร้างแบบจำลอง logit หนึ่งค่าต่อจุดตัดระหว่างระดับแต่ละจุด ส่วน conditional logistic regression สร้างแบบจำลอง logit โดยให้ intercept แยกกันสำหรับชุดที่จับคู่แต่ละชุด แล้วกำจัด intercept เหล่านั้นออก

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

birth cohort จำลอง

cohort มีทารกแรกเกิด 500 คน โดย 200 คนเกิดจากมารดาที่ฝากครรภ์ไม่เพียงพอ และ 120 คนคลอดก่อนกำหนด (preterm birth) ซึ่งบางคนอยู่ในทั้งสองกลุ่ม ทารกแต่ละคนมีระดับ $Y$ ภายใน 28 วัน คือ 0 ไม่มี 1 เล็กน้อย 2 ปานกลาง 3 รุนแรง โดยรวมแล้วทารก 256 คนไม่มีภาวะแทรกซ้อน 101 คนมีภาวะเล็กน้อย 71 คนปานกลาง และ 72 คนรุนแรง

จำนวนเหล่านี้สร้างจากแบบจำลองที่ทราบค่า จึงตรวจสอบค่าประมาณทุกค่ากับแบบจำลองนั้นได้ สำหรับจุดตัดแต่ละจุดระหว่างระดับ แบบจำลองกำหนด log odds ของการอยู่ที่ระดับ $j$ หรือแย่กว่า:

$$\operatorname{logit} P(Y \ge j) = \alpha_j + \beta x + \gamma_j z$$

ในที่นี้ $x$ เท่ากับ 1 เมื่อฝากครรภ์ไม่เพียงพอ $z$ เท่ากับ 1 เมื่อคลอดก่อนกำหนด และ $j = 1, 2, 3$ ระบุจุดตัดทั้งสามจุด intercept คือ $\alpha_1 = -0.4$, $\alpha_2 = -1.4$ และ $\alpha_3 = -2.5$ ความชันของ ANC คือ $\beta = \ln 2$ ซึ่งเป็น odds ratio 2.00 ที่ทุกจุดตัด ความชันของการคลอดก่อนกำหนดคือ $\gamma_1 = 0.3$, $\gamma_2 = 0.7$ และ $\gamma_3 = 1.2$ ซึ่งเป็น odds ratio 1.35, 2.01 และ 3.32 เพิ่มขึ้นตามความรุนแรง

ความจริงจำลองนี้ให้ ANC มี odds ratio ค่าเดียวที่ทุกจุดตัด และให้การคลอดก่อนกำหนดมี odds ratio ต่างกันที่แต่ละจุด จำนวนที่ใช้เป็นจำนวนคาดหวังภายใต้แบบจำลอง ปัดเป็นทารกจำนวนเต็ม ไม่ใช่การสุ่ม ค่า odds ratio ที่ประมาณได้ด้านล่างจึงใกล้ค่าที่ใช้สร้างข้อมูล

Stata: สร้างจำนวนของ cohort

โค้ด Stata c7_sim.do (บรรทัด 53-61 จาก 265)
* expected number at or above each cut, rounded to whole newborns
forvalues j = 1/3 {
    generate int atleast`j' = floor(newborns*invlogit(`alpha`j'' + ln(2)*inadequate_anc + `gamma`j''*preterm) + 0.5)
}
generate int sev0 = newborns - atleast1
generate int sev1 = atleast1 - atleast2
generate int sev2 = atleast2 - atleast3
generate int sev3 = atleast3
list inadequate_anc preterm newborns sev0 sev1 sev2 sev3, noobs
ผลลัพธ์จากการรัน c7_sim.log
. * expected number at or above each cut, rounded to whole newborns
. forvalues j = 1/3 {
  2.     generate int atleast`j' = floor(newborns*invlogit(`alpha`j'' + ln(2)*inadequate_anc + `gamma`j''*preterm) + 0.5)
  3. }

. generate int sev0 = newborns - atleast1

. generate int sev1 = atleast1 - atleast2

. generate int sev2 = atleast2 - atleast3

. generate int sev3 = atleast3

. list inadequate_anc preterm newborns sev0 sev1 sev2 sev3, noobs

  +-----------------------------------------------------------+
  | inadeq~c   preterm   newborns   sev0   sev1   sev2   sev3 |
  |-----------------------------------------------------------|
  |        0         0        240    144     49     29     18 |
  |        0         1         60     31      9      7     13 |
  |        1         0        140     60     34     26     20 |
  |        1         1         60     21      9      9     21 |
  +-----------------------------------------------------------+
ข้อมูลจำลอง ทั้งสองส่วนเป็นข้อความที่ตัดมาจากสคริปต์จำลองเบื้องหลังชุดบทความนี้และผลลัพธ์ที่สคริปต์พิมพ์ออกมา ชื่อไฟล์และเลขบรรทัดบอกเพียงตำแหน่งของข้อความที่ตัดมาในสคริปต์นั้น และทุกอย่างที่ต้องใช้เพื่ออ่านก็แสดงไว้ครบแล้ว ลูปคำนวณจำนวนทารกที่คาดว่าจะอยู่ที่ระดับนั้นหรือสูงกว่าในสี่กลุ่มย่อยของ ANC และการคลอดก่อนกำหนด ปัดเป็นทารกจำนวนเต็ม คำสั่ง list พิมพ์จำนวนที่ได้ตามระดับ คือ sev0 ถึง sev3

แบบจำลองแรก: ผลลัพธ์สองค่าในทารกแรกเกิดที่เป็นอิสระต่อกัน

ชนิดของผลลัพธ์: สองค่า คือมีภาวะแทรกซ้อนใด ๆ (ระดับ 1 หรือแย่กว่า) เทียบกับไม่มี แบบแผนการสุ่มตัวอย่าง: ทารกแรกเกิดที่เป็นอิสระต่อกันจาก birth cohort เดียว แบบจำลอง: logistic regression ทั่วไป โดยมีการฝากครรภ์ไม่เพียงพอและการคลอดก่อนกำหนดเป็นตัวแปรร่วม การคลอดก่อนกำหนดอาจอยู่บนเส้นทางจาก ANC ไปสู่ภาวะแทรกซ้อนเอง ซึ่งจะกลับมากล่าวถึงใกล้ท้ายบทความ

จากทารก 500 คน มี 244 คนที่มีภาวะแทรกซ้อนระดับใดระดับหนึ่ง สัมประสิทธิ์ของ ANC ที่ประมาณได้เท่ากับ 0.691 บนมาตรวัด log odds เมื่อยกกำลังด้วยเลขชี้กำลังจะได้ odds ratio 2.00 (ช่วงเชื่อมั่น 95% 1.38 ถึง 2.88) ปรับด้วยการคลอดก่อนกำหนดแล้ว odds ratio ของการคลอดก่อนกำหนดคือ 1.40 (0.92 ถึง 2.13)

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

Stata: logit และ logistic กับภาวะแทรกซ้อนใด ๆ

โค้ด Stata c7_sim.do (บรรทัด 97-107 จาก 265)
generate byte any_morbidity = severity >= 1
quietly count if any_morbidity == 1
canon bin.n_any 0 r(N)
logit any_morbidity i.inadequate_anc i.preterm, nolog
canon logit.b_anc 4 _b[1.inadequate_anc]
canon logit.se_anc 4 _se[1.inadequate_anc]
canon logit.b_preterm 4 _b[1.preterm]
canon logit.se_preterm 4 _se[1.preterm]
canon logit.b_cons 4 _b[_cons]
canon logit.ll 4 e(ll)
logistic any_morbidity i.inadequate_anc i.preterm, nolog
ผลลัพธ์จากการรัน c7_sim.log
. generate byte any_morbidity = severity >= 1

. quietly count if any_morbidity == 1

. canon bin.n_any 0 r(N)
CANON c7.bin.n_any 244

. logit any_morbidity i.inadequate_anc i.preterm, nolog

Logistic regression                                     Number of obs =    500
                                                        LR chi2(2)    =  17.79
                                                        Prob > chi2   = 0.0001
Log likelihood = -337.53237                             Pseudo R2     = 0.0257

----------------------------------------------------------------------------------
   any_morbidity | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc |   .6912976   .1868748     3.70   0.000     .3250298    1.057565
       1.preterm |   .3354302   .2149308     1.56   0.119    -.0858265     .756687
           _cons |   -.404775   .1254744    -3.23   0.001    -.6507003   -.1588497
----------------------------------------------------------------------------------

. canon logit.b_anc 4 _b[1.inadequate_anc]
CANON c7.logit.b_anc 0.6913

. canon logit.se_anc 4 _se[1.inadequate_anc]
CANON c7.logit.se_anc 0.1869

. canon logit.b_preterm 4 _b[1.preterm]
CANON c7.logit.b_preterm 0.3354

. canon logit.se_preterm 4 _se[1.preterm]
CANON c7.logit.se_preterm 0.2149

. canon logit.b_cons 4 _b[_cons]
CANON c7.logit.b_cons -0.4048

. canon logit.ll 4 e(ll)
CANON c7.logit.ll -337.5324

. logistic any_morbidity i.inadequate_anc i.preterm, nolog

Logistic regression                                     Number of obs =    500
                                                        LR chi2(2)    =  17.79
                                                        Prob > chi2   = 0.0001
Log likelihood = -337.53237                             Pseudo R2     = 0.0257

----------------------------------------------------------------------------------
   any_morbidity | Odds ratio   Std. err.      z    P>|z|     [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc |   1.996304   .3730589     3.70   0.000     1.384072    2.879352
       1.preterm |   1.398542   .3005898     1.56   0.119     .9177535    2.131204
           _cons |   .6671269   .0837074    -3.23   0.001     .5216803    .8531246
----------------------------------------------------------------------------------
Note: _cons estimates baseline odds.
ข้อมูลจำลอง ข้อความที่ตัดมาจากสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ คำสั่ง logit พิมพ์สัมประสิทธิ์บนมาตรวัด log odds ส่วน logistic ประมาณแบบจำลองเดียวกันและพิมพ์ odds ratio หลังปัดเศษ แถว inadequate_anc อ่านได้ 0.691 ในตารางแรกและ 2.00 ในตารางที่สอง บรรทัดที่ขึ้นต้นด้วย CANON พิมพ์โดยฟังก์ชันช่วยเล็ก ๆ ในสคริปต์ที่สะท้อนผลลัพธ์แต่ละค่าออกมา ให้อ่านตัวเลขท้ายบรรทัดและไม่ต้องสนใจคำนำหน้า

R: glm() กับภาวะแทรกซ้อนใด ๆ

โค้ด R c7_sim_r.R (บรรทัด 77-89 จาก 222)
o$any_morbidity <- as.integer(o$severity >= 1)
canon("bin.n_any", sum(o$any_morbidity), 0)
fl <- glm(any_morbidity ~ inadequate_anc + preterm, family = binomial, data = o)
print(summary(fl)$coefficients)
se_l <- sqrt(diag(vcov(fl)))
canon("logit.b_anc", coef(fl)[["inadequate_anc"]])
canon("logit.se_anc", se_l[["inadequate_anc"]])
canon("logit.b_preterm", coef(fl)[["preterm"]])
canon("logit.se_preterm", se_l[["preterm"]])
canon("logit.b_cons", coef(fl)[["(Intercept)"]])
canon("logit.ll", as.numeric(logLik(fl)))
# the same fit shown as odds ratios
print(exp(cbind(OR = coef(fl), confint.default(fl))))
ผลลัพธ์จากการรัน c7_sim_r.log
> o$any_morbidity <- as.integer(o$severity >= 1)

> canon("bin.n_any", sum(o$any_morbidity), 0)
CANON c7.bin.n_any 244

> fl <- glm(any_morbidity ~ inadequate_anc + preterm,
+     family = binomial, data = o)

> print(summary(fl)$coefficients)
                 Estimate Std. Error   z value     Pr(>|z|)
(Intercept)    -0.4047750  0.1254744 -3.225956 0.0012555244
inadequate_anc  0.6912976  0.1868748  3.699256 0.0002162321
preterm         0.3354302  0.2149308  1.560643 0.1186080947

> se_l <- sqrt(diag(vcov(fl)))

> canon("logit.b_anc", coef(fl)[["inadequate_anc"]])
CANON c7.logit.b_anc 0.6913

> canon("logit.se_anc", se_l[["inadequate_anc"]])
CANON c7.logit.se_anc 0.1869

> canon("logit.b_preterm", coef(fl)[["preterm"]])
CANON c7.logit.b_preterm 0.3354

> canon("logit.se_preterm", se_l[["preterm"]])
CANON c7.logit.se_preterm 0.2149

> canon("logit.b_cons", coef(fl)[["(Intercept)"]])
CANON c7.logit.b_cons -0.4048

> canon("logit.ll", as.numeric(logLik(fl)))
CANON c7.logit.ll -337.5324

> print(exp(cbind(OR = coef(fl), confint.default(fl))))
                      OR     2.5 %    97.5 %
(Intercept)    0.6671269 0.5216803 0.8531246
inadequate_anc 1.9963043 1.3840720 2.8793525
preterm        1.3985420 0.9177535 2.1312037
ข้อมูลจำลอง ข้อความที่ตัดมาจากสคริปต์ R ของสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ glm() ที่ใช้ family = binomial ใช้ logit link ตารางสัมประสิทธิ์ของมันจึงอยู่บนมาตรวัด log odds การยกกำลังสัมประสิทธิ์และช่วง 95% ของมันด้วยเลขชี้กำลังให้ odds ratio ที่ logistic ของ Stata พิมพ์ บรรทัดที่ขึ้นต้นด้วย CANON พิมพ์โดยฟังก์ชันช่วยเล็ก ๆ ในสคริปต์ที่สะท้อนผลลัพธ์แต่ละค่าออกมา ให้อ่านตัวเลขท้ายบรรทัดและไม่ต้องสนใจคำนำหน้า

แบบจำลองที่สอง: ระดับที่มีลำดับในทารกกลุ่มเดียวกัน

ชนิดของผลลัพธ์: หมวดที่มีลำดับ คือระดับภาวะแทรกซ้อนสี่ระดับ แบบแผนการสุ่มตัวอย่าง: ทารกแรกเกิดอิสระต่อกันกลุ่มเดิม แบบจำลอง: การถดถอยโลจิสติกสำหรับผลลัพธ์แบบลำดับ (ordinal logistic regression) ซึ่งใช้ทุกระดับพร้อมกันและรักษาลำดับของมันไว้ [2]

ผลลัพธ์สี่ระดับมี จุดตัด (cut point) สามจุด คือเส้นแบ่งระหว่างระดับที่อยู่ติดกัน จุดตัดแต่ละจุดแบ่งทารกเป็นสองกลุ่ม คือระดับ $j$ หรือแย่กว่า กับดีกว่าระดับ $j$ log odds ของการอยู่ที่ระดับ $j$ หรือแย่กว่าคือ cumulative logit ของจุดตัดนั้น

ordinal logistic regression ประมาณ cumulative logit หนึ่งค่าต่อจุดตัด แต่ละค่ามี intercept ของตัวเอง ทว่าใช้ชุดความชันร่วมกันชุดเดียว:

$$\operatorname{logit} P(Y \ge j) = \alpha_j + \beta x + \gamma z$$

intercept $\alpha_j$ ลดลงเมื่อ $j$ เพิ่ม เพราะทารกที่ถึงระดับสูงขึ้นมีน้อยลง ความชัน $\beta$ และ $\gamma$ ไม่มี $j$ กำกับ odds ratio ของ ANC คือ $e^{\beta}$ จึงเท่ากันทุกจุดตัด นี่คือ ข้อสมมติ proportional odds (proportional odds assumption) ต่างจากแบบจำลองที่ใช้สร้างข้อมูล แบบจำลองนี้บังคับให้การคลอดก่อนกำหนดมีความชันค่าเดียวที่ทุกจุดตัด

คำสั่ง ologit ของ Stata และ polr() ของ R พิมพ์ intercept ด้วยเครื่องหมายตรงข้าม เป็นจุดตัด $\kappa_j = -\alpha_j$ ที่เพิ่มขึ้นตาม $j$ ความชันคงเครื่องหมายเดิม สัมประสิทธิ์บวกทำให้ระดับสูงมีโอกาสมากขึ้น

ใน cohort จำลอง odds ratio ของ ANC ที่ปรับด้วยการคลอดก่อนกำหนดคือ 2.01 (ช่วงเชื่อมั่น 95% 1.43 ถึง 2.82) ใกล้เคียง 2.00 ที่ใช้สร้างข้อมูล เป็น odds ratio ของการอยู่ในระดับที่สูงขึ้น ไม่ว่าจะลากเส้นแบ่งระหว่างระดับตรงไหน ช่วงของมันแคบกว่า 1.38 ถึง 2.88 ของแบบจำลองสองค่า เพราะแบบจำลองแบบลำดับใช้ระดับที่อยู่เหนือจุดตัดแรกด้วย

อ่านให้เป็นรูปธรรม odds ratio นี้ใช้กับทุกวิธีแบ่งระดับพร้อมกัน ในทารกที่มีสถานะคลอดก่อนกำหนดเหมือนกัน การฝากครรภ์ไม่เพียงพอเพิ่ม odds ของภาวะแทรกซ้อนใด ๆ ประมาณสองเท่า ของภาวะปานกลางหรือรุนแรง และของภาวะรุนแรง เช่นเดียวกัน การยึดสถานะการคลอดก่อนกำหนดไว้ให้เท่ากันเช่นนี้อาจละส่วนหนึ่งของผลของ ANC ไป ดังที่หัวข้อ สิ่งที่ odds ratio ทั้งสามมีร่วมกัน อธิบาย

เมื่อผลลัพธ์พบบ่อย การเพิ่ม odds เป็นสองเท่าไม่ได้เพิ่มความเสี่ยงเป็นสองเท่า ดังที่ความน่าจะเป็นด้านล่างแสดง จุดตัดคือ 0.458, 1.361 และ 2.267 พิจารณาทารกที่คลอดครบกำหนด ไม่ได้คลอดก่อนกำหนด และฝากครรภ์เพียงพอ แบบจำลองให้ทารกคนนี้มีความน่าจะเป็น 0.613 ที่จะไม่มีภาวะแทรกซ้อน และ 0.094 ที่จะมีภาวะรุนแรง

เมื่อฝากครรภ์ไม่เพียงพอ ค่าเหล่านี้กลายเป็น 0.441 และ 0.172 odds ratio ค่าเดียวจึงเลื่อนการกระจายของระดับทั้งหมดไปทางด้านรุนแรง การเลื่อนเพียงครั้งเดียวนี้คือสิ่งที่ข้อสมมติ proportional odds ให้ และเป็นสิ่งที่มันเสี่ยง

odds ratio ของการคลอดก่อนกำหนดคือ 1.84 (1.24 ถึง 2.74) ต่างจาก odds ratio ของ ANC ตรงที่มันสรุป odds ratio ซึ่งในความจริงจำลองเพิ่มจาก 1.35 ที่จุดตัดแรกเป็น 3.32 ที่จุดตัดที่สาม คำถามว่าตัวเลขเดียวสรุปได้เหมาะสมหรือไม่คือสิ่งที่การตรวจ proportional odds ถาม

Stata: ologit กับระดับทั้งสี่

โค้ด Stata c7_sim.do (บรรทัด 117-117 จาก 265)
ologit severity i.inadequate_anc i.preterm, nolog
ผลลัพธ์จากการรัน c7_sim.log
. ologit severity i.inadequate_anc i.preterm, nolog

Ordered logistic regression                             Number of obs =    500
                                                        LR chi2(2)    =  27.91
                                                        Prob > chi2   = 0.0000
Log likelihood = -597.08706                             Pseudo R2     = 0.0228

----------------------------------------------------------------------------------
        severity | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc |   .6966207   .1732433     4.02   0.000     .3570702    1.036171
       1.preterm |   .6124173   .2012158     3.04   0.002     .2180416    1.006793
-----------------+----------------------------------------------------------------
           /cut1 |   .4581421   .1220414                      .2189454    .6973388
           /cut2 |   1.360579   .1350555                      1.095875    1.625283
           /cut3 |   2.266819   .1627571                      1.947821    2.585817
----------------------------------------------------------------------------------
ข้อมูลจำลอง ข้อความที่ตัดมาจากสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ สัมประสิทธิ์ทั้งสองเป็น log odds ratio ที่ใช้ร่วมกันทั้งสามจุดตัด เมื่อยกกำลังแถว inadequate_anc คือ 0.697 จะได้ odds ratio 2.01 แถว /cut1 ถึง /cut3 คือจุดตัด $\kappa_j$

R: polr() กับระดับทั้งสี่

โค้ด R c7_sim_r.R (บรรทัด 95-97 จาก 222)
fo <- polr(factor(severity, levels = 0:3, ordered = TRUE) ~ inadequate_anc + preterm, data = o,
           method = "logistic", Hess = TRUE)
print(summary(fo))
ผลลัพธ์จากการรัน c7_sim_r.log
> fo <- polr(factor(severity, levels = 0:3, ordered = TRUE) ~
+     inadequate_anc + preterm, data = o, method = "logistic",
+     Hess = TRUE)

> print(summary(fo))
Call:
polr(formula = factor(severity, levels = 0:3, ordered = TRUE) ~
    inadequate_anc + preterm, data = o, Hess = TRUE, method = "logistic")

Coefficients:
                Value Std. Error t value
inadequate_anc 0.6966     0.1732   4.021
preterm        0.6124     0.2012   3.044

Intercepts:
    Value   Std. Error t value
0|1  0.4581  0.1220     3.7540
1|2  1.3606  0.1351    10.0742
2|3  2.2668  0.1628    13.9276

Residual Deviance: 1194.174
AIC: 1204.174
ข้อมูลจำลอง ข้อความที่ตัดมาจากสคริปต์ R ของสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ polr() จากแพ็กเกจ MASS รับผลลัพธ์เป็น ordered factor สัมประสิทธิ์ของมันตรงกับ Stata และ intercept 0|1, 1|2 และ 2|3 คือจุดตัดสามจุดเดียวกัน มันพิมพ์ค่า t และไม่พิมพ์ค่า P

ตรวจ proportional odds ทีละจุดตัด

ตรวจข้อสมมติ proportional odds ได้โดยตรง ประมาณ logistic regression ทั่วไปที่จุดตัดแต่ละจุดทั้งสามจุด สามแบบจำลองเปรียบเทียบระดับ 1 หรือแย่กว่ากับไม่มี ระดับ 2 หรือแย่กว่ากับต่ำกว่านั้น และระดับ 3 กับต่ำกว่านั้น ถ้าข้อสมมติเป็นจริง odds ratio ของตัวแปรร่วมแต่ละตัวควรใกล้เคียงกันทั้งสามจุด

ตารางด้านล่างแสดงว่า ANC เป็นไปตามที่ข้อสมมติต้องการ โดยมี odds ratio 2.00, 2.01 และ 2.01 ที่สามจุดตัด การคลอดก่อนกำหนดไม่เป็นเช่นนั้น คือ 1.40, 2.05 และ 3.31 ใกล้ความจริงจำลองที่ 1.35, 2.01 และ 3.32 แบบจำลองที่จุดตัดแรกคือแบบจำลองสองค่าข้างต้น ตัวเลขจึงตรงกับแบบจำลองนั้น

ข้อมูลจำลอง ทารก 500 คน สามแถวแรกแต่ละแถวเป็น logistic regression แยกกันที่มีตัวแปรร่วมทั้งสองตัว แถวสุดท้ายคือแบบจำลองแบบลำดับ ความจริงจำลองคือแบบจำลองที่ใช้สร้างจำนวนเหล่านี้
จุดตัดodds ratio ของ ANC ค่าที่ประมาณได้odds ratio ของ ANC ความจริงจำลองodds ratio ของการคลอดก่อนกำหนด ค่าที่ประมาณได้odds ratio ของการคลอดก่อนกำหนด ความจริงจำลอง
ระดับ 1 หรือแย่กว่า เทียบกับไม่มี2.002.001.401.35
ระดับ 2 หรือแย่กว่า เทียบกับระดับ 0 ถึง 12.012.002.052.01
ระดับ 3 เทียบกับระดับ 0 ถึง 22.012.003.313.32
แบบจำลองแบบลำดับ odds ratio ค่าเดียวสำหรับทุกจุดตัด2.012.001.84ไม่มีค่าเดียว

Stata: logistic หนึ่งแบบจำลองต่อจุดตัด

โค้ด Stata c7_sim.do (บรรทัด 138-149 จาก 265)
* the proportional-odds check, step 1: a separate binary logit at each cut (severity >= 1, >= 2, >= 3);
* plain 0/1 covariates, so that e(b) holds inadequate_anc, preterm and the constant in that order
forvalues j = 1/3 {
    generate byte ge`j' = severity >= `j'
    logit ge`j' inadequate_anc preterm, nolog
    matrix bcut`j' = e(b)
    predict double pi`j', pr
    canon cut`j'.b_anc 4 _b[inadequate_anc]
    canon cut`j'.or_anc 4 exp(_b[inadequate_anc])
    canon cut`j'.b_preterm 4 _b[preterm]
    canon cut`j'.or_preterm 4 exp(_b[preterm])
}
ผลลัพธ์จากการรัน c7_sim.log
. * the proportional-odds check, step 1: a separate binary logit at each cut (severity >= 1, >= 2, >= 3);
. * plain 0/1 covariates, so that e(b) holds inadequate_anc, preterm and the constant in that order
. forvalues j = 1/3 {
  2.     generate byte ge`j' = severity >= `j'
  3.     logit ge`j' inadequate_anc preterm, nolog
  4.     matrix bcut`j' = e(b)
  5.     predict double pi`j', pr
  6.     canon cut`j'.b_anc 4 _b[inadequate_anc]
  7.     canon cut`j'.or_anc 4 exp(_b[inadequate_anc])
  8.     canon cut`j'.b_preterm 4 _b[preterm]
  9.     canon cut`j'.or_preterm 4 exp(_b[preterm])
 10. }

Logistic regression                                     Number of obs =    500
                                                        LR chi2(2)    =  17.79
                                                        Prob > chi2   = 0.0001
Log likelihood = -337.53237                             Pseudo R2     = 0.0257

--------------------------------------------------------------------------------
           ge1 | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc |   .6912976   .1868748     3.70   0.000     .3250298    1.057565
       preterm |   .3354302   .2149308     1.56   0.119    -.0858265     .756687
         _cons |   -.404775   .1254744    -3.23   0.001    -.6507003   -.1588497
--------------------------------------------------------------------------------
CANON c7.cut1.b_anc 0.6913
CANON c7.cut1.or_anc 1.9963
CANON c7.cut1.b_preterm 0.3354
CANON c7.cut1.or_preterm 1.3985

Logistic regression                                     Number of obs =    500
                                                        LR chi2(2)    =  24.29
                                                        Prob > chi2   = 0.0000
Log likelihood = -287.11966                             Pseudo R2     = 0.0406

--------------------------------------------------------------------------------
           ge2 | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc |   .6964978   .2038964     3.42   0.001     .2968681    1.096127
       preterm |   .7169957   .2238947     3.20   0.001     .2781702    1.155821
         _cons |  -1.411917   .1515253    -9.32   0.000    -1.708901   -1.114933
--------------------------------------------------------------------------------
CANON c7.cut2.b_anc 0.6965
CANON c7.cut2.or_anc 2.0067
CANON c7.cut2.b_preterm 0.7170
CANON c7.cut2.or_preterm 2.0483

Logistic regression                                     Number of obs =    500
                                                        LR chi2(2)    =  29.04
                                                        Prob > chi2   = 0.0000
Log likelihood = -191.55985                             Pseudo R2     = 0.0705

--------------------------------------------------------------------------------
           ge3 | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc |   .6983656   .2648017     2.64   0.008     .1793638    1.217367
       preterm |     1.1974   .2683201     4.46   0.000     .6715027    1.723298
         _cons |  -2.501102   .2175031   -11.50   0.000      -2.9274   -2.074803
--------------------------------------------------------------------------------
CANON c7.cut3.b_anc 0.6984
CANON c7.cut3.or_anc 2.0105
CANON c7.cut3.b_preterm 1.1974
CANON c7.cut3.or_preterm 3.3115
ข้อมูลจำลอง ข้อความที่ตัดมาจากสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ ลูปประมาณ logistic regression หนึ่งแบบจำลองต่อจุดตัด คือ ge1 ถึง ge3 และเก็บสัมประสิทธิ์กับความน่าจะเป็นที่ทำนายของแต่ละแบบจำลองไว้ใช้กับ Brant test ถัดไป บรรทัดที่ขึ้นต้นด้วย CANON พิมพ์โดยฟังก์ชันช่วยเล็ก ๆ ในสคริปต์ที่สะท้อนผลลัพธ์แต่ละค่าออกมา ให้อ่านตัวเลขท้ายบรรทัดและไม่ต้องสนใจคำนำหน้า

Brant test: ความชันเท่ากันทุกจุดตัดหรือไม่

การเปรียบเทียบ odds ratio สามค่าด้วยตาไม่ใช่การทดสอบ เพราะค่าที่ประมาณแต่ละค่ามีความคลาดเคลื่อนจากการสุ่มตัวอย่าง Brant test ถามว่าความชันของแบบจำลองแยกตามจุดตัดต่างกันเกินกว่าความคลาดเคลื่อนนั้นจะอธิบายได้หรือไม่ [3] เป็น การทดสอบ Wald (Wald test) คือชั่งผลต่างของความชันเทียบกับความแปรปรวนจากการสุ่มตัวอย่างที่ประมาณได้ และให้สถิติไคสแควร์

แบบจำลองแยกใช้ทารกกลุ่มเดียวกัน ความชันของมันจึงสัมพันธ์กัน วิธีของ Brant ประมาณความสัมพันธ์นั้นจากแบบจำลองทั้งสามจุดตัด ขณะที่การเปรียบเทียบแบบตรงไปตรงมาจะมองข้ามมัน [3] ทดสอบได้ทั้งตัวแปรร่วมทั้งหมดพร้อมกันหรือทีละตัว องศาอิสระ (degrees of freedom) นับผลต่างของความชันที่ทดสอบ ในที่นี้คือสองค่าต่อตัวแปรร่วมหนึ่งตัว

สำหรับตัวแปรร่วมทั้งหมดรวมกัน สถิติเท่ากับ 11.72 ที่องศาอิสระ 4 (P = 0.020) การทดสอบพบหลักฐานต่อข้อสมมติสำหรับแบบจำลองโดยรวม สำหรับการคลอดก่อนกำหนดเพียงอย่างเดียว สถิติเท่ากับ 11.65 ที่องศาอิสระ 2 (P = 0.003)

สำหรับ ANC เพียงอย่างเดียว สถิติเท่ากับ 0.001 ที่องศาอิสระ 2 (P = 0.999) เพราะจำนวนถูกสร้างโดยไม่มีสัญญาณรบกวนจากการสุ่มตัวอย่าง สถิติของ ANC จึงอยู่เกือบที่ศูนย์ ในตัวอย่างจริงที่ข้อสมมติเป็นจริง สถิติแปรผันไปตามแต่ละตัวอย่าง และค่า P อาจตกที่ใดก็ได้ระหว่าง 0 ถึง 1

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

Stata: Brant test ที่คำนวณใน Mata

โค้ด Stata c7_sim.do (บรรทัด 150-184 จาก 265)
* step 2, the Brant (1990) Wald test that the slopes are equal across the three cuts. For cuts j and l the
* covariance of the two slope vectors is A_j^-1 [X' diag(pi_hi (1 - pi_lo)) X] A_l^-1, with
* A_j = X' diag(pi_j (1 - pi_j)) X, lo = min(j, l) and hi = max(j, l); the contrasts compare cut 1 with
* cuts 2 and 3 (4 df for both covariates, 2 df for each covariate)
mata:
X = st_data(., ("inadequate_anc", "preterm")), J(st_nobs(), 1, 1)
P = st_data(., ("pi1", "pi2", "pi3"))
B1 = st_matrix("bcut1")
B2 = st_matrix("bcut2")
B3 = st_matrix("bcut3")
b = (B1[|1, 1 \ 1, 2|], B2[|1, 1 \ 1, 2|], B3[|1, 1 \ 1, 2|])'
V = J(6, 6, 0)
for (j = 1; j <= 3; j++) {
    for (l = 1; l <= 3; l++) {
        lo = min((j, l))
        hi = max((j, l))
        Ajl = invsym(cross(X, P[., j]:*(1 :- P[., j]), X)) * cross(X, P[., hi]:*(1 :- P[., lo]), X) * invsym(cross(X, P[., l]:*(1 :- P[., l]), X))
        V[|2*j - 1, 2*l - 1 \ 2*j, 2*l|] = Ajl[|1, 1 \ 2, 2|]
    }
}
D = (1, 0, -1, 0, 0, 0 \ 0, 1, 0, -1, 0, 0 \ 1, 0, 0, 0, -1, 0 \ 0, 1, 0, 0, 0, -1)
Da = D[(1 \ 3), .]
Dp = D[(2 \ 4), .]
st_numscalar("sc_brant_all", (D*b)' * invsym(D*V*D') * (D*b))
st_numscalar("sc_brant_anc", (Da*b)' * invsym(Da*V*Da') * (Da*b))
st_numscalar("sc_brant_pre", (Dp*b)' * invsym(Dp*V*Dp') * (Dp*b))
end
canon brant.chi2_all 4 scalar(sc_brant_all)
canon brant.df_all 0 4
canon brant.p_all 4 chi2tail(4, scalar(sc_brant_all))
canon brant.chi2_anc 4 scalar(sc_brant_anc)
canon brant.p_anc 4 chi2tail(2, scalar(sc_brant_anc))
canon brant.chi2_preterm 4 scalar(sc_brant_pre)
canon brant.p_preterm 4 chi2tail(2, scalar(sc_brant_pre))
canon brant.df_each 0 2
ผลลัพธ์จากการรัน c7_sim.log
. * step 2, the Brant (1990) Wald test that the slopes are equal across the three cuts. For cuts j and l the
. * covariance of the two slope vectors is A_j^-1 [X' diag(pi_hi (1 - pi_lo)) X] A_l^-1, with
. * A_j = X' diag(pi_j (1 - pi_j)) X, lo = min(j, l) and hi = max(j, l); the contrasts compare cut 1 with
. * cuts 2 and 3 (4 df for both covariates, 2 df for each covariate)
. mata:
------------------------------------------------- mata (type end to exit) -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
: X = st_data(., ("inadequate_anc", "preterm")), J(st_nobs(), 1, 1)

: P = st_data(., ("pi1", "pi2", "pi3"))

: B1 = st_matrix("bcut1")

: B2 = st_matrix("bcut2")

: B3 = st_matrix("bcut3")

: b = (B1[|1, 1 \ 1, 2|], B2[|1, 1 \ 1, 2|], B3[|1, 1 \ 1, 2|])'

: V = J(6, 6, 0)

: for (j = 1; j <= 3; j++) {
>     for (l = 1; l <= 3; l++) {
>         lo = min((j, l))
>         hi = max((j, l))
>         Ajl = invsym(cross(X, P[., j]:*(1 :- P[., j]), X)) * cross(X, P[., hi]:*(1 :- P[., lo]), X) * invsym(cross(X, P[., l]:*(1 :- P[., l]), X))
>         V[|2*j - 1, 2*l - 1 \ 2*j, 2*l|] = Ajl[|1, 1 \ 2, 2|]
>     }
> }

: D = (1, 0, -1, 0, 0, 0 \ 0, 1, 0, -1, 0, 0 \ 1, 0, 0, 0, -1, 0 \ 0, 1, 0, 0, 0, -1)

: Da = D[(1 \ 3), .]

: Dp = D[(2 \ 4), .]

: st_numscalar("sc_brant_all", (D*b)' * invsym(D*V*D') * (D*b))

: st_numscalar("sc_brant_anc", (Da*b)' * invsym(Da*V*Da') * (Da*b))

: st_numscalar("sc_brant_pre", (Dp*b)' * invsym(Dp*V*Dp') * (Dp*b))

: end
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

. canon brant.chi2_all 4 scalar(sc_brant_all)
CANON c7.brant.chi2_all 11.7179

. canon brant.df_all 0 4
CANON c7.brant.df_all 4

. canon brant.p_all 4 chi2tail(4, scalar(sc_brant_all))
CANON c7.brant.p_all 0.0196

. canon brant.chi2_anc 4 scalar(sc_brant_anc)
CANON c7.brant.chi2_anc 0.0011

. canon brant.p_anc 4 chi2tail(2, scalar(sc_brant_anc))
CANON c7.brant.p_anc 0.9994

. canon brant.chi2_preterm 4 scalar(sc_brant_pre)
CANON c7.brant.chi2_preterm 11.6476

. canon brant.p_preterm 4 chi2tail(2, scalar(sc_brant_pre))
CANON c7.brant.p_preterm 0.0030

. canon brant.df_each 0 2
CANON c7.brant.df_each 2
ข้อมูลจำลอง ข้อความที่ตัดมาจากสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ Mata คือภาษาเมทริกซ์ของ Stata บล็อกนี้ซ้อนความชันหกค่าของแบบจำลองสามจุดตัด สร้างความแปรปรวนร่วมข้ามจุดตัดด้วยสูตรของ Brant และเปรียบเทียบจุดตัดแรกกับจุดตัดที่สอง และกับจุดตัดที่สาม ได้สถิติ Wald สำหรับความชันทั้งหมด สำหรับ ANC เพียงอย่างเดียว และสำหรับการคลอดก่อนกำหนดเพียงอย่างเดียว บรรทัดที่ขึ้นต้นด้วย CANON พิมพ์โดยฟังก์ชันช่วยเล็ก ๆ ในสคริปต์ที่สะท้อนผลลัพธ์แต่ละค่าออกมา ให้อ่านตัวเลขท้ายบรรทัดและไม่ต้องสนใจคำนำหน้า

R: แบบจำลองแยกตามจุดตัดและ Brant test

โค้ด R c7_sim_r.R (บรรทัด 117-162 จาก 222)
# the proportional-odds check, step 1: a separate binary logit at each cut (severity >= 1, >= 2, >= 3)
fits <- lapply(1:3, function(j) {
  o$ge <- as.integer(o$severity >= j)
  glm(ge ~ inadequate_anc + preterm, family = binomial, data = o)
})
for (j in 1:3) {
  f <- fits[[j]]
  canon(sprintf("cut%d.b_anc", j), coef(f)[["inadequate_anc"]])
  canon(sprintf("cut%d.or_anc", j), exp(coef(f)[["inadequate_anc"]]))
  canon(sprintf("cut%d.b_preterm", j), coef(f)[["preterm"]])
  canon(sprintf("cut%d.or_preterm", j), exp(coef(f)[["preterm"]]))
}
# step 2, the Brant (1990) Wald test that the slopes are equal across the three cuts. For cuts j and l the
# covariance of the two slope vectors is A_j^-1 [X' diag(pi_hi (1 - pi_lo)) X] A_l^-1, with
# A_j = X' diag(pi_j (1 - pi_j)) X, lo = min(j, l) and hi = max(j, l); the contrasts compare cut 1 with
# cuts 2 and 3 (4 df for both covariates, 2 df for each covariate)
X <- cbind(inadequate_anc = o$inadequate_anc, preterm = o$preterm, cons = 1)
P <- sapply(fits, fitted)
b <- unlist(lapply(fits, function(f) coef(f)[c("inadequate_anc", "preterm")]))
V <- matrix(0, 6, 6)
for (j in 1:3) {
  for (l in 1:3) {
    lo <- min(j, l)
    hi <- max(j, l)
    Ajl <- solve(crossprod(X, X * (P[, j] * (1 - P[, j])))) %*%
      crossprod(X, X * (P[, hi] * (1 - P[, lo]))) %*%
      solve(crossprod(X, X * (P[, l] * (1 - P[, l]))))
    V[(2 * j - 1):(2 * j), (2 * l - 1):(2 * l)] <- Ajl[1:2, 1:2]
  }
}
D <- rbind(c(1, 0, -1, 0, 0, 0), c(0, 1, 0, -1, 0, 0), c(1, 0, 0, 0, -1, 0), c(0, 1, 0, 0, 0, -1))
wald <- function(Dm) {
  db <- Dm %*% b
  as.numeric(t(db) %*% solve(Dm %*% V %*% t(Dm)) %*% db)
}
chi_all <- wald(D)
chi_anc <- wald(D[c(1, 3), ])
chi_pre <- wald(D[c(2, 4), ])
canon("brant.chi2_all", chi_all)
canon("brant.df_all", 4, 0)
canon("brant.p_all", pchisq(chi_all, 4, lower.tail = FALSE))
canon("brant.chi2_anc", chi_anc)
canon("brant.p_anc", pchisq(chi_anc, 2, lower.tail = FALSE))
canon("brant.chi2_preterm", chi_pre)
canon("brant.p_preterm", pchisq(chi_pre, 2, lower.tail = FALSE))
canon("brant.df_each", 2, 0)
ผลลัพธ์จากการรัน c7_sim_r.log
> fits <- lapply(1:3, function(j) {
+     o$ge <- as.integer(o$severity >= j)
+     glm(ge ~ inadequate_anc + preterm, family = binomial, data = o)
+ })

> for (j in 1:3) {
+     f <- fits[[j]]
+     canon(sprintf("cut%d.b_anc", j), coef(f)[["inadequate_anc"]])
+     canon(sprintf("cut%d.or_anc", j), exp(coef(f)[["inadequate_anc"]]))
+     canon(sprintf("cut%d.b_preterm", j), coef(f)[["preterm"]])
+     canon(sprintf("cut%d.or_preterm", j), exp(coef(f)[["preterm"]]))
+ }
CANON c7.cut1.b_anc 0.6913
CANON c7.cut1.or_anc 1.9963
CANON c7.cut1.b_preterm 0.3354
CANON c7.cut1.or_preterm 1.3985
CANON c7.cut2.b_anc 0.6965
CANON c7.cut2.or_anc 2.0067
CANON c7.cut2.b_preterm 0.7170
CANON c7.cut2.or_preterm 2.0483
CANON c7.cut3.b_anc 0.6984
CANON c7.cut3.or_anc 2.0105
CANON c7.cut3.b_preterm 1.1974
CANON c7.cut3.or_preterm 3.3115

> X <- cbind(inadequate_anc = o$inadequate_anc, preterm = o$preterm,
+     cons = 1)

> P <- sapply(fits, fitted)

> b <- unlist(lapply(fits, function(f) coef(f)[c("inadequate_anc",
+     "preterm")]))

> V <- matrix(0, 6, 6)

> for (j in 1:3) {
+     for (l in 1:3) {
+         lo <- min(j, l)
+         hi <- max(j, l)
+         Ajl <- solve(crossprod(X, X * (P[, j] * (1 - P[, j])))) %*%
+             crossprod(X, X * (P[, hi] * (1 - P[, lo]))) %*% solve(crossprod(X,
+             X * (P[, l] * (1 - P[, l]))))
+         V[(2 * j - 1):(2 * j), (2 * l - 1):(2 * l)] <- Ajl[1:2,
+             1:2]
+     }
+ }

> D <- rbind(c(1, 0, -1, 0, 0, 0), c(0, 1, 0, -1, 0,
+     0), c(1, 0, 0, 0, -1, 0), c(0, 1, 0, 0, 0, -1))

> wald <- function(Dm) {
+     db <- Dm %*% b
+     as.numeric(t(db) %*% solve(Dm %*% V %*% t(Dm)) %*% db)
+ }

> chi_all <- wald(D)

> chi_anc <- wald(D[c(1, 3), ])

> chi_pre <- wald(D[c(2, 4), ])

> canon("brant.chi2_all", chi_all)
CANON c7.brant.chi2_all 11.7179

> canon("brant.df_all", 4, 0)
CANON c7.brant.df_all 4

> canon("brant.p_all", pchisq(chi_all, 4, lower.tail = FALSE))
CANON c7.brant.p_all 0.0196

> canon("brant.chi2_anc", chi_anc)
CANON c7.brant.chi2_anc 0.0011

> canon("brant.p_anc", pchisq(chi_anc, 2, lower.tail = FALSE))
CANON c7.brant.p_anc 0.9994

> canon("brant.chi2_preterm", chi_pre)
CANON c7.brant.chi2_preterm 11.6476

> canon("brant.p_preterm", pchisq(chi_pre, 2, lower.tail = FALSE))
CANON c7.brant.p_preterm 0.0030

> canon("brant.df_each", 2, 0)
CANON c7.brant.df_each 2
ข้อมูลจำลอง ข้อความที่ตัดมาจากสคริปต์ R ของสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ glm() สามแบบจำลองคือแบบจำลองแยกตามจุดตัด และขั้นตอนเมทริกซ์ทำซ้ำการคำนวณของ Stata สถิติและค่า P ตรงกับของ Stata บรรทัดที่ขึ้นต้นด้วย CANON พิมพ์โดยฟังก์ชันช่วยเล็ก ๆ ในสคริปต์ที่สะท้อนผลลัพธ์แต่ละค่าออกมา ให้อ่านตัวเลขท้ายบรรทัดและไม่ต้องสนใจคำนำหน้า

เมื่อ proportional odds ไม่เป็นจริง

ตัวแปรร่วมที่ถูกชี้ไม่ได้ทำให้แบบจำลองแบบลำดับใช้ไม่ได้ ในที่นี้การตรวจชี้เฉพาะการคลอดก่อนกำหนด odds ratio ของ ANC ที่ 2.01 ยังสรุป odds ratio ของจุดตัดสามค่าที่เกือบเท่ากัน ส่วน odds ratio ของการคลอดก่อนกำหนดที่ 1.84 ไม่เป็นเช่นนั้น เพราะอยู่ระหว่าง odds ratio ของจุดตัดที่ไล่จาก 1.40 ถึง 3.31

มีทางออกสามแบบที่ผ่อนข้อสมมตินี้ แต่ละแบบถูกระบุชื่อไว้ที่นี่และอธิบายในบทความที่ลิงก์ไว้ด้านล่าง partial proportional odds model ให้ตัวแปรร่วมที่ถูกชี้ ในที่นี้คือการคลอดก่อนกำหนด มี odds ratio ของตัวเองที่แต่ละจุดตัด generalized ordered logit model ปล่อยตัวแปรร่วมทุกตัวในแบบเดียวกัน multinomial logistic regression ทิ้งลำดับและเปรียบเทียบแต่ละระดับกับระดับอ้างอิง

ใน cohort นี้ รายงานที่สมเหตุสมผลแบบหนึ่งคือคง odds ratio แบบลำดับของ ANC ไว้ และให้การคลอดก่อนกำหนดมี odds ratio เฉพาะจุดตัดสามค่า partial proportional odds model ให้ได้ทั้งสองในการประมาณครั้งเดียว ไม่ว่าทางใด รายงานควรระบุว่าได้ทำการตรวจและพบอะไร

บทความแยกเรื่อง ordinal logistic regression และเรื่อง multinomial เทียบกับ ordinal regression อธิบายทางออกเหล่านี้ ไม่ว่าเลือกแบบใด การแสดง odds ratio เฉพาะจุดตัดควบคู่กับค่าสรุปอาจช่วยผู้อ่าน ค่า P ของ Brant ยังขึ้นกับขนาดตัวอย่างด้วย การศึกษาขนาดใหญ่อาจชี้ความเบี่ยงเบนเล็กน้อย และการศึกษาขนาดเล็กอาจมองข้ามความเบี่ยงเบนขนาดใหญ่

แบบจำลองที่สาม: คู่ case-control ที่จับคู่

ชนิดของผลลัพธ์: สองค่าอีกครั้ง คือ case หรือ control แบบแผนการสุ่มตัวอย่าง: ชุดที่จับคู่ แต่ละชุดคือ case หนึ่งคนกับ control ที่ถูกสุ่มมาให้ แบบจำลอง: conditional logistic regression ซึ่งเปรียบเทียบ case แต่ละคนกับ control ของมันเท่านั้น [4]

case ทั้ง 20 คนมีภาวะแทรกซ้อนภายใน 28 วันหลังเกิด control ของแต่ละคนคือทารกที่ยังไม่มีภาวะแทรกซ้อนในวันที่ case มี คู่จึงถูกจับคู่ตามอายุเป็นวัน และ control อาจกลายเป็น case ในภายหลังได้ การเลือก control จากผู้ที่ยังอยู่ในกลุ่มเสี่ยงเรียกว่า การสุ่มกลุ่มควบคุมจากกลุ่มเสี่ยง (risk-set sampling) อธิบายไว้ในบทความแยกเรื่องการเลือกกลุ่มควบคุม

case กับ control ของมันก่อให้เกิด ชุดที่จับคู่ (matched set) ในที่นี้คือหนึ่งคู่ ซึ่งมีรหัสระบุของตัวเองในข้อมูล คู่ 20 คู่ให้ 40 แถว หนึ่งแถวต่อสมาชิกหนึ่งคน สิ่งสัมผัสคือการฝากครรภ์ไม่เพียงพอ ซึ่งบันทึกไว้สำหรับสมาชิกทั้งสองคนของทุกคู่

ประเภทของการจับคู่อธิบายไว้ในบทความเรื่อง conditional logistic regression เช่นเดียวกับ overmatching คือการจับคู่ตามปัจจัยที่สัมพันธ์กับสิ่งสัมผัสแต่ไม่ใช่ตัวกวน ซึ่งอาจเสียความแม่นยำ

ตัวอย่างคำนวณด้วยมือ คู่ที่จับคู่ 20 คู่ จำแนกตาม ANC ของ case และของ control ของมัน
ANC ของ caseANC ของ controlจำนวนคู่ชนิดของคู่
ไม่เพียงพอไม่เพียงพอ8สอดคล้องกัน
ไม่เพียงพอเพียงพอ4ไม่สอดคล้องกัน
เพียงพอไม่เพียงพอ1ไม่สอดคล้องกัน
เพียงพอเพียงพอ7สอดคล้องกัน
ทุกคู่ทุกแบบ20สอดคล้องกัน 15 คู่ ไม่สอดคล้องกัน 5 คู่

คู่ที่สอดคล้องกันและคู่ที่ไม่สอดคล้องกัน

คู่ที่สอดคล้องกัน (concordant pair) คือคู่ที่ case และ control มีสถานะสิ่งสัมผัสเหมือนกัน มี 15 คู่ โดย 8 คู่สัมผัสทั้งสองคน และ 7 คู่ไม่สัมผัสทั้งสองคน คู่ที่ไม่สอดคล้องกัน (discordant pair) คือคู่ที่ทั้งสองต่างกัน มี 4 คู่ที่เฉพาะ case สัมผัส และ 1 คู่ที่เฉพาะ control สัมผัส การแบ่งนี้เป็นตัวขับเคลื่อนทุกอย่างที่ตามมา

ตัวอย่างคำนวณด้วยมือ คู่ 20 คู่ที่ยุบรวมเป็นตารางเดียว โดยไม่สนใจการจับคู่
ANCCaseControl
ไม่เพียงพอ129
เพียงพอ811
รวม2020

ตัวอย่างคำนวณด้วยมือ: รวมคู่เข้าด้วยกัน

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

logistic regression ทั่วไปไม่สนใจว่า control คนไหนอยู่กับ case คนไหน มองเห็น case 20 คนและ control 20 คน และรวมเป็นตารางข้างบน คือสิ่งสัมผัสจำแนกตามสถานะ case

  1. Case ที่สัมผัส

    \[ 8 + 4 = 12 \]

    คู่ที่สัมผัสทั้งสองคน บวกคู่ที่เฉพาะ case สัมผัส case อีก 7 + 1 = 8 คนฝากครรภ์เพียงพอ

  2. Control ที่สัมผัส

    \[ 8 + 1 = 9 \]

    คู่ที่สัมผัสทั้งสองคน บวกคู่ที่เฉพาะ control สัมผัส control อีก 7 + 4 = 11 คนฝากครรภ์เพียงพอ

  3. odds ของการสัมผัสใน case

    \[ \frac{12}{8} = 1.50 \]

    ฝากครรภ์ไม่เพียงพอเทียบกับเพียงพอใน case 20 คน

  4. odds ของการสัมผัสใน control

    \[ \frac{9}{11} = 0.818 \]

    odds เดียวกันใน control 20 คน

  5. odds ratio แบบไม่จับคู่

    \[ \mathrm{OR} = \frac{1.50}{0.818} = 1.83 \]

    อัตราส่วนของ odds การสัมผัสสองค่า เช่นเดียวกับตาราง case-control แบบไม่จับคู่ทั่วไป

ผลลัพธ์: เมื่อรวมคู่เข้าด้วยกัน odds ของการฝากครรภ์ไม่เพียงพอใน case สูงเป็น 1.83 เท่าของ control

Conditional logistic regression: เปรียบเทียบภายในคู่

conditional logistic regression ถามคำถามที่แคบกว่า เมื่อกำหนดว่าสมาชิกของคู่เป็น case เพียงคนเดียว สิ่งสัมผัสเปลี่ยนว่าสมาชิกคนใดเป็น case อย่างไร มีเพียงคู่ที่สมาชิกต่างกันในสิ่งสัมผัสเท่านั้นที่ตอบได้

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

คู่ที่ไม่สอดคล้องกัน 5 คู่มีข้อมูลทั้งหมดเกี่ยวกับสิ่งสัมผัส ใน 4 คู่ในนั้น สมาชิกที่สัมผัสเป็น case ส่วนอีก 1 คู่ สมาชิกที่ไม่สัมผัสเป็น case พิจารณาคู่แบบ 1:1 ที่มีสิ่งสัมผัสแบบสองค่าหนึ่งตัวและไม่มีตัวแปรร่วมอื่น เมื่อนั้น odds ratio แบบมีเงื่อนไข (conditional OR) หรือ odds ratio จับคู่ $\mathrm{OR}_c$ ก็คืออัตราส่วนของสองจำนวนนี้:

$$\mathrm{OR}_c = \frac{4}{1} = 4.00$$

ตัวเศษนับคู่ที่เฉพาะ case สัมผัส และตัวส่วนนับคู่ที่เฉพาะ control สัมผัส ตัวอย่างคำนวณด้วยมือถัดไปแสดงว่าอัตราส่วนนี้มาจากไหน [5]

ตัวอย่างคำนวณด้วยมือ: ทำไม intercept ของคู่จึงหักล้างกัน

ตัวอย่างคำนวณด้วยมือ ให้คู่ที่ $k$ มี intercept ของตัวเองคือ $\omega_k$ แทนทุกสิ่งที่สมาชิกของคู่มีร่วมกัน และมีความชันของสิ่งสัมผัสร่วมกันคือ $\beta$ เขียน $x_1$ แทนสิ่งสัมผัสของสมาชิกที่กลายเป็น case และ $x_2$ แทนของสมาชิกอีกคน ภาวะน่าจะเป็น (likelihood) คือความน่าจะเป็นของข้อมูลที่สังเกต อ่านเป็นฟังก์ชันของสัมประสิทธิ์ conditional likelihood ถามว่าสมาชิกคนที่ 1 ไม่ใช่คนที่ 2 เป็น case น่าจะเป็นเท่าใด เมื่อกำหนดว่าคู่มี case เพียงคนเดียว [5]

  1. ความน่าจะเป็นของหนึ่งคู่

    \[ P_k = \frac{e^{\omega_k + \beta x_1}}{e^{\omega_k + \beta x_1} + e^{\omega_k + \beta x_2}} \]

    ภายใต้แบบจำลองโลจิสติก $e^{\omega_k + \beta x}$ คือ odds ที่สมาชิกจะเป็น case เมื่อกำหนดว่าคู่มี case หนึ่งคน ความน่าจะเป็นของ case คือสัดส่วนของ odds ของสมาชิกทั้งสอง

  2. intercept ของคู่หักล้างกัน

    \[ P_k = \frac{e^{\beta x_1}}{e^{\beta x_1} + e^{\beta x_2}} \]

    การหารทั้งเศษและส่วนด้วย $e^{\omega_k}$ ทำให้ intercept หายไป จึงไม่ต้องประมาณสิ่งที่สมาชิกของคู่มีร่วมกัน

  3. คู่ที่สอดคล้องกัน

    \[ x_1 = x_2: \quad P_k = \frac{e^{\beta x_1}}{2\, e^{\beta x_1}} = 0.5 \]

    ไม่ว่า $\beta$ จะเป็นเท่าใด ความน่าจะเป็นคือ 0.5 คู่ที่สอดคล้องกัน 15 คู่จึงไม่มีข้อมูลเกี่ยวกับ $\beta$

  4. คู่ที่ไม่สอดคล้องกัน

    \[ q = \frac{e^{\beta}}{1 + e^{\beta}} \]

    เขียน $q$ แทนความน่าจะเป็นที่สมาชิกที่สัมผัสเป็น case สมาชิกที่ไม่สัมผัสเป็น case ด้วยความน่าจะเป็น $1 - q$

  5. conditional likelihood

    \[ L(\beta) \propto q^{4}\,(1 - q)^{1} \]

    ใน 4 คู่ที่ไม่สอดคล้องกัน สมาชิกที่สัมผัสเป็น case และใน 1 คู่ สมาชิกที่ไม่สัมผัสเป็น case คู่ที่สอดคล้องกันเพียงคูณ $L(\beta)$ ด้วยค่าคงที่

  6. ค่าสูงสุด

    \[ q = \frac{4}{5} = 0.80 \]

    ผลคูณ $q^{4}(1 - q)$ มีค่ามากที่สุดเมื่อ $q$ เท่ากับสัดส่วนของคู่ที่ไม่สอดคล้องกันซึ่งสมาชิกที่สัมผัสเป็น case เมื่อนั้น $1 - q = 0.20$

  7. กลับสู่ odds ratio

    \[ e^{\beta} = \frac{q}{1 - q} = \frac{0.80}{0.20} = 4.00 \]

    แก้ขั้นของคู่ที่ไม่สอดคล้องกันหา $e^{\beta}$ ได้ odds $q/(1 - q)$ ซึ่งคือ odds ratio จับคู่ $\mathrm{OR}_c$

ผลลัพธ์: conditional likelihood ให้ odds ratio จับคู่ 4.00 คืออัตราส่วน 4/1 ของคู่ที่ไม่สอดคล้องกันสองชนิด ส่วนคู่ที่สอดคล้องกันไม่มีบทบาท

คู่ทั้ง 20 คู่เดียวกันใน Stata และ R

สคริปต์จำลองพิมพ์จำนวนคู่สี่แบบ ขยายเป็น 40 แถว และประมาณทั้งสองแบบจำลอง logistic regression ทั่วไปให้ 1.83 ช่วงเชื่อมั่น 95% 0.52 ถึง 6.43 conditional logistic regression ให้ 4.00 ช่วงเชื่อมั่น 95% 0.45 ถึง 35.79 โค้ด Stata ยังประมาณแบบจำลองทั่วไปซ้ำโดยจัดกลุ่ม standard error ตามคู่ ซึ่งจะกลับมากล่าวถึงในความเข้าใจผิดที่พบบ่อยใกล้ท้ายบทความ

ค่าประมาณจากการจับคู่ไม่แม่นยำ สัมประสิทธิ์ของมันคือ 1.386 ซึ่งเป็น log ของ 4.00 บนมาตรวัด log นั้น standard error ของมันคือ $\sqrt{1/4 + 1/1} = 1.118$ สูตรนี้ใช้ได้กับคู่ 1:1 ที่มีสิ่งสัมผัสแบบสองค่าหนึ่งตัวและไม่มีตัวแปรร่วมอื่น

standard error ยกกำลังสอง คือความแปรปรวน เท่ากับผลบวก $1/4 + 1/1$ ส่วนใหญ่มาจากคู่เดียวที่เฉพาะ control สัมผัส เมื่อมีคู่ที่ให้ข้อมูล 5 คู่ ช่วงจึงกว้างจากต่ำกว่า 1 ไปไกลเกิน 4.00 ข้อมูลจึงตัดความเป็นไปได้ที่ไม่มีความสัมพันธ์ไม่ได้

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

Stata เรียกผลลัพธ์ของ clogit ว่า conditional (fixed-effects) logistic regression ผลคงที่ (fixed effects) คือ intercept ของคู่ $\omega_k$ ซึ่ง conditional likelihood กำจัดออกแทนการประมาณ ผลลัพธ์จึงไม่มีพจน์คงที่

clogit() ของ R พิมพ์การเรียก coxph() ซึ่งคือการถดถอย Cox ของ R สำหรับข้อมูลเวลาถึงเหตุการณ์ เพราะประมาณ conditional likelihood ผ่านรูทีนนั้น บทความเรื่องการเลือกกลุ่มควบคุมตามเรื่องลิงก์นั้นและความหมายต่อ odds ratio จับคู่

Stata: ป้อนคู่ทั้ง 20 คู่

โค้ด Stata c7_sim.do (บรรทัด 197-225 จาก 265)
clear
input byte anc_case byte anc_ctrl int n_pairs
1 1 8
0 0 7
1 0 4
0 1 1
end
quietly summarize n_pairs
canon pairs.n 0 r(sum)
quietly summarize n_pairs if anc_case == anc_ctrl
canon pairs.concordant 0 r(sum)
quietly summarize n_pairs if anc_case != anc_ctrl
canon pairs.discordant 0 r(sum)
quietly summarize n_pairs if anc_case == 1 & anc_ctrl == 0
canon pairs.case_only_inadequate 0 r(sum)
quietly summarize n_pairs if anc_case == 0 & anc_ctrl == 1
canon pairs.control_only_inadequate 0 r(sum)
expand n_pairs
drop n_pairs
generate int set_id = _n
reshape long anc_, i(set_id) j(role) string
rename anc_ inadequate_anc
generate byte case = role == "case"
bysort set_id: assert _N == 2
bysort set_id (case): assert case[1] == 0 & case[2] == 1
quietly count
canon pairs.rows 0 r(N)
* the pooled 2 x 2 table that ordinary logistic regression sees
tabulate inadequate_anc case
ผลลัพธ์จากการรัน c7_sim.log
. clear

. input byte anc_case byte anc_ctrl int n_pairs

     anc_case  anc_ctrl   n_pairs
  1. 1 1 8
  2. 0 0 7
  3. 1 0 4
  4. 0 1 1
  5. end

. quietly summarize n_pairs

. canon pairs.n 0 r(sum)
CANON c7.pairs.n 20

. quietly summarize n_pairs if anc_case == anc_ctrl

. canon pairs.concordant 0 r(sum)
CANON c7.pairs.concordant 15

. quietly summarize n_pairs if anc_case != anc_ctrl

. canon pairs.discordant 0 r(sum)
CANON c7.pairs.discordant 5

. quietly summarize n_pairs if anc_case == 1 & anc_ctrl == 0

. canon pairs.case_only_inadequate 0 r(sum)
CANON c7.pairs.case_only_inadequate 4

. quietly summarize n_pairs if anc_case == 0 & anc_ctrl == 1

. canon pairs.control_only_inadequate 0 r(sum)
CANON c7.pairs.control_only_inadequate 1

. expand n_pairs
(16 observations created)

. drop n_pairs

. generate int set_id = _n

. reshape long anc_, i(set_id) j(role) string
(j = case ctrl)

Data                               Wide   ->   Long
-----------------------------------------------------------------------------
Number of observations               20   ->   40
Number of variables                   3   ->   3
j variable (2 values)                     ->   role
xij variables:
                      anc_case anc_ctrl   ->   anc_
-----------------------------------------------------------------------------

. rename anc_ inadequate_anc

. generate byte case = role == "case"

. bysort set_id: assert _N == 2

. bysort set_id (case): assert case[1] == 0 & case[2] == 1

. quietly count

. canon pairs.rows 0 r(N)
CANON c7.pairs.rows 40

. * the pooled 2 x 2 table that ordinary logistic regression sees
. tabulate inadequate_anc case

inadequate |         case
      _anc |         0          1 |     Total
-----------+----------------------+----------
         0 |        11          8 |        19
         1 |         9         12 |        21
-----------+----------------------+----------
     Total |        20         20 |        40
ข้อมูลตัวอย่างคำนวณด้วยมือ ป้อนในสคริปต์จำลองของชุดบทความนี้ เป็นข้อความที่ตัดมาจากสคริปต์นั้นและผลลัพธ์ของสคริปต์ บล็อก input พิมพ์รูปแบบของคู่สี่แบบพร้อมจำนวน ขั้น expand และ reshape ทำให้มีหนึ่งแถวต่อสมาชิกของคู่ โดยมี set_id ของแต่ละคู่ และ assert ตรวจว่าแต่ละคู่มี case หนึ่งคนและ control หนึ่งคน คำสั่ง tabulate คือตารางรวมที่ logistic regression ทั่วไปมองเห็น บรรทัดที่ขึ้นต้นด้วย CANON พิมพ์โดยฟังก์ชันช่วยเล็ก ๆ ในสคริปต์ที่สะท้อนผลลัพธ์แต่ละค่าออกมา ให้อ่านตัวเลขท้ายบรรทัดและไม่ต้องสนใจคำนำหน้า

Stata: logistic regression ทั่วไปกับคู่ที่จับคู่

โค้ด Stata c7_sim.do (บรรทัด 240-241 จาก 265)
* ordinary logistic regression ignores the matched sets
logistic case i.inadequate_anc, nolog
ผลลัพธ์จากการรัน c7_sim.log
. * ordinary logistic regression ignores the matched sets
. logistic case i.inadequate_anc, nolog

Logistic regression                                     Number of obs =     40
                                                        LR chi2(1)    =   0.91
                                                        Prob > chi2   = 0.3413
Log likelihood = -27.27303                              Pseudo R2     = 0.0163

----------------------------------------------------------------------------------
            case | Odds ratio   Std. err.      z    P>|z|     [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc |   1.833333   1.174412     0.95   0.344      .522362    6.434448
           _cons |   .7272727   .3379347    -0.69   0.493      .292535    1.808076
----------------------------------------------------------------------------------
Note: _cons estimates baseline odds.
ข้อมูลตัวอย่างคำนวณด้วยมือ ข้อความที่ตัดมาจากสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ เมื่อไม่สนใจชุดที่จับคู่ odds ratio เท่ากับ 1.83 ช่วง 95% คือ 0.52 ถึง 6.43 ในข้อมูล case-control ค่า _cons ขึ้นกับจำนวน control ที่สุ่มต่อ case หนึ่งคน ไม่ใช่ odds พื้นฐานของภาวะแทรกซ้อน

Stata: การประมาณแบบรวมด้วย cluster-robust standard error

โค้ด Stata c7_sim.do (บรรทัด 246-247 จาก 265)
* clustering the standard error on the matched set changes the interval, not the odds ratio
logistic case i.inadequate_anc, vce(cluster set_id) nolog
ผลลัพธ์จากการรัน c7_sim.log
. * clustering the standard error on the matched set changes the interval, not the odds ratio
. logistic case i.inadequate_anc, vce(cluster set_id) nolog

Logistic regression                                     Number of obs =     40
                                                        Wald chi2(1)  =   1.82
                                                        Prob > chi2   = 0.1770
Log pseudolikelihood = -27.27303                        Pseudo R2     = 0.0163

                                    (Std. err. adjusted for 20 clusters in set_id)
----------------------------------------------------------------------------------
                 |               Robust
            case | Odds ratio   std. err.      z    P>|z|     [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc |   1.833333    .823118     1.35   0.177     .7604568    4.419858
           _cons |   .7272727   .1778604    -1.30   0.193     .4503282    1.174534
----------------------------------------------------------------------------------
Note: _cons estimates baseline odds.
ข้อมูลตัวอย่างคำนวณด้วยมือ ข้อความที่ตัดมาจากสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ นี่คือ logistic regression ทั่วไปแบบเดียวกัน แต่ใช้ vce(cluster set_id) เพื่อขอ cluster-robust standard error ซึ่งยอมให้สมาชิกสองคนในคู่เดียวกันสัมพันธ์กันได้ odds ratio ยังเป็น 1.83 มีเพียง standard error และช่วง 95% ที่เปลี่ยน ในที่นี้เป็น 0.76 ถึง 4.42 หัวข้อความเข้าใจผิดที่พบบ่อยใกล้ท้ายบทความอธิบายว่าทำไมสิ่งนี้ไม่ได้ทำให้การประมาณแบบรวมกลายเป็นการวิเคราะห์แบบจับคู่

R: คู่ที่จับคู่และการประมาณ logistic regression ทั่วไป

โค้ด R c7_sim_r.R (บรรทัด 165-198 จาก 222)
pairs <- data.frame(anc_case = c(1, 0, 1, 0), anc_ctrl = c(1, 0, 0, 1), n_pairs = c(8, 7, 4, 1))
canon("pairs.n", sum(pairs$n_pairs), 0)
canon("pairs.concordant", sum(pairs$n_pairs[pairs$anc_case == pairs$anc_ctrl]), 0)
canon("pairs.discordant", sum(pairs$n_pairs[pairs$anc_case != pairs$anc_ctrl]), 0)
canon("pairs.case_only_inadequate", sum(pairs$n_pairs[pairs$anc_case == 1 & pairs$anc_ctrl == 0]), 0)
canon("pairs.control_only_inadequate", sum(pairs$n_pairs[pairs$anc_case == 0 & pairs$anc_ctrl == 1]), 0)
pp <- pairs[rep(seq_len(nrow(pairs)), pairs$n_pairs), ]
pp$set_id <- seq_len(nrow(pp))
m <- rbind(data.frame(set_id = pp$set_id, case = 1L, inadequate_anc = pp$anc_case),
           data.frame(set_id = pp$set_id, case = 0L, inadequate_anc = pp$anc_ctrl))
m <- m[order(m$set_id, m$case), ]
stopifnot(all(table(m$set_id) == 2), all(tapply(m$case, m$set_id, sum) == 1))
canon("pairs.rows", nrow(m), 0)
# the pooled 2 x 2 table that ordinary logistic regression sees
print(table(inadequate_anc = m$inadequate_anc, case = m$case))
a <- sum(m$case == 1 & m$inadequate_anc == 1)
bb <- sum(m$case == 1 & m$inadequate_anc == 0)
cc <- sum(m$case == 0 & m$inadequate_anc == 1)
dd <- sum(m$case == 0 & m$inadequate_anc == 0)
canon("pooled.cases_inadequate", a, 0)
canon("pooled.cases_adequate", bb, 0)
canon("pooled.controls_inadequate", cc, 0)
canon("pooled.controls_adequate", dd, 0)
canon("pooled.odds_cases", a / bb)
canon("pooled.odds_controls", cc / dd)
# ordinary logistic regression ignores the matched sets
fp <- glm(case ~ inadequate_anc, family = binomial, data = m)
print(summary(fp)$coefficients)
bp <- coef(fp)[["inadequate_anc"]]
sp <- sqrt(vcov(fp)["inadequate_anc", "inadequate_anc"])
canon("pooled.or", exp(bp))
canon("pooled.se", sp)
canon("pooled.or_lb", exp(bp - z * sp))
canon("pooled.or_ub", exp(bp + z * sp))
ผลลัพธ์จากการรัน c7_sim_r.log
> pairs <- data.frame(anc_case = c(1, 0, 1, 0), anc_ctrl = c(1,
+     0, 0, 1), n_pairs = c(8, 7, 4, 1))

> canon("pairs.n", sum(pairs$n_pairs), 0)
CANON c7.pairs.n 20

> canon("pairs.concordant", sum(pairs$n_pairs[pairs$anc_case ==
+     pairs$anc_ctrl]), 0)
CANON c7.pairs.concordant 15

> canon("pairs.discordant", sum(pairs$n_pairs[pairs$anc_case !=
+     pairs$anc_ctrl]), 0)
CANON c7.pairs.discordant 5

> canon("pairs.case_only_inadequate", sum(pairs$n_pairs[pairs$anc_case ==
+     1 & pairs$anc_ctrl == 0]), 0)
CANON c7.pairs.case_only_inadequate 4

> canon("pairs.control_only_inadequate", sum(pairs$n_pairs[pairs$anc_case ==
+     0 & pairs$anc_ctrl == 1]), 0)
CANON c7.pairs.control_only_inadequate 1

> pp <- pairs[rep(seq_len(nrow(pairs)), pairs$n_pairs),
+     ]

> pp$set_id <- seq_len(nrow(pp))

> m <- rbind(data.frame(set_id = pp$set_id, case = 1,
+     inadequate_anc = pp$anc_case), data.frame(set_id = pp$set_id,
+     case = 0, inadequate_anc = pp$anc_ctrl))

> m <- m[order(m$set_id, m$case), ]

> stopifnot(all(table(m$set_id) == 2), all(tapply(m$case,
+     m$set_id, sum) == 1))

> canon("pairs.rows", nrow(m), 0)
CANON c7.pairs.rows 40

> print(table(inadequate_anc = m$inadequate_anc, case = m$case))
              case
inadequate_anc  0  1
             0 11  8
             1  9 12

> a <- sum(m$case == 1 & m$inadequate_anc == 1)

> bb <- sum(m$case == 1 & m$inadequate_anc == 0)

> cc <- sum(m$case == 0 & m$inadequate_anc == 1)

> dd <- sum(m$case == 0 & m$inadequate_anc == 0)

> canon("pooled.cases_inadequate", a, 0)
CANON c7.pooled.cases_inadequate 12

> canon("pooled.cases_adequate", bb, 0)
CANON c7.pooled.cases_adequate 8

> canon("pooled.controls_inadequate", cc, 0)
CANON c7.pooled.controls_inadequate 9

> canon("pooled.controls_adequate", dd, 0)
CANON c7.pooled.controls_adequate 11

> canon("pooled.odds_cases", a/bb)
CANON c7.pooled.odds_cases 1.5000

> canon("pooled.odds_controls", cc/dd)
CANON c7.pooled.odds_controls 0.8182

> fp <- glm(case ~ inadequate_anc, family = binomial,
+     data = m)

> print(summary(fp)$coefficients)
                 Estimate Std. Error    z value  Pr(>|z|)
(Intercept)    -0.3184537  0.4646602 -0.6853476 0.4931246
inadequate_anc  0.6061358  0.6405884  0.9462172 0.3440378

> bp <- coef(fp)[["inadequate_anc"]]

> sp <- sqrt(vcov(fp)["inadequate_anc", "inadequate_anc"])

> canon("pooled.or", exp(bp))
CANON c7.pooled.or 1.8333

> canon("pooled.se", sp)
CANON c7.pooled.se 0.6406

> canon("pooled.or_lb", exp(bp - z * sp))
CANON c7.pooled.or_lb 0.5224

> canon("pooled.or_ub", exp(bp + z * sp))
CANON c7.pooled.or_ub 6.4344
ข้อมูลตัวอย่างคำนวณด้วยมือ ข้อความที่ตัดมาจากสคริปต์ R ของสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ data frame m เก็บหนึ่งแถวต่อสมาชิกของคู่ glm() พิมพ์สัมประสิทธิ์ inadequate_anc บนมาตรวัด log odds และบรรทัด CANON ให้ odds ratio ของมัน 1.83 พร้อมขอบเขต บรรทัดที่ขึ้นต้นด้วย CANON พิมพ์โดยฟังก์ชันช่วยเล็ก ๆ ในสคริปต์ที่สะท้อนผลลัพธ์แต่ละค่าออกมา ให้อ่านตัวเลขท้ายบรรทัดและไม่ต้องสนใจคำนำหน้า

Stata: conditional logistic regression กับคู่ที่จับคู่

โค้ด Stata c7_sim.do (บรรทัด 252-253 จาก 265)
* conditional logistic regression compares the case with its own control inside each matched set
clogit case i.inadequate_anc, group(set_id) or nolog
ผลลัพธ์จากการรัน c7_sim.log
. * conditional logistic regression compares the case with its own control inside each matched set
. clogit case i.inadequate_anc, group(set_id) or nolog

Conditional (fixed-effects) logistic regression         Number of obs =     40
                                                        LR chi2(1)    =   1.93
                                                        Prob > chi2   = 0.1650
Log likelihood = -12.89922                              Pseudo R2     = 0.0695

----------------------------------------------------------------------------------
            case | Odds ratio   Std. err.      z    P>|z|     [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc |          4   4.472136     1.24   0.215     .4470826    35.78757
----------------------------------------------------------------------------------
ข้อมูลตัวอย่างคำนวณด้วยมือ ข้อความที่ตัดมาจากสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ option group(set_id) บอก clogit ว่าแถวใดเป็นชุดที่จับคู่แต่ละชุด และ option or พิมพ์ odds ratio Stata แสดงเป็น 4 คือ 4.00 ช่วง 95% คือ 0.45 ถึง 35.79

R: clogit() กับคู่ที่จับคู่

โค้ด R c7_sim_r.R (บรรทัด 206-208 จาก 222)
# conditional logistic regression compares the case with its own control inside each matched set
fc <- clogit(case ~ inadequate_anc + strata(set_id), data = m)
print(summary(fc))
ผลลัพธ์จากการรัน c7_sim_r.log
> fc <- clogit(case ~ inadequate_anc + strata(set_id),
+     data = m)

> print(summary(fc))
Call:
coxph(formula = Surv(rep(1, 40L), case) ~ inadequate_anc + strata(set_id),
    data = m, method = "exact")

  n= 40, number of events= 20

                coef exp(coef) se(coef)    z Pr(>|z|)
inadequate_anc 1.386     4.000    1.118 1.24    0.215

               exp(coef) exp(-coef) lower .95 upper .95
inadequate_anc         4       0.25    0.4471     35.79

Concordance= 0.575  (se = 0.075 )
Likelihood ratio test= 1.93  on 1 df,   p=0.2
Wald test            = 1.54  on 1 df,   p=0.2
Score (logrank) test = 1.8  on 1 df,   p=0.2
ข้อมูลตัวอย่างคำนวณด้วยมือ ข้อความที่ตัดมาจากสคริปต์ R ของสคริปต์จำลองของชุดบทความนี้และผลลัพธ์ของสคริปต์ strata(set_id) ระบุชุดที่จับคู่ การเรียกที่พิมพ์ออกมาเป็นแบบจำลอง Cox เพราะ clogit() ประมาณ conditional likelihood ด้วยรูทีน Cox ของ R คอลัมน์ exp(coef) คือ odds ratio จับคู่ 4.00 และสองคอลัมน์สุดท้ายของตารางที่สองให้ช่วงเชื่อมั่น 95% ของมัน

ทำไม 1.83 และ 4.00 จึงต่างกัน

ตัวเลขสองค่านี้ตอบคำถามต่างกัน odds ratio แบบไม่จับคู่เปรียบเทียบ odds ของการสัมผัสของ case ทั้งหมดกับ control ทั้งหมด odds ratio แบบจับคู่เปรียบเทียบการสัมผัสภายในคู่ ระหว่างทารกที่มีสิ่งที่สมาชิกของคู่มีร่วมกันเหมือนกัน

การรวมทำให้การเปรียบเทียบภายในคู่เลือนไป คู่ที่สอดคล้องกันทุกคู่เพิ่มสิ่งสัมผัสเดียวกันเข้าทั้งสองคอลัมน์ของตารางรวม คู่ที่สัมผัสทั้งสองคนเพิ่มสมาชิกที่สัมผัส 8 คนเข้าแต่ละคอลัมน์ คู่ที่ไม่สัมผัสทั้งสองคนเพิ่มสมาชิกที่ไม่สัมผัส 7 คนเข้าแต่ละคอลัมน์ ในที่นี้ 15 จาก 20 คู่สอดคล้องกัน

เมื่อ case 12 จาก 20 คน และ control 9 จาก 20 คนฝากครรภ์ไม่เพียงพอ ถ้าสิ่งสัมผัสของสมาชิกสองคนในคู่ไม่สัมพันธ์กัน จะมีคู่ที่สอดคล้องกันราว 10 คู่ (0.60 x 0.45 + 0.40 x 0.55 = 0.49 ดังนั้น 0.49 x 20 = 9.8) เทียบกับ 15 คู่ที่พบจริง ผลนี้เข้ากันได้กับการที่ตัวแปรที่ใช้จับคู่สัมพันธ์กับสิ่งสัมผัส แม้คู่ 20 คู่จะแสดงเรื่องนี้ไม่ได้

เมื่อตัวแปรที่ใช้จับคู่ ในที่นี้คืออายุเป็นวัน สัมพันธ์กับสิ่งสัมผัส การจับคู่จะดึงสิ่งสัมผัสของ control เข้าหาของ case การวิเคราะห์ที่ไม่สนใจชุดจึงมักทำให้ odds ratio ลำเอียงเข้าหา 1 [5] ส่วนที่เพิ่มร่วมกันข้างต้นคือวิธีที่ความลำเอียงนั้นปรากฏในตารางรวม

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

non-collapsibility ซึ่งนิยามไว้ด้านล่างสำหรับแบบจำลอง cohort คือช่องว่างระหว่าง odds ratio ภายในกลุ่มย่อยกับ odds ratio ที่เฉลี่ยข้ามกลุ่มย่อย ซึ่งเกิดขึ้นเมื่อกลุ่มย่อยทำนายผลลัพธ์ได้แม้ไม่มีตัวกวน ในที่นี้ทุกคู่มี case หนึ่งคนและ control หนึ่งคนที่อายุเป็นวันเท่ากัน ตัวคู่เองและอายุจึงทำนายไม่ได้ว่าใครเป็น case ช่องว่างระหว่าง 1.83 กับ 4.00 จึงไม่ได้มาจาก non-collapsibility ช่องว่างนี้สะท้อนแบบแผนการศึกษา คือการเปรียบเทียบภายในคู่กับการเปรียบเทียบแบบรวม ซึ่งการจับคู่อาจดึงสิ่งสัมผัสของ control เข้าหาของ case ร่วมกับความบังเอิญเมื่อคู่ที่ให้ข้อมูลมีน้อยเช่นนี้

แล้วตัวเลขไหนควรอยู่ในรายงาน การศึกษาออกแบบมาเพื่อเปรียบเทียบภายในคู่ odds ratio จับคู่จึงมักเป็นค่าที่ควรรายงาน ควรมาพร้อมช่วงเชื่อมั่นของมัน 0.45 ถึง 35.79 ซึ่งเข้ากันได้ทั้งกับการไม่มีความสัมพันธ์และกับความสัมพันธ์ที่แรง

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

ตัวอย่างคำนวณด้วยมือ เปลี่ยนจำนวนคู่ในแต่ละรูปแบบ คือสัมผัสทั้งสองคน เฉพาะ case สัมผัส เฉพาะ control สัมผัส และไม่สัมผัสทั้งสองคน แผงจะสร้างตารางรวม odds ratio แบบไม่จับคู่ odds ratio แบบจับคู่ และจำนวนคู่ที่ให้ข้อมูล (คู่ที่ไม่สอดคล้องกัน) ใหม่ แผงเริ่มที่ตัวอย่างคำนวณด้วยมือ คู่ 8, 4, 1 และ 7 คู่ให้ 1.83 และ 4.00 จากคู่ที่ให้ข้อมูล 5 คู่ ลองเพิ่มคู่ที่สอดคล้องกันแล้วดู odds ratio แบบจับคู่คงที่ ขณะที่แบบไม่จับคู่เปลี่ยน

link เดียว แบบจำลองสามแบบ: ตารางสำหรับตัดสินใจ

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

แถวสุดท้ายใช้ อัตราส่วนฮาซาร์ด (hazard ratio) ซึ่งเปรียบเทียบอัตราการเกิดเหตุการณ์ ณ ขณะใดขณะหนึ่งในผู้ที่ยังไม่เกิดเหตุการณ์ ภายใต้ ข้อสมมติฮาซาร์ดเป็นสัดส่วน (proportional hazards) อัตราส่วนนั้นคงที่ตลอดช่วงติดตาม

แบบจำลองทั้งสามในบทความนี้ คำสั่งเป็นแม่แบบ y คือผลลัพธ์ (case ในแถวที่จับคู่) x คือสิ่งสัมผัส และ set_id คือรหัสของชุดที่จับคู่
ชนิดของผลลัพธ์และแบบแผนการสุ่มตัวอย่างแบบจำลองStataRสัมประสิทธิ์ที่ยกกำลัง
ผลลัพธ์สองค่า คนอิสระต่อกัน เช่น cohort หรือการศึกษา case-control ที่ไม่จับคู่Logistic regressionlogistic y i.xglm(y ~ x, family = binomial)odds ratio ของผู้สัมผัสเทียบกับผู้ไม่สัมผัส โดยคงตัวแปรร่วมอื่นไว้ ในการศึกษา case-control สิ่งที่ค่านี้ประมาณขึ้นกับวิธีสุ่มกลุ่มควบคุม (เป็น odds ratio เฉพาะเมื่อกลุ่มควบคุมมาจากผู้ที่ยังไม่เกิดผลลัพธ์เมื่อสิ้นสุดการติดตาม)
หมวดที่มีลำดับ คนอิสระต่อกันOrdinal logistic regression (proportional odds)ologit y i.xMASS::polr(factor(y, ordered = TRUE) ~ x)odds ratio ค่าเดียวสำหรับหมวดที่สูงขึ้น สมมติว่าเท่ากันทุกจุดตัด ตรวจด้วย Brant test
ผลลัพธ์สองค่า (case หรือ control) ชุดที่จับคู่ รวมถึง control ที่สุ่มจากกลุ่มเสี่ยงConditional logistic regressionclogit case i.x, group(set_id) orsurvival::clogit(case ~ x + strata(set_id))odds ratio ภายในชุดที่จับคู่ เมื่อสุ่มจากกลุ่มเสี่ยงจะประมาณ hazard ratio ภายใต้ proportional hazards ตัวแปรที่ใช้จับคู่ไม่มีสัมประสิทธิ์

สิ่งที่ odds ratio ทั้งสามมีร่วมกัน

ในแบบจำลองทั้งสาม สัมประสิทธิ์คือผลต่างของ log odds และเลขชี้กำลังของมันคือ odds ratio สิ่งที่ odds ratio นั้นมีเงื่อนไขอยู่ต่างกัน logistic regression และ ordinal logistic regression มีเงื่อนไขตามตัวแปรร่วมอื่น ส่วน conditional logistic regression มีเงื่อนไขตามชุดที่จับคู่ด้วย ชุดเงื่อนไขแต่ละชุดอาจนิยาม ปริมาณเป้าหมายของการประมาณ (estimand) ของตัวเอง คือปริมาณที่แน่ชัดซึ่งการศึกษาตั้งใจประมาณ

ใน birth cohort odds ratio ของ ANC ที่ 2.00 และ 2.01 มีเงื่อนไข (conditional) ตามการคลอดก่อนกำหนด การคลอดก่อนกำหนดเป็น ปัจจัยพยากรณ์โรค (prognostic factor) คือลักษณะที่ทำนายผลลัพธ์ แบบจำลองที่ไม่มีมันจะให้ odds ratio ระดับประชากร (marginal) ซึ่งเฉลี่ยตามสถานะการคลอดก่อนกำหนด odds ratio ไม่ยุบรวมได้ (non-collapsible) ดังนั้นทั้งสองค่าไม่จำเป็นต้องเท่ากัน แม้การคลอดก่อนกำหนดจะไม่สัมพันธ์กับ ANC [6]

odds ratio แบบมีเงื่อนไขกับแบบระดับประชากร รวมทั้ง hazard ratio ทั้งสองแบบ อาจต่างกันได้แม้ไม่มีตัวกวนเลย ในการทดลองแบบสุ่ม การปรับด้วยปัจจัยพยากรณ์โรคที่มีผลแรงมักทำให้ odds ratio ขยับห่างจาก 1 มากขึ้น นั่นคือการเปลี่ยน estimand ไม่ใช่การกำจัดความลำเอียง กรณีการทดลองแบบสุ่ม [7] คำนวณไว้ในตอนของชุดบทความนี้ว่าด้วย odds ratio แบบมีเงื่อนไขและระดับประชากร

รายงานจึงควรระบุว่า odds ratio แต่ละค่ามีเงื่อนไขตามอะไร ควบคู่กับแบบจำลอง การฝากครรภ์ไม่เพียงพออาจนำไปสู่การคลอดก่อนกำหนด ทำให้การคลอดก่อนกำหนดเป็น ตัวแปรสื่อกลาง (mediator) คือขั้นบนเส้นทางไปสู่ภาวะแทรกซ้อน ถ้าเป็นเช่นนั้น odds ratio 2.00 และ 2.01 ละเส้นทางนั้นไป ต่างจาก ผลรวม (total effect) ของ ANC ที่รวมมันไว้ ตัวแปรร่วมใดควรอยู่ในเงื่อนไขเป็นคำถามเชิงเหตุผล ซึ่งควรตัดสินด้วยแผนภาพเชิงเหตุผลก่อนการวิเคราะห์

ความเข้าใจผิดที่พบบ่อยและวิธีแก้

  • "clogit ให้ odds ratio ที่ใหญ่กว่า จึงแม่นยำกว่า"

    ขนาดของค่าประมาณไม่ได้บอกอะไรเกี่ยวกับความแม่นยำ ในที่นี้ 4.00 จากการจับคู่อยู่ห่างจาก 1 มากกว่า 1.83 ที่ไม่จับคู่ แต่ช่วงของมัน 0.45 ถึง 35.79 กว้างกว่า 0.52 ถึง 6.43 มาก ในข้อมูลอื่น odds ratio จับคู่อาจเป็นค่าที่เล็กกว่า

    วิธีแก้: ตัวเลขสองค่านี้ตอบคำถามต่างกัน: odds ratio แบบจับคู่เปรียบเทียบการได้รับสิ่งสัมผัสภายในคู่ ส่วน odds ratio แบบไม่คำนึงถึงการจับคู่รวมทุกคู่เข้าด้วยกัน ควรเลือกแบบจำลองจากรูปแบบการศึกษาและ estimand แล้วรายงานช่วงความเชื่อมั่น

  • "ologit ใช้ได้กับผลลัพธ์แบบลำดับทุกชนิด"

    ใน cohort จำลอง แบบจำลองแบบลำดับให้การคลอดก่อนกำหนดมี odds ratio ค่าเดียวคือ 1.84 ขณะที่ odds ratio ตามจุดตัดอยู่ที่ 1.40 ถึง 3.31 Brant test ชี้การคลอดก่อนกำหนด (P = 0.003) และไม่ชี้ ANC (P = 0.999)

    วิธีแก้: ologit สมมติ odds ratio ค่าเดียวที่ทุกจุดตัด ให้ตรวจข้อสมมตินั้น และบอกว่าทำอะไรเมื่อข้อสมมตินั้นไม่เป็นจริง

  • "clogit ทิ้งคู่ที่สอดคล้องกัน"

    ทุกคู่เข้าสู่ conditional likelihood เมื่อมีสิ่งสัมผัสแบบสองค่าหนึ่งตัวและไม่มีตัวแปรร่วมอื่น คู่ที่สอดคล้องกันให้ความน่าจะเป็น 0.5 ไม่ว่าสัมประสิทธิ์เป็นเท่าใด จึงไม่มีข้อมูลเกี่ยวกับสัมประสิทธิ์นั้น เมื่อมีตัวแปรร่วมมากขึ้น มันยังให้ข้อมูลแก่ตัวแปรที่ต่างกันภายในคู่ได้

    วิธีแก้: เก็บทุกชุดที่จับคู่ไว้ในข้อมูล และรายงานว่ามีกี่ชุดที่ไม่สอดคล้องกันในสิ่งสัมผัส เพราะชุดเหล่านั้นกำหนดความแม่นยำของ odds ratio ของสิ่งสัมผัส

  • "cluster-robust standard error แก้การวิเคราะห์ที่ไม่สนใจการจับคู่ได้"

    cluster-robust standard error ยอมให้หน่วยสังเกตภายในกลุ่ม ในที่นี้คือคู่ มีความสัมพันธ์กันได้ การจัดกลุ่มตามคู่ใน logistic regression แบบรวม ดังผลลัพธ์ Stata ที่แสดงพร้อมโค้ดของคู่ที่จับคู่ ทำให้ช่วงเชื่อมั่น 95% แคบลงเป็น 0.76 ถึง 4.42 จาก 0.52 ถึง 6.43 ช่วงแคบลงในที่นี้เพราะแต่ละคู่มี case เพียงหนึ่งคนโดยการออกแบบ odds ratio ยังเป็น 1.83 ซึ่งยังเป็นการเปรียบเทียบแบบรวม ช่วงทั้งสองจึงไม่ตอบคำถามภายในคู่

    วิธีแก้: วิเคราะห์ชุดที่จับคู่ด้วย conditional logistic regression standard error ที่ต่างไปเปลี่ยนช่วง ไม่เปลี่ยนค่าประมาณหรือคำถามที่มันตอบ

  • "logistic regression ทั่วไปที่ใส่ตัวบ่งชี้หนึ่งตัวต่อคู่ให้ odds ratio จับคู่"

    ตัวบ่งชี้ของแต่ละคู่ถูกประมาณจากทารกเพียงสองคน วิธี maximum likelihood ทั่วไป ซึ่งเป็นวิธีประมาณค่าที่ใช้กันบ่อย จึงผลัก odds ratio ให้ห่างจาก 1 ยิ่งขึ้น สำหรับคู่ 1:1 ที่มีสิ่งสัมผัสแบบสองค่าหนึ่งตัว odds ratio ของแบบจำลองนั้นคือกำลังสองของค่าจับคู่ [5]

    วิธีแก้: ใช้ clogit วางเงื่อนไขตามชุดที่จับคู่ conditional likelihood กำจัด intercept ของคู่แทนการประมาณมัน

  • "odds ratio 4.00 หมายถึงความเสี่ยงของภาวะแทรกซ้อนเป็นสี่เท่า"

    odds ratio จับคู่คืออัตราส่วนของ odds ภายในคู่ ไม่ใช่อัตราส่วนของความเสี่ยง ในที่นี้ control แต่ละคนมาจากทารกที่ยังอยู่ในกลุ่มเสี่ยง ณ เวลาที่ case เกิดเหตุการณ์ ภายใต้ proportional hazards แบบแผนนี้ทำให้ค่านั้นประมาณ hazard ratio ไม่ใช่ risk ratio ที่ 28 วัน บทความเรื่องการเลือกกลุ่มควบคุมคำนวณเรื่องนี้ไว้

    วิธีแก้: ระบุมาตรวัดที่แบบแผนการศึกษาประมาณ และรายงาน 4.00 พร้อมช่วงเชื่อมั่น 95% ของมัน 0.45 ถึง 35.79

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

อภิธานศัพท์

logit link
link ที่เป็น ln{p/(1 - p)} หรือ log odds ซึ่งผูกตัวทำนายเชิงเส้นเข้ากับความน่าจะเป็น สัมประสิทธิ์ของมันคือ log odds ratio
cut point (จุดตัดระหว่างหมวดของผลลัพธ์แบบลำดับ)
เส้นแบ่งระหว่างหมวดสองหมวดที่อยู่ติดกันของผลลัพธ์แบบลำดับ จุดตัดแต่ละจุดมี intercept ของตัวเองใน ordinal logistic regression
cumulative logit
log odds ของการอยู่ที่หมวดหนึ่งหรือสูงกว่าของผลลัพธ์แบบลำดับ มีหนึ่งค่าต่อจุดตัดแต่ละจุด
ordinal logistic regression
แบบจำลองสำหรับผลลัพธ์แบบลำดับที่ประมาณ cumulative logit ที่จุดตัดทุกจุดโดยใช้ชุดความชันร่วมกันชุดเดียว
proportional odds assumption (ข้อสมมติ proportional odds)
ข้อสมมติที่ว่า odds ratio ของตัวแปรร่วมเท่ากันที่ทุกจุดตัดของผลลัพธ์แบบลำดับ
Brant test
Wald test ของข้อสมมติ proportional odds ที่เปรียบเทียบความชันของแบบจำลองแยกตามจุดตัด ทั้งโดยรวมและทีละตัวแปรร่วม
matched set (ชุดที่จับคู่)
case หนึ่งคนพร้อม control ที่สุ่มหรือจับคู่ให้มัน คือหนึ่งคู่ในการจับคู่ 1:1
concordant pair (คู่ที่สอดคล้องกัน)
คู่ที่จับคู่ซึ่ง case และ control มีสิ่งสัมผัสเหมือนกัน เมื่อมีสิ่งสัมผัสแบบสองค่าหนึ่งตัว คู่นี้ไม่ให้ข้อมูลเกี่ยวกับสิ่งสัมผัสใน conditional logistic regression
discordant pair (คู่ที่ไม่สอดคล้องกัน)
คู่ที่จับคู่ซึ่ง case และ control ต่างกันในสิ่งสัมผัส เมื่อจับคู่ 1:1 มีสิ่งสัมผัสแบบสองค่าหนึ่งตัวและไม่มีตัวแปรร่วมอื่น odds ratio จับคู่คืออัตราส่วนของคู่สองชนิด
conditional likelihood
ความน่าจะเป็นที่สมาชิกคนใดของชุดที่จับคู่เป็น case เมื่อกำหนดว่าชุดมี case เพียงหนึ่งคน intercept ของชุดหักล้างกันหมดไป
conditional logistic regression
แบบจำลองโลจิสติกที่ให้ intercept แยกกันสำหรับชุดที่จับคู่แต่ละชุด ประมาณผ่าน conditional likelihood เพื่อให้ case แต่ละคนถูกเปรียบเทียบกับ control ของมันเท่านั้น
risk-set sampling (การสุ่มกลุ่มควบคุมจากกลุ่มเสี่ยง)
การเลือก control ของ case แต่ละคนจากผู้ที่ยังอยู่ในกลุ่มเสี่ยง ณ เวลาที่ case เกิดเหตุการณ์
non-collapsibility (ความไม่ยุบรวมได้)
คุณสมบัติของ odds ratio และ hazard ratio ที่ค่าระดับประชากรอาจต่างจากค่ามีเงื่อนไขร่วมโดยไม่มีตัวกวนเลย
mediator (ตัวแปรสื่อกลาง)
ตัวแปรบนเส้นทางเหตุและผลจากสิ่งสัมผัสไปสู่ผลลัพธ์ ซึ่งผลส่วนหนึ่งอาจผ่านไป
cluster-robust standard error
standard error ที่ยอมให้หน่วยสังเกตภายในกลุ่ม เช่นคู่ที่จับคู่ สัมพันธ์กันได้ เปลี่ยนช่วงแต่ไม่เปลี่ยนค่าประมาณ

เอกสารอ้างอิง

  1. StataCorp. Stata Base Reference Manual: [R] logit. Stata documentation (PDF) at stata.com. Accessed 11 October 2026. https://www.stata.com/manuals/rlogit.pdf
  2. StataCorp. Stata Base Reference Manual: [R] ologit. Stata documentation (PDF) at stata.com. Accessed 11 October 2026. https://www.stata.com/manuals/rologit.pdf
  3. Brant R. Assessing proportionality in the proportional odds model for ordinal logistic regression. Biometrics. 1990;46(4):1171 ff. doi:10.2307/2532457 https://doi.org/10.2307/2532457
  4. StataCorp. Stata Base Reference Manual: [R] clogit. Stata documentation (PDF) at stata.com. Accessed 11 October 2026. https://www.stata.com/manuals/rclogit.pdf
  5. Breslow NE, Day NE. Statistical Methods in Cancer Research, Volume I: The Analysis of Case-Control Studies. IARC Scientific Publications No. 32. Lyon: International Agency for Research on Cancer; 1980. https://publications.iarc.who.int/Book-And-Report-Series/Iarc-Scientific-Publications/Statistical-Methods-In-Cancer-Research-Volume-I-The-Analysis-Of-Case-Control-Studies-1980
  6. Greenland S, Pearl J, Robins JM. Confounding and collapsibility in causal inference. Stat Sci. 1999;14(1):29-46. doi:10.1214/ss/1009211805 https://doi.org/10.1214/ss/1009211805
  7. Gail MH, Wieand S, Piantadosi S. Biased estimates of treatment effect in randomized experiments with nonlinear regressions and omitted covariates. Biometrika. 1984;71(3):431-444. doi:10.1093/biomet/71.3.431 https://doi.org/10.1093/biomet/71.3.431

ประเด็นสำคัญ

  • logistic regression, ordinal logistic regression และ conditional logistic regression ใช้ logit link ร่วมกัน สัมประสิทธิ์ที่ยกกำลังแต่ละค่าจึงเป็น odds ratio
  • ชนิดของผลลัพธ์และแบบแผนการสุ่มตัวอย่าง ที่ระบุก่อนพิมพ์คำสั่งใด ๆ เป็นตัวกำหนดว่าจะใช้แบบจำลองใดในสามแบบ
  • ordinal logistic regression สมมติ odds ratio ค่าเดียวที่ทุกจุดตัด ให้ตรวจทั้งโดยรวมและรายตัวแปรร่วมด้วย Brant test
  • ในคู่ที่จับคู่ 1:1 ที่มีสิ่งสัมผัสแบบสองค่าหนึ่งตัว มีเพียงคู่ที่ไม่สอดคล้องกันที่ให้ข้อมูลแก่ odds ratio ของสิ่งสัมผัส 4 ต่อ 1 ให้ 4.00 ขณะที่การรวมคู่ให้ 1.83
  • ค่าประมาณจากการจับคู่ตอบคำถามของแบบแผนการศึกษา แต่เมื่อมีคู่ที่ไม่สอดคล้องกัน 5 คู่ ช่วง 0.45 ถึง 35.79 แสดงว่าไม่แม่นยำ

อ่านต่อในวิกิ: [[conditional-vs-marginal-or-th]] [[case-control-sampling-schemes-stata-and-r-th]] [[ordinal-logistic-regression]] [[conditional-logistic-regression]] [[multinomial-vs-ordinal-regression]]

0
ถึงนักอ่านชาวไทยและต่างชาติทำความเข้าใจบริบททางการแพทย์ของผมอ่านต่อ →ถึงนักอ่านชาวไทยและต่างชาติทำความเข้าใจเนื้อหาของผมที่นอกเหนือจากการแพทย์อ่านต่อ →

ความคิดเห็น

ยังไม่มีความคิดเห็น มาเป็นคนแรกกันเลย

เข้าสู่ระบบเพื่อแสดงความคิดเห็น