Random Slope: เมื่อแต่ละคนพก "เส้นทาง" ของตัวเองมาด้วย

บทคัดย่อ

Random intercept ยกโปรไฟล์ทั้งเส้นขึ้นลงโดยไม่เปลี่ยนรูปทรง จึงทำให้ความกระจายที่ครั้งวัดหนึ่งต่างจากอีกครั้งไม่ได้ บทความนี้เติมคอลัมน์ที่สองให้ Z ให้แต่ละคนพกทั้งค่าชดเชยและอัตราการเปลี่ยนแปลงของตัวเอง และทำให้ G ขยายเป็น matrix ขนาด 2×2 ในการศึกษาย่อยแบบสุ่มที่จำลองขึ้น marginal variance เป็นฟังก์ชันกำลังสองของเวลา คือ 8.25, 7.81, 7.69 และ 8.41 ที่สัปดาห์ 0, 2, 4 และ 8 ส่วน correlation ทั้งหกไล่จาก 0.61 ถึง 0.74 แทนที่จะเท่ากัน ซึ่งปิดฉาก compound symmetry ผู้อ่านจะเข้าใจว่าทำไม correlation ของ intercept กับ slope จึงพูดถึงจุดอ้างอิงของเวลาเท่านั้น singular fit บอกอะไร และควรระบุ random structure ล่วงหน้าแค่ไหน

---

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

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

ตอนที่ 4 เป็นตอนที่ยุบ linear mixed model ทุกตัวลงเหลือบรรทัดเดียวนั้น และทิ้งแผนที่ไว้ให้ว่ามีอยู่แค่สองที่เท่านั้นที่เราเปลี่ยนอะไรได้ โพสต์นี้เปลี่ยนที่แรก คือเติมคอลัมน์ที่สองให้ $Z$ เพื่อให้แต่ละคนไม่ได้พกแค่ค่าชดเชย แต่พก "อัตราการเปลี่ยนแปลง" ของตัวเองมาด้วย ส่วนตอนถัดไป เมื่อ Compound Symmetry ไม่จริง จะปล่อย $Z$ ไว้เฉย ๆ แล้วหันไปจัดการกับ $R$ แทน

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

แถวตัวเลขสี่ค่าที่ random intercept ซึ่งมี residual variance ค่าเดียวสร้างไม่ได้

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

นี่คือหลักฐานชิ้นแรกของโพสต์นี้ มันคือเส้นทแยงมุมของ marginal covariance matrix ซึ่งก็คือ variance ของ outcome ที่การนัดแต่ละครั้ง ในการศึกษาย่อยแบบสุ่มที่จำลองขึ้น ซึ่งคุณจะได้รู้จักในอีกไม่กี่บรรทัดข้างหน้า

สัปดาห์	0	2	4	8
Var(TEWL) (g·m⁻²·h⁻¹)²	8.25	7.81	7.69	8.41

สี่การนัด สี่ค่า variance ที่ไม่เท่ากันเลย ใน ตอนที่ 2 เส้นทแยงมุมเป็น 10.00 เท่ากันหมดทั้งสี่ตำแหน่ง และค่านอกเส้นทแยงมุมเป็น 6.25 ทุกช่อง ความ "แบน" นั้นไม่ได้เกิดจากความบังเอิญของการจำลอง แต่ถูกบังคับด้วยโครงสร้างของโมเดล เพราะ random intercept เติมตัวเลขค่าเดียวกันลงในทุกแถวของคนคนนั้น มันจึงยกทั้งโปรไฟล์ขึ้นหรือลง โดยไม่แตะรูปร่างของโปรไฟล์เลย specification แบบนั้นมีส่วนประกอบอยู่เพียงสองอย่าง คือค่าชดเชยหนึ่งค่าต่อหนึ่งคน กับ residual variance ค่าเดียวที่ใช้ร่วมกันทุกการนัด และไม่มีอย่างใดเลยที่จะทำให้การกระจายที่สัปดาห์ที่ 8 ต่างจากการกระจายที่สัปดาห์ที่ 4 ได้ เงื่อนไขต่อท้ายนี้สำคัญ และหัวข้อท้ายโพสต์จะกลับมาที่มันอีกครั้ง เพราะ random intercept ที่ปล่อยให้ residual variance ต่างกันไปในแต่ละการนัด ก็ดัดเส้นทแยงมุมให้ไม่คงที่ได้เหมือนกัน ซึ่งเป็นเนื้อหาของตอนที่ 6

แถวตัวเลขข้างบนจึงประกาศตั้งแต่ก่อนจะมีพีชคณิตสักบรรทัดว่า มีบางอย่างถูกเติมเข้าไปใน $Z$ แล้ว

ลองอ่าน "พัด" นี้ดู ผู้เข้าร่วมจำลองแปดคนเริ่มต้นที่ความสูงต่างกันอย่างเห็นได้ชัด เพราะ SD ของ intercept ที่ 2.5 g·m⁻²·h⁻¹ กำหนดไว้อย่างนั้น จึงไม่มีจุดตั้งต้นร่วมกัน เส้นของแต่ละคนเอียงด้วยมุมที่ต่างกัน กลุ่มเส้นหุบเข้าหากันเล็กน้อยราวกลางช่วงติดตาม แล้วบานออกอีกครั้งเมื่อถึงสัปดาห์ที่ 8 เส้นที่วาดคือเส้นทางส่วนบุคคลที่ยังไม่มี noise คือ $b_{0i} + b_{1i}t$ รอบเส้นของประชากร ความกว้าง ของกลุ่มเส้นที่ตำแหน่งใด ๆ บนแกนนอนจึงเป็นส่วนเบี่ยงเบนมาตรฐาน ไม่ใช่ variance คือ 2.50, 2.41, 2.39 และ 2.53 g·m⁻²·h⁻¹ ที่สัปดาห์ที่ 0, 2, 4 และ 8 ทีนี้ยกกำลังสองสี่ค่านั้น แล้วบวก $\sigma^2 = 2.00$ ซึ่งเป็นสิ่งที่การวัดครั้งเดียวแบกไว้บนเส้นทางนั้น แถวในตารางข้างบนจะกลับมาพอดี คือ $6.25 + 2.00 = 8.25$, $5.81 + 2.00 = 7.81$, $5.69 + 2.00 = 7.69$ และ $6.41 + 2.00 = 8.41$ ทุกข้อความในโพสต์นี้อ่านออกมาจากภาพเดียวภาพนี้

