ค่าชดเชยส่วนตัวที่ถูกแชร์: Random Intercept คืออะไรกันแน่

บทคัดย่อ

Random intercept ไม่ใช่จุดตัดแกน Y ที่สุ่มมา แต่คือค่าเดียวต่อหนึ่งคน ที่ถูกใส่ซ้ำไม่เปลี่ยนแปลงในทุกแถวของคนคนนั้น บทความนี้ใช้ข้อมูลจำลองการวัด TEWL ที่ผิวหนังสี่ตำแหน่งในผู้ใหญ่ 80 คน พาไปเห็นค่าชดเชยประจำตัวในข้อมูลดิบก่อน fit โมเดล แล้วพิสูจน์ในสามบรรทัดว่าการแชร์ตัวเลขค่าเดียวกันเท่ากับการทำให้เกิด covariance นั่นคือการวัดสองครั้งจากคนเดียวกันมี covariance เท่ากับ between-person variance 6.25 จาก variance ทั้งหมด 10.00 ของหนึ่งแถว ผู้อ่านจะเข้าใจว่า correlation คือผลของโครงสร้างข้อมูล ไม่ใช่ตัวเลือกในโปรแกรม ได้อ่าน random = ~ 1 | id ทีละ token และเห็นว่าทำไม offset จึงหักล้างกันใน within-person contrast จนได้ SE 0.43 เทียบกับ 0.50

---

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

ซีรีส์นี้ติดตามปริมาณเพียงตัวเดียวตั้งแต่สเปรดชีตแผ่นแรกจนถึงคอลัมน์ standard error สุดท้าย นั่นคือ $b_i$ หรือ "ตัวเลขที่คนหนึ่งคนพกติดตัวเข้าไปในทุกการวัดที่เขามีส่วนร่วม" โพสต์แรกนี้จะพาไปดูตัวเลขนั้นในข้อมูลดิบก่อนที่จะ fit โมเดลใด ๆ แล้วจึงพิสูจน์ด้วยพีชคณิตสามบรรทัดว่า การแชร์ตัวเลขเดียวกันนั้นมีความหมายเท่ากับการทำให้เกิด correlation ขึ้นมาเอง โพสต์ถัดไปจะขยาย covariance ค่าเดียวที่ได้จากที่นี่ให้กลายเป็น matrix ขนาด 4×4 เต็มรูป — Compound Symmetry กับ ICC: ออฟเซ็ตเดียว correlation เดียว

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

ชุดข้อมูล และคำถามหนึ่งข้อก่อนจะเริ่มสร้างโมเดล

การศึกษาจำลองที่อยู่เบื้องหลังซีรีส์นี้ตั้งใจให้ธรรมดาที่สุด ผู้เข้าร่วม 80 คน — 40 คนเป็น moderate atopic dermatitis (AD) และอีก 40 คนไม่มีโรคภูมิแพ้ใด ๆ (control) — มาพบแพทย์เพียงครั้งเดียว และวัด transepidermal water loss (TEWL, g·m⁻²·h⁻¹) ด้วยเครื่อง closed-chamber evaporimeter ที่ผิวหนังซึ่งไม่มีรอยโรค 4 ตำแหน่งมาตรฐาน โดยวัดเรียงลำดับเดิม ภายใต้อุณหภูมิและความชื้นห้องเดียวกัน ค่า TEWL ที่สูงขึ้นหมายถึง skin barrier ที่รั่วมากขึ้น

องค์ประกอบ	ค่า
ผู้เข้าร่วม	80 คน (AD 40 คน, control 40 คน)
หน่วยที่วัดซ้ำ	site, $k = 4$ ตำแหน่งต่อคน
ชื่อหน่วย	forearm · hand · shin · back (reference = forearm)
จำนวน observation	320 แถว สมดุลเต็มรูปแบบ
ตัวแปรกลุ่ม	group ∈ {AD, control}, reference = control
Outcome	tewl, ต่อเนื่อง, g·m⁻²·h⁻¹
ตัวแปรระดับบุคคล	age (ปี), sex, phototype (Fitzpatrick I–VI)

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

id	group	forearm	hand	shin	back
S01	control	14.1	20.1	11.4	17.0
S02	control	9.7	10.9	8.3	10.4
S03	AD	17.6	27.2	14.2	18.9
S04	AD	16.9	20.6	12.7	17.3

ลองใช้เวลาสักสามสิบวินาทีกับตารางนี้ก่อนอ่านต่อ เปรียบเทียบ S01 กับ S02 ทั้งคู่เป็น control ทั้งคู่ถูกวัดที่ตำแหน่งเดียวกันสี่ตำแหน่งด้วยเครื่องเดียวกัน แต่ S01 สูงกว่า S02 ทุกตำแหน่ง ทั้งที่ forearm ที่ hand ที่ shin และที่ back ระยะห่างไม่ได้เท่ากันทั้งสี่ครั้ง แต่ไม่เคยเปลี่ยนเครื่องหมายเลย

เฉลย: ทุกคนกำลังถูกเทียบกับตัวเลขผิดตัว

