ทางเลือกแบบ Population-Averaged: GLS, GEE และจุดที่โชคของ linear model หมดลง

บทคัดย่อ

Marginal model ไม่อธิบายว่า correlation ภายในคนมาจากไหน มันระบุ covariance ออกมาตรง ๆ แล้วเดินหน้าต่อ บทความนี้เอาการปฏิเสธนั้นมาคิดจริงจัง แยกคำถามคลินิกสองข้อหลัง output ชุดเดียวกัน คือระดับบุคคลกับระดับประชากร แล้วไล่ผ่าน GLS และ GEE ซึ่งประกอบด้วย mean model, working correlation structure และ sandwich variance estimator โดยค่าประมาณจุดยังถูกต้องแม้ working correlation จะผิด เมื่อใช้ identity link ทั้งสองตระกูลให้คำตอบตรงกันพอดี จึงเป็นเหตุผลที่ outcome ต่อเนื่องยกโทษให้ความต่างนี้ แต่เมื่อใช้ logit link ไม่ตรงกัน โมเดลจำลองตัวเดียวกันให้ conditional odds ratio 6.05 และ population-averaged odds ratio 3.72 จึงควรเลือกด้วย estimand และสมมติฐานเรื่องข้อมูลขาดหาย ไม่ใช่ความเคยชิน

Infographic summary of population-averaged approaches: GLS and GEE — ภาพสรุป

เรามาถึงไหนแล้ว. ซีรีส์นี้ว่าด้วยการศึกษาที่ผู้เข้าร่วมหนึ่งคนให้ข้อมูลหลายแถว — ในที่นี้คือผิวสี่ตำแหน่ง ในงานอื่นอาจเป็นการติดตามหลายครั้ง — ทำให้แถวเหล่านั้นเป็นอิสระต่อกันไม่ได้ วิธีแก้ที่ใช้กันคือ mixed model ซึ่งให้ค่าออฟเซ็ตประจำตัวคนละหนึ่งค่า ใส่ซ้ำเหมือนกันในทุกแถวของคนคนนั้น และค่าตัวเดียวนี้เองที่เป็นตัวสร้าง correlation ระหว่างค่าวัดของคนคนเดียวกัน ตอนที่ 6 จบลงด้วยโมเดลที่ปฏิเสธจะทำแบบนั้นทั้งหมด: gls ระบุ marginal covariance matrix ออกมาตรง ๆ และไม่เคยถามว่า correlation มาจากไหน ตอนนี้เราจะเอาการ "ปฏิเสธไม่อธิบาย" นั้นมาคิดอย่างจริงจัง เพราะมันคือทางเลือกที่สามที่ใช้ได้จริง ไม่ใช่ทางเลือกที่ด้อยกว่า สำหรับ outcome แบบต่อเนื่อง การปฏิเสธนี้แทบไม่ต้องจ่ายอะไรเลย แต่สำหรับ outcome แบบ binary มันเปลี่ยนค่าสัมประสิทธิ์เลยทีเดียว

เนื้อหาข้างล่างไม่ได้ตั้งสมมติว่าคุณอ่านตอนที่ 1 ถึง 6 มาแล้ว และคำศัพท์แต่ละคำจะถูกแนะนำใหม่ตรงจุดที่ใช้ครั้งแรก ถ้าอยากเริ่มจากภาพรวมที่ค่อย ๆ ไล่ระดับก่อน แผนที่ของทั้งซีรีส์อยู่ที่นี่: คู่มือซีรีส์

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

ที่ราวด์วอร์ด

แพทย์สองคนยืนดูผลลัพธ์ชุดเดียวกันจากข้อมูลจำลองชุดเดิม คนแรกถามถึงผู้ป่วยที่อยู่ตรงหน้า คือผู้หญิงที่เป็น moderate atopic dermatitis ซึ่งจากค่าที่เคยวัดไว้ก่อนหน้านี้ เรารู้อยู่แล้วว่าผิวของเธอเป็นคนที่ barrier รั่วมากกว่าคนทั่วไป เธออยากรู้ว่าค่าที่วัดจากผิวของผู้ป่วยคนนี้แห้งกว่าคนที่มี barrier รั่วพอ ๆ กันแต่ไม่ได้เป็นโรคนี้อยู่เท่าไร แพทย์คนที่สองถามคำถามที่กว้างกว่า: ถ้าเอาประชากร atopic dermatitis ทั้งหมดในโครงการไปเทียบกับประชากรที่ไม่มีโรคภูมิแพ้เลย โดยเฉลี่ยแล้วผิวจะแห้งกว่ากันเท่าไร ทั้งสองเป็นคำถามทางคลินิกที่ชอบธรรมทั้งคู่ และไม่มีคำถามไหนเป็น counterfactual ระดับบุคคล คำถามแรกเทียบคนสองคนที่ถูกตรึงไว้ที่จุดตั้งต้นระดับบุคคลเดียวกัน ไม่ได้เทียบผู้ป่วยคนนี้กับตัวเธอเองในเวอร์ชันที่ไม่เคยเป็นโรค ทั้งสองคำถามตอบได้จากข้อมูลจำลองชุดเดียวกัน

แต่คำตอบไม่ได้เท่ากันเสมอไป

Outcome (จำลอง)	Conditional (subject-specific)	Population-averaged (marginal)
TEWL ที่ hand, AD − control (g·m⁻²·h⁻¹)	+8.0	+8.0
Xerosis ที่ตำแหน่งหนึ่ง, AD เทียบ control (odds ratio)	6.05	3.72

ทุกค่าในตารางนี้เป็นคุณสมบัติที่แน่นอนของโมเดลจำลองที่สร้างมันขึ้นมา รวมถึง odds ratio ทั้งสองค่า สูตรประมาณที่คุ้นเคยกัน ซึ่งจะเดินให้ดูทีละขั้นในช่วงหลังของบทความนี้ ให้ค่า population-averaged odds ratio ประมาณ 4.0 ส่วนการอินทิเกรตโมเดลเดียวกันนั้นตรง ๆ ให้ค่าที่แน่นอนคือ 3.72 และส่วนต่างระหว่างสองค่านี้เป็นของสูตรประมาณ ไม่ใช่ของโมเดล