ทำไมต้องเปลี่ยน design

ตำแหน่งผิวหนังสี่จุดในตัวอย่างหลักของซีรีส์เป็นชุดตำแหน่งวัดที่ตายตัวและไม่มีลำดับ (unordered) โดยแต่ละจุดมีค่าเฉลี่ยเป็นของตัวเอง มันไม่ได้เรียงกันเป็นลำดับ back ไม่ได้ "มาทีหลัง" shin และการเอา label ของตำแหน่งไปคูณกับ slope ก็ไม่ได้ให้ความหมายอะไร Random slope ต้องการมิติที่คนคนหนึ่ง "เคลื่อนที่" ไปตามนั้นได้ ซึ่งตำแหน่งผิวหนังสี่จุดไม่ได้ให้มิติแบบนั้น ตำแหน่งทั้งสี่จะกลับมาในตอนที่ 6 ซึ่งเป็นโครงสร้างที่เหมาะกับคำถามของตอนนั้นพอดี

โพสต์นี้จึงเปลี่ยนไปใช้ การศึกษาย่อยแบบนัดติดตาม (visit sub-study) อย่างชัดเจน และใช้เฉพาะในโพสต์นี้เท่านั้น นี่คือการศึกษาย่อยแบบสุ่มที่จำลองขึ้นภายในโครงการสมมติเดียวกัน โดยผู้เข้าร่วม 80 คนที่มีความบกพร่องของ skin barrier ในระดับ AD ถูกสุ่มแบบ 1:1 ไปยังกลุ่ม intensive emollient เทียบกับ standard care และวัด transepidermal water loss ที่ ท้องแขน (volar forearm) จุดเดียว ใน สี่การนัด คือสัปดาห์ที่ 0, 2, 4 และ 8 รวมเป็น 80 × 4 = 320 observations สมดุลเต็มรูปแบบ และเป็นภาพสะท้อนกลับด้านของ design เดิม คือวัดหนึ่งตำแหน่งสี่ครั้ง แทนที่จะวัดสี่ตำแหน่งครั้งเดียว

โมเดลที่ใช้สร้างข้อมูลคือ

\[ \text{tewl}_{ij} \;=\; \beta_0 + \beta_1 t_{ij} + \beta_{\text{arm}}\,\text{arm}_i + \beta_2\,(\text{arm}_i \times t_{ij}) \;+\; b_{0i} + b_{1i} t_{ij} \;+\; e_{ij} \]

โดยมีค่าจำลองดังนี้

พารามิเตอร์	ค่า	อ่านว่า
$\beta_0$ (ค่าเฉลี่ยสัปดาห์ที่ 0 ทั้งสองกลุ่ม)	18.0	จุดที่ทุกคนเริ่มต้น
$\beta_1$ (standard care ต่อสัปดาห์)	−0.60	เส้นทางของประชากรภายใต้ standard care
$\beta_{\text{arm}}$ (main effect ของ arm ที่สัปดาห์ที่ 0)	0	ถูกตั้งเป็นศูนย์ตอนสร้างข้อมูล ไม่ใช่ถูกตัดออกจากโมเดล
$\beta_2$ (ลดเพิ่มในกลุ่ม intensive ต่อสัปดาห์)	−0.25	treatment effect ซึ่งอยู่ใน interaction ทั้งหมด
$\tau_0^2$ (variance ของ intercept)	6.25 (SD 2.5)	คนแต่ละคนห่างกันแค่ไหนที่สัปดาห์ที่ 0
$\tau_1^2$ (variance ของ slope)	0.04 (SD 0.20 ต่อสัปดาห์)	อัตราการเปลี่ยนแปลงของแต่ละคนห่างกันแค่ไหน
$\text{corr}(b_{0i}, b_{1i})$	−0.30 → $\tau_{01} = -0.15$	ระดับตั้งต้นกับอัตราเปลี่ยนแปลงเดินทางไปด้วยกันอย่างไร
$\sigma^2$ (residual)	2.00	สิ่งที่เหลืออยู่ในการวัดครั้งเดียว

ส่วน fixed effect มีสองจุดที่ควรพูดถึงอย่างละหนึ่งประโยคก่อนจะผ่านไป จุดแรก main effect ของ arm ถูกตั้งไว้ที่ ศูนย์ ตอนสร้างข้อมูล ซึ่งไม่เหมือนกับการตัดมันออกจากโมเดล เหตุผลคือการสุ่มทำให้สองกลุ่มเท่ากันที่สัปดาห์ที่ 0 โดยค่าคาดหมาย ไม่ใช่ในกลุ่มตัวอย่างที่เกิดขึ้นจริง สูตร week * arm จึงขยายออกเป็น week + arm + week:arm และโมเดลที่ fit ก็ประมาณค่าพจน์นั้น แทนที่จะสมมติทิ้งไป ส่วน treatment effect เองอยู่ใน interaction กับเวลาทั้งหมด จุดที่สอง ภายใต้ standard care ประชากรจำลองลดลงจาก 18.0 เหลือ $18.0 - 8 \times 0.60 = 13.2$ ในแปดสัปดาห์ ขณะที่กลุ่ม intensive ลดด้วยอัตรา $-0.60 - 0.25 = -0.85$ ต่อสัปดาห์ และไปจบที่ $18.0 - 8 \times 0.85 = 11.2$

