เมื่อ Compound Symmetry ไม่จริง: variance ไม่เท่ากัน, unstructured correlation และเพดาน identifiability

บทคัดย่อ

หลายคนเชื่อว่า random intercept แปลว่าทุกคู่ของการวัดซ้ำจากคนเดียวกัน correlate เท่ากัน ประโยคนี้ไม่ผิดแต่ถูกตัดครึ่ง สิ่งที่มันรับประกันคือการเติม between-person variance ค่าเดียวกันให้ covariance ทุกคู่ ส่วน correlation จะเท่ากันต้องอาศัยข้อสองที่มักไม่มีใครพูด คือ residual ต้องไม่ correlate กันและมี variance เท่ากัน บทความนี้ผ่อนครึ่งหลังนั้น เมื่อให้ residual variance ต่างกันตามตำแหน่ง covariance ยังเป็น 6.25 นอกเส้นทแยงมุมทุกช่อง แต่ correlation ทั้งหกไล่จาก 0.44 ถึง 0.74 แล้วนับต้นทุน covariance matrix สมมาตร 4×4 มีปริมาณอิสระสิบค่า unstructured correlation บวก variance ที่ไม่เท่ากันใช้ครบสิบ ถ้าใส่ variance ระดับบุคคลทับอีกก็เป็นสิบเอ็ดพารามิเตอร์ต่อสิบปริมาณ

---

หนึ่งค่าต่อหนึ่งคน: เข้าใจ Linear Mixed Model ตั้งแต่รากฐาน — ตอนที่ 6 จาก 8

เพิ่งมาเจอตอนนี้เป็นตอนแรกใช่ไหม ขอปูพื้นสั้น ๆ ก่อน เมื่อเราวัดคนไข้คนเดิมหลายครั้ง ค่าที่วัดได้จะไม่เป็นอิสระต่อกัน และ mixed model จัดการเรื่องนี้ด้วยการแจกตัวเลขส่วนตัวให้คนไข้แต่ละคนหนึ่งค่า เป็นค่าชดเชยที่ยกโปรไฟล์ของคนคนนั้นขึ้นหรือลงทั้งเส้น ทุกอย่างที่โมเดลสมมติเกี่ยวกับการวัดซ้ำถูกอัดอยู่ในบรรทัดเดียวคือ $V_i = Z_i G Z_i' + R_i$ ซึ่งมีของที่เปลี่ยนได้เพียงสองชิ้น คือ $Z$ ที่กำหนดว่าแต่ละคนได้พกอะไรติดตัวมาบ้าง และ $R$ ที่กำหนดว่าการวัดครั้งเดียวแบกอะไรไว้เอง โพสต์นี้เป็นเรื่องของ $R$ และเป็นเรื่องของประโยคหนึ่งที่พวกเราส่วนใหญ่ถูกสอนมาในเวอร์ชันที่จริงแค่ครึ่งเดียว ถ้าอยากได้ภาพรวมทั้งแปดตอนแบบไม่ต้องลุยสมการ อ่านได้ที่ คู่มือซีรีส์

ตอนที่ 5 เราเติมความรวยให้พจน์ด้านซ้ายของ $V_i = Z_i G Z_i' + R_i$ ด้วยการเพิ่มคอลัมน์ที่สองให้ $Z$ และแต่ละคนก็เลิกพกแค่ "ค่าชดเชย" เปลี่ยนเป็นพก "เส้นทาง" ของตัวเองมาด้วย ตอนที่ 6 นี้เราจะปล่อย $Z$ ไว้เฉย ๆ แล้วหันไปเล่นงานพจน์ด้านขวาแทน คือ $R_i$ หรือ residual covariance ซึ่งเป็นก้อนที่งานวิเคราะห์เกือบทุกชิ้นตรึงไว้ที่ $\sigma^2 I$ อย่างเงียบ ๆ โดยไม่เคยเขียนสมมติฐานนี้ลงในเมธอดเลย เราจะไล่พจน์นี้ขึ้นบันได covariance structure ไปจนชนเพดานแข็ง ๆ แล้ว ตอนที่ 7 จะพาเดินออกทางประตูบานสุดท้ายที่เหลืออยู่ในห้อง

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

ประโยคที่เรากำลังจะเอามาส่องดู

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

ผู้อ่านเกือบทุกคนมาถึงตรงนี้พร้อมประโยคหนึ่งที่ถูกติดตั้งไว้ในหัวแล้วว่า "random intercept แปลว่าการวัดซ้ำทุกคู่ในคนเดียวกัน correlate กันเท่ากันหมด" ประโยคนี้ไม่ได้ผิด แต่มันถูกตัดครึ่งหลังทิ้งไป และครึ่งหลังนั้นคือส่วนที่ทำงานหนักที่สุด สิ่งที่ random intercept รับประกันจริง ๆ คือมัน สมทบ ปริมาณเท่ากันเป๊ะเข้าไปใน covariance ของทุกคู่การวัดในคนเดียวกัน นั่นคือ $\tau^2$ ส่วนที่ว่าปริมาณที่สมทบเข้าไปนั้นจะเป็น covariance ทั้งก้อนหรือไม่ เป็นคนละเรื่อง และเป็นเรื่องของ residual โดย covariance ของหนึ่งคู่เท่ากับ $\tau^2 + R_{jk}$ ซึ่งจะยุบเหลือแค่ $\tau^2$ ก็ต่อเมื่อ residual ของการวัดทั้งสองไม่มี correlation ต่อกัน ส่วนการที่ correlation เท่ากันด้วยนั้นเป็นข้อกล่าวอ้างที่ไกลออกไปอีกขั้นหนึ่ง และมันจะจริงก็ต่อเมื่อสมมติฐานอีกข้อหนึ่ง ซึ่งไม่มีใครในห้องพูดออกมาดัง ๆ บังเอิญเป็นจริงไปด้วย