จริง ๆ แล้วสายตาของเราไม่ได้เปรียบเทียบ S01 กับ S02 แต่กำลังเปรียบเทียบแต่ละคนกับ "โปรไฟล์ที่กลุ่มของเขาควรจะเป็น" และเพราะชุดข้อมูลนี้เป็นข้อมูลจำลอง เราจึงรู้โปรไฟล์นั้นแบบเป๊ะ ๆ เนื่องจากเราเขียนมันขึ้นมาก่อนจะสร้างข้อมูลแม้แต่แถวเดียว

site	control	AD	ผลต่าง (AD − control)
forearm	12.0	18.0	+6.0
hand	16.0	24.0	+8.0
shin	10.0	13.0	+3.0
back	14.0	19.0	+5.0

ผลต่างไม่เท่ากันในแต่ละตำแหน่ง แปลว่า interaction ระหว่าง site กับ group มีอยู่จริงในข้อมูลจำลองชุดนี้ และต้องอยู่ในโมเดลทุกตัวที่เรา fit ต่อจากนี้ ทีนี้ลองเอาค่าที่วัดได้ทั้งสี่ค่าของแต่ละคน ลบด้วยโปรไฟล์ของกลุ่มตัวเอง

id	group	forearm	hand	shin	back
S01	control	+2.1	+4.1	+1.4	+3.0
S02	control	−2.3	−5.1	−1.7	−3.6
S03	AD	−0.4	+3.2	+1.2	−0.1
S04	AD	−1.1	−3.4	−0.3	−1.7

S01 อยู่เหนือโปรไฟล์ของกลุ่มตัวเองทั้งสี่ตำแหน่ง S02 อยู่ต่ำกว่าทั้งสี่ตำแหน่ง และ S04 ก็อยู่ต่ำกว่าทั้งสี่ตำแหน่งเช่นกัน นี่คือรูปแบบที่เห็น และคำอธิบายแบบ measurement error ที่เป็นอิสระต่อกันเข้ากับรูปแบบนี้ได้ไม่ค่อยสนิทนัก ที่น่าสังเกตคือ S03 เป็นตัวอย่างแย้งที่ซื่อสัตย์ — บวกสอง ลบสอง — ซึ่งเราจะกลับมาอธิบายในอีกสักครู่ เพราะมันคือสิ่งที่กันไม่ให้เราเชื่อการ "ดูด้วยตา" มากเกินไป

ลืมคำว่า correlation ไปก่อน

ลองคิดว่า ordinary linear regression ที่ fit บนทั้ง 320 แถว กำลังอ้างอะไรอยู่ มันบอกว่า control ทุกคนที่วัดที่ forearm มีค่าคาดหมายเท่ากันหมดคือ 12.0 บอกว่า AD ทุกคนที่วัดที่ hand มีค่าคาดหมายเท่ากันหมดคือ 24.0 และบอกว่าสิ่งที่เหลือในแต่ละแถวคือ $e_{ij}$ ซึ่งเป็น error ที่ถูกสุ่มขึ้นมาใหม่และเป็นอิสระต่อกันในทั้ง 320 แถว ภายใต้โมเดลแบบนั้น สองแถวที่มาจาก S01 มี "อะไรร่วมกัน" ไม่มากไปกว่าแถวหนึ่งจาก S01 กับอีกแถวหนึ่งจาก S02 เลย

ตารางข้างบนเข้ากับคำอ้างนี้ได้ไม่ค่อยสนิท ค่าเบี่ยงเบนทั้งสี่ของ S01 คือ +2.1, +4.1, +1.4 และ +3.0 และกลไกที่สุ่ม error สี่ค่าอย่างเป็นอิสระจากการแจกแจงที่สมมาตร จะต้องอาศัยความบังเอิญจึงจะได้ค่าทั้งสี่ตกอยู่ฝั่งเดียวกันของศูนย์ ซึ่งเกิดขึ้นด้วยความน่าจะเป็น $2 \times (1/2)^4 = 1/8$ สำหรับคนหนึ่งคน การเจอคนแบบนี้หนึ่งคนจึงไม่ใช่เรื่องน่าตื่นเต้น แต่ตารางที่คนส่วนใหญ่เป็นแบบนี้เป็นอีกเรื่องหนึ่ง ดังนั้นคำถามจึงไม่ใช่ว่าควรเผื่อที่ให้ "อะไรบางอย่างที่เป็นของคนนั้น" หรือไม่ แต่คือจะเติมอะไรเข้าไป และมีต้นทุนอะไรตามมา

ทำให้ offset ปรากฏชัด

ก่อนลุยของยาก: สิ่งที่ตามมาข้างล่างคือความคิดทางคลินิกอย่างเดียวที่ถูกเขียนเป็นสัญลักษณ์ นั่นคือผู้ป่วยบางคน "ค่าสูงโดยพื้นฐาน" และเขาจะสูงในทุกครั้งที่เราวัดเขา เราเจอเรื่องนี้บนวอร์ดอยู่แล้ว เช่นผู้ป่วยที่ความดันสูงกว่าคนอื่นเล็กน้อยทั้งรอบเช้า รอบบ่าย และรอบดึก ไม่ใช่เพราะวัดผิดไปทางเดียวกันสามครั้ง แต่เพราะคนคนนั้น "ตั้งอยู่" ตรงนั้นของเขาเอง สมการข้างล่างเพียงตั้งชื่อให้ "จุดที่ผู้ป่วยตั้งอยู่" ว่า $b_i$ และตั้งชื่อให้ความแกว่งของการวัดครั้งใดครั้งหนึ่งว่า $e_{ij}$ เท่านั้น ไม่ได้อ้างอะไรมากไปกว่านี้