# random intercept + random slope on time (Post 5, visit sub-study)
fitS <- lmer(tewl ~ week * arm + (week | id), data = visits)

$G$ กลายเป็น matrix ขนาด 2×2

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

เมื่อมี random effect ตัวเดียว $G$ เป็นตัวเลขค่าเดียว พอมีสองตัว มันกลายเป็น symmetric matrix ขนาด 2×2 ที่บรรจุ variance สองค่าและ covariance หนึ่งค่า

\[ \begin{pmatrix} b_{0i}\\ b_{1i}\end{pmatrix} \sim N\!\left(\begin{pmatrix}0\\0\end{pmatrix},\; G=\begin{pmatrix} 6.25 & -0.15\\ -0.15 & 0.04 \end{pmatrix}\right) \]

เส้นทแยงมุมต้องอ่านในหน่วยของสิ่งที่มันไปคูณด้วย ค่า SD ของ intercept ที่ 2.5 g·m⁻²·h⁻¹ บอกว่าผู้เข้าร่วมที่อยู่สูงกว่าค่าเฉลี่ยหนึ่งส่วนเบี่ยงเบนมาตรฐาน เริ่มต้นด้วยผิวที่รั่วมากกว่า 2.5 หน่วยตั้งแต่วันสุ่ม ส่วน SD ของ slope ที่ 0.20 g·m⁻²·h⁻¹ ต่อสัปดาห์ บอกว่าผู้เข้าร่วมที่อยู่สูงกว่าค่าเฉลี่ยของ $b_{1i}$ หนึ่งส่วนเบี่ยงเบนมาตรฐาน มี อัตรา การเปลี่ยนแปลง สูงกว่า เส้นของประชากร 0.20 หน่วยต่อสัปดาห์ ภายใต้ standard care นั่นหมายถึงลดลง 0.40 แทนที่จะเป็น 0.60 หน่วยต่อสัปดาห์ คือดีขึ้น ช้ากว่า ไม่ใช่เร็วกว่า ส่วนต่างนั้นสะสมเป็น $8 \times 0.20 = 1.6$ หน่วยเมื่อถึงการนัดครั้งสุดท้าย ซึ่งมีขนาดพอ ๆ กับ treatment effect เองที่ $8 \times 0.25 = 2.0$ หน่วย การระวังเครื่องหมายตรงนี้คุ้มกับประโยคที่เพิ่มเข้ามา เพราะ slope เฉลี่ยเป็นค่าลบ คำว่า "เร็วกว่า" กับ "สูงกว่า" จึงชี้ไปคนละทาง

ค่านอกเส้นทแยงมุมคือของใหม่ มันคือ covariance และ covariance ที่มีหน่วยผสมกันแบบนี้รู้สึกตามได้ยาก จึงควรแปลงมันเสียก่อน คือ $-0.15 = -0.30 \times 2.5 \times 0.20$ หรือ correlation เท่ากับ −0.30 ระหว่างจุดที่คนคนหนึ่งเริ่มต้น กับความเร็วที่คนคนนั้นเปลี่ยนแปลง ในชุดข้อมูล จำลอง ชุดนี้ ผู้ที่เริ่มต้นด้วยผิวที่รั่วมากกว่ามีแนวโน้มจะดีขึ้นเร็วกว่าเล็กน้อย บนภาพพัด เส้นที่เริ่มสูงที่สุดคือเส้นที่ชันที่สุดด้วย ซึ่งเป็นเหตุผลที่กลุ่มเส้นหุบเข้าก่อนที่จะบานออกอีกครั้ง

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

ผลลัพธ์ข้อแรก: variance กลายเป็นฟังก์ชันของเวลา

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

เพราะตอนนี้ $Y$ มี $b_{0i} + b_{1i}t$ อยู่ข้างใน marginal variance จึงรับมรดกทั้ง random effect ทั้งสองตัวและ covariance ระหว่างกันมาด้วย

\[ \operatorname{Var}(Y_{ij}\mid t) \;=\; \tau_0^2 \;+\; 2t\,\tau_{01} \;+\; t^2\tau_1^2 \;+\; \sigma^2 \]

นี่คือฟังก์ชันกำลังสองใน $t$ และควรลงมือคำนวณให้ครบอย่างน้อยสักครั้ง ที่สัปดาห์ที่ 0 พจน์ของ slope ทั้งสองหายไป เหลือเพียง variance ของ intercept กับ residual คือ $6.25 + 0 + 0 + 2.00 = 8.25$ ที่สัปดาห์ที่ 2 พจน์ผลคูณไขว้เริ่มติดลบขณะที่พจน์กำลังสองยังเล็กมาก คือ $6.25 + 2(2)(-0.15) + (2^2)(0.04) + 2.00 = 6.25 - 0.60 + 0.16 + 2.00 = 7.81$ ที่สัปดาห์ที่ 4 ทั้งสองพจน์ยังหักล้างกันต่อไป คือ $6.25 - 1.20 + 0.64 + 2.00 = 7.69$ และที่สัปดาห์ที่ 8 พจน์กำลังสองเอาชนะได้ในที่สุด เพราะมันโตตาม $t^2$ ขณะที่พจน์ที่ติดลบโตตาม $t$ เท่านั้น คือ $6.25 - 2.40 + 2.56 + 2.00 = 8.41$

เส้นโค้งนี้ต่ำสุดตรงที่อนุพันธ์ $2\tau_{01} + 2t\tau_1^2$ เท่ากับศูนย์ นั่นคือที่ $t^{*} = -\tau_{01}/\tau_1^2 = 0.15/0.04 = 3.75$ สัปดาห์ ซึ่ง variance เท่ากับ 7.69 พัดจึงแคบที่สุดก่อนถึงการนัดครั้งที่สามเล็กน้อย และกว้างที่สุดที่การนัดครั้งสุดท้าย