ประโยคที่ถูกต้องคือ random intercept บวกกับ residual ที่เป็นอิสระต่อกันและมี variance เท่ากัน จึงจะได้ compound symmetry ถ้าตัดครึ่งหลังทิ้ง ครึ่งแรกยังอยู่ครบไม่สะเทือน แต่ compound symmetry ตายเงียบ ๆ ไปแล้ว ในโพสต์นี้เราจะเอาแว่นขยายไปส่องครึ่งหลัง นับดูว่าการผ่อนคลายสมมติฐานนั้นมีราคาเท่าไร แล้วไปหาจุดที่ covariance matrix ไม่เหลือที่ว่างให้ parameter ตัวถัดไปอีกแล้ว

Spec B: offset ที่แชร์ค่าเดิม แต่ residual variance สี่ค่าไม่เท่ากัน

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

Spec A ซึ่งพาเรามาตั้งแต่ตอนที่ 1 ถึงตอนที่ 4 ให้ทั้งสี่ตำแหน่งวัดใช้ residual variance ค่าเดียวกัน คือ $\sigma^2 = 3.75$ ทุกตำแหน่ง วางอยู่ใต้ between-person variance $\tau^2 = 6.25$ รวมเป็น marginal variance เท่ากับ 10.00 ที่ทุกตำแหน่ง Spec B เป็นการจำลองคนละชุด ไม่ใช่การวิเคราะห์ข้อมูล Spec A ซ้ำ โดยใช้ดีไซน์เดียวกันบนผู้เข้าร่วม 80 คน ตำแหน่งวัดสี่ตำแหน่งเรียงลำดับเดิม ค่าเฉลี่ยประชากรชุดเดิม และ $\tau^2 = 6.25$ ค่าเดิม แต่แต่ละตำแหน่งได้รับ residual variance เป็นของตัวเอง

ตำแหน่ง (site)	residual variance $\sigma_j^2$ (g²·m⁻⁴·h⁻²)	marginal variance $\tau^2+\sigma_j^2$ (g²·m⁻⁴·h⁻²)
forearm	1.75	8.00
hand	2.75	9.00
shin	13.75	20.00
back	3.75	10.00

ข้อสังเกตคือ back ถูกกำหนดให้มี residual variance เท่ากับ 3.75 พอดี ซึ่งเป็นตัวเลขเดียวกับที่ Spec A ใช้กับทั้งสี่ตำแหน่ง ดังนั้นจึงมีเพียงตำแหน่งเดียวเท่านั้นที่ได้ผลรวมคุ้นตาที่ 10.00 ส่วน shin ถูกกำหนดให้มี residual variance สูงที่สุดโดยการสร้างข้อมูล นั่นเป็นคุณสมบัติของการจำลองล้วน ๆ ไม่ใช่ข้อกล่าวอ้างว่าผิวหน้าแข้งมีพฤติกรรมอย่างไรในความเป็นจริง ทีนี้ลองสมมติว่านำโมเดลแบบ pooled variance ไป fit กับข้อมูล Spec B ค่า $\widehat{\sigma}^{2}$ ค่าเดียวที่รายงานออกมาจะตกอยู่ใกล้ ๆ ค่าเฉลี่ยของ variance ทั้งสี่ตำแหน่ง ซึ่งในดีไซน์ที่สมดุลก็คือค่าเฉลี่ยเลขคณิตธรรมดา $22/4 = 5.5$ ไม่ใช่ 3.75 ค่าประมาณนั้นไม่ได้คำนวณผิด เพียงแต่มันกำลังตอบคำถามเกี่ยวกับค่าเฉลี่ยที่ไม่มีตำแหน่งไหนเป็นเจ้าของจริง ๆ

จุดหักมุม: covariance เท่ากันเป๊ะ แต่ correlation ต่างกันหกค่า

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

ทีนี้เอาแว่นขยายไปส่องช่องนอกเส้นทแยงมุม การพิสูจน์จากตอนที่ 1 ไม่ได้เปลี่ยนไปเลยแม้แต่บรรทัดเดียว

\[ \operatorname{Cov}(Y_{ij}, Y_{ik}) = \operatorname{Cov}(b_i + e_{ij},\ b_i + e_{ik}) = \operatorname{Var}(b_i) = \tau^2 = 6.25 \]

อ่านฝั่งขวาให้ดี ไม่มี $\sigma_j^2$ โผล่มาเลย และไม่มีดัชนีตำแหน่งอยู่ในนั้นด้วย residual variance ถูกตัดออกไปจาก covariance ตั้งแต่ก่อนเริ่มคำนวณ เพราะ $e_{ij}$ กับ $e_{ik}$ เป็นอิสระต่อกัน และทั้งคู่เป็นอิสระจาก $b_i$ ช่องนอกเส้นทแยงมุมทุกช่องของ covariance matrix ใน Spec B จึงอ่านได้ว่า 6.25 เหมือนกับใน Spec A ทุกประการ offset ที่ถูกแชร์ยังคงถูกแชร์อยู่ และยังแชร์ด้วยขนาดเท่าเดิม