เติมเทอมเดียวเข้าไปในโมเดล — ตัวเลขหนึ่งค่าที่ผูกกับ "คน" ไม่ได้ผูกกับ "แถว"

\[ Y_{ij} = \mu_{j,\,g(i)} + b_i + e_{ij}, \qquad b_i \sim N(0,\ \tau^2), \qquad e_{ij} \sim N(0,\ \sigma^2), \qquad b_i \perp e_{ij}, \quad e_{ij} \perp e_{ik}\ (j \neq k) \]

อ่านสามก้อนนี้เป็นสามหน้าที่ที่แยกกัน $\mu_{j,g(i)}$ คือสัญญาณระดับประชากรของตำแหน่ง $j$ ในกลุ่มที่คนที่ $i$ สังกัดอยู่ ส่วน $b_i$ คือระยะที่คนคนนี้อยู่ห่างจากสัญญาณนั้น และ $e_{ij}$ คือสิ่งที่เหลือ ณ ตำแหน่งนั้นโดยเฉพาะ ข้อความเรื่องความเป็นอิสระสองข้อที่ต่อท้ายบรรทัดไม่ใช่ของประดับ — offset เป็นอิสระจาก noise และ noise ที่ตำแหน่งหนึ่งเป็นอิสระจาก noise ที่อีกตำแหน่งหนึ่ง — และพีชคณิตข้างล่างจะใช้ทั้งสองข้อนี้จนหมด ในข้อมูลจำลองชุดนี้ $\tau^2 = 6.25$ (SD 2.5) และ $\sigma^2 = 3.75$ (SD 1.94) และสาระทั้งหมดของโพสต์นี้คือข้อสังเกตเดียวว่า $b_i$ ไม่มี $j$ อยู่ในตัวมันเลย

เมื่อเขียนกระจายออกมา ชุดข้อมูลของผู้เข้าร่วมทั้งสี่คนคือผลบวกของสามคอลัมน์

id	site	ค่าเฉลี่ยประชากร $\mu_{j,g}$	+ offset ของคน $b_i$	+ noise $e_{ij}$	= ค่าที่วัดได้
S01	forearm	12.0	+2.5	−0.4	14.1
S01	hand	16.0	+2.5	+1.6	20.1
S01	shin	10.0	+2.5	−1.1	11.4
S01	back	14.0	+2.5	+0.5	17.0
S02	forearm	12.0	−3.0	+0.7	9.7
S02	hand	16.0	−3.0	−2.1	10.9
S02	shin	10.0	−3.0	+1.3	8.3
S02	back	14.0	−3.0	−0.6	10.4
S03	forearm	18.0	+0.8	−1.2	17.6
S03	hand	24.0	+0.8	+2.4	27.2
S03	shin	13.0	+0.8	+0.4	14.2
S03	back	19.0	+0.8	−0.9	18.9
S04	forearm	18.0	−2.0	+0.9	16.9
S04	hand	24.0	−2.0	−1.4	20.6
S04	shin	13.0	−2.0	+1.7	12.7
S04	back	19.0	−2.0	+0.3	17.3

ภาพที่ 1. กายวิภาคของหนึ่งแถวข้อมูล — หนึ่ง observation ประกอบด้วยสามส่วนที่บวกกัน: สัญญาณระดับประชากร 12.0, offset ของคน +2.5 และ measurement noise −0.4 รวมเป็นค่าที่วัดได้ 14.1 ของ S01 ที่ forearm

คอลัมน์ offset ซ้ำเดิมทุกแถว ส่วนคอลัมน์ noise ไม่เคยซ้ำเลย S01 พา +2.5 เข้าไปในสี่แถว S02 พา −3.0 เข้าไปในสี่แถว และ noise เปลี่ยนค่าทุกแถวจริง ๆ ซึ่งเป็นเหตุผลพอดีว่าทำไม S03 ที่มี offset เพียง +0.8 เมื่อเทียบกับ residual SD 1.94 จึงดูเหมือนคนที่มี error เป็นอิสระธรรมดา ทั้งที่ถูกสร้างจากกลไกเดียวกันเป๊ะ นี่คือคำเตือนแรกของโพสต์นี้ offset เป็นคุณสมบัติของโมเดล แต่การที่มัน "มองเห็นได้" ในผู้เข้าร่วมคนใดคนหนึ่งเป็นเรื่องของโชค

ทำไมจึงเรียกว่า intercept

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