ไม่มีบรรทัดไหนในตารางนี้ที่ผิด และไม่มีบรรทัดไหนที่เป็น "ค่าที่แก้ให้ถูกแล้ว" ของอีกบรรทัดหนึ่ง บรรทัดแรกคือความบังเอิญที่ทำให้ linear mixed model ดูใจดีกับเรา ส่วนบรรทัดที่สองคือจุดที่ความบังเอิญนั้นหมดอายุ

ภาพประกอบ: สองปรัชญาวางเทียบกัน คือเงาผู้ป่วยคนเดียวกับภาพฝูงชน โดยแต่ละฝั่งมีคำถามของตัวเอง

สองปรัชญา หนึ่งตาราง

ก่อนลุยของยาก: ตารางข้างล่างวางวิธีมองข้อมูลชุดเดียวกันสองแบบไว้ข้างกัน แบบแรกสร้าง correlation ภายในคนขึ้นมาจากสิ่งที่เราเรียกชื่อและประมาณค่าได้ คือ person effect หนึ่งค่าต่อผู้ป่วยหนึ่งคน ส่วนแบบที่สองแค่เขียนลงไปตรง ๆ ว่าค่าวัดของคนคนหนึ่งเกาะกันแค่ไหน แล้วเดินหน้าต่อ พูดง่าย ๆ คือความต่างระหว่าง "อธิบายว่าทำไมค่าสองค่าจากผู้ป่วยคนเดียวกันถึงคล้ายกัน" กับ "แค่เผื่อไว้ว่ามันคล้ายกัน" มีเพียงแบบแรกที่ให้ค่าความแตกต่างระหว่างบุคคลออกมาให้รายงาน และมีเพียงแบบที่สองที่ไม่สนใจว่าความแตกต่างนั้นมาจากไหน

โมเดลทุกตัวในซีรีส์นี้จนถึงตอนนี้ สร้าง marginal covariance ขึ้นมาจากชิ้นส่วนที่มองไม่เห็น: หนึ่งค่าต่อหนึ่งคน, variance ของค่านั้น, และ residual variance ตอนที่ 4 เขียนสูตรนี้ไว้ครั้งเดียวจบเป็น $V_i = Z_i G Z_i' + R_i$ ส่วน marginal model ปฏิเสธที่จะให้สูตร มันเขียน $V_i$ ลงไปเลย

	Random-effects (conditional)	Marginal (population-averaged)
สิ่งที่คุณระบุ	$Z_i$, $G$, $R_i$	$V_i$ โดยตรง หรือเวอร์ชัน working ของมัน
correlation มาจากไหน	จากค่าแฝงที่ถูกแชร์ $b_i$	ไม่มีที่มา — ประกาศออกมาเฉย ๆ
$\beta$ หมายถึงอะไร	การเปลี่ยนแปลงสำหรับ คนคนหนึ่ง โดยตรึง $b_i$ ไว้	การเปลี่ยนแปลงของ ค่าเฉลี่ยประชากร
ผลลัพธ์ระดับบุคคล	$\widehat{\tau}^{2}$, ICC, BLUPs, subject-specific prediction	ไม่มี
วิธีประมาณค่า	likelihood (REML/ML)	least squares หรือ estimating equations

สองคอลัมน์นี้ไม่ได้แย่งงานกันทำ มันคือสองคำตอบต่อคำถามว่า correlation มีไว้เพื่ออะไร — เป็นปริมาณทางวิทยาศาสตร์ที่เราอยากประมาณค่า หรือเป็น nuisance ที่ต้องรองรับให้ถูกต้องเพื่อให้ mean structure ถูกประมาณอย่างซื่อสัตย์

`gls`: พูด covariance ออกมาดัง ๆ

ก่อนลุยของยาก: โค้ดข้างล่างทำอยู่อย่างเดียว คือบอกโปรแกรมว่าตารางขนาดสี่คูณสี่ของ variance และ correlation ระหว่างผิวสี่ตำแหน่งของคนคนหนึ่งหน้าตาเป็นอย่างไรได้บ้าง โดยไม่พูดถึง person effect เลยสักคำ ทุกอย่างที่ความไม่เป็นอิสระต้องการถูกป้อนด้วยมือผ่านอาร์กิวเมนต์สองตัวเท่านั้น

Marginal specification ที่ตรงไปตรงมาที่สุดคือตัวที่ตอนที่ 6 แนะนำไปแล้ว และไม่ต้องใช้เครื่องมือใหม่เลย

# unstructured marginal covariance, no latent person effect (Posts 6, 7)
fitU <- gls(tewl ~ site * group + age + sex + phototype,
            correlation = corSymm(form = ~ siteN | id),
            weights     = varIdent(form = ~ 1 | site), data = skin)

อ่านสองอาร์กิวเมนต์นี้เป็นประโยคเดียว: corSymm ให้ค่า correlation อิสระหกค่าระหว่างสี่ตำแหน่ง และ varIdent ให้ variance เฉพาะของแต่ละตำแหน่งอีกสี่ค่า รวมกันคือใช้ครบพอดีสิบค่าอิสระของ symmetric covariance matrix ขนาด $4\times4$ ซึ่งเป็นเพดานที่เรานับไว้ในตอนที่ 6 สิ่งที่ fit นี้ประมาณค่าคือเวกเตอร์ fixed effect $\beta$ กับปริมาณ covariance สิบค่า ส่วนสิ่งที่มันไม่ได้ประมาณคือ $\tau^2$ เพราะในโมเดลไม่มี $b_i$ ที่จะเป็นเจ้าของค่านั้น ความไม่เป็นอิสระภายในคนถูกจัดการครบถ้วน เพียงแต่จัดการด้วยการประกาศ ไม่ใช่ด้วยการแยกส่วน