จากตรงนี้มีสองข้อสรุปที่คนละชนิดกัน ในทางคณิตศาสตร์ เส้นทแยงมุมของ $V$ ไม่คงที่อีกต่อไป และไม่มีความประณีตในการประมาณค่าใด ๆ ที่จะทำให้มันกลับมาคงที่ได้ตราบใดที่ $\tau_1^2 > 0$ ในทางปฏิบัติ variance function รูปแบบนี้อาจเป็นเหตุผลให้คุณลองดูเส้นทแยงมุมที่โมเดลของคุณเองบอกเป็นนัย แทนที่จะดูแต่ตาราง fixed effect เพราะโมเดลที่ถูกบังคับให้ variance คงที่ต้องไปประนีประนอมที่ใดที่หนึ่ง และจุดที่มันประนีประนอมมักอยู่ที่ปลายทั้งสองข้างของช่วงติดตาม

ถ้าปิดสวิตช์ variance ของ slope ทิ้ง คือให้ $\tau_1^2 = 0$ และดังนั้น $\tau_{01} = 0$ สมการเดิมจะคืนค่า $6.25 + 2.00 = 8.25$ ที่ทุกการนัด นั่นคือคำตอบของ random intercept และมันคือค่าที่สัปดาห์ที่ 0 พอดี ซึ่งอธิบายว่าทำไมโมเดลที่มีแต่ intercept จึงมักดูมีเหตุผลที่สุดตรงจุดที่ข้อมูลบอกอะไรเกี่ยวกับ slope ได้น้อยที่สุด

ผลลัพธ์ข้อที่สอง: ค่านอกเส้นทแยงมุมก็ขยับด้วย

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

covariance ระหว่างสองการนัดในคนคนเดียวกันได้มาด้วยวิธีเดียวกัน คือรวบรวมพจน์ที่การวัดทั้งสองครั้งใช้ร่วมกัน

\[ \operatorname{Cov}(Y_{ij}, Y_{ik}) \;=\; \tau_0^2 \;+\; (t_j+t_k)\,\tau_{01} \;+\; t_j t_k\,\tau_1^2 \]

ลองคิดคู่ที่ห่างที่สุดก่อน สำหรับสัปดาห์ที่ 0 กับ 8 ค่า $t_j t_k = 0$ พจน์กำลังสองจึงหายไปทั้งหมด เหลือเพียงพจน์ไขว้ที่ทำงาน คือ $6.25 + (0+8)(-0.15) + 0 = 5.05$ หารด้วยส่วนเบี่ยงเบนมาตรฐานสองตัวที่เกี่ยวข้องเพื่อให้ได้ correlation คือ $5.05/\sqrt{8.25 \times 8.41} = 5.05/8.33 = 0.61$ จากนั้นลองคู่ที่ติดกัน คือสัปดาห์ที่ 2 กับ 4 ได้ $6.25 + (6)(-0.15) + (8)(0.04) = 6.25 - 0.90 + 0.32 = 5.67$ และ $5.67/\sqrt{7.81 \times 7.69} = 5.67/7.75 = 0.73$ ครบทั้งหกคู่พร้อม correlation ที่ตามมาเป็นดังนี้

คู่ (สัปดาห์)	0,2	0,4	0,8	2,4	2,8	4,8
covariance	5.95	5.65	5.05	5.67	5.39	5.73
correlation	0.74	0.71	0.61	0.73	0.67	0.71

เอาไปวางเทียบกับ matrix ของตอนที่ 2 ที่ทุกช่องนอกเส้นทแยงมุมเป็น 6.25 และทุก correlation เป็น 0.625 จะเห็นว่าไม่ว่าเส้นทแยงมุมหรือค่านอกเส้นทแยงมุมก็ไม่คงที่อีกต่อไปแล้ว compound symmetry ตรงนี้ไม่ได้ถูกละเมิดแบบคร่าว ๆ แต่เป็นเท็จอย่างเป็นทางการโดยโครงสร้าง ตั้งแต่วินาทีที่ random slope เข้ามา ถ้าปิดพจน์ของ slope ทิ้ง correlation ทั้งหกจะยุบลงเหลือค่าเดียวคือ $6.25/8.25 = 0.76$ ซึ่งก็คือโครงสร้างแบนราบที่ตอนก่อน ๆ อธิบายไว้

และควรสังเกตด้วยว่ารูปแบบนี้ ไม่ใช่ อะไร ในการจำลองชุดนี้ correlation ทั้งหกค่าลดลงจริงเมื่อการนัดสองครั้งห่างกันมากขึ้น คือ 0.74 ที่ระยะห่างสั้นที่สุด และ 0.61 ที่ระยะห่างยาวที่สุด อย่างไรก็ตาม มันไม่ได้เป็น ฟังก์ชัน ของระยะห่างนั้น และคุณสมบัตินี้ต่างหากที่สำคัญ สัปดาห์ที่ 0 กับ 2 ให้ correlation 0.74 ขณะที่สัปดาห์ที่ 2 กับ 4 ซึ่งห่างกันสองสัปดาห์เท่ากันเป๊ะ ให้ 0.73 ส่วนคู่ที่ห่างกันสี่สัปดาห์ทั้งสองคู่แยกจากกันที่ทศนิยมตำแหน่งที่สาม คือ 0.709 สำหรับสัปดาห์ที่ 0 กับ 4 เทียบกับ 0.713 สำหรับสัปดาห์ที่ 4 กับ 8 และ covariance ของสองคู่นั้นแยกกันให้เห็นชัดที่ 5.65 กับ 5.73 โครงสร้างแบบ autoregressive ทำอย่างนั้นไม่ได้ เพราะ correlation ของมันคือ $\rho^{|t_j - t_k|}$ ซึ่งรู้จักแค่ระยะห่างระหว่างการนัดสองครั้งเท่านั้น คู่ที่ห่างเท่ากันทุกคู่จึงต้องให้ค่าเดียวกันเสมอ Random slope จึงไม่ใช่โครงสร้างแบบ autoregressive ที่เปลี่ยนชื่อมาใส่ มันสร้างรูปแบบกำลังสองที่รูปร่างถูกกำหนดโดย $G$ และโดย ตำแหน่งของการนัดบนแกนเวลา ไม่ใช่โดยระยะห่างเพียงอย่างเดียว

