SE เท่ากันสี่ค่า: วิธีปกป้อง mixed model ต่อหน้าอาจารย์

บทคัดย่อ

Mixed model ที่มี repeated factor พิมพ์ SE ค่าเดียวกันทุกระดับ แล้วมีคนถามว่าโปรแกรมพังไหม คำตอบคือไม่ ภายใต้ดีไซน์สมดุลและ residual variance ค่าเดียวที่ pooled ไว้ SE ของ contrast แต่ละหน่วยเท่ากับรากที่สองของ 2 เท่าของ variance รวมหารด้วย n ซึ่งไม่มี index ของหน่วย ค่าที่เท่ากันจึงเป็นสิ่งที่พีชคณิตบังคับ ตอนปิดซีรีส์พิสูจน์ประโยคนั้น ระบุส่วนประกอบสี่อย่างที่พึ่งพา แล้วจงใจทำลายทีละข้อ เช่น เมื่อค่าวัดตำแหน่งหนึ่งหาย 12 คน SE ขยับจาก 0.71 เป็นราว 0.74 และปฏิเสธความผิดพลาดด้านตรงข้าม เพราะ output ที่ถูกต้องไม่เคยยืนยันสมมติฐานตัวเอง ผู้อ่านจะได้สคริปต์สี่ขั้นสำหรับพูดปกป้อง และ checklist รายงานผล

---

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

เรามาถึงไหนแล้ว. ซีรีส์นี้ว่าด้วยการศึกษาที่ผู้เข้าร่วมหนึ่งคนให้ข้อมูลหลายแถว ในที่นี้คือผิวสี่ตำแหน่งที่วัดในคนเดียวกันทั้ง 80 คน ทำให้แถวเหล่านั้นถูกมองว่าเป็นอิสระต่อกันไม่ได้ วิธีแก้ที่ใช้มาตลอดคือ mixed model ซึ่งให้ค่าออฟเซ็ตประจำตัวคนละหนึ่งค่าคือ $b_i$ ใส่ซ้ำเหมือนกันในทุกแถวของคนคนนั้น และค่าตัวเดียวนี้เองที่ทำให้ค่าวัดของคนคนเดียวกันเดินทางไปด้วยกัน เจ็ดตอนที่ผ่านมาเราตามค่านี้ไปเรื่อย ๆ เข้าไปใน covariance matrix เข้าไปเป็น intraclass correlation เข้าไปอยู่ในสมการทั่วไป $V_i = Z_i G Z_i' + R_i$ และเข้าไปอยู่ในสองจุดที่ mixed model ทุกตัวถูกขยายได้ ตอนที่แล้ว ทางเลือกแบบ Population-Averaged: GLS, GEE และจุดที่โชคของ linear model หมดลง แสดงให้เห็นว่าเกิดอะไรขึ้นเมื่อโมเดลปฏิเสธที่จะอธิบายว่า correlation มาจากไหน แล้วบอกตรง ๆ ไปเลยว่ามันเป็นเท่าไร ตอนปิดซีรีส์นี้จะแปลงทั้งซีรีส์ให้เหลือหนึ่งย่อหน้าที่คุณพูดปกป้องได้จริงต่อหน้าคนที่ยังไม่เชื่อ และมันเริ่มต้นแบบที่มักจะเริ่มเสมอ คือด้วยคิ้วที่ยกขึ้นข้างหนึ่ง

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

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

ฉากเปิด

ผลลัพธ์อยู่บนจอแล้ว โมเดลคือตัวเดียวกับที่ซีรีส์นี้ fit มาตั้งแต่ตอนที่ 1 คือค่า transepidermal water loss (TEWL, g·m⁻²·h⁻¹) ที่วัดจากผิวหนังบริเวณที่ไม่มีรอยโรค 4 ตำแหน่งมาตรฐาน ในผู้ใหญ่จำลอง 80 คน แบ่งเป็นผู้ที่มี moderate atopic dermatitis 40 คนและผู้ที่ไม่มีโรคภูมิแพ้ 40 คน โดยมี interaction แบบ site × group มี covariate ระดับบุคคลสามตัว และมี random intercept หนึ่งค่าต่อผู้เข้าร่วมหนึ่งคน สี่ตำแหน่ง สี่ค่าประมาณของผลต่าง AD − control และคอลัมน์ standard error ที่อ่านได้เป็นตัวเลขเดียวกันสี่ครั้ง

site	AD − control (g·m⁻²·h⁻¹)	SE (g·m⁻²·h⁻¹)
forearm	6.0	0.71
hand	8.0	0.71
shin	3.0	0.71
back	5.0	0.71

ค่าประมาณของผลต่างไม่ได้เท่ากัน ข้อมูลจำลองชุดนี้ถูกสร้างขึ้นให้มี site × group interaction จริง และตัวเลข 6.0, 8.0, 3.0 และ 5.0 ก็สะท้อนสิ่งนั้น มีเพียงคอลัมน์ของความไม่แน่นอนเท่านั้นที่ซ้ำกัน คำถามที่ตามมามักเป็นคำถามเดียวกันเสมอ และเป็นคำถามที่ยุติธรรมมาก คือ โปรแกรมก๊อปตัวเลขไล่ลงมาทั้งคอลัมน์หรือเปล่า

คำตอบสั้น ๆ ก่อนจะมีพีชคณิต

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

เริ่มจากเหตุผลว่าทำไมคอลัมน์นี้ถึงดูผิด ในตารางเกือบทุกตารางที่แพทย์อ่าน แต่ละแถวมี standard error ของตัวเอง เพราะแต่ละแถวตั้งอยู่บนจำนวนผู้ป่วยที่ต่างกันและการกระจายที่ต่างกัน ค่าประมาณจากกลุ่มย่อย 30 คนย่อมแม่นยำน้อยกว่าค่าประมาณจากกลุ่มย่อย 300 คน ตัวเลขสี่ค่าที่เท่ากันเป๊ะในคอลัมน์ความไม่แน่นอนจึงดูเหมือนอุบัติเหตุจากการก๊อปวาง

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

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