correlation เป็นคนละปริมาณ เพราะมันเอา covariance นั้นไปหารด้วย standard deviation ของสองตำแหน่งที่ต่างกัน

\[ \operatorname{Corr}(Y_{ij}, Y_{ik}) = \frac{\tau^2}{\sqrt{(\tau^2+\sigma_j^2)(\tau^2+\sigma_k^2)}} \]

ตัวเศษเป็นค่าคงที่ ตัวส่วนไม่ใช่ ลองคำนวณด้วยมือสักสองช่องจากหกช่อง แล้วกลไกจะชัดจนปฏิเสธไม่ได้

\[ \operatorname{Corr}(\text{forearm},\text{hand}) = \frac{6.25}{\sqrt{8.00\times 9.00}} = \frac{6.25}{\sqrt{72}} = \frac{6.25}{8.485} = 0.74 \]

\[ \operatorname{Corr}(\text{forearm},\text{shin}) = \frac{6.25}{\sqrt{8.00\times 20.00}} = \frac{6.25}{\sqrt{160}} = \frac{6.25}{12.649} = 0.49 \]

เลข 6.25 ตัวเดียวกันเข้าไปในทั้งสองบรรทัด แล้วออกมาเป็นคนละค่ากันคนละโลก ทำแบบเดียวกันให้ครบทั้งหกคู่

คู่ตำแหน่ง	ตัวส่วน $\sqrt{(\tau^2+\sigma_j^2)(\tau^2+\sigma_k^2)}$	correlation
forearm–hand	$\sqrt{8\times9}$	0.74
forearm–shin	$\sqrt{8\times20}$	0.49
forearm–back	$\sqrt{8\times10}$	0.70
hand–shin	$\sqrt{9\times20}$	0.47
hand–back	$\sqrt{9\times10}$	0.66
shin–back	$\sqrt{20\times10}$	0.44

offset ที่แชร์กันค่าเดียว covariance ที่เท่ากันเป๊ะทุกช่อง และ correlation ที่ต่างกันหกค่าตั้งแต่ 0.44 ถึง 0.74 compound symmetry เป็นเท็จอย่างเป็นทางการในชุดข้อมูลจำลองนี้ และ random intercept ไม่ได้มีส่วนในการตายของมันเลย คนร้ายอยู่ที่ตัวส่วนล้วน ๆ

varIdent() ทำอะไรกันแน่

ก่อนลุยของยาก: หัวข้อนี้เป็นเรื่องการจดบัญชี ไม่ใช่ทฤษฎี แต่การจดบัญชีนี่แหละที่ทำให้คนอ่าน output แล้วสะดุด nlme ไม่ได้พิมพ์ residual standard deviation ออกมาสี่ค่า มันพิมพ์ออกมาค่าเดียว คือของ level ที่บังเอิญมาเป็นอันดับแรกของ factor แล้วพิมพ์ตัวคูณให้ตัวที่เหลือ residual standard deviation ของหน่วยวัดอื่นจึงเท่ากับค่าที่พิมพ์ออกมานั้น คูณด้วยตัวคูณของหน่วยนั้นเอง

โมเดลที่สร้าง Spec B ถูก fit ด้วยการผูก residual variance หนึ่งค่าเข้ากับแต่ละ level ของ site

# Spec B — site-specific residual variance (Post 6)
fitB <- lme(tewl ~ site * group + age + sex + phototype,
            random = ~ 1 | id, weights = varIdent(form = ~ 1 | site), data = skin)

ใน nlme นั้น varIdent() ไม่ได้ประมาณค่า variance อิสระสี่ค่าโดยตรง แต่ประมาณ residual standard deviation ของ stratum อ้างอิงหนึ่งค่า ซึ่งคือ forearm อันเป็น factor level แรก แล้วประมาณตัวคูณ (multiplier) ให้ stratum ที่เหลือ โดยตรึงตัวคูณของ stratum อ้างอิงไว้ที่ 1 เพื่อให้ระบุค่าได้ ให้อ่าน output ตามนั้น กล่าวคือ ค่า residual standard error ที่พิมพ์ออกมาเป็นของตำแหน่งอ้างอิง ตารางของ variance function เก็บตัวคูณอีกสามค่า และ residual SD ของแต่ละตำแหน่งคือค่าอ้างอิงคูณด้วยตัวคูณของตำแหน่งนั้น เมื่อนับเป็น parameter ก็คือ $1 + (k-1) = 4$ ตัว ในขณะที่โมเดล pooled ใช้เพียง 1 ตัว ส่วน $\tau^2$ ของ random intercept ยังถูกประมาณอยู่บนนั้นเหมือนเดิม และยังรายงานอยู่ในตาราง VarCorr เหมือนเดิม ค่าที่ใช้สร้างข้อมูล Spec B คือ 6.25 และการ fit ครั้งหนึ่ง ๆ จะคืนค่าประมาณที่ตกอยู่ใกล้ค่านั้น ไม่ใช่ตรงกับค่านั้นเป๊ะ