GEE: mean model, working correlation และแซนด์วิช

ก่อนลุยของยาก: GEE ง่ายกว่าสัญลักษณ์ของมันมาก ขั้นแรกคุณเขียนว่าค่าเฉลี่ยของ outcome ควรหน้าตาเป็นอย่างไร ซึ่งเป็นส่วนที่คุณสนใจจริง ๆ ขั้นที่สองคุณเดาว่าค่าวัดซ้ำ ๆ ของผู้ป่วยคนหนึ่งเกาะกันอย่างไร ขั้นที่สามคุณปล่อยให้การกระจายจริงของข้อมูลมาแก้ standard error ให้ เผื่อว่าที่เดาไว้นั้นไม่ดี สมการข้างล่างเป็นเพียงวิธีพูดสามข้อนี้ให้เป็นทางการเท่านั้น ถ้าพีชคณิตเริ่มเบลอ ให้จำประโยคเดียวไว้ก่อน คือการเดา correlation ผิดทำให้เสียความแม่นยำ ไม่ได้ทำให้ผลผิด ตราบใดที่ค่าเฉลี่ยถูกระบุไว้ถูกต้อง

Generalized estimating equations ขยายแนวคิดเดียวกันนี้ไปยัง outcome ที่ไม่ต่อเนื่อง GEE มีองค์ประกอบสามอย่างพอดี และการแยกพูดทีละอย่างก็คือทั้งหมดของวิธีนี้แล้ว

หนึ่ง mean model: link function กับ linear predictor เหมือน generalized linear model ธรรมดาทุกประการ สอง working correlation structure ซึ่งอธิบายว่าค่าที่วัดภายใน cluster เดียวกันเกาะกันอย่างไร สาม sandwich (robust) variance estimator ซึ่งเป็นตัวที่ทำให้องค์ประกอบที่สองมีความยืดหยุ่น ค่าประมาณของพารามิเตอร์ได้มาจากการแก้สมการ

\[ \sum_{i=1}^{N} D_i' V_i^{-1}\left(Y_i - \mu_i\right) = 0, \qquad D_i = \frac{\partial\mu_i}{\partial\beta'} \]

ซึ่งไม่ใช่ likelihood แต่เป็นชุดสมการที่คำตอบของมัน consistent สำหรับ $\beta$ ตราบใดที่ mean model ถูกระบุถูกต้อง ลองอ่านทีละชิ้นจากขวาไปซ้าย: $Y_i - \mu_i$ คือระยะที่ cluster นี้เบี่ยงจากจุดที่ mean model วางไว้, $V_i^{-1}$ ตัดสินว่า residual ของ cluster นี้แต่ละตัวควรมีน้ำหนักเท่าไร และ $D_i$ แปลผลลัพธ์กลับมาอยู่บนสเกลของสัมประสิทธิ์ ดัชนี $i$ วิ่งไปตาม cluster ที่เป็นอิสระต่อกันจำนวน $N$ ก้อน ซึ่งในตัวอย่างที่เราใช้อยู่คือผู้เข้าร่วม 80 คน ไม่ใช่ 320 แถว คุณสมบัติเชิงกลุ่มตัวอย่างขนาดใหญ่ทุกข้อของวิธีนี้อ้างถึงจำนวนนี้ ไม่ได้อ้างถึงจำนวนการวัด ส่วน working correlation อาศัยอยู่ใน $V_i$ เท่านั้น และ $V_i$ ปรากฏในฐานะน้ำหนักเท่านั้น

นี่คือเหตุผลว่าทำไม working correlation ที่ผิดจึงไม่ทำให้ point estimate พัง มันทำให้เสีย efficiency และทำให้ standard error แบบ model-based ผิดไป ซึ่งเป็นจุดที่องค์ประกอบที่สามได้ชื่อของมันมา

\[ \widehat{\operatorname{Var}}(\widehat{\beta}) = \underbrace{\left(\textstyle\sum D_i'V_i^{-1}D_i\right)^{-1}}_{\text{ขนมปัง}} \underbrace{\left(\textstyle\sum D_i'V_i^{-1}\widehat{r}_{i}\widehat{r}_{i}' V_i^{-1}D_i\right)}_{\text{ไส้}} \underbrace{\left(\textstyle\sum D_i'V_i^{-1}D_i\right)^{-1}}_{\text{ขนมปัง}} \]

สองก้อนนอกคือสิ่งที่คุณจะใช้ถ้า working correlation ถูกต้อง ส่วนก้อนกลางแทนที่ covariance ที่สมมติไว้ด้วยผลคูณไขว้เชิงประจักษ์ของ residual ที่สังเกตได้จริง $\widehat{r}_{i} = Y_i - \widehat{\mu}_{i}$ ทีละ cluster พูดอีกอย่างคือข้อมูลได้รับอนุญาตให้ล้มข้อสมมติเรื่อง variance ส่วนคำถามที่ว่าความ robust นี้คุ้มหรือไม่ในงานหนึ่ง ๆ เป็นการตัดสินใจเรื่องดีไซน์ ไม่ใช่คุณสมบัติของตัวประมาณค่า — ที่ผู้เข้าร่วมจำลอง 80 คนอาจพึ่งพาได้อย่างสมเหตุสมผล แต่ถ้ามี cluster เพียงไม่กี่ก้อนก็อาจพึ่งไม่ได้

เลือก working correlation อย่างไร

ก่อนลุยของยาก: working correlation คือการเดาตั้งต้นของคุณว่าค่าวัดซ้ำ ๆ ของผู้ป่วยคนหนึ่งคล้ายกันแบบไหน จะเป็นไม่คล้ายกันเลย คล้ายกันเท่ากันทุกคู่ คล้ายกันน้อยลงเมื่อห่างกันตามเวลา หรือแต่ละคู่มีค่าของตัวเอง สี่ตัวเลือกข้างล่างก็คือการเดาสี่แบบนี้ ส่วนข้อควรระวังที่ตามมาเป็นเรื่องว่าต้องมีผู้ป่วยกี่คนแซนด์วิชถึงจะไว้ใจได้ว่าจะช่วยกู้การเดาที่ไม่ดีให้

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