นี่ยังเป็นเหตุผลว่าทำไมการซ้ำกันนี้จึงไม่หายไปเมื่อเปลี่ยนโปรแกรม ถ้า fit โมเดลเดียวกันด้วยสมมติฐานชุดเดียวกันใน R หรือใน Stata คอลัมน์ที่ซ้ำกันแบบเดิมก็จะกลับมา เพราะทั้งสองโปรแกรมกำลังคำนวณสมการที่ไม่มีอะไรจำเพาะกับตำแหน่งอยู่เลย มีเพียง between-person variance, residual variance และจำนวนผู้เข้าร่วม การซ้ำกันจึงเป็นเรื่องเลขคณิต ไม่ใช่บั๊กของการแสดงผล และมันยังเท่ากันอยู่ไม่ว่าจะพิมพ์ทศนิยมกี่ตำแหน่งก็ตาม

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

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

สามปริมาณที่มักถูกสับสน

ก่อนลุยของยาก: ปริมาณสามตัวข้างล่างจะไม่ปนกัน ถ้าผูกแต่ละตัวไว้กับคำถามทางคลินิกหนึ่งข้อ ผู้ป่วยแต่ละคนต่างกันแค่ไหน — นั่นคือ raw SD ค่าที่วัดครั้งหนึ่งอยู่ห่างจากระดับปกติของผู้ป่วยคนนั้นเองแค่ไหน — นั่นคือ residual SD แล้วเรามั่นใจแค่ไหนกับค่าเฉลี่ยที่เพิ่งประมาณได้ — นั่นคือ standard error มีเพียงตัวที่สามที่เล็กลงเมื่อเก็บผู้ป่วยเพิ่ม ซึ่งเป็นเหตุผลตรง ๆ ว่าทำไมมันจึงเท่ากันได้ในสี่ตำแหน่งที่การกระจายดิบไม่เท่ากัน

คำว่า "ความแปรปรวน" ในผลลัพธ์ของ mixed model ทำหน้าที่อย่างน้อยสามอย่าง และคำอธิบายที่สับสนระหว่างสามอย่างนี้จะพังทันทีที่มีคำถามที่สอง

หนึ่ง raw standard deviation ของแต่ละตำแหน่ง อธิบายว่าค่าที่วัดได้ที่ตำแหน่งนั้นกระจายรอบค่าเฉลี่ยของตำแหน่งมากแค่ไหน เมื่อดูข้ามผู้เข้าร่วมทุกคนในกลุ่ม ในข้อมูลจำลองชุดนี้ ค่าที่โมเดลบอกไว้เมื่อคุมค่า covariate ระดับบุคคลให้คงที่ คือ marginal SD เท่ากับ $\sqrt{\tau^2 + \sigma^2} = \sqrt{10.00} \approx 3.16$ g·m⁻²·h⁻¹ และขอให้สังเกตว่าค่านี้รวมการกระจายระหว่างบุคคลไว้ด้วย เพราะคนสองคนที่วัดตำแหน่งเดียวกันก็ต่างกันอยู่แล้วจริง ๆ ส่วน SD ที่คำนวณจากคอลัมน์ข้อมูลดิบตรง ๆ จะไม่ต่ำกว่าค่านี้ เพราะมันแบกความแปรปรวนที่ age, sex และ phototype อธิบายได้เอาไว้ด้วย

สอง residual standard deviation อธิบายสิ่งที่เหลืออยู่หลังจากตัด mean structure และ person effect ออกไปแล้ว คือตำแหน่งหนึ่งอยู่ห่างจากโปรไฟล์ของคนคนนั้นเองแค่ไหน ในที่นี้เท่ากับ $\sqrt{\sigma^2} = \sqrt{3.75} \approx 1.94$ g·m⁻²·h⁻¹ ส่วน between-person SD ที่ถูกตัดออกไปคือ $\sqrt{\tau^2} = 2.50$

สาม standard error ของค่าเฉลี่ยที่ประมาณได้ของแต่ละตำแหน่ง อธิบายว่าเราระบุตำแหน่งของค่าเฉลี่ยระดับ ประชากร ได้แม่นยำเพียงใดเมื่อมีขนาดตัวอย่างเท่านี้ ในดีไซน์นี้เท่ากับ 0.50 g·m⁻²·h⁻¹ ค่านี้ไม่ได้อธิบายตัวข้อมูลเลย แต่อธิบายตัวค่าประมาณ

ข้อความที่ไม่เป็นเหตุเป็นผลต่อกันควรพูดออกมาตรง ๆ เพราะเป็นขั้นตอนที่รับน้ำหนักของคำอธิบายทั้งหมด คือ $\text{SD}_{\text{shin}} \neq \text{SD}_{\text{forearm}}$ ไม่ได้ แปลว่า $\operatorname{SE}(\widehat{\mu}_{\text{shin}}) \neq \operatorname{SE}(\widehat{\mu}_{\text{forearm}})$ สองข้อความนี้พูดถึงคนละวัตถุ และข้อความแรกจะเดินทางไปหาข้อความหลังไม่ได้ถ้าไม่มีสมมติฐานที่คุณตรวจสอบได้มาช่วย

สมการที่ไม่มี index ของตำแหน่งอยู่เลย

ก่อนลุยของยาก: สมการถัดไปสั้นมาก และประเด็นของมันเป็นเรื่องที่ "มองเห็น" มากกว่าเรื่องพีชคณิต ลองดูฝั่งขวาแล้วหาว่ามีอะไรในนั้นที่บอกว่า forearm, hand หรือ shin บ้าง จะพบว่ามีเพียง between-person variance, residual variance และจำนวนผู้เข้าร่วม ไม่มีตำแหน่งอยู่ตรงไหนเลย ปริมาณที่สูตรของมันไม่เคยเอ่ยถึงตำแหน่ง ย่อมมีค่าต่างกันตามตำแหน่งไม่ได้

ภายใต้โครงสร้าง Spec A ของข้อมูลจำลอง คือ random intercept หนึ่งตัว pooled residual variance ค่าเดียว วัดครบสี่ตำแหน่งในทุกคน และมีผู้เข้าร่วมกลุ่มละ 40 คน standard error ของค่าเฉลี่ยที่ประมาณได้ของแต่ละตำแหน่งภายในกลุ่มหนึ่งมีรูปปิดดังนี้

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

ลองอ่านฝั่งขวาแล้วมองหาตำแหน่งที่วัด มันไม่อยู่ตรงนั้น สมการนี้มีเพียง between-person variance, residual variance และขนาดตัวอย่าง ส่วน index ของตำแหน่ง $j$ ไม่ปรากฏที่ใดเลย ส่วน contrast ระหว่าง AD กับ control ที่ตำแหน่งหนึ่งเป็นการเปรียบเทียบสองกลุ่มที่เป็นอิสระต่อกันกลุ่มละ 40 คน variance ของมันจึงเป็นผลบวกของสองปริมาณข้างต้น

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

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