ตัวแปรร่วมระดับบุคคลในคำสั่งนั้น ได้แก่ age, sex และ phototype เข้าไปอยู่ในโครงสร้างค่าเฉลี่ยเท่านั้น ทั้งสามตัวมีค่าคงที่ภายในคนหนึ่งคน จึงอาจดูดซับความแปรปรวนระหว่างบุคคลไปบางส่วนจากที่ $\tau^2$ ต้องแบกไว้ แต่ไม่ได้ไปแตะพีชคณิตของ covariance ในหัวข้อก่อนหน้าเลย

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

อีกทิศทางหนึ่ง: บันไดของ correlation

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

การผ่อนคลายครึ่ง variance ของ compound symmetry เป็นการเดินหมากหนึ่งตา การผ่อนคลายครึ่ง correlation เป็นอีกตาหนึ่ง และโครงสร้างเหล่านี้เรียงตัวเป็นบันไดสั้น ๆ ที่จัดลำดับตามจำนวน parameter ที่ใช้อธิบายบล็อก marginal ขนาด $k \times k$ ของหนึ่งคน คือ $V_i$ การนับข้างล่างนี้เป็นการนับ parameter ของบล็อกที่ประกอบเสร็จแล้วก้อนนั้น ทุกตัวที่อยู่ในนั้น ไม่ว่าจะเข้ามาในรูปของ $\tau^2$ หรือในรูปของส่วนหนึ่งของ $R_i$ เพราะบล็อกก้อนนี้เองคือเจ้าของ "สิบช่อง" ที่เรากำลังจะไปนับกัน สำหรับ $k = 4$ ได้ดังนี้

โครงสร้าง	บอกอะไรเกี่ยวกับการวัดสี่ครั้งของคนหนึ่งคน	จำนวน parameter ของ $V_i$ แบบ marginal, $k=4$
Independence, $V_i = \sigma^2 I_4$	ไม่มีความสัมพันธ์ภายในคนเลย ทุกแถวถูกมองว่าไม่เกี่ยวข้องกัน	1
Compound symmetry (exchangeable)	variance เท่ากันทุกตำแหน่ง และ correlation ค่าเดียวสำหรับทุกคู่ คือสิ่งที่ Spec A ผลิตออกมา	2
First-order autoregressive, AR(1)	variance ค่าเดียวร่วมกัน และ correlation ที่ลดลงตามระยะห่างระหว่างครั้งที่วัด	2
Heteroscedastic compound symmetry (Spec B)	variance หนึ่งค่าต่อหนึ่งหน่วยวัด และ covariance ร่วม $\tau^2$ ค่าเดียวสำหรับทุกคู่	5
Unstructured (corSymm + varIdent)	ทุก variance และทุก correlation รายคู่ อิสระต่อกันหมด	10

AR(1) มีข้อควรระวังเฉพาะตัวหนึ่งข้อ มันมีความหมายก็ต่อเมื่อมิติของการวัดซ้ำถูกเรียงลำดับจริง ๆ และระยะห่างระหว่างการวัดมีความหมายจริง ๆ ซึ่งเป็นเหตุผลที่มันไปอยู่กับ visit sub-study ใน ตอนที่ 5 ไม่ใช่ที่นี่ ตำแหน่งวัดสี่ตำแหน่งของตัวอย่างเดินเรื่องนี้เป็นเซตตายตัวที่ไม่มีลำดับ เป็นระดับเชิงนาม (nominal) ไม่ใช่ตำแหน่งบนสเกลใด ๆ ลำดับ forearm, hand, shin, back เป็นแค่ธรรมเนียมการติดป้าย การ fit โครงสร้างที่ถือว่า hand "อยู่ใกล้" shin มากกว่า back จึงเป็นการยัดเยียดมาตรวัดระยะที่ดีไซน์ไม่ได้มีอยู่จริง ในทางกลับกัน compound symmetry ไม่เปลี่ยนค่าไม่ว่าจะเรียงลำดับสี่หน่วยวัดแบบไหน

ถ้าอยากเห็นบันไดขยับจริง ๆ widget ICC explorer ที่สร้างไว้ใน ตอนที่ 2 มีตัวเลือก $k$ ให้ปรับ compound symmetry อยู่ที่ 2 parameter ตลอดไม่ว่าจะดัน $k$ ไปไกลแค่ไหน ในขณะที่จำนวน parameter ของ unstructured คือ $k(k+1)/2$ ซึ่งเท่ากับ 3 เมื่อ $k=2$ และเท่ากับ 21 เมื่อ $k=6$

เพดาน identifiability นับกันเป็นเหรียญ

ก่อนลุยของยาก: ชุดข้อมูลตอบคำถามได้มากที่สุดเท่ากับจำนวนคำตอบที่มันบรรจุอยู่จริง ๆ ฟังดูชัดเจนจนไม่น่าต้องพูด แต่เราเขียนโมเดลที่แอบถามมากกว่าที่ข้อมูลมีอยู่หนึ่งคำถามได้จริง ๆ และการเก็บคนไข้เพิ่มก็ไม่มีวันแก้ได้ เพราะคำตอบที่ขาดไปนั้นไม่เคยมีอยู่ให้เก็บตั้งแต่แรก ไม่กี่บรรทัดข้างล่างนี้จะนับจำนวนคำถามก่อน แล้วค่อยนับจำนวนคนถาม ถ้าตัวเลขที่สองมากกว่าตัวเลขแรก โมเดลนั้นก็ over-parameterised