เขียนส่วน random ได้สองแบบ

สูตร (week | id) ขอทั้ง random effect สองตัว และ covariance ระหว่างกัน รวมเป็น variance สองค่ากับ covariance หนึ่งค่า คือสามพารามิเตอร์ใน $G$ ส่วนทางเลือก (1 | id) + (0 + week | id) ขอ variance สองค่าเหมือนกัน แต่ตรึง covariance ไว้ที่ศูนย์ จึงใช้สองพารามิเตอร์แทนที่จะเป็นสาม

# correlated random intercept and slope — three parameters in G
fitS <- lmer(tewl ~ week * arm + (week | id), data = visits)

# independent random intercept and slope — two parameters in G
fitS_ind <- lmer(tewl ~ week * arm + (1 | id) + (0 + week | id), data = visits)

specification แบบที่สองไม่ใช่แค่รุ่นย่อส่วนของแบบแรก แต่เป็นการอ้างคนละอย่าง มันยืนยันว่า ณ จุดกำเนิดที่คุณเลือกไว้ ระดับตั้งต้นของคนคนหนึ่งไม่ได้บอกอะไรเกี่ยวกับอัตราการเปลี่ยนแปลงของคนคนนั้นเลย ในการจำลองชุดนี้ ข้ออ้างดังกล่าวเป็นเท็จที่สัปดาห์ที่ 0 เพราะ correlation ถูกกำหนดไว้ที่ −0.30 โดยโครงสร้าง อย่างไรก็ตาม รูปแบบ independent อาจเป็นทางเลือกเชิงปฏิบัติที่สมเหตุสมผลเมื่อรูปแบบ correlated ไม่ converge โดยมีเงื่อนไขว่าการลดรูปนั้นต้องถูกรายงาน ไม่ใช่ถูกใช้อย่างเงียบ ๆ โครงสร้าง random แบบสองคอลัมน์นี้เขียนได้เช่นกันผ่าน random-effects formula แบบสองคอลัมน์ใน nlme และผ่าน random-effects specification ที่สอดคล้องกันใน mixed ของ Stata syntax ต่างกัน แต่จำนวนพารามิเตอร์เท่ากัน

จุดที่คุณวางให้ week เป็นศูนย์ เปลี่ยนความหมายของ $G$

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

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

random intercept คือระยะที่คนคนหนึ่งเบี่ยงจากเส้นของประชากร ที่ $t = 0$ ไม่ว่าคุณจะวาง $t = 0$ ไว้ตรงไหนก็ตาม พอย้ายจุดกำเนิด คุณก็กำลังถามถึงปริมาณคนละตัว ทั้ง variance ของ intercept และ covariance ระหว่าง intercept กับ slope จึงเปลี่ยนตามไปด้วย การย้ายจุดอ้างอิงไปที่เวลา $c$ ให้ intercept ตัวใหม่คือ $b_{0i}^{*} = b_{0i} + c\,b_{1i}$ และดังนั้น

\[ \operatorname{Var}(b_{0i}^{*}) = \tau_0^2 + 2c\,\tau_{01} + c^2\tau_1^2, \qquad \operatorname{Cov}(b_{0i}^{*}, b_{1i}) = \tau_{01} + c\,\tau_1^2 \]

ถ้าย้ายไปอ้างอิงที่สัปดาห์ที่ 4 ซึ่งเป็นกลางช่วงติดตาม covariance จะกลายเป็น $-0.15 + 4(0.04) = +0.01$ โดยมี variance ของ intercept เท่ากับ $6.25 - 1.20 + 0.64 = 5.69$ และได้ correlation เท่ากับ $0.01/\sqrt{5.69 \times 0.04} = 0.02$ ข้อมูลจำลองชุดเดิม โมเดลตัวเดิม ค่าประมาณ fixed effect ของ interaction กับการรักษาตัวเดิม แต่ correlation ระหว่าง intercept กับ slope เดินจาก −0.30 ไปเป็นแทบเท่ากับศูนย์ เพียงเพราะจุดกำเนิดขยับไปสี่สัปดาห์

covariance จะเป็นศูนย์พอดีที่ $c = -\tau_{01}/\tau_1^2 = 3.75$ สัปดาห์ ซึ่งเป็น 3.75 ตัวเดียวกับที่ทำให้ variance function ต่ำสุด นี่ไม่ใช่ความบังเอิญ และพีชคณิตสองบรรทัดที่อยู่บนหน้านี้แล้วก็พิสูจน์มันได้ คืออนุพันธ์ของ variance function ที่เป็น $2\tau_{01} + 2t\,\tau_1^2$ เท่ากับสองเท่าของ covariance หลังย้ายจุดอ้างอิงที่เป็น $\tau_{01} + c\,\tau_1^2$ พอดี ค่าหนึ่งจึงเป็นศูนย์ก็ต่อเมื่ออีกค่าหนึ่งเป็นศูนย์ เอกลักษณ์นี้เป็นจริงกับ $G$ ใด ๆ ที่มี $\tau_1^2 > 0$ ไม่ใช่เฉพาะกับตัวเลขชุดนี้ จุดกำเนิดที่ทำให้ระดับตั้งต้นกับอัตราการเปลี่ยนแปลงไม่สัมพันธ์กัน คือจุดกำเนิดที่พัดแคบที่สุดพอดี