อย่างไรก็ตาม ย่อหน้าข้างต้นกล่าวอ้างสองเรื่องพร้อมกัน และควรแยกออกจากกันก่อนที่จะมีคนมาแยกให้ เรื่องแรกคือ ความเท่ากัน ของ standard error ทั้งสี่ค่า ซึ่งเป็นความเท่ากันแบบเป๊ะ มันตามมาจากความสมมาตรของ design matrix ดังที่หัวข้อถัดไปจะแสดง และไม่เปลี่ยนแม้จะใส่ covariate ระดับบุคคลเข้าไปด้วย เรื่องที่สองคือ ค่ากลางร่วม ที่เท่ากับ 0.71 ซึ่งเป็นเพียงค่าประมาณแบบแทนค่าลงสูตร คำนวณจาก variance components ที่ใช้สร้างข้อมูลในโมเดลค่าเฉลี่ยรายเซลล์ที่ยังไม่ปรับ covariate ดังนั้นโมเดลที่ fit จริง ซึ่งประมาณ $\tau^2$ และ $\sigma^2$ จากข้อมูล และมี age, sex, phototype อยู่ใน mean structure จะพิมพ์ตัวเลขสี่ตัวที่ยังเท่ากันเองทุกประการ แต่ขยับออกจาก 0.71 ไปเล็กน้อย ข้อยกเว้นเดียวคือค่า 0.43 ของ contrast ภายในคนเดียวกันที่จะพิสูจน์ในหัวข้อถัดไป เพราะ covariate ที่คงที่ภายในคนหักล้างกันไปในผลต่างภายในคน การปรับ covariate จึงขยับค่านี้ไม่ได้เลย

รูปทั่วไป: เรขาคณิตของ contrast

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

รูปปิดข้างต้นเป็นกรณีเฉพาะเชิงเลขคณิต ส่วนข้อความรูปทั่วไปที่ควรใช้ตอบนักระเบียบวิธีคือเรื่องเรขาคณิตของ contrast เพราะ contrast เชิงเส้นใด ๆ ที่เขียนเป็น $\widehat{\Delta}_{s} = c_s'\widehat{\beta}$ มี

\[ \operatorname{SE}(\widehat{\Delta}_{s}) \;=\; \sqrt{c_s'\,\operatorname{Var}(\widehat{\beta})\,c_s}, \qquad \operatorname{Var}(\widehat{\beta}) = \left(\textstyle\sum_i X_i' V_i^{-1} X_i\right)^{-1} \]

contrast แบบ AD − control ทั้งสี่ตำแหน่งใช้โมเดลเดียวกันและใช้ $\operatorname{Var}(\widehat{\beta})$ ก้อนเดียวกัน เวกเตอร์ $c_{\text{forearm}}, c_{\text{hand}}, c_{\text{shin}}, c_{\text{back}}$ ต่างกันแค่ว่าเลือกคอลัมน์ indicator และ interaction คอลัมน์ไหนเท่านั้น เมื่อผู้เข้าร่วมทุกคนให้ข้อมูลหนึ่งค่าที่ทุกตำแหน่ง เมื่อสองกลุ่มมีขนาดเท่ากัน และเมื่อ $V_i$ เป็น compound-symmetric matrix ตัวเดียวกันสำหรับทุกคน ตำแหน่งที่ถูกเลือกเหล่านั้นจะอยู่ในจุดที่สมมาตรกันเชิงโครงสร้างภายใน $\operatorname{Var}(\widehat{\beta})$ quadratic form ทั้งสี่จึงเท่ากัน

\[ c_{\text{forearm}}'\operatorname{Var}(\widehat{\beta})\,c_{\text{forearm}} \;=\; c_{\text{hand}}'\operatorname{Var}(\widehat{\beta})\,c_{\text{hand}} \;=\; c_{\text{shin}}'\operatorname{Var}(\widehat{\beta})\,c_{\text{shin}} \;=\; c_{\text{back}}'\operatorname{Var}(\widehat{\beta})\,c_{\text{back}} \]

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

ส่วนประกอบสี่อย่าง

ความเท่ากันเกิดจากเงื่อนไขสี่อย่างที่ทำงานพร้อมกัน และแต่ละอย่างถอดออกได้ หนึ่ง pooled residual variance ค่าเดียว คือเราสั่งให้โมเดลประมาณ $\sigma^2$ เพียงค่าเดียวที่ทั้งสี่ตำแหน่งใช้ร่วมกัน จึงไม่มีตำแหน่งใดได้รับอนุญาตให้มีระดับ noise ของตัวเอง สอง ผู้เข้าร่วมกลุ่มเดียวกันถูกวัดครบทุกตำแหน่ง คือแต่ละคนให้ข้อมูลครบสี่แถว การเปรียบเทียบระหว่างตำแหน่งจึงตั้งอยู่บนคนกลุ่มเดียวกัน สาม ขนาดกลุ่มเท่าเดิมในทุกตำแหน่ง คือมีผู้เข้าร่วม AD 40 คนและ control 40 คนที่ทุกตำแหน่ง ค่า $n$ จึงเท่ากันทุกจุดที่มันปรากฏในสมการ สี่ covariate ที่คงที่ภายในคนคนเดียวกัน คือ age, sex และ phototype มีค่าเดียวกันในทั้งสี่แถวของผู้เข้าร่วมคนหนึ่ง

ส่วนประกอบข้อที่สามมักถูกพูดแรงเกินจริง รวมถึงในประโยคข้างต้นด้วย สิ่งที่ความเท่ากันต้องการคือขนาดของสองกลุ่มต้องไม่เปลี่ยนไปตามตำแหน่ง ไม่ได้ต้องการให้สองกลุ่มมีขนาดเท่ากันเอง ถ้ามี AD 30 คนและ control 50 คนที่ทุกตำแหน่ง standard error ทั้งสี่ค่าก็ยังเท่ากันอยู่ดี เพียงแต่ค่ากลางร่วมของทั้งสี่จะไม่ใช่ตัวเลขที่สูตรรูปปิดแบบ $n$ เดียวข้างต้นให้ไว้

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