ข้อโต้แย้งทั้งหมดอยู่ตรงนี้ และมันเป็นเลขคณิต ไม่ใช่ความเห็น covariance matrix ระดับบุคคลเป็นเมทริกซ์สมมาตรขนาด $4 \times 4$ จึงบรรจุปริมาณอิสระได้

\[ \frac{k(k+1)}{2} = \frac{4 \times 5}{2} = 10 \]

ปริมาณ นั่นคือ variance สี่ค่าบนเส้นทแยงมุม และ covariance หกค่าเหนือเส้นทแยงมุม มีสิบช่อง ไม่มีมากกว่านั้นแม้แต่ช่องเดียว ทีนี้มานับว่า unstructured ใช้ไปเท่าไร varIdent ซื้อ variance ไปสี่ค่า corSymm ซื้อ correlation ไปหกค่า รวมเป็นสิบ parameter สำหรับสิบปริมาณ งบประมาณถูกใช้จนหมดเกลี้ยงพอดี

\[ \underbrace{4}_{\texttt{varIdent}} \;+\; \underbrace{6}_{\texttt{corSymm}} \;+\; \underbrace{1}_{\tau^2} \;=\; 11 \quad \text{parameter สำหรับ} \quad \frac{k(k+1)}{2} = 10 \quad \text{ปริมาณ} \]

การเติม random intercept เข้าไปในโมเดลนั้นคือการขอให้เหรียญสิบเอ็ดเหรียญลงไปในช่องสิบช่อง เหรียญที่สิบเอ็ดไม่ได้โลภ แต่มันไม่มีบ้านอยู่ต่างหาก นี่คือความซ้ำซ้อน (redundancy) ในความหมายที่เคร่งครัด ไม่ใช่แค่พอดีตัวเกินไป และไม่ว่าจะเก็บข้อมูลเพิ่มอีกเท่าไรก็แก้ไม่ได้

สามวิธีแบ่ง เมทริกซ์เดียวกัน

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

สัญชาตญาณเบื้องหลังเพดานนี้ควรถูกเล่าอย่างช้า ๆ เพราะมันอธิบายว่าทำไมความล้มเหลวที่เกิดขึ้นจึงไม่ใช่อุบัติเหตุเชิงตัวเลข สมมติว่าเราปล่อยให้ residual covariance เป็นค่าใดก็ได้เรียกว่า $c$ แทนที่จะบังคับให้เป็นศูนย์ สิ่งที่ marginal matrix ของหนึ่งคนรายงานออกมามีสองอย่าง

\[ \operatorname{Cov}(Y_{ij}, Y_{ik}) = \tau^2 + c, \qquad \operatorname{Var}(Y_{ij}) = \tau^2 + \sigma_j^2 \]

ตรึงเป้าหมายไว้ที่เมทริกซ์ของ Spec B ซึ่งมี 6.25 ในทุกช่องนอกเส้นทแยงมุม และมี 8.00, 9.00, 20.00, 10.00 บนเส้นทแยงมุม คู่ $(\tau^2, c)$ ใดก็ตามที่ทำให้ $\tau^2 + c = 6.25$ เมื่อจับคู่กับ residual variance $\sigma_j^2 = V_{jj} - \tau^2$ จะสร้างเมทริกซ์นั้นกลับมาได้ครบทุกช่อง ข้อจำกัดมีเพียงข้อเดียว คือเมทริกซ์ residual ที่ตามมาต้องยังเป็น covariance matrix ที่ถูกต้องด้วย ซึ่งบีบให้ $\tau^2$ อยู่ในช่วงตั้งแต่ 0 ถึงประมาณ 7.03 ยกตัวอย่างสามชุดที่ล้วนอยู่ในช่วงนั้นอย่างสบาย ๆ

การแบ่ง	$\tau^2$	residual covariance $c$	residual variance (forearm, hand, shin, back)	สร้าง $V$ เดิมได้ไหม
A	6.25	0.00	1.75, 2.75, 13.75, 3.75	ได้
B	3.25	3.00	4.75, 5.75, 16.75, 6.75	ได้
C	0.00	6.25	8.00, 9.00, 20.00, 10.00	ได้

ลองตรวจการแบ่งแบบ B ด้วยมือ: $3.25 + 3.00 = 6.25$ ในทุกช่องนอกเส้นทแยงมุม และ $3.25 + 4.75 = 8.00$, $3.25 + 5.75 = 9.00$, $3.25 + 16.75 = 20.00$, $3.25 + 6.75 = 10.00$ บนเส้นทแยงมุม เมทริกซ์ที่ fit ได้เหมือนของการแบ่งแบบ A ทุกประการ และของแบบ C ก็เหมือนกัน ทั้ง ๆ ที่แบบ C ไม่มี person effect อยู่เลย เพราะความสัมพันธ์ทั้งหมดถูกดูดเข้าไปอยู่ใน residual covariance หมดแล้ว