Independence — ทำงานราวกับว่าแถวภายใน cluster เดียวกันไม่เกี่ยวข้องกัน ซึ่งตัวมันเองก็เป็นข้อสมมติแบบ working อยู่ดี ไม่ใช่การไม่สมมติอะไรเลย แล้วปล่อยให้แซนด์วิชรับหน้าที่ซ่อมทั้งหมด ใช้ได้จริงและมักใช้การได้ดีกว่าที่คิด แต่สิ้นเปลืองเมื่อ correlation สูง
Exchangeable — correlation ค่าเดียวสำหรับทุกคู่ นี่คือคู่ขนานเชิง working ของ compound symmetry matrix จากตอนที่ 2 ที่ทุกช่องนอกเส้นทแยงมุมเท่ากับ 0.625
AR-1 — correlation ที่ลดลงตามระยะห่าง ต้องมีมิติที่เรียงลำดับได้ ซึ่งสี่ตำแหน่งบนร่างกายที่สลับกันได้ไม่ได้ให้มา จึงเป็นเรื่องของ visit sub-study ในตอนที่ 5 ไม่ใช่ที่นี่
Unstructured — correlation อิสระหกค่าเมื่อ $k=4$ ประมาณค่าโดยไม่บังคับโครงสร้างใด ๆ เพิ่มเติม แต่ยังต้องถูกจำกัดให้เมทริกซ์ working ยังคง positive definite แลกกับพารามิเตอร์หกตัวและความต้องการจำนวน cluster ที่มากขึ้น

ข้อควรระวังที่ซื่อสัตย์ที่สุดผูกอยู่กับแซนด์วิช ไม่ใช่กับตัวเลือกข้างต้น Robust variance estimator เป็นคำสัญญาที่ตั้งอยู่บนจำนวน cluster ที่มาก: มัน consistent เมื่อจำนวน cluster อิสระเพิ่มขึ้น และเป็นที่ทราบกันว่าเมื่อ cluster น้อย มันจะ anti-conservative คือให้ช่วงความเชื่อมั่นที่แคบกว่าที่ควรจะเป็น ดังนั้นเมื่อจำนวน cluster น้อย ควรพิจารณา robust variance แบบที่มีการแก้ bias ซึ่งมีให้ใช้ในโปรแกรมส่วนใหญ่ และ mixed model แบบ likelihood-based อาจเป็นทางเลือกที่ปกป้องได้ง่ายกว่า

ความบังเอิญของ linear model

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

ตอนที่ 3 สัญญาไว้ว่า random intercept ขยับ covariance ไม่ได้ขยับ mean คำสัญญานั้นมีเงื่อนไขห้อยท้ายอยู่ — สำหรับ linear model — และตอนนี้คือตอนที่เงื่อนไขนั้นถึงกำหนดชำระ

เมื่อ link เป็น identity ค่าเฉลี่ยแบบ conditional ของการวัดหนึ่งครั้งคือ $\mu_{j,g} + b_i$ และค่าเฉลี่ยแบบ marginal ได้จากการเฉลี่ยข้าม person effect เนื่องจาก $E(b_i)=0$ โดยการสร้าง

\[ E(Y_{ij}) = E_b\!\left[\mu_{j,g} + b_i\right] = \mu_{j,g} + E(b_i) = \mu_{j,g} \]

สัมประสิทธิ์ fixed effect ของโมเดล conditional และ marginal จึงอธิบายปริมาณเดียวกันบนสเกลเดียวกัน นี่คือทั้งหมดของความบังเอิญนี้ และมันเป็นเรื่องพีชคณิต ไม่ใช่โชคของชุดข้อมูล: การเฉลี่ยเส้นตรงข้ามค่ารบกวนที่มีค่าเฉลี่ยเป็นศูนย์ ย่อมได้เส้นตรงเส้นเดิมกลับมา ส่วนตัวแปรร่วมระดับบุคคลที่โค้ดในซีรีส์นี้ใส่ไว้ — age, sex และ phototype — มีค่าคงที่ภายในคนคนเดียวกัน อยู่ใน mean structure ล้วน ๆ และไม่กระทบข้อโต้แย้งข้างต้นเลย บนข้อมูลจำลองที่ outcome เป็นค่าต่อเนื่อง ผลต่าง AD − control ที่ hand เท่ากับ +8.0 g·m⁻²·h⁻¹ ไม่ว่าคุณจะ fit lme ด้วย random = ~ 1 | id, gls ด้วย unstructured covariance หรือ GEE ด้วย exchangeable working correlation ทั้งสามวิธีเห็นตรงกันเรื่อง estimand และต่างกันแค่วิธีลงบัญชี correlation เท่านั้น แต่ทั้งสามไม่จำเป็นต้องให้ standard error เท่ากัน เพราะกลไกด้าน variance ต่างกันจริง

นี่คือเหตุผลที่นักวิเคราะห์จำนวนมากไม่เคยต้องคิดเรื่องความแตกต่างนี้เลย outcome แบบต่อเนื่องให้อภัยเรา

จุดที่โชคหมดลง