ยังมีเงื่อนไขอีกหนึ่งข้อที่รองอยู่ใต้ทั้งสี่ข้อและมองข้ามได้ง่าย เพราะ Spec A ทำให้มันเป็นจริงอยู่แล้วโดยไม่ต้องประกาศ คือ marginal covariance ต้องเป็น compound-symmetric ตั้งแต่ต้น ซึ่งเป็นสิ่งที่ random intercept ตัวเดียวบวกกับ pooled residual variance ค่าเดียวให้มา ตอนที่ 5 คือตัวอย่างค้าน ถ้าเพิ่ม random slope เข้าไป เส้นทแยงมุมของ $V_i$ จะกลายเป็นฟังก์ชันของเวลา standard error ของแต่ละครั้งที่วัดจึงต่างกัน แม้ส่วนประกอบทั้งสี่ข้อข้างต้นจะเป็นจริงครบก็ตาม

ทำลายความสมดุลแล้วดูว่าเกิดอะไรขึ้น

ก่อนลุยของยาก: วิธีพิสูจน์ที่สะอาดที่สุดว่าความเท่ากันนั้นเป็นของจริง คือทำลายมันทิ้งโดยตั้งใจ ข้างล่างนี้ผู้เข้าร่วม 12 คนจะไม่มีค่าวัดที่ shin แล้ว standard error ของ shin ก็แยกตัวออกจากอีกสามตำแหน่งทันที จากนั้นเราจะเทียบวิธีคำนวณสองแบบ แบบไร้เดียงสาปฏิบัติกับ shin ราวกับว่ามันเป็นการศึกษาเดี่ยว ๆ ที่มีผู้ป่วย 28 คน ส่วนโมเดลทำได้ดีกว่านั้น เพราะผู้เข้าร่วม 12 คนนั้นยังให้ข้อมูลอีกสามตำแหน่งไว้ และ mixed model ใช้ค่าวัดอื่น ๆ ของคนคนหนึ่งมาบอกอะไรบางอย่างเกี่ยวกับค่าที่หายไป คล้ายกับที่คุณดูแขนขาอีกสามข้างของผู้ป่วยก่อนจะประเมินว่าข้างที่ไม่ได้บันทึกไว้น่าจะเป็นอย่างไร การยืมข้อมูลแบบนี้เองคือสิ่งที่ออฟเซ็ตที่ถูกแชร์ซื้อมาให้

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

แบบแรกคือการคำนวณแบบไร้เดียงสา ซึ่งปฏิบัติกับ shin ราวกับว่ามันอยู่ลำพัง แล้วผสมขนาดตัวอย่างสองค่าเข้าไปในสูตรของดีไซน์ที่สมดุล

\[ \operatorname{SE}_{\text{naive}}(\widehat{\Delta}_{\text{shin}}) \;\approx\; \sqrt{\frac{\tau^2+\sigma^2}{28} + \frac{\tau^2+\sigma^2}{40}} \;=\; \sqrt{\frac{10.00}{28} + \frac{10.00}{40}} \;=\; \sqrt{0.357+0.250} \;=\; \sqrt{0.607} \;\approx\; 0.78 \]

การ fit ด้วยวิธี likelihood ทำได้ดีกว่านั้น และดีกว่าด้วยเหตุผลที่ซีรีส์นี้สร้างมาทั้งเรื่อง เพราะ variance components ถูกใช้ร่วมกันข้ามตำแหน่ง ผู้เข้าร่วม 12 คนที่ไม่มีค่าที่ shin จึงไม่ได้หายไปจากค่าเฉลี่ยของ shin ทั้งหมด สามตำแหน่งที่เขาให้ข้อมูลไว้ผูกกับ shin ผ่าน $\tau^2$ และดึงสารสนเทศส่วนหนึ่งที่แถวที่หายไปพรากไปกลับคืนมาได้ ถ้าให้ $n_1 = 28$ คือผู้เข้าร่วมที่ข้อมูลครบ $n_2 = 12$ คือผู้ที่ข้อมูลไม่ครบ และ $N = n_1 + n_2 = 40$ วิธี generalised least squares ภายใต้ $V_i$ ที่สมมติไว้จะให้

\[ \operatorname{Var}(\widehat{\mu}_{\text{shin},\text{AD}}) \;=\; \frac{\tau^2+\sigma^2}{n_1} \;-\; \frac{3\,n_2\,\tau^4}{n_1 N (\sigma^2+3\tau^2)} \;=\; 0.3571 \;-\; 0.0558 \;=\; 0.3013 \]

พจน์ที่ถูกลบออกไปคือสารสนเทศที่ดึงกลับมาได้จากอีกสามตำแหน่งพอดี ส่วนเลข 3 ในตัวเศษก็คือจำนวนตำแหน่งที่เหลืออยู่นั่นเอง กลุ่ม control มีข้อมูลครบและให้ $10.00/40 = 0.25$ เหมือนเดิม standard error ของ contrast ที่ shin จึงเป็น

\[ \operatorname{SE}(\widehat{\Delta}_{\text{shin}}) \;=\; \sqrt{0.3013 + 0.25} \;=\; \sqrt{0.5513} \;\approx\; 0.74 \]

site	ผู้เข้าร่วมต่อกลุ่ม	contrast SE จากโมเดล	contrast SE จากสูตรไร้เดียงสา
forearm	40 / 40	0.71	0.71
hand	40 / 40	0.71	0.71
shin	28 / 40	0.74	0.78
back	40 / 40	0.71	0.71

ทั้งสองคอลัมน์เป็นค่ารูปปิดที่คำนวณจาก variance components ที่ใช้สร้างข้อมูล หน่วยเป็น g·m⁻²·h⁻¹ ทั้งคู่ให้ค่าตรงกันที่ทุกตำแหน่งที่ดีไซน์ยังสมดุล และแยกออกจากกันเฉพาะตำแหน่งที่ไม่สมดุลแล้วเท่านั้น

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