ดังนั้นคำแนะนำเชิงปฏิบัติจึงไม่ใช่แค่ "ให้ centring เวลา" แต่แคบกว่านั้น คือให้วาง $t = 0$ ไว้ที่ช่วงเวลาที่คุณต้องการให้ intercept อธิบายจริง ๆ เช่น วันสุ่ม วันเริ่มยา หรือวันเริ่มมีอาการ แล้วตีความ correlation ระหว่าง intercept กับ slope ในฐานะข้อความเกี่ยวกับ ช่วงเวลานั้น เท่านั้น การรายงาน correlation ระหว่าง intercept กับ slope โดยไม่ระบุจุดอ้างอิง แทบจะตีความไม่ได้เลย

เมื่อ random slope ประมาณค่าไม่ได้

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

Random slope คือ variance component และ variance component ต้องการข้อมูลจึงจะถูก identify ได้ มีสามสถานการณ์ที่พบซ้ำ ๆ

ข้อแรก จำนวนการวัดซ้ำต่อคนน้อยเกินไป การประมาณ slope เฉพาะบุคคลต้องมีจุดวัดภายในคนมากพอที่จะแยกแนวโน้มส่วนตัวออกจาก residual noise เมื่อมีการนัดเพียงสองครั้งที่เวลาเดียวกันสำหรับทุกคน อุปสรรคไม่ใช่ความยาก แต่เป็นการนับ matrix ของ covariance ขนาด $2 \times 2$ บรรจุปริมาณที่ต่างกันเพียงสามค่า คือ variance สองค่ากับ covariance หนึ่งค่า ขณะที่โมเดล random slope แบบ correlated ขอถึงสี่ค่า คือ $\tau_0^2$, $\tau_{01}$, $\tau_1^2$ และ $\sigma^2$ สี่พารามิเตอร์ดึงกลับมาจากสามปริมาณไม่ได้ โมเดลจึงไม่ถูก identify เลย และไม่มีขนาดตัวอย่างเท่าใดที่ซ่อมได้ นี่คือการนับแบบเดียวกับที่ตอนที่ 6 จะทำอย่างละเอียดภายใต้หัวข้อเพดานของ identifiability การนัดตั้งแต่สามครั้งขึ้นไปทำให้ identify ได้ในหลักการ และเมื่อจำนวนผู้เข้าร่วมมีไม่มาก ค่าประมาณก็ยังอาจไม่เสถียรแม้ในทางเทคนิคจะ identify ได้แล้วก็ตาม

ข้อที่สอง ตัวแปรที่จะให้ slope ไม่มีความแปรผันภายในคน กลุ่มการรักษาที่สุ่มในการศึกษาย่อยนี้คงที่ภายในคนคนเดียวกัน การใส่ random effect ของ arm ภายใน id จึงประมาณค่าไม่ได้เลย มีเพียงตัวแปรที่ขยับอยู่ภายในคน ซึ่งในที่นี้คือเวลา เท่านั้นที่แบก random slope เหนือคนคนนั้นได้

ข้อที่สาม boundary fit และ singular fit ค่า likelihood อาจสูงสุดที่ขอบของ parameter space และ output ก็ประกาศขอบนั้นได้มากกว่าหนึ่งแบบ คุณอาจเห็น variance ของ slope ที่ถูกประมาณเป็นศูนย์หรือใกล้ศูนย์มาก เห็น correlation ระหว่าง intercept กับ slope ที่เป็น −1.00 หรือ +1.00 พอดี หรือเห็นข้อความ boundary (singular) fit ของ lme4 (Bates และคณะ, 2015) สองอย่างแรกไม่ใช่สถานการณ์เดียวกัน variance ของ slope ที่เป็นศูนย์คือการล้างมิติของ slope ทิ้งทั้งมิติ เพราะ design ไม่มีข้อมูลเรื่องความต่างของอัตราระหว่างคนให้ดึงกลับมาได้เลย ส่วน correlation ที่เป็น ±1.00 พอดีบอกอีกอย่างหนึ่ง คือ variance ทั้งสองค่าถูกประมาณเป็นบวก แต่ $G$ ยุบลงเหลือ rank เดียว random effect สองตัวจึงกลายเป็นมิติเดียวที่ปลอมตัวมา singular fit ไม่ใช่ความล้มเหลวของซอฟต์แวร์ และไม่ใช่เหตุผลอัตโนมัติที่จะเก็บโมเดลนั้นไว้ ไม่ว่าจะเป็นแบบใด มันคือข้อมูลกำลังรายงานว่าโครงสร้าง random ที่คุณขอไปนั้นซับซ้อนเกินกว่าที่ design จะรองรับได้

กฎเชิงปฏิบัติ ซึ่งขอกล่าวอย่างระมัดระวังตามที่มันควรเป็น คือเมื่อจำนวนการนัดต่อคนมีน้อย โครงสร้าง random slope แบบ correlated อาจมากเกินกว่าที่ข้อมูลจะ identify ได้ และการประกาศเส้นทางการลดรูปไว้ล่วงหน้า คือจาก correlated slope ไปเป็น independent slope และไปเป็น intercept อย่างเดียว น่าจะเป็นสิ่งที่ควรระบุไว้ในแผนการวิเคราะห์ก่อนที่จะ fit โมเดลตัวแรก

ควรใส่ random structure มากแค่ไหน

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