ก่อนลุยของยาก: ทีนี้เปลี่ยนจากค่าที่วัดได้เป็น outcome แบบมี/ไม่มี ความบังเอิญเมื่อครู่จะพังทันที เพราะความเสี่ยงและ odds ไม่ได้เฉลี่ยแบบเดียวกับค่าเฉลี่ยทั่วไป ลองนึกถึงผู้ป่วยสี่คนที่มีความเสี่ยงต่อผิวแห้งเฉพาะตัวเท่ากับ 5%, 20%, 60% และ 95% ค่าเฉลี่ยของความเสี่ยงสี่ค่านี้ไม่เท่ากับความเสี่ยงของผู้ป่วยที่อยู่ตรงกลาง เพราะสเกลความเสี่ยงถูกบีบที่ปลายทั้งสองข้าง การเฉลี่ยความสัมพันธ์ที่เป็นเส้นโค้งข้ามผู้ป่วยหลายคนจึงทำให้เส้นแบนลง และเส้นที่แบนลงก็คือ odds ratio ที่เล็กลง ทุกอย่างข้างล่างคือการคิดบัญชีอย่างละเอียดว่าเล็กลงเท่าไร

ทีนี้ลองเปลี่ยน outcome ดู Extension C ของการจำลองใช้ผู้เข้าร่วม 80 คนเดิมและสี่ตำแหน่งเดิม แต่บันทึกว่ามีผิวแห้งที่มองเห็นได้ — xerosis — ที่ตำแหน่งนั้นหรือไม่ โมเดลที่สร้างข้อมูลคือ logistic random-intercept model

\[ \operatorname{logit}\big[P(\text{xerosis}_{ij}=1 \mid b_i)\big] = \beta_0 + \beta_1\,\text{AD}_i + b_i, \qquad b_i \sim N(0,\ \tau^2) \]

โดยในการจำลองกำหนด $\beta_0 = -1.00$, $\beta_1 = +1.80$ และ $\tau^2 = 2.00$ โครงสร้างแบบ subject-specific นี้คือโครงสร้างเดียวกับตอนที่ 1 เป๊ะ ๆ — หนึ่งค่าต่อหนึ่งคน ใส่ซ้ำแบบไม่เปลี่ยนแปลงในทุกแถวของคนคนนั้น — เพียงแต่ตอนนี้ใส่บนสเกล logit เมื่อ exponentiate $\beta_1$ จะได้ conditional odds ratio เท่ากับ 6.05

ขั้นตอนการเฉลี่ยที่เมื่อครู่ยังไม่มีพิษภัย ตอนนี้มีแล้ว เพราะฟังก์ชัน inverse-logit ไม่เป็นเชิงเส้น

\[ E_b\!\left[\operatorname{expit}(\eta + b)\right] \;\neq\; \operatorname{expit}\!\left(\eta + E_b[b]\right) \]

การเฉลี่ยเส้นโค้งรูปตัว S ข้ามการกระจายของ person effect ไม่ได้คืนเส้นโค้งรูปเดิมกลับมา แต่คืนเส้นที่แบนราบกว่า เส้นโค้งแบบ population-averaged จึงชันน้อยกว่าเส้นโค้งแบบ subject-specific และสัมประสิทธิ์แบบ population-averaged จึงมีขนาดเล็กกว่าในเชิงสัมบูรณ์ สูตรประมาณมาตรฐานของตัวคูณที่ทำให้หดลง ซึ่งเสนอโดย Zeger, Liang และ Albert คือ

\[ \beta^{M} \approx \beta^{C}\left(1 + c^2\tau^2\right)^{-1/2}, \qquad c = \frac{16\sqrt{3}}{15\pi} \approx 0.588 \]

ซึ่งแปลงสัมประสิทธิ์ conditional เป็นคู่ขนานแบบ population-averaged โดยใช้เพียง between-person variance ค่าเดียว ลองเดินเลขคณิตทีละขั้นด้วยค่าจากการจำลอง

\[ c^2 = 0.588^2 = 0.346, \qquad 1 + 0.346\times 2.00 = 1.692, \qquad \sqrt{1.692} = 1.30 \]

\[ \beta^{M} \approx \frac{1.80}{1.30} \approx 1.38 \qquad\Rightarrow\qquad \text{marginal OR} = e^{1.38} \approx 4.0 \]

เมื่อใช้ตัวคูณเดียวกันนี้กับ intercept จะได้ความเสี่ยงแบบ marginal ที่สูตรประมาณบอกเป็นนัย: $-1.00/1.30 = -0.77$ สำหรับกลุ่ม control และ $(-1.00+1.80)/1.30 = 0.62$ สำหรับกลุ่ม AD ซึ่งเมื่อแปลงเป็นสเกลความน่าจะเป็นได้ประมาณ 32% และ 65% การประกอบความเสี่ยงสองค่านี้กลับเป็น odds ratio ให้ค่า 4.0 อีกครั้ง ต้องพูดตรง ๆ ว่านี่ไม่ใช่การตรวจทานอะไรทั้งสิ้น ทั้งสองเส้นทางหาร linear predictor ตัวเดียวกันด้วย $\sqrt{1+c^2\tau^2}$ ตัวเดียวกัน เส้นทางที่สองจึงยุบลงทางพีชคณิตเหลือ $e^{\beta_1/\sqrt{1+c^2\tau^2}}$ ซึ่งก็คือเส้นทางแรกที่ถูกเขียนออกมาบนสเกลความน่าจะเป็นแล้วอ่านกลับ ความสอดคล้องกันนี้จึงถูกบังคับด้วยการสร้าง ไม่ใช่สิ่งที่ได้มาจากการตรวจสอบ และซีรีส์ที่ยืนยันตลอดว่าผลลัพธ์ที่โปรแกรมพิมพ์ออกมายืนยันข้อสมมติที่สร้างมันไม่ได้ ก็ไม่ควรรับข้อความที่เป็นจริงโดยนิยามมาเป็นหลักฐานเช่นกัน