เมื่อทั้งสาม specification ผลิต $V$ ตัวเดียวกัน มันก็ผลิต multivariate normal density ตัวเดียวกัน และผลิตค่า likelihood เท่ากันที่ทุกจุดข้อมูล likelihood จึงแบนราบไปตามทิศทางนั้น ข้อมูลไม่สามารถเลือกได้ว่าการแบ่งแบบไหนถูก ไม่ใช่เพราะตัวอย่างน้อย แต่เพราะคำถามนี้ไม่มีคำตอบอยู่ในข้อมูลตั้งแต่ต้น สิ่งที่ข้อมูลตรึงไว้ได้ และตรึงได้แม่นยำขึ้นเรื่อย ๆ เมื่อขนาดตัวอย่างโตขึ้น คือผลรวม $\tau^2 + c = 6.25$ ส่วนคำถามที่ว่าผลรวมนั้นจะถูกแบ่งตรงจุดไหนภายในช่วงที่ยอมรับได้ ไม่ใช่คำถามที่การเก็บข้อมูลเพิ่มจะตอบให้ได้

หน้าตาของ over-parameterisation ใน output

โปรแกรมจะไม่พิมพ์ออกมาว่า "โมเดลนี้ระบุค่าไม่ได้" มันพิมพ์อาการออกมาแทน อาการส่วนใหญ่เป็นอาการร่วมของ specification ที่ over-parameterised ไม่ว่าจะออกมาจากแพ็กเกจไหนก็ตาม และควรจำให้ได้ตั้งแต่แรกเห็น ประการแรก การ fit อาจไม่ converge หรือ converge ก็ต่อเมื่อเพิ่มจำนวนรอบการวนซ้ำ พร้อมคำเตือนเกี่ยวกับ covariance matrix หรือเกี่ยวกับ step ของ optimiser ประการที่สอง variance component อาจไปนั่งอยู่บนขอบเขตพอดี เช่น $\widehat{\tau}^{2}$ ยุบลงเป็นศูนย์ หรือ correlation ของ random effects ถูกประมาณได้เท่ากับ $\pm 1$ พอดี ซึ่ง lme4 จะรายงานว่าเป็น singular fit ประการที่สาม การ fit ใหม่ด้วย starting value คนละชุดอาจทำให้ variance component ขยับไปมากพอสมควร ในขณะที่ค่าประมาณของ fixed effect และ covariance matrix ที่ fit ได้แทบไม่เปลี่ยน ซึ่งเป็นการที่ likelihood แบน ๆ เผยตัวออกมาตรง ๆ สุดท้าย standard error ของ variance component อาจใหญ่มหาศาล หรืออาจคำนวณไม่ออกเลย

อาการเหล่านี้ไม่ควรถูกอ่านว่าเป็นหลักฐานว่าความแปรปรวนระหว่างบุคคลไม่มีอยู่จริง ค่าประมาณที่ขอบเขตอย่าง $\widehat{\tau}^{2} = 0$ ในโมเดลที่มี residual covariance อิสระอยู่ด้วย คือเลขคณิตของการแบ่งแบบ C ไม่ใช่ข้อค้นพบเกี่ยวกับตัวบุคคล สิ่งที่ควรทำคือถอดความซ้ำซ้อนออกจาก specification ไม่ใช่สรุปว่า person effect ไม่เคยมีอยู่

ถ้า gls ไม่มี person effect แล้วมันเมินการวัดซ้ำหรือเปล่า

ไม่ นี่คือความเข้าใจผิดที่พบบ่อยที่สุดเกี่ยวกับ marginal specification และควรตอบให้ชัดถ้อยชัดคำ การ fit ด้วย gls พร้อม corSymm และ varIdent คือการระบุ covariance matrix ภายในคนออกมาตรง ๆ แทนที่จะอนุมานมันจากปริมาณแฝง

# 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)

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

สิ่งที่เสียไปเป็นเรื่องจริงและควรพูดอย่างซื่อสัตย์ คุณเสีย $\tau^2$ ในฐานะ estimand ที่ตีความได้ และเสีย ICC ในความหมายของการแบ่งสัดส่วน variance ที่ ตอนที่ 2 สร้างไว้ทั้งตอน ส่วน correlation แบบ marginal ทั้งหกค่ายังถูกประมาณอยู่ และประมาณอย่างอิสระด้วย เพียงแต่ correlation ที่อ่านออกมาจากเมทริกซ์ตรง ๆ เป็นการบรรยายความสัมพันธ์ ไม่ใช่การแบ่งสัดส่วนของความสัมพันธ์นั้นระหว่างตัวบุคคลกับครั้งที่วัด นอกจากนี้คุณเสีย BLUPs ไปด้วย จึงไม่มีการทำนายระดับบุคคล และไม่มีทางถามว่าผู้เข้าร่วมคนไหนอยู่เหนือโปรไฟล์ของกลุ่มตัวเอง ส่วนสิ่งที่ยังเหลืออยู่ก็เป็นเรื่องจริงเช่นกัน ได้แก่ ค่าประมาณ fixed effect ที่ยัง consistent, standard error ที่เคารพความสัมพันธ์ภายในคน และ contrast ที่ถูกต้องภายใต้ covariance ที่ระบุไว้ การแลกเปลี่ยนนี้จะยอมรับได้หรือไม่ ขึ้นอยู่กับว่าคำถามวิจัยของคุณเป็นคำถามเกี่ยวกับ "คน" หรือเกี่ยวกับ "ประชากร" ซึ่งเป็นเนื้อหาทั้งหมดของตอนถัดไป

เลือก covariance model ด้วยหลักฐาน ไม่ใช่ด้วยความสวยงาม

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

