V = ZGZ′ + R: สมการเดียวที่อยู่เบื้องหลัง mixed model ทุกตัว

บทคัดย่อ

เลขคณิตของสามตอนก่อนยุบเหลือ matrix notation บรรทัดเดียว คือ marginal covariance ของการวัดจากคนหนึ่งคนเท่ากับ Z G Z ทรานสโพส บวก R บทความนี้หยุดบวกเลขแล้ววาดแผนที่ อธิบายแต่ละก้อนด้วยหนึ่งประโยค แทนค่าตัวอย่างจำลองประจำซีรีส์ และคูณ matrix ให้ดูจริงว่าคอลัมน์ของเลขหนึ่งสี่ตัว คูณ between-person variance 6.25 คูณแถวของเลขหนึ่งสี่ตัว ให้ค่า 6.25 ทุกช่อง และเมื่อบวก 3.75 บนเส้นทแยงมุมก็ได้ compound symmetry matrix ตัวเดิมที่สร้างด้วยมือในตอนที่ 1 ถึง 3 สมการนี้มีก้อนที่แก้ได้สองก้อน การตัดสินใจเรื่อง covariance ที่เหลือจึงแก้ได้ที่ Z กับ G หรือที่ R เท่านั้น และการเขียนสมการออกมาได้ไม่เท่ากับมีเหตุผลรองรับ

---

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

ตอนนี้เรามาถึงไหนแล้ว โพสต์ที่ 3 — Conditional หรือ Marginal? สองคำถามในโมเดลเดียว — ปิดงานเลขคณิตด้วยการเขียน covariance matrix ทั้งสองแบบของผู้เข้าร่วมจำลองคนเดียวกันออกมาทีละช่อง: เป็นอิสระต่อกันเมื่อรู้ว่าเป็นคนไหน และ correlate กันที่ 0.625 เมื่อไม่รู้ว่าเป็นคนไหน โพสต์นี้จะหยุดบวกเลขแล้วเริ่มวาดแผนที่ เพราะทุกอย่างที่คุณอ่านมาสามตอนคือกรณีพิเศษหนึ่งกรณีของ matrix notation บรรทัดเดียว ส่วนโพสต์ที่ 5 — Random Slope: เมื่อแต่ละคนพก "เส้นทาง" ของตัวเองมาด้วย — คือครั้งแรกที่เราจะจงใจไปแก้อะไรบางอย่างข้างในสมการนี้

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

Matrix ที่คุณสร้างมาแล้วสามรอบ

นี่คือสิ่งเดิมอีกครั้ง — marginal covariance matrix ของค่า transepidermal water loss (TEWL) สี่ค่าที่ผู้เข้าร่วมหนึ่งคนให้กับเรา เรียงตามลำดับตำแหน่งที่ตรึงไว้ตลอดซีรีส์คือ forearm, hand, shin, back:

\[ V=\begin{pmatrix} 10.00 & 6.25 & 6.25 & 6.25\\ 6.25 & 10.00 & 6.25 & 6.25\\ 6.25 & 6.25 & 10.00 & 6.25\\ 6.25 & 6.25 & 6.25 & 10.00 \end{pmatrix} \]

คุณเดินมาถึง matrix ตัวนี้ด้วยเส้นทางที่ต่างกันสามเส้นแล้ว ในโพสต์ที่ 1 คุณได้มาหนึ่งช่องจากการตัดทุกอย่างที่ไม่ได้แชร์กันระหว่างสองแถวของคนเดียวกันทิ้ง เหลือ $\operatorname{Cov}(Y_{ij},Y_{ik}) = \operatorname{Var}(b_i) = \tau^2 = 6.25$ ในโพสต์ที่ 2 คุณเติมอีกสิบห้าช่องที่เหลือด้วยข้อสังเกตว่าไม่มีตำแหน่งไหนพิเศษกว่ากัน จึงได้ compound symmetry และ intraclass correlation เท่ากับ $6.25/10.00 = 0.625$ และในโพสต์ที่ 3 คุณสร้างตัวเลขชุดเดิมขึ้นมาใหม่เป็นรอบที่สามจาก law of total variance และ law of total covariance แล้วเห็นค่านอกเส้นทแยงมุมโผล่ขึ้นมาทันทีที่ offset ของคนถูกซ่อนกลับไป