จุดยืนแบบ maximal (Barr และคณะ, 2013) เห็นว่าสำหรับการทดสอบสมมติฐานเชิงยืนยัน โครงสร้าง random effect ควรบรรจุ random slope ทุกตัวที่ design รองรับ เหตุผลคือการละ random slope ของพจน์ที่แปรผันจริงระหว่างหน่วย อาจทำให้ Type I error ของพจน์นั้นสูงเกินจริงอย่างรุนแรง ส่วนจุดยืนแบบ parsimony (Matuschek และคณะ, 2017) ยอมรับข้อโต้แย้งนั้นในหลักการ แต่แสดงให้เห็นว่าโครงสร้างแบบ maximal ที่ถูก fit กับข้อมูลที่รองรับไม่ไหว จะสูญเสีย power และให้ผลลัพธ์ที่ไม่เสถียร และโครงสร้างที่เลือกโดยคำนึงถึงดุลระหว่าง Type I error กับ power ทำงานได้ดีกว่าในสถานการณ์จริง Pinheiro และ Bates (2000) ให้กรอบพื้นฐานของทั้งสองจุดยืนไว้ว่า covariance model ของ random effect คือชุดพารามิเตอร์ที่ต้องถูก identify ไม่ใช่รายการความปรารถนา

ข้อเสนอแนะที่เคารพทั้งสองฝ่าย ซึ่งขอเสนอในฐานะความเห็นมากกว่ากฎ คือโครงสร้าง random น่าจะถูกกำหนดล่วงหน้าจาก design และคำถามวิจัย โดยคงไว้ซึ่ง random slope ของตัวแปรที่สมมติฐานหลักถามถึง เมื่อจำนวนการวัดซ้ำต่อคนพอจะรองรับได้ และลดรูปตามเส้นทางที่ประกาศไว้เมื่อรองรับไม่ไหว สิ่งที่ไม่ควรเกิดขึ้นคือการเลือกโครงสร้าง random จากการดูค่า p ของ fixed effect เพราะนั่นทำให้ผลการอนุมานที่รายงานขึ้นกับสิ่งที่มันควรจะทดสอบเสียเอง

เปรียบเทียบ random structure โดยไม่หลอกตัวเอง

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

การเปรียบเทียบนี้ถูกกำกับด้วยข้อเท็จจริงสามข้อ ข้อแรก การทำ likelihood-ratio test ระหว่างโครงสร้าง random ที่ซ้อนกัน ต้อง fit ด้วย fixed effect ชุดเดียวกัน เพราะค่า likelihood แบบ REML เปรียบเทียบข้าม specification ของ fixed effect ที่ต่างกันไม่ได้ ข้อที่สอง การทดสอบว่า variance component ตัวเดียว เป็นศูนย์หรือไม่ วางค่าภายใต้สมมติฐานหลักไว้ที่ขอบของ parameter space การแจกแจงอ้างอิงที่ถูกต้องจึงเป็นส่วนผสม 50:50 ของ $\chi^2_0$ กับ $\chi^2_1$ ไม่ใช่ $\chi^2_1$ เพียว ๆ การเทียบค่าสถิติกับการแจกแจงไคสแควร์ตามแบบเดิมจึงเป็นการทดสอบที่อนุรักษ์นิยม และค่า p ที่ได้จะใหญ่เกินจริงมากกว่าเล็กเกินจริง ผลลัพธ์ที่สะอาดนี้เป็นของ variance component ทีละตัวเท่านั้น และไม่ควรถูกเหมารวม เพราะการตัด random slope ออกพร้อมกับ covariance ของมันคือการตัดพารามิเตอร์สองตัวที่มีเพียงตัวเดียวติดขอบ การแจกแจงอ้างอิงที่ถูกต้องจึงต้องคำนวณออกมา ไม่ใช่สมมติเอา ข้อที่สาม AIC และ BIC เป็นหลักฐานสนับสนุนเรื่องโครงสร้าง ไม่ใช่หลักฐานชี้ขาด และไม่มีตัวใดใช้แทนการดู residual และเส้นทางที่สังเกตได้จริง

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

Random slope อธิบายว่า correlation มาจากไหน แต่ไม่ได้รับรองว่าคำอธิบายนั้นถูกต้อง

ข้อแรก variance ที่ไม่คงที่ในข้อมูลจริงไม่ได้พิสูจน์ว่าต้องมี random slope โมเดลที่มี random intercept คงที่แต่ให้ residual variance ต่างกันตามครั้งที่วัด สามารถสร้างเส้นทแยงมุมที่หน้าตาใกล้เคียงกันมากได้ และ specification ทั้งสองแยกจากกันไม่ได้ด้วยการมองกราฟ variance เฉย ๆ ซึ่งตอนที่ 6 จะหยิบคำถามนี้ขึ้นมาตรง ๆ

ข้อที่สอง correlation ระหว่าง intercept กับ slope ขึ้นกับจุดกำเนิด ดังที่การคำนวณเรื่อง centring ข้างต้นแสดงให้เห็น จึงไม่ควรรายงานมันเป็นข้อค้นพบทางชีววิทยาเรื่องการลู่เข้าหรือ regression to the mean โดยไม่ระบุจุดอ้างอิง

ข้อที่สาม ไม่มีอะไรในโพสต์นี้ที่เปลี่ยนความสัมพันธ์ให้กลายเป็นข้อความเชิงสาเหตุ treatment effect ในการศึกษาย่อยจำลองนี้ตีความเป็นการเปรียบเทียบแบบสุ่มได้ ก็เพราะการศึกษาย่อยนี้ถูกจำลองให้เป็นการสุ่มเท่านั้น ในข้อมูล cohort เชิงสังเกต โมเดล random slope ตัวเดียวกันจะยังแบกปัญหา confounding ทุกอย่างที่มันมีอยู่ตั้งแต่ต้น

ข้อที่สี่ ค่า slope เฉพาะบุคคลที่ mixed model ทำนายออกมาเป็นปริมาณที่ถูก shrink แล้ว ไม่ใช่ค่าที่วัดได้ การดึงมันออกมาวิเคราะห์ต่อเสมือนเป็นอัตราการเปลี่ยนแปลงที่สังเกตได้จริง คือการวิเคราะห์ขั้นที่สองที่ละเลยความไม่แน่นอนของมัน และอาจให้ข้อสรุปที่มั่นใจเกินจริงอย่างมาก ในกรณีที่ intercept และ slope เฉพาะบุคคลจำเป็นต้องถูกนำไปใช้ที่อื่นจริง ๆ เช่นใช้เป็นตัวทำนายเหตุการณ์ทางคลินิกในภายหลัง โครงสร้างที่เหมาะสมคือ joint model ที่ประมาณทั้งสองส่วนไปพร้อมกัน ซึ่งเป็นเนื้อหาของ Joint model สำหรับข้อมูล longitudinal และ survival

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