ถ้าทุกคนใช้โปรไฟล์ของตำแหน่งชุดเดียวกัน และต่างกันแค่ $b_i$ เส้นโปรไฟล์ ค่าเฉลี่ย ของแต่ละคนก็คือเส้นโปรไฟล์ประชากรที่ถูกเลื่อนขึ้นลงในแนวดิ่ง ขนานกันเสมอ ไม่เอียง ไม่ตัดกัน เส้นค่าเฉลี่ยของ S01 อยู่สูงกว่าโปรไฟล์ control อยู่ 2.5 หน่วยทั้งสี่ตำแหน่ง ส่วนของ S02 อยู่ต่ำกว่า 3.0 หน่วยทั้งสี่ตำแหน่ง ทั้งนี้ค่าที่วัดได้จริงไม่ได้วางอยู่บนเส้นเหล่านั้น และโมเดลก็ไม่ได้เรียกร้องให้เป็นเช่นนั้นด้วย ค่าเบี่ยงเบนที่เกิดขึ้นจริงของ S01 คือ +2.1, +4.1, +1.4 และ +3.0 ซึ่งก็คือเส้นขนานหนึ่งเส้นบวกกับ $e_{ij}$ ที่สุ่มมาสี่ครั้งแยกกัน ความขนานจึงเป็นคุณสมบัติของค่าเฉลี่ยแบบมีเงื่อนไข ไม่ใช่คุณสมบัติของข้อมูลที่วัดได้

ภาพที่ 2. สองโปรไฟล์ของค่าเฉลี่ยแบบมีเงื่อนไข — เส้นค่าเฉลี่ยประชากรทั้งสี่ตำแหน่งเป็นสีเทา เส้นค่าเฉลี่ยของ S01 ขนานและอยู่เหนือกว่า (+2.5) เส้นของ S02 ขนานและอยู่ต่ำกว่า (−3.0) ระยะห่างแนวดิ่งคงที่โดยการสร้างข้อมูล ส่วนค่าที่วัดได้จริงกระจายรอบแต่ละเส้นด้วย $e_{ij}$

ชื่อ "intercept" จึงพูดถึง ตำแหน่งที่เส้นวางอยู่ ไม่ได้พูดถึงว่าแกน Y เป็นจุดที่น่าสนใจ ส่วนคำว่า "random" ก็เข้าใจผิดกันง่ายพอ ๆ กัน มันไม่ได้แปลว่ามั่ว และไม่ได้แปลว่าโปรแกรมเป็นคนเลือกตัวเลขให้ แต่แปลว่า $b_i$ ถูกมองเป็นค่าที่สุ่มมาจากการแจกแจงหนึ่ง — ในที่นี้คือ $N(0,\ 6.25)$ — แทนที่จะเป็นพารามิเตอร์แยกอีก 80 ตัวที่ต้องประมาณค่าให้ผู้เข้าร่วมทีละคน และเพราะการแจกแจงนั้นมีศูนย์กลางที่ศูนย์ การเติม $b_i$ จึงไม่ไปแตะค่าเฉลี่ยระดับประชากรเลย สิ่งที่มันเปลี่ยนคืออย่างอื่น ซึ่งเป็นเนื้อหาของสองหัวข้อถัดไป

จังหวะที่แถวข้อมูลเลิกเป็นอิสระต่อกัน

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

เขียนสมการทั้งสี่ของ S01 ออกมา แล้วมองที่คอลัมน์กลาง

\[ \begin{aligned} Y_{\text{S01},\,\text{forearm}} &= 12.0 + 2.5 - 0.4 = 14.1 \\ Y_{\text{S01},\,\text{hand}} &= 16.0 + 2.5 + 1.6 = 20.1 \\ Y_{\text{S01},\,\text{shin}} &= 10.0 + 2.5 - 1.1 = 11.4 \\ Y_{\text{S01},\,\text{back}} &= 14.0 + 2.5 + 0.5 = 17.0 \end{aligned} \]

ค่า $+2.5$ ตัวเดียวกันปรากฏครบทั้งสี่แถวโดยไม่เปลี่ยนแปลง ไม่มีใครสั่งให้โมเดลสร้าง correlation ไม่มีพารามิเตอร์ correlation ถูกพิมพ์ลงไปที่ไหน และไม่มีการเลือก working correlation structure ใด ๆ ทั้งสิ้น แต่ถึงอย่างนั้น สี่แถวนี้ก็มีของบางอย่างร่วมกัน ซึ่งสี่แถวจากคนสี่คนไม่มี

ภาพที่ 3. พัด — โหนด $b_i$ หนึ่งจุดอยู่ด้านบน มีลูกศรสี่เส้นชี้ลงไปที่ forearm, hand, shin และ back พ่อแม่หนึ่ง ลูกสี่ ดังนั้นพี่น้องจึงสัมพันธ์กัน

พีชคณิตสามบรรทัด

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

หยิบสองตำแหน่งที่ต่างกัน $j \neq k$ จากคนเดียวกัน แล้วกระจาย covariance ออกมา

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