การสาธิตนี้มีข้อควรระวังสองข้อที่ต้องพูดตามตรง ข้อแรก ตัวเลขทั้งสองเป็นค่ารูปปิดที่คำนวณจาก variance components ที่ใช้สร้างข้อมูล ไม่ใช่จากค่าที่ถูกประมาณ ดังนั้นการ fit จริงซึ่งมี $\widehat{\tau}^{2}$, $\widehat{\sigma}^{2}$ และ covariate ระดับบุคคล จะให้ค่าที่ใกล้เคียงแต่ไม่ตรงเป๊ะ ข้อสอง ทั้งสองค่าตั้งอยู่บนสมมติฐานว่าค่าวัดที่ shin หายไปด้วยเหตุผลที่โมเดลมองเห็นได้ ถ้าการหายไปขึ้นกับข้อมูลที่บันทึกไว้เท่านั้น เช่น กลุ่มของผู้เข้าร่วม covariate ของเขา และค่าที่อีกสามตำแหน่ง การ fit ด้วยวิธี likelihood ก็ยังใช้ได้อยู่ นั่นคือกรณี missing at random แต่ถ้าการหายไปขึ้นกับค่าที่ shin ที่ไม่ได้ถูกบันทึกนั้นเอง คือคนที่ผิวแห้งที่สุดเป็นคนที่ไม่ได้ถูกวัดพอดี ก็ไม่มี covariance structure ใดซ่อมมันได้ กรณีหลังนี้ต้องอาศัยสมมติฐานที่ระบุออกมาชัด ๆ พร้อม sensitivity analysis ของมันเอง

วิธีที่สองในการทำลายความเท่ากันคือคงความสมดุลไว้อย่างสมบูรณ์ แล้วไปโจมตี pooled variance แทน ในข้อมูลจำลอง Spec B ของตอนที่ 6 ทั้งสี่ตำแหน่งมี residual variance ต่างกัน marginal variance $\tau^2+\sigma_j^2$ จึงกลายเป็น 8.00, 9.00, 20.00 และ 10.00 แทนที่จะเป็น 10.00 เท่ากันทุกตำแหน่ง เมื่อใช้สูตร contrast เดิมแต่ใส่ variance ของแต่ละตำแหน่งลงไป ค่าทั้งสี่ก็แยกออกจากกันทันที ที่ตำแหน่ง shin ซึ่งมี marginal variance เท่ากับ 20.00

\[ \operatorname{SE}(\widehat{\Delta}_{\text{shin}}) \;\approx\; \sqrt{\frac{2(\tau^2+\sigma^2_{\text{shin}})}{n}} \;=\; \sqrt{\frac{2\times 20.00}{40}} \;=\; \sqrt{1.00} \;=\; 1.00 \]

ส่วนที่ตำแหน่ง back ซึ่งมี marginal variance เท่ากับ 10.00 การคำนวณแบบเดียวกันให้ $\sqrt{2\times 10.00/40}=\sqrt{0.50}\approx 0.71$ และ forearm กับ hand ซึ่งเป็นสองตำแหน่งที่เงียบที่สุดในข้อมูลจำลองชุดนั้นจะอยู่ต่ำกว่าทั้งคู่ ดีไซน์เดียว ข้อมูลสมดุลชุดเดียว แต่ standard error ออกมาสี่ค่าไม่เท่ากัน

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

SE ที่เท่ากันไม่ได้แปลว่าเป็นค่าประมาณสี่ค่าที่อิสระต่อกัน

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

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

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

ค่าออฟเซ็ตของคนคนนั้นปรากฏในทั้งสองพจน์ด้วยเครื่องหมายเดียวกันและหักล้างกันพอดี เหลือไว้เพียง residual สองตัว standard error ของผลต่างนี้จึงตัด between-person variance ทิ้งไปทั้งก้อน

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

ตัวเลขเดียวกันนี้เดินทางมาถึงได้ด้วยเส้นทางที่สอง ซึ่งควรแสดงไว้เพราะมันทำให้ความไม่อิสระมองเห็นได้ ค่าเฉลี่ยของแต่ละตำแหน่งมี variance เท่ากับ $0.25$ และค่าเฉลี่ยสองตำแหน่งภายในกลุ่มเดียวกันมี covariance เท่ากับ $\tau^2/n = 6.25/40 = 0.15625$ เพราะทั้งคู่คำนวณจากผู้เข้าร่วมชุดเดียวกัน ดังนั้น

\[ \operatorname{Var}(\widehat{\mu}_{\text{forearm}} - \widehat{\mu}_{\text{back}}) \;=\; 0.25 + 0.25 - 2(0.15625) \;=\; 0.1875, \qquad \sqrt{0.1875} \approx 0.43 \]

ทีนี้ลองเทียบตัวเลขสองตัวที่โมเดลรายงาน คือ 0.43 สำหรับผลต่างระหว่างตำแหน่งภายในคนเดียวกัน กับ 0.50 สำหรับค่าเฉลี่ยของแต่ละตำแหน่งเดี่ยว ๆ ผลต่างของสองปริมาณถูกประมาณได้ แม่นยำกว่า แต่ละปริมาณเสียอีก ซึ่งเป็นไปไม่ได้ทางเลขคณิตหากค่าเฉลี่ยสองตัวนั้นอิสระต่อกัน และเป็นลายนิ้วมือโดยตรงของออฟเซ็ตที่ถูกแชร์ โมเดลที่ปฏิบัติกับค่าวัดทั้งสี่ราวกับว่าเป็นอิสระต่อกันผลิตตัวเลขนี้ไม่ได้ และการเทียบคู่นี้มักเป็นวิธีที่เร็วที่สุดในการแสดงว่า random intercept ทำงานจริง

ความผิดพลาดด้านตรงข้าม ซึ่งเป็นส่วนที่คำอธิบายส่วนใหญ่พลาด

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

1. "SE ที่เท่ากันสี่ค่าแปลว่ามีบั๊ก" ผิด เพราะมันคือผลทางเลขคณิตของ balanced design ที่มี pooled residual variance ค่าเดียว
2. "ในเมื่อ SE เป็นผลลัพธ์ที่ถูกต้องตามกระบวนการ สมมติฐาน pooled variance ก็ต้องถูกต้อง" ผิดเท่า ๆ กัน เพราะโปรแกรมคำนวณความไม่แน่นอนโดยมีเงื่อนไขบน covariance structure ที่เราส่งให้มัน และมันจะทำหน้าที่นั้นอย่างว่าง่ายเหมือนเดิมแม้โครงสร้างนั้นจะผิด

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

เมื่อไรที่ SE เท่ากันควรทำให้กังวลจริง ๆ