บันไดนี้ชวนให้เกิดนิสัยเสียอย่างหนึ่ง คือปีนขึ้นไปเรื่อย ๆ จนกว่า output จะดูน่าประทับใจ ขั้นตอนการเลือกจึงควรถูกเขียนเป็นโปรโตคอลสั้น ๆ มากกว่าจะเป็นรสนิยม ประการแรก พล็อต residual เทียบกับ fitted value แล้วดูการกระจายแยกตามหน่วยวัดแต่ละหน่วย ถ้าเห็นการบานออกอย่างมีโครงสร้าง หรือเห็นแถบที่หนากว่าเพื่อนอย่างชัดเจนที่ตำแหน่งใดตำแหน่งหนึ่ง นั่นคือ diagnostic ที่เป็นเหตุผลให้ใช้ residual แบบ heteroscedastic ประการที่สอง เปรียบเทียบการกระจายของ residual รายหน่วยเป็นตัวเลขด้วย เช่นสรุป standardised residual รายตำแหน่ง เพื่อให้สิ่งที่ตาเห็นจากกราฟถูกวัดออกมาเป็นค่า ประการที่สาม เปรียบเทียบโครงสร้าง variance ที่ซ้อนกัน (nested) ด้วย likelihood-ratio test โดยตรึง fixed effect ให้เหมือนกันเป๊ะและ fit ทั้งสองโมเดลด้วย REML เพราะ REML likelihood เปรียบเทียบข้าม mean structure ที่ต่างกันไม่ได้ ประการที่สี่ ใช้ AIC และ BIC เป็นหลักฐานสนับสนุน ไม่ใช่ผู้ตัดสิน และพึงผูกข้อควรระวังเรื่องขอบเขตไว้กับการเปรียบเทียบที่มันเป็นเจ้าของจริง ๆ กล่าวคือ การทดสอบว่า variance ของ random effect เท่ากับศูนย์หรือไม่นั้นอยู่บนขอบของปริภูมิพารามิเตอร์และเป็นการทดสอบแบบ conservative ในขณะที่การเปรียบเทียบระหว่างแบบ pooled กับ varIdent เป็นการทดสอบอัตราส่วนของ variance ซึ่งค่าภายใต้สมมติฐานหลักอยู่ภายในปริภูมิ จึงใช้การแจกแจงอ้างอิงตามปกติได้ สุดท้าย ควรระบุการเปรียบเทียบนี้ไว้ล่วงหน้า แทนที่จะรันหลังจากเห็นผลของ fixed effect แล้ว

มีสองอย่างที่ไม่ควรอยู่ในโปรโตคอลนี้ อย่างแรก การเลือก covariance structure โดยดูจาก p-value ของ fixed effect ที่มันผลิตออกมา ไม่ใช่การตรวจสอบโมเดล แต่คือการเลือกคำตอบก่อนแล้วค่อยเลือกวิธี อย่างที่สอง ความเรียบร้อยสวยงามของ output ไม่ใช่หลักฐานของอะไรทั้งสิ้น กฎที่ ตอนที่ 8 จะยืนอยู่บนนั้นควรถูกรับมาใช้ตั้งแต่ตอนนี้ เปลี่ยน covariance model เพราะ diagnostic ขัดแย้งกับสมมติฐาน ไม่ใช่เพราะ output ดูเรียบร้อยเกินไป

สิ่งที่โพสต์นี้ไม่ได้อนุญาตให้สรุป

ข้อโต้แย้งข้างต้นถูกอ่านเกินเลยได้ง่ายอยู่สามทิศทาง และควรปิดประตูทั้งสามทางไว้ให้ชัด

ทางแรก นี่ไม่ใช่เหตุผลให้ fit โครงสร้างที่รวยที่สุดเป็นค่าตั้งต้น parameter สิบตัวอาจจะพอไหวกับผู้เข้าร่วมจำลอง 80 คนที่วัดคนละสี่ครั้ง อย่างไรก็ตาม specification แบบ unstructured เดียวกันที่ $k = 10$ ครั้งจะขอไปถึง 55 ตัว และเมื่อจำนวน cluster ไม่มากนัก ค่าประมาณอาจไม่เสถียรและ standard error อาจแคบเกินจริง ความประหยัด parameter ไม่ใช่ความขี้ขลาด มันอาจเป็นเส้นแบ่งระหว่างโมเดลที่ประมาณค่าได้กับโมเดลที่เป็นแค่เครื่องประดับ

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

ทางที่สาม การอนุมานย้อนกลับก็ทำไม่ได้เช่นกัน โมเดลที่ fit ออกมาแล้วรายงาน correlation เท่ากันทุกคู่ ไม่ได้พิสูจน์ว่า correlation เท่ากันจริง มันแค่ บังคับ ให้เท่ากัน ความเท่ากันใน output ของ compound-symmetry fit คือการพูดซ้ำสมมติฐานที่เรายื่นให้โปรแกรมไป และตรวจสอบได้ทางเดียวคือผ่อนคลายสมมติฐานนั้นแล้วดูว่าอะไรเปลี่ยนไปบ้าง

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