สามในสี่เทอมหายไปด้วยข้อสมมติเรื่องความเป็นอิสระสองข้อที่เราเขียนไว้ในโมเดลเอง ดังนั้น covariance ระหว่างการวัดสองครั้งจากคนเดียวกันจึงเท่ากับ variance ของสิ่งที่ทั้งสองแถวแชร์กันพอดี บรรทัดนี้ไม่ได้อ้างการแจกแจงใด ๆ และไม่ได้ประมาณอะไรเลย มันเป็น identity แต่เป็น identity ที่ใช้ข้อสมมติทั้งสองข้อ ไม่ใช่ข้อเดียว ถ้า residual มี covariance ที่ไม่เป็นศูนย์ ค่านั้นจะถูก บวกเพิ่ม เข้าไปกับ $\tau^2$ ไม่ใช่มาแทนที่ ซึ่งเป็นสถานการณ์ที่โพสต์ที่ 6 จะรื้อออกมาดูทีละชิ้น และขอให้สังเกตด้วยว่าคำตอบไม่มีทั้ง $j$ และ $k$ อยู่ในตัวมัน ซึ่งเป็นข้อเท็จจริงที่โพสต์ที่ 2 จะใช้เวลาทั้งโพสต์อยู่กับมัน

Variance ของ observation เดียวก็ได้มาจากการกระจายชุดเดียวกัน

\[ \operatorname{Var}(Y_{ij}) = \operatorname{Var}(b_i) + \operatorname{Var}(e_{ij}) = \tau^2 + \sigma^2 = 6.25 + 3.75 = 10.00 \]

แปลว่าจาก variance ทั้งหมด 10.00 หน่วย TEWL กำลังสอง ที่การวัดหนึ่งครั้งแบกไว้ มี 6.25 หน่วยที่ถูกแชร์กับการวัดครั้งอื่น ๆ ของคนเดียวกัน และอีก 3.75 หน่วยเป็นของแถวนั้นแถวเดียว อัตราส่วนระหว่างสองตัวเลขนี้มีชื่อเรียกของมัน และมันคือเนื้อหาทั้งหมดของโพสต์ถัดไป

สองตัวอย่างใกล้ตัวที่เป็นเรื่องเดียวกัน

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

ทั้งสองอุปมาบอกได้แค่ โครงสร้าง และไม่มีอันไหนอธิบาย กลไก ในตัวอย่างที่เราใช้อยู่ $b_i$ เป็นเพียงที่รองรับทุกอย่างในระดับบุคคลที่โมเดลไม่ได้ใส่ไว้ และไม่ควรถูกอ่านเป็นปริมาณทางชีววิทยา

อ่าน random = ~ 1 | id ทีละ token

ทุกอย่างข้างบนถูกสั่งให้โปรแกรมทำได้ด้วยบรรทัดสั้น ๆ บรรทัดเดียว และบรรทัดนี้ควรค่าแก่การอ่านช้า ๆ สักครั้ง มากกว่าการก๊อปวางไปเรื่อย ๆ มันมีแค่สี่ token คือเครื่องหมาย tilde เลขหนึ่ง ขีดตั้ง และชื่อคอลัมน์ที่บอกว่าใครเป็นใคร โดยแต่ละ token ทำหน้าที่เดียวเป๊ะ ๆ

คำสั่งสำหรับชุดข้อมูลจำลองนี้สั้นมาก

# Spec A — random intercept, pooled residual variance
library(nlme)
fitA <- lme(tewl ~ site * group + age + sex + phototype,
            random = ~ 1 | id, data = skin, method = "REML")

# equivalent in lme4
library(lme4)
fitA2 <- lmer(tewl ~ site * group + age + sex + phototype + (1 | id), data = skin)

* Stata equivalent of Spec A
* mixed estimates by ML unless reml is requested; lme() and lmer() above are REML
mixed tewl i.site##i.group age i.sex i.phototype || id:, reml

อ่าน random = ~ 1 | id จากในออกนอก เลข 1 คือคอลัมน์ของเลขหนึ่ง ซึ่งแปลว่า "หนึ่งค่า" ไม่ใช่ slope ของอะไร เครื่องหมาย | แปลว่า "สำหรับแต่ละระดับของ" และ id คือชื่อของระดับนั้น หนึ่งค่าต่อผู้เข้าร่วมหนึ่งคน รวม 80 ค่า ทุกอย่างที่อยู่ทางซ้ายของขีดตั้งคือสิ่งที่จะถูกใช้ซ้ำ ส่วนสิ่งที่อยู่ทางขวาเป็นตัวกำหนดว่าจะสุ่มใหม่บ่อยแค่ไหน คำสั่ง (1 | id) ใน lme4 และ || id: ใน Stata พูดสามเรื่องเดียวกันนี้ด้วยลำดับที่ต่างกันเท่านั้น

ที่สำคัญคือ ส่วน fixed effect ของสูตร — site * group + age + sex + phototype — ไม่ได้ถูกแตะเลย คำสั่ง random เปลี่ยนสิ่งที่โมเดลสมมติเกี่ยวกับ covariance ระหว่างแถว แต่ไม่ได้เปลี่ยนสิ่งที่โมเดลสมมติเกี่ยวกับค่าเฉลี่ย

ถ้าไม่ใส่ คุณกำลังอ้างอะไรอยู่

มีคำถามหนึ่งที่ถูกถามในห้องสอบเสมอ คือถ้าไม่สนใจเรื่องทั้งหมดนี้แล้ว fit regression ธรรมดาไปเลย จะเสียหายตรงไหน คำตอบที่ซื่อตรงไม่ใช่ "p-value จะผิด" แต่คือ "ขึ้นกับว่าคุณกำลังเปรียบเทียบอะไรอยู่" สามกรณีข้างล่างชี้ไปคนละทิศทาง และการแยกสามกรณีนี้ให้ออกก็คือทักษะทั้งหมดของเรื่องนี้