ความเท่ากันเป็นสิ่งที่คาดหวังได้ภายใต้ส่วนประกอบสี่อย่างข้างต้น มันจะกลายเป็นเหตุผลให้ต้องตรวจสอบให้หนักขึ้นก็ต่อเมื่อข้อมูลขัดแย้งกับสมมติฐานที่ผลิตมันขึ้นมา มีสัญญาณห้าข้อที่ควรใส่ใจ หนึ่ง การกระจายของ residual แตกต่างกันชัดเจนระหว่างหน่วยที่วัดใน diagnostics เช่น box plot ของ residual แยกตามตำแหน่งที่ความกว้างไม่อาจถือว่าใช้แทนกันได้ สอง กราฟ residual เทียบกับ fitted แสดง heteroscedasticity ที่มีรูปแบบ เช่น รูปกรวยที่บานออกตามค่า fitted สาม ข้อมูลขาดหายมากและไม่เท่ากันระหว่างหน่วยที่วัด ในขณะที่ standard error ยังเท่ากันเป๊ะ ซึ่งบ่งชี้ว่าความเท่ากันถูกค้ำไว้ด้วยสมมติฐานมากกว่าด้วยข้อมูล สี่ ค่าประมาณ variance ที่ชนขอบหรือเป็น singular ซึ่ง variance component ยุบลงเป็นศูนย์ และโครงสร้างที่ถูกรายงานไม่ใช่โครงสร้างที่ข้อมูลรองรับ [1] ห้า มีทางเลือกแบบ heteroscedastic ที่มีเหตุผลรองรับและทำให้ confidence interval หรือข้อสรุปเชิงเนื้อหาเปลี่ยนไปอย่างมีนัยสำคัญ ซึ่งเป็นข้อเดียวในห้าข้อที่ชี้ขาด เพราะมันแสดงว่าสมมติฐานไม่ได้แค่น่าสงสัย แต่ส่งผลจริง

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

วิธีที่ถูกต้องในการยุติเรื่องนี้: sensitivity analysis ที่ระบุไว้ล่วงหน้า

มีคำตอบเดียวที่ยุติบทสนทนานี้ได้อย่างเหมาะสม และมันไม่ใช่การโต้เถียง คือ refit โมเดลด้วย residual variance แยกตามหน่วยที่วัด โดยคง fixed-effect structure ไว้เหมือนเดิมทุกประการ แล้วรายงานว่าอะไรเปลี่ยนไปบ้าง [2]

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

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

มีสามข้อเชิงกระบวนการที่ทำให้สิ่งนี้เป็น sensitivity analysis ไม่ใช่การตกปลาหาผลลัพธ์ หนึ่ง การเปรียบเทียบต้องใช้ fixed effects ชุดเดียวกัน เพราะค่า likelihood แบบ REML จาก mean structure ที่ต่างกันเทียบกันไม่ได้ สอง การวิเคราะห์ควรถูกระบุไว้ล่วงหน้าใน protocol หรือใน statistical analysis plan เพื่อไม่ให้ถูกอ่านว่าเป็นปฏิกิริยาต่อค่า p ที่ไม่พึงประสงค์ สาม การเปลี่ยนโครงสร้างที่ไม่ได้ระบุไว้ล่วงหน้าไม่ได้แปลว่าไม่ชอบธรรม แต่ต้องติดป้ายให้ตรงกับสิ่งที่มันเป็น โครงสร้างที่เปลี่ยนเพราะ residual diagnostics ชี้ไปทางนั้น ให้รายงานว่าเป็นการวิเคราะห์แบบ post hoc ที่ถูกกระตุ้นด้วย diagnostics ส่วนโครงสร้างที่เปลี่ยนเพราะค่า p ของ fixed effects ขยับ ไม่ใช่ sensitivity analysis ในรูปแบบใดเลย ไม่ว่าจะเรียกมันว่าอะไรในต้นฉบับก็ตาม

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

As a pre-specified sensitivity analysis, the model was refitted with site-specific residual variances and an identical fixed-effect structure. The estimated site × group differences changed by at most [___] g·m⁻²·h⁻¹, the width of every 95% confidence interval changed by less than [___], and the direction and statistical significance of every reported contrast were preserved.

ถ้ามีบางอย่างเปลี่ยนไปจริง ให้พูดถึงมันในประโยคเดียวกับที่คุณอ้างความคงทนของผลลัพธ์

As a pre-specified sensitivity analysis, the model was refitted with site-specific residual variances and an identical fixed-effect structure. The principal comparison at [___] was preserved in both magnitude and direction, whereas the comparison at [___] became less precise (95% CI widened from [___] to [___]) and no longer met the conventional threshold of p < 0.05. This secondary comparison is therefore reported as imprecise and hypothesis-generating.

อะไรที่ไม่ใช่เหตุผลทางสถิติ

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

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

สคริปต์สี่ขั้น

นี่คือคำตอบที่ใช้พูดออกมาได้จริง ใช้เวลาประมาณหนึ่งนาที

1. เป็นบั๊กหรือเปล่า ไม่ใช่ ภายใต้ balanced design ที่มี pooled residual variance ค่าเดียว standard error ของ contrast ของแต่ละตำแหน่งคือ $\sqrt{2(\tau^2+\sigma^2)/n}$ และไม่มี index ของตำแหน่งอยู่ในสมการนั้นเลย ค่าทั้งสี่เท่ากันเพราะพีชคณิตบังคับให้เป็นเช่นนั้น

2. ทำไมถึงเท่ากันเป๊ะ เพราะเงื่อนไขสี่ข้อเป็นจริงพร้อมกัน คือ pooled residual variance ค่าเดียว ผู้เข้าร่วมกลุ่มเดียวกันถูกวัดครบทุกตำแหน่ง ขนาดกลุ่มไม่เปลี่ยนไปตามตำแหน่ง และ covariate ที่คงที่ภายในคนคนเดียวกัน ถอดข้อใดข้อหนึ่งออก ความเท่ากันก็ไม่ถูกรับประกันอีกต่อไป เช่น ถ้าค่าวัดที่ shin หายไป 12 คน contrast SE ของตำแหน่งนั้นจะขยับจาก 0.71 ไปเป็นประมาณ 0.74 เมื่อโมเดลยืมอีกสามตำแหน่งกลับมา และจะเป็นประมาณ 0.78 ก็ต่อเมื่อคิดตำแหน่งนั้นแยกลำพังด้วยสูตรไร้เดียงสาเท่านั้น

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

4. ถ้าผ่อนสมมติฐานนั้นแล้วผลจะกลับด้านหรือไม่ นี่คือ sensitivity analysis ที่ระบุไว้ล่วงหน้า ใช้ residual variance แยกตามตำแหน่งบน fixed effects ชุดเดิม และนี่คือสิ่งที่เปลี่ยนและไม่เปลี่ยนอย่างชัดเจน