การตรวจทานของจริงคือหยุดใช้สูตรประมาณแล้วอินทิเกรตออกมาตรง ๆ การเฉลี่ย $\operatorname{expit}(\beta_0 + \beta_1\text{AD}_i + b)$ ข้าม $b \sim N(0,\ 2.00)$ ด้วยวิธีเชิงตัวเลข ให้ความเสี่ยงแบบ marginal ที่แน่นอนเท่ากับ 32.5% และ 64.2% ให้ odds แบบ marginal เท่ากับ 0.481 และ 1.792 และให้ population-averaged odds ratio ที่แน่นอนเท่ากับ 3.72 สูตรประมาณจึงให้ค่า marginal odds ratio สูงเกินไปราว 7% ที่ค่า $\tau^2$ ระดับนี้ — น้อยพอที่ประเด็นเชิงการสอนยังยืนอยู่ได้ และมากพอที่ 4.0 กับ 3.72 ต้องไม่ถูกวางในตารางราวกับเป็นปริมาณเดียวกัน

ความจริงจำลองชุดเดียว สัมประสิทธิ์สองค่า: 6.05 กับ 3.72

ภาพประกอบ: เส้นโค้ง conditional ที่ชันกว่าและเส้นโค้ง marginal ที่แบนกว่าบนแกนเดียวกัน แสดงการ attenuation

อ่านออกเสียงทั้งสองค่า

ตัวเลขสองค่านี้เลิกสับสนทันทีที่พูดออกมาเป็นประโยคทางคลินิก

Conditional odds ratio พูดว่า สำหรับผู้ป่วยสองคนที่มีแนวโน้มพื้นฐานต่อผิวแห้งเท่ากัน คนที่เป็น atopic dermatitis มี odds ของการเกิด xerosis ที่ตำแหน่งหนึ่ง ๆ ประมาณหกเท่า การเปรียบเทียบถูกตรึงไว้ภายในชั้นของระดับบุคคล ซึ่งเป็นเหตุผลที่มันถูกเรียกว่า subject-specific

Population-averaged odds ratio พูดว่า เมื่อเฉลี่ยข้ามประชากรจำลองทั้งหมด กลุ่ม atopic dermatitis มี odds ของการเกิด xerosis ประมาณ 3.7 เท่า การเปรียบเทียบเป็นการเทียบระหว่างคนสองกลุ่ม และความแตกต่างระหว่างบุคคลถูกเฉลี่ยกลืนหายไป ไม่ได้ถูกตรึงไว้

มีข้อเท็จจริงเชิงเลขคณิตสองข้อที่ควรพูดตรง ๆ เพราะมันจะสลายความลึกลับส่วนใหญ่ทิ้งไป ข้อแรก สัมประสิทธิ์แบบ marginal อยู่ใกล้ค่า null มากกว่าสัมประสิทธิ์แบบ conditional เสมอ และภายในกรอบของสูตรประมาณ ระยะห่างระหว่างสองค่าถูกกำหนดโดย $\tau^2$ ทั้งหมด เมื่อ $\tau^2 = 0$ ตัวคูณกลายเป็น 1 และคำตอบทั้งสองตรงกันพอดี ซึ่งการคำนวณแบบแน่นอนก็ยืนยันข้อนี้ ข้อที่สอง เมื่อคำนวณแบบแน่นอน ขนาดของการหดลงไม่ได้ขึ้นกับ $\tau^2$ อย่างเดียว แต่ขึ้นกับตำแหน่งบนเส้นโค้ง logistic ที่การเปรียบเทียบนั้นตั้งอยู่ด้วย นั่นคือขึ้นกับ $\beta_0$ ด้วย ความขึ้นต่อกันข้อนี้เองที่ทำให้ 4.0 กับ 3.72 ไม่ใช่ค่าเดียวกัน การหดลงตรงนี้ไม่ใช่ bias ไม่ใช่ confounding และไม่ใช่ผลจากกลุ่มตัวอย่างเล็ก — มันคือสิ่งที่เกิดขึ้นเมื่อเราเฉลี่ยฟังก์ชันที่ไม่เป็นเชิงเส้น

คำถามวิจัยของคุณต้องการ estimand ตัวไหน

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

การเปรียบเทียบระดับประชากรและระดับนโยบายเอนไปทาง marginal ถ้าข้อสรุปที่คุณตั้งใจจะตีพิมพ์อยู่ในรูป กลุ่มนี้มีความชุกสูงกว่ากลุ่มนั้น หรือ ถ้าประชากรทั้งคลินิกได้รับการดูแลต่างออกไป ค่าเฉลี่ยของ outcome จะขยับเท่านี้ estimand ที่คุณต้องการคือแบบ population-averaged และ marginal model เล็งไปที่มันโดยตรง ในทางกลับกัน การทำนายระดับบุคคล การแยกส่วน variance การรายงาน ICC การใช้ BLUPs และคำถามใด ๆ ที่ความแตกต่างระหว่างบุคคลคือตัววัตถุทางวิทยาศาสตร์เอง ล้วนต้องใช้รูปแบบ random-effects เพราะปริมาณเหล่านั้นไม่มีอยู่เลยใน marginal model

Estimand แบบ subject-specific อาจเป็นตัวเลือกที่เหมาะกว่าเมื่อประโยคทางคลินิกเป็นเรื่องภายในคนจริง ๆ เช่น การตัดสินใจ step-down สำหรับผู้ป่วยที่อยู่ในการติดตามอยู่แล้ว มากกว่าจะเป็นการเทียบประชากรสองกลุ่ม เมื่อ estimand สองตัวใกล้กัน ซึ่งจะเป็นเช่นนั้นเสมอเมื่อ $\tau^2$ เล็ก การเลือกอาจตัดสินด้วยความสะดวกได้อย่างสมเหตุสมผล แต่เมื่อสองค่าห่างกันมาก การเลือกต้องตัดสินด้วยคำถามวิจัย และควรกำหนดไว้ล่วงหน้า ไม่ใช่มาตัดสินหลังจากเห็นแล้วว่าตัวไหนใหญ่กว่า

สิ่งที่ GEE ให้คุณไม่ได้