ถ้า fit ordinary least squares ลงบนทั้ง 320 แถว คุณกำลังอ้าง — ในภาษาคลินิกล้วน ๆ — ว่าการวัดครั้งที่สองของผู้เข้าร่วมที่อยู่ในการศึกษาแล้ว ให้ข้อมูลใหม่มากเท่ากับการวัดครั้งแรกของผู้เข้าร่วมที่ยังไม่เคยถูกตรวจเลย ภายใต้ข้อมูลจำลองชุดนี้ คำอ้างดังกล่าวผิดในลักษณะที่เจาะจงมาก นั่นคือ สองแถวจากคนเดียวกันแชร์ variance กันอยู่แล้ว 6.25 หน่วย จากทั้งหมด 10.00 หน่วย TEWL กำลังสอง ที่แต่ละแถวแบกไว้ ขณะที่สองแถวจากคนละคนไม่ได้แชร์อะไรกันเลย

สิ่งที่ตามมาไม่ใช่ความผิดพลาดเดียวที่มีทิศทางเดียว ทิศทางของมันขึ้นกับว่ากำลังดู contrast ตัวไหน จึงควรพูดให้ตรงมากกว่าพูดเป็นคำขวัญ ลองดู contrast สามแบบที่ดีไซน์นี้รองรับ แบบแรก ผลต่าง AD กับ control ที่ เฉลี่ยข้ามทั้งสี่ตำแหน่ง ผู้เข้าร่วมแต่ละคนส่งแถวที่สัมพันธ์กันเข้าไปคนละสี่แถว การวิเคราะห์ที่ละเลย offset จึงนับ 160 แถวต่อกลุ่มเป็น 160 แถวอิสระ standard error ที่ได้จึงเล็กเกินจริงไปมาก คือประมาณ 0.35 เทียบกับค่าจริง 0.60 ในข้อมูลจำลองชุดนี้ แบบที่สอง ผลต่าง AD กับ control ณ ตำแหน่งเดียว ตรงนี้ผู้เข้าร่วมแต่ละคนส่งเข้าไปคนละหนึ่งแถวพอดี ค่าที่นำมาเทียบกันเป็นอิสระต่อกันจริง ๆ การละเลย offset จึงแทบไม่มีต้นทุนอะไร standard error ที่รายงานออกมาราว 0.71 ถือว่าใกล้เคียงของจริง แบบที่สาม contrast ระหว่าง ตำแหน่งภายในคนเดียวกัน ทิศทางกลับด้าน เพราะการวิเคราะห์ที่ไม่มี person effect จะปล่อยให้ $b_i$ ค้างอยู่ในเทอม error ทำให้ residual variance พองขึ้น และทำให้การเปรียบเทียบที่แม่นยำจริง ๆ ดูไม่แม่นยำ

ที่น่าสนใจคือ แบบที่สองมักเป็นแบบที่คนคิดว่าพังแล้ว ทั้งที่ไม่ได้พัง ส่วนแบบแรกต่างหากที่พังอยู่เงียบ ๆ ดังนั้นประโยคเดียวแบบ "เราปรับสำหรับการวัดซ้ำแล้ว" แทบไม่ได้บอกอะไรผู้อ่านเลย สิ่งที่ต้องระบุคือตัว contrast เอง เพราะในสามค่านี้ standard error สองค่าผิดคนละทิศทาง และอีกค่าหนึ่งไม่ได้ผิดเลย

ผลตอบแทนแรก: ผลต่างที่วัดได้แม่นกว่าค่าเฉลี่ยที่ประกอบมันขึ้นมา

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

ภายใต้ข้อมูลจำลองนี้ ที่มีผู้เข้าร่วม $n = 40$ คนต่อกลุ่มและดีไซน์สมดุล standard error ของค่าเฉลี่ยที่ตำแหน่งหนึ่งภายในกลุ่มหนึ่งคือ

\[ \operatorname{SE}(\widehat{\mu}_{j}) = \sqrt{\frac{\tau^2 + \sigma^2}{n}} = \sqrt{\frac{10.00}{40}} = 0.50 \]

เพราะสิ่งที่คนหนึ่งคนใส่เข้าไปในค่าเฉลี่ยนั้นแบกความแปรปรวนมาทั้งสองแหล่ง ทั้ง offset และ noise ทีนี้ลองเปลี่ยนไปดู contrast ภายในคนเดียวกัน คือ forearm ลบด้วย back ของคนคนเดียวกัน แล้วดูว่าพีชคณิตทำอะไร

\[ Y_{i,\text{forearm}} - Y_{i,\text{back}} = (\mu_{\text{forearm}} - \mu_{\text{back}}) + (e_{i,\text{forearm}} - e_{i,\text{back}}), \qquad b_i - b_i = 0 \]