เมื่อของสิ่งเดียวเดินมาถึงเราได้จากสามเส้นทางที่เป็นอิสระต่อกัน โดยทั่วไปมันมักเป็นกรณีพิเศษของอะไรที่ใหญ่กว่า และในที่นี้ก็เป็นอย่างนั้นจริง ๆ linear mixed model ทุกตัวที่เคย fit กันมา — ไม่ว่าจะมี random intercept ตัวเดียวหรือ random effect ห้าตัว วัดสี่ตำแหน่งบนร่างกายหรือสี่สิบ visit ข้อมูลสมดุลหรือหลุดเป็นหย่อม ๆ — ล้วนมี marginal covariance matrix ในรูป $Z_i G Z_i' + R_i$ และ matrix ที่พิมพ์อยู่ข้างบนคือสิ่งที่นิพจน์นี้คืนออกมาเมื่อโมเดลมี offset ร่วมเพียงค่าเดียว บวกกับ residual ที่เป็นอิสระต่อกันและมี variance เท่ากันทุกตัว มีข้อกำกับหนึ่งข้อที่ควรพูดไว้ตรงนี้มากกว่าจะซุกไว้เป็นเชิงอรรถ นั่นคือการเขียน covariance เป็นก้อนต่อคนหนึ่งก้อนอย่างที่โพสต์นี้ทำตลอด ยังสมมติเพิ่มด้วยว่าแถวข้อมูลแบ่งออกเป็นคลัสเตอร์รายคนที่เป็นอิสระต่อกัน รูปทั่วไปนั้นไม่สะเทือนเมื่อเจอ crossed random effect แต่ก้อนรายคน $V_i$ สะเทือน

ผลตอบแทนที่ได้ไม่ใช่ความสวยงาม แต่คือสมการนี้มีที่ให้วางการตัดสินใจเรื่อง covariance ได้ สองที่เท่านั้น คนที่เป็นเจ้าของสมการนี้จึงจัดหมวดหมู่ตัวเลือกเรื่อง covariance ของ mixed model ทุกตัวที่จะเจอในชีวิต — random slope, residual variance ที่ไม่เท่ากันในแต่ละหน่วยวัด, autoregressive error, unstructured covariance — ได้ด้วยการชี้ไปที่ก้อนใดก้อนหนึ่งในสองก้อนนั้น

ก้อนต่าง ๆ ในสมการ ก้อนละหนึ่งประโยค

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

Linear mixed model รูปทั่วไป ในรูปแบบสองขั้นตอนตามที่ Laird และ Ware เสนอไว้ เขียนเวกเตอร์การวัดทั้งหมดของคนหนึ่งคนพร้อมกันทีเดียว:

\[ Y_i = X_i\beta + Z_i b_i + e_i, \qquad b_i \sim N(0, G), \quad e_i \sim N(0, R_i), \quad b_i \perp e_i \]

อ่านจากซ้ายไปขวา ทีละก้อน

1. $Y_i$ คือกองข้อมูลของคนที่ $i$ — ในตัวอย่างของเราคือเวกเตอร์ขนาด $4\times1$ ที่บรรจุค่าที่ forearm, hand, shin และ back ของผู้เข้าร่วมคนนั้น
2. $X_i\beta$ คือสิ่งที่ ประชากร ทำ นั่นคือ design matrix ของ fixed effect คูณกับสัมประสิทธิ์ชุดที่ทุกคนใช้ร่วมกัน
3. $Z_i b_i$ คือการที่ คนคนนี้ เบี่ยงออกจากประชากร ประกอบด้วย design matrix $Z_i$ ซึ่งบอกว่าอนุญาตให้เบี่ยงในทิศทางใดได้บ้าง และเวกเตอร์ $b_i$ ซึ่งบอกว่าสำหรับคนนี้เบี่ยงไปเท่าไร
4. $e_i$ คือสิ่งที่เหลือหลังจากหักสองก้อนแรกออกแล้ว — สัญญาณรบกวนระหว่างการวัดแต่ละครั้งที่ทั้งโครงสร้างประชากรและ offset ส่วนตัวของคนอธิบายไม่ได้
5. $G$ คือ covariance matrix ของการเบี่ยงส่วนบุคคล ส่วน $R_i$ คือ covariance matrix ของส่วนที่เหลือ