1. พล็อตเส้นทางของผู้เข้าร่วมก่อน fit อะไรทั้งสิ้น ถ้าเส้นที่สังเกตได้ไม่ขนานกันอย่างเห็นได้ชัด random slope ก็ควรค่าแก่การพิจารณา แต่ถ้าเส้นขนานกันและสั้น มันอาจ identify ไม่ได้ ไม่ว่าในทางทฤษฎีจะน่าสนใจแค่ไหน
2. วางจุดกำเนิดของเวลาอย่างตั้งใจ ให้ $t = 0$ อยู่ที่ช่วงเวลาที่คุณต้องการให้ intercept อธิบาย และระบุจุดอ้างอิงนั้นทุกครั้งที่รายงาน $G$
3. กำหนดโครงสร้าง random และเส้นทางการลดรูปไว้ล่วงหน้า การประกาศทางถอยไว้ก่อน fit ตัวแรกอาจช่วยกันไม่ให้โครงสร้างถูกเลือกจากเวอร์ชันที่ให้ค่า p ที่ถูกใจกว่า
4. รายงาน $G$ ไม่ใช่แค่ $\beta$ variance สองค่า covariance และ correlation ที่ตามมาคือผลการศึกษา ในงาน longitudinal มันอาจให้ข้อมูลกับผู้อ่านมากกว่าทศนิยมตำแหน่งที่สามของ fixed effect
5. ถือว่า singular fit คือข้อมูล รายงานมัน รายงานว่าคุณลดรูปไปเป็นอะไร และหลีกเลี่ยงการนำเสนอโมเดลที่ลดรูปแล้วราวกับว่ามันเป็นแผนตั้งแต่แรก
6. เปรียบเทียบโครงสร้าง random ที่ซ้อนกันด้วย REML บน fixed effect ชุดเดียวกัน และอย่าลืมว่าการทดสอบ variance component ตัวเดียว ที่ขอบเป็นการทดสอบแบบอนุรักษ์นิยม ก่อนจะสรุปว่า variance ของ slope "ไม่มีนัยสำคัญ"

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

• Random slope เติมคอลัมน์ที่สองให้ $Z$ ทำให้แต่ละคนพกอัตราการเปลี่ยนแปลงส่วนตัวมาเพิ่มจากค่าชดเชยส่วนตัว และ $G$ โตจากตัวเลขค่าเดียวเป็น matrix 2×2 ที่บรรจุ 6.25, 0.04 และ −0.15 ในการจำลองนี้
• marginal variance กลายเป็นฟังก์ชันกำลังสองของเวลา คือ 8.25, 7.81, 7.69 และ 8.41 ที่สัปดาห์ที่ 0, 2, 4 และ 8 โดยต่ำสุดที่สัปดาห์ที่ 3.75 เส้นทแยงมุมของ $V$ จึงไม่คงที่อีกต่อไป
• ค่านอกเส้นทแยงมุมก็ไม่คงที่เช่นกัน correlation จำลองทั้งหกค่าอยู่ระหว่าง 0.61 ถึง 0.74 และไม่ได้เป็นฟังก์ชันของระยะห่างเพียงอย่างเดียว การนัดสองครั้งที่ห่างกันสองสัปดาห์ให้ 0.74 ที่ต้นช่วงติดตาม แต่ให้ 0.73 ที่ช่วงถัดมา ซึ่งฆ่า compound symmetry อย่างเป็นทางการ และแยก random slope ออกจากโครงสร้างแบบ autoregressive
• correlation ระหว่าง intercept กับ slope ที่ −0.30 เป็นข้อความเกี่ยวกับสัปดาห์ที่ 0 เท่านั้น การย้ายจุดอ้างอิงไปที่สัปดาห์ที่ 4 เปลี่ยนมันเป็น 0.02 โดยไม่มีข้อมูลสักค่าเดียวเปลี่ยนไป
• Random slope ต้องการการวัดซ้ำต่อคนที่มากพอ และตัวแปรที่แปรผันจริงภายในคน ส่วน singular fit คือข้อมูลที่กำลังปฏิเสธคำขอ ไม่ใช่ซอฟต์แวร์ที่กำลังพัง
• การตัดสินว่าจะใส่ random structure มากแค่ไหน เป็นดุลยพินิจที่ปกป้องได้ระหว่างจุดยืน maximal กับ parsimony และควรถูกกำหนดล่วงหน้าจาก design มากกว่าจะถูกเลือกจากค่า p ของ fixed effect

ตอนต่อไป ในตอนที่ 6 เราจะปล่อย $Z$ ไว้อย่างนั้น แล้วหันไปจัดการกับ $R$ และจะพบว่า random intercept ให้ covariance ค่าเดียวกันกับทุกคู่ของการวัด แต่ไม่ได้บังคับให้ correlation เท่ากันด้วยตัวมันเอง

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

1. Pinheiro JC, Bates DM. Mixed-Effects Models in S and S-PLUS. Springer; 2000.
2. Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. J Stat Softw. 2015;67(1):1–48.
3. Barr DJ, Levy R, Scheepers C, Tily HJ. Random effects structure for confirmatory hypothesis testing: keep it maximal. J Mem Lang. 2013;68(3):255–278.
4. Matuschek H, Kliegl R, Vasishth S, Baayen H, Bates D. Balancing Type I error and power in linear mixed models. J Mem Lang. 2017;94:305–315.