Offset ตัดกันหมดพอดีในทุกคน เพราะมันเป็นตัวเลขเดียวกันในทั้งสองแถว สิ่งที่เหลือคือผลต่างของ noise สองตัวที่เป็นอิสระต่อกัน standard error ของผลต่างภายในคนที่ประมาณได้จึงเป็น

\[ \operatorname{SE}(\widehat{\mu}_{\text{forearm}} - \widehat{\mu}_{\text{back}}) = \sqrt{\frac{2\sigma^2}{n}} = \sqrt{\frac{7.50}{40}} = 0.43 \]

ซึ่ง เล็กกว่า ค่า 0.50 ของค่าเฉลี่ยแต่ละตัวเสียอีก นี่ไม่ใช่ความย้อนแย้ง และไม่ใช่กฎธรรมชาติทั่วไป มันจริงเมื่อ $\sigma^2 < \tau^2$ เท่านั้น และในข้อมูลจำลองชุดนี้ 3.75 เล็กกว่า 6.25 จริง โมเดลที่ปฏิบัติต่อการวัดทั้งสี่ครั้งว่าเป็นอิสระต่อกันไม่มีทางผลิตเลข 0.43 ออกมาได้ และตัวเลขตัวเดียวนี้เองคือเหตุผลที่โพสต์สุดท้ายของซีรีส์จะปกป้องคอลัมน์ standard error ต่อหน้ากรรมการสอบได้

มีข้อกำกับหนึ่งข้อที่ต้องติดไปกับตัวเลขทั้งสองค่านี้ ทั้งคู่เป็นค่าที่ได้จากการแทนค่า (plug-in) นั่นคือสูตรปิดที่แทนด้วย variance component ชุดที่ใช้สร้างข้อมูลจำลอง สำหรับโมเดลค่าเฉลี่ยรายเซลล์แบบสมดุลที่ยังไม่ได้ปรับตัวแปรใด ๆ ไม่ใช่ผลลัพธ์ที่ได้จากโมเดลในบล็อกโค้ดข้างบน เมื่อ age, sex และ phototype เข้าไปอยู่ในส่วน fixed effect ค่าเฉลี่ยรายตำแหน่งและ contrast ระหว่างคนจะขยับไปเล็กน้อย ขณะที่ผลต่างภายในคนจะไม่ขยับเลย เพราะตัวแปรที่คงที่ภายในคนถูกตัดทิ้งไปพร้อมกับ $b_i$ ในการลบครั้งเดียวกันนั่นเอง

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

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

ข้อสรุปที่โพสต์นี้ปกป้องได้นั้นแคบมาก และข้อสรุปข้างเคียงที่ฟังดูคล้ายกันยังไม่ได้ถูกปกป้อง ข้อแรก การเห็น offset ชัด ๆ ในผู้เข้าร่วมสองคนไม่ใช่หลักฐานว่า random intercept เป็นโครงสร้าง covariance ที่ถูกต้อง offset ของ S03 ที่ +0.8 ไม่เคยแสดงเครื่องหมายที่คงเส้นคงวาเลย และการวัดสี่ครั้งต่อคนก็ไม่ได้เยอะอะไร ข้อสอง $b_i$ ไม่ใช่คำอธิบายของอะไรทั้งสิ้น มันดูดซับทุกอย่างที่เป็นระดับบุคคลและไม่ได้อยู่ในโมเดล การอ่านมันเป็นคุณสมบัติทางสรีรวิทยาของผิวหนังผู้ป่วยจะกลายเป็นการอ้างเชิงสาเหตุที่ดีไซน์นี้รองรับไม่ได้ ข้อสาม ประโยค "random intercept จัดการเรื่องการวัดซ้ำให้แล้ว" เป็นได้อย่างมากแค่ครึ่งประโยค เพราะเทอมนี้บังคับรูปแบบ covariance แบบเฉพาะเจาะจง ซึ่งความสมเหตุสมผลของมันต้องถูกตรวจสอบ ไม่ใช่ถูกสมมติ และโพสต์ที่ 2 กับที่ 6 จะตรวจสอบให้ ข้อสุดท้าย ไม่มี random effect ตัวไหนซ่อม mean structure ที่ผิดได้ ถ้า interaction ระหว่าง site กับ group ควรอยู่ในโมเดลแล้วเราไม่ใส่ ไม่มีอะไรในส่วน random ที่จะกู้มันกลับมาได้

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