มิติของแต่ละก้อนสำคัญกว่าที่ตาเห็น เพราะมันคือสิ่งที่ทำให้การคูณ matrix ในหัวข้อถัดไปทำงานได้:

ก้อน	มิติ (ตัวอย่างในซีรีส์)	บรรจุอะไร
$Y_i$	4 × 1	ค่า TEWL สี่ค่า ตำแหน่งละหนึ่งค่า
$X_i$	4 × p	คอลัมน์ fixed effect ของสี่แถวของคนนี้
$\beta$	p × 1	สัมประสิทธิ์ระดับประชากร
$Z_i$	4 × 1	คอลัมน์ของเลขหนึ่ง — offset ร่วมค่าเดียว ใส่ให้ทุกแถว
$b_i$	1 × 1	offset ของคนนี้
$G$	1 × 1	$\tau^2$ ความแปรปรวนของ offset เหล่านั้น
$e_i$	4 × 1	residual สี่ค่า
$R_i$	4 × 4	covariance matrix ของ residual

สองประโยคจากโมเดลเดียว

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

เมื่อเรียกชื่อทุกก้อนได้แล้ว โมเดลนี้พูดข้อความเชิงความน่าจะเป็นสองประโยคที่ต่างกัน และเนื้อหาทั้งโพสต์ที่ 3 อยู่ในช่องว่างระหว่างสองประโยคนี้

เมื่อ condition บน random effect ของคนคนนั้น — เปิดซองจดหมายออกแล้ว — การวัดทั้งสี่ครั้งกระจายแบบ normal รอบค่าเฉลี่ยที่รวม offset ของคนนี้ไว้ด้วย โดยมี covariance เป็น $R_i$ เท่านั้น ไม่มีอะไรอื่น:

\[ Y_i \mid b_i \sim N\!\left(X_i\beta + Z_i b_i,\ R_i\right) \]

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