ขั้นที่ 1 และ 2 เป็นพีชคณิตและพูดแบบยืนยันได้เลย ส่วนขั้นที่ 3 และ 4 คือจุดที่ความน่าเชื่อถือถูกสร้างขึ้นจริง

checklist ขั้นต่ำสำหรับการรายงาน mixed model

ผู้อ่านควรประกอบ covariance structure ของคุณขึ้นมาใหม่ได้จากตัวบทความเพียงอย่างเดียว โดยปกติเจ็ดข้อก็เพียงพอ หนึ่ง mean structure รวมทุก interaction ที่ fit ไว้ ไม่ใช่เฉพาะตัวที่มีนัยสำคัญ สอง โครงสร้างของ random effects พร้อมระบุให้ชัดว่า grouping factor คืออะไร สาม สมมติฐานของ residual covariance ว่าเป็นแบบ pooled แยกตามหน่วย หรือ unstructured โดยเรียกชื่อออกมาตรง ๆ ไม่ใช่ปล่อยให้เดา สี่ วิธีประมาณค่า (REML หรือ ML) และวิธีประมาณ degrees of freedom ที่ใช้ในการทดสอบ ห้า ค่าประมาณของ variance components และ ICC ที่ตามมา โดยค่าที่ใช้สร้างตัวอย่างจำลองชุดนี้คือ $\tau^2 = 6.25$, $\sigma^2 = 3.75$ และ $\rho = 0.625$ ส่วนการ fit ด้วย REML จริงจะให้ค่าประมาณที่ใกล้เคียงค่าเหล่านี้ ไม่ใช่เท่ากันเป๊ะ ซึ่งเป็นเหตุผลตรง ๆ ว่าทำไมสิ่งที่ควรอยู่ในบทความคือค่าประมาณ ไม่ใช่ค่าที่ใช้ออกแบบข้อมูล หก โปรแกรมและเวอร์ชัน พร้อมคำสั่งที่ใช้ fit หากรูปแบบของวารสารเปิดให้ทำได้ และสุดท้าย sensitivity analysis ของ covariance structure ที่ระบุไว้ล่วงหน้าพร้อมผลของมัน

หมายเหตุหนึ่งย่อหน้าเรื่อง degrees of freedom

denominator degrees of freedom เป็นคำถามที่ยังเปิดอยู่จริงใน mixed model ไม่ใช่เรื่องที่ถูกมองข้าม และซีรีส์นี้จงใจกันมันออกไปนอกขอบเขต ในดีไซน์ที่สมดุลและมี grouping factor เพียงชั้นเดียว วิธีประมาณที่ใช้กันทั่วไป ได้แก่ containment, Satterthwaite และ Kenward–Roger มักให้ผลใกล้เคียงกันมาก และการเลือกวิธีแทบไม่เปลี่ยนข้อสรุป แต่เมื่อข้อมูลไม่สมดุล จำนวน cluster น้อย หรือโครงสร้าง random ซับซ้อน ทั้งสามวิธีอาจแยกออกจากกันมากพอที่จะมีผล และ Kenward–Roger อาจเหมาะกว่าเป็นพิเศษเพราะปรับค่าประมาณ variance ให้รองรับความไม่แน่นอนของ covariance parameters ด้วย กฎการรายงานง่ายกว่าทฤษฎีมาก คือให้ระบุว่าใช้วิธีประมาณแบบใด เพราะค่า t ที่ไม่มี denominator degrees of freedom กำกับนั้นทำซ้ำไม่ได้

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

มีการอ่านเกินสามแบบที่ควรปฏิเสธอย่างชัดเจน รวมถึงรูปที่แรงเกินไปของข้อสรุปของบทความนี้เองด้วย

บทความนี้ไม่ได้อนุญาตให้สรุปว่า compound symmetry เพียงพอโดยทั่วไป ซีรีส์นี้แสดงไปแล้วสองครั้งว่ามันไม่เพียงพอ คือ random slope ทำลายมันด้วยการทำให้เส้นทแยงมุมกลายเป็นฟังก์ชันของเวลา และ residual variance ที่ไม่เท่ากันทำลายความเท่ากันของ correlation ทั้งที่ covariance ยังเท่ากันทุกช่อง คำถามที่ว่า pooled residual variance มีเหตุผลรองรับหรือไม่ เป็นคำถามเกี่ยวกับ diagnostics ดีไซน์ และคำถามวิจัยของคุณ และต้องตอบใหม่ทุกครั้ง

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

สุดท้าย บทความนี้ไม่ได้อนุญาตให้เชื่อว่า covariance structure เลือกได้จากการดูว่าค่า p ของ fixed effects เปลี่ยนไปอย่างไร กระบวนการแบบนั้นทำลายความถูกต้องของ inference ที่มันถูกใช้ปกป้อง โครงสร้างต้องถูกเลือกจาก residual diagnostics จากดีไซน์ และจากการเปรียบเทียบโมเดล variance ที่ซ้อนกันด้วยวิธี likelihood บน fixed effects ชุดเดียวกัน

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

ห้าอย่าง เรียงตามลำดับที่มักต้องใช้จริง

1. ก่อนอื่นให้ plot residual แยกตามหน่วยที่วัด ถ้าการกระจายดูใช้แทนกันได้ การคง pooled residual variance ไว้ก็อาจมีเหตุผลเพียงพอ แต่ถ้าหน่วยใดหน่วยหนึ่ง noisy กว่าอย่างเห็นได้ชัด ข้อสังเกตนั้นควรอยู่ในแผนการวิเคราะห์ ไม่ใช่ในส่วนอภิปราย
2. เขียน covariance structure ลงใน protocol พร้อมกับ sensitivity analysis ที่จะใช้ตรวจสอบมัน โครงสร้างที่เลือกไว้ล่วงหน้าคือการตัดสินใจ ส่วนโครงสร้างเดียวกันที่เลือกหลังเห็นค่า p แล้วไม่ใช่
3. รายงาน SD และ SE แยกกันโดยติดป้ายกำกับทั้งคู่ และอย่าปล่อยให้คอลัมน์ใดในตารางกำกวมระหว่างสองค่านี้ ผู้อ่านกู้ความแตกต่างคืนจากตัวเลขเปล่า ๆ ไม่ได้
4. รายงาน variance components และ ICC ที่ตามมา เพราะมันทำให้ผู้อ่านประกอบ design effect กลับขึ้นมาได้ และตัดสินได้ว่าชุดข้อมูลนี้มีสารสนเทศที่เป็นอิสระอยู่จริงมากแค่ไหน
5. เมื่อโครงสร้างสองแบบให้ผลไม่ตรงกัน ให้รายงานทั้งสองแบบ ข้อสรุปที่รอดจากการระบุโมเดลทางเลือกที่มีเหตุผลรองรับมีค่ามากกว่าข้อสรุปที่ไม่เคยถูกทดสอบ ส่วนข้อสรุปที่ไม่รอดก็ควรถูกอธิบายตามที่มันเป็น

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