1. เขียนสมมติฐานของ residual covariance ลงในเมธอดให้ชัด ในประโยคเดียวกับที่เขียนโครงสร้าง random effects ประโยคว่า "random intercept ระดับผู้เข้าร่วม ร่วมกับ residual ที่เป็นอิสระต่อกันและ homoscedastic" คือ specification ที่ครบถ้วน ส่วน "ใช้ mixed model เพื่อจัดการกับการวัดซ้ำ" ยังไม่ครบ
2. ดูการกระจายของ residual รายหน่วยวัดก่อนเลือกโครงสร้าง และให้ diagnostic เป็นตัวขับการตัดสินใจ แทนที่จะให้สถิติความกลมกลืนของโมเดลเป็นตัวขับ
3. ถ้า diagnostic ขัดแย้งกับการ pool ก็ควร fit เวอร์ชัน heteroscedastic แล้วเปรียบเทียบบน fixed effect ชุดเดียวกัน และอาจรายงานการเปรียบเทียบนั้นเป็น sensitivity analysis ที่ระบุไว้ล่วงหน้า มากกว่าจะยกให้เป็นโมเดลหลัก
4. อย่าวาง random intercept ทับลงบน unstructured residual covariance ให้เลือกอย่างใดอย่างหนึ่ง คืออธิบายความสัมพันธ์ด้วย person effect แฝง หรือระบุความสัมพันธ์นั้นออกมาตรง ๆ เพราะการใช้ทั้งคู่คือการขอ parameter ที่ covariance matrix ไม่มีที่ว่างให้เก็บ
5. ให้ estimand เป็นตัวตัดสินเมื่อสอง specification สะดวกไม่เท่ากัน ถ้าคำถามต้องการ $\tau^2$, ICC, การแบ่งสัดส่วน variance หรือการทำนายระดับบุคคล รูปแบบ random effects อาจเป็นทางเดียวที่ตอบได้ ไม่ว่าทางเลือกแบบ marginal จะดูยืดหยุ่นกว่าเพียงใด
6. รายงานว่าอะไรเปลี่ยน ไม่ใช่รายงานว่าไม่มีอะไรเปลี่ยน ถ้า sensitivity analysis ทำให้ช่วงความเชื่อมั่นหรือ contrast รองเปลี่ยนไป นั่นอาจคุ้มค่าที่จะเขียนเป็นประโยคของมันเองในผลการศึกษา

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

• random intercept สมทบค่า $\tau^2 = 6.25$ เท่ากันเข้าไปใน covariance ของการวัดซ้ำทุกคู่ในคนเดียวกัน และเมื่อ residual ไม่มี correlation ต่อกัน ค่าที่สมทบเข้าไปนั้นก็คือ covariance ทั้งก้อน แต่สิ่งที่มันไม่เคยบังคับเลยคือการให้แชร์ correlation ค่าเดียวกัน
• compound symmetry ต้องการสมมติฐานสองข้อ คือ offset ร่วมหนึ่งค่า และ residual ที่เป็นอิสระต่อกันโดยมี variance เท่ากัน และข้อที่สองคือครึ่งที่มักถูกละไว้ไม่พูดถึง
• ใน Spec B ค่า covariance เท่ากับ 6.25 ในทุกช่องนอกเส้นทแยงมุม ในขณะที่ correlation ทั้งหกค่าไล่ตั้งแต่ 0.44 ถึง 0.74 เพราะแต่ละค่าถูกหารด้วยคู่ของ marginal standard deviation ที่ต่างกัน
• covariance matrix สมมาตรขนาด $4 \times 4$ บรรจุปริมาณอิสระได้ $k(k+1)/2 = 10$ ค่า corSymm บวก varIdent ใช้ไปครบ 10 พอดี และการเติม $\tau^2$ ทำให้กลายเป็น 11 parameter สำหรับ 10 ปริมาณ ซึ่งเป็นความซ้ำซ้อนที่ขนาดตัวอย่างเท่าไรก็แก้ไม่ได้
• เนื่องจากการแบ่ง $(\tau^2, c)$ ได้หลายแบบให้เมทริกซ์ที่ fit ออกมาเหมือนกันเป๊ะ ในที่นี้คือทุกค่า $\tau^2$ ตั้งแต่ 0 ถึงประมาณ 7.03 likelihood จึงแบนไปตามทิศทางนั้น และคำเตือนเรื่อง convergence, ค่าประมาณที่ขอบเขต และความไวต่อ starting value คือหน้าตาของความแบนนั้นเมื่อมองจากภายนอก
• การ fit ด้วย gls ที่ระบุ covariance matrix ออกมาตรง ๆ นับรวมการวัดซ้ำไว้ครบถ้วน สิ่งที่มันยอมสละคือคำอธิบายว่าความสัมพันธ์มาจากไหน พร้อมกับ ICC ในความหมายของการแบ่งสัดส่วน variance และการทำนายระดับบุคคล

ตอนต่อไป

ตอนที่ 7 จะเดินถนนสาย marginal ไปจนสุด ทั้ง GLS และ GEE ที่ปฏิเสธจะอธิบาย correlation เลยแม้แต่น้อย แล้วแสดงให้เห็นว่าสำหรับ outcome ต่อเนื่องราคาที่จ่ายแทบเป็นศูนย์ ในขณะที่สำหรับ outcome แบบ binary มันเปลี่ยนค่า odds ratio เลยทีเดียว

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

1. Verbeke G, Molenberghs G. Linear Mixed Models for Longitudinal Data. New York: Springer; 2000.
2. Diggle PJ, Heagerty P, Liang KY, Zeger SL. Analysis of Longitudinal Data. 2nd ed. Oxford: Oxford University Press; 2002.
3. Pinheiro JC, Bates DM. Mixed-Effects Models in S and S-PLUS. New York: Springer; 2000.