มีสี่อย่างที่ไม่มีอยู่เลย และรู้ไว้ก่อนย่อมดีกว่าไปพบตอนกำลังเขียนผลการศึกษา หนึ่ง ไม่มี $\tau^2$: working correlation parameter ไม่ใช่การแยกส่วน variance จึงไม่มีการแบ่งระหว่างบุคคลกับภายในบุคคลให้ดู สอง ด้วยเหตุนั้นจึงไม่มี ICC ในความหมายของการแบ่งสัดส่วน variance แบบตอนที่ 2 สาม ไม่มี BLUPs จึงไม่มีการทำนายระดับบุคคล สุดท้าย ไม่มี likelihood ซึ่งแปลว่าไม่มี likelihood-ratio test และไม่มี AIC การเปรียบเทียบโมเดลต้องใช้ QIC แทน

ถ้า ICC คือเป้าหมายทางวิทยาศาสตร์ ไม่ใช่ nuisance parameter marginal model ก็เป็นเครื่องมือที่ผิด — และควรจำไว้ด้วยว่า ICC ของ reliability ในการศึกษาด้านการวัด เป็นคนละปริมาณอีกตัวหนึ่ง ที่ใช้สูตรร่วมกันกับ ICC ของ clustering แต่ตอบคนละคำถาม

ข้อควรระวังเรื่องข้อมูลขาดหาย

ก่อนลุยของยาก: มีความต่างเชิงปฏิบัติข้อหนึ่งที่มักตัดสินการเลือกก่อนข้อถกเถียงเรื่อง estimand เสียอีก และเป็นเรื่องของคนที่หลุดจากการติดตาม ลองถามเป็นคำถามทางคลินิกตรง ๆ ว่า คนที่ขาดนัดไปนั้นเราดูออกจากสิ่งที่มีอยู่ในเวชระเบียนแล้วหรือไม่ เช่น ค่าครั้งก่อนแย่กว่า อยู่ในกลุ่มที่โรครุนแรงกว่า หรืออายุมากกว่า ถ้าดูออก mixed model แบบ likelihood-based รับมือได้ด้วยตัวเอง ในขณะที่ GEE แบบธรรมดาต้องการให้การขาดหายนั้นแทบไม่เกี่ยวกับอะไรเลย ย่อหน้าข้างล่างพูดเรื่องนี้อย่างเป็นทางการ

ยังมีความแตกต่างอีกข้อหนึ่งที่บ่อยครั้งเป็นตัวตัดสิน และไม่เกี่ยวกับ estimand เลย

Mixed model แบบ likelihood-based ใช้ได้ภายใต้ข้อสมมติ missing at random กล่าวคือการขาดหายอาจขึ้นกับข้อมูลที่สังเกตได้ ส่วน GEE ที่ใช้ sandwich variance แบบปกติต้องการข้อสมมติที่แข็งกว่าคือ missing completely at random เว้นแต่จะเพิ่ม inverse-probability weight เข้าไป ในการศึกษาทางคลินิกแบบ longitudinal การหลุดหายจากการติดตามมักสัมพันธ์กับประวัติที่สังเกตได้ คือคนที่ขาดนัดมักเป็นคนที่ค่าที่วัดครั้งก่อนแย่ที่สุด รูปแบบนี้เข้าข่าย missing at random ไม่ใช่ missing completely at random ข้อสมมติที่อ่อนกว่าจึงมักเป็นข้อสมมติที่สมเหตุสมผลกว่า ดังนั้นเมื่อคาดว่าดีไซน์จะสูญเสียผู้เข้าร่วมในลักษณะที่ขึ้นกับข้อมูลที่บันทึกไว้แล้ว mixed model แบบ likelihood-based หรือ GEE ที่ถ่วงน้ำหนักอย่างเหมาะสมอาจเป็นการระบุโมเดลที่ปลอดภัยกว่า และเหตุผลนี้ควรถูกเขียนไว้ในแผนการวิเคราะห์ ไม่ใช่ไปแก้ต่างเอาทีหลัง

ส่วนการหลุดหายที่ขึ้นกับค่าที่ไม่ได้ถูกสังเกตเอง — missing not at random ซึ่งเป็นความหมายที่คำว่า informative ถูกสงวนไว้ให้ตามธรรมเนียม — ไม่ได้ถูกซ่อมด้วยโมเดลใดในสองแบบนี้ กรณีนั้นต้องอาศัยการระบุข้อสมมติเกี่ยวกับกลไกการขาดหายอย่างชัดเจน พร้อม sensitivity analysis รอบข้อสมมตินั้น และการเลือกระหว่าง mixed model กับ GEE ก็ไม่ใช่การตัดสินใจหลักอีกต่อไป

สิ่งที่บทความนี้ไม่ได้อนุญาตให้สรุป

มีการอ่านเกินสามแบบที่เอื้อมถึงได้ง่ายจากตรงนี้ และทั้งสามแบบผิด

หนึ่ง ความบังเอิญของ identity link ไม่ใช่กฎทั่วไปของโมเดลที่ไม่เป็นเชิงเส้น มันจำเพาะกับ identity link และสิ่งที่ใกล้เคียงข้อยกเว้นที่สุด — log link ซึ่ง person effect ถูกแยกออกมาเป็นตัวคูณ และมีเพียง intercept ที่ดูดซับการเลื่อนไว้ — เป็นการยืนยันประเด็นนี้ ไม่ใช่ทำให้อ่อนลง สำหรับ logit link สัมประสิทธิ์ต่างกันจริง และทั้ง 6.05 และ 3.72 ถูกต้องทั้งคู่

สอง สูตรการหดลงเป็นสูตรประมาณ และเขียนด้วย $\approx$ ด้วยเหตุผลนั้น บนค่าที่ใช้จำลอง มันให้ population-averaged odds ratio สูงเกินไปราว 7% ซึ่งแม่นพอสำหรับ $\tau^2$ ขนาดปานกลางที่จะสื่อประเด็นเชิงการสอนและใช้ตรวจทานโมเดลที่ fit ได้อย่างคร่าว ๆ แต่ไม่ควรใช้แปลง conditional odds ratio ที่ตีพิมพ์แล้วให้เป็น marginal แล้วรายงานราวกับว่าเป็นค่าที่ประมาณมาเอง