• ภายใต้ balanced design ที่มี pooled residual variance ค่าเดียว standard error ของ contrast ของทุกหน่วยคือปริมาณเดียวกันที่ถูกคำนวณซ้ำ ๆ เพราะสมการ $\sqrt{2(\tau^2+\sigma^2)/n}$ ไม่มี index ของหน่วยอยู่เลย
• raw SD, residual SD และ SE ของค่าเฉลี่ยที่ประมาณได้ เป็นสามปริมาณที่ต่างกัน และการที่สองหน่วยมีการกระจายไม่เท่ากันไม่ได้แปลว่า standard error ของค่าเฉลี่ยทั้งสองจะไม่เท่ากัน
• การถอดส่วนประกอบข้อใดข้อหนึ่งในสี่ข้อออก ได้แก่ pooled residual variance, การวัดครบทุกหน่วยในทุกคน, ขนาดกลุ่มที่ไม่เปลี่ยนไปตามหน่วยที่วัด และ covariate ที่คงที่ภายในคนคนเดียวกัน ทำให้ความเท่ากันไม่ถูกรับประกันอีกต่อไป ดังที่ contrast ของ shin แสดงให้เห็น คือเมื่อมีค่าวัดหายไป 12 คน มันขยับจาก 0.71 ไปเป็นประมาณ 0.74 ภายใต้การ fit ด้วยวิธี likelihood และเป็นประมาณ 0.78 ภายใต้สูตรไร้เดียงสาที่คิดตำแหน่งนั้นแยกลำพังโดยไม่นับการยืมข้อมูล
• SE ที่เท่ากันไม่ได้แปลว่าเป็นค่าประมาณสี่ค่าที่อิสระต่อกัน เพราะ contrast ภายในคนเดียวกันมี SE เท่ากับ 0.43 เทียบกับ 0.50 ของค่าเฉลี่ยเดี่ยว ๆ เนื่องจากออฟเซ็ตที่ถูกแชร์หักล้างกันไปในผลต่าง
• ผลลัพธ์ที่ถูกต้องตามกระบวนการไม่เคยยืนยันสมมติฐานที่ผลิตมันขึ้นมา คำตอบจริงของคำถามนี้จึงมีเพียงอย่างเดียว คือ sensitivity analysis ที่ระบุไว้ล่วงหน้าด้วย residual variance แยกตามหน่วย พร้อมรายงานว่าอะไรเปลี่ยนและอะไรไม่เปลี่ยน

ที่ที่ซีรีส์นี้จบลง: หนึ่งค่า จากตารางแรกถึงคอลัมน์สุดท้าย

ซีรีส์นี้ตามปริมาณเดียวมาตั้งแต่ต้นจนจบ ในตอนที่ 1 มันมองเห็นได้ด้วยตาเปล่า คือผู้เข้าร่วมจำลองคนหนึ่งอยู่เหนือโปรไฟล์ของประชากรในทั้งสี่ตำแหน่ง ส่วนอีกคนอยู่ต่ำกว่าในทั้งสี่ตำแหน่ง เป็นออฟเซ็ตที่คุณเห็นได้ก่อนจะ fit โมเดลใด ๆ ในตอนที่ 2 ออฟเซ็ตนั้นเติม covariance matrix ให้มีค่านอกเส้นทแยงมุมเท่ากันหมด และได้ชื่อของตัวเองว่า intraclass correlation ที่ 0.625 ในตอนที่ 3 มันแยกโมเดลตัวเดียวกันออกเป็นสองการอ่านที่ถูกต้องทั้งคู่ คือเป็นอิสระต่อกันเมื่อรู้ตัวคน และ correlate กันเมื่อมองข้ามประชากร ในตอนที่ 4 มันกลายเป็นผลคูณเมทริกซ์ $Z_i G Z_i'$ พร้อมกับแผนที่ที่บอกว่ามีเพียงสองสิ่งเท่านั้นที่เปลี่ยนได้ ตอนที่ 5 และ 6 เปลี่ยนทีละอย่าง และตอนที่ 7 ปฏิเสธการแยกส่วนทั้งหมด ส่วนในตอนสุดท้ายนี้ กลไกเดิมตัวนั้น คือออฟเซ็ต $b_i$ ที่คนหนึ่งพกติดตัว ความแปรปรวน $\tau^2$ ของมันเมื่อมองข้ามประชากร และ correlation $\rho$ ที่มันสร้างขึ้นระหว่างสองแถวใดก็ตามของคนคนเดียวกัน คือสิ่งที่ทำให้ standard error สี่ค่าเท่ากัน และทำให้ contrast ภายในคนเดียวกันเล็กกว่าค่าเฉลี่ยที่ประกอบมันขึ้นมา

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

ซีรีส์นี้ตามลำดับ:

1. ค่าชดเชยส่วนตัวที่ถูกแชร์: Random Intercept คืออะไรกันแน่
2. Compound Symmetry กับ ICC: ออฟเซ็ตเดียว correlation เดียว
3. Conditional หรือ Marginal? สองคำถามในโมเดลเดียว
4. V = ZGZ′ + R: สมการเดียวที่อยู่เบื้องหลัง mixed model ทุกตัว
5. Random Slope: เมื่อแต่ละคนพก "เส้นทาง" ของตัวเองมาด้วย
6. เมื่อ Compound Symmetry ไม่จริง: variance ไม่เท่ากัน, unstructured correlation และเพดาน identifiability
7. ทางเลือกแบบ Population-Averaged: GLS, GEE และจุดที่โชคของ linear model หมดลง
8. SE เท่ากันสี่ค่า: วิธีปกป้อง mixed model ต่อหน้าอาจารย์ — คุณอยู่ที่นี่

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

1. Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. J Stat Softw. 2015;67(1):1–48.
2. Verbeke G, Molenberghs G. Linear Mixed Models for Longitudinal Data. New York: Springer; 2000.