\[ Y_i \sim N\!\left(X_i\beta,\ Z_i G Z_i' + R_i\right) \]

บรรทัดที่สองนี้ไม่ใช่การพิสูจน์ใหม่ มันคือการอินทิเกรตที่คุณทำด้วยมือไปแล้วในโพสต์ที่ 3 สำหรับ $k = 4$ เพียงแต่เขียนครั้งเดียวให้ครอบคลุมทุก $Z$ ทุก $G$ และทุก $R$ พร้อมกัน และเป็นรูปแบบที่ Verbeke และ Molenberghs ใช้นำเสนอ marginal covariance ของ linear mixed model สังเกตว่าอะไรขยับและอะไรไม่ขยับ: random effect ไปเปลี่ยน covariance แต่ปล่อยค่าเฉลี่ยไว้เฉย ๆ นี่เป็นคุณสมบัติของ identity link ไม่ใช่กฎสากล และโพสต์ที่ 7 จะมาเก็บหนี้ก้อนนี้เมื่อ outcome เปลี่ยนเป็น binary

แทนค่าตัวอย่างที่เราใช้ทั้งซีรีส์

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

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

การศึกษาจำลองนี้ให้ $k = 4$ ตำแหน่งต่อผู้เข้าร่วมหนึ่งคน และ offset ร่วมหนึ่งค่าต่อหนึ่งคน ดังนั้น $Z_i$ คือคอลัมน์ของเลขหนึ่งขนาด $4\times1$ และ $G$ คือ matrix ขนาด $1\times1$ ที่บรรจุ $\tau^2 = 6.25$ ส่วน residual เป็นอิสระต่อกันและใช้ variance ค่าเดียวร่วมกัน ดังนั้น $R_i = \sigma^2 I_4 = 3.75\,I_4$ ก้อนซ้ายมือคูณออกมาเป็น $4\times1$ คูณสเกลาร์ คูณ $1\times4$:

\[ Z_i G Z_i' = \begin{pmatrix}1\\1\\1\\1\end{pmatrix} (6.25) \begin{pmatrix}1&1&1&1\end{pmatrix} = 6.25\begin{pmatrix}1&1&1&1\\1&1&1&1\\1&1&1&1\\1&1&1&1\end{pmatrix} = 6.25\,J_4 \]

ทุกช่องของผลคูณนี้เท่ากับ $6.25$ — รวมถึงช่องบนเส้นทแยงมุมด้วย — เพราะทุกช่องคือ $1 \times 6.25 \times 1$ เหมือนกันหมด เมื่อบวกก้อน residual เข้าไป $3.75$ จะไปลงเฉพาะเส้นทแยงมุม:

\[ V_i = Z_i G Z_i' + R_i = 6.25\,J_4 + 3.75\,I_4 = \begin{pmatrix} 10.00 & 6.25 & 6.25 & 6.25\\ 6.25 & 10.00 & 6.25 & 6.25\\ 6.25 & 6.25 & 10.00 & 6.25\\ 6.25 & 6.25 & 6.25 & 10.00 \end{pmatrix} \]

Matrix จากโพสต์ที่ 2 กลับมาแล้ว ตรงกันทุกหลัก โดยไม่มีใครไปตัดสินใจว่า correlation ควรเป็นเท่าไร ประโยคที่ควรหิ้วกลับบ้านจากโพสต์นี้คือ compound symmetry ไม่ใช่ตัวเลือกที่ซอฟต์แวร์ยื่นให้แล้วเรารับไว้ แต่มันคือสิ่งที่คอลัมน์ของเลขหนึ่ง บวกกับ residual variance ค่าเดียวที่ใช้ร่วมกัน ผลิตออกมา ประโยคนี้มีสองท่อน และทั้งสองท่อนรับน้ำหนักคนละส่วน โครงสร้าง "หนึ่งทั้งกระดาน" ของ $J_4$ มาจากการที่ $Z_i$ เป็นเลขหนึ่งทั้งคอลัมน์ ซึ่งเป็นสิ่งที่ตรึงทุกช่องนอกเส้นทแยงมุมไว้ที่ $6.25$ ส่วนเส้นทแยงมุมที่เท่ากันมาจากการที่ $R_i$ เป็น variance ค่าเดียวคูณ identity และ intraclass correlation ที่ $6.25/10.00 = 0.625$ คือผลลัพธ์ทางเลขคณิตของทั้งสองอย่างรวมกัน ถ้าเก็บคอลัมน์ของเลขหนึ่งไว้เหมือนเดิมแต่ปล่อยให้ residual variance ต่างกันไปในแต่ละตำแหน่ง ค่านอกเส้นทแยงมุมจะยังราบเท่ากันที่ $6.25$ อยู่ ในขณะที่ compound symmetry ตายไปแล้วอย่างเป็นทางการ ซึ่งนั่นคือการสาธิตที่โพสต์ที่ 6 สร้างขึ้นมาทั้งโพสต์ เปลี่ยนอินพุตตัวใดตัวหนึ่ง เอาต์พุตก็เปลี่ยนตาม ซึ่งนั่นคือแผนงานของโพสต์ที่ 5 และ 6 พอดี

อ่านเปเปอร์แล้วแปลว่าอะไร: หัวข้อวิธีการที่เขียนเพียงว่า "ใส่ random intercept ของผู้เข้าร่วมไว้ด้วย" ได้บอกโครงสร้าง covariance ที่โมเดลสมมติไว้ให้คุณครบแล้ว คือ compound symmetry ต่อให้ไม่มีคำนี้ปรากฏอยู่เลยสักคำ

อ่าน design matrix $X_i$

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

ฝั่งขวาของโมเดลยังมี matrix อีกตัวหนึ่ง และเป็นตัวที่ผู้อ่านส่วนใหญ่ข้ามไป สำหรับโมเดลที่เราใช้ — TEWL กับ site, group, interaction ของทั้งสอง และตัวแปรระดับบุคคล — $X_i$ มีสี่แถว ตำแหน่งละหนึ่งแถว และคอลัมน์ของมันแบ่งได้เป็นสองตระกูล:

คอลัมน์ของ $X_i$	จำนวนคอลัมน์	เปลี่ยนค่าภายในคนเดียวกันหรือไม่
Intercept	1	ไม่
ตัวบ่งชี้ site (hand, shin, back โดยมี forearm เป็น reference)	3	ใช่
ตัวบ่งชี้ group (AD เทียบ control)	1	ไม่
Interaction ระหว่าง site กับ group	3	ใช่
age	1	ไม่
sex	1	ไม่
phototype (Fitzpatrick ในฐานะ factor)	น้อยกว่าจำนวนระดับที่ปรากฏอยู่หนึ่ง	ไม่

ไล่สายตาลงมาตามสี่แถวของผู้เข้าร่วมหนึ่งคน จะมีแค่สองตระกูลที่ขยับ คือตัวบ่งชี้ site และเทอม interaction ส่วน age, sex, phototype และป้าย group ถูกพิมพ์ค่าเดิมซ้ำครบทั้งสี่แถว เพราะมันเป็นคุณสมบัติของ คน ไม่ใช่ของ การวัด นี่ไม่ใช่รายละเอียดทางบัญชี แต่คือเมล็ดพันธุ์ของข้อโต้แย้งในโพสต์ที่ 8 ที่ standard error สี่ค่าออกมาเท่ากันเป๊ะ: contrast ที่อยู่ในคอลัมน์ที่ขยับทั้งหมดสามารถประมาณค่าได้ภายในคนคนเดียว จึงไม่ต้องจ่ายค่า between-person variance เลย ในขณะที่ contrast ที่ต้องใช้คอลัมน์ที่ไม่ขยับทำแบบนั้นไม่ได้

ตัวเลขจำลองทำให้ประเด็นนี้จับต้องได้ ถ้าเราผูกตัวเลขแต่ละตัวไว้กับปริมาณที่มันประมาณจริง ๆ ค่าเฉลี่ยของกลุ่มหนึ่งที่ตำแหน่งหนึ่ง ซึ่งเป็นค่าเฉลี่ยจากผู้เข้าร่วม 40 คนของกลุ่มนั้น คนละหนึ่งแถว มี standard error เท่ากับ $\sqrt{10.00/40} = 0.50$ ส่วน contrast AD − control ที่ตำแหน่งเดียวกันนั้น ซึ่งเอาคนสองชุดชุดละ 40 คนที่ไม่ทับกันมาเทียบกัน มี standard error เท่ากับ $\sqrt{2\times10.00/40} = 0.71$ และผลต่างระหว่างสองตำแหน่งภายในคนเดียวกันอย่าง forearm − back มี standard error เท่ากับ $\sqrt{2\times3.75/40} = 0.43$ ซึ่งเล็กกว่า standard error ของค่าเฉลี่ยรายตำแหน่งแต่ละตัวเสียอีก ทั้งสามค่านี้เป็นค่ารูปปิดของโมเดลค่าเฉลี่ยรายเซลล์แบบสมดุลที่ยังไม่ได้ปรับตัวแปรร่วม การใส่คอลัมน์ตัวแปรร่วมของ $X_i$ เข้าไปใน fit จะขยับสองค่าแรกออกจากตัวเลขข้างต้นไปเล็กน้อย ในขณะที่ผลต่างภายในคนไม่ขยับตาม เพราะทุกคอลัมน์ที่คงที่ภายในคนเดียวกันตัดทิ้งกันหมด — เป็นการตัดกันชุดเดียวกับที่โพสต์ที่ 8 จะใช้เวลาทั้งโพสต์ปกป้อง เพียงแต่มองจากมุมของ design matrix

อ่านเปเปอร์แล้วแปลว่าอะไร: การเปรียบเทียบภายในคนกับการเปรียบเทียบระหว่างกลุ่มที่พิมพ์อยู่ในตารางเดียวกัน ไม่ได้มีความแม่นยำเท่ากัน และ design matrix คือที่ที่ความต่างนั้นมองเห็นได้

แผนที่: มีให้แก้แค่สองที่เท่านั้น

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

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

คุณเปลี่ยน…	ได้อะไรมา	จ่ายอะไรไป
$Z_i$ กับ $G$ — เติมคอลัมน์ให้ $Z$ ขยาย $G$ ให้ใหญ่ขึ้น	ให้แต่ละคนเบี่ยงได้มากกว่าค่าคงที่ค่าเดียว เช่น มีเส้นทางของตัวเอง มีแนวโน้มเฉพาะตำแหน่ง หรือมีการตอบสนองต่อการรักษาของตัวเอง marginal variance จะเลิกคงที่ และ correlation จะเลิกเท่ากันทุกคู่	พารามิเตอร์ variance และ covariance เพิ่มขึ้น พื้นผิว likelihood ยากขึ้น และค่าประมาณอาจตกไปที่ขอบเขตเมื่อข้อมูลรองรับไม่ไหว — โพสต์ที่ 5
$R_i$ — ผ่อนรูปทรงของ residual covariance	residual variance ที่ต่างกันได้ในแต่ละหน่วยวัด residual correlation ที่ลดลงตามระยะห่าง หรือ matrix แบบ unstructured ที่ไม่สมมติอะไรเลย	พารามิเตอร์หมดเร็ว เพดาน identifiability มาถึงเร็วกว่าที่คนส่วนใหญ่คิด และเสีย $\tau^2$ ในฐานะปริมาณที่แยกออกมาได้เมื่อชนเพดานนั้น — โพสต์ที่ 6

แผนที่นี้เป็นข้อความเรื่อง notation ไม่ใช่ใบอนุญาตให้สร้างโมเดล โครงสร้าง covariance ใด ๆ ที่คุณนึกออกล้วนเขียนในรูปนี้ได้ แต่การที่โครงสร้างหนึ่ง เขียนได้ ไม่ใช่หลักฐานแม้แต่น้อยว่ามัน มีเหตุผลรองรับ สำหรับข้อมูลของคุณ การเลือกว่าจะ fit $G$ แบบไหนและ $R$ แบบไหนยังคงเป็นการตัดสินใจเชิงเนื้อหา ที่ต้องเถียงจากดีไซน์ซึ่งผลิตข้อมูลชุดนั้นขึ้นมา และจาก diagnostics ของโมเดลที่ fit แล้ว ไม่ใช่เถียงจากว่าเอาต์พุตอันไหนดูเรียบร้อยกว่ากัน

ซอฟต์แวร์พิมพ์ $G$ กับ $R$ ไว้ตรงไหน

ทั้งสองก้อนถูกพิมพ์ออกมาโดยทุกโปรแกรมมาตรฐาน เพียงแต่อยู่คนละหัวข้อ นี่คือโมเดล Spec A ที่ fit สองวิธีใน R:

# 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

ใน R คำสั่ง VarCorr() จะรายงานตาราง random effect ซึ่งก็คือ $G$ — ในที่นี้คือตัวเลขตัวเดียวคือ $\tau^2$ — พร้อมกับเทอม residual ของ $R$ ให้อ่านชื่อคอลัมน์ก่อนอ่านตัวเลข เพราะสิ่งที่อยู่ข้างใน $R$ คือ residual variance $\sigma^2 = 3.75$ ในขณะที่เอาต์พุตหลายแบบพิมพ์ออกมาเฉพาะรากที่สองของมัน คือ residual standard deviation ประมาณ $1.94$ และสองค่านี้ใช้แทนกันไม่ได้ ส่วนฟังก์ชัน getVarCov() จะส่ง matrix ที่ประกอบแล้วมาให้ก็จริง แต่ค่าเริ่มต้นไม่ใช่ตัวที่เราต้องการในที่นี้ เพราะค่าเริ่มต้นคือ type = "random.effects" ซึ่งคืน $G$ ออกมา ดังนั้น $V_i$ ที่พิมพ์ไว้ข้างบนต้องสั่งขอตรง ๆ ด้วย getVarCov(fitA, individuals = "S01", type = "marginal") สำหรับ Stata ตาราง random effect และแถบ residual ของคำสั่ง mixed บรรจุของสองสิ่งเดียวกันในลำดับเดียวกัน โดยรายงานเป็น variance เว้นแต่จะสั่ง stddeviations การลองประกอบ $V_i$ ด้วยมือจากสองแถบนั้นสักครั้งตอนเริ่มโครงการ เป็นนิสัยเล็ก ๆ ที่อาจช่วยประหยัดความสับสนได้มากในภายหลัง เพราะหลังจากนั้นทุกบรรทัดของเอาต์พุตจะมีที่อยู่ที่แน่นอนอยู่ในสมการ

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

ความเป็นรูปทั่วไปของ $V_i = Z_i G Z_i' + R_i$ ถูกอ่านเกินเลยได้ง่ายมาก จึงควรตีเส้นไว้สี่เส้นให้ชัด

ข้อแรก สมการนี้เป็นภาชนะ ไม่ใช่คำแนะนำ มันรองรับ compound symmetry ได้สบายพอ ๆ กับที่รองรับโครงสร้างที่ขัดแย้งกับ compound symmetry และมันไม่มีความเห็นใด ๆ ว่าข้อมูลของคุณสมควรได้แบบไหน

ข้อสอง $G$ กับ $R$ ไม่ได้แยกกู้คืนจากข้อมูลได้เสมอไป — $R$ ที่อิสระมากพอสามารถดูดซับทุกอย่างที่ $G$ เคยแบกไว้ได้หมด และการนับพารามิเตอร์ที่แสดงว่าเหตุการณ์นี้เกิดขึ้นเมื่อไรคือธุระของโพสต์ที่ 6

ข้อสาม $b_i$ ไม่ใช่พารามิเตอร์รายคนที่โมเดลประมาณค่าให้ ค่าที่ซอฟต์แวร์คืนมาคือค่าทำนายแบบหดตัวของปริมาณสุ่มตัวหนึ่ง และควรถูกอ่านในฐานะนั้น

ข้อสุดท้าย การแบ่งงานอย่างเรียบร้อยระหว่าง "random effect ไปเปลี่ยน covariance" กับ "fixed effect แบกค่าเฉลี่ย" เป็นคุณสมบัติของ linear model ที่ใช้ identity link และมันไม่รอดข้ามไปถึง outcome แบบ binary

สิ่งที่จงใจไม่พูดถึงในซีรีส์นี้

ซีรีส์นี้ละเว้นหกเรื่องโดยตั้งใจ และขอเรียกชื่อทุกเรื่องหนึ่งครั้ง แทนที่จะเงียบ ๆ ข้ามไป

1. กลไกการประมาณค่าแบบ REML เทียบกับ ML — สองวิธีต่างกันที่การประมาณ variance component และที่ว่า likelihood แบบใดเทียบกันได้บ้าง ผลกระทบเชิงปฏิบัติจะพูดตรงจุดที่มันมีผลจริง คือโพสต์ที่ 5 ไม่ใช่พิสูจน์ที่นี่
2. Denominator degrees of freedom (Satterthwaite, Kenward–Roger) — เป็นปัญหาการประมาณค่าเรื่อง การทดสอบ ไม่ใช่เรื่องโครงสร้าง covariance ที่ซีรีส์นี้สร้างขึ้นมารอบ ๆ โพสต์ที่ 8 จะให้พื้นที่หนึ่งย่อหน้าในฐานะรายการที่ต้องรายงาน
3. Bayesian mixed model — $Z$, $G$ และ $R$ ชุดเดิมปรากฏอยู่ครบ แต่ prior และการคำนวณ posterior เป็นหลักสูตรอีกวิชาหนึ่ง
4. Crossed และ partially crossed random effects — รูปทั่วไป $V = ZGZ' + R$ ขยายไปรองรับได้โดยไม่ต้องแก้อะไร แต่ก้อนรายคน $V_i$ ขยายตามไปไม่ได้ เพราะแถวข้อมูลไม่ได้แบ่งเป็นคลัสเตอร์รายคนที่เป็นอิสระต่อกันอีกต่อไป และสัญชาตญาณแบบ "หนึ่งค่าต่อหนึ่งคน" ก็ไปต่อไม่ได้เช่นกัน จึงสมควรมีที่ทางของมันเอง
5. Generalised และ non-linear mixed model ที่เกินกว่ากรณี logistic กรณีเดียวในโพสต์ที่ 7 — ซีรีส์นี้ต้องการเพียงแง่มุมที่ conditional กับ marginal แยกทางกันเท่านั้น
6. กลไกของข้อมูลสูญหาย ที่เกินกว่าข้อโต้แย้งเรื่องความสมดุลในโพสต์ที่ 8 — การอธิบายให้ครบต้องยกอนุกรมวิธานของ missingness มาทั้งชุด

ทำอะไรได้บ้างในงานวิเคราะห์ของคุณ

ข้อแรก คราวหน้าที่คุณอ่านเอาต์พุตของ mixed model ให้หาตำแหน่งของ $G$ และ $R$ ให้เจอก่อนอ่านสัมประสิทธิ์แม้แต่ตัวเดียว เพราะสมมติฐานเรื่อง covariance เป็นตัวกำหนด standard error ทุกค่าที่พิมพ์อยู่ข้างล่างนั้น

ข้อสอง ลองเขียน $Z_i$ ของผู้เข้าร่วมหนึ่งคนออกมาด้วยมือ — สี่แถว กับจำนวนคอลัมน์เท่าที่ random effect ของคุณกำหนด — โดยทั่วไปแค่นี้ก็มักพอจะเผยให้เห็นว่าคุณกำลังขอโครงสร้างที่ดีไซน์รองรับไหวจริงหรือไม่

ข้อสาม เมื่อใดที่เจอตัวเลือกของ mixed model ที่ไม่คุ้นเคยในคู่มือหรือในคอมเมนต์ของ reviewer ให้ถามว่ามันไปแก้ก้อนไหนในสองก้อน คำถามเดียวนี้อาจใช้แทนความสับสนเรื่องศัพท์ covariance structure ได้เกือบทั้งหมด

ข้อสี่ ระบุไว้ในแผนการวิเคราะห์ว่าคอลัมน์ใดของ $X_i$ คงที่ภายในคนเดียวกัน เพราะ contrast ที่สร้างจากคอลัมน์เหล่านั้นกับ contrast ที่สร้างจากคอลัมน์ที่ขยับภายในคน จะมีความแม่นยำไม่เท่ากัน และอาจคุ้มค่าที่จะเขียนบอกไว้ในต้นฉบับก่อนที่ reviewer จะพูดแทนคุณ

ข้อสุดท้าย ระบุโครงสร้างที่สมมติของ $G$ และ $R$ ไว้ในหัวข้อวิธีการอย่างชัดเจน แทนที่จะบอกแค่ชื่อคำสั่งของซอฟต์แวร์ เพราะชื่อคำสั่งไม่ใช่สมมติฐานเรื่อง covariance และผู้อ่านย่อมสร้างอย่างหลังขึ้นมาจากอย่างแรกไม่ได้

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

• Linear mixed model ทุกตัว ไม่ว่าซับซ้อนเพียงใด มี marginal covariance เท่ากับ $Z_i G Z_i' + R_i$ และ matrix แบบ compound symmetry ที่เราสร้างด้วยมือมาตั้งแต่โพสต์ที่ 1 ถึง 3 ก็คือสิ่งที่นิพจน์นี้คืนออกมาเมื่อ $Z_i$ เป็นคอลัมน์ของเลขหนึ่ง และ $R_i$ เป็น variance ค่าเดียวคูณ identity
• การคูณควรทำให้เห็นสักครั้ง: คอลัมน์เลขหนึ่งขนาด $4\times1$ คูณ $\tau^2 = 6.25$ คูณแถวเลขหนึ่งขนาด $1\times4$ ได้ $6.25\,J_4$ และเมื่อบวก $3.75\,I_4$ ก็ได้เส้นทแยงมุม 10.00 กับค่านอกเส้นทแยงมุม 6.25 กลับมาตรงเป๊ะ
• ดังนั้น compound symmetry จึงไม่ใช่ตัวเลือกที่ซอฟต์แวร์เลือกให้แทนคุณ แต่เป็นผลลัพธ์ทางเลขคณิตของ offset ร่วมหนึ่งค่า บวกกับ residual ที่เป็นอิสระและมี variance เท่ากัน
• สมการนี้มีก้อนที่แก้ได้เพียงสองก้อน ทุกการตัดสินใจเรื่อง covariance ที่รออยู่ข้างหน้าจึงเป็นการแก้ $Z$ กับ $G$ (โพสต์ที่ 5) หรือแก้ $R$ (โพสต์ที่ 6) อย่างใดอย่างหนึ่ง
• ใน $X_i$ คอลัมน์ของ site ขยับภายในคนเดียวกัน ส่วน age, sex, phototype และ group ไม่ขยับ ความไม่สมมาตรนี้คือสิ่งที่จะอธิบาย standard error สี่ค่าที่เท่ากันในโพสต์ที่ 8
• เขียนได้ไม่เท่ากับมีเหตุผลรองรับ รูปทั่วไปนี้ยินดีจดโครงสร้าง covariance ที่ดีไซน์ของคุณรองรับไม่ไหวให้เสมอ และมีเพียงดีไซน์กับ diagnostics เท่านั้นที่บอกได้ว่าควร fit แบบไหน

ตอนต่อไป ใน Random Slope: เมื่อแต่ละคนพก "เส้นทาง" ของตัวเองมาด้วย เราจะเติมคอลัมน์ที่สองให้ $Z$ ดู $G$ โตขึ้นเป็น matrix ขนาด 2×2 และเห็น variance เลิกคงที่

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

1. Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982;38(4):963–974.
2. Verbeke G, Molenberghs G. Linear Mixed Models for Longitudinal Data. Springer; 2000.