1. พล็อตโปรไฟล์รายคนก่อน fit อะไรทั้งนั้น spaghetti plot ของแต่ละคนข้ามหน่วยที่วัดซ้ำอาจทำให้คุณเห็น offset ได้ตรง ๆ และมีต้นทุนแค่โค้ดบรรทัดเดียว
2. ระบุตัวแปรกลุ่มอย่างตั้งใจ ตรวจว่าแถวที่คุณเชื่อว่าเป็นของคนเดียวกันใช้ identifier เดียวกันจริง เพราะ id ที่ระบุผิดจะแปลง person effect ให้กลายเป็น noise แบบเงียบ ๆ
3. รายงาน $\widehat{\tau}^{2}$ และ $\widehat{\sigma}^{2}$ เป็นตัวเลข ประโยค "ใส่ random intercept ของผู้เข้าร่วมไว้แล้ว" ยังไม่ถือว่ารายงานครบ variance component ทั้งสองตัวควรถูกรายงานคู่ไปกับ fixed effect เพื่อให้ผู้อ่านประกอบภาพได้ว่าโมเดลสมมติอะไรไว้
4. ติดป้ายให้ทุก contrast ก่อนจะตีความ standard error ของมัน ว่าเป็นแบบภายในคน แบบระหว่างคนที่ผู้เข้าร่วมส่งข้อมูลคนละหนึ่งแถว หรือแบบระหว่างคนที่รวมแถวซ้ำของแต่ละคนเข้าด้วยกัน offset ตัดทิ้งในแบบแรก ไม่กัดในแบบที่สอง และเป็นตัวชี้ขาดในแบบที่สาม ทั้งสามแบบจึงพกปริมาณข้อมูลไม่เท่ากันแม้จะมาจากโมเดลเดียวกัน
5. เก็บตัวแปรต่อเนื่องระดับบุคคลให้ต่อเนื่องต่อไป การแบ่ง age เป็นกลุ่มเพื่อให้ตารางดูเรียบร้อยเป็นการทิ้งข้อมูลและอาจบิดเบือนค่าประมาณ ถ้าความสัมพันธ์ไม่เป็นเส้นตรง fractional polynomial หรือ spline เป็นทางเลือกแรกที่ดีกว่า
6. อย่าเพิ่งอ่าน offset ที่ประมาณได้ของคนหนึ่งคนเป็นค่าวัดทางคลินิก ว่า $\widehat{b}_{i}$ มีไว้ทำอะไรและไม่ได้มีไว้ทำอะไร จะถูกจัดการอย่างเป็นเรื่องเป็นราวในโพสต์ที่ 3

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

• Random intercept คือค่าเดียวต่อหนึ่งคน ที่ถูกบวกเข้าไปแบบไม่เปลี่ยนแปลงในทุกแถวของคนนั้น ไม่ใช่จุดตัดแกน Y ที่สุ่มมา และคำว่า "random" หมายถึงการที่ $b_i$ เป็นค่าที่สุ่มจากการแจกแจง ไม่ใช่พารามิเตอร์อิสระที่ต้องประมาณทีละคน
• เพราะตัวเลขเดียวกันปรากฏในทุกแถวของคนหนึ่งคน covariance ระหว่างการวัดสองครั้งใด ๆ ของคนนั้นจึงเท่ากับ $\operatorname{Var}(b_i) = \tau^2 = 6.25$ ซึ่งเป็น identity ไม่ใช่ค่าประมาณ และในนั้นไม่มี index ของตำแหน่งอยู่เลย
• การวัด TEWL จำลองหนึ่งครั้งแบก variance ไว้ $\tau^2 + \sigma^2 = 6.25 + 3.75 = 10.00$ หน่วย TEWL กำลังสอง โดย 6.25 หน่วยถูกแชร์กับการวัดครั้งอื่นของคนเดียวกัน และ 3.75 หน่วยเป็นของแถวนั้นเท่านั้น
• correlation ระหว่างการวัดซ้ำเป็นผลลัพธ์ของโครงสร้างข้อมูล ไม่ใช่ออปชันของโมเดล และการไม่ใส่ offset ไม่ใช่ความผิดพลาดเดียวที่มีทิศทางเดียว เพราะมันทำให้ standard error ของ contrast ที่รวมแถวซ้ำของแต่ละคนเล็กเกินจริง ปล่อย contrast ที่ผู้เข้าร่วมส่งคนละหนึ่งแถวไว้ใกล้เคียงของเดิม และทำให้การเปรียบเทียบภายในคนไม่แม่นเท่าที่ควรจะเป็น
• offset ตัดกันพอดีใน contrast ภายในคนเดียวกัน จึงเป็นเหตุผลที่ SE จำลองของ forearm − back เท่ากับ 0.43 เทียบกับ 0.50 ของค่าเฉลี่ยเดี่ยว ซึ่งจริงในที่นี้เพราะ $\sigma^2 < \tau^2$ ไม่ใช่จริงเป็นกฎสากล และตัวเลขทั้งสองเป็นค่าแทนค่าของโมเดลสมดุลที่ยังไม่ปรับตัวแปร ไม่ใช่ผลลัพธ์จากการ fit จริง
• การมองเห็น offset ด้วยตาไม่ได้ยืนยันว่าโครงสร้าง covariance ที่ random intercept บังคับใช้นั้นถูกต้อง และโครงสร้างนั้นคือเนื้อหาของซีรีส์ที่เหลือทั้งหมด

โพสต์ถัดไป

โพสต์ที่ 2 จะเติมช่องที่เหลืออีกสิบห้าช่อง: offset ที่แชร์กันหนึ่งค่าบวกกับ residual ที่มี variance เท่ากันให้ผลเป็น compound symmetry แบบเป๊ะ ๆ และอัตราส่วนระหว่าง 6.25 กับ 10.00 ก็มีชื่อที่นักระบาดวิทยาทุกคนใช้อยู่แล้ว — Compound Symmetry กับ ICC: ออฟเซ็ตเดียว correlation เดียว

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

1. Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982;38(4):963–974.