สาม คำว่า population-averaged เป็นคำที่พูดถึงการเฉลี่ยข้าม random effect เท่านั้น ไม่ได้แปลว่า marginal ในความหมายเชิงสาเหตุ ไม่ได้แปลว่าจัดการ confounding แล้ว และความสัมพันธ์แบบ population-averaged จากดีไซน์เชิงสังเกตก็ยังคงเป็นความสัมพันธ์อยู่ดี การอ่านเชิงสาเหตุใด ๆ ยังต้องอาศัยดีไซน์และชุดตัวแปรที่ปรับมารองรับอยู่เช่นเดิม

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

หนึ่ง เขียน estimand ออกมาเป็นประโยคทางคลินิกก่อนเลือกโปรแกรม แล้วตรวจว่าประโยคนั้นมีคำว่า สำหรับผู้ป่วยคนหนึ่ง หรือคำว่า ในประชากร อยู่ สอง สำหรับ outcome ต่อเนื่องจากดีไซน์ที่สมดุล ทั้งสองรูปแบบปกป้องได้ทั้งคู่ จึงควรระบุว่าคุณ fit ตัวไหนและเพราะอะไร แทนที่จะปล่อยให้ผู้อ่านเดาเอง สาม สำหรับ outcome แบบ binary หรือ count ให้กำกับ odds ratio ทุกค่าในตารางผลว่าเป็น conditional หรือ population-averaged เพราะทั้งสองแบบถูกวางเรียงข้างกันข้ามงานวิจัยราวกับเป็นปริมาณเดียวกันอยู่เป็นประจำ สี่ ถ้ารายงาน conditional odds ratio ให้รายงาน $\widehat{\tau}^{2}$ ควบคู่ไปด้วย เพราะระยะห่างระหว่างสองสเกลเป็นฟังก์ชันของค่านั้น และผู้อ่านสร้างมันขึ้นใหม่เองไม่ได้ ห้า เมื่อคาดว่าการหลุดหายจากการติดตามจะขึ้นกับข้อมูลที่สังเกตไว้แล้ว ให้เลือก mixed model แบบ likelihood-based หรือ GEE ที่ถ่วงน้ำหนัก และระบุไว้ในโปรโตคอล ส่วนกรณีที่การหลุดหายอาจขึ้นกับค่าที่ไม่ได้ถูกสังเกตเอง ให้วางแผน sensitivity analysis ไว้ตั้งแต่ตอนนั้นเลย สุดท้าย เมื่อจำนวน cluster น้อย ให้พิจารณา robust variance ที่แก้ bias แล้ว และอ่านค่า coverage ของช่วงความเชื่อมั่นจากแซนด์วิชที่ไม่ได้แก้ด้วยความระมัดระวัง

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

Marginal model ไม่อธิบายว่า correlation ภายในคนมาจากไหน มันระบุ covariance ออกมาตรง ๆ และไม่ว่าจะเลือกทางไหน ความไม่เป็นอิสระก็ถูกจัดการครบถ้วนเท่ากัน
GEE มีองค์ประกอบสามอย่างพอดี — mean model, working correlation structure และ sandwich variance estimator — และ point estimate ยังคง consistent เมื่อสิ่งที่ผิดมีเพียง working correlation
เมื่อ link เป็น identity สัมประสิทธิ์ fixed effect แบบ conditional และ marginal ตรงกัน เพราะ $E(b_i)=0$ ซึ่งเป็นเหตุผลที่ outcome ต่อเนื่องให้อภัยความแตกต่างนี้ทั้งหมด
เมื่อ link เป็น logit สัมประสิทธิ์ไม่ตรงกัน: ในตัวอย่างจำลอง โมเดลที่สร้างข้อมูลตัวเดียวกันให้ conditional OR เท่ากับ 6.05 และ population-averaged OR ที่แน่นอนเท่ากับ 3.72 พร้อมความเสี่ยงแบบ marginal ประมาณ 32% และ 64% ส่วนสูตรประมาณแบบ Zeger–Liang–Albert ให้ค่าหลังนี้ประมาณ 4.0 ซึ่งสูงเกินไปราว 7% และต้องกำกับให้ชัดว่าค่าไหนเป็นค่าไหน
ภายในกรอบของสูตรประมาณ การหดลงถูกกำหนดโดย $\tau^2$ เพียงตัวเดียว มันดึงค่าประมาณแบบ marginal เข้าหา null เสมอ และหายไปพอดีเมื่อ $\tau^2 = 0$
GEE ไม่ให้ $\tau^2$ ไม่ให้ ICC ไม่ให้ BLUPs และไม่ให้ likelihood การเลือกระหว่างสองตระกูลนี้จึงควรตัดสินด้วย estimand และข้อสมมติเรื่องข้อมูลขาดหาย ไม่ใช่ด้วยความเคยชิน

ตอนต่อไป. ตอนที่ 8 จะเอาทุกอย่างที่ซีรีส์นี้สร้างมาเล็งไปที่คอลัมน์ผลลัพธ์เพียงคอลัมน์เดียว: standard error สี่ค่าที่เท่ากันเป๊ะ กับอาจารย์ที่อยากรู้ว่าโปรแกรมเสียหรือเปล่า

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

Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982;38(4):963–974.
Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13–22.
Zeger SL, Liang KY, Albert PS. Models for longitudinal data: a generalized estimating equation approach. Biometrics. 1988;44(4):1049–1060.
Verbeke G, Molenberghs G. Linear Mixed Models for Longitudinal Data. Springer; 2000.
Diggle PJ, Heagerty P, Liang KY, Zeger SL. Analysis of Longitudinal Data. 2nd ed. Oxford University Press; 2002.
Hubbard AE, Ahern J, Fleischer NL, et al. To GEE or not to GEE: comparing population average and mixed models for estimating the associations between neighborhood risk factors and health. Epidemiology. 2010;21(4